Film icing dynamics simulation method based on multilayer particle representation

Abstract

The invention discloses a film icing dynamics simulation method based on multilayer particle representation, and belongs to the technical field of crossing of computer graphics and physical simulation. The invention provides a novel phase field mapping method, ice crystal phase field updating is carried out on a layer of particles with fixed positions, the ice crystal center is used as an anchor point to establish a mapping relation with the surface of a thin film, and a structured icing pattern is kept on the moving particles. Meanwhile, according to the method, a temperature updating equation is introduced into a film dynamic model, a boundary particle film surfactant concentration calculation formula based on the moving least square method is deduced, and the stability of a simulation system is enhanced. According to the method, the defect of related research on film icing phenomenon simulation in the field of computer graphics is filled, and the method has the advantages of being high in stability under complex conditions and supporting a user to control the ice crystal shape.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of computer graphics and physical simulation, and particularly relates to a thin film icing dynamic simulation method based on a multi-layer particle representation. BACKGROUND

[0002] Thin film icing is a natural physical phenomenon, and its dynamic process involves phase transition, heat conduction, and coupling of fluid and solid, which has important scientific research and visual effect application value. Traditional thin film fluid simulation methods mainly achieve this by solving Navier-Stokes, while crystal simulation relies on the diffusion equation of temperature and phase field. In recent years, particle methods have been widely used in fluid simulation, greatly improving the ability to capture convection and vortex details in simulation, but there is still no method to model and simulate thin film icing phenomena, and the following shortcomings still need to be overcome:

[0003] 1. Difficulty in maintaining ice crystal pattern: Traditional phase field methods are difficult to maintain the structured pattern of ice crystals under dynamic particle representation, especially under the action of complex surface flow, the detailed branches of ice crystals are easy to distort due to uneven particle distribution or movement.

[0004] 2. Insufficient coupling of temperature field and phase field: Existing methods usually ignore the influence of temperature change on surface tension, and cannot accurately simulate the Marangoni flow driven by the temperature gradient caused by the release of latent heat during phase transition, making it difficult to reproduce the real ice crystal dynamic behavior.

[0005] 3. Unstable boundary condition processing: Under complex boundary conditions (such as ice crystal and fluid coupling, multi-phase interaction), the simulation system of existing methods is easy to collapse or produce non-physical phenomena, lacking a robust boundary processing mechanism.

[0006] In the field of computer graphics, there have been attempts to simulate ice crystal growth and thin film dynamics through phase field theory and particle methods. For example, Kobayashi (1993) proposed an ice crystal growth model based on phase field, and Ren et al. (2018) achieved artistic control of ice crystal shape by introducing a direction field. However, these methods are mainly for static scenes and cannot handle ice crystal movement on dynamic thin film surfaces. The MELP (Moving Eulerian-Lagrangian Particles) method performs well in thin film and foam simulation, but lacks support for temperature field and phase transition, and cannot be directly applied to icing simulation.

[0007] In addition, the existing fluid-rigid coupling method (such as the SPH boundary processing technology) performs poorly in complex phase change scenarios, especially when ice crystals interact with the thin film boundary, which is prone to numerical instability. Therefore, there is an urgent need for a thin film icing simulation method that can maintain the structured pattern of ice crystals under dynamic particle representation, accurately simulate the coupling of temperature field and phase field, and have robust boundary processing capability. SUMMARY

[0008] The purpose of the present application is to propose a thin film icing simulation method under multi-layer particle representation, which can maintain the complex structured shape of ice crystals growing under the action of complex surface flow, and solve the stability problem caused by ice crystal movement, temperature change and complex boundary conditions in the process of soap bubble icing. The present application fills the gap in the prior art by introducing a phase field mapping method and multi-layer particle representation, combining the diffusion equations of temperature and phase field and the moving least squares method boundary processing technology, and provides an efficient, stable and controllable solution for physical simulation of thin film icing.

[0009] To achieve the above purpose, the technical scheme adopted by the present application is:

[0010] A thin film icing dynamics simulation method based on multi-layer particle representation, comprising the following steps:

[0011] a. A theoretical framework for capturing potential thermodynamic mechanisms in two-dimensional and three-dimensional thin film icing processes is proposed;

[0012] The potential thermodynamic mechanisms in the two-dimensional and three-dimensional thin film icing processes are simulated using three layers of particles, namely Euler particles, Lagrangian particles and graph particles. The Euler particles and Lagrangian particles together represent the two-dimensional and three-dimensional thin film surface, and the main dynamics calculation process is performed on the Euler particles. The Lagrangian particles and graph particles together represent the ice crystals, and the icing calculation process is performed on the graph particles. In the three-layer particle, the Lagrangian particles are used to represent the material state and visualization, and are connected to the graph particles and Euler particles to transfer the current phase field properties and other physical quantities to the graph particles, and collect thermodynamic physical quantities from the graph particles. The dynamics physical quantities and thermodynamic physical quantities are then transferred to the Euler particles, and the momentum is collected from the Euler particles. The three types of particles define local coordinate systems according to their normal vectors, which are used to project neighboring particles onto their tangent planes;

[0013] b. An advanced phase field mapping method is proposed to maintain the structured pattern of ice crystals under moving particle representation, and to support user control of the pattern during simulation;

[0014] The user can directly control the pattern during simulation by pre-designed gray map editing and fixed local particle direction field control of the pattern particle, or by directly updating the phase field according to the normalized time-varying directed distance field, and calculating the update of the temperature field according to the update of the phase field;

[0015] c. A boundary particle thin film surfactant concentration calculation formula based on the moving least square method is derived under the multi-layer particle representation to improve the stability of the simulation system;

[0016] The boundary particle thin film surfactant concentration calculation formula based on the moving least square method is that, during initialization, the boundary particle sampling is performed at the position where the solid boundary is close to the Euler particle, and the particle spacing is consistent with the initial spacing of the Euler particle; before each iteration of each relaxation Jacobi iteration of the implicit incompressible smoothed particle hydrodynamics (IISPH), the sampled boundary particle estimates its surfactant concentration using the moving least square method according to the surfactant concentration of the neighboring Euler particle.

[0017] Further, in the three-layer particle, the Lagrangian particle is used to represent the material state and visualization, and connects the pattern particle and the Euler particle, and transmits the current phase field attribute and other physical quantities to the pattern particle, and collects the thermodynamic physical quantity from the pattern particle, and then transmits the dynamic physical quantity and the thermodynamic physical quantity to the Euler particle, and collects the momentum from the Euler particle.

[0018] Further, the Euler particles are distributed relatively uniformly, and in the calculation of each time step, the redistribution is performed according to the number density of the Euler particles to ensure that the number of the Euler particles per unit area is relatively consistent in the thin film; and in the rendering process, the Euler particles are used for surface reconstruction of the thin film mesh body.

[0019] Further, the temperature equation is introduced in the dynamic calculation process, and the update formula of the thin film surfactant concentration based on the implicit incompressible smoothed particle hydrodynamics (IISPH) is re-derived, and the relaxation Jacobi iteration method is used for calculation in each time step; the tangential component and the normal component of the velocity are separated, and the update of the two components is calculated separately.

[0020] Further, the Lagrangian particles are sampled according to the Poisson sampling method and the maximum number of ice crystals allowed in the current simulation scene during initialization to obtain an initial set of crystal center particles, each of which represents an ice crystal and contains a plurality of Lagrangian particles, the initial values of the phase field attributes of the Lagrangian particles are set to specific values between 0 and 1 that allow ice crystal growth, and the initial values of the phase field attributes of the Lagrangian particles that do not belong to the set of crystal center particles are set to 0; the thickness attribute is updated at each time step according to the surrounding Euler particles, and the thickness attribute is interpolated onto the surface of the thin film mesh reconstructed in claim 3 during rendering; the ice crystal mask is calculated according to the phase field attribute to determine the fluid part and the ice crystal part, and the ice crystal to which different Lagrangian particles belong is determined according to the connectivity of the ice crystal part; during the position attribute update, the fluid part is directly updated according to its own speed attribute, and the ice crystal part is first calculated to obtain the total momentum, total angular momentum, center of mass, total mass, and total moment of inertia of the ice crystal according to the Lagrangian particles contained in the ice crystal, then the linear velocity and angular velocity of the ice crystal are calculated, and finally the particle velocity is calculated according to the linear velocity, angular velocity, center of mass, and current position attribute of the ice crystal, and the position attribute is updated according to the particle velocity.

