CN117094239A - High-viscosity fluid rheological simulation method based on dynamic pressure correction - Google Patents
High-viscosity fluid rheological simulation method based on dynamic pressure correction Download PDFInfo
- Publication number
- CN117094239A CN117094239A CN202310504931.4A CN202310504931A CN117094239A CN 117094239 A CN117094239 A CN 117094239A CN 202310504931 A CN202310504931 A CN 202310504931A CN 117094239 A CN117094239 A CN 117094239A
- Authority
- CN
- China
- Prior art keywords
- density
- fluid
- correction
- constraint
- representing
- 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.)
- Pending
Links
- 239000012530 fluid Substances 0.000 title claims abstract description 146
- 238000012937 correction Methods 0.000 title claims abstract description 107
- 238000004088 simulation Methods 0.000 title claims abstract description 66
- 238000000034 method Methods 0.000 title claims abstract description 39
- 230000003068 static effect Effects 0.000 claims abstract description 14
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 9
- 239000002245 particle Substances 0.000 claims description 138
- 239000011159 matrix material Substances 0.000 claims description 28
- 239000013598 vector Substances 0.000 claims description 9
- 238000000518 rheometry Methods 0.000 claims description 5
- 239000011164 primary particle Substances 0.000 claims description 3
- 238000004364 calculation method Methods 0.000 abstract description 4
- 230000010355 oscillation Effects 0.000 abstract description 3
- 230000002238 attenuated effect Effects 0.000 abstract 1
- 239000008280 blood Substances 0.000 description 11
- 210000004369 blood Anatomy 0.000 description 11
- 230000008878 coupling Effects 0.000 description 4
- 238000010168 coupling process Methods 0.000 description 4
- 238000005859 coupling reaction Methods 0.000 description 4
- 230000005484 gravity Effects 0.000 description 4
- 238000001356 surgical procedure Methods 0.000 description 4
- 230000001133 acceleration Effects 0.000 description 3
- 230000001276 controlling effect Effects 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 230000010354 integration Effects 0.000 description 3
- 238000009825 accumulation Methods 0.000 description 2
- 230000009471 action Effects 0.000 description 2
- 238000013016 damping Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 239000000463 material Substances 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 230000001105 regulatory effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 238000005507 spraying Methods 0.000 description 2
- -1 batter Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000003292 glue Substances 0.000 description 1
- 235000012907 honey Nutrition 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000006188 syrup Substances 0.000 description 1
- 235000020357 syrup Nutrition 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/25—Design optimisation, verification or simulation using particle-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/12—Timing analysis or timing optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Computing Systems (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Fluid Mechanics (AREA)
- Operations Research (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention provides a high-viscosity fluid rheological simulation method based on dynamic pressure correction, which comprises the steps of firstly calculating the implicit viscosity force speed, taking a follow-up derivative of density as a source item, defining a correction force as an inverse gradient form of pressure, completing the calculation of a density dynamic correction coefficient, secondly taking the correction force as a prediction iteration quantity of a next speed field, finally taking a residual error of the actual density and the static density as the source item, expanding a Lagrange continuous equation by Euler integral, carrying out the calculation of solving a new divergence correction coefficient, and then taking the correction force as an iteration quantity of the last prediction speed. The invention can solve the problems of the traditional SPH algorithm that the truncation error on the speed change is caused, the volume error is accumulated along with the time, the volume is attenuated, the density correction steady-state error is overlarge, the occurrence of density oscillation and the instability of high-viscosity simulation are caused.
Description
Technical Field
The invention relates to the technical field of high-viscosity fluid rheology, in particular to a high-viscosity fluid rheology simulation method based on dynamic pressure correction.
Background
SPH (Smoothed Particle Hydrodynamics) A fluid simulation method is widely applied to various fields such as game interaction, film special effects, CFD numerical simulation, astronomy research, marine research and the like. Through many years of development history, the method can realize the characteristic simulation of different fluids and the simulation of multiphase flow coupling interaction, and can well meet the real-time performance and stability required by simulation scenes.
However, since the invention in 1977, SPH methods have all been problematic in that they suffer from density oscillations and volume decay over time due to four factors: (1) The coordinate updating of particles in the SPH method is realized based on a Lagrangian method, and the Euler integral discretization expression exists in the acting force items among the particles, so that the particles with inverse pressure gradient can be more attractive along with the reduction of the mutual distance due to the introduction of a kernel function, and the density of the current flow field can be increased; (2) The density error introduced in the SPH method is used for correcting the pressure, and the static correction coefficient can cause overshoot of the correction quantity each time, so that larger steady-state error is caused; (3) The Lagrangian continuous equation of the SPH method only limits the error of density, and no divergence constraint and no related formula, which would lead to problems of simulation and volume decay over time. (4) In the traditional SPH method, the angular momentum loss of boundary conditions in fluid-solid coupling of a complex model and a deformation model cannot be overcome: when the boundary particles are not enough to generate the pressure and the viscosity force are asymmetric, the loss of angular momentum is caused, and when part of particles are firstly contacted with the boundary, the particles close to the boundary are extruded by the upper layer particles due to the asymmetry of the pressure and can never stay at the boundary to rebound, and the speed loss at the boundary is caused by the asymmetry of the viscosity force, which is generally called false visual effect in the field.
Disclosure of Invention
Based on the method, the invention aims to provide a high-viscosity fluid rheological simulation method based on dynamic pressure correction, and the traditional SPH method is improved by setting a dynamic correction coefficient under the combined action of constant density constraint and non-divergence constraint so as to solve the defects of the traditional SPH method, thereby remarkably improving the stability of fluid simulation.
According to the invention, a high-viscosity fluid rheological simulation method based on dynamic pressure correction is provided, and comprises the following steps:
initializing the attribute of the high-viscosity fluid, wherein the attribute comprises speed, coordinates, a neighborhood matrix, density and a correction factor, and carrying out solution of implicit viscosity force according to a combined algorithm of Laplace approximation and implicit linear matrix equation set of a speed field and Jacobian iteration to obtain a speed simulation result under non-pressure;
before the pressure correction term is introduced, calculating a global time step of the high-viscosity fluid according to CFL conditions;
taking the residual error of the predicted density and the static density as a first source item, and feeding back and forth a correction coefficient for controlling a density constant constraint condition according to the first source item so as to establish a density constant constraint model under the density constant constraint condition, so that a speed field after the pressure correction item of the density constant constraint is obtained according to the density constant constraint model;
updating the position of the primary particles according to the velocity field after the pressure correction term of constant density constraint, and updating velocity field attributes based on the updated position, wherein the velocity field attributes comprise a neighborhood matrix, density and correction factors;
taking the updated random derivative of the density as a second source item, and utilizing the second source item to feedback control a correction coefficient without a divergence constraint condition, thereby establishing a divergence constraint-free model, and predicting a velocity field after passing through a pressure correction item without a divergence constraint according to the divergence constraint-free model;
when the simulation time is greater than a preset time threshold, outputting a speed field, and updating the flow field particle coordinates of the next time according to the speed field and the global time step.
Preferably, the step of obtaining the velocity simulation result under non-pressure according to the solution of the implicit viscous force by the combined algorithm of the Laplace approximation and the implicit linear matrix equation set of the velocity field and the Jacobian iteration comprises the following steps:
calculating to obtain a speed prediction corresponding to the viscous force term according to the following formula:
wherein,indicating the predicted speed of fluid particle i at time n+1,/->Indicating the velocity of fluid particle i at time n ρ i Represents the density of the fluid particles i, μ represents the viscosity coefficient,/->A laplace representing the velocity of fluid particle i at time n+1;
the laplace velocity of the fluid particle i at time n+1 is calculated according to the following formula:
wherein d is the spatial dimension, m j Representing the mass, ρ, of the fluid particles j j Indicating the density of the fluid particles j,a velocity vector indicating that fluid particle i points to fluid particle j,/->Coordinate vector indicating that fluid particle i points to fluid particle j,/->Indicating the spatial distance from fluid particle i to fluid particle j, h indicating the kernel function W ij Nuclear width of->Representing a kernel function W ij Is a gradient of (a).
Preferably, the step of obtaining the velocity simulation result under non-pressure by performing the implicit viscous force calculation according to the combined algorithm of the laplace approximation of the velocity field and the implicit linear matrix equation set and the jacobian iteration further comprises:
and carrying out Jacobian iterative solution according to the following formula to converge out the speed caused by the viscous force term at the next moment:
wherein the method comprises the steps of
Wherein Δt represents the global time step, A represents the neighbor matrix of the velocity field Laplacian, A ij Matrix elements representing the ith row and jth column of the neighborhood matrix,representing the mass difference between fluid particles i and j, μ being the viscosity coefficient, ++>Transpose of distance vector representing fluid particle i and fluid particle j, A ii Matrix elements representing the ith row and ith column of the neighborhood matrix.
Preferably, the step of calculating a global time step of the high viscosity fluid according to CFL conditions comprises:
the global time step is calculated according to the following formula:
wherein v is max Represents the maximum velocity of the particles in space, and d represents the particle diameter.
Preferably, the step of taking the residual error of the predicted density and the static density as a first source term and feeding back a correction coefficient for controlling the density constant constraint condition according to the first source term to establish a density constant constraint model under the density constant constraint condition comprises the following steps:
the predicted density is calculated according to the following formula:
wherein,the predicted density ρ of the fluid particle i at time n+1 is represented by 0 Representing the static density of the initialization setting, +.>Represents the predicted density of fluid particles i at time n, < >>Represents the predicted density of fluid particles j at time n, < >>Correction parameters representing fluid particles i under a constant density constraint model, +.>Representing the correction coefficient of the fluid particles j under the constant density constraint model;
the solving equation of the correction parameters under the density constant constraint model is constructed according to the following formula:
wherein alpha is i An inverse ratio representing a weighted sum of the neighborhood particles;
the solving equation for the satellite derivative of density is constructed according to the following formula:
wherein,the satellite derivative of the density of the fluid particles i representing the moment n+1,/o>Representing the pressure of the fluid particles i>Representing the pressure of fluid particle i against fluid particle j;
determining a pressure correction term subject to a density constant constraint according to the following formula:
wherein,a gradient representing the pressure of the fluid particles i;
bringing the pressure term into a solution equation of the correction parameter under the density constant constraint model to obtain the correction parameter under the density constant constraint model:
representing the mass weight sum,/of all fluid particles j within the kernel range>The sum of weights representing the square of the masses of all fluid particles j within the kernel range.
Preferably, the step of predicting the velocity field after the pressure correction term of the constant density constraint according to the constant density constraint model includes:
wherein,representing the predicted velocity of fluid particles i at time n+1, output via the density-invariant constraint model,/->The velocity at time n of the fluid particle i input into the density constant constraint model is represented.
Preferably, the updated satellite derivative of the density is taken as a second source term, and the second source term is utilized to feedback control a correction coefficient without a divergence constraint condition, so as to build a model without the divergence constraint:
discretizing a solving equation of the satellite derivative of the density to obtain:
and calculating a correction coefficient without a divergence constraint condition according to the following formula:
wherein,representing the correction parameters of the fluid particles i in the model without the divergence constraint.
Preferably, the step of predicting the velocity field after passing through the non-divergence-constraint pressure correction term according to the non-divergence-constraint model includes:
the velocity field is calculated according to the following formula:
wherein,the predicted speed of the fluid particle i at the time n+1, which is output by the non-divergence constraint model, is represented,representing the velocity of the fluid particle i at time n, input into the non-divergence-constrained model, +.>Representing the correction factor of the fluid particle j in the model without the divergence constraint.
Compared with the prior art, the invention has the following advantages:
(1) In the invention, the implicit jacobian iteration is adopted for the resolving of the viscous force, and the viscous force item does not carry the strain rate tensor attribute, so that the problem that the velocity gradient of the viscous force of boundary particles is not zero is avoided, the angular momentum loss and damping movement of the boundary are eliminated, and the boundary stability of the particles under the high-viscosity movement is ensured;
(2) According to the invention, the stability of simulation is ensured through dynamic correction of density and divergence, and the accumulation of divergence errors along with time can not occur when the time step is increased on the premise of ensuring conservation of momentum according to CFL conditions, so that the simulation of larger time step can be realized, and the performance of a computer consumed by the simulation is reduced;
(3) The residual errors of the predicted density and the static density are selected as source items, and the overshoot of the flow field density is regulated through the dynamic density correction coefficient, so that the final steady-state error range is ensured to be limited to 1.2%;
(4) The invention selects the follow-up derivative of the density as a source term, and limits the shrinkage of the volume of the flow field through the dynamic divergence correction coefficient, so that the volume of the flow field cannot be reduced along with the time, thereby ensuring the accuracy of the flow field;
(5) The correction factors contained in the dynamic correction coefficients selected by the invention can be used in density constraint and divergence constraint at the same time, so that the simulation consumption is reduced, and the correction factors are used as single particle attributes, thereby being convenient for local control;
(6) The particles storing the flow field attribute are uniformly distributed in the simulation, so that the convergence rate in the simulation can be effectively improved, and the simulation stability is improved;
(7) The density and divergence constraints of the invention can be processed in parallel by the GPU hardware equipment, thereby further reducing the time overhead in the rheological simulation process of the high-viscosity fluid.
Additional aspects and advantages of the invention will be set forth in part in the description which follows and, in part, will be obvious from the description, or may be learned by practice of the invention.
Drawings
FIG. 1 is a flow chart of a method for rheological simulation of a high viscosity fluid based on dynamic pressure correction in accordance with an embodiment of the present invention;
FIG. 2 is a comparison of haemorheology simulations at different viscosity coefficients;
FIG. 3 is a graph showing the comparison of flow field average density results;
FIG. 4 is a schematic diagram of angular momentum contrast of a flow field approaching a boundary;
fig. 5 is a diagram of an example of application in a dental surgical system.
The invention will be further described in the following detailed description in conjunction with the above-described figures.
Detailed Description
In order that the invention may be readily understood, a more complete description of the invention will be rendered by reference to the appended drawings. Several embodiments of the invention are presented in the figures. This invention may, however, be embodied in many different forms and should not be construed as limited to the embodiments set forth herein. Rather, these embodiments are provided so that this disclosure will be thorough and complete.
Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention belongs. The terminology used herein in the description of the invention is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The term "and/or" as used herein includes any and all combinations of one or more of the associated listed items.
Referring to fig. 1, a flow chart of a method for simulating rheological of a high viscosity fluid based on dynamic pressure correction according to an embodiment of the invention is shown, the method includes steps S01 to S06, wherein:
step S01: initializing the attribute of the high-viscosity fluid, wherein the attribute comprises speed, coordinates, a neighborhood matrix, density and a correction factor, and carrying out solution of implicit viscosity force according to a combined algorithm of Laplace approximation and implicit linear matrix equation set of a speed field and Jacobian iteration to obtain a speed simulation result under non-pressure;
it should be noted that the high viscosity fluid in this embodiment may be honey, blood, batter, syrup, glue, etc., and the rheological simulation of blood is firstly incompressible under normal atmospheric pressure conditions, and the flow has high viscosity, so that the Navier-Stokes equation needs to be used as the basis of the incompressible Newtonian fluid simulation, and the construction of the physical model of fluid motion is based on three assumptions: the first is that the movement of the fluid in space is continuous, i.e. the inflow and outflow per unit time is the same. Second is that the terms in the fluid velocity field are differentiable, such as density, pressure, etc. Third, satisfy Newton's law of three conservation. The N-S equation may also be referred to as a control equation for the fluid, including the continuous equation and the momentum equation:
the first term is a continuous equation, also called a mass conservation equation, describing that the fluid mass increase of the finite control element per unit time is equal to the fluid mass exiting the finite control element face. And the second term of the equation is a momentum equation, describing the relationship between the acceleration and the force of the fluid unit.
In some optional embodiments of the present invention, the non-pressure term in the N-S equation, that is, the gravity term and the viscous force term are solved, the gravity term is the gravity acceleration g, and the viscous force term, in the simulation of the high viscous force blood system, is necessarily required to use the solution of the implicit viscous force, the combined algorithm of the laplace approximation of the given velocity field and the implicit linear matrix equation set is used to realize the solution of the implicit viscous force, and in the simulation of the high viscosity blood, the velocity prediction corresponding to the viscous force term is given first:
wherein,indicating the predicted speed of fluid particle i at time n+1,/->Indicating the velocity of fluid particle i at time n ρ i Represents the density of the fluid particles i, μ represents the viscosity coefficient,/->A laplace representing the velocity of fluid particle i at time n+1;
and then taking the Laplace approximation of the velocity field into velocity prediction corresponding to the viscous force term:
wherein d is the spatial dimension, m j Representing the mass, ρ, of the fluid particles j j Indicating the density of the fluid particles j,a velocity vector indicating that fluid particle i points to fluid particle j,/->Coordinate vector indicating that fluid particle i points to fluid particle j,/->Indicating the spatial distance from fluid particle i to fluid particle j, h indicating the kernel function W ij Nuclear width of->Representing a kernel function w ij Is a gradient of (2);
the equation is written into a linear system, wherein the predicted speed at the next moment is set as an unknown quantity, and the speed caused by the viscous force term at the next moment can be converged by iteratively solving the linear system:
wherein the method comprises the steps of
Wherein Δt represents the global time step, A represents the neighbor matrix of the velocity field Laplacian, A ij Matrix elements representing the ith row and jth column of the matrix,representing the mass difference between fluid particles i and j, μ being the viscosity coefficient, ++>The transpose of the distance vectors of fluid particle i and fluid particle j is shown.
Finally, through jacobian iteration, in the simulation of a high-viscosity-force blood system, a linear equation set is solved through implicit integration, iteration is continued until the weighted residual error between the speed at the next moment and the speed of the particles in the previous moment is smaller than a certain threshold value, a loop is jumped out, and the speed caused by viscosity force at the next moment is predicted, for example, the following formula is an example of coordinate iteration:
x i [k+1] ＝x i [k] +α·diag(A) -1 (b-Ax [k] )
wherein x is i [k+1] For the next moment of coordinates of the fluid particles i, x i [k] For the last moment of the fluid particle i, α is the relaxation coefficient of the jacobian iteration, customized by the user, diag (A) -1 And b is an observation value of coordinates, namely coordinate values obtained by multiplying the CFL time step by the speed, and A is a neighbor matrix of the Laplacian of the speed field at the previous moment.
Referring to fig. 2, in order to simulate the improved implicit viscosity force, when μ increases, the particles slow down due to the increase in viscosity coefficient, resulting in stacking of particles in the third and fourth graphs, and when μ increases to 0.75, the velocity field does not collapse instantaneously. It follows that it is suggested to simulate blood ejection using an implicit viscous force of μ=0.5 in dental surgery.
When μ=0.5 shown in fig. 3 is the same, the blood splash simulation applied to the gingival cutting in the dental surgery is performed, the first-spraying blood is firstly contacted with air under the action of the implicit viscosity force, and then the sprayed blood is smoothed by laplace because of the high outlet speed, so that the stack effect is formed with the first-spraying blood, and the expected simulation effect of high viscosity force in the dental surgery blood simulation is achieved
Step S02: before the pressure correction term is introduced, calculating a global time step of the high-viscosity fluid according to CFL conditions;
updated once during each simulation cycle, and used together by two pressure constraints, which greatly reduces the computational effort of the iteration. All non-pressure supplied accelerations such as implicit viscosity force, gravity force, etc. are calculated prior to the introduction of the pressure correction term. And then adjusting the time step in the simulation in real time according to the CFL condition:
wherein v is max The maximum speed of the particles in the space is represented, d the particle diameter, and the time step is self-adaptive according to the speed attribute of the speed field, and no artificial setting is needed.
Step S03: taking the residual error of the predicted density and the static density as a first source item, and feeding back and forth a correction coefficient for controlling a density constant constraint condition according to the first source item so as to establish a density constant constraint model under the density constant constraint condition, so that a speed field after the pressure correction item of the density constant constraint is obtained according to the density constant constraint model;
in order to achieve density constant pressure correction, the goal is to achieve residual error of actual density and static densityIs determined by explicit euler integration of the lagrangian continuity equation:
wherein,represents n +Predicted density ρ of fluid particle i at time 1 0 Representing the static density of the initialization setting, +.>Represents the predicted density of fluid particles i at time n, < >>Represents the predicted density of fluid particles j at time n, < >>Correction parameters representing fluid particles i under a constant density constraint model, +.>Representing the correction coefficient of the fluid particles j under the constant density constraint model;
wherein correction parameters are corrected prior to initializationPre-calculation similar to no divergence pressure correction is used to improve convergence, +.>Expressed as:
wherein alpha is i An inverse ratio representing the weighted sum of the neighborhood particles, as determined by the lagrangian continuous equation;
and then constructing a solving equation of the satellite derivative of the density according to the following formula:
wherein,the satellite derivative of the density of the fluid particles i representing the moment n+1,/o>Representing the pressure of the fluid particles i>Representing the pressure of fluid particle i against fluid particle j;
determining a pressure correction term subject to a density constant constraint according to the following formula:
wherein,a gradient representing the pressure of the fluid particles i;
bringing the pressure term into the Lagrangian continuous equation and correcting the parametersIt is proposed that correction factor alpha is included i ：
Representing the mass weight sum,/of all fluid particles j within the kernel range>The sum of weights representing the square of the masses of all fluid particles j within the kernel range.
After obtaining the pressure correction term of constant density constraint, predicting the speed field of the next time step:
wherein,representing the predicted velocity of fluid particles i at time n+1, output via the density-invariant constraint model,/->The velocity at time n of the fluid particle i input into the density constant constraint model is represented.
Step S04: updating the position of the primary particles according to the velocity field after the pressure correction term of constant density constraint, and updating velocity field attributes based on the updated position, wherein the velocity field attributes comprise a neighborhood matrix, density and correction factors;
wherein the factor alpha is pre-calculated i Is used by a density constant pressure correction solver to calculate a correction pressure to correct the density error of each neighborhood particleThe time integration in line 13 calculates the new position of the particle. Thus in the new position, all neighborhoods N i Density ρ i Sum factor alpha i Real-time updates are required.
Step S05: taking the updated random derivative of the density as a second source item, and utilizing the second source item to feedback control a correction coefficient without a divergence constraint condition, thereby establishing a divergence constraint-free model, and predicting a velocity field after passing through a pressure correction item without a divergence constraint according to the divergence constraint-free model;
it should be noted that in some alternative embodiments of the present invention, the velocity field is not generally non-diverged after the density bias correction is completed. Finally, intervention from the non-divergence constraint is required, in order to achieve non-divergence pressure correction, the objective is to achieve a minimization of the satellite derivative of density, which is determined by discretizing the lagrangian continuous equation:
the non-divergence-constraint pressure correction is realized by a non-divergence correction parameter calculated according to the following formula:
wherein,representing the correction parameters of the fluid particles i under the non-divergence constraint model;
and (3) bringing the satellite derivative of the density into a non-divergence correction parameter to obtain a non-divergence constraint pressure correction term, so as to predict the speed field of the next time step:
i.e. the velocity field is calculated according to the following formula:
wherein,the predicted speed of the fluid particle i at the time n+1, which is output by the non-divergence constraint model, is represented,representing the velocity of the fluid particle i at time n, input into the non-divergence-constrained model, +.>Representing the correction factor of the fluid particle j in the model without the divergence constraint.
Referring to fig. 4, in fig. 4 (a), the blue line is a static density curve, the red line is an average density curve of the conventional SPH method, the average density is in a continuously increasing state before 100 frames, the average density of the fluid particles after 100 frames shows a continuously correcting overshoot oscillation condition, the overshoot amount is 62% around 100 frames, and the steady state error is maintained at a level of 22% when the fluid is run to 300 frames, which directly causes the volume of the fluid to be unequal in different time periods. FIG. 4 (b) is a graph showing the average density of the constant density non-scattering SPH at a rate of 0.013kg/m over time, 100 frames ago 3 The average density after 175 frames will be maintained at an average of 1.06 ρ per frame 0 Constant level (6% steady state error). From this, it can be seen that a constant total volume of fluid can be maintained at different time points by a constant density non-divergence SPH-based hemorheology simulation method, this experiment; the accuracy, rapidity and stability of the pressure correction of the dynamic coefficient proposed by the invention are demonstrated.
Referring to fig. 5, the haemorheological simulation of the virtual dental procedure is not independent of the coupling with the gingival tissue, in the conventional SPH method, the darker the color, the closer to zero the F-norm representing the strain rate tensor, because the pressure and viscous forces cannot overcome the loss of angular momentum at the fluid-solid coupling boundary conditions, as in fig. 5 (a). The F-norm of the strain rate tensor is not zero due to the absence of particles at the boundary. The haemorheology simulation method based on density-constant non-divergence SPH largely overcomes the problem of angular momentum loss at the boundary, as shown in fig. 5 (b), and the experiment proves the accuracy, robustness and expandability of the invention applied in the field of virtual dental surgery.
Step S06: when the simulation time is greater than a preset time threshold, outputting a speed field, and updating the flow field particle coordinates of the next time according to the speed field and the global time step.
When the simulation time (program running time) is greater than a time threshold which is considered to be set, the speed field output by the non-divergence constraint model is the final speed field, the final speed field is multiplied by the time step which is equal to the final coordinate, after the coordinate is updated according to the final coordinate, the program running time is cleared again, the steps S01 to S06 are repeatedly executed, and if the simulation time is less than or equal to the time threshold, the steps S01 to S05 are repeated until the simulation time is greater than the time threshold, so that the final speed field is output.
Compared with the prior art, the invention has the following advantages:
(1) In the invention, the implicit jacobian iteration is adopted for the resolving of the viscous force, and the viscous force item does not carry the strain rate tensor attribute, so that the problem that the velocity gradient of the viscous force of boundary particles is not zero is avoided, the angular momentum loss and damping movement of the boundary are eliminated, and the boundary stability of the particles under the high-viscosity movement is ensured;
(2) According to the invention, the stability of simulation is ensured through dynamic correction of density and divergence, and the accumulation of divergence errors along with time can not occur when the time step is increased on the premise of ensuring conservation of momentum according to CFL conditions, so that the simulation of larger time step can be realized, and the performance of a computer consumed by the simulation is reduced;
(3) The residual errors of the predicted density and the static density are selected as source items, and the overshoot of the flow field density is regulated through the dynamic density correction coefficient, so that the final steady-state error range is ensured to be limited to 1.2%;
(4) The invention selects the follow-up derivative of the density as a source term, and limits the shrinkage of the volume of the flow field through the dynamic divergence correction coefficient, so that the volume of the flow field cannot be reduced along with the time, thereby ensuring the accuracy of the flow field;
(5) The correction factors contained in the dynamic correction coefficients selected by the invention can be used in density constraint and divergence constraint at the same time, so that the simulation consumption is reduced, and the correction factors are used as single particle attributes, thereby being convenient for local control;
(6) The particles storing the flow field attribute are uniformly distributed in the simulation, so that the convergence rate in the simulation can be effectively improved, and the simulation stability is improved;
(7) The density and divergence constraints of the invention can be processed in parallel by the GPU hardware equipment, thereby further reducing the time overhead in the rheological simulation process of the high-viscosity fluid.
In the description of the present specification, a description referring to terms "one embodiment," "some embodiments," "examples," "specific examples," or "some examples," etc., means that a particular feature, structure, material, or characteristic described in connection with the embodiment or example is included in at least one embodiment or example of the present invention. In this specification, schematic representations of the above terms do not necessarily refer to the same embodiments or examples. Furthermore, the particular features, structures, materials, or characteristics described may be combined in any suitable manner in any one or more embodiments or examples.
The foregoing examples illustrate only a few embodiments of the invention and are described in detail herein without thereby limiting the scope of the invention. It should be noted that it will be apparent to those skilled in the art that several variations and modifications can be made without departing from the spirit of the invention, which are all within the scope of the invention. Accordingly, the scope of protection of the present invention is to be determined by the appended claims.
Claims (8)
1. A method of rheological simulation of a high viscosity fluid based on dynamic pressure correction, the method comprising:
initializing the attribute of the high-viscosity fluid, wherein the attribute comprises speed, coordinates, a neighborhood matrix, density and a correction factor, and carrying out solution of implicit viscosity force according to a combined algorithm of Laplace approximation and implicit linear matrix equation set of a speed field and Jacobian iteration to obtain a speed simulation result under non-pressure;
before the pressure correction term is introduced, calculating a global time step of the high-viscosity fluid according to CFL conditions;
taking the residual error of the predicted density and the static density as a first source item, and feeding back and forth a correction coefficient for controlling a density constant constraint condition according to the first source item so as to establish a density constant constraint model under the density constant constraint condition, so that a speed field after the pressure correction item of the density constant constraint is obtained according to the density constant constraint model;
updating the position of the primary particles according to the velocity field after the pressure correction term of constant density constraint, and updating velocity field attributes based on the updated position, wherein the velocity field attributes comprise a neighborhood matrix, density and correction factors;
taking the updated random derivative of the density as a second source item, and utilizing the second source item to feedback control a correction coefficient without a divergence constraint condition, thereby establishing a divergence constraint-free model, and predicting a velocity field after passing through a pressure correction item without a divergence constraint according to the divergence constraint-free model;
when the simulation time is greater than a preset time threshold, outputting a speed field, and updating the flow field particle coordinates of the next time according to the speed field and the global time step.
2. The method for rheological simulation of high viscosity fluid based on dynamic pressure correction according to claim 1, wherein the step of obtaining the velocity simulation result under non-pressure by solving the implicit viscosity force according to the combined algorithm of the laplace approximation of the velocity field and the implicit linear matrix equation set and jacobian iteration comprises:
calculating to obtain a speed prediction corresponding to the viscous force term according to the following formula:
wherein,indicating the predicted speed of fluid particle i at time n+1,/->Indicating the velocity of fluid particle i at time n ρ i Represents the density of the fluid particles i, μ represents the viscosity coefficient,/->A laplace representing the velocity of fluid particle i at time n+1;
the laplace velocity of the fluid particle i at time n+1 is calculated according to the following formula:
wherein d is the spatial dimension, m j Representing the mass, ρ, of the fluid particles j j Indicating the density of the fluid particles j,a velocity vector indicating that fluid particle i points to fluid particle j,/->Coordinate vector indicating that fluid particle i points to fluid particle j,/->Indicating the spatial distance from fluid particle i to fluid particle j, h indicating the kernel function W ij Nuclear width of->Representing a kernel function W ij Is a gradient of (a).
3. The method for rheological simulation of a high viscosity fluid based on dynamic pressure correction according to claim 2, wherein the step of obtaining the velocity simulation result under non-pressure by performing the solution of the implicit viscosity force according to the combined algorithm of the laplace approximation of the velocity field and the implicit linear matrix equation set and jacobian iteration further comprises:
and carrying out Jacobian iterative solution according to the following formula to converge out the speed caused by the viscous force term at the next moment:
wherein the method comprises the steps of
Wherein Δt represents the global time step, A represents the neighbor matrix of the velocity field Laplacian, A ij Matrix elements representing the ith row and jth column of the neighborhood matrix,representing the mass difference between fluid particles i and j, μ being the viscosity coefficient, ++>Transpose of distance vector representing fluid particle i and fluid particle j, A ii Matrix elements representing the ith row and ith column of the neighborhood matrix.
4. A method of rheological simulation of a highly viscous fluid based on dynamic pressure correction as set forth in claim 3 wherein the step of calculating a global time step for the highly viscous fluid based on CFL conditions comprises:
the global time step is calculated according to the following formula:
wherein v is max Represents the maximum velocity of the particles in space, and d represents the particle diameter.
5. The method for rheological simulation of a high viscosity fluid based on dynamic pressure correction according to claim 4, wherein the step of taking the residual error of the predicted density and the static density as a first source term and feeding back the correction coefficient for controlling the constant density constraint condition according to the first source term to build the constant density constraint model under the constant density constraint condition comprises:
the predicted density is calculated according to the following formula:
wherein,the predicted density ρ of the fluid particle i at time n+1 is represented by 0 Representing the static density of the initialization settings,represents the predicted density of fluid particles i at time n, < >>Represents the predicted density of fluid particles j at time n, < >>Correction parameters representing fluid particles i under a constant density constraint model, +.>Representing the correction coefficient of the fluid particles j under the constant density constraint model;
the solving equation of the correction parameters under the density constant constraint model is constructed according to the following formula:
wherein alpha is i An inverse ratio representing a weighted sum of the neighborhood particles;
the solving equation for the satellite derivative of density is constructed according to the following formula:
wherein,the satellite derivative of the density of the fluid particles i representing the moment n+1,/o>The pressure of the fluid particles i is indicated,representing the pressure of fluid particle i against fluid particle j;
determining a pressure correction term subject to a density constant constraint according to the following formula:
wherein,a gradient representing the pressure of the fluid particles i;
bringing the pressure term into a solution equation of the correction parameter under the density constant constraint model to obtain the correction parameter under the density constant constraint model:
representing the mass weight sum,/of all fluid particles j within the kernel range>The sum of weights representing the square of the masses of all fluid particles j within the kernel range.
6. The method of dynamic pressure correction-based high viscosity fluid rheology simulation according to claim 5, wherein predicting a velocity field after a pressure correction term of a density constant constraint according to a density constant constraint model includes:
wherein,the predicted velocity of the fluid particle i at the time n+1 output by the density constant constraint model is represented,the velocity at time n of the fluid particle i input into the density constant constraint model is represented.
7. The dynamic pressure correction-based high viscosity fluid rheology simulation method according to claim 6, wherein the updated satellite derivative of density is used as a second source term, and the second source term is used to feedback control the correction coefficient without the divergence constraint condition, so as to build a model without the divergence constraint:
discretizing a solving equation of the satellite derivative of the density to obtain:
and calculating a correction coefficient without a divergence constraint condition according to the following formula:
wherein,representing the correction parameters of the fluid particles i in the model without the divergence constraint.
8. The dynamic pressure correction-based high viscosity fluid rheology simulation method according to claim 7, wherein the step of predicting a velocity field after passing through a non-divergence-constrained pressure correction term according to a non-divergence-constrained model includes:
the velocity field is calculated according to the following formula:
wherein,indicating the predicted velocity of fluid particle i at time n+1, without the output of the divergence-constrained model,/->Representing the velocity of the fluid particle i at time n, input into the non-divergence-constrained model, +.>Representing the correction factor of the fluid particle j in the model without the divergence constraint.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310504931.4A CN117094239A (en) | 2023-05-06 | 2023-05-06 | High-viscosity fluid rheological simulation method based on dynamic pressure correction |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310504931.4A CN117094239A (en) | 2023-05-06 | 2023-05-06 | High-viscosity fluid rheological simulation method based on dynamic pressure correction |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN117094239A true CN117094239A (en) | 2023-11-21 |
Family
ID=88777735
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310504931.4A Pending CN117094239A (en) | 2023-05-06 | 2023-05-06 | High-viscosity fluid rheological simulation method based on dynamic pressure correction |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117094239A (en) |
-
2023
- 2023-05-06 CN CN202310504931.4A patent/CN117094239A/en active Pending
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Müller et al. | Position based dynamics | |
| Baaijens | A fictitious domain/mortar element method for fluid–structure interaction | |
| Breuer et al. | Fluid–structure interaction using a partitioned semi-implicit predictor–corrector coupling scheme for the application of large-eddy simulation | |
| Kadapa et al. | A fictitious domain/distributed Lagrange multiplier based fluid–structure interaction scheme with hierarchical B-spline grids | |
| US10210287B2 (en) | Augmented material point method for simulating phase changes and varied materials | |
| US20150161810A1 (en) | Position based fluid dynamics simulation | |
| CN109002630B (en) | Rapid simulation method for super-elastic material | |
| Ganzenmüller et al. | Hourglass control for Smooth Particle Hydrodynamics removes tensile and rank-deficiency instabilities: Hourglass control for SPH | |
| Liu et al. | Automatic scaled boundary finite element method for three-dimensional elastoplastic analysis | |
| Loppi et al. | Locally adaptive pseudo-time stepping for high-order flux reconstruction | |
| US11270041B2 (en) | Position-based dynamics simulation | |
| US20230061175A1 (en) | Real-Time Simulation of Elastic Body | |
| Su et al. | Energy Conservation for the Simulation of Deformable Bodies. | |
| Chen et al. | An alternative updated Lagrangian formulation for finite particle method | |
| Thürey et al. | Interactive free surface fluids with the lattice Boltzmann method | |
| Raessi et al. | A semi‐implicit finite volume implementation of the CSF method for treating surface tension in interfacial flows | |
| JP2004219389A (en) | Instantaneous buckling model for woven fabric, history model, woven fabric simulation method based on the model, and computer-readable recording medium stored with program performing the same | |
| Bukač et al. | Time-adaptive partitioned method for fluid-structure interaction problems with thick structures | |
| Li et al. | A partitioned framework for coupling LBM and FEM through an implicit IBM allowing non-conforming time-steps: Application to fluid-structure interaction in biomechanics | |
| Dagenais et al. | A prediction-correction approach for stable sph fluid simulation from liquid to rigid | |
| Dalmon et al. | Fluids-membrane interaction with a full Eulerian approach based on the level set method | |
| Joneidi et al. | Isogeometric boundary integral analysis of drops and inextensible membranes in isoviscous flow | |
| US20060129351A1 (en) | Information processing apparatus and method for solving simultaneous linear equations | |
| CN117094239A (en) | High-viscosity fluid rheological simulation method based on dynamic pressure correction | |
| Zhang et al. | A moving mesh finite difference method for non-monotone solutions of non-equilibrium equations in porous media |
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 |