A simulation method for phase transition of elastic-viscoplastic material based on material point method

Abstract

The application discloses a simulation method of elastic-viscoplastic material phase transition based on a material point method, which is characterized by the following steps: 1) a unified numerical model for describing the elastic-viscoplastic material and non-Newtonian fluid phenomenon is established for simulating the material in the phase transition process; 2) Allen-Cahn and Cahn-Hilliard phase field equations are discretely modeled and introduced into the material point method; and 3) an elastic-viscoplastic material control strategy driven by phase field evolution is introduced to eliminate the discontinuity in performance caused by the change of physical properties of the mixed contact surface in the process. Compared with the prior art, the application has the advantages that the simulation of the three states of solid, intermediate and liquid of the complex elastic-viscoplastic material is simulated, the details generated by the evolution between phases are simulated, a new flexible phase field evolution based on the material point method and the corresponding material performance can be seamlessly and controllably displayed, and thus a realistic phase transition simulation effect is realized, and the application has a unique advantage in the generation of complex phase transition animation.

Description

Technical Field

[0001] This invention relates to the field of animation simulation technology, specifically a simulation method for the phase transition process of elastoviscoplastic materials based on the material point method and phase field model. Background Technology

[0002] Realistic animations depicting the coexistence of fluids and solids, and their subtle interactions, play an indispensable role in creating stunning visual effects in today's video games, digital entertainment, and virtual reality. Despite recent successes in simulating multiphase flows, viscoelastic / viscoplastic materials, and elastic and plastic deformation, the fundamental theory and numerical description of complex multiphase interactions, as well as the temporal evolution of intermediate states in fluid-solid mixtures during phase transitions—e.g., the transformation of a boiled egg from a fluid to a solid, or cement mixtures—remain largely unexplored. Rapid and efficient simulation and its adoption for mass production in the film and television industry remain major challenges for computer graphics.

[0003] To seamlessly handle the elastic-viscoplastic (EVP) and fluid behavior during phase transitions and to simulate their interactions in a unified manner, an EVP physical model from continuum mechanics is needed. This model can represent particles / elastic / viscoplastic and even fluids in a unified way and supports more complex phenomena than existing simulation frameworks, such as multiphase and multimaterial interactions. Simultaneously, to correctly address interface evolution in multiphase interactions, phase field theory models (Allen-Cahn and Cahn-Hilliard phase fields) are introduced. The basic idea is to use an auxiliary scalar field where specific phase field equations, such as diffusion equations, can be defined and discretized. The introduction of phase field theory and the application of phase field models provide a stable evolution method from an energy perspective for phase diffusion and separation phenomena at solid-liquid mixture interfaces, as well as other possible phase transition phenomena.

[0004] Phase-field models often require higher-order stability and accuracy, and how to discretize and solve them more effectively and stably has always been a key research problem. The Material Point Method (MPM), as a particle-mesh hybrid method, possesses numerical capabilities from both Lagrange and Euler perspectives, providing stable explicit or implicit discretization solutions for higher-order partial differential equations, such as phase-field models including the Cahn-Hilliard and Allen-Cahn equations. Furthermore, due to its simplicity and stability, the material point-based abstraction method of the MPM is beneficial for the numerical realization of physical constitutive models such as the elastoviscoplastic model.

[0005] In conclusion, combining the matter point method with the phase field model can provide a strong foundation for realistic simulation of phase transition processes, and will greatly enhance the visual experience of related simulations. Summary of the Invention

[0006] The purpose of this invention is to address the shortcomings of existing technologies by proposing a simulation method for the phase transition process of elastoviscoplastic materials. This method employs a realistic material phase transition simulation based on the matter point method and incorporating a phase field model. It realizes the simulation of multiphase, multi-material interaction processes existing in real life, as well as the elastoviscoplastic behavior of intermediate products. This method introduces phase field theory, improving the phase transition phenomena generated around the interface during multiphase material interaction from both a theoretical and practical perspective. This results in more realistic simulation results, a simpler method, and better performance. Applied to various phase transition phenomena, i.e., multi-material scenarios, it can seamlessly and controllably demonstrate a new flexible phase field evolution based on the matter point method and the corresponding material behavior, thereby achieving realistic phase transition simulation effects and showing promising application prospects.

[0007] The specific technical solution to achieve the purpose of this invention is: a simulation method for the phase transition process of elastoviscoplastic materials based on the matter point method. Its characteristic is that this method employs a unified simulation method for elastoviscoplastic materials and non-Newtonian fluids based on the matter point method, and the discretization of the Allen-Cahn and Cahn-Hilliard phase field equations and their matter point method. It utilizes a phase field-driven dynamic control strategy for the properties of elastoviscoplastic materials to achieve a realistic simulation of the phase transition process of elastoviscoplastic materials. Specifically, it includes the following steps:

[0008] 1) A unified simulation method for elasto-viscoplastic materials and non-Newtonian fluids based on the matter point method, specifically including:

[0009] a) Particle properties are transferred to the background mesh

[0010] First, the simulated target object is discretized into a component carrying a certain mass m. p and volume The simulation involves particles and a background mesh G that covers the simulation area. The background mesh is then reinitialized (all carried values ​​are set to zero), and the mass and momentum carried by the particles are considered according to the APIC form and the spline interpolation function w. ip The transfer to a similar background grid is specifically expressed by the following formulas (a1) to (a2):

[0011]

[0012]

[0013] in, The local velocity affine used in APIC ensures that the system possesses better momentum and angular momentum conservation properties; m n Represents the mass at time step n (subscript i represents the attribute carried by the mesh, and subscript p represents the particle attribute); This represents the velocity at time step n.

[0014] b) Update background mesh point momentum

[0015] This step requires considering the strain (deformation gradient) of the particles surrounding the grid points. The stress is then calculated according to the elastic constitutive model of the system, and the corresponding stress is solved by equation (b) based on the rate at which the corresponding mesh points are updated according to the stress.

[0016]

[0017] The velocity of the corresponding grid point is solved by the following equation (c):

[0018]

[0019] Where Ψ is the energy density function defined by the elastic constitutive model; Represents the initial volume of the particle; The deformation gradient of particle p; w ip For the spline interpolation function between grid point i and particle p; This represents the velocity of the grid point in the next time step; Δt is the time step size. The impact of collisions on the grid velocity also needs to be addressed; different normal and tangential velocities at the collision point should be handled differently depending on the type of collision.

[0020] c) Background mesh properties are passed back to particles and deformation gradients are updated.

[0021] Based on the updated background mesh, the updated physical properties are transmitted back to the surrounding particles using a spline interpolation function. Simultaneously, the local velocity affine matrix is ​​updated for APIC, as expressed by equations (d) to (f) below:

[0022]

[0023]

[0024]

[0025] in, This represents the velocity of the particle at the next time step; For the local velocity affine used in APIC in the next time step; Represents the velocity of the grid point at the next time step; Δt represents the position of particle p at the next time step (updated from the position at the current time step); Δt represents the time step size.

[0026] The particle then updates its strain (deformation gradient) based on the velocity gradient of the surrounding grid points, resulting in the test strain F defined by the following equation (g). tr :

[0027]

[0028] Among them, F tr The deformation gradient is a prediction of the degree of deformation at the next time step, which is used for subsequent corrections. The velocity gradient at particle p passes through The deformation gradient is estimated using this method; other parameters are explained in the formula above. This is a deformation gradient predicted based on velocity changes. Further correction is needed based on the viscoplastic constitutive model of the material, using a return mapping algorithm for deformation gradients that violate the yield criterion.

[0029] To broadly describe materials of various phases, we combined the Drucker-Prager and Cam-Clay models with the von Mises model, respectively, and established two yield criteria (VMCC and VMDP) to describe the solid / intermediate state. These criteria are defined by equations (h) to (i) below:

[0030] y vmdp =C f tr(τ)+||s||-C c (h);

[0031]

[0032] Among them, y vmdp y vmcc These are the forms obtained by combining the yield criteria of von Mises with those of Drucker Prager and Cam Clay, respectively; C f C c These represent the friction angle and degree of polymerization of the material, respectively; τ and s = dev(τ) are the Kirchhoff stress and shear stress obtained from the deformation gradient using the elastic constitutive model, i.e. When y < 0, the material point violates the yield criterion, and a back-mapping algorithm is needed to correct it using the backward Euler form expressed by the following equation (j):

[0033] b n+1 -b tr =-2δγG(τ) n+1 )b n+1 (j).

[0034] in, It is a non-associative return mapping rule with volume-preserving properties; δγ=γΔt is the plastic flow distance; b=FF T Let b be the left Cauchy Green tensor. tr =F tr F tr,T .

[0035] c-4: According to the Herschel-Bulkley model, the viscosity coefficient η and shear rate coefficient h are used to control the non-Newtonian behavior of the material points during the return mapping process. When h < 1, it exhibits shear thinning, and when h > 1, it exhibits shear thickening. The specific formula is expressed by the following equation (k):

[0036]

[0037] in, d represents the number of dimensions; y vm This is the corresponding yield criterion function.

[0038] c-4: Solve for s using the bisection method n+1 Then, the corrected F can be finally obtained through the following equation (l). n+1 That is, the particle deformation gradient at the next time step:

[0039]

[0040] 2) Discrete modeling of the Allen-Cahn and Cahn-Hilliard phase-field equations is performed, and the matter point method is introduced. The phase transition caused by material contact and external influences is simulated through the evolution of the phase field. Specifically, this includes:

[0041] a) Introducing the phase field as an additional field into the computational flow of the matter point method, in step a) above, the transfer of particle properties to the background mesh increases the transmission of the phase field between particles and the mesh. The specific formula is expressed by the following equation (m):

[0042]

[0043] Among them, c n This represents the phase field value at time step n (the subscript i represents the attribute carried by the mesh, and the subscript p represents the particle attribute).

[0044] b) Based on the weak form of the Allen-Cahn (AC) and Cahn-Hilliard (CH) equations, the phase field solution is added during the grid update of the background grid point momentum in step b) above, to update the phase field value of the next time step. The matrix form of the specific solution for the CH phase field is expressed by the following equation (n):

[0045]

[0046] The matrix form for solving the AC phase field is expressed by the following equation (o):

[0047]

[0048] Where D is the diffusion rate of the phase field; ∈ is the width coefficient of the phase field interface; M, r, L are the discrete matrices of the mass moment, potential moment, and Laplace operator of the phase field in the form of the matter point method, respectively; c n+1 Let be the column vector formed by the phase field values ​​of all grid points. Let be the column vector of the free energies of all grid points;

[0049] The mass moment M, potential moment r, and discrete L of the phase field in the form of the matter point method are matrices represented by the following equations (p) to (r):

[0050]

[0051]

[0052]

[0053] Here, H is the user-defined phase field potential energy function, and the meanings of the other parameters have been explained above.

[0054] 3) Introduce a phase-field evolution-driven dynamic property control strategy for elastoviscoplastic materials to eliminate the discontinuities in performance caused by changes in the physical properties of the mixed contact surface. Specifically, this includes:

[0055] a) Based on the phase field value c carried by each particle, the range of phase field values ​​is divided into solid phase: [0, c a ], intermediate phase: [c a ,c b ], liquid phase [c b [,1], where 0 represents a completely dry solid, 1 represents a completely liquid, and c is set. m To represent the state where the mesophase has the highest viscosity; c a ,c b The parameters set by the user are used to control the phase range, and in experiments they are generally set to 0.3 and 0.7 respectively.

[0056] b) Smooth interpolation between e1 and e2 is performed using Hermite polynomials, the function S1 of which is defined by the following equation (s):

[0057]

[0058] Where e1 and e2 are the starting and ending points of the gradual change, respectively; x is the independent variable.

[0059] c) Define the smooth function S2, which is the opposite of S1, by the following equation (t):

[0060] S2(e1,e2,x)=S1(-e2,-e1,-x)(t).

[0061] d) For VMDP, the parameters as a function of the phase field are defined by the following equations (u) to (v):

[0062]

[0063]

[0064] Where, φ c C is the friction angle; c This is the degree of polymerization coefficient; These are the user-defined initial friction angle of the solid and the maximum polymerization coefficient during the phase transformation process, respectively, which represent the magnitude of the yield stress; c is the phase field value of the corresponding particle.

[0065] e) Set the polymerization coefficient and hardening coefficient according to the following equations (w) to (x):

[0066] C c =κsinh(ξmax(-α,0))+C c (w);

[0067] ξ=S2(0,c m ,c)ξ 0 (x).

[0068] Where, ξ 0 c is the maximum value of the hardness coefficient. m This is the location of maximum viscosity; as the phase field value c approaches c... m The material's hardening ability gradually decreases to 0 while its viscosity η gradually increases.

[0069] f) The change in viscosity during a phase transition as defined by the following equation (y):

[0070] η=min(S1(c a ,c m ,c),S2(c m ,c b ,c))(y).

[0071] Where η is the dynamic viscosity coefficient during the phase transition process, and the definitions of other parameters can be explained above.

[0072] Compared with the prior art, the present invention has the following beneficial technical effects and significant technical progress:

[0073] 1) This invention explores and expands state-of-the-art physics-driven simulation and phase-field modeling techniques to establish a completely unified framework for simulating multiphase interactions with fine interface effects, and can fully demonstrate the phase transition process of intermediate phases, providing a stable and realistic simulation framework for seamless phase transition processes.