[0021] Further, the graph particles are initialized to be uniformly distributed in a square in two-dimensional space, and each ice crystal composed of Lagrangian particles is divided into a sub-region, which is bound to the set of crystal center particles, and the position and binding relationship remain unchanged during the entire simulation process; during initialization, for each sub-region, the part with a phase field value of 0 is the template part, and the part with a phase field value of 0 is the environment part.

[0022] Further, the ice calculation process is to update the template part and its neighborhood environment part based on the phase field, temperature field, and local molecular direction field, and for the environment part, no physical quantity is updated, and after the physical quantity is updated, the template part and the environment part are not divided again according to the phase field threshold, but only the environment part connected to the template part and having a phase field value greater than a certain threshold is reset to the template part, and the remaining part remains unchanged.

[0023] Further, the physical quantity transfer process between the Lagrangian particles and the Euler particles is to use a three-dimensional smoothing kernel function to calculate the interpolation weight.

[0024] Further, the physical quantity transfer process between the Lagrangian particles and the graph particles is to use a two-dimensional smoothing kernel function to calculate the interpolation weight on the tangent plane defined by the local coordinate system defined according to the normal vector of the graph particle; when the Lagrangian particles transfer physical quantities to the graph particles, only the environment part of the graph particles can obtain the physical quantities; when the graph particles transfer physical quantities to the Lagrangian particles, the Lagrangian particles only obtain the physical quantities from the template part of the graph particles.

[0025] The thin film icing dynamics simulation method based on the multi-layer particle representation can more realistically reproduce the dynamics process of two-dimensional and three-dimensional thin films in cold and room temperature environments, maintains the structured pattern of ice crystals, restores the movement mode of ice crystals, and obtains high-quality thin film icing simulation results.

[0026] Advantages and beneficial effects of the present application:

[0027] The present application realizes the thin film icing dynamics simulation method based on the multi-layer particle representation. The method can effectively maintain the structured pattern of icing in the ice crystal movement process; the method can generate the Marangoni flow generated by the temperature gradient; the method can maintain the stability of the simulation system under the conditions of ice crystal movement, temperature change and complex boundary conditions; the method supports direct user control of the pattern of ice crystals. BRIEF DESCRIPTION OF DRAWINGS

[0028] Figure 1 The figure is a schematic diagram of the thin film icing dynamics simulation method based on the multi-layer particle representation in the present application.

[0029] Figure 2 The figure is a schematic diagram of the effect of the three-dimensional spherical bubble thin film icing process simulation at-40 DEG C in the present application.

[0030] Figure 3 The figure is a schematic diagram of the effect of the three-dimensional spherical bubble thin film icing process simulation at-40 DEG C in the present application.

[0031] Figure 4 The figure is a schematic diagram of the effect of the two-dimensional five-star-shaped thin film icing process simulation at-40 DEG C in the present application.

[0032] Figure 5 The figure is a schematic diagram of the effect of the three-dimensional spherical bubble thin film icing process simulation at-40 DEG C in the present application. DETAILED DESCRIPTION

[0033] In order to make the purpose, technical scheme and advantages of the present application clearer, the technical scheme in the present application will be described clearly and completely below in combination with the drawings in the present application. Obviously, the described embodiments are part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present application. The following embodiments are used to illustrate the present application, but cannot be used to limit the scope of the present application.

[0034] The application provides a thin film icing dynamics simulation method based on a multi-layer particle representation.

[0035] In some embodiments of the application, the thin film icing dynamics simulation method captures the potential thermodynamic mechanism in the two-dimensional and three-dimensional thin film icing process by using three layers of particles, namely Euler particles, Lagrange particles and graph particles. The Euler particles and the Lagrange particles jointly represent the two-dimensional and three-dimensional thin film surface, and the main dynamics calculation process is performed on the Euler particles. The Lagrange particles and the graph particles jointly represent the ice crystals, and the main icing calculation process is performed on the graph particles. In the three layers of particles, the Lagrange particles are used to represent the material state and visualization, and are connected with the graph particles and the Euler particles. The current phase field attributes and other physical quantities are transmitted to the graph particles, and the thermodynamic physical quantities are collected from the graph particles. Then, the dynamics physical quantities and the thermodynamic physical quantities are transmitted to the Euler particles, and the momentum is collected from the Euler particles. The three kinds of particles define local coordinate systems according to their normal vectors, and are used to project the neighborhood particles onto the tangent plane of the particles.

[0036] The Euler particles are distributed relatively uniformly, and are redistributed according to the number density of the Euler particles in each time step calculation to ensure that the number of the Euler particles per unit area is relatively consistent at different positions of the thin film. The Euler particles are used for thin film mesh surface reconstruction in the rendering process.

[0037] The Lagrangian particle is sampled according to a Poisson sampling method and a maximum ice crystal number allowed in a current simulation scene to obtain an initial crystal center particle group, each crystal center particle group represents an ice crystal and contains a plurality of Lagrangian particles, the initial value of the phase field attribute of the Lagrangian particles is set to a specific value between 0 and 1 that can allow ice crystal growth, and the initial value of the phase field attribute of the Lagrangian particles that do not belong to the crystal center particle group is set to 0; the thickness attribute is updated according to surrounding Euler particles at each time step, and the thickness attribute is interpolated to the surface of the film mesh body reconstructed in claim 3 during rendering; the ice crystal mask is calculated according to the phase field attribute to determine the fluid part and the ice crystal part, and the ice crystal to which different Lagrangian particles belong is determined according to the connectivity of the ice crystal part; during the position attribute update, the fluid part is directly updated according to the speed attribute, the total momentum, the total angular momentum, the center of mass, the total mass and the total rotational inertia of the ice crystal are calculated according to the Lagrangian particles contained in the ice crystal, the linear velocity and the angular velocity of the ice crystal are calculated, finally, the particle speed is calculated according to the linear velocity, the angular velocity, the center of mass and the current position attribute of the ice crystal, and the position attribute is updated according to the particle speed.

[0038] The graph particle is initialized as a square uniform distribution in a two-dimensional space, a sub-region is divided for each ice crystal composed of the Lagrangian particle, and the sub-region is bound with the crystal center particle group, and the position and the binding relationship remain unchanged during the whole simulation process; during the initialization, for each sub-region, the part with the phase field value of 0 is the template part, and the part with the phase field value of 0 is the environment part.

[0039] In some embodiments of the present application, the dynamics calculation process in the thin-film icing dynamics simulation method introduces a temperature equation and re-derives an update formula of the thin-film surfactant concentration based on an implicit incompressible smoothed particle hydrodynamics (hereinafter referred to as IISPH), and the relaxation Jacobi iteration method is used for calculation at each time step; the tangential component and the normal component of the velocity are separated and the update of the two components is calculated separately.

[0040] In some embodiments of the present application, the icing calculation process in the thin-film icing dynamics simulation method, for the template part and the environment part of the neighborhood of the template part during the initialization of the graph particle, the update is based on the phase field, the temperature field and the local molecular direction field, for the environment part, the update of the physical quantity is not performed, after the update of the physical quantity, the template part and the environment part are not divided again according to the phase field threshold value, but only the environment part connected with the template part and having a phase field value greater than a certain threshold value is reset to the template part, and the remaining part remains unchanged.

[0041] In some embodiments of the present application, the process of physical quantity transmission between Lagrangian particles and Eulerian particles in the thin-film icing dynamics simulation method uses three-dimensional smooth kernel function for interpolation weight calculation.

[0042] In some embodiments of the present application, the process of physical quantity transmission between Lagrangian particles and Eulerian particles in the thin-film icing dynamics simulation method uses three-dimensional smooth kernel function for interpolation weight calculation.

[0043] In some embodiments of the present application, the boundary particle thin-film surfactant concentration calculation formula based on the moving least square method is that, during initialization, boundary particles are sampled at the position where the solid boundary is close to the Eulerian particle, and the particle spacing is consistent with the initial spacing of the Eulerian particle. Before each iteration of each relaxation Jacobi iteration of the IISPH, the sampled boundary particle estimates its surfactant concentration using the moving least square method according to the surfactant concentration of the neighboring Eulerian particles.

[0044] Reference Figure 1 A thin-film icing dynamics simulation method based on multi-layer particle representation provided by the present application is described in detail, including the following steps:

[0045] First, the ice crystal growth and shape maintaining mechanism based on the graph particle. Only the phase field, temperature field and local molecular direction field on the Lagrangian particle are transmitted to the environment area of the graph particle, and only the physical quantities in the template area of the graph particle and the environment area in the neighborhood thereof are updated, and then only the physical quantities in the template area of the graph particle are transmitted to the ice crystal bound by each sub-area and the Lagrangian particles in the neighborhood thereof. The specific implementation includes the following steps:

[0046] a1. Phase field mapping. For each ice crystal C, according to each Lagrangian particle L belonging to the ice crystal and the two-dimensional rotation angle β L of the current state of L relative to the stored state on the graph particle M, the two-dimensional rotation angle β C of the ice crystal C is calculated in the following manner:

[0047]

[0048] Then, for each Lagrangian particle L that may need to map the position (i.e. has been iced or is about to be iced), if its position is close enough to the ice crystal C, the rotation matrix R(β C of β C is used to rotate the phase field of the Lagrangian particle L to the phase field of the ice crystal C.) and the center position of the sub-region on the graph particle M bound to ice crystal C. Calculate the position on the sub-region after mapping. The calculation formula is as follows:

[0049]

[0050] Where T C (v) represents projecting vector v onto the two-dimensional plane containing ice crystal C, x L and x C These are the positions of the Lagrange particle L and the ice crystal C, respectively.

[0051] The a2.L2M mapping, for thermodynamic physical quantities q (i.e., phase field ζ and temperature field T), defines the following transfer process:

[0052]

[0053] in Represents the set of Lagrange particles. This represents the graph particle M from the set of Lagrange particles in its neighborhood. Physical quantities are obtained at the graph particle M. W(M,L) represents the smooth kernel function between the graph particle M and the Lagrange particle L. According to the above definition, the relevant physical quantities on the graph particle M are the phase field ζ, the temperature field T, and the local molecular orientation field θ. ori You can obtain it in the following way:

[0054]

[0055] Where L n It is the Lagrange particle closest to the graph particle M.

[0056] a3. Icing Calculation. For each graph particle M in the template part of the graph particle and its neighborhood environment part, the phase field ζ is updated as follows (the subscript M is omitted for ease of reading):

[0057]

[0058] Among them, M ζ It is a constant controlling the freezing rate, θ is the local freezing direction, (u,v) are the coordinates of the graph particle in the two-dimensional plane, and the anisotropy function ε, the double potential well function g, the weighting function p, and the temperature-dependent pure solid free energy density f are all defined. s and the free energy density of pure liquid phase f l Direction-dependent free energy density f ori The definition of an isofunction is given by the following formula:

[0059]

[0060] wherein, is a constant controlling the width of the fluid- solid interface, δ is a constant controlling the strength of anisotropy, J is a constant controlling the number of main branches of ice crystals, α is a supercooling coefficient in the range of (0, 1), γ is a thermal scaling coefficient, T e is the melting temperature.

[0061] Similarly, the update formula of the local molecular orientation field θ ori is:

[0062]

[0063] wherein M ori is a constant describing the variability of the local molecular orientation field.

[0064] According to the thermal diffusion coefficient a, latent heat coefficient K, ambient temperature T env and time constant τ, the update of temperature T can be performed in the following way:

[0065]

[0066] a4. M2L transfer, similar to L2M transfer, for each Lagrangian particle L, define its way of obtaining physical quantities from graph particles as:

[0067]

[0068] wherein C is the ice crystal to which L belongs, is the set of graph particles, denotes the set of neighborhood graph particles that L obtains from the mapped position of L, denotes that L obtains physical quantities from , M n denotes the graph particle closest to the position . In particular, for the un-frozen Lagrangian particles, obtain physical quantities from the sub-area of graph particles corresponding to the nearest ice crystal according to the above formula.

[0069] The above steps a1-a4 represent the calculation of the natural icing process according to the phase field theory. If it is desired to design the ice crystal pattern in a user-controlled manner, there are two ways: control by pre-designed grayscale map editing and fixing the local molecular orientation field of the graph particles, or replace the phase field update formula with a normalized time-varying directed distance field, and calculate the update of the temperature field according to the update of the phase field. It should be noted that, in order to prevent visual flicker, the directed distance field should be non-decreasing over time.

