Ice crystal nucleation-growth-crushing cross-scale particle swarm efficient calculation method

Abstract

The invention discloses an efficient calculation method for an ice crystal nucleation-growth-crushing cross-scale particle swarm, and belongs to the technical field of ice slurry preparation and cold storage engineering. The method comprises the following steps: constructing a microscopic-mesoscopic-macroscopic three-layer cascaded cross-scale particle swarm calculation framework; implementing an adaptive particle splitting-merging algorithm; establishing a coarse-fine particle bidirectional mapping-correction dynamic coupling mechanism; respectively mapping smoothed particle fluid dynamic kernel function summation, dissipative particle dynamic force calculation and nucleation rate update to a tensor core; and outputting a cross-scale evolution result of the ice crystal scale from nanometer to millimeter and the time from microsecond to hour. According to the invention, full-link tracking from nanometer nucleation to millimeter agglomeration is realized in a single particle framework for the first time, and scale splitting in a traditional method is overcome; cNT, interfacial tension and nucleation rate prediction error lt are corrected through anisotropic hydrogen bond potential and shearing; and the method is obviously superior to an empirical model.

Description

Technical Field

[0001] This invention belongs to the field of ice slurry preparation and cold storage engineering technology, specifically involving an efficient and high-precision calculation method for the dynamic process of ice crystal nucleation-growth-fracture based on cross-scale particle swarms. Background Technology

[0002] Dynamic ice storage technology utilizes off-peak electricity at night to produce ice slurry, which is then released during the day for cooling, making it a key energy storage method. During the ice-making, transportation, storage, and cooling processes, ice crystals undergo complex dynamic events such as nucleation, growth, aggregation, and fragmentation. The spatiotemporal distribution of these events directly determines the system's heat transfer efficiency, flow resistance, and cold storage density. However, ice crystal sizes span six orders of magnitude from nanometers to millimeters, and timescales range from microsecond-level molecular collisions to hourly cold storage and release cycles. Traditional numerical methods struggle to simultaneously capture both microscopic nucleation mechanisms and macroscopic rheological behavior, leading to a significant trade-off between simulation accuracy and computational efficiency.

[0003] Existing technologies mainly employ three approaches: (1). Macroscopic homogeneous model: Treating ice slurry as a continuous medium, the energy equation and viscosity model are modified by adding source terms. However, it is impossible to analyze the evolution of individual ice crystals and it is difficult to predict non-Newtonian abrupt changes caused by aggregation and breakup. The error is often greater than 30%.

[0004] (2). Mesogroup equilibrium equation (PBE): The number density function is introduced to couple the flow and particle size distribution evolution. However, the nucleation and fragmentation models depend on empirical parameters, and the grid method faces the "curse of dimensionality" in solving high-dimensional PBE. The computational cost increases exponentially with the piecewise particle size.

[0005] (3). Molecular-scale simulation: Molecular dynamics or Monte Carlo simulation can accurately describe the formation of hydrogen bond networks and critical nuclei. However, the time step is limited to the femtosecond level, and the maximum simulation volume is less than micrometer cubic, which is more than 109 times smaller than the engineering scale. Therefore, it cannot be directly used for system-level prediction.

[0006] In recent years, meshless Lagrangian particle methods (SPH and DPD) have been explored for ice slurry simulation due to their inherent parallelism, adaptability, and ability to handle large deformations. However, current research remains limited to a single scale: SPH only tracks the macroscopic ice volume fraction, ignoring nucleation randomness; while DPD can characterize interfacial tension, it lacks a realistic nucleation thermodynamic closed loop; and microscopic CNT models cannot be coupled with the macroscopic flow field in real time. More importantly, when the particle size exceeds one million, kernel function summation and nucleation rate updates become bottlenecks, causing a sharp drop in CPU parallel efficiency, making it difficult to meet the needs of online digital twins. Therefore, there is an urgent need for a cross-scale, efficient, high-precision, and engineering-applicable method for calculating ice crystal dynamics. Summary of the Invention

[0007] The purpose of this invention is to overcome the shortcomings of existing technologies and provide an efficient computational method for cross-scale particle swarm optimization of ice crystal nucleation, growth, and fragmentation. This method couples smoothed particle hydrodynamics (SPH), improved dissipative particle dynamics (DPD), classical nucleation theory (CNT), and shear correction models through a three-layer cascaded particle framework of "micro-meta-macro," achieving high-fidelity, real-time numerical simulation of the full-scale evolution of ice crystals from nanoscale nucleation to millimeter-scale agglomeration and fragmentation. This provides core algorithmic support for the optimized design, operation control, and digital twin of dynamic ice storage systems.

[0008] The specific technical solution adopted in this invention is as follows: This invention provides an efficient computational method for cross-scale particle swarm optimization involving ice crystal nucleation, growth, and fragmentation, comprising the following steps: S1: Construct a cross-scale particle swarm computing framework with a three-tiered cascade of micro-meso-macro levels; After establishing the static hierarchical architecture, in order to accurately characterize the unsteady morphological evolution of ice crystal particles in complex flow fields, it is necessary to further introduce a discrete phase morphological evolution algorithm, as shown in S2.

[0009] S2, Implement the adaptive particle split-merge algorithm: Based on the aforementioned cross-scale particle swarm optimization framework, the ice crystal particle scale evolves according to a discrete phase distribution function, which includes aggregation and fragmentation terms. Simultaneously, Brownian motion and fluid turbulent kinetic energy are coupled to calculate the aggregation rates between the carrier fluid and ice crystals, and between ice crystals themselves. Since single-scale evolution models cannot provide feedback on nonlinear interactions between multiple levels, it is necessary to construct a cross-scale bidirectional information transmission and dynamic correction mechanism that can guarantee the conservation law, as shown in S3.