[0074] 2) Based on the powerful material point method, this invention realizes a unified elastic-viscoplastic constitutive model, and introduces a discretized phase field model into it, which ensures the accuracy of the phase transition process and the realism of the material performance from the perspective of physical theory.

[0075] 3) When applied to various phase transition phenomena, i.e. multiple material scenarios, this framework can seamlessly and controllably demonstrate a new flexible phase field evolution based on the material point method and the corresponding material behavior, thereby achieving realistic phase transition simulation effects. Attached Figure Description

[0076] Figure 1 This is a schematic diagram illustrating the process by which water and powder are mixed to produce a viscous intermediate substance, as described in this invention.

[0077] Figure 2 This is a schematic diagram illustrating the simulation of multiphase materials and the phase transition process between them in this invention;

[0078] Figure 3 This is a schematic diagram illustrating the diffusion separation evolution of the AC and CH phase field models under random conditions to verify the present invention.

[0079] Figure 4 This is a schematic diagram simulating the melting process of candy by heat according to the present invention;

[0080] Figure 5 This is a schematic diagram illustrating the cooling process of melted candy after it comes into contact with the ground, as described in this invention.

[0081] Figure 6 This is a schematic diagram illustrating the volume change process of an egg after being heated and undergoing a phase change, as described in this invention.

[0082] Specific implementation methods

[0083] The present invention specifically includes the following steps:

[0084] 1) A unified simulation method for elastic-viscoplastic materials and non-Newtonian fluids based on the material point method.

[0085] By combining the Drucker-Prager and Cam-Clay models with the von Mises yield criterion and introducing non-Newtonian flow theory using the material point method, a unified elastic-viscoplastic model is achieved to describe materials such as particulate solids, elastomers, viscous bodies, and fluids. This allows for the description of material properties that may occur during phase transitions using a unified model, and provides model support for the seamless changes in material properties during phase transitions.

[0086] 2) Discretization of the Allen-Cahn and Cahn-Hilliard phase-field equations and their material point method implementation

[0087] By discretizing the Allen-Cahn and Cahn-Hilliard phase field models and introducing the matter point method, we provide underlying theoretical support for the phase transition process in multiphase contact, realize accurate modeling of the interface evolution and ensure continuity during the phase transition process, and control the evolution trend through a user-defined potential energy function.

[0088] 3) Dynamic property control strategy for elastoviscoplastic materials driven by phase field

[0089] By utilizing the evolution of phase field values ​​to dynamically control the parameters of the constitutive model of elastoviscoplastic materials, it is possible to effectively describe changes in properties such as elasticity, plasticity, and viscosity that may occur during phase transitions. This invention can effectively describe the changes in material properties during the common phase transition process from solid to intermediate to liquid, thereby achieving rich and realistic simulations of phase transition processes.

[0090] The material point method simulation process of the elastic-viscoplastic constitutive model in step 1) is described in detail below:

[0091] First, the simulated target object is discretized into a component carrying a certain mass m. p and volume Particles, and a background mesh G that can cover the simulation area. The background mesh is then reinitialized (deformation gradient F). p Set as the identity matrix), and set the time step for this iteration according to the CFL conditions. Consider assigning the mass and momentum carried by the particle to the APIC form and using the spline interpolation function w. ip The transfer to a nearby background grid is expressed by the following formulas (a1) to (a2):

[0092]

[0093]

[0094] in, The local velocity affine used in APIC ensures that the system has better momentum and angular momentum conservation properties. ip For spline interpolation function; m n Represents the mass at time step n (subscript i represents the attribute carried by the mesh, and subscript p represents the particle attribute); This represents the velocity at time step n.

[0095] The spline interpolation function w ip It provides a smooth interpolation near a spatial location, typically set by the following equation (a3):

[0096]

[0097] After transferring particle properties to the mesh, momentum and velocity need to be updated on the mesh. This step requires considering the strain (deformation gradient) of particles around the mesh points. The system's elastic constitutive model is used to solve for the corresponding stress release, and then the velocity of the corresponding mesh points is updated according to the stress, specifically in the form of equations (b) to (c) below:

[0098]

[0099]

[0100] Where Ψ is the energy density function defined by the elastic constitutive model; Represents the initial volume of the particle; The deformation gradient of particle p; w ip For the spline interpolation function between grid point i and particle p; Δt represents the velocity of the grid point at the next time step; Δt is the time step size.

