Interfacial order parameter asymmetry correction based simulation method for seawater freezing desalination across scales

Abstract

The application provides a seawater freezing desalination cross-scale simulation method based on interface order parameter asymmetric correction, and belongs to the technical field of seawater desalination based on computer data processing; first, a sea ice-seawater two-phase microscopic molecular dynamics model is constructed, interface microscopic structure features are extracted, then a quantitative asymmetric correction function is constructed, a mesoscopic phase field model with asymmetric correction is established based on the quantitative asymmetric correction function, multi-physical field dynamic evolution is coupled and solved, finally, a cross-scale simulation model is generated and a prediction result is output; the application effectively solves the morphology distortion problem of a traditional phase field model when simulating high salinity crystallization, and provides a very reliable digital design tool for supercooling degree adjustment and refrigeration rate optimization in an actual seawater freezing desalination process.

Description

Technical Field

[0001] This invention belongs to the field of seawater desalination technology based on computer data processing, and particularly relates to a cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters. Background Technology

[0002] Seawater cryo-desalination technology has become an important development direction in seawater desalination due to its low energy consumption and environmental friendliness. Seawater cryo-desalination utilizes the phase equilibrium relationship between water and salt under freezing point conditions to achieve separation. Under moderately supercooled conditions, water crystallizes out of seawater, while salt ions are repelled to the solid-liquid interface and aggregate. These high-concentration salts easily become trapped inside ice crystals, forming "salt cells," thus significantly reducing seawater desalination efficiency. To explore the underlying mechanisms for improving cryo-desalination rates, researchers typically employ numerical simulation techniques such as macroscopic computational fluid dynamics (CFD), mesoscopic phase field method (PF), and microscopic molecular dynamics (MD). These numerical simulation techniques can reveal ice crystal growth dynamics and solute redistribution mechanisms at different spatiotemporal scales.

[0003] Currently, researchers mostly use phase-field models to simulate the evolution of mesoscopic ice dendrites during seawater crystallization, including mesoscopic-scale simulations and microscopic-scale simulations. Existing simulation methods have the following main problems when predicting the efficiency of seawater freezing and desalination:

[0004] (1) Inaccurate physical description: For the sake of mathematical simplicity, the traditional phase-field method usually assumes that the order parameter is symmetrically distributed on both sides of the solid-liquid interface. However, molecular dynamics reveals that the molecular order at the interface has significant asymmetric characteristics. This "mathematical simplification assumption" leads to insufficient accuracy of the traditional model in describing the nonlinear driving force under high salinity seawater.

[0005] (2) Model prediction distortion: Due to the neglect of the asymmetric physical properties of the interface, traditional models are difficult to accurately predict how salt in seawater is captured by ice crystals, especially in the later stage of dendrite differentiation in seawater. The calculated Pelet number and other key dynamic indicators have large deviations from the experimental results, with relative errors as high as 240% or more.

[0006] (3) Lack of cross-scale mechanism: At present, there is no coupling mechanism in the industry that can effectively map the interface asymmetric features at the atomic scale to the mesoscopic dendrite evolution, which hinders the synergistic prediction from molecular evolution behavior to macroscopic dilution effect.

[0007] The main reason for these shortcomings is that seawater crystallization is a complex, multi-scale physical process. Traditional methods fail to fully utilize information at the microscale to correct the free energy function at the mesoscale. Furthermore, atomic-scale simulations are limited in their spatiotemporal scale and lack quantitative asymmetric correction methods. These factors prevent existing numerical simulation platforms from accurately reflecting the coupling between salt repulsion and interface development during seawater freezing, ultimately hindering researchers from accurately evaluating and optimizing seawater freezing and desalination processes. Summary of the Invention

[0008] To address the above problems, this invention proposes a cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters, comprising the following steps:

[0009] S1, Select the water molecule force field and electrolyte force field that can characterize the phase transition characteristics to construct a two-phase model of sea ice-seawater, perform relaxation calculations in the equilibrium state at the preset melting point temperature, and output the motion trajectory data of water molecules and salt ions in the equilibrium state.

[0010] S2, calculate the bond order parameters of each water molecule in equilibrium, fit the spatial distribution of the bond order parameter values, and statistically output the spatial distribution curve of water molecule number density and the spatial distribution data of salt ion concentration in the interface transition region.

[0011] S3, based on the output of S2, uses the hyperbolic tangent function to fit and obtain the asymmetric normalized order parameter distribution equation, and establishes the exponential correlation of the salt ion concentration gradient with the order parameter evolution, generating the exponential asymmetric correction function of the salt ion concentration with the order parameter evolution.

[0012] S4. The exponential asymmetric correction function is used as a correction term and embedded into the free energy density functional of the phase field model. After partial derivative calculation, the phase field evolution control equation containing the asymmetric chemical driving force term is generated.

[0013] S5, based on the exponential asymmetric correction function and the phase field evolution control equation, generates a cross-scale simulation model for seawater freezing and desalination, and performs simulation calculations of the seawater freezing and desalination process.

[0014] Preferably, the specific process for simulating the seawater freezing and desalination process is as follows:

[0015] Based on Fick's second law, the concentration field equation is derived, and the temperature field equation is established based on Fourier's law. The grid step size and time step size of the mesoscopic simulation are obtained. The control equations of the phase field, concentration field, and temperature field are spatiotemporally discretized, random thermal perturbation terms are added, and numerical iteration is performed to output the dynamically evolving phase field, concentration field, and temperature field distribution matrices. Based on the output distribution matrices, dendritic morphology evolution characteristic data and crystallization desalination rate data of the seawater freezing desalination process are extracted and output.