[0010] S3, establish a dynamic coupling mechanism of "coarse-fine particle bidirectional mapping-correction": Macroscopic smooth particle hydrodynamics dynamically solves the mass, momentum, and energy equations and splits on demand to generate mesoscopic dissipative particle dynamic clusters. These dissipative particle dynamic clusters further trigger microscopic nucleation events. In this process, nucleation, growth, and fragmentation information is back-weighted to the macroscopic smooth particle hydrodynamics particles, updating their ice crystal integral number, momentum, and energy. This ensures that the mass, momentum, and energy are conserved across scales without losing the details of the supercooled water flow field, establishing a dynamic coupling mechanism of coarse-fine particle bidirectional mapping and correction. To address the challenge of the exponential increase in computational load caused by high-precision coupling mechanisms, a high-performance heterogeneous parallel computing architecture based on tensor operations needs to be designed to meet the timeliness requirements of engineering-level simulations, as shown in S4.

[0011] S4, building a heterogeneous parallel architecture for tensor processing units: The smooth particle hydrodynamic kernel function summation, dissipative particle dynamic force calculation and nucleation rate update are mapped to tensor cores respectively, and the tensor parallel processing unit is used to realize the real-time calculation of the entire process of nucleation-growth-fragmentation at the scale of tens of millions of particles. Relying on the aforementioned high-throughput computing power to achieve real-time solution throughout the entire process, the final step is to complete the digital twin closed loop from numerical simulation to engineering application through the output and verification of multi-dimensional spatiotemporal data, as shown in S5.

[0012] S5: Outputs cross-scale evolution results of ice crystals from nanometers to millimeters and time from microseconds to hours. Through virtual-real interaction with experiments, it achieves high-fidelity prediction and provides digital twin support for dynamic ice storage systems.

[0013] Preferably, in S1, the microscopic layer of the cross-scale particle swarm computing framework embeds classical nucleation theory and shear-corrected nucleation model within dissipative particle dynamics clusters to predict the critical radius, contact angle state, and nucleation rate of non-uniform nucleation in real time, and to track critical nucleus growth and breakup events; the mesoscopic layer embeds dissipative particle dynamics clusters within smooth particle hydrodynamics particles, which simulate the hydrogen bond network, dipole orientation, aggregate chain conformation, and hydrogen-oxygen bond arrangement of water molecules through conservative forces, dissipative forces, and random forces, thereby calculating local interfacial tension, interfacial energy release rate, and dissipation intensity; the macroscopic layer discretizes the ice slurry flow field using smooth particle hydrodynamics, with each smooth particle hydrodynamic particle carrying macroscopic density, velocity, temperature, and ice crystal integral number to describe the macroscopic flow, heat transfer, and aggregation and breakup dynamics of ice slurry particles.

[0014] Preferably, the method for constructing the microlayer involves adding a shear-induced energy term to the free energy barrier in classical nucleation theory to increase the energy barrier, thereby reducing the critical radius in the shear flow and increasing the nucleation rate; the system free energy difference in the shear flow for: ; in, r The nucleation radius, It is the ice-water interface energy per unit area. It is the change in volume free energy per unit volume. The effective viscosity coefficient of the shear flow field is... The flow field shear rate, Shear modulus; The heterogeneous nucleation coefficient is given by the following formula: ; in, c Young's contact angle The cosine value; , , Let be the radius of the ice crystal particle at the current moment; The critical radius is given by the following formula: .

[0015] Preferably, in the method for constructing the mesoscopic layer, the force on the dissipative particle dynamics particle is proportional to its mass, and the momentum equation of the dissipative particle dynamics particle is: ; in, and Particles mass and velocity vectors External force; For particles For particles The forces acting upon it consist of conservative forces, dissipative forces, and random forces; Introducing anisotropic hydrogen bonding potential into the conservative force term Using a tetrahedral hydrogen bond network analogous to water, the expression is as follows: ; in, For particles mass; particle and particles The position vectors are respectively and , For particles and Relative distance: , ; Define from particles Pointing to particles The unit vector is ,Right now ; For unit mass particles and Conservative force coefficients between them, based on the equation of state of water and particle number density The relationship between water and its dimensionless compressibility factor Sure: ; ; in, p The pressure of water. T Thermodynamic temperature Boltzmann's constant; The conservative force coefficient between particles; It is a constant value, take ; Orientation function For a second-order Legendre polynomial: ; in, and For dissipative particle dynamics particles i and j The dipole vector is a unit vector used to reproduce the hydrogen bond orientation dependence of water molecules and the anisotropy of interfacial tension; weight function Based on the Lennard-Jones potential function, the following is determined: ; in, This is the critical particle distance; beyond this distance, the interparticle forces are negligible. The dipole vector is expressed by an order parameter tensor used to describe the degree of local orientational order in an anisotropic particle system. Q -The tensor is determined as follows: First, the spatial distribution of coarse-grained neighbor points in dissipative particle dynamics is analogous to the arrangement of microscopic molecular positions to calculate the centroid. : ; Wherein, weight function Using the Lucy function: ; Subsequently, construct the second-order moment matrix. : ; in, The charge distribution of water molecules; Finally, diagonalization yields the order parameter tensor. : ; in, For unit tensors, Maximum eigenvalue λ max The corresponding eigenvector is the dipole vector.

[0016] Preferably, the random force acts along the line connecting the particle's centers of mass and is proportional to the particle's mass; Acting on particles random forces on for: ; in, The intensity of the random force; Distance between particles Relevant kernel functions; These are random numbers with a mean of 0 and a variance of 1. The discrete time step; For particles Pointing to particles Unit vector; random force intensity The determination is based on the dimensional criterion: ; The proportionality constant is based on simulation data from the Standard Particle Model: , , , The calibration is as follows: .

[0017] in, For particle mass; The particle velocity; It represents the distance between two particles.

[0018] Preferably, the macroscopic layer uses a smoothed particle hydrodynamics method to discretize the ice slurry flow field. Based on the smoothed particle hydrodynamics method, a modeling factor and a diffusion term are added to the mass, momentum, and energy equations respectively. A weakly compressible scheme is used to stabilize the density and pressure fields and avoid spurious high-frequency noise. The specific calculation format is as follows: ; ; ; ; in: , , , Particles i Its density, velocity, pressure, and internal energy; For particles j Volume; , These are the reference density and the speed of sound, respectively. The size of the supporting domain radius; For velocity modeling factor, For pressure modeling factor, The heat flux modeling factor is determined as follows: ; ; ; in: , , ; , , Particles i and j The average speed of sound, average density, and average support domain radius are determined as follows: , , ; , , The density, velocity, and internal energy diffusion terms are defined as follows: ; ; ; Among them, parameters =1.0, =2.0, =0.1, =0.01, =0.1, =0.2, =0.02, =0.01, =0.01.

[0019] Preferably, S2 is as follows: S21: The ice crystal particle size is described using a discrete phase distribution function, and the variation in ice crystal particle size originates from the combined contribution of two dynamic processes: agglomeration and fragmentation caused by interparticle collisions. ; in, The scale is represented as The discrete phase distribution function; For aggregate terms, This is a fragmentation term; during evolution, it is calculated by pairwise collisions based on the number of neighboring nodes of the current central node, as follows: ; ; in, The scale is represented as and The aggregation rate of ice crystal particles; A weighting function that is related to the distance between particles; Let be the probability of particle breakage. For the broken kernel function; Particle breakage probability model Related to the size of the parent and daughter particles, the size is The particle size is broken down to the scale of The particle probability function is as follows: ; in, For ice crystal shape factor; The scale size is Discrete phase distribution function The range of particle sizes covered Adaptive updates are as follows: ; in, To minimize the calculated particle size, For the current moment, ice crystal particles i Particle size, ice crystal particles i The particle size of the neighboring point; After one time step of evolution, the particle i The particle size follows the discrete phase distribution function The particle size range covered has been updated as follows: ; In discrete phase distribution function Within the covered particle size range, the particle size distribution function is determined based on smooth spline interpolation: Cubic interp ( ); It should be noted that the particle breakage probability model With particle breakage probability Referring to the same variable. Among them, It only indicates the probability of breakage and does not involve specific calculations, while Involving specific calculations, writing This is to indicate the particle size before crushing. Larger than the particle size after crushing .

[0020] S22: Coupled Brownian motion with fluid turbulent kinetic energy, calculate the aggregation rate between the carrier fluid and ice crystals, and between ice crystals, as follows: Brown reunion rate for: ; in, The crystallization constant of the selected carrier solution depends on the mass and velocity of the crystallizing particles and the dynamic viscosity of the solution. By relating the kinetic energy of a single molecule to its thermodynamic temperature The calibration relationship extends to the kinetic energy at the molecular cluster scale: ,Will Revised to: ; Overall aggregation rate of ice crystal particles considering the effect of fluid turbulent kinetic energy for: ; in, The turbulent dissipation rate is calculated using a classical turbulence model.

[0021] Preferably, in S3, during the establishment of the dynamic coupling mechanism of coarse-fine particle bidirectional mapping-correction, the macroscopic scale smooth particle hydrodynamics computational domain inlet adopts experimental synchronous data-driven boundary conditions, the outlet adopts a non-reflective smooth particle hydrodynamics-dissipative particle dynamics coupling buffer layer, and the sidewall adopts a contact angle slip boundary corrected by the Young equation, ensuring that the experimental error with the real ice slurry circuit is <3%.

[0022] Preferably, in S4, the tensor parallel processing unit adopts a particle-background mesh dual decomposition, that is, the spatial domain is decomposed into the tensor processing unit core array, and each core maintains a local particle list; the smooth particle hydrodynamic kernel function summation adopts a tensor sparse accumulator, the dissipative particle dynamic force and heat calculation uses tensor dot product instructions, and the ice crystal particle size update adopts a batch sigmoid unit.

[0023] Preferably, in S5, the output results include: macroscopic subcooled water background flow field, mesoscopic molecular cluster scale dynamic reconstruction quantity, ice crystal nucleation, aggregation and breakup rates, ice crystal particle distribution function, cold storage rate, spatiotemporal hotspot map of aggregation and breakup events, which are used for real-time optimization of operating parameters and fault diagnosis of dynamic ice storage system.

[0024] Compared with the prior art, the present invention has the following advantages: (1). Cross-scale unification: For the first time, the entire link from nanonucleation to millimeter aggregation was tracked within a single particle framework, overcoming the scale fragmentation of traditional methods; (2). High fidelity: By correcting CNTs with anisotropic hydrogen bonding potential and shear, the prediction error of interfacial tension and nucleation rate is <3%, which is significantly better than the empirical model; (3) High efficiency: The combination of adaptive particle management and tensor processing units enables real-time computation of tens of millions of particles, which is two orders of magnitude faster than the traditional grid method; (4) Easy to deploy: Modular design, supports plug-and-play use of existing ice storage SCADA systems, and realizes online digital twin and closed-loop optimization. Attached Figure Description

[0025] Figure 1 Ice storage tanks for large-scale ice slurry energy storage systems; Figure 2 A method for tracking and calculating the entire process of ice crystal nucleation, growth, and fragmentation in a three-tiered cascade of "micro-meso-macro" levels; Figure 3 A schematic diagram of a multi-scale particle swarm optimization framework; Figure 4 This is a schematic diagram of the adaptive split-merge algorithm; Figure 5 This is a flowchart of the bidirectional mapping-correction process; Figure 6 Diagram of heterogeneous parallel architecture for tensor processing units; Figure 7 For experimental verification and digital twin interface diagram; This invention revolves around the entire lifecycle of ice crystal nucleation, growth, and fragmentation, proposing systematic innovations in both algorithms and software to form a cross-scale digital twin solution that can be practically implemented in engineering. The following detailed description, through specific implementation methods, elaborates on five aspects: theoretical model, numerical methods, parallel architecture, experimental verification, and software interface. Detailed Implementation

[0026] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below. Technical features in various embodiments of the present invention can be combined accordingly without mutual conflict.

[0027] This embodiment is described based on a large-scale ice slurry energy storage tank, such as... Figure 1 As shown in the diagram, the ice storage tank is externally wrapped with an insulation layer. Inlet and outlet pipes, respectively, are installed on both side walls to connect to the outside. Internally, it contains a stirring device for agitating the solution, a cooling coil for cooling the solution, and a temperature sensor for measuring the solution temperature. The ice storage tank's liquid level and parameters such as temperature and flow rate during operation can be monitored and adjusted via a control panel.

[0028] based on Figure 1 The apparatus shown in this embodiment provides an efficient computational method for cross-scale particle swarm optimization, such as the method for calculating ice crystal nucleation, growth, and fragmentation. Figure 2As shown, this method addresses the accuracy and efficiency bottlenecks of traditional ice crystal dynamics models in cross-scale coupling, tracking of the entire nucleation-growth-fragmentation process, and efficient computation. It innovatively constructs a three-tiered cross-scale particle swarm computing framework of "macro-meso-micro": (1) The macro layer uses smooth particle hydrodynamics (SPH) to describe the macroscopic flow and heat transfer of ice slurry, as well as the dynamic behaviors of aggregation and fragmentation between ice crystal particles; (2) The meso-scale layer introduces improved dissipative particle dynamics (DPD) to simulate the hydrogen bonds, dipole orientation, aggregate chain conformation, and hydrogen-oxygen bond arrangement of water molecules to characterize interfacial tension, interfacial energy release, and dissipation intensity; (3) The micro-scale layer embeds classical nucleation theory (CNT) and shear-corrected nucleation model to achieve real-time prediction of non-uniform nucleation, contact angle state, critical radius, and nucleation rate. The three layers are dynamically coupled through a "coarse-fine particle bidirectional mapping-correction" mechanism. Macroscopic SPH particles generate / annihilate mesoscopic DPD clusters in real time, and the DPD clusters trigger microscopic nucleation events. The nucleation, growth, and fragmentation information are then back-weighted to the macroscopic field to ensure the conservation of mass, momentum, and energy.

[0029] To further improve computational efficiency, this invention proposes an "adaptive particle splitting-merging algorithm": when the local ice crystal number density gradient or shear rate exceeds a threshold, particle splitting is automatically triggered to resolve the fine structure; when the gradient is below the threshold, particle aggregation is implemented to reduce the degrees of freedom. Simultaneously, a heterogeneous parallel architecture of the TPU tensor processing unit is designed, mapping SPH kernel function summation, DPD force and thermal calculations, and nucleation rate updates to the TPU tensor core, enabling real-time computation of the entire nucleation-growth-fragmentation process at a particle scale of tens of millions. While maintaining an error of less than 3% compared to experimental data, this invention improves computational speed by two orders of magnitude compared to traditional grid methods, enabling simultaneous simulation of cross-scale evolution of ice crystals from nanometers to millimeters and time periods from microseconds to hours. This provides a high-fidelity digital twin tool for the optimized design and operational control of dynamic ice storage systems.

[0030] The efficient computational method for cross-scale particle swarm optimization of ice crystal nucleation, growth, and fragmentation described above in this invention specifically includes the following steps: S1: Construct a three-tiered cascaded cross-scale particle swarm computing framework of "micro-meso-macro" level, such as... Figure 3 As shown.

[0031] The micro-layer incorporates classical nucleation theory (CNT) and a shear-modified nucleation model within the DPD (Dissipative Particle Dynamics) cluster. This allows for real-time prediction of the critical radius, contact angle state, and nucleation rate in non-uniform nucleation, while also tracking critical nucleus growth and fragmentation events. The shear-modified nucleation model adds a shear-induced energy term to the CNT free energy barrier, thereby increasing the energy barrier, which in turn reduces the critical radius in the shear flow and improves the nucleation rate. The simplified steps are as follows: Step 1-1: Construct the system free energy difference in shear flow: ; in, r The nucleation radius, It is the ice-water interface energy per unit area. It is the change in volume free energy per unit volume. The effective viscosity coefficient of the shear flow field is... The flow field shear rate, Shear modulus; The heterogeneous nucleation coefficient depends on the contact angle, the radius of the nucleating particle, and the curvature of the nucleation site, as shown in the following formula: ; in, c Young's contact angle The cosine value; , , Let be the radius of the ice crystal particle at the current moment.

[0032] Steps 1-2: Determine the critical radius and critical free energy difference The formula is as follows: ; ; Steps 1-3 determine the number of ice nuclei per unit volume, using the following formula: ; It is the number density of water molecules in the liquid phase. T Thermodynamic temperature Boltzmann's constant; Steps 1-4 determine the possible rate at which water molecules can enter the ice nucleus, using the following formula: ; It is the number of molecules within the jump distance range on the nuclear surface. It is the activation energy for the diffusion of water molecules across the nucleus boundary. h It is Planck's constant.

[0033] Steps 1-5 determine the Zeldovich factor, using the following formula: ; in, The volume of liquid water: , The volume of a single water molecule It is the critical Gibbs free energy barrier.

[0034] Steps 1-6 determine the steady-state nucleation rate of supercooled water in the heterogeneous nucleation process: .

[0035] The mesoscopic layer consists of modified dissipative particle dynamics (DPD) clusters embedded within SPH particles. DPD particles simulate the hydrogen bond network, dipole orientation, polymer chain conformation, and hydroxyl bond arrangement of water molecules through conservative, dissipative, and random forces. This allows for the calculation of local interfacial tension, interfacial energy release rate, and dissipation intensity, including the following steps: Step 2-1: Construct the momentum equation for the DPD particle, where the force acting on the particle is proportional to its mass: ; in, and Each of them is a particle mass and velocity vectors For particles For particles The force, It is an external force. It consists of conservative forces, dissipative forces, and random forces, and specifically includes the following steps: Step 2-2, in the conservative force term Introducing anisotropic hydrogen bonding potentials, analogous to the tetrahedral hydrogen bond network of water: ; in, For particles The quality; For particles and Relative distance: , , .

[0036] Steps 2-3, based on the equation of state of water and particle number density The relationship between water and its dimensionless compressibility factor Determine the unit mass of particles and The gravitational force generated between : ; ; p The pressure of water. T Thermodynamic temperature Boltzmann's constant; It is a constant value. .

[0037] Steps 2-4: Determine the orientation function For a second-order Legendre polynomial: ; Among them, among them, and For dissipative particle dynamics particles i and j The dipole vector is a unit vector used to reproduce the hydrogen bond orientation dependence of water molecules and the anisotropy of interfacial tension.

[0038] Steps 2-5: Determine the weighting function in the conservative force based on the Lennard-Jones potential function. as follows: ; in, This is the critical particle distance; beyond this distance, the interparticle force is negligible.

[0039] Steps 2-6 employ the order parameter tensor, which describes the degree of local orientation order in anisotropic particle systems. Q - The tensor determines the dipole vector, specifically including the following steps: (1). The centroid is calculated by analogy between the spatial distribution of coarse-grained neighbor points in the DPD and the positional arrangement of microscopic molecules: ; Wherein, weight function Using the Lucy function: ; (2). Construct the second-order moment matrix: ; in, This represents the charge distribution of water molecules. According to the classical water potential energy model TIP4P / 2005, O = -1.1128e, H = +0.5564e, which allows us to determine... = [0.5564, 0.5564, -1.1128].

[0040] (3). Diagonalization yields Q -tensor: ; Maximum eigenvalue λ max The corresponding eigenvector is the dipole vector.

[0041] Steps 2-7: Construct particles along the line connecting the particle's center of mass. Effect on The random force on is: ; in, The intensity of the random force; Distance between particles Relevant kernel functions; These are random numbers with a mean of 0 and a variance of 1.

[0042] Steps 2-8: Determine the intensity of random force based on dimensional criteria. : ; Steps 2-9, based on simulation data from the standard particle model: , , , Calibrate random force strength The proportionality constants are as follows: ; Steps 2-10: Construct particles along the line connecting the particle's center of mass. Effect on The dissipation force on it is: ; in: It is the intensity of dissipation force; Distance between particles The relevant weighting function; each pair of particles is independent for each time step. Relative velocity: .

[0043] Steps 2-11: Based on the dissipative fluctuation theorem and the law of conservation of energy, construct the equilibrium conditions between dissipative forces and stochastic forces: ; Steps 2-12 determine the functional relationship between the weighting function and the inter-particle distance in dissipative and random forces: ; Steps 2-13, the heat transfer equation for DPD particles is constructed as follows: ; in, and Particles Specific heat capacity and temperature; , and These represent the interparticle mesoscopic thermal conductivity, viscous thermal flux, and random thermal flux, respectively, induced by particles. The specific steps include: Steps 2-14: Determine the mesoscopic heat transfer flux: ; Steps 2-15: Adjust the intensity of the contact heat flow. as follows: ; in, Specific heat capacity, It is the mesoscopic thermal conductivity.

[0044] Steps 2-16: Determine the viscous heat flow: ; Steps 2-17: Determine random heat flow: ; Steps 2-18: Determine the random heat flux intensity based on the dissipative fluctuation theorem. and contact heat flow intensity Relationship: ; Steps 2-19: Construct random numbers that conform to a Gaussian distribution. and as follows: ; ; Ensure that the random forces and random heat flows between different particle pairs at different times are independent of each other, and guarantee the conservation of momentum and energy of the particle system.

[0045] Steps 2-20 modify the Velocity-Verlet algorithm by using the current position, velocity, and force of the DPD particle to calculate the position and velocity at the next moment, and then using the new position and velocity to calculate the new force. The iterative format is as follows: ; ; ; ; ; ; ; The macroscopic layer employs the Smooth Particle Hydrodynamics (SPH) method to discretize the ice slurry flow field. Each SPH particle carries macroscopic density, velocity, temperature, and ice crystal integral number to describe the macroscopic flow, heat transfer, and aggregation and breakup dynamics of the ice slurry particles. Based on the macroscopic SPH particles, modeling factors and diffusion terms are added to the mass, momentum, and energy equations, respectively. A weakly compressible scheme is used to stabilize the density / pressure field and avoid spurious high-frequency noise. The main steps are as follows: Step 3-1, construct the SPH numerical calculation format: ; ; ; ; in: , , , Particles i Its density, velocity, pressure, and internal energy; For particles j The volume. , For reference density and speed of sound, The radius of the supporting domain.

[0046] Step 3-2: Determine the velocity modeling factor. Pressure modeling factor Heat flux modeling factor ,as follows: Velocity scaling factor: ; Pressure modeling factor: ; Heat flux modeling factor: ; in: , , .

[0047] Step 3-3: Identify the particles i and j average speed of sound Average density and average support domain radius for: , , .

[0048] Steps 3-4: Determine the density, velocity, and internal energy diffusion terms. , , as follows: ; ; ; Steps 3-5: Determine the series of parameters in the macroscopic SPH calculation format. , , , , , , , , as follows:

[0049] S2: Implement the adaptive particle split-merge algorithm; The ice crystal particle size evolves according to a discrete phase distribution function, which includes aggregation and fragmentation terms. The aggregation rate between the carrier fluid and ice crystals, and between ice crystals themselves, is influenced by a combination of Brownian aggregation and fluid turbulent kinetic energy. The particle fragmentation function is decomposed into two parts: a fragmentation kernel model characterizing the particle fragmentation rate coefficient, and a particle fragmentation probability model. For example... Figure 4 The specific steps are as follows: Step 4-1: Construct a distribution function to describe the ice crystal particle size, and the variation in ice crystal size originates from the combined contribution of two dynamic processes: agglomeration and fragmentation caused by inter-particle collisions. ; in, The scale is represented as The discrete phase distribution function; For aggregate terms, This is a broken term.

[0050] Step 4-2: Calculate the aggregation term by pairwise collisions based on the number of neighboring nodes of the central node. and broken items ,as follows: ; ; in, The scale is represented as and The aggregation rate of ice crystal particles; A weighting function that is related to the distance between particles; Let be the probability of particle breakage. This is the broken kernel function.

[0051] Step 4-3: Construct the aggregation rates between the carrier fluid and ice crystals, and between ice crystals themselves, under the combined effects of Brownian motion and fluid turbulent kinetic energy. The Brownian aggregation rate is: ; in, The crystallization constant of the selected carrier solution depends on the mass and velocity of the crystallizing particles and the dynamic viscosity of the solution. .

[0052] Step 4-4: Relationship between the kinetic energy of a single molecule and the thermodynamic temperature The calibration relationship extends to the kinetic energy at the molecular cluster scale: , correct for: ; Steps 4-5 further consider the effect of fluid turbulence kinetic energy to determine the overall aggregation rate of ice crystal particles as follows: ; in, The turbulent dissipation rate is calculated using a classical turbulence model.

[0053] Steps 4-6 decompose the particle breakage function into two parts: the breakage core model of the particle breakage rate coefficient. Particle breakage probability model .

[0054] Steps 4-7 involve associating the particle breakage probability model with the scales of the parent and daughter particles to determine the scale. The particle size is broken down to the scale of The particle probability function is as follows: ; in, For ice crystal shape factor.

[0055] Steps 4-8: Remove ice crystal particles i Particle size distribution function The range of particle sizes covered Adaptive updates are as follows: ; in, To minimize the calculated particle size, For the current moment, ice crystal particles i Particle size, ice crystal particles i The particle size of the neighboring point.

[0056] Steps 4-9, after one time step of evolution, the particle... iThe particle size is updated according to the particle size range covered by the particle size distribution function as follows: ; Steps 4-10: Within the particle size range covered by the particle size distribution function, the particle size distribution function is determined based on smooth spline interpolation. Cubic interp ( ); S3: Establish a dynamic coupling mechanism of "coarse-fine particle bidirectional mapping-correction", the specific steps of which are as follows: Figure 5 ,include: Step 5-1: Macroscopic SPH particles dynamically solve the mass, momentum, and energy equations and split as needed to generate mesoscopic DPD clusters. The macroscopic field distribution provides background parameters for the dynamic evolution of mesoscopic DPD clusters, and the DPD clusters further trigger microscopic nucleation events. Step 5-2: Nucleation, growth, and fragmentation information are back-weighted to macroscopic SPH particles to update their ice crystal integral number, momentum, and energy, ensuring cross-scale mass, momentum, and energy conservation without losing details of the supercooled water flow field. Step 5-3: During the "coarse-fine particle" bidirectional coupling process, the macroscopic-scale SPH computational domain inlet uses experimentally synchronized data to drive boundary conditions, the outlet uses a non-reflective SPH-DPD coupling buffer layer, and the sidewall uses a contact angle slip boundary corrected by the Young equation to ensure that the experimental error with the real ice slurry circuit is <3%.

[0057] S4: Building a heterogeneous parallel architecture for tensor processing units, the main steps are as follows: Figure 6 ,include: Step 6-1: The SPH kernel function summation, DPD force calculation and nucleation rate update are mapped to tensor cores respectively. The tensor parallel processing unit is used to realize the real-time calculation of the entire process of nucleation-growth-fragmentation at the scale of tens of millions of particles. Step 6-2: The parallel architecture of the tensor processing unit adopts "particle-background mesh dual decomposition": the spatial domain is decomposed into the core array of the tensor processing unit, and each core maintains a local particle list; the SPH kernel function summation adopts the tensor sparse accumulator, the DPD force and heat calculation uses the tensor dot product instruction, and the ice crystal particle size update adopts the batch Sigmoid unit.

[0058] S5: Outputs cross-scale evolution results of ice crystals from nanometers to millimeters and time from microseconds to hours. Output results include: macroscopic-scale supercooled water background flow field, mesoscopic molecular cluster-scale dynamic reconstruction, ice crystal nucleation, aggregation, and fragmentation rates, ice crystal particle distribution function, cold storage rate, and spatiotemporal hotspot maps of aggregation / fragmentation events. These results are used for real-time optimization of operating parameters and fault diagnosis of dynamic ice storage systems. Through virtual-real interaction with experiments, high-fidelity predictions are achieved, providing digital twin support for dynamic ice storage systems, such as... Figure 7 .

[0059] S6: All algorithms in steps S1 to S5 are encapsulated in Python / Matlab hybrid code and directly embedded into the existing ice storage digital twin platform through seamless switching between dynamic and static graphs.

[0060] The core contribution of this invention in the "model optimization" section is: replacing the "static supercooled water" assumption in the traditional nucleation theory (CNT) with the "shear-supercooled coupling" assumption, and using eight closed-loop formulas from steps 1-1 to 1-6 above to incorporate engineering-measurable parameters such as shear rate, contact angle, curvature, and viscous dissipation into the nucleation rate formula in one go, enabling the model to operate within 0-2000 s. -1 Within the shear rate range, the experimental error is < 3%. Specifically, the 6-step optimization steps of the shear-corrected nucleation model are as follows: Step 1: Free Energy Reconfiguration Adding a shear-induced energy term to the free energy of classical nucleation theory yields a new system free energy, as shown in the following equation: ; Step 2: Further correction of heterogeneity coefficients: The original homogeneity coefficient only considered the static contact angle. By introducing the nucleation curvature and the dynamic Young's contact angle, the final heterogeneity coefficient is obtained, as shown in the following formula: ; Step 3 Analytical solution for critical radius: right Differentiate and set the derivative to zero to give an explicit solution: ; This formula changes the critical radius from the original... The determination becomes simultaneously determined by the effective viscosity coefficient, shear rate, and shear modulus of the shear flow field, thus realizing "the greater the shear, the better." The physical trend is "the smaller".

[0061] Step 4 Nucleation Free Energy Barrier and factor As shown in the following formula: ; ; Step 5 Steady-state nucleation rate J The closed-form formula is shown below: ; Performing steps 1-5 above yields an explicit algebraic formula containing only local macroscopic quantities, which can be updated within each step without iteration.

[0062] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the invention. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention. Therefore, all technical solutions obtained through equivalent substitution or transformation fall within the protection scope of the present invention.

Claims

1. A high-efficiency computing method for ice crystal nucleation-growth-fragmentation across-scale particle population, characterized in that, Includes the following steps: S1: Construct a cross-scale particle swarm computing framework with a three-tiered cascade of micro-meso-macro levels; S2: Based on the aforementioned cross-scale particle swarm optimization framework, the ice crystal particle scale evolves according to a discrete phase distribution function, which includes an aggregation term and a fragmentation term; simultaneously, Brownian motion and fluid turbulent kinetic energy are coupled to calculate the aggregation rate between the carrier fluid and the ice crystal, and between ice crystals themselves. S3: Macroscopic smooth particle fluid dynamics particles dynamically solve the mass, momentum, and energy equations and split on demand to generate mesoscopic dissipative particle dynamic clusters. The dissipative particle dynamic clusters further trigger microscopic nucleation events. In this process, the nucleation, growth, and fragmentation information is back-weighted to the macroscopic smooth particle fluid dynamics particles, updating their ice crystal integral number, momentum, and energy. This ensures that the cross-scale mass, momentum, and energy are conserved without losing the details of the supercooled water flow field, establishing a dynamic coupling mechanism of coarse-fine particle bidirectional mapping-correction. S4: The smooth particle hydrodynamic kernel function summation, dissipative particle dynamic force calculation and nucleation rate update are mapped to tensor cores respectively, and the tensor parallel processing unit is used to realize the real-time calculation of the entire process of nucleation-growth-fragmentation at the scale of tens of millions of particles. S5: Outputs cross-scale evolution results of ice crystals from nanometers to millimeters and time from microseconds to hours. Through virtual-real interaction with experiments, it achieves high-fidelity prediction and provides digital twin support for dynamic ice storage systems.

2. The ice crystal nucleation-growth-fragmentation multiscale particle population high-efficiency computational method of claim 1, wherein, In S1, the microscopic layer of the cross-scale particle swarm computing framework embeds classical nucleation theory and shear-corrected nucleation models within dissipative particle dynamics clusters to predict the critical radius, contact angle state, and nucleation rate of non-uniform nucleation in real time, and to track critical nucleus growth and breakup events. The mesoscopic layer embeds dissipative particle dynamics clusters within smooth particle hydrodynamics particles. These particles simulate the hydrogen bond network, dipole orientation, aggregate chain conformation, and hydrogen-oxygen bond arrangement of water molecules through conservative forces, dissipative forces, and random forces, thereby calculating local interfacial tension, interfacial energy release rate, and dissipation intensity. The macroscopic layer discretizes the ice slurry flow field using smooth particle hydrodynamics. Each smooth particle hydrodynamic particle carries macroscopic density, velocity, temperature, and ice crystal integral number to describe the macroscopic flow, heat transfer, and aggregation and breakup dynamics of ice slurry particles.

3. The efficient computational method for cross-scale particle swarm optimization based on ice crystal nucleation-growth-fragmentation according to claim 2, characterized in that, The construction method of the micro layer is to increase a shear induction energy item in a free energy barrier of a classical nucleation theory to increase an energy barrier, thereby reducing a critical radius in the shear flow and increasing a nucleation rate; and a system free energy difference in the shear flow is : ； wherein, r is the nucleation radius, is the ice-water interfacial energy per unit area, is the change in volume free energy per unit volume, is the effective viscosity coefficient of the shear flow field, is the shear rate of the flow field, is the shear modulus; is the heterogeneous nucleation coefficient, which is given by the formula: ； in, c Young's contact angle The cosine value; , , Let be the radius of the ice crystal particle at the current moment; The critical radius is given.