[0070] Second, the momentum, mass, volume, temperature, and phase field of the Lagrangian particles are transferred to the Eulerian particles, and the geometry and dynamics are calculated according to the distribution of the Eulerian particles. In the calculation process, the moving least square method is used to estimate the surfactant concentration on the boundary particles. The specific implementation includes the following steps:

[0071] b. L2E transfer, for each Eulerian particle E, from the set of Lagrangian particles in its neighborhood The way to obtain an arbitrary physical quantity q at

[0072]

[0073] If ε is defined as the set of Eulerian particles, then the calculation of a in the above formula is as follows: L

[0074]

[0075] According to the above formula, the mass m, surfactant c, volume V, and momentum p are transferred from the Lagrangian particles to the Eulerian particles. In particular, the temperature is transferred in the following way: In addition, the affine momentum

[0076]

[0077] The calculation of the matrix B L and D L will be given in the E2L transfer step of step d below. The velocity of the Eulerian particle E

[0078] c. Film calculation, the film calculation is divided into geometry and dynamics.

[0079] In the geometry calculation, first update the thickness η of each Lagrangian particle L in the following way: L

[0080]

[0081] For each Eulerian particle E, the calculation of its thickness is η E = V E / a E , and the calculation of the surfactant concentration is Γ E = c E / a E , where a E is the control area of E, and the calculation is as follows:

[0082]

[0083] E​​' E denotes other Euler particles in the neighborhood of E. If we define E ' In the local coordinate system of E, we have coordinates (u, v, z), is a two-dimensional differential operator on the tangent plane of E, then the mean curvature at E can be calculated in the following way:

[0084]

[0085] Before the dynamic calculation, we first divide the velocity of Euler particles into normal component and tangential component and update them respectively.

[0086] In the normal direction, we first define the external pressure p out of the film as atmospheric pressure. If the film is not closed, we also set the internal pressure p in as atmospheric pressure, and if the film is closed, we calculate it in the following way:

[0087]

[0088] where n0 is the molar mass of the internal gas, is the ideal gas constant, O is an arbitrary point inside the film, and n E is the normal of the film at E.

[0089] Then, we can update the normal velocity component according to the following formula:

[0090]

[0091] where p is the density, σ0 is the surface tension of pure water, is the component of external force in the tangential direction of E.

[0092] In the tangential direction, we first need to solve the new surfactant concentration Γ E using the relaxation Jacobi iteration method, the formula is as follows (the subscript E has been omitted for easy reading):

[0093]

[0094] where T * , Γ * and η * represent the temperature, surfactant concentration and thickness at E after the L2E transfer in step b respectively. After solving, we can update the tangential component of the velocity E

[0095]

[0096] ​In particular, when there is a solid boundary in the scene, the particles need to be sampled on the boundary, and the surfactant concentration on the boundary particles B is calculated where and are x B are the first two components of d B , while d B is calculated as follows:

[0097]

[0098] where denotes the coordinates of an Euler particle E in the B neighborhood in the local coordinate system R of B.

[0099] α B , β B , and γ B are calculated as follows:

[0100]

[0101] where

[0102] Third, the momentum on the Euler particles is transferred to the Lagrangian particles, and the Lagrangian particles are divided into fluid and ice crystal parts according to the phase field value, and the two are updated in different ways. The specific implementation includes the following steps:

[0103] d. E2L transfer, for each Lagrangian particle L, the velocity u L and the matrix B L and the matrix D L are obtained from the neighborhood Euler particles as follows:

[0104]

[0105] e. Fluid-structure coupling calculation. For Euler particles, there is no clear distinction between solid and liquid phases, and no collision between ice crystals and boundaries is considered. For Lagrangian particles, the solid and liquid phases are distinguished by a specific threshold value according to the phase field value, and once the phase field of a liquid particle exceeds the threshold value, it is converted from liquid to solid.

[0106] In solid processing and fluid-structure coupling, it is assumed that each ice crystal is a rigid body, and if the relative speed of adjacent ice crystals is small enough, they will be merged. For each ice crystal, the total mass, total momentum, and center of mass are first calculated, and the total angular momentum and total moment of inertia are further obtained through the center of mass, and finally the linear velocity, angular velocity, and final velocity of the Lagrangian particle of the rigid body are calculated. In addition, before the solid processing starts, a number of XSPH calculations are first performed to apply a virtual viscosity to ensure that there is no gap between the fluid and the solid.