[0016] Preferably, in step S1, a water molecule force field and an electrolyte force field characterizing the phase transition properties of water molecules are selected, and salt ion pairs are inserted into the liquid phase region to configure a seawater liquid phase with a preset salinity. Then, a solid ice substrate is spliced ​​with the seawater liquid phase to generate a sea ice-seawater two-phase model with a preset spatial size. Subsequently, the equilibrium relaxation process controls the system temperature at a preset melting point and maintains the pressure at a preset pressure. At the same time, periodic boundary conditions are applied to the microscopic model, and the monitoring program monitors the convergence state of the total energy and density of the system in real time. When the solid and liquid phases reach a dynamic coexistence equilibrium, the data on the motion trajectories of generated water molecules and salt ions are output.

[0017] Preferably, the specific process of S2 is as follows:

[0018] Multiple candidate cutoff radii are set, and the bond order parameter values ​​of water molecule coordinate data at the corresponding cutoff radii are calculated. Probability density distribution curves at different cutoff radii are generated. By comparing multiple sets of probability density distribution curves, the cutoff radius with the best discrimination is selected as the statistical cutoff radius. The intersection point of the solid-liquid phase distribution curves at the optimal statistical cutoff radius is extracted, and the critical threshold for phase differentiation is set accordingly. Water molecules with bond order parameters greater than the critical threshold are marked as solid-phase ice molecules, and water molecules with bond order parameters less than the critical threshold are marked as liquid-phase water molecules.

[0019] Next, the number of water molecules in each tiny statistical region is counted, and the spatial distribution data of the average bond order parameter along the freezing direction is calculated. The spatial distribution data of the average bond order parameter is fitted using the hyperbolic tangent function to extract and generate the solid-liquid interface position coordinate data. Subsequently, based on the solid-liquid interface position coordinate data, the spatial distribution matrix of water molecule and salt ion number density in the two-phase model is calculated. The number density baseline value of the liquid phase region is extracted, and the number density gradient transition region between the solid and liquid phases is extracted, thereby calibrating the physical width of the interface transition region. Finally, the spatial distribution data of water molecule number density and salt ion concentration in the interface transition region are exported.

[0020] Preferably, in step S3, based on the spatial distribution curve of water molecule number density output from S2, a hyperbolic tangent function combined with a nonlinear quadratic term is used for fitting to generate an asymmetric normalized order parameter distribution equation; simultaneously, based on the spatial distribution data of salt ion concentration output from S2, an exponential fitting is performed to generate an exponential asymmetric correction function for the evolution of salt ion concentration with order parameters; specifically:

[0021] The average bond sequence parameter spatial distribution data is used, and the average bond sequence parameter spatial distribution data is normalized to generate normalized sequence parameter distribution data. The normalized sequence parameter distribution data at multiple time points is extracted, and the time-averaged bond sequence parameter distribution data is calculated and generated. The time-averaged bond sequence parameter distribution data is spatially fitted using the hyperbolic tangent function combined with nonlinear quadratic terms to generate the normalized sequence parameter asymmetric spatial distribution equation.

[0022] Subsequently, spatial distribution data of salt ion concentration is extracted, and the concentration data of the included anions and cations are extracted separately. The average value of the two is calculated to generate overall local concentration distribution data of the water-salt binary solution. The concentration distribution data is mapped to the generated normalized order parameter asymmetric distribution equation to generate a mapping dataset of salt ion concentration changing with the order parameter. Based on the mapping dataset, exponential fitting is performed and an exponential asymmetric correction equation is output.

[0023] Preferably, the specific implementation process of S4 is as follows:

[0024] The exponential asymmetric correction equation is used as a microscopic asymmetric correction term and embedded into the chemical free energy density functional of the basic phase field model. Subsequently, the equation is derived, and the partial derivative of the free energy density functional equation is calculated during the derivation process to extract and generate asymmetric chemical driving force partial derivative data. Anisotropic parameters between the solid and liquid interfaces are set and introduced to obtain the asymmetric phase field evolution control equation. Random thermal perturbation matrix data related to the interface position are superimposed on the equation iteration terms.

[0025] Preferably, Fick's second law and the condition of equal chemical potential are introduced to construct the concentration field evolution control equation. At the same time, based on the law of thermal conduction and combined with the latent heat parameter of crystallization phase transformation, the temperature field evolution control equation is established. The two evolution control equations contain specific physical parameter characteristics, including solute diffusion coefficient, thermal diffusion coefficient, latent heat release constant of solidification, solid-liquid phase volume fraction derivative, and solution specific heat capacity parameter. Subsequently, the two control equations are merged with the output asymmetric phase field evolution control equation to finally generate a multiphysics field coupled dynamic equation set.

[0026] The phase field evolution control equation, concentration field evolution control equation, and temperature field evolution control equation are discretized by time difference and space. Numerical solutions and data updates are performed sequentially according to the phase field, concentration field, and temperature field. The multiphysics solution outputs the dynamically evolving ice crystal phase field distribution matrix, concentration field distribution matrix, and temperature field distribution matrix.

[0027] Preferably, the process also includes the following:

[0028] The kinetic parameters of the dendrite tips are extracted by extracting the phase field distribution matrix, concentration field distribution matrix, and temperature field distribution matrix, and real-time kinetic characteristic data of the dendrite tips, including the radius of curvature and growth rate of the dendrite tips, are obtained. Combined with the solute diffusion coefficient in liquid seawater and the real-time kinetic characteristic data, the simulated predicted Peckle number is calculated and generated. Continuous image data of a real seawater directional freezing experimental system are acquired, and the actual dendrite tip growth rate and actual radius of curvature are extracted using image processing algorithms to calculate and generate the actual experimental Peckle number. The simulated predicted Peckle number is compared with the experimental Peckle number to calculate and obtain relative error data. The system outputs a cross-scale simulation model after error correction, and uses this model to directly output the predicted results of seawater dendrite morphology evolution and solute redistribution, realizing information mapping from the microscopic to the mesoscopic.

[0029] Compared with the prior art, the present invention has the following beneficial effects:

[0030] This invention effectively solves the morphology distortion problem of traditional phase-field models in simulating high-salinity crystallization. By introducing a microscale asymmetric driving force into the free energy functional, the phase-field evolution equation can more accurately guide the fine development of the solid-liquid interface. In simulating the directional competitive growth of ice crystals, the modified columnar dendritic morphology of the ice crystals becomes more slender, with its diameter distribution significantly shrinking from 16-45 grid steps before correction to 22-34 grid steps, and the average diameter shrinking to 89% of the original. This slender dendritic channel promotes the faster expulsion of salt ions from the ice phase, effectively reducing the retention of salt cells inside the ice crystal. Quantitative simulation data show that the expelled salt significantly increases the liquid phase salt concentration in the surrounding seawater channels from 0.068-0.079 kg / kg to 0.071-0.083 kg / kg, with an average increase in liquid phase concentration of approximately 14.5%. At the same time, the narrowing of the seawater channels provides more space for crystal expansion, resulting in an overall increase in the crystallization desalination rate of the system by approximately 3.67%.

[0031] This invention's cross-scale simulation method significantly improves the theoretical prediction accuracy of key kinetic indicators for cryogenic crystallization. The Peclet number is a core physical quantity measuring the competition between ice crystal growth and solute diffusion. In real cryogenic desalination directional crystallization experiments, the Peclet number during the dendrite differentiation stage stabilizes at around 0.28. The average Peclet number calculated by the traditional symmetric phase-field model is only 0.098, with a relative error as high as 65%, exhibiting severe interfacial diffusion distortion. After introducing an asymmetric correction term, the average Peclet number output by the simulation increases to 0.21, and the relative error with the experimentally measured value is significantly reduced to 25%. Furthermore, the high porosity resulting from slender dendrites enhances the latent heat conduction efficiency, enabling an accelerated temperature decrease of 0.2–0.4 K in local regions of the simulation system. These precise quantitative improvements demonstrate that this method enhances the physical reliability of mesoscopic simulations. Researchers can accurately predict interfacial evolution behavior at different temperatures and salinities without relying entirely on costly physical experiments, thus providing a highly reliable digital design tool for supercooling adjustment and refrigeration rate optimization in actual seawater cryogenic desalination processes. Attached Figure Description

[0032] To more clearly illustrate the technical solutions of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the following description is only one embodiment of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0033] Figure 1 This is a schematic diagram of the overall process of the cross-scale simulation method of the present invention.

[0034] Figure 2 This is a configuration diagram of the microscopic sea ice-seawater two-phase model in the embodiment.

[0035] Figure 3 The following diagrams are provided for the solid-liquid interface location and water molecule number density distribution in the embodiments: (a) is the fitting diagram of the spatial distribution of the average order parameter along the x-axis and the interface location; (b) is the spatial distribution of water molecule number density and the characteristic calibration diagram of the interface transition region.

[0036] Figure 4 The figures shown are fitting curves of the concentration at the interface as a function of sequential parameters in the embodiments, where (a) is a number density distribution of water molecules and salt ions along spatial coordinates; and (b) is an exponential fitting curve of the salt ion concentration as a function of sequential parameters.

[0037] Figure 5 The diagram shows a comparison of the phase field morphology evolution of the traditional model and the modified model at different time steps in the embodiment.

[0038] Figure 6This is a comparison of the salt concentration field distribution of the traditional model and the modified model at different time steps in the embodiment.

[0039] Figure 7 This is a comparison of the temperature field distribution of the traditional model and the modified model at different time steps in the embodiment.

[0040] Figure 8 This is a comparison diagram of the multiphysics cross-sectional spatial distribution of the traditional model and the modified model in the embodiments of the present invention;

[0041] (a) and (b) are comparison diagrams of the spatial distribution of phase field sequence parameters at the Y=100 and Y=200 sections, respectively.

[0042] (c) and (d) are comparison diagrams of the spatial distribution of solute concentration at the Y=100 and Y=200 sections, respectively;

[0043] (e) and (f) are comparison diagrams of the spatial distribution of the temperature field at the Y=100 and Y=200 sections, respectively.

[0044] Figure 9 The diagram shows the seawater directional crystallization physical experiment system and results in the embodiments of the present invention; (a) is a physical diagram of the seawater directional freezing physical experiment system; (b) is a diagram of the seawater directional crystallization results under continuous time series.

[0045] Figure 10 The figures show a comparison of the Peckle number output from physical experiments and cross-scale simulations in the embodiments of the present invention over time; (a) is a time evolution diagram of the experimental Peckle number; and (b) is a comparison diagram of the time evolution of the Peckle number calculated by the traditional model and the modified model. Detailed Implementation

[0046] The overall process of this invention is as follows Figure 1 As shown:

[0047] S1, Constructing a two-phase molecular dynamics model of sea ice-seawater:

[0048] A two-phase model of sea ice and seawater is constructed by selecting water molecule force field and electrolyte force field that can characterize phase transition characteristics; relaxation calculation is performed in equilibrium state at preset melting point temperature, and water molecule coordinate data and salt ion position trajectory data in equilibrium state are output.

[0049] S2, Extracting interface microstructure features:

[0050] Based on the water molecule coordinate data and salt ion position trajectory data output by S1, the bond order parameters of each water molecule are calculated; the spatial distribution of the bond order parameter values ​​is fitted using the hyperbolic tangent function to extract the position coordinates of the solid-liquid interface, and the spatial distribution curve of water molecule number density and salt ion concentration in the interface transition region are statistically output.

[0051] S3, Constructing a quantitative asymmetric correction function:

[0052] Based on the spatial distribution curve of water molecule number density output by S2, the hyperbolic tangent function combined with nonlinear quadratic terms is used for fitting to generate an asymmetric normalized order parameter distribution equation; at the same time, based on the spatial distribution data of salt ion concentration output by S2, an exponential fitting is performed to generate an exponential asymmetric correction function for the evolution of salt ion concentration with order parameters.