[0101] The model used in this invention is the Neo-Hookean model modified as shown in equation (z) below:

[0102]

[0103] Where b = FF T Let be the left Cauchy-Green tensor; J = det(F) represents the rate of change of volume; κ is the bulk modulus of the material; μ is the shear modulus of the material; and d is the simulation dimension. The energy is divided into a shear component Ψ. dev (F) and volume part Ψ vol (F) makes subsequent viscoplastic processing easier. After that, the impact of collisions on mesh velocity needs to be addressed; different normal and tangential velocities at the collision points should be handled differently depending on the type of collision.

[0104] Based on the updated background mesh, the updated physical properties need to be passed back to the surrounding particles using the spline interpolation function. At the same time, the local velocity affine matrix is ​​updated for APIC. The specific formulas are expressed by equations (d) to (f) below:

[0105]

[0106]

[0107]

[0108] in, This represents the velocity of the particle at the next time step; For the local velocity affine used in APIC in the next time step; Represents the velocity of the grid point at the next time step; Δt represents the position of particle p at the next time step (updated from the position at the current time step); Δt is the time step size.

[0109] The particle then updates its strain (deformation gradient) based on the velocity gradient of the surrounding grid points, resulting in the test strain F defined by the following equation (g). tr :

[0110]

[0111] Among them, F tr The deformation gradient is a prediction of the degree of deformation at the next time step, which is used for subsequent corrections. The velocity gradient at particle p passes through The method of estimation is as follows; other parameters are explained in the formula above. Trial strain F tr The deformation gradient is predicted based on the velocity change. It needs to be corrected by the return mapping algorithm based on the viscoplastic constitutive model of the material for deformation gradients that violate the yield criterion.

[0112] To broadly describe materials of various phases, the Drucker-Prager and Cam-Clay models were combined with the von Mises model, respectively, to establish two yield criteria (VMCC and VMDP) for describing solid / intermediate states, specifically defined by equations (h) to (i) below:

[0113] y vmdp =C f tr(τ)+||s||-C c (h);

[0114]

[0115] Among them, y vmdp yvmcc These are the forms obtained by combining the yield criteria of von Mises with those of Drucker Prager and Cam Clay, respectively; C f C c These represent the friction angle and degree of polymerization of the material, respectively; τ and s = dev(τ) are the Kirchhoff stress and shear stress obtained from the deformation gradient using the elastic constitutive model, i.e. When y < 0, the material point violates the yield criterion, and a back-mapping algorithm is needed to correct it using the backward Euler form expressed by the following equation (j):

[0116] b n+1 -b tr =-2δγG(τ) n+1 )b n+1 (j).

[0117] in, It is a non-associative return mapping rule with volume-preserving properties; δγ=γΔt is the plastic flow distance; b=FF T Let be the left Cauchygreen tensor. To realize non-Newtonian flow in the intermediate phase, based on the Herschel-Bulkley model, we use the viscosity coefficient η and shear rate coefficient h during the return mapping process to control the non-Newtonian behavior of the material points. When h < 1, it exhibits shear thinning; when h > 1, it exhibits shear thickening. The specific formula is defined by the following equation (k):

[0118]

[0119] in, d represents the number of dimensions; y vm This is the corresponding yield criterion function. s is solved using the bisection method. n+1 Then, the corrected F can be finally obtained through the following equation (l). n+1 That is, the deformation gradient of the particle in the next time step:

[0120]

[0121] This completes the iteration of the current time step.

[0122] The specific process of discretization and evolution of the phase-field model in this invention is as follows:

[0123] Introducing the phase field as an additional field into the material point method's calculation process, the specific formula for increasing the transmission of the phase field between particles and the mesh in the above steps is expressed by the following equation (m):

[0124]

[0125] Where c is the phase field value carried on each particle, which evolves synchronously with the update of the matter point method as an additional scalar field; (subscript i represents the attribute carried by the mesh, and subscript p represents the particle attribute).

[0126] The weak form of the Allen-Cahn (AC) and Cahn-Hilliard (CH) equations is derived, ultimately manifested in the addition of a phase field solution step during mesh update to update the phase field value for the next time step. The specific matrix form of the solution is defined by the following equations (n) to (o):

[0127] For the CH phase field:

[0128] For phase field AC:

[0129] Where D is the diffusion rate of the phase field, and ∈ represents the width of the phase field interface. n+1 Let be the column vector formed by the phase field values ​​of all grid points. Let be the column vector of the free energies of all grid points; in addition, three auxiliary matrices are constructed for convenient calculation, M, r, and L being the mass moment, potential moment, and discretization of the Laplace operator in the form of the matter point method, respectively, as defined by the following equations (p) to (r):

[0130]

[0131]

[0132]

[0133] H is a user-defined phase field potential energy function. By setting the potential energy function, the phase field can be made to evolve in a directional potential well. The meanings of the other parameters have been explained above.

[0134] Solving large linear equation systems using the conjugate gradient method with Jacobi preconditions typically achieves convergence within 10–20 iterations. Based on the different properties of the AC and CH phase fields, different phase fields can be selected according to the requirements of the simulation scenario. At this point, the phase field value c for the next time step can be obtained. n+1 Based on the updated phase field values, we can control the subsequent material properties.

[0135] The phase-field driven control process for the elastic-viscoplastic material properties in this invention is specifically as follows:

[0136] Based on the phase field value c carried by each particle, the range of phase field values ​​is divided into the solid phase: [0, c a ], intermediate phase: [c a ,c b], liquid phase [c b [,1], where 0 represents a completely dry solid, 1 represents a completely liquid, and c is set. m To represent the state where the mesophase has the highest viscosity; c a ,c b The parameters set by the user are used to control the phase range, and in experiments they are generally set to 0.3 and 0.7 respectively.

[0137] b) Smooth interpolation between e1 and e2 is performed using Hermite polynomials, the function S1 of which is defined by the following equation (s):

[0138]

[0139] Where e1 and e2 are the starting and ending points of the gradual change, respectively; x is the independent variable.

[0140] c) Define the smooth function S2, which is the opposite of S1, by the following equation (t):

[0141] S2(e1,e2,x)=S1(-e2,-e1,-x)(t).

[0142] d) For VMDP, the parameters as a function of the phase field are defined by the following equations (u) to (v):

[0143]

[0144]

[0145] Where, φ c C is the friction angle; c This is the degree of polymerization coefficient; These are the user-defined initial friction angle of the solid and the maximum polymerization coefficient during the phase transformation process, respectively, which represent the magnitude of the yield stress; c is the phase field value of the corresponding particle.

[0146] e) Set the polymerization coefficient and hardening coefficient according to the following equations (w) to (x):

[0147] C c =κsinh(ξmax(-α,0))+C c (w);

[0148] ξ=S2(0,c m ,c)ξ 0 (x);

[0149] Where, ξ 0 c is the maximum value of the hardness coefficient. m This is the location of maximum viscosity; as the phase field value c approaches c... m The material's hardening ability gradually decreases to 0 while its viscosity η gradually increases.

[0150] f) The change in viscosity during a phase transition as defined by the following equation (y):

[0151] η=min(S1(c a ,c m ,c),S2(c m ,c b ,c)) (y);

[0152] Where η is the dynamic viscosity coefficient during the phase transition process, and the definitions of other parameters can be explained above. By controlling the material properties through the phase field value, the VMDP and VMCC models will automatically degenerate into the von Mises model when the phase field value enters the intermediate phase stage. Combined with formula (11), we can achieve the corresponding viscoplastic effect in the intermediate phase, and the entire process from the solid phase to the intermediate phase is seamless.

[0153] The specific implementation process of this invention is as follows: Simulated objects in the scene are discretized into particles, and a background mesh covering the simulation area is created. Each particle is then assigned initial velocity, mass, volume, deformation gradient, phase field value, and other attributes. This invention uses a semi-implicit Euler method to update the time step and restricts the time step size through CFL conditions. The specific flow of the phase transition simulation method based on the matter point method and phase field model is shown in Algorithm 1 below:

[0154]

[0155] The material phase transition simulation method involved in this invention is illustrated below:

[0156] See Figure 1 The image illustrates a scenario where powder and water are mixed. Water is poured into the soluble powder, and the mixture is stirred with a stirring rod. When the stirring rod is lifted, the viscous intermediate product rises along with the powder, creating a stringy, viscoplastic effect. This scenario uses a VMDP model as the constitutive model for the powder and water mixture, and employs a CH phase field (potential well set to 0.5). It can be observed that upon contact between water and particles, a smooth, gradually changing contact surface is created due to the application of the phase field model. Furthermore, the material's behavior transitions from a liquid and solid state to an intermediate viscous state, consistent with the entire process of a contact phase transition. During the stirring rod's ascent, the stringy phenomenon caused by the viscosity of the intermediate substance mimics the viscoplastic deformation that occurs during the preparation of rice or noodles.