4. The efficient computational method for cross-scale particle swarm optimization based on ice crystal nucleation-growth-fragmentation according to claim 2, characterized in that, In the method for constructing the mesoscopic layer, the force on the dissipative particle dynamics particle is proportional to its mass, and the momentum equation of the dissipative particle dynamics particle is: ； in, and Particles mass and velocity vectors External force; For particles For particles The forces acting upon it consist of conservative forces, dissipative forces, and random forces; Introducing anisotropic hydrogen bonding potential into the conservative force term Using a tetrahedral hydrogen bond network analogous to water, the expression is as follows: ； in, For particles mass; particle and particles The position vectors are respectively and , For particles and Relative distance: , ; Define from particles Pointing to particles The unit vector is ,Right now ; For unit mass particles and The gravitational force generated between them; Orientation function For a second-order Legendre polynomial: ； in, and For particles i and j The dipole vector is a unit vector used to reproduce the hydrogen bond orientation dependence of water molecules and the anisotropy of interfacial tension; weight function Based on the Lennard-Jones potential function, the following is determined: ； in, This is the critical particle distance; beyond this distance, the interparticle force is negligible.

5. The efficient computational method for cross-scale particle swarm optimization of ice crystal nucleation-growth-fragmentation according to claim 4, characterized in that, The random force acts along the line connecting the particle's centers of mass and is proportional to the particle's mass; Acting on particles random forces on for: ； in, For particles and particles Random force intensity between; Distance between particles Relevant kernel functions; These are random numbers with a mean of 0 and a variance of 1. The discrete time step; For particles Pointing to particles Unit vector; random force intensity The determination is based on the dimensional criterion: ； The proportionality constant is based on simulation data from the Standard Particle Model: , , , The calibration is as follows: ； in, For particle mass; The particle velocity; It represents the distance between two particles.

6. The efficient computational method for cross-scale particle swarm optimization based on ice crystal nucleation-growth-fragmentation according to claim 2, characterized in that, The macroscopic layer uses a smoothed particle hydrodynamics method to discretize the ice slurry flow field. Based on smoothed particle hydrodynamics, a modeling factor and a diffusion term are added to the mass, momentum, and energy equations respectively. A weakly compressible scheme is used to stabilize the density and pressure fields and avoid spurious high-frequency noise. The specific calculation format is as follows: ； ； ； ； in: , , , Particles i Its density, velocity, pressure, and internal energy; For particles j Volume; , These are the reference density and the speed of sound, respectively. The size of the supporting domain radius; For velocity modeling factor, For pressure modeling factor, The heat flux modeling factor is determined as follows: ； ； ； in: , , ; , , Particles i and j The average speed of sound, average density, and average support domain radius are determined as follows: , , ; , , The density, velocity, and internal energy diffusion terms are defined as follows: ； ； ； Among them, model parameters =1.0, =2.0, =0.1, =0.01, =0.1, =0.2, =0.02, =0.01, =0.

01.

7. The efficient computational method for cross-scale particle swarm optimization of ice crystal nucleation-growth-fragmentation according to claim 1, characterized in that, S2 is specifically as follows: S21: The ice crystal particle size is described using a discrete phase distribution function, and the variation in ice crystal particle size originates from the combined contribution of two dynamic processes: agglomeration and fragmentation caused by interparticle collisions. ； in, The scale is represented as The discrete phase distribution function; For aggregate terms, This is a fragmentation term; during evolution, it is calculated by pairwise collisions based on the number of neighboring nodes of the current central node, as follows: ； ； in, The scale is represented as and The aggregation rate of ice crystal particles; A weighting function that is related to the distance between particles; Let be the probability of particle breakage. For the broken kernel function; S22: Ice crystal aggregation rate due to the combined effects of Brownian motion and fluid turbulent kinetic energy. for: ； in, The turbulent dissipation rate is calculated using a classical turbulence model.

8. The efficient computational method for cross-scale particle swarm optimization of ice crystal nucleation-growth-fragmentation according to claim 1, characterized in that, In S3, during the establishment of the dynamic coupling mechanism of coarse-fine particle bidirectional mapping-correction, the macroscopic scale smooth particle hydrodynamics computational domain inlet adopts experimental synchronous data-driven boundary conditions, the outlet adopts a non-reflective smooth particle hydrodynamics-dissipative particle dynamics coupling buffer layer, and the sidewall adopts a contact angle slip boundary corrected by the Young equation to ensure that the experimental error with the real ice slurry circuit is <3%.

9. The efficient computational method for cross-scale particle swarm optimization based on ice crystal nucleation-growth-fragmentation according to claim 1, characterized in that, In S4, the tensor parallel processing unit adopts a particle-background mesh dual decomposition, that is, the spatial domain is decomposed into the core array of the tensor processing unit, and each core maintains a local particle list; the smooth particle hydrodynamic kernel function summation adopts a tensor sparse accumulator, the dissipative particle dynamic force and heat calculation uses tensor dot product instructions, and the ice crystal particle size update adopts a batch sigmoid unit.

10. The efficient computational method for cross-scale particle swarm optimization based on ice crystal nucleation-growth-fragmentation according to claim 1, characterized in that, In S5, the output results include: macroscopic supercooled water background flow field, mesoscopic molecular cluster scale dynamic reconstruction quantity, ice crystal nucleation, aggregation and breakup rates, ice crystal particle distribution function, cold storage rate, spatiotemporal hotspot map of aggregation and breakup events, which are used for real-time optimization of operating parameters and fault diagnosis of dynamic ice storage system.

Images