[0053] S4, Establish an asymmetric modified mesoscopic phase-field model:

[0054] The exponential asymmetric correction function generated by S3 is used as a correction term and embedded into the free energy density functional of the phase field model. After partial derivative calculation, the phase field evolution control equation containing the asymmetric chemical driving force term is generated.

[0055] S5, Coupled solution of multiphysics dynamic evolution:

[0056] The grid step size and time step size of the mesoscopic simulation are obtained. The phase field evolution control equation obtained from S4, the concentration field equation derived based on Fick's second law, and the temperature field equation established based on Fourier's law are spatiotemporally discretized. A random thermal perturbation term is added to the phase field equation and numerical iterative solution is performed to output the dynamically evolving phase field distribution matrix, concentration field distribution matrix and temperature field distribution matrix.

[0057] S6 generates a multi-scale simulation model and outputs prediction results:

[0058] Based on the exponential asymmetric correction function and the phase field evolution control equation, a cross-scale simulation model of seawater freezing desalination based on interface order parameter asymmetric correction is generated. The simulation calculation of the seawater freezing desalination process is performed using the cross-scale simulation model. Based on the phase field distribution matrix and concentration field distribution matrix output by S5, dendrite morphology evolution characteristic data and crystallization desalination rate data of the seawater freezing desalination process are extracted and output.

[0059] The present invention will be further described below with reference to embodiments. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0060] Example 1:

[0061] S1. Constructing a microscopic simulation model of molecular dynamics

[0062] A force field model capable of characterizing the phase transition properties of water molecules was selected to construct a solid-phase ice crystal model. Relaxation calculations were performed on this model under an isothermal and isobaric system to obtain a stable initial ice crystal structure. A liquid phase was then constructed based on the optimized ice crystal units, and salt ions of a predetermined concentration were introduced into the liquid phase to simulate a seawater environment. The solid and liquid phases were then joined to construct a two-phase micromechanical model of sea ice and seawater, as shown below. Figure 2 As shown. In this embodiment, the simulation system is set to operate at a preset melting point temperature and standard atmospheric pressure, and three-dimensional periodic boundary conditions are applied.

[0063] Model Construction: The TIP4P / ice water molecule model and the Madrid 2019 force field were selected. A hexagonal ice cell containing 896 water molecules was constructed, and the ice substrate was replicated along the x-axis. Simultaneously, Na+ was randomly inserted into the liquid phase region. + and Cl - Ion pairs were used to achieve a seawater molar concentration of 0.6 M to construct a complete simulation system with dimensions of 271.2 Å × 33.2 Å × 31.8 Å.

[0064] Equilibrium relaxation: The Nose-Hoover method was used to control the system temperature near the melting point, and the pressure was maintained at 1 atm. Periodic boundary conditions were used to eliminate size effects, and the Velocity-Verlet algorithm was employed for integrating the equations of motion with a time step of 1 fs. The SHAKE algorithm was used to constrain the bond lengths and bond angles of water molecules to maintain molecular rigidity. Short-range van der Waals interactions were calculated using the Lennard-Jones potential function with a cutoff radius of 10 Å; long-range electrostatic interactions were calculated using the PPPM method with a cutoff radius of 8.5 Å. During relaxation, the convergence of the total energy and density of the system was monitored until the solid and liquid phases reached a dynamic coexistence equilibrium, obtaining a stable configuration for subsequent analysis, and simultaneously outputting atomic motion trajectory files.

[0065] S2. Extracting the microstructural features of the sea ice-seawater interface.

[0066] The main purpose of this step is to identify solid-liquid phase molecules based on the simulated trajectory of molecular dynamics equilibrium by introducing bond order parameters, accurately locate the interface transition region, and provide a microscopic physical basis for subsequent mesoscopic-scale asymmetric corrections.

[0067] Bond order parameters are used as an indicator to measure the degree of local ordering in the sea ice-seawater interface. The data analysis process involves comparing probability density distribution curves, optimizing the statistical cutoff radius, and determining the critical threshold for order parameters to distinguish phases. The identification algorithm uses this threshold to differentiate between solid sea ice and liquid seawater. Subsequently, a statistical algorithm divides the simulated system into multiple small statistical regions and uses a smooth transition function to fit the spatial distribution of the average bond order parameter. This fitting step precisely locates the solid-liquid interface. Simultaneously, the statistical program extracts the spatial distribution of water molecule and salt ion number densities within the interface transition region and determines the width of the interface transition region. The specific process includes:

[0068] Order parameter calculation and cutoff radius optimization: Extract the equilibrium water molecule trajectory data output by S1; set the cutoff radii to 3.0Å, 3.2Å, 3.5Å and 3.7Å respectively, and calculate the bond order parameter values ​​of the water molecule trajectory at the corresponding cutoff radius, generating probability density distribution curves under different cutoff radii; compare multiple sets of probability density distribution curves, and select 3.5Å as the optimal statistical cutoff radius; extract the intersection point of the solid-liquid phase distribution curves under the optimal statistical cutoff radius, and set 0.36 as the phase distinction critical threshold; based on the critical threshold, traverse and determine the bond order parameter values ​​of water molecules: water molecules with bond order parameters greater than 0.36 are marked as solid phase ice molecules, and water molecules with bond order parameters less than 0.36 are marked as liquid phase water molecules, and output solid-liquid phase molecule labeling data.

[0069] Solid-liquid interface localization and number density statistics:

[0070] Extract the output solid-liquid phase molecular labeling data and water molecule coordinate data; divide the sea ice-seawater two-phase model into 271×10 small statistical regions on the xz plane, count the number of water molecules in each of these small statistical regions, and calculate and generate the spatial distribution data of the average bond order parameter along the x-axis, such as... Figure 3 As shown in (a), the spatial distribution data of the average bond order parameter is fitted using the hyperbolic tangent function to extract and generate the solid-liquid interface position coordinate data. The fitting function is as follows:

[0071] ;

[0072] in To calculate the generated average bond order parameter distribution data along the x-axis, and is the average bond order parameter between hexagonal ice and liquid water, and w is the width of the interface region. SL and I LS The position of the solid-liquid interface is represented by x, which is the extracted spatial coordinate data along the x-axis, and tanh is the hyperbolic tangent function.

[0073] Based on the solid-liquid interface coordinate data, the spatial distribution matrix of water molecule and ion number densities within the sea ice-seawater two-phase model was calculated and generated; the intermediate liquid water number density baseline value in the liquid water region was extracted to be 32 nm. 3 The number density gradient transition region between the solid and liquid phases was extracted, and the width of the interface transition region was calibrated to 15 small statistical regions with an absolute physical size of 16.3 Å. Figure 3 As shown in (b); the final output is the spatial distribution data of water molecule number density and salt ion concentration in the interface transition region, which are used as input parameters for constructing the correction function in S3.

[0074] S3. Constructing an asymmetric correction function for order parameters and concentration.

[0075] The main purpose of this step is to transform the spatial distribution law of physical quantities obtained at the microscale into a continuous mathematical expression, and to construct a quantitative correction function containing asymmetric features, so as to realize the parameter mapping and cross-scale feedback of microscopic information to mesoscopic models.

[0076] At the microscale, the distribution of order parameters at the interface was statistically analyzed, revealing a significant asymmetric distribution characteristic: a sharp decrease in order parameters on the solid side and a gradual change on the liquid side. Time averaging was performed on evolution trajectory files from multiple time points. The fitting program utilized a hyperbolic tangent function combined with a nonlinear quadratic term to construct a normalized equation for the asymmetric spatial distribution of the order parameters. The spatial distribution data of water molecule number density and salt ion concentration within the transition region of the statistical interface were also analyzed. Figure 4 As shown in (a), the relationship between the order parameter and salt ion concentration in the fitted interface transition region is as follows. Figure 4 As shown in (b), an exponential asymmetric correction correlation equation is established between the order parameter and salt ion concentration in the interface transition region. This correlation equation serves as the core correction basis for subsequent mesoscopic seawater freezing simulations. The specific process includes:

[0077] Normalization of the asymmetric distribution of the order parameter and construction of the fitting equation:

[0078] Extract the average bond order parameter spatial distribution data output by S2; normalize the average bond order parameter spatial distribution data to generate normalized order parameter distribution data; extract the normalized order parameter distribution data at multiple time points, calculate and generate time-averaged bond order parameter distribution data; fit the time-averaged bond order parameter distribution data using the hyperbolic tangent function combined with a nonlinear quadratic term to generate and output the normalized order parameter asymmetric spatial distribution equation as follows:

[0079] ;

[0080] In the formula, The normalized order parameter asymmetric distribution data generated for calculation is given by x, which is the input spatial coordinate data along the x-axis, tanh is the hyperbolic tangent mathematical function, and the constant term coefficient is 0.00065. The asymmetric feature parameters extracted for calculation are given by the existence of the quadratic term coefficient, which physically proves that the interface exhibits an inherent asymmetry that varies with the crystallization driving force on both sides.

[0081] Extract the spatial distribution data of salt ion concentration from the output of S2; extract the Na+ ion concentration from it. + and Cl - Concentration data are used to calculate the average of the two values, generating overall local concentration distribution data for the water-salt binary solution. Concentration distribution data within the interface transition region is extracted, and this data is directly mapped to the normalized order parameter asymmetric distribution data to generate a mapping dataset showing the salt ion concentration as a function of the order parameter. Based on this mapping dataset, exponential fitting is performed to generate and output the following exponential asymmetric correction equation:

[0082] ;

[0083] In the formula, To calculate the asymmetric correction data for the generated salt ion concentration, The input is the normalized order parameter asymmetric distribution data, and exp is the natural exponential mathematical function. Finally, the exponential asymmetric correction equation is output as the input data for the asymmetric correction term of the free energy density functional of the embedded mesoscopic phase field model in S4.

[0084] S4. Establish an asymmetric modified mesoscopic phase-field model.

[0085] The main purpose of this step is to break the symmetry assumption of traditional phase-field models. This invention embeds asymmetric correction terms extracted at the microscale into the mesoscale phase-field governing equations to construct a cross-scale phase-field dynamics model, accurately simulating interface evolution and solute repulsion during seawater freezing and desalination. The specific process includes:

[0086] Construction of the asymmetric modified free energy density function:

[0087] Extracting the exponential asymmetric correction equation from the S3 output The exponential asymmetric correction equation is explicitly embedded as a microscopic asymmetric correction term into the chemical free energy density functional of the basic KKS phase-field model. By fusing solid-phase free energy data, liquid-phase free energy data, and the double-well potential function, a free energy density functional equation containing the asymmetric correction term is calculated and generated, as follows:

[0088] ;

