CN109726431B - Self-adaptive SPH fluid simulation method based on average kernel function and iterative density change rate - Google Patents

Self-adaptive SPH fluid simulation method based on average kernel function and iterative density change rate Download PDF

Info

Publication number
CN109726431B
CN109726431B CN201811409064.1A CN201811409064A CN109726431B CN 109726431 B CN109726431 B CN 109726431B CN 201811409064 A CN201811409064 A CN 201811409064A CN 109726431 B CN109726431 B CN 109726431B
Authority
CN
China
Prior art keywords
particle
sph
support domain
radius
change rate
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
Application number
CN201811409064.1A
Other languages
Chinese (zh)
Other versions
CN109726431A (en
Inventor
朱鉴
李泓澍
陈炳丰
蔡瑞初
郝志峰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Guangdong University of Technology
Original Assignee
Guangdong University of Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Guangdong University of Technology filed Critical Guangdong University of Technology
Priority to CN201811409064.1A priority Critical patent/CN109726431B/en
Publication of CN109726431A publication Critical patent/CN109726431A/en
Application granted granted Critical
Publication of CN109726431B publication Critical patent/CN109726431B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention provides a self-adaptive SPH fluid simulation method based on an average kernel function and an iterative density change rate, under an SPH frame, the particle density and a support domain are solved in an iterative mode, the physical quantity of a target particle is interpolated by using a neighbor particle group with equal total mass, the stress of the target particle is calculated, so that the displacement of the target particle is obtained, the state of the target particle is updated, and aiming at the self-adaptive particle support domain, a symmetrical object based on the particle stress is obtained; and finally, by means of the average kernel function, the problem of asymmetric interaction force of particles introduced by the variable support domain is solved, so that the simulation system is more stable.

Description

Self-adaptive SPH fluid simulation method based on average kernel function and iterative density change rate
Technical Field
The invention relates to the technical field of computer simulation, in particular to a self-adaptive SPH fluid simulation method based on an average kernel function and an iterative density change rate.
Background
In the field of computer graphics, in recent years, the Particle-based lagrangian method is gradually becoming a main tool for implementing fluid simulation, and particularly, the Smooth Particle Hydrodynamics (SPH) method, which can describe the process of medium motion more naturally, simulate more fluid details such as foam, water bloom and the like, and process more complicated fluid surfaces, and is small in calculation amount and easy to implement.
However, due to the discretization of particles and the limitation of a simulation area, the SPH method inevitably increases the numerical dissipation, which will result in the loss of many tiny details in the fluid simulation process, and the SPH needs a large number of particles to simulate a high-precision large-scale scene, while a large number of high-density particles will result in a massive calculation requirement, because the conventional SPH particle support area is fixed, in areas with different particle distribution densities, the difference of the number of neighboring particles in the support area will occur, it is difficult to obtain a kernel approximation progress with a consistent calculation area, resulting in an excessive particle interpolation error; it is always compelled for users to make a trade-off between simulation effect and system performance, and therefore it is always very meaningful to improve the computational accuracy of the SPH method.
Disclosure of Invention
Aiming at the defects of the prior art, the invention provides a self-adaptive SPH fluid simulation method based on an average kernel function and an iteration density change rate, which can solve the calculation error caused by the influence of a fixed particle support domain in the traditional SPH method; and by designing a set of particle neighbor search scheme and average interpolation kernel function average kernel function, the problem of asymmetric particle interaction force caused by a variable support domain is solved. According to the system experiment result, the simulation effect is more vivid, and the fluid movement is more solid.
The technical scheme of the invention is as follows: an adaptive SPH fluid simulation method based on an average kernel function and an iterative density change rate comprises the following steps:
s1), obtaining the density change rate of the SPH particle fluid and the change of the radius r of the particle support domain through iteration;
s2), setting a particle support domain radius filtering function;
s3) obtaining the radius r of the particle support domain of each SPH particle according to the step S1)iRadius of the particle support domain r for each SPH particleiLet r beijGet riAnd rjOf the obtained R value is obtained by a filter functionijWhen the Euclidean distance d between two particlesijLess than RijThat is, the particle i and the particle j are adjacent particles;
and S4) carrying out physical quantity interpolation solving on each SPH particle according to the average gradient kernel function.
In the method, in step S1), the change of the SPH particle support domain radius r is obtained from the fluid density change rate, and the calculation formula is a change rate iterative formula:
Figure BDA0001878063610000021
Figure BDA0001878063610000022
where dt is the fluid simulation time step, ρiIs the density of the particle i, mjIs the mass of particle j, vijIs the difference in velocity between particle i and particle j, i.e. vij＝vi-vjAnd W is a kernel function,
Figure BDA00018780636100000216
is W pairs (x)i-xj) Computing a gradient, where xiIs the position of the particle i, xjIs the position of the particle j and,
Figure BDA0001878063610000023
for W to symmetric particle support domain rijCalculating partial derivative, D is dimension of fluid simulation system, riIs the radius of the particle i and,
Figure BDA0001878063610000024
is the rate of change of the radius of the particle support domain for particle i,
Figure BDA0001878063610000025
is the density change rate of the particle i;
using an iterative approach to solve, i.e. ri≥rj，rij＝ri，
Figure BDA0001878063610000026
Otherwise, rij＝rj，
Figure BDA0001878063610000027
Preset of
Figure BDA0001878063610000028
Figure BDA0001878063610000029
The change rate of the radius of the particle support domain for the particle j at the previous time step is substituted into the above formula (1) to obtain
Figure BDA00018780636100000210
Corrected by the formula (2)
Figure BDA00018780636100000211
Then, the method is iterated through the formulas (1) and (2) until
Figure BDA00018780636100000212
Is stabilized to obtain
Figure BDA00018780636100000213
By using
Figure BDA00018780636100000214
And a time step dt, updating the density of the particle i and the particle support domain radius.
In the above method, in step S2), the particle support domain radius filter function is set as:
Figure BDA00018780636100000215
wherein R isminAnd RmaxLower and upper limits for the radius of the particle support domain, typically 2 rparticle≤Rmin≤Rmax≤8·rparticle，rparticleIs the SPH particle radius; the problem of too large a support domain caused by too small a density of isolated particles is prevented by a filtering function.
In the above method, in step S4), the average kernel function is represented as a mean value of a conventional kernel function, including its gradient and laplacian form:
Figure BDA0001878063610000031
Figure BDA0001878063610000032
Figure BDA0001878063610000033
wherein i, j is the particle index, xiIs the position of the particle i, RiThe particle support domain for the filtered particle i.
The method also comprises the step of searching particle neighbors according to the particle support domain, wherein the neighbor searching mode is a full-pairing searching method, a linked list searching method or a tree searching method, and a derivative method of the method.
The beneficial effects of the invention are as follows: according to the method, through iteration of a group of equation sets, the density change rate and the support domain change rate of the particles with smaller errors than those of the traditional method are obtained, and the calculation result of the interaction force of the particles is effectively and accurately improved through more accurate density; secondly, the adaptive particle support domain enables the interpolation effect of the force to be more excellent and the calculation result to be closer to the physical fact; and finally, by means of the average kernel function, the problem of asymmetric interaction force of particles introduced by the variable support domain is solved, so that the simulation system is more stable.
Drawings
FIG. 1 is a flow chart of the SPH simulation process of the method of the present invention;
FIG. 2 is a schematic diagram of a conventional SPH method particle support domain;
FIG. 3 is a diagram illustrating the difference in the number of neighboring particles in the particle support domain caused by the density difference in the conventional SPT method;
FIG. 4 is a schematic diagram showing the result of adjusting the radius of the particle support domain by the method of the present invention;
FIG. 5 is a comparison graph of the simulation process of the SPH according to the conventional SPH and the method of the present invention, wherein FIG. a is a diagram of the simulation process of the conventional SPH, and FIG. b is a diagram of the simulation process of the SPH according to the method of the present invention.
Detailed Description
The following further describes embodiments of the present invention with reference to the accompanying drawings:
as shown in fig. 1, an adaptive SPH fluid simulation method based on an average kernel function and an iterative density change rate includes the following steps:
s1), creating a fluid model, initializing parameters, and building a particle model for fluid simulation, wherein the particle model attributes comprise but are not limited to position, density, speed, mass and radius of a support domain;
s2) searching for particle neighbors according to the particle support domain, wherein the neighbor searching mode is a full-pairing searching method, a linked list searching method or a tree searching method, and a derivative method of the methods, and the tree searching method is recommended to be used, because the tree searching method is more excellent in a system with non-uniform particle support domains;
particle support domain radius r for each SPH particleiLet rijGet riAnd rjAccording to the filter function RijWhen the Euclidean distance d of two particlesijLess than RijThat is, particle i and particle j are neighboring particles to each other, that is:
Figure BDA0001878063610000041
wherein R isminAnd RmaxLower and upper limits for the radius of the particle support domain, typically 2 rparticle≤Rmin≤Rmax≤8·rparticle，rparticleIs the SPH particle radius; the problem that the support domain is too large due to too small isolated particle density is prevented from being generated through a filtering function;
s3), iterative formulation solving the particle density and support domain, SPH is an interpolation method of particle systems, using which the field quantities defined only at discrete particle positions can be computed at any position in space, for which SPH uses a radially symmetric smoothing kernel to distribute the quantities in a local neighborhood of each particle, according to which a scalar a is weighted and interpolated with the contributions of all particles at the x position:
Figure BDA0001878063610000042
wherein j is the index of the neighboring particle of the particle i, and mjIs the mass of the particle j, xiDenotes the position of the particle i, pjIs the density of particle j, W is the kernel function, AjIs the field magnitude of the discrete particle j;
then there are: ρ (x)i)＝∑jmjW, deriving it to yield:
Figure BDA0001878063610000043
where dt is the fluid simulation time step, ρiIs the density of the particles i, mjIs the mass of particle j, vijIs the difference in velocity between particle i and particle j, i.e. vij＝vi-vjAnd W is a function of a kernel,
Figure BDA0001878063610000047
is W pair (x)i-xj) Computing a gradient where xiIs the position of the particle i and,
Figure BDA0001878063610000044
for W pairs of symmetric particle support domains rijCalculating a partial derivative, riIs the radius of the particle i and,
Figure BDA0001878063610000045
is the rate of change of the radius of the particle support domain for particle i,
Figure BDA0001878063610000046
is the density change rate of the particle i;
in SPH, the radius of the particle support domain is selected as much as possibleThe number of the neighbors of the particle is kept unchanged in calculation so as to obtain the kernel approximation precision consistent with the full field, and the radius r of the support domain of the particle i is setiThe region containing a constant mass MsumI.e. by
Figure BDA0001878063610000051
And (5) obtaining the following results by differentiating the left side and the right side of the step (5) with respect to time:
Figure BDA0001878063610000052
in the formula, V is the volume of a particle support domain under a 3-dimensional condition, S is the area of the particle support domain under a 2-dimensional condition, and D is the dimension of a fluid simulation system;
using an iterative approach to solve, i.e. ri≥rj，rij＝ri，
Figure BDA0001878063610000053
Otherwise, rij＝rj，
Figure BDA0001878063610000054
Preset of
Figure BDA0001878063610000055
Figure BDA0001878063610000056
Substituting the radius change rate of the particle support domain for the particle j at the previous time step into the formula (1) to obtain
Figure BDA0001878063610000057
Corrected by the formula (2)
Figure BDA0001878063610000058
Then, the method is iterated through the formulas (1) and (2) until
Figure BDA0001878063610000059
Is stabilized to obtain
Figure BDA00018780636100000510
By using
Figure BDA00018780636100000511
And a time step dt, updating the density of the particle i and the radius of the particle support domain;
s4), calculating the stress of the particles, and expressing momentum conservation by using a Navier-Stokes equation:
Figure BDA00018780636100000512
wherein g is an external force density field, mu is the viscosity of the fluid,
Figure BDA00018780636100000513
represents a simulated pressure, ρ g represents a simulated external force,
Figure BDA00018780636100000514
to simulate viscous forces.
So resultant force density field
Figure BDA00018780636100000515
This makes it possible to obtain:
Figure BDA00018780636100000516
in the formula, aiIs the acceleration of particle i;
according to the relevant physical formula, the following can be known:
Figure BDA00018780636100000517
Figure BDA00018780636100000518
the stability, accuracy and speed of the SPH method depend to a large extent on the choice of the smoothing kernel W, so that the conventional kernel is chosen to be:
Figure BDA0001878063610000061
due to the asymmetric particle interaction force problem introduced by the variable support domain, the problem is solved using an average kernel function, WijThe mean, expressed as the mean of a conventional kernel function, including its gradient and laplacian form:
Figure BDA0001878063610000062
Figure BDA0001878063610000063
Figure BDA0001878063610000064
wherein i, j is the particle index, xiIs the position of the particle i, RiA particle support domain for filtered particles i;
s5) updating the particle state according to the state obtained in the step S4)
Figure BDA0001878063610000065
Updating the speed and position information of the particles;
s6), judging whether the simulation is finished or not, and if not, continuing to execute the steps from S2) to S4) of the next time step.
FIG. 2 is a schematic diagram of a particle support domain of a conventional SPH method, and it can be seen that the conventional SPH has a fixed particle support domain; as can be seen from fig. 3, the conventional SPH uses a fixed particle support domain, which causes a difference in the number of neighboring particles in the support domain when the density is different; as can be seen in fig. 4: by the method, the particle support domain is adjusted, so that the number of the neighbor particles in the support domain tends to be the same; and as can be seen from fig. 5, the simulation results are closer to the physical reality, and the fluid motion is more concise.
The foregoing embodiments and description have been presented only to illustrate the principles and preferred embodiments of the invention, and various changes and modifications may be made therein without departing from the spirit and scope of the invention as hereinafter claimed.

Claims (5)

1. An adaptive SPH fluid simulation method based on an average kernel function and an iterative density change rate is characterized by comprising the following steps:
s1), obtaining the density change rate of the SPH particle fluid and the change of the radius r of the particle support domain through iteration; the change of the SPH particle support domain radius r is obtained according to the fluid density change rate, and the calculation formula is a change rate iterative formula:
Figure FDA0003856655010000011
Figure FDA0003856655010000012
where dt is the fluid simulation time step, ρiIs the density of the particle i, mjIs the mass of the particle j, vijIs the difference in velocity between particle i and particle j, i.e. vij＝vi-vjAnd W is a kernel function,
Figure FDA0003856655010000013
is W pairs (x)i-xj) Computing a gradient, where xiIs the position of the particle i, xjIs the position of particle j;
Figure FDA0003856655010000014
for W to symmetric particle support domain rijPartial derivative is calculated, D is the dimension of the fluid simulation system, riIs the radius of the particle i and,
Figure FDA0003856655010000015
is the rate of change of the radius of the particle support domain for particle i,
Figure FDA0003856655010000016
is the density change rate of the particle i;
using an iterative approach to solve, i.e. ri≥rj，rij＝ri，
Figure FDA0003856655010000017
Otherwise, rij＝rj，
Figure FDA0003856655010000018
Preset of
Figure FDA0003856655010000019
Figure FDA00038566550100000110
The change rate of the radius of the particle support domain for the particle j at the previous time step is substituted into the above formula (1) to obtain
Figure FDA00038566550100000111
Corrected by the formula (2)
Figure FDA00038566550100000112
Then, the method is iterated through the formulas (1) and (2) until the method is finished
Figure FDA00038566550100000113
Is stabilized to obtain
Figure FDA00038566550100000114
By using
Figure FDA00038566550100000115
And a time step dt, updating the density of the particle i and the radius of the particle support domain;
s2), setting a particle support domain radius filtering function;
Figure FDA00038566550100000116
wherein R isminAnd RmaxThe lower limit and the upper limit of the radius of the particle support domain are taken as 2 rparticle≤Rmin≤Rmax≤8·rparticle，rparticleIs the SPH particle radius; the problem that the support domain is too large due to too small isolated particle density is prevented from being generated through a filtering function;
s3) obtaining the radius r of the particle support domain of each SPH particle according to the step S1)iRadius of the particle support domain r for each SPH particleiLet rijGet riAnd rjOf the obtained R value is obtained by a filter functionijWhen the Euclidean distance d of two particlesijLess than RijThat is, the particle i and the particle j are adjacent particles;
s4), carrying out physical quantity interpolation solving on each SPH particle according to an average gradient kernel function, wherein the average kernel function is expressed as the average value of a conventional kernel function and comprises the gradient and Laplace operator form:
Figure FDA0003856655010000021
Figure FDA0003856655010000022
Figure FDA0003856655010000023
wherein i, j is the particle index, xiIs the position of the particle i, RiThe particle support domain for the filtered particle i.
2. The adaptive SPH fluid simulation method based on the average kernel function and the iterative density change rate according to claim 1, characterized in that: the simulation method further includes creating a fluid model and initializing parameters, building a particle model for fluid simulation, the particle model attributes including, but not limited to, position, density, velocity, mass, support domain radius.
3. The adaptive SPH fluid simulation method based on average kernel function and iterative density change rate of claim 1, wherein: in the step S3), a neighbor searching mode is adopted for searching the particle neighbors according to the particle support domain, and the neighbor searching mode is a full-pairing searching method, a linked list searching method or a tree searching method and a derivative method of the method.
4. The adaptive SPH fluid simulation method based on average kernel function and iterative density change rate of claim 1, wherein: step S4) also comprises the step of calculating the stress of the particles by using a Navier-Stokes momentum conservation equation, namely:
Figure FDA0003856655010000031
wherein g is an external force density field, mu is a viscosity coefficient of the fluid,
Figure FDA0003856655010000032
representing a simulated pressure, pg representing a simulated external force,
Figure FDA0003856655010000033
in order to simulate the viscous forces,
Figure FDA0003856655010000034
is the rate of change of the particle velocity;
the resultant density field is therefore:
Figure FDA0003856655010000035
this gives:
Figure FDA0003856655010000036
in the formula, aiIs the acceleration of particle i;
according to the relevant physical formula, the following can be known:
Figure FDA0003856655010000037
Figure FDA0003856655010000038
the stability, accuracy and speed of the SPH method depend on the choice of the smoothing kernel W, so that the conventional kernel is chosen to be:
Figure FDA0003856655010000039
5. the adaptive SPH fluid simulation method based on average kernel function and iterative density change rate of claim 4, wherein: the simulation method further comprises the step of obtaining
Figure FDA00038566550100000310
UpdatingAnd (4) judging whether the simulation is finished or not according to the speed and the position information of the particles, and if not, continuing to execute the steps from S2) to S4) of the next time step.
CN201811409064.1A 2018-11-23 2018-11-23 Self-adaptive SPH fluid simulation method based on average kernel function and iterative density change rate Active CN109726431B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811409064.1A CN109726431B (en) 2018-11-23 2018-11-23 Self-adaptive SPH fluid simulation method based on average kernel function and iterative density change rate

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811409064.1A CN109726431B (en) 2018-11-23 2018-11-23 Self-adaptive SPH fluid simulation method based on average kernel function and iterative density change rate

Publications (2)

Publication Number Publication Date
CN109726431A CN109726431A (en) 2019-05-07
CN109726431B true CN109726431B (en) 2022-11-01

Family

ID=66295792

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811409064.1A Active CN109726431B (en) 2018-11-23 2018-11-23 Self-adaptive SPH fluid simulation method based on average kernel function and iterative density change rate

Country Status (1)

Country Link
CN (1) CN109726431B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113011075B (en) * 2021-03-01 2024-06-11 南京师范大学 Free surface identification method suitable for multi-resolution particle method

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108269299A (en) * 2017-01-04 2018-07-10 北京航空航天大学 A kind of viscous fluid modeling method based on SPH method approximate solutions

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR101267627B1 (en) * 2011-06-13 2013-05-27 한국과학기술원 SPH Fluid simulation method and system for Multi-Level Vorticity, recording medium for the same
US20150161810A1 (en) * 2013-12-10 2015-06-11 Nvidia Corporation Position based fluid dynamics simulation
CN105320782B (en) * 2014-06-16 2019-06-21 复旦大学 A kind of characteristic size grade CMP process emulation mode for considering polishing fluid and influencing
CN104200015B (en) * 2014-08-20 2017-06-16 清华大学 A kind of fluid simulation method and device
CN106650064B (en) * 2016-12-09 2019-07-26 华东师范大学 A kind of condensation phenomenon emulation mode based on particle model
CN107025332B (en) * 2017-03-07 2021-05-14 华南理工大学 Visualization method for microscopic water diffusion process on fabric surface based on SPH

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108269299A (en) * 2017-01-04 2018-07-10 北京航空航天大学 A kind of viscous fluid modeling method based on SPH method approximate solutions

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"Study of droplet deformation, heat-conduction and solidification using incompressible smoothed particle hydrodynamics method ";Hong-bing Xiong等;《9th International Conference on Hydrodynamics》;20101015;全文 *

Also Published As

Publication number Publication date
CN109726431A (en) 2019-05-07

Similar Documents

Publication Publication Date Title
Umetani et al. Learning three-dimensional flow for interactive aerodynamic design
Lian et al. Implementation of regularized isogeometric boundary element methods for gradient‐based shape optimization in two‐dimensional linear elasticity
Losasso et al. Simulating water and smoke with an octree data structure
McNamara et al. Fluid control using the adjoint method
US20080275677A1 (en) System, methods, and computer readable media, for product design using coupled computer aided engineering models
US20080162090A1 (en) System, methods, and computer readable media, for product design using t-spline deformation
CN108763683B (en) New WENO format construction method under trigonometric function framework
US20080319722A1 (en) Water Particle Manipulation
CN108197072B (en) High-precision intermittent Galerkin artificial viscous shock wave capturing method based on weighted conservative variable step
CN113850008A (en) Self-adaptive grid disturbance domain updating acceleration method for aircraft aerodynamic characteristic prediction
CN107451383B (en) Calibration method for initial bed sand gradation of plane two-dimensional water sand mathematical model
CN109783935B (en) Implementation method for improving splash fluid stability based on ISPH
CN110717269A (en) Fluid surface detail protection method based on grid and particle coupling
CN109726431B (en) Self-adaptive SPH fluid simulation method based on average kernel function and iterative density change rate
Arami Fadafan et al. Moving particle semi-implicit method with improved pressures stability properties
KR100588000B1 (en) Apparatus and method for capturing free surface of fluid in computer animation
CN110598283B (en) Fluid simulation method based on SPH kernel function
KR101562863B1 (en) Method for simulating fluid flow by using the lattice Boltzmann theory and recording medium for performing the method
JP5750091B2 (en) Fluid simulation method
KR100779993B1 (en) Method for simulating fluid with applying contol value to pressure field, recording medium and apparatus thereof
KR101514806B1 (en) Method and apparatus of fluid simulation
Cotela Dalmau et al. Simulation of two-and three-dimensional viscoplastic flows using adaptive mesh refinement
JP6167554B2 (en) Channel shape optimization method and channel shape optimization apparatus
JP7395456B2 (en) Simulation device and program
Ghoneim A new technique for numerical simulation of dendritic solidification using a meshfree interface finite element method

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