[0107] f. Position update, for all Euler particles, first position update is performed, then according to the number density of the particles, a virtual pressure is calculated (the calculation process refers to the iterative calculation method of the surfactant concentration on the Euler particles), and the particle position is updated again through the virtual pressure, that is, the redistribution step, which ensures the relatively uniform particle distribution. For all Lagrange particles, first position update is performed according to the calculated velocity, and then the position of the Lagrange particle is projected onto the surface represented by the Euler particles.

[0108] Based on the above embodiments, the film icing dynamics simulation method based on multi-layer particle representation of the present application realizes the simulation of the film icing phenomenon, has the advantages of high stability in complex situations and supporting user control of ice crystal shape.

[0109] Through verification, the method of the present application can more realistically reproduce the dynamics process of two-dimensional and three-dimensional films icing in cold and room temperature environments in a computer, maintain the structured pattern of ice crystals, and restore the movement mode of ice crystals, obtaining high-quality film icing simulation results. For details, see the following simulation results: Figures 2-5 .

[0110] Referring to Figure 2 , according to the reference image, the simulation method proposed by the present application can obtain high-quality simulation results of three-dimensional spherical film (soap bubble) icing, and can realize asymmetric ice crystal growth.

[0111] Referring to Figure 3 , according to the reference image, the simulation method proposed by the present application can fix the local particle direction field according to the gray scale image (the rightmost image), and realize the user-controlled spiral icing pattern.

[0112] Referring to Figure 4 , according to the reference image, the simulation method proposed by the present application can perform stable two-dimensional film simulation in the case of complex five-point star-shaped boundary.

[0113] Referring to Figure 5 , according to the reference image, the simulation method proposed by the present application can also be applied to the simulation of the icing process of a three-dimensional hemispherical film at room temperature (20℃), and obtain simulation results consistent with the laws of thermodynamics.

[0114] It should be further pointed out that the above embodiments are only used to understand the technical solutions of the present application, and are not used to limit the protection scope of the present application. Any obvious adjustment and modification of the technical solutions of the present application which belongs to the technical concept of the present application should also belong to the protection scope of the present application.

Claims

1. A method for thin film icing dynamics simulation based on multi-layer particle representation, characterized in that The method comprises the following steps: a. A theoretical framework for capturing the potential thermodynamic mechanisms in the process of two-dimensional and three-dimensional thin film icing is proposed; The potential thermodynamic mechanisms in the process of two-dimensional and three-dimensional thin film icing are simulated by using three layers of particles, namely Euler particles, Lagrange particles and graph particles, wherein the Euler particles and the Lagrange particles jointly represent the two-dimensional and three-dimensional thin film surfaces, and the main dynamic calculation process is performed on the Euler particles; the Lagrange particles and the graph particles jointly represent ice crystals, and the icing calculation process is performed on the graph particles; in the three layers of particles, the Lagrange particles are used to represent the material state and visualization, and are connected with the graph particles and the Euler particles, so as to transmit the current phase field attributes and other physical quantities to the graph particles, collect thermodynamic physical quantities from the graph particles, and then transmit the dynamic physical quantities and the thermodynamic physical quantities to the Euler particles and collect momentum from the Euler particles; the three kinds of particles all define local coordinate systems according to their normal vectors, so as to project the neighboring particles onto the tangent plane of the particles; b. An advanced phase field mapping method is proposed to maintain the pattern of ice crystal structuring under the mobile particle representation, and to support user direct control over the pattern during simulation; The support for the user direct control over the pattern during simulation is realized by pre-designed grayscale map editing and fixed local particle direction field control of the graph particles, or by directly updating the phase field through the normalized time-varying directed distance field, and calculating the update of the temperature field according to the update of the phase field; c. A boundary particle thin film surfactant concentration calculation formula based on the moving least square method under the multi-layer particle representation is derived to improve the stability of the simulation system; The boundary particle thin film surfactant concentration calculation formula based on the moving least square method is that, during initialization, the boundary particles are sampled at the position where the solid boundary is close to the Euler particles, and the particle spacing is consistent with the initial spacing of the Euler particles; before each iteration of each relaxation Jacobi iteration of the implicit incompressible smoothed particle hydrodynamics (IISPH), the sampled boundary particles estimate their own surfactant concentration by using the moving least square method according to the surfactant concentration of the neighboring Euler particles.

2. The simulation method of claim 1, wherein: In the three layers of particles, the Lagrange particles are used to represent the material state and visualization, and are connected with the graph particles and the Euler particles, so as to transmit the current phase field attributes and other physical quantities to the graph particles, collect thermodynamic physical quantities from the graph particles, and then transmit the dynamic physical quantities and the thermodynamic physical quantities to the Euler particles and collect momentum from the Euler particles.

3. The simulation method of claim 1, wherein: The Euler particles are distributed relatively uniformly, and are redistributed according to the number density of the Euler particles in each time step to ensure that the number of the Euler particles per unit area is relatively consistent at different positions of the thin film; and the Euler particles are used for surface reconstruction of the thin film mesh body in the rendering process.

4. The simulation method of claim 1, wherein: The dynamic calculation process introduces a temperature equation, and re-derives an update formula of the thin film surfactant concentration based on the implicit incompressible smoothed particle hydrodynamics (IISPH), and the formula is calculated by using the relaxation Jacobi iteration method in each time step; the tangential component and the normal component of the velocity are separated, and the update of the two components is calculated separately.

5. The simulation method of claim 1, wherein: The Lagrangian particles are sampled according to a Poisson sampling method and a maximum number of ice crystals allowed in a current simulation scene to obtain an initial set of crystal center particles, each set of crystal center particles representing an ice crystal and containing a plurality of Lagrangian particles, the phase field attribute of the Lagrangian particles being initially set to a specific value between 0 and 1 to allow ice crystal growth, and the phase field attribute of the Lagrangian particles not belonging to the set of crystal center particles being initially set to 0; the thickness attribute is updated according to surrounding Euler particles at each time step, and the thickness attribute is interpolated to the surface of the film mesh reconstructed in claim 3 during rendering; the ice crystal mask is calculated according to the phase field attribute to determine the fluid part and the ice crystal part, and the ice crystal to which different Lagrangian particles belong is determined according to the connectivity of the ice crystal part; During the position attribute update, the fluid part is directly updated according to the speed attribute, and the ice crystal part is first calculated according to the total momentum, total angular momentum, center of mass, total mass and total moment of inertia of the ice crystal containing the Lagrangian particles, and then the linear velocity and angular velocity of the ice crystal are calculated, and finally the particle velocity is calculated according to the linear velocity, angular velocity, center of mass and current position attribute of the ice crystal, and the position attribute is updated according to the particle velocity.

6. The simulation method of claim 1, wherein: The graph particles are initialized as a square uniform distribution in two-dimensional space, and each ice crystal composed of Lagrangian particles in claim 5 is divided into a sub-region, which is bound to the set of crystal center particles in claim 5, and the position and binding relationship remain unchanged during the entire simulation process. During initialization, for each sub-region, the part with a phase field value of 0 is the template part, and the part with a phase field value of 0 is the environment part.

7. The simulation method of claim 1, wherein: The ice calculation process is to update the template part and its neighborhood environment part based on the phase field, temperature field and local molecular direction field, and for the environment part, no physical quantity is updated, and after the physical quantity is updated, the template part and the environment part are not divided again according to the phase field threshold, but only the environment part connected with the template part and having a phase field value greater than a certain threshold is reset to the template part, and the remaining part remains unchanged.

8. The simulation method of claim 1, wherein: The physical quantity transmission process between the Lagrangian particles and the Euler particles is to use a three-dimensional smoothing kernel function to calculate the interpolation weight.

9. The simulation method of claim 1, wherein: The physical quantity transmission process between the Lagrangian particles and the graph particles is to use a two-dimensional smoothing kernel function to calculate the interpolation weight on the tangent plane defined according to the normal vector of the graph particles in claim 2; only the environment part of the graph particles in claim 6 can obtain the physical quantity when the Lagrangian particles transmit the physical quantity to the graph particles; and only the template part of the graph particles in claim 6 can obtain the physical quantity when the graph particles transmit the physical quantity to the Lagrangian particles.

Images