[0089] in To calculate the generated corrected chemical free energy density distribution data, c is the input concentration distribution matrix, T is the input temperature field distribution matrix, and f s (c s ) and f L (c L The data are the calculated and extracted solid and liquid phase free energy densities, respectively. For the volume fraction of the solid-liquid phase, is a monotonic interpolation function. Let w be the set double-well potential function, w be the set barrier height parameter, and k be the set coefficient of the asymmetric term. The exponential asymmetric correction data input for step S3;

[0090] Derivation and generation of the governing equations for phase field evolution:

[0091] Based on the generated free energy density functional equation containing asymmetric correction terms, and setting the dilute solution approximation condition and the constraint condition of equal chemical potential in both solid and liquid phases, partial derivatives of the free energy density functional equation are calculated to extract and generate asymmetric chemical driving force partial derivative terms. Set and introduce anisotropic parameters between the solid and liquid interfaces. The energy gradient coefficients are corrected; a random thermal perturbation matrix related to the interface position is superimposed on the equation iteration terms to generate and output the phase field evolution control equations with anisotropic and asymmetric corrections, as follows:

[0092] ;

[0093] In the formula, To calculate the evolution rate matrix of the output phase field order parameters over time, M is the phase field mobility, R is the gas constant, T is the temperature, and V is the velocity. m For molar volume, This refers to the partial derivative of the asymmetric chemical driving force in the microscopic feedback. Furthermore, this invention superimposes a random thermal perturbation term related to the interface position into the phase-field control equation to induce the growth of higher-order lateral dendrite arms at the dendrite growth interface, making the morphological evolution in the mesoscopic numerical simulation more consistent with the physical process of real seawater crystallization.

[0094] S5. Coupled solution of multiphysics dynamic evolution

[0095] The main purpose of this step is to establish continuous governing equations describing solute diffusion and latent heat transport, and to solve the phase field, concentration field, and temperature field through spatiotemporal discretization and coupling to realistically reproduce the dendrite competitive growth morphology during seawater cryogenic desalination. The specific process includes:

[0096] Equations for the evolution of solute concentration and temperature fields are constructed: By introducing Fick's second law and the equal chemical potential condition of the KKS model, the governing equations for the concentration field are derived as follows:

[0097] ;

[0098] Meanwhile, based on Fourier's law and considering the latent heat released during the seawater crystallization phase transition, the governing equations for the temperature field of the binary solution are established as follows:

[0099] ;

[0100] In the formula D is the solute diffusion coefficient. T Where is the thermal diffusivity, and L is the latent heat released upon solidification. C is the derivative of the solid-liquid phase volume fraction. P The specific heat capacity of the solution. Integrating the concentration field evolution control equation, the temperature field evolution control equation, and the phase field evolution control equation output in step S4, a multiphysics coupled dynamic equation set is generated;

[0101] Multiphysics Coupled Discrete Solution and Boundary Condition Setting: The mesoscopic computational domain is divided into a two-dimensional grid with a grid step size of 0.1 μm, and the boundary conditions are set according to the numerical stability conditions. t≤ x 2 / (5D T Calculate and set a safe time step. The forward Euler algorithm was used to discretize the phase field and concentration field equations using time-difference and spatial discretization, and the alternating explicit-implicit difference scheme was used to solve the temperature field equation. For simulation boundary treatment, Zero-Neumann boundary conditions were applied to the phase and concentration fields; Neumann boundary conditions were applied to the temperature field to accurately simulate the dynamics of constant heat flux transport. In the iterative calculation at each time step, the computational domain was divided into interface and non-interface regions, and numerical solutions were updated sequentially according to the phase field, solute field, and temperature field. Furthermore, at the interface... By introducing a stochastic thermal perturbation term with decreasing intensity, anisotropic competitive growth of ice crystals and differentiation of higher-order dendrites are induced in the supercooled environment. The final output is a dynamically evolving phase field result, as shown below. Figure 5 As shown, the concentration distribution is as follows: Figure 6 As shown, the temperature gradient distribution results are as follows: Figure 7 As shown. Simultaneously, a comparative analysis of the lateral results of the phase field, concentration field, and temperature field at different ice crystal heights before and after correction is performed, as shown below. Figure 8 As shown.

[0102] S6. Generate a cross-scale simulation model and output the prediction results.

[0103] The main purpose of this step is to quantitatively compare the numerical simulation results of the asymmetric modified phase-field model with real directional freezing experimental data by introducing dimensionless physical quantities, thereby verifying the accuracy of this invention in predicting dendrite morphology evolution and solute redistribution. The specific process includes:

[0104] Extraction and calculation of key dynamic indicators:

[0105] Extract the ice crystal phase field matrix output by S5; extract the radius of curvature and growth rate data of dendrite tips from the data matrix in real time; extract the diffusion coefficient parameter of the solute in the liquid seawater; substitute the above extracted data into the Peckley number calculation formula:

[0106] ;

[0107] In the formula, R is the radius of curvature of the dendrite tip, V is the growth rate of the dendrite tip, and D... l Let be the diffusion coefficient of the solute in liquid seawater. The simulation-predicted Pelet data for the time series is calculated and output iteratively over time steps.

[0108] Relative error calculation and model output:

[0109] Acquire real seawater directional crystallization image feature data generated by a physical experimental system based on a semiconductor cooling stage, a constant temperature circulating bath, and a CCD microscopic observation module; extract the actual tip growth rate data and actual radius of curvature data from the image feature data using image processing algorithms, and calculate and generate the real experimental Pelet number Pe. exp Extract the simulated prediction Pelet data and substitute it into the relative error calculation formula:

[0110] ;

[0111] In the formula, E represents the calculated relative error data between the simulation and the experiment; the characteristic coefficients of the exponential asymmetric correction equation in S3 are iteratively updated based on the relative error data E until the relative error data converges to within the preset accuracy threshold; finally, a cross-scale simulation model of seawater freezing desalination based on interface order parameter asymmetric correction is output after closed-loop error correction; numerical calculations are performed using the simulation model to directly output the predicted data of seawater crystallization dendrite morphology evolution and solute redistribution concentration data.

[0112] Example 2:

[0113] This embodiment verifies the above process through experimental simulation.

[0114] Setting initial parameters and boundary conditions for the mesoscopic simulation: Seawater is considered as a water-salt binary solution with a salinity of 30 ppt. The parameter configuration module calculates and sets the spatial and time steps of the mesoscopic grid based on numerical stability conditions. This step also configures thermophysical parameters such as latent heat of solidification, initial undercooling, thermal diffusivity, and solute diffusivity for the simulation system. Zero-Neumann boundary conditions are used in the phase and concentration field solutions, while Neumann boundary conditions are used in the temperature field solution. These explicit boundary conditions accurately simulate the dynamics of heat and mass transport during the actual crystallization process. Detailed simulation parameters are shown in Table 1.

[0115] Table 1. Seawater thermal property parameter settings

[0116]

[0117] A finite difference algorithm is employed to perform temporal and spatial discretization of the asymmetric phase field equations, concentration field diffusion equations, and temperature field. The computational domain is divided into interface and non-interface regions, and numerical solutions are iterated sequentially in the order of phase field, concentration field, and temperature field. During the solution process, the asymmetric correction coefficients are updated in real time, and a random thermal perturbation term with decreasing intensity towards both sides is added at the phase interface to induce anisotropic lateral branching and competitive growth of ice crystals in a supercooled environment. The dynamically evolving dendrite morphology, concentration distribution, and temperature gradient of the ice crystals are output.

[0118] At the same time, a physical experimental system for directional crystallization of seawater was established, such as Figure 9 As shown in Figure (a), the dendritic morphology of seawater solution during crystallization was observed and analyzed, and compared with the results of ice crystal simulation. The main component of the seawater solution is sodium chloride, so analytical grade sodium chloride reagent and deionized water were used to prepare a 3% sodium chloride solution. The main experimental equipment used included: a TF-HXD3012 low-temperature constant-temperature circulating bath (parameters shown in Table 2), an AST-JX101D metallurgical microscope, an AM200RZA industrial camera, and an XD-6048 semiconductor cooling plate (parameters shown in Table 3).

[0119] Table 2 Parameters of Low Temperature Constant Temperature Bath

[0120]

[0121] Table 3 Parameters of Semiconductor Cooling Plate

[0122]

[0123] The experiment first involved constructing the experimental setup according to the experimental design. A low-temperature constant-temperature circulating tank was connected to the inlet and outlet of the refrigerant. This circulating tank was responsible for controlling the temperature and flow rate of the refrigerant. Simultaneously, the required sodium chloride solution was prepared in advance. A metal substrate was installed at the bottom of the seawater crystallization tank, and the entire crystallization tank was placed on a heat-conducting substrate. This allowed the low-temperature refrigerant to exchange heat fully with the crystallization tank through the heat-conducting substrate. After the heat exchange stabilized, nitrogen gas was introduced from the inlet of the crystallization tank, flowing to expel water vapor from the device. In the subsequent data acquisition phase, a laser was placed on one side of the crystallization tank to provide a stable light source, while a CCD industrial camera was placed on the other side. The camera was directly connected to a computer to acquire and record the entire directional crystallization growth process in real time, such as... Figure 9 As shown in (b).

[0124] The actual tip growth rate data was calculated from the captured images, and the radius of curvature was measured using the camera's built-in scale function. The results were plotted as follows: Figure 10 As shown in (a), the statistical results are shown in Table 4.

[0125] Table 4. Numerical values ​​of tip growth rate and radius of curvature

[0126]

[0127] Depend on Figure 10 The comparison results in (b) show that the modified model effectively suppressed the fluctuation of Pelet number in the later stage of simulation, and also more closely approximates the experimental Pelet number. The comparison results of the three Pelet numbers in the middle stage of ice crystal growth are shown in Table 5.

[0128] Table 5 Comparison of Pecklet numbers between simulation and experimental methods

[0129]

[0130] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

[0131] While the specific embodiments of the present invention have been described above, they are not intended to limit the scope of protection of the present invention. Those skilled in the art should understand that various modifications or variations that can be made by those skilled in the art without creative effort based on the technical solutions of the present invention are still within the scope of protection of the present invention.

Claims

1. A cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters, characterized in that, The process includes the following: S1, Select the water molecule force field and electrolyte force field that can characterize the phase transition characteristics to construct a two-phase model of sea ice-seawater, perform relaxation calculations in the equilibrium state at the preset melting point temperature, and output the motion trajectory data of water molecules and salt ions in the equilibrium state. S2, calculate the bond order parameters of each water molecule in equilibrium, fit the spatial distribution of the bond order parameter values, and statistically output the spatial distribution curve of water molecule number density and the spatial distribution data of salt ion concentration in the interface transition region. S3, based on the output of S2, uses the hyperbolic tangent function to fit and obtain the asymmetric normalized order parameter distribution equation, and establishes the exponential correlation of the salt ion concentration gradient with the order parameter evolution, generating the exponential asymmetric correction function of the salt ion concentration with the order parameter evolution. S4. The exponential asymmetric correction function is used as a correction term and embedded into the free energy density functional of the phase field model. After partial derivative calculation, the phase field evolution control equation containing the asymmetric chemical driving force term is generated. S5, based on the exponential asymmetric correction function and the phase field evolution control equation, generates a cross-scale simulation model for seawater freezing and desalination, and performs simulation calculations of the seawater freezing and desalination process.

2. The cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters as described in claim 1, characterized in that: The specific process of simulation calculation for performing seawater freezing and desalination is as follows: Based on Fick's second law, the concentration field equation is derived, and the temperature field equation is established based on Fourier's law. The grid step size and time step size of the mesoscopic simulation are obtained. The control equations of the phase field, concentration field, and temperature field are spatiotemporally discretized, random thermal perturbation terms are added, and numerical iteration is performed to output the dynamically evolving phase field, concentration field, and temperature field distribution matrices. Based on the output distribution matrices, dendritic morphology evolution characteristic data and crystallization desalination rate data of the seawater freezing desalination process are extracted and output.

3. The cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters as described in claim 1, characterized in that: In step S1, a water molecule force field and an electrolyte force field characterizing the phase transition properties of water molecules are selected, and salt ion pairs are inserted into the liquid phase region to configure a seawater liquid phase with a preset salinity. Then, a solid ice substrate is spliced ​​with the seawater liquid phase to generate a sea ice-seawater two-phase model with a preset spatial size. Subsequently, the equilibrium relaxation process controls the system temperature at a preset melting point and maintains the pressure at a preset pressure. At the same time, periodic boundary conditions are applied to the microscopic model, and the monitoring program monitors the convergence state of the total energy and density of the system in real time. When the solid and liquid phases reach a dynamic coexistence equilibrium, the data on the motion trajectories of generated water molecules and salt ions are output.

4. The cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters as described in claim 1, characterized in that: The specific process of S2 is as follows: Multiple candidate cutoff radii are set, and the bond order parameter values ​​of water molecule coordinate data at the corresponding cutoff radii are calculated. Probability density distribution curves at different cutoff radii are generated. By comparing multiple sets of probability density distribution curves, the cutoff radius with the best discrimination is selected as the statistical cutoff radius. The intersection point of the solid-liquid phase distribution curves at the optimal statistical cutoff radius is extracted, and the critical threshold for phase differentiation is set accordingly. Water molecules with bond order parameters greater than the critical threshold are marked as solid-phase ice molecules, and water molecules with bond order parameters less than the critical threshold are marked as liquid-phase water molecules. Next, the number of water molecules in each tiny statistical region is counted, and the spatial distribution data of the average bond order parameter along the freezing direction is calculated. The spatial distribution data of the average bond order parameter is fitted using the hyperbolic tangent function to extract and generate the solid-liquid interface position coordinate data. Subsequently, based on the solid-liquid interface position coordinate data, the spatial distribution matrix of water molecule and salt ion number density in the two-phase model is calculated. The number density baseline value of the liquid phase region is extracted, and the number density gradient transition region between the solid and liquid phases is extracted, thereby calibrating the physical width of the interface transition region. Finally, the spatial distribution data of water molecule number density and salt ion concentration in the interface transition region are exported.

5. The cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters as described in claim 4, characterized in that: In S3, based on the spatial distribution curve of water molecule number density output from S2, a hyperbolic tangent function combined with a nonlinear quadratic term is used for fitting to generate an asymmetric normalized order parameter distribution equation; simultaneously, based on the spatial distribution data of salt ion concentration output from S2, an exponential fitting is performed to generate an exponential asymmetric correction function for the evolution of salt ion concentration with order parameters; specifically: The average bond sequence parameter spatial distribution data is used, and the average bond sequence parameter spatial distribution data is normalized to generate normalized sequence parameter distribution data. The normalized sequence parameter distribution data at multiple time points is extracted, and the time-averaged bond sequence parameter distribution data is calculated and generated. The time-averaged bond sequence parameter distribution data is spatially fitted using the hyperbolic tangent function combined with nonlinear quadratic terms to generate the normalized sequence parameter asymmetric spatial distribution equation. Subsequently, spatial distribution data of salt ion concentration is extracted, and the concentration data of the included anions and cations are extracted separately. The average value of the two is calculated to generate overall local concentration distribution data of the water-salt binary solution. The concentration distribution data is mapped to the generated normalized order parameter asymmetric distribution equation to generate a mapping dataset of salt ion concentration changing with the order parameter. Based on the mapping dataset, exponential fitting is performed and an exponential asymmetric correction equation is output.

6. The cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters as described in claim 4, characterized in that: The specific implementation process of S4 is as follows: The exponential asymmetric correction equation is used as a microscopic asymmetric correction term and embedded into the chemical free energy density functional of the basic phase-field model. Subsequently, the equation is derived, and the partial derivative of the free energy density functional equation is calculated during the derivation process to extract and generate asymmetric chemical driving force partial derivative data. Anisotropic parameters between the solid and liquid interfaces are set and introduced to obtain the asymmetric phase field evolution control equation. Random thermal perturbation matrix data related to the interface position are superimposed in the iterative terms of the equation.

7. The cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters as described in claim 2, characterized in that: Fick's second law and the condition of equal chemical potential are introduced to construct the concentration field evolution control equation. At the same time, based on the law of thermal conduction and combined with the latent heat parameter of crystallization phase transformation, the temperature field evolution control equation is established. The two evolution control equations contain specific physical parameter characteristics, including solute diffusion coefficient, thermal diffusion coefficient, latent heat release constant of solidification, solid-liquid phase volume fraction derivative, and solution specific heat capacity parameter. Subsequently, the two governing equations are merged with the output asymmetric phase field evolution governing equations to finally generate a set of multiphysics coupled dynamic equations. The phase field evolution control equation, concentration field evolution control equation, and temperature field evolution control equation are discretized by time difference and space, and numerical solutions and data updates are performed sequentially according to the phase field, concentration field, and temperature field. The multiphysics solution outputs the dynamically evolving ice crystal phase field distribution matrix, concentration field distribution matrix, and temperature field distribution matrix.

8. The cross-scale simulation method for seawater freezing and desalination based on asymmetric correction of interface order parameters as described in claim 1, characterized in that: It also includes the following processes: The kinetic parameters of the dendrite tip are extracted by extracting the phase field distribution matrix, concentration field distribution matrix, and temperature field distribution matrix, and real-time kinetic characteristic data of the dendrite tip, including the radius of curvature and growth rate of the dendrite tip, are obtained. Combined with the solute diffusion coefficient in liquid seawater and the real-time kinetic characteristic data, the simulated predicted Peckle number is calculated and generated. Continuous image data of the real seawater directional freezing experimental system are obtained, and the actual dendrite tip growth rate and actual radius of curvature are extracted using image processing algorithms to calculate and generate the real experimental Peckle number. The simulated predicted Pelet number is compared with the experimental Pelet number to calculate the relative error data; the system outputs a cross-scale simulation model after error correction, and uses this model to directly output the prediction results of seawater dendrite morphology evolution and solute redistribution, realizing information mapping from the microscopic to the mesoscopic.

Images