[0157] See Figure 2The figure further illustrates the simulation of complex multiphase, multimaterial contacts within the framework. An armadillo-shaped, fragile object falls onto a beach, breaks into numerous fragments, and eventually melts into a viscous intermediate mixture due to water erosion. This experiment demonstrates that, in handling complex multiphase contacts, the framework of this invention can provide simultaneous simulation of various materials, such as particulates, fragile bodies, and fluids, while simultaneously simulating the phase transition processes between them.

[0158] See Figure 3 The figure shows the phase separation effect produced by AC and CH phase field models under a random phase field condition. It can be observed that as the phase field evolves, the red and blue phases gradually separate and produce a clear and smooth interface, demonstrating the correctness of the discretization and solution of the phase field model in this invention. At the same time, the CH phase field can produce a smoother and more active separation phenomenon compared with the AC phase field.

[0159] See Figure 4 and Figure 5 The image shows a rabbit-shaped candy being placed into a funnel, where it breaks due to the fragility of the material. Upon heating and melting, the candy melts, and the heating process can be observed to propagate smoothly in a wave-like pattern. The molten, viscous candy then flows out of the funnel and cools upon contact with the ground (color changes indicate temperature changes). As it cools and its viscosity increases, the candy creates a swirling visual effect.

[0160] See Figure 6 The figure shows the visual effect of the present invention simulating the phase transition of an egg when heated. As the temperature rises, the volume of the egg gradually expands, demonstrating that the present invention can not only handle material properties, but also make very rich controls on other properties such as volume based on the phase field.

[0161] The above examples are merely specific embodiments of the present invention. Obviously, the present invention is not limited to the above embodiments, and many other phase transition simulation scenarios are possible. All phase transition simulation scenarios that can be directly derived or conceived by those skilled in the art from the content disclosed in this invention should be considered within the scope of protection of this invention.

Claims

