CN111104753B - Viscous incompressible fluid simulation method based on SPH - Google Patents
Viscous incompressible fluid simulation method based on SPH Download PDFInfo
- Publication number
- CN111104753B CN111104753B CN201911390333.9A CN201911390333A CN111104753B CN 111104753 B CN111104753 B CN 111104753B CN 201911390333 A CN201911390333 A CN 201911390333A CN 111104753 B CN111104753 B CN 111104753B
- Authority
- CN
- China
- Prior art keywords
- fluid
- particle
- particles
- equation
- pressure
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 239000012530 fluid Substances 0.000 title claims abstract description 297
- 238000000034 method Methods 0.000 title claims abstract description 77
- 238000004088 simulation Methods 0.000 title claims abstract description 30
- 230000008569 process Effects 0.000 claims abstract description 23
- 239000002245 particle Substances 0.000 claims description 267
- 238000004364 calculation method Methods 0.000 claims description 34
- 239000011159 matrix material Substances 0.000 claims description 17
- 230000009471 action Effects 0.000 claims description 16
- 238000012937 correction Methods 0.000 claims description 5
- 239000006185 dispersion Substances 0.000 claims description 4
- 230000005484 gravity Effects 0.000 claims description 3
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 230000000694 effects Effects 0.000 abstract description 3
- 238000004804 winding Methods 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 7
- 239000000463 material Substances 0.000 description 5
- 230000008859 change Effects 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 239000010408 film Substances 0.000 description 2
- 235000012907 honey Nutrition 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 description 1
- 235000007688 Lycopersicon esculentum Nutrition 0.000 description 1
- 240000003768 Solanum lycopersicum Species 0.000 description 1
- 239000010426 asphalt Substances 0.000 description 1
- 230000001808 coupling effect Effects 0.000 description 1
- 238000012217 deletion Methods 0.000 description 1
- 230000037430 deletion Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000004519 grease Substances 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 235000015067 sauces Nutrition 0.000 description 1
- 238000010008 shearing Methods 0.000 description 1
- 230000008719 thickening Effects 0.000 description 1
- 239000010409 thin film Substances 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention relates to a viscous incompressible fluid simulation method based on SPH, which is a viscous incompressible fluid simulation method based on SPH with high precision and high robustness. The invention can stably simulate Newtonian fluid with different viscosity coefficient values and non-Newtonian fluid with variable viscosity coefficient, can keep the surface morphology of the viscous fluid for a long time, and can accurately simulate a plurality of fine motion processes such as the stretching and thinning process, the rope winding effect, the film morphology and the like of the viscous fluid.
Description
Technical Field
The invention belongs to the field of computer graphics, and particularly relates to a viscous incompressible fluid simulation method based on SPH.
Background
In production and life of people, viscous fluid is seen everywhere, such as petroleum, asphalt, honey, tomato sauce and the like, and therefore the viscous fluid is an important simulation object in the field of graphics. Although the meshless particle method/smooth particle fluid dynamics (SPH) method is a well-established and widely used fluid simulation method, the system method still has difficulty in achieving stable, fast, and accurate simulation of viscous fluids. Particularly, the traditional meshless particle method/smooth particle fluid dynamics method performs poorly when simulating the draw down characteristics of viscous fluids, simulating viscous fluids in the form of thin films, maintaining the details of the surface morphology of highly viscous fluids, and simulating non-newtonian fluid materials with widely varying viscosity coefficients.
The reason why the conventional meshless particle method/smooth particle fluid dynamics (SPH) method is poor in simulating viscous fluid is mainly as follows. (1) In the traditional grid-free particle method, in order to ensure the convenience of a numerical calculation process, the solution of the pressure force (pressure generated by density change in a fluid material) and the viscous force (shearing force) in fluid particles is divided into two independent calculation steps, so the solution result of the viscous force and the solution result of the pressure force are often in conflict with each other, the solution results of the viscous force and the solution results of the pressure force are mutually damaged, and finally, a plurality of complex and fine motion characteristics of viscous fluid are lost; (2) when a meshless particle method is used for simulating fluid, due to the fact that neighborhood particles of particles near the surface of the fluid are lost, physical quantity calculation errors near the surface of the fluid can be caused, and simulation quality is further influenced; (3) because the stretching instability problem exists in the gridless particle method, in order to ensure the stability (robustness) of the simulation process, the conventional gridless particle method often eliminates the situation that the pressure is less than 0, and therefore the stretching state of the viscous fluid cannot be accurately simulated.
Disclosure of Invention
The invention solves the problems: aiming at the problem that the traditional grid-free particle method is difficult to accurately and stably simulate the motion of viscous fluid, a viscous incompressible fluid simulation method based on SPH with high precision and high robustness is provided, a unified calculation scheme of pressure/viscous force is provided, the problem of neighborhood loss of fluid surface particles is solved by introducing a boundary condition without constructing a gshost particle, negative pressure is kept when a pressure Poisson equation is calculated, and finally, the rapid, accurate and high-robustness simulation of the viscous fluid can be realized.
The invention provides a high-precision and high-robustness viscous fluid rapid simulation method aiming at the problem that the traditional non-grid particle method is difficult to accurately and stably simulate the motion of viscous fluid. According to the method, the non-physical conflict and drift error generated in the solving process are avoided through the unified solving scheme of incompressible constraint and viscous constraint, the problem of surface particle loss is rapidly and efficiently solved in a mode of not distributing Ghost particles, meanwhile, the influence of negative pressure at partial positions of the fluid in the stretching process is considered, and stable, accurate and rapid simulation of Newtonian fluid and non-Newtonian fluid with different viscosity values is finally achieved. The invention can keep the surface shape of the viscous fluid for a long time; the method can simulate a plurality of fine motion processes such as the stretching and thinning process, the rope winding effect, the film form and the like of the viscous fluid; and the simulation of the mutual coupling effect of different viscous fluids in the same simulation scene can be realized.
The specific technical scheme of the invention is as follows: an SPH-based viscous incompressible fluid simulation method, as shown in FIG. 2, comprising the steps of:
step 1: based on the Navistokes equation and the condition of no velocity divergence, the discretization of fluid particles is realized by adopting smooth particle fluid dynamics (SPH) to obtain a dynamic equation;
step 2: updating a velocity field and a position field of the fluid particles under the action of mass force according to the kinetic equation obtained in the step 1;
and step 3: performing integral iteration of incompressibility and viscosity constraint;
and 4, step 4: establishing a Poisson equation of the pressure of each fluid particle according to the kinetic equation obtained in the step 1, discretizing the Poisson equation of the pressure and converting the pressure into a linear system, obtaining the pressure value of each fluid particle by solving the linear system, and correcting corresponding elements in a coefficient matrix of the Poisson equation linear system of the pressure of the fluid particles near a boundary according to the fluid particles which are not lost in the neighborhood inside the fluid when a part of the fluid particles are positioned near the fluid boundary;
and 5: calculating the pressure gradient of each fluid particle based on the pressure values of the fluid particles acquired in the step 4, and correcting the calculation result of the pressure gradient of the fluid particles near the boundary according to the fluid particles which are not lost and are located in the neighborhood inside the fluid when the part of the fluid particles are located near the fluid boundary;
step 6: constructing a velocity viscosity constraint implicit differential equation containing the pressure gradient action according to the pressure gradient of each fluid particle in the step 5, discretizing and converting the velocity viscosity constraint implicit equation containing the pressure gradient action into a linear system, and solving a fluid particle velocity field after applying the viscosity constraint;
and 7: after step 6, finishing 1 iteration, and accumulating the iteration times;
and 8: if the iteration times are less than the set value or the iteration is not converged, repeating the operations from the step (4) to the step (7);
and step 9: taking the fluid particle velocity field obtained by calculation after the iteration is finished as the velocity field of the fluid particles at the current moment, and updating the position field of each fluid particle based on the velocity field;
step 10: updating the simulation time, and repeating the operation between the steps 2-9 until the termination time or the set simulation termination condition is reached, and finishing the whole simulation process.
In step 4, the process of the dispersion and conversion of the poisson equation of the pressure into the linear system is as follows:
the left term of the pressure poisson equation, i.e., the left end discretization/equation form of the equation, is as follows:
if the fluid particles corresponding to a certain row of the coefficient matrix are located in the fluid, the left term of the equation is:
If the fluid particles corresponding to a certain row of the coefficient matrix are positioned near the surface of the fluid, the left term of the equation is as follows:
coefficient of fluid particles near surface A0The calculating method of (2):
each fluid particle needs to be calculated as Ai＝∑jaijTaking A of the particles which are filled in the neighborhood inside the fluidiAs A0(ii) a Or in the initial stage of simulation, each particle calculates Ai＝∑jaijFinally A0The value of (A) is as follows: a. the0＝max(Ai)；
In the above formulas, i represents the current fluid particle to be solved, j represents the neighborhood fluid particle of the fluid particle i, and pi、piRespectively representing the pressure value of the fluid particle i to be solved and the pressure value of the neighborhood particle j of the particle i, aijTo pressDiscretized strong Poisson equation coefficient, VjIs the volume of the fluid particle j, rijIs the relative distance, w, between the fluid particle i and the fluid particle jijAnd eta is a kernel function in the SPH method and is a constant with a small value.
In step 4, the poisson equation form of the pressure is as follows:
the poisson equation of the pressure is specifically in the form:
secondly, the source item discretization scheme of the Poisson equation of the pressure intensity is as follows:
third, when the Poisson equation of each pressure is calculated, the left term of the equation
Density correction is needed to compensate the density loss in the fluid motion process, and the compensation amount is calculated in a specific form:
in the above formulas, i represents the current fluid particle to be solved, j represents the neighborhood fluid particle of the fluid particle i, and piFor the pressure value, v, of the fluid particle to be solvediIs the velocity field, v, of the fluid particle iijRefers to the relative velocity field, V, between particle i and particle jjVolume of particle j, wijIs a kernel function in SPH, alpha is a constant coefficient, rho0Is the initial density of the fluid, piΔ t is the simulated time step for the current density of the fluid particle i.
If the fluid particles corresponding to a row of the coefficient matrix are located near the fluid surface, the method for determining whether the fluid particles are located near the fluid surface comprises: each particle stores a surface judgment mark AiWherein the calculation formula is Ai＝∑jaijIf A isi<A0If the particle is located near the surface of the fluid, the particle is considered to be located at the surface of the fluid, otherwise, the particle is located inside the fluid;
wherein a isijThe calculation formula of (A) is as follows:
i represents the fluid particle currently being solved for, j represents the fluid particle in the neighborhood of the fluid particle i, aijFor the discretized coefficient of the pressure Poisson equation, VjIs the volume of the fluid particle j, rijIs the relative distance, w, between the fluid particle i and the fluid particle jijFor the kernel function in the SPH method, η is a constant whose value is less than one percent of the distance between adjacent fluid particles.
In step 5, when a part of the fluid particles is located near the fluid boundary, the pressure gradient calculation result of the fluid particles near the boundary needs to be corrected, and the specific correction method is as follows:
when the fluid particles are inside the fluid:
when the fluid particles are near the surface of the fluid:
in the above formulas, i represents the current fluid particle to be solved, j represents the neighborhood fluid particle of the fluid particle i,
is the pressure value of the fluid particles i,
is the pressure value, V, of the fluid particle jjIs the volume of particle j, rijIs the relative distance, w, between the fluid particle i and the fluid particle jijIs the kernel function in SPH.
In the step 6, a velocity viscosity constraint implicit differential equation containing a pressure gradient action is constructed, and a velocity viscosity constraint implicit square differential equation containing the pressure gradient action is discretized and converted into a linear system in the specific form:
the velocity viscosity constraint implicit differential equation right term of the pressure gradient action is as follows, and the discretization scheme of the second order differential of the velocity field and the viscosity coefficient is as follows:
in the above
The unknown quantity to be solved is the velocity viscosity constraint implicit differential equation under the action of the pressure gradient, namely the fluid particle velocity field after the viscosity constraint is applied;
i represents the fluid particle currently being solved for, j represents the fluid particle in the neighborhood of the fluid particle i,
is the pressure value of the fluid particles i,
to apply the velocity field of the fluid particles i before viscous confinement,
to apply a velocity field of the fluid particles i after viscous confinement,
means a relative velocity field, V, between particles i and j after application of a viscous constraintjIs the volume of particle j, rijIs the relative distance, x, between fluid particle i and fluid particle jijIs the relative position field between particle i and particle j, wijAs kernel functions in SPH，μiAnd mujViscosity coefficients, ρ, for a fluid particle i and a neighboring fluid particle j, respectively0At is the initial density of the fluid, Δ t is the simulated time step.
In the step 1, the dispersion of the fluid particles is realized by adopting a Smooth Particle Hydrodynamics (SPH) under the navistokes equation and the condition of no velocity divergence, and the specific form is as follows:
the NaviStokes equation:
no velocity divergence condition:
the physical field discretization calculation scheme is as follows: phi is ai＝-∑jVjφjwij；
The first order differential calculation scheme of the physical field is as follows:
in the above formula, v is a velocity field, p is a pressure, τ is a viscous stress, ρ is a density, f is a volume force (gravity), μ is a viscosity coefficient of the fluid, i represents a fluid particle to be solved currently, j represents a neighborhood fluid particle of the fluid particle i, and φj、φiAre any physical quantities, V, of particles j and i, respectivelyjVolume of particle j, wijIs the kernel function in SPH, phiijIs the difference in physical quantity between particle i and particle j.
Further, the convergence characterization parameter of the ensemble iteration of step (12) may be a pressure value or a fluid particle velocity field in the poisson equation and other approximate forms. And (5) when the difference value between the convergence characterization parameters of two continuous iterations is less than the convergence condition in the step (8), terminating the iteration and carrying out the next calculation.
Further, all the fluid particles need to be processed from step (5) to step (13) synchronously, so that the whole processing process is completed.
Compared with the prior art, the method has the following advantages:
(1) the physical significance of the Poisson equation, the calculation mode of viscosity constraint and the pressure-viscosity unified iteration scheme provided by the invention is clear and is closer to the action process of real viscous fluid, so that the surface details of the viscous fluid can be more effectively maintained;
(2) the method for processing the particles near the surface of the fluid can effectively avoid calculation errors and errors caused by the neighborhood loss of the particles on the surface, so that the method can better simulate the process of stretching and thinning the viscous fluid;
(3) the overall scheme provided by the invention has very high stability, and can simulate the range from 0 Pa.s to very high, namely 8X 108A Newtonian fluid with a viscosity coefficient of Pa · s can also simulate a shear-thinning and shear-thickening fluid with a viscosity coefficient varying over a wide range.
(4) When the surface tension effect is introduced in the implementation steps of the method, high-viscosity and high-surface-tension fluids such as honey, grease and the like can be simulated more truly, and the film form of the fluids can be simulated stably;
(5) each step in the invention has high parallelism, thus being suitable for the prior GPU architecture.
Drawings
FIG. 1 is a unified solver diagram of incompressible and viscous constraints for viscous fluids;
FIG. 2 is a flow chart of an implementation of the present invention;
fig. 3 is a schematic view of the position of the fluid particles.
Detailed Description
In order to make the aforementioned objects, features and advantages of the present invention comprehensible, embodiments accompanied with figures are described in detail below, without limiting the present invention thereto.
1. The hardware platform of the method adopts a six-core CPU with the model number of Intel i7-8700K, the main frequency of 3.7GHz, an NVIDIA GeForce GTX 2080 display card and a video memory of 8 GB. The system program is written in C + +, a CUDA language is adopted for accelerating the parallel computing part, the program is compiled and executed by means of Microsoft Visual Studio 2015, and open source libraries such as OpenGL and Freeplus are adopted in the development process.
2. The method takes a NaviStokes equation (momentum equation) and a non-velocity-divergence and non-compressible condition as the dynamic basis of the simulation process, and the specific form is as follows.
The navistokes equation/momentum equation:
where f is the mass force (vector) experienced by the fluid particle, ρ is the fluid density value, p is the pressure field (scalar) of the fluid particle, and τ is the viscous force field (vector) within the fluid.
Calculation formula of viscous force field:
wherein
The velocity field gradient of the fluid particles (third order matrix in 3 dimensions and second order matrix in 2 dimensions) is the viscosity coefficient (scalar).
Non-velocity divergence incompressible condition:
where v is the fluid mass point velocity field.
In the above formula, v is a velocity field, p is a pressure, τ is a viscous stress, ρ is a density, f is a volume force (gravity), and μ is a viscosity coefficient.
3. The invention completes the discretization of various physical fields based on the fluid dynamics of smooth particles.
The general physical field discretization formula is: phi is ai＝-∑jVjφjwijWherein phiiIs a physical field of some kind (e.g. density field, etc.), the j fluid particles are the neighborhood particles of the i fluid particles, wijIs a kernel function whose value is related to the distance between the i fluid particles and the j fluid particles;
the first order differential discretization formula is:
wherein phiiA certain type of physical field (such as a speed field, a pressure field and the like);
the discretization formula of the second-order velocity field is as follows:
the second-order pressure field discretization formula is as follows:
middle phi of the above formulaj、φiThe physical quantities are arbitrary types of the particles j and i.
4. Fig. 1 is a schematic diagram of a unified solver for incompressible constraint and viscous constraint of viscous fluid, which is implemented by calculating a poisson equation and an implicit viscous constraint equation and performing overall iteration until an obtained pressure field or velocity field converges, so as to finally obtain a fluid particle velocity field satisfying both an incompressible condition and a viscous constraint condition, thereby avoiding a conflict between incompressible calculation and viscous calculation in a conventional SPH simulation method.
In order to avoid reduction of calculation accuracy caused by drift errors generated during iterative calculation of incompressible constraint and viscous constraint, a similar discretization scheme is adopted for a Poisson equation and an implicit viscous constraint equation.
In order to realize the simulation of the stretching state of the viscous fluid, the negative pressure generated in the process of solving the poisson equation needs to be reserved.
In order to ensure the calculation accuracy of the high-order differential physical quantity of the particles near the fluid boundary, the problem of neighborhood deletion of the particles near the boundary needs to be solved. As shown in fig. 3, the particle a is a particle inside the fluid, and there are enough particles in the vicinity of the fluid particle, so that the accuracy of calculating each physical quantity of the particle a can be ensured, but the fluid particle located near the fluid boundary, such as the particle B in fig. 3, has fewer adjacent particles, so that the accuracy of calculating the physical quantity of the fluid particle near the boundary is seriously affected. The conventional SPH method generates the gshost particles outside the free surface of the fluid, but it takes up a lot of memory space and causes extra computation time consumption. The present invention assumes that there are ghost particles, but no ghost particles are generated. Instead, when assembling the main diagonal elements of the coefficient matrix of the poisson equation linear system, the main diagonal elements of the coefficient matrix corresponding to the boundary particles are directly set to be uniform values, and the values are the same as the main diagonal elements of the coefficient matrix corresponding to the particles with the neighborhoods arranged inside the fluid. By the method, the generation of the ghost particles can be avoided, and the accuracy of second order differential calculation of the boundary particle pressure field can be ensured. The specific calculation formula is as follows.
When the fluid particles are located at a position near the surface of the fluid:
when the fluid particles are located at a position inside the fluid:
where p is the pressure field (scalar) of the fluid particles.
After the pressure Poisson equation is converted into a linear system, the main diagonal element of the coefficient matrix of the fluid particle i pressure is AiThe calculation method comprises the following steps: the method specifically comprises the following steps: a. thei＝∑jaij(ii) a Wherein a isijFor the j coefficient of the ith equation, the calculation formula is:
when the fluid particles i are surface particles, the main diagonal coefficient of the fluid particles i is A0A, which is numerically equal to the diagonal elements of the coefficient matrix for the fluid particle i when the particle neighborhood within the fluid is fulli。
Calculating the pressure field gradient of each fluid particle is equivalent to calculating the pressure to which the fluid particle is subjected. The particles on the surface of the fluid should take the action of air into consideration, and in order to ensure the accuracy of the gradient of the particle pressure field near the surface of the fluid, the traditional SPH method needs to be realized by arranging the ghost particles outside the free surface. The method avoids extra space consumption and time consumption caused by setting the ghost particles by the following calculation.
When the fluid particles are inside the fluid:
when the fluid particles are near the surface of the fluid:
in the viscous constraint solving process of the velocity, correction of fluid particles near the boundary does not need to be considered, and therefore, in the overall iteration process, the pressure solving (poisson equation/incompressible constraint) needs to be ensured to be before the viscous constraint, otherwise, more serious calculation errors can be caused.
The implicit viscous constraint equation for velocity is in the form:
wherein
The discretization scheme of the term is:
in the equation
Is the fluid particle velocity field prior to application of the viscous constraint,
the fluid particle velocity field after applying the viscosity constraint is also the unknown to be solved for this equation.
The Poisson equation and the implicit viscosity equation can be converted into a linear equation set, the coefficient matrix order of the linear equation set is the same as the number of fluid particles, and the linear equation set can be converged under GC and Jacobi method iterative solution.
5. Whether the fluid particles are located near the free surface of the fluid can be calculated by AiJudging (namely the main diagonal elements of the coefficient matrix in the linear system of the Poisson equation corresponding to the fluid particles i), when A isiIs less than A0When the fluid particles i are considered to be located near the free surface of the fluid.
6. Simulation of shear-thinning or shear-thickening fluids was achieved using the following physical model.
The viscosity coefficient of the fluid particles i is calculated as:
wherein,
is a scalar quantity related to the shear rate of the fluid particle:
when the value of m is positive, the fluid material is the shear thickening fluid; when m is negative, the fluid material is shear thinning fluid.
In the above formulas, DiIs the shear deformation tensor of the fluid particle i,
for characterizing the extent of material deformation at the i-position of the fluid particle, μ0、μ∞Respectively, a lower limit value and an upper limit value of the change of the coefficient of viscosity of the fluid simulated at present, k is a constant coefficient larger than zero, and m is an integer which can be negative.
7. In the invention, after the uniform solution of the incompressible constraint and the viscous constraint of the viscous fluid is completed, the surface tension constraint is applied to the fluid surface particles, and the viscous fluid with large surface tension can be simulated.
8. The physical model of the incompressible condition is true for the no velocity divergence condition (fluid particle velocity divergence is zero), i.e.:
The above embodiments are merely illustrative and not restrictive, and those skilled in the art may modify the technical solution of the present invention without departing from the spirit and scope of the present invention, which is defined by the claims.
Claims (7)
1. An SPH-based viscous incompressible fluid simulation method comprising the steps of:
step 1: based on the Navistokes equation and the condition of no velocity divergence, the discretization of fluid particles is realized by adopting smooth particle fluid dynamics (SPH) to obtain a dynamic equation;
step 2: updating a velocity field and a position field of the fluid particles under the action of mass force according to the kinetic equation obtained in the step 1;
and step 3: performing integral iteration of incompressibility and viscosity constraint;
and 4, step 4: establishing a Poisson equation of the pressure of each fluid particle according to the kinetic equation obtained in the step 1, discretizing the Poisson equation of the pressure and converting the pressure into a linear system, obtaining the pressure value of each fluid particle by solving the linear system, and correcting corresponding elements in a coefficient matrix of the Poisson equation linear system of the pressure of the fluid particles near a boundary according to the fluid particles which are not lost in the neighborhood inside the fluid when a part of the fluid particles are positioned near the fluid boundary;
and 5: calculating the pressure gradient of each fluid particle based on the pressure values of the fluid particles acquired in the step 4, and correcting the calculation result of the pressure gradient of the fluid particles near the boundary according to the fluid particles which are not lost and are located in the neighborhood inside the fluid when the part of the fluid particles are located near the fluid boundary;
step 6: constructing a velocity viscosity constraint implicit differential equation containing the pressure gradient action according to the pressure gradient of each fluid particle in the step 5, discretizing and converting the velocity viscosity constraint implicit equation containing the pressure gradient action into a linear system, and solving a fluid particle velocity field after applying the viscosity constraint;
and 7: after step 6, finishing 1 iteration, and accumulating the iteration count;
and 8: if the iteration times are less than the set value or the iteration is not converged, repeating the operations from the step (4) to the step (7);
and step 9: taking the fluid particle velocity field obtained by calculation after the iteration is finished as the velocity field of the fluid particles at the current moment, and updating the position field of each fluid particle based on the velocity field;
step 10: updating the simulation time, and repeating the operation between the steps 2-9 until the termination time or the set simulation termination condition is reached, and finishing the whole simulation process.
2. The method of claim 1, wherein: in step 4, the process of the dispersion and conversion of the poisson equation of the pressure into the linear system is as follows:
the left term of the pressure poisson equation, i.e., the left end discretization/equation form of the equation, is as follows:
if the fluid particles corresponding to a certain row of the coefficient matrix are located in the fluid, the left term of the equation is:
If the fluid particles corresponding to a certain row of the coefficient matrix are positioned near the surface of the fluid, the left term of the equation is as follows:
coefficient of fluid particles near surface A0The calculating method of (2):
each fluid particle needs to be calculated as Ai＝∑jaijTaking A of the particles which are filled in the neighborhood inside the fluidiAs A0(ii) a Or in the initial stage of simulation, each particle calculates Ai＝∑jaijFinally A0The values of (A) are as follows: a. the0＝max(Ai)；
In the above formulas, i represents the current fluid particle to be solved, j represents the neighborhood fluid particle of the fluid particle i, and pi、pjRespectively representing the pressure value of the fluid particle i to be solved and the pressure value of the neighborhood particle j of the particle i, aijFor the discretized coefficient of the pressure Poisson equation, VjIs the volume of the fluid particle j, rijIs the relative distance, w, between the fluid particle i and the fluid particle jijFor the kernel function in the SPH method, η is a constant whose value is less than one percent of the distance between the majority of adjacent fluid particles.
3. The method of claim 1, wherein: in step 4, the poisson equation form of the pressure is as follows:
the poisson equation of the pressure is specifically in the form:
secondly, the source item discretization scheme of the Poisson equation of the pressure intensity is as follows:
third, during each calculation of the Poisson equation of the pressure intensityLeft term of equation
Density correction is needed to compensate the density loss in the fluid motion process, and the compensation amount is calculated in a specific form:
in the above formulas, i represents the current fluid particle to be solved, j represents the neighborhood fluid particle of the fluid particle i, and piFor the pressure value, v, of the fluid particle to be solvediIs the velocity field, v, of the fluid particle iijRefers to the relative velocity field, V, between particle i and particle jjVolume of particle j, wijIs a kernel function in SPH, alpha is a constant coefficient, rho0Is the initial density of the fluid, piΔ t is the simulated time step for the current density of the fluid particle i.
4. The method of claim 2, wherein: if the fluid particles corresponding to a row of the coefficient matrix are located near the fluid surface, the method for determining whether the fluid particles are located near the fluid surface comprises: each particle stores a surface judgment mark AiWherein the calculation formula is Ai＝∑jaijIf A isi<A0If the particle is located near the surface of the fluid, the particle is considered to be located at the surface of the fluid, otherwise, the particle is located inside the fluid;
wherein a isijThe calculation formula of (A) is as follows:
i represents the fluid particle currently being solved for, j represents the fluid particle in the neighborhood of the fluid particle i, aijFor the discretized coefficient of the pressure Poisson equation, VjIs the volume of the fluid particle j, rijIs the relative distance, w, between the fluid particle i and the fluid particle jijFor the kernel function in the SPH method, η is a constant whose value is less than one percent of the distance between adjacent fluid particles。
5. The method of claim 1, further comprising: in step 5, when a part of the fluid particles is located near the fluid boundary, the pressure gradient calculation result of the fluid particles near the boundary needs to be corrected, and the specific correction method is as follows:
when the fluid particles are inside the fluid:
when the fluid particles are near the surface of the fluid:
in the above formulas, i represents the current fluid particle to be solved, j represents the neighborhood fluid particle of the fluid particle i,
is the pressure value of the fluid particles i,
is the pressure value, V, of the fluid particle jjIs the volume of particle j, rijIs the relative distance, w, between the fluid particle i and the fluid particle jijIs the kernel function in SPH.
6. The method of claim 1, wherein: in the step 6, a velocity viscosity constraint implicit differential equation containing a pressure gradient action is constructed, and a velocity viscosity constraint implicit square differential equation containing the pressure gradient action is discretized and converted into a linear system in the specific form:
the velocity viscosity constraint implicit differential equation right term of the pressure gradient action is as follows, and the discretization scheme of the second order differential of the velocity field and the viscosity coefficient is as follows:
in the above
The unknown quantity to be solved is the velocity viscosity constraint implicit differential equation under the action of the pressure gradient, namely the fluid particle velocity field after the viscosity constraint is applied;
i represents the fluid particle currently being solved for, j represents the fluid particle in the neighborhood of the fluid particle i,
is the pressure value of the fluid particles i,
to apply the velocity field of the fluid particles i before viscous confinement,
to apply a velocity field of the fluid particles i after viscous confinement,
means a relative velocity field, V, between particles i and j after application of a viscous constraintjIs the volume of particle j, rijIs the relative distance, x, between fluid particle i and fluid particle jijIs the relative position field between particle i and particle j, wijIs a kernel function in SPH, muiAnd mujViscosity coefficients, ρ, for a fluid particle i and a neighboring fluid particle j, respectively0At is the initial density of the fluid, Δ t is the simulated time step.
7. The method of claim 1, wherein: in the step 1, the dispersion of the fluid particles is realized by adopting a Smooth Particle Hydrodynamics (SPH) under the navistokes equation and the condition of no velocity divergence, and the specific form is as follows:
the NaviStokes equation:
no velocity divergence condition:
the physical field discretization calculation scheme is as follows: phi is ai＝-∑jVjφjwij；
The first order differential calculation scheme of the physical field is as follows:
in the above formula, v is a velocity field, p is a pressure, τ is a viscous stress, ρ is a density, f is a volume force (gravity), μ is a viscosity coefficient of the fluid, i represents a fluid particle to be solved currently, j represents a neighborhood fluid particle of the fluid particle i, and φj、φiAny physical quantities, V, of particles j and i, respectivelyjVolume of particle j, wijIs the kernel function in SPH, phiijIs the difference in physical quantity between particle i and particle j.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911390333.9A CN111104753B (en) | 2019-12-30 | 2019-12-30 | Viscous incompressible fluid simulation method based on SPH |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911390333.9A CN111104753B (en) | 2019-12-30 | 2019-12-30 | Viscous incompressible fluid simulation method based on SPH |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111104753A CN111104753A (en) | 2020-05-05 |
| CN111104753B true CN111104753B (en) | 2022-03-25 |
Family
ID=70424046
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911390333.9A Active CN111104753B (en) | 2019-12-30 | 2019-12-30 | Viscous incompressible fluid simulation method based on SPH |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111104753B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113158531B (en) * | 2021-02-07 | 2022-06-21 | 南开大学 | Single-component and multi-component incompressible fluid simulation method utilizing deformation gradient |
| CN115248989A (en) * | 2021-04-25 | 2022-10-28 | 北京航空航天大学 | Viscous fluid simulation method based on yield criterion constraint |
| CN115270651B (en) * | 2022-06-20 | 2024-03-15 | 北京科技大学 | Monocular video-oriented non-Newtonian fluid simulation reconstruction method |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109063375A (en) * | 2018-09-07 | 2018-12-21 | 中山大学 | The analogy method and system of incompressible fluid based on secrecy and without divergence |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR102205845B1 (en) * | 2013-10-28 | 2021-01-22 | 삼성전자주식회사 | Method and apparatus for modeling based on particles |
-
2019
- 2019-12-30 CN CN201911390333.9A patent/CN111104753B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109063375A (en) * | 2018-09-07 | 2018-12-21 | 中山大学 | The analogy method and system of incompressible fluid based on secrecy and without divergence |
Non-Patent Citations (3)
| Title |
|---|
| 一种新的有效模拟不可压缩流体的SPH方法;张中兴等;《计算机应用研究》;20141029(第03期);全文 * |
| 基于N-S方程的真实感流体仿真;武斌博等;《微计算机信息》;20090405(第10期);全文 * |
| 核设施退役虚拟仿真中烟尘输运过程建模及算法研究;刘中坤等;《原子能科学技术》;20100920;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111104753A (en) | 2020-05-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111104753B (en) | Viscous incompressible fluid simulation method based on SPH | |
| Feppon et al. | Topology optimization of thermal fluid–structure systems using body-fitted meshes and parallel computing | |
| US10114911B2 (en) | Fluid structure interaction simulation method and apparatus, and computer-readable storage medium | |
| Wang et al. | Delaunay graph and radial basis function for fast quality mesh deformation | |
| Yu et al. | A high-order spectral difference method for unstructured dynamic grids | |
| Ren et al. | Dual-support smoothed particle hydrodynamics in solid: variational principle and implicit formulation | |
| Zhang et al. | High-fidelity aerostructural optimization with integrated geometry parameterization and mesh movement | |
| Da Ronch et al. | Adaptive design of experiments for efficient and accurate estimation of aerodynamic loads | |
| CN108717722A (en) | Fluid animation generation method and device based on deep learning and SPH frames | |
| Liu et al. | High-order particle method for solving incompressible Navier–Stokes equations within a mixed Lagrangian–Eulerian framework | |
| Greco et al. | A stabilized formulation with maximum entropy meshfree approximants for viscoplastic flow simulation in metal forming | |
| Martínez et al. | Updated Lagrangian free surface flow simulations with natural neighbour Galerkin methods | |
| De Corato et al. | Finite element formulation of fluctuating hydrodynamics for fluids filled with rigid particles using boundary fitted meshes | |
| Sverdrup et al. | An embedded boundary approach for efficient simulations of viscoplastic fluids in three dimensions | |
| Li et al. | Analysis of the accuracy and pressure oscillation of the lattice Boltzmann method for fluid–solid interactions | |
| Wang et al. | Implicit large eddy simulation of the NASA CRM high-lift configuration near stall | |
| Cui et al. | A high-order edge-based smoothed finite element (ES-FEM) method with four-node triangular element for solid mechanics problems | |
| Rossinelli et al. | Vortex methods for incompressible flow simulations on the GPU | |
| JP7042562B2 (en) | Optimal pressure projection for uncompressed transient and steady-state Navier-Stokes equations | |
| Wang et al. | Simulating fluid-structure interactions with a hybrid immersed smoothed point interpolation method | |
| Fan et al. | The point collocation method with a local maximum entropy approach | |
| Obrecht et al. | High-performance implementations and large-scale validation of the link-wise artificial compressibility method | |
| Ouyang et al. | A hybrid smoothed particle hydrodynamics coupled to a fictitious domain method for particulate flows and its application in a three-dimensional printing process | |
| Gao et al. | Predicting fluid–structure interaction with graph neural networks | |
| Cetinaslan | Position‐based simulation of elastic models on the GPU with energy aware gauss‐seidel algorithm |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |