Summary of the invention
The object of the present invention is to provide the construction method of Free Surface flow model in a kind of improved semi implicit algorithm.
For achieving the above object, the present invention has adopted following technical scheme:
1) physical parameter of measurement flow field and fluid, physical parameter comprises dimension of flow field, fluid density and viscosity;
2) initial distribution of setting particle;
3) through step 2) after, the interim velocity field of explicit calculating;
4) through after the step 3), the estimation particle position;
5) through after the step 4), calculate population density, then build-up pressure Poisson equation;
6) through after the step 5), find the solution the pressure gradient term of pressure Poisson equation according to pressure field, then revise particle rapidity and position.
In the described step 3), adopt the interim velocity field of particle method large eddy simulation turbulence model (subparticle stress model) explicit calculating of introducing, the large eddy simulation turbulence model comprises the dynamic Smagorinsky model of static Smagorinsky model or Lagrangian Form.
The Large eddy simulation method of Eulerian mesh method simulation turbulent flow is coupled in the particle method, the governing equation of large eddy simulation will have more than the governing equation of routine a subparticle stress item, and adopting model coefficient is the subparticle stress item that the dynamic Smagorinsky model of the static Smagorinsky model of constant or Lagrangian Form is constructed the large eddy simulation turbulence model.
The interim velocity field of the described explicit calculating of particle method large eddy simulation turbulence model may further comprise the steps:
1) when adopting the Smagorinsky model, the subparticle stress model is as follows:
Wherein,
The deformation rate tensor that can separate yardstick, δ
IjExpression Crow Buddhist nun restrains tensor, v
tRepresent respectively eddy viscosity and tubulence energy with k; Tubulence energy item and pressure term are coupled and find the solution, and eddy viscosity calculates by following formula:
Wherein, Δ represents to filter yardstick, C
sIt is model coefficient;
2) adopt C
s=0.15 model coefficient as the subparticle stress model is called static Smagorinsky model; Perhaps adopt the method for double filtration, according to Germano equation calculating model coefficient everywhere, be called dynamic Smagorinsky model; The Germano equation is as follows:
Wherein
Subparticle stress after the expression large scale is filtered,
Represent the subparticle stress after the double filtration of large small scale;
3) for dynamic Smagorinsky model, suppose that dynamic Smagorinsky model coefficient is irrelevant with the filtration yardstick, with formula (3), (4) described Smagorinsky model substitution Germano equation, can obtain following relation:
Definition
And
Then the dynamic Smagorinsky model coefficient of every bit is
Wherein,
Strain rate tensor after the expression large scale is filtered, L
KkExpression subparticle stress increment L
IjThe diagonal line component;
4) adopt the time averaging method of following the trail of particle to determine dynamic Smagorinsky model coefficient, be defined in the time averaging M on the particle trajectory
IjL
IjAnd M
IjM
IjBe respectively I
LMAnd I
MMThen dynamically the Smagorinsky model coefficient is as follows:
5) obtain after the model coefficient, adopt the combination of gradient former in the MPS method and the dispersion models subparticle stress model that disperses, finish behind subparticle stress model discrete, adopt explicit algorithm to find the solution subparticle stress item, the first order derivative of the speed of using in the subparticle stress model all is based on the velocity distribution of layer in a period of time.
In the interim velocity field of described explicit calculating, adopt to become viscosity Newtonian fluid model and process constitutive equation and be the non-Newtonian fluid of μ=f (| γ |), μ is kinetic viscosity, and γ is shearing rate.
It is the Free Surface of the non-Newtonian fluid of μ=f (| the γ |) problem that flows for constitutive equation, non-Newtonian fluid is considered as becoming the viscosity Newtonian fluid to be similar to, formed the MPS method that is applicable to the Flows of Non-newtonian Fluids problem, in each time step, utilize the velocity field computational flow shearing rate γ everywhere of a upper time step, recycling constitutive equation μ=f (| γ |) calculates flow field kinetic viscosity μ everywhere, this is equivalent to non-Newtonian fluid is changed into the Newtonian fluid that becomes viscosity, the Fluid Control Equation of Newtonian fluid stands good for non-Newtonian fluid, and difference is that this moment, μ was parameter.
It is the Free Surface of the non-Newtonian fluid of μ=f (| the γ |) problem that flows for constitutive equation, can abandon using MPS method kernel function commonly used, and employing can not only embody the rule of particle interaction, its space derivative is the cubic spline kernel function of correct SPH method also, and then the shear stress of the discrete non-Newtonian fluid of divergence discrete solution of employing SPH, formed the MPS method that is applicable to the Flows of Non-newtonian Fluids problem.
The present invention is mainly for improved semi implicit algorithm (Moving Particle Semi-implicit, MPS), turbulence model and two kinds of non-Newtonian models of subparticle stress have been incorporated in the particle method, so that the fine model engineering class of particle class methods energy is mobile.The present invention proposes based on the MPS method of large eddy simulation and for the MPS method of the non-Newtonian fluid of equation shape such as μ=f (| γ |), expanded the scope of application of MPS method.
Embodiment
The invention will be further described below in conjunction with accompanying drawing.
The overall calculation process of algorithm of the present invention comprises following a few step (referring to shown in Figure 1):
1. measure the physical parameter of flow field and fluid, comprise dimension of flow field, fluid density and viscosity;
2. set the particle initial placement according to initial flow-field;
3. the turbulence model of the movement-based particle semi implicit algorithm of employing the present invention proposition and the construction method fluid flow of non-Newtonian models are simulated;
With particle model at the State-output of each time step to computer screen, show the real-time flow process in flow field.
The below improves the each several part in the flow process and does detailed introduction.Discrete solution and the implementation step of basic improved semi implicit algorithm can be with reference to Koshizuka S.and Oka Y.Moving-particle semi-implicit method for fragmentation of incompressible fluid.Nuclear Science and Engineering1996; 123:421-434.
1 is used for the subparticle stress model of particle method large eddy simulation
Fig. 1 step 8,9, the 10th, the step of Considering Turbulence effect has wherein been used Static and dynamic subparticle stress model, next two kinds of models and discrete way is described in detail.For the N-S equation of incompressibility fluid, filter the LES governing equation that obtains in physical space as follows:
In the formula,
It is the large scale speed component after filtering.The right of equation has increased after filtering
This is called as subparticle stress item (Sub-Particle Stress, SPS) in particle method.
Filter operation is carried out in mesh scale in the grid class methods, should carry out at particle size in the particle class methods.In the present invention, still adopt the Smagorinsky model of widespread use in the grid class methods to come the modeling of subparticle stress model, and adopt static model and two kinds of methods of dynamic model to calculate the model coefficient of Smagorinsky model.The key that the subparticle stress model is introduced in the particle method is: how (1) adopts the Smagorinsky model that disperses of the discrete model in the particle class methods; (2) the subparticle model after how will dispersing is coupled in the algorithm of MPS.
When adopting the Smagorinsky model, the subparticle stress model is as follows:
Wherein,
The deformation rate tensor that can separate yardstick, δ
IjExpression Crow Buddhist nun restrains tensor, v
tRepresent respectively eddy viscosity and tubulence energy with k, they can characterize the intensity of turbulent flow.Tubulence energy item and pressure term are coupled and find the solution, and eddy viscosity can calculate by following formula:
Wherein, Δ represents to filter yardstick, C
sIt is model coefficient.
The present invention adopts two kinds of methods to determine the coefficient of Smagorinsky model: the empirical constant C that at first can adopt widespread use in the grid class methods
s=0.15 coefficient as the subparticle stress model, such scheme are called static Smagorinsky model; Another kind method is to adopt the method for double filtration according to Germano equation calculating model coefficient everywhere, and such scheme is called dynamic Smagorinsky model.For for simplicity, the small scale Δ
1Filter with subscript "-" expression, large scale Δ
2Filter with subscript "~" expression.Use Δ
1, Δ
2The subparticle stress that double filtration obtains can be than separately with the yardstick Δ
2Filter subparticle stress and have more a part; The increment L of this part subparticle stress
IjExpression, its expression formula is as follows:
The Germano equation is as follows:
Suppose dynamic Smagorinsky model coefficient C
D(annotate: in order to narrate conveniently, dynamically there are following relation in Smagorinsky model coefficient and static Smagorinsky model coefficient
Irrelevant with the filtration yardstick, with the above-mentioned Germano equation of above-mentioned Smagorinsky model substitution, can obtain following relation:
Definition
Then the dynamic Smagorinsky model coefficient of every bit is
Can determine the dynamic Smagorinsky model coefficient of any time, any position according to following formula, but the C that calculates like this
DMay alter a great deal, what cause calculating is unstable, therefore adopts average method to obtain model coefficient.Because the motion that the MPS method can be followed the trail of particle easily, the present invention adopts Meneveau C, Lund TS, Cabot WH.A Lagrangian dynamic subgrid-scale model of turbulence.Journal of Fluid Mechanics 1996; The time averaging method of the tracking particle that 319:353-85 proposes is determined dynamic Smagorinsky model coefficient, is defined in the time averaging M on the particle trajectory
IjL
IjAnd M
IjM
IjAs follows:
Wherein, and W (t-t ') be the weight function of an exponential form, so that I
LMAnd I
MMCan find the solution by following discrete method:
Wherein, T is the slack time of current time, and n represents current iterative steps, and Meveneau advises that its computing formula is as follows:
T=1.5Δ
2(I
LMI
MM)
-1/8(13)
Through above-mentioned calculating, can obtain the I at the different moment, diverse location place
MMAnd I
LM, dynamically the Smagorinsky model coefficient is as follows:
Obtain after the dynamic Smagorinsky model coefficient, be updated to the Smagorinsky model (formula (3) and formula (4)) of front, just can calculate the subparticle stress of dynamic Smagorinsky model.
Differential term in the SPS model all is first order derivative, can adopt gradient former in the MPS method and the combination of dispersion models to disperse.Specifically, subparticle stress
Model in, the first order derivative item of the speed that only relates on different coordinate directions; And gradient former can provide the single order partial derivative of a speed component on each coordinate direction, therefore, can be fully calculates subparticle stress by the gradient discrete model of MPS itself
After obtaining subparticle stress, just can calculate again the acceleration of the particle that is produced by subparticle stress item according to the divergence operator in the MPS method.In dynamic Smagorinsky model, also use large scale in the filtration that can separate on the yardstick, according to the characteristics of MPS discrete model, the filter operation of large scale can adopt average weighted method to disperse, and is as follows:
Wherein, the expression formula of Gaussian function G (r, Δ) is as follows under two-dimensional case:
Wherein x represents the position coordinates of particle, and r is two distances between the particle.
Finish behind subparticle stress model discrete, the present invention adopts explicit algorithm to find the solution subparticle stress item, and the first order derivative of the speed of namely using in the model all is based on the velocity distribution of layer in a period of time.
The present invention adopts as shown in Figure 2 dam break model to calculate the effect of the example check turbulence model of parallel efficiency, highly is being 0.5m, and length is in the tank of 1.0m a fluid column to be arranged, liquid-column height is H=0.5m, the width of fluid column is L=0.25m, and liquid is Newtonian fluid, and density is 1000kg/m
3, kinematic viscosity is 0.00001m
2/ s.Be configured at particle to have calculated the Dam Break Problems based on the Static and dynamic subparticle stress of Smagorinsky model under 50 * 100 the condition, for the experimental result with Martin compares, the position Z in dam break forward position carries out nondimensionalization with L, and time t adopts
Carry out nondimensionalization, its result and experimental result and the result contrast of only considering the molecule stickiness are as shown in Figure 3.As can be seen from Figure 3: only consider that the result of molecular viscosity and experimental result have larger difference, the position that shows as the dam break forward position of simulation will be ahead of the experimental result of synchronization; Consider that analog result and experimental result behind the subparticle stress effect are identical better; Simultaneously, can also see, in the later stage of calculating, static Smagorinsky model than dynamic Smagorinsky model more near experimental result.This has proved that employing subparticle stress model can improve analog result really.
2 to become the approximate non-Newtonian models of viscosity Newtonian fluid
Fig. 1 step 4,5,7 has provided the first disposal route of non-Newtonian fluid, is about to non-Newtonian fluid and is considered as becoming the viscosity Newtonian fluid and is similar to.In each time step, utilize the velocity field computational flow shearing rate γ everywhere of a upper time step, recycling constitutive equation μ=f (| γ |) calculates flow field kinetic viscosity μ everywhere.Through above-mentioned processing, be equivalent to non-Newtonian fluid is changed into the Newtonian fluid that becomes viscosity, the Fluid Control Equation of Newtonian fluid stands good for the non-Newtonian fluid of constitutive equation shape such as μ=f (| γ |), and difference is that this moment, μ was parameter.
For the two-dimensional flow problem, the expression formula of shearing rate γ is
Wherein, x, y represent respectively two coordinate directions of two-dimensional coordinate, and u, v represent respectively the speed of velocity on two coordinate directions;
In the MPS method, the partial derivative of speed can be obtained with following equation:
Wherein, x
i, x
j, y
i, y
jRepresent respectively the position of particle on coordinate direction, subscript i represents the center particle, and subscript j represents ambient particles, r
i, r
jThe coordinate vector of expression particle, w represents kernel function, n
0Represent initial population density,
The expression Hamiltonian;
After obtaining shearing rate γ, recycling constitutive equation μ=f (| γ |) can calculate flow field kinetic viscosity μ everywhere on the current time level.According to the Newtonian fluid equation of momentum:
Wherein,
The random derivative of expression speed, ρ represents particle density, and p represents pressure, and u represents velocity, and F represents body force;
At this moment, as long as non-Newtonian fluid is considered as becoming the Newtonian fluid of viscosity, formula (22) stands good, and μ has only wherein become variable from constant.Therefore, the explicit part calculating formula of MPS method calculating Flows of Non-newtonian Fluids is:
Wherein,
The interim speed of expression particle i,
The interim speed of n the time step of expression particle i, Δ t represents time step,
The temporary position of expression particle i,
The temporary position of n the time step of expression particle i;
3 adopt the non-Newtonian model of cubic spline kernel function
Fig. 1 step 4,5,6 have provided the second disposal route of non-Newtonian fluid, because the kernel function form of MPS method is fairly simple coarse, cause the scope of application of discrete scheme narrower, expansion is relatively poor, for the Free Surface of the non-Newtonian fluid of constitutive equation shape such as μ=f (| the γ |) problem that flows, the present invention also can not use MPS method kernel function commonly used, and employing can not only embody the rule of particle interaction, its space derivative is correct SPH method cubic spline kernel function commonly used also, and then the divergence of the shear stress tensor τ of the discrete non-Newtonian fluid of divergence discrete solution of employing SPH, the Free Surface of simulation non-Newtonian fluid flows.
For non-Newtonian fluid, equation of momentum writing:
Wherein, τ is shear stress tensor, f=F/ ρ;
Interparticle interaction is all finished by kernel function in MPS method and the SPH method, but the kernel function of two kinds of methods requires to some extent difference.In the MPS method, kernel function only is a weight function, so this kernel function needs only energy reaction particle spacing to the strong and weak rule that interacts.The most frequently used MPS method kernel function is shown below:
Wherein, r
eBe the effective radius of influence of particle, r is interparticle distance.
Because the SPH method acts directly on the discrete of various operators on the kernel function, so kernel function not only will embody the rule of particle interaction, will guarantee that also its space derivative is also correct.In order to adopt the divergence operator discrete solution of SPH method, the present invention adopts the cubic spline kernel function commonly used suc as formula the SPH shown in (27).
Wherein, q=r/h, h are smooth length.Because the fluid in the MPS method is real incompressible fluid, so h is definite value, h=r
e/ 2.
Gradient operator of the present invention and Laplace operator still adopt traditional MPS method discrete scheme, adopt the divergence of the discrete shear stress tensor τ of following SPH divergence operator discrete scheme commonly used: (τ represents shear stress tensor, also can be referred to as shear stress when not emphasizing tensor character)
Wherein, m represents mass particle, w
IjExpression particle j is to the kernel function of particle i, and N represents the sum of particle i ambient particles.The expression formula of particle density ρ is
Similarly, because the MPS method can embody the characteristic of real incompressible fluid, the m in the following formula, ρ is constant, particle density and fluid macroscopic view equal density.With formula (29) substitution formula (28), obtain the divergence operator discrete scheme of this paper:
Wherein:
Wherein, r
IjDistance between expression particle i and particle j.
It should be noted that τ is a second-order tensor, it and one-dimensional vector
Dot product result be one-dimensional vector.
Like this, calculating the mobile the second MPS algorithm arrangement of non-Newtonian fluid free surface organizes as follows: in each time step, utilize the velocity field computational flow shearing rate γ everywhere of a upper time step, recycling constitutive equation μ=f (| γ |) calculates flow field kinetic viscosity μ everywhere, then calculate shear stress tensor τ according to formula (32), and utilize formula (30) that the divergence of shear stress tensor τ is dispersed.
τ=2μ(γ)D (32)
μ is viscosity, and D is strain rate tensor.In the two-dimensional problems, strain rate tensor D is expressed as:
Adopt and document [1] (Shao, S.D.and E.Lo, Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface.Advances in water resources2003.26 (7): 787-800.) and document [2] (Komatina, D.and M.Jovanovic, Experimental study of steady and unsteady free surface flows with water-clay mixtures.Journal of Hydraulic research 1997.35 (5): 579-590.) identical mud-aqueous mixtures Dam Break Problems is verified two kinds of non-Newtonian models that the present invention proposes.
The example moulded dimension is 2m a length as shown in Figure 4, and degree of tilt is S
0Mud-aqueous mixtures that the degree of depth is 0.1m is housed in=0.1% the tank.Be consistent with document [1], it is considered as the Cross fluid, its constitutive equation is:
Simplifying method in the list of references [1], formula (34) can be write following form:
Wherein, μ
BBe Ben-Hur surrender viscosity, τ
BBe the Ben-Hur yield stress.
Fluid density is 1200kg/m
3, yield stress is taken as 25.0Pa, and Ben-Hur surrender viscosity is 0.07Ns/m
2
The employing of the present invention primary consistent with document [1] arranged, the primary spacing is 0.01m, be evenly arranged altogether 2000 fluid particles and 1256 wall particles, result of calculation when t=0.1s, 0.3s and 0.6s as shown in Figure 5 and Figure 6, all coincide better with document [1], proved the correctness of two kinds of non-Newtonian fluid algorithms that the present invention proposes.Be further contrast, the distribution of particles when having compared t=0.37s and 0.6s, as shown in Figure 7.Can see, two kinds of non-Newtonian fluid algorithm simulation results that the present invention proposes all coincide well with experimental result and SPH method analog result.By contrast, the dam break front end Free Surface and the experimental result that adopt the non-Newtonian fluid algorithm of cubic spline kernel function to obtain are best in consistent manner.