1. A method of simulation of phase transformation of elastic-viscoplastic material based on the material point method, characterized in that, The simulation method specifically comprises the following steps: 1) based on the material point method, Drucker-Prager and Cam-Clay are combined with von Mises yield criterion respectively, a unified numerical model for describing elastic-viscoplastic material and non-Newtonian fluid phenomena is established, which is used for simulating the material in the phase transition process, the simulation of the unified elastic-viscoplastic material and non-Newtonian fluid based on the material point method specifically comprises: a) particle properties are transmitted to the background grid The simulation target object is discretized into particles carrying a certain mass m p and volume The particles, and a background grid G that can cover the simulation region, are reinitialized, i.e. all the carried values are set to zero, and the mass and momentum carried by the particles are transferred to the nearby background grid G according to the APIC form and the spline interpolation function w ip The specific formula is represented by the following (a1)~(a2) formula: where m n represents the quality at the n-th time step, the index i represents the attribute carried by the mesh, and the index p represents the particle attribute; represents the velocity at the n-th time step; is the local velocity affine used in the APIC; w ip is the spline interpolation function; b) the momentum of the background grid point is updated According to the deformation gradient of the particles around the grid point and the elastic constitutive of the system is solved by the stress of the following (b) formula The stress is then updated according to The velocity of the corresponding grid point is updated by the following (c) formula where Ψ is the energy density function defined by the elastic constitutive relation; V p 0 is the initial volume of the particle; is the deformation gradient of particle p; w ip is the spline interpolation function between grid point i and particle p; represents the velocity of the next time step grid point; Δt is the time step length; c) the background grid properties are transmitted back to the particle and the deformation gradient is updated c-1: the surrounding particles are transmitted back according to the updated physical properties of the background grid and the spline interpolation function, and the local velocity affine matrix is updated for APIC, and the specific formula is represented by the following (d)-(f) formulas: where, represents the velocity of the particle at the next time step; is the local velocity affine used in the APIC for the next time step; represents the velocity of the grid point at the next time step; is the position of the particle p at the next time step; Δt is the time step length; c-2: velocity of surrounding grid points deformation gradient is updated to obtain a deformation gradient F predicted by the following equation (g) tr : where F tr is the shape deformation gradient, which is a prediction of the shape deformation at the next time step for subsequent correction; is estimated by the velocity gradient at particle p. c-3: according to the viscoplastic constitutive model of the material, the deformation gradient that violates the yield criterion is corrected according to the return mapping algorithm, and the two yield criteria y of the solid / intermediate state are described by the following (h)-(i) formulas: y vmdp = C f tr(τ) + ||s|| - C c (h); where y vmdp , y vmcc are the forms of von Mises and Drucker Prager and Cam Clay yield criteria mixed respectively; C f , C c are the friction angle and the degree of aggregation of the material respectively; τ and s = dev(τ) are the Kirchhoff stress and the shear stress obtained from the elastic constitutive according to the deformation gradient, i.e. When y < 0, the material point violates the yield criterion, then the return mapping algorithm is needed to correct using the backward Euler form expressed by the following (j) formula: b n+1 -b tr = -2δγG(τ n+1 )b n+1 (j); wherein, is a non-associative return mapping rule with the property of volume preservation; δγ = γΔt is the plastic flow distance; b = FF T is the left Cauchy-Green tensor, i.e. b tr = F tr F tr,T ; c-4: according to the Herschel-Bulkley model, the viscosity coefficient η and the shear rate coefficient h are used to control the performance of the material point non-Newtonian behavior in the return mapping process, when h<1, it shows shear thinning, and when h>1, it shows shear thickening, and the specific formula is represented by the following (k) formula: wherein, d is the number of dimensions; y vm is the corresponding yield criterion function; c-5: Solve s using bisection n+1 Then, the modified F can be finally obtained by the following equation (1) n+1 i.e. the next time step particle p deformation gradient: 2) the Allen-Cahn and Cahn-Hilliard phase field equations are discretely modeled and introduced into the material point method, and the phase transition caused by the contact and external influence of the material is simulated through the evolution of the phase field, specifically comprising: a) the phase field is introduced into the calculation process of the material point method as an additional field, and the transmission of the phase field between the particle and the grid is increased in the above step a) particle properties are transmitted to the background grid, and the specific formula is represented by the following (m) formula: where c n is the phase field value at the nth time step, with subscript i representing the property carried by the grid, and subscript p representing the particle property; b) according to the weak form derivation of the Allen-Cahn (AC) and Cahn-Hilliard (CH) equations, the phase field solution is added in the grid update in the above step b) the momentum of the background grid point is updated, which is used to update the phase field value of the next time step, and the matrix form of the specific solution of the CH phase field is represented by the following (n) formula: The matrix form of the specific solution of the AC phase field is represented by the following (o) formula: where D is the diffusion rate of the phase field; ∈ is the width coefficient of the phase field interface; M, r, L are the mass moment, potential energy moment and Laplace operator of the phase field respectively; c n+1 is a column vector composed of the phase field values of all grid points; is a column vector composed of the free energies of all grid points; The mass moment M, potential energy moment r and Laplacian operator of the phase field are discretely L in the form of the material point method, which are represented by the following (p)-(r) formulas: Wherein, H is a user-defined phase field potential energy function; 3) the attribute dynamic control strategy of the elastic-viscoplastic material driven by the phase field evolution is introduced, and the discontinuity in the performance caused by the physical property change process of the mixed contact surface is eliminated, specifically comprising: a) According to the phase field value c carried by each particle, the range of phase field value is divided into solid phase: [0, c a ], mesophase: [c a , c b ], liquid phase [c b , 1], wherein 0 is a completely dry solid, 1 is a complete liquid, and c m is set to represent the state of maximum mesophase viscosity; c a , c b are parameters set by the user to control the phase range, c a , c b are set to 0.3 and 0.7, respectively; b) Hermite polynomials are used for smooth interpolation between e1 and e2, and the function S1 is defined by the following (s) formula: Wherein, e1 and e2 are the starting point and the end point of the gradient respectively; x is the independent variable; c) the smooth function S2 opposite to S1 is defined by the following (t) formula: S2(e1,e2,x)=S1(-e2,-e1,-x) (t); d) For VMDP, the parameter variation with phase field is defined by the following (u)~(v) equations: where φ c is the friction angle; C c is the degree of polymerization coefficient; are the user-defined initial friction angle of the solid and the maximum degree of polymerization during the phase transformation, respectively, i.e., representing the size of the yield stress; c is the phase field value corresponding to the particle. e) The aggregation coefficient and the hardening coefficient are set by the following (w)~(x) equations, respectively: C c = K sinh(ξ max(-a, 0)) + C c (w); ξ = S2(0, c m ,c)ξ 0 (x) where ξ 0 is the maximum value of the hardness coefficient; c m is the maximum viscosity position, as the phase field value c m approaches c , the hardening ability of the material gradually decreases to 0 and the viscosity η gradually increases; f) The viscosity variation in the phase transition process is defined by the following (y) equation: η = min(Si(c a ,c m ,c), S2(c m ,c b ,c)) (y); wherein η is the dynamic viscosity coefficient in the phase transition process.

Images