CN102867094A - Establishment method for free surface flow model in moving particle semi-implicit algorithm - Google Patents

Abstract

The invention provides an establishment method for a free surface flow model in a moving particle semi-implicit algorithm, comprising the following steps of: introducing a turbulence model which comprises a static Smagorinsky model and a dynamic Smagorinsky model in the form of Lagrange; treating a non-Newtonian fluid having a constitutive equation like mu = f(the absolute value ofgamma) by adopting a variable-viscosity Newtonian fluid model; and introducing a cubic spline kernel function, discretizing the shear stress of the non-Newtonian fluid by adopting the divergence discretization scheme of a smooth particle fluid dynamic method, and treating the free surface flow of the non-Newtonian fluid. According to the method provided by the invention, the sub-particle stress model of turbulence can be considered, and the method provided by the invention can be used for calculation for the non-Newtonian fluid.

Description

The construction method of Free Surface flow model in a kind of improved semi implicit algorithm

Technical field

The present invention relates to the numerical simulation field that turbulent flow and non-Newtonian fluid free surface flow, relate to the optimization of turbulent flow and non-Newtonian fluid numerical model.

Background technology

In ship hydrodynamics, space flight, hydraulic engineering and petroleum engineering, all there is a large amount of Free Surfaces problem that flows.For example: large-scale carrier fluid boats and ships must be considered the impact for ship running stability of rocking of internal liquid in when design, and liquid rocking belowdecks is exactly the typical Free Surface problem that flows.The rubble flow disaster that causes of dam break and for example, the formation development of rubble flow also are the typical Free Surface problems that flows.The processes such as common evaporation and condensation, inkjet printing, chemical reaction all are accompanied by the Free Surface phenomenon that flows in the daily life.So common and for example this is important for phenomenon just because of Free Surface flows, and is one of focus of Fluid Flow in A research field to its research always.

Along with the development of computer technology, numerical simulation is widely used in the fields such as physics, chemistry, biology, Aero-Space, has been considered to since last century the third scientific research means except theoretical analysis and experiment are observed.Observe than experiment, numerical simulation has that efficient is high, cost less, the easy advantage such as reproduction.Mostly relate to expensive, large scale or microscale research object because Free Surface flows, experimental study is difficult to carry out, so numerical simulation also becomes research freedom surface one of the most strong instrument that flows.

Free Surface flows, and conventional mesh class methods that the numerical simulation study field adopts are foremost PIC (Particle-In-Cell) method, MAC (Marker And Cell) method, VOF (Volume of Fluid) method and Level Set method.Although these class methods flow and obtained many progress on the problem processing Free Surface, they arrange and accurately follow the trail of at the interface or the aspects such as reconstruct all are faced with the difficulties such as algorithm complexity, calculated amount are large at the irregular interface grid.

The particle class methods are fast-developing in recent years a kind of algorithms, because their employings is the motion of following the trail of particle under Lagrangian framework, therefore the position that can follow the trail of Free Surface is easy to the emulation Free Surface and flows.Smoothed Particle Hydrodynamics (SPH) and improved semi implicit algorithm (MPS) are two kinds of common particle class methods, the SPH method adopts the explicit processing pressure item of state equation, therefore fluid is approximate incompressible, and MPS implicit expression processing pressure item, fluid is really incompressible; The differentiating operator discrete solution of SPH method acts directly on the kernel function, therefore kernel function is had higher requirement, and the differentiating operator of MPS is then discrete by the particle weighted mean, and lower to the kernel function requirement, form is fairly simple coarse.

In addition, the common Free Surface of nature and engineering field flows and is mostly the turbulent flow of non-Newtonian fluid, and such as Dam Break Problems, the current of Dam Break Problems have carried a large amount of mud more, needs to adopt non-Newtonian fluid to be used for modeling, and is complicated turbulent flow.Therefore the particle class methods are applied to above-mentioned mobile in, need to introduce turbulence model and non-Newtonian models at particle method.The MPS method is mainly used in the flow field problem of the less Reynolds number of Newtonian fluid at present, studies less for non newtonian, turbulence model.

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:

${\stackrel{\‾}{\mathrm{\τ}}}_{\mathrm{ij}}={2v}_t{\stackrel{\‾}{S}}_{\mathrm{ij}}-\frac{2}{3}{\mathrm{k\δ}}_{\mathrm{ij}}---\left(3\right)$

Wherein, ${\stackrel{\‾}{S}}_{\mathrm{ij}}=\frac{1}{2}(\frac{\∂{\stackrel{\‾}{u}}_i}{{\∂x}_j}+\frac{\∂{\stackrel{\‾}{u}}_j}{{\∂x}_i})$ 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:

$v_t={\left(C_s\mathrm{\Δ}\right)}^2{\left(2{\stackrel{\‾}{S}}_{\mathrm{ij}}{\stackrel{\‾}{S}}_{\mathrm{ij}}\right)}^{1/2}---\left(4\right)$

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:

$L_{\mathrm{ij}}={\widetilde{\mathrm{\τ}}}_{\mathrm{ij}}-\stackrel{~}{{\stackrel{\‾}{\mathrm{\τ}}}_{\mathrm{ij}}}---\left(6\right)$

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

$C_D=\frac{M_{\mathrm{ij}}{L}_{\mathrm{ij}}}{M_{\mathrm{ij}}{M}_{\mathrm{ij}}}---\left(8\right)$

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:

$C_D=C_s^2=\frac{I_{\mathrm{LM}}}{I_{\mathrm{MM}}}---\left(14\right)$

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.

Description of drawings

Fig. 1 is the calculation process of algorithm of the present invention;

Fig. 2 is the physical dimension of Dam Break Problems;

Fig. 3 is the contrast of Dam Break Problems analog result and experimental result;

Fig. 4 is Cross fluid dam break example schematic diagram;

Fig. 5 is to become the approximate non-Newtonian fluid MPS method simulation Cross fluid dam break result of viscosity Newtonian fluid; (a) t=0.1s, (b) t=0.3s, (c) t=0.6s; T represents the dam break development time;

Fig. 6 is for adopting the non-Newtonian model MPS method simulation Cross fluid dam break result of cubic spline kernel function; (a) t=0.1s, (b) t=0.3s, (c) t=0.6s; T represents the dam break development time;

Fig. 7 be t=0.37s(Fig. 7 a) with t=0.6s(Fig. 7 b) time Cross fluid dam break Free Surface distribute; T represents the dam break development time.

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:

$\frac{\∂{\stackrel{\‾}{u}}_i}{{\∂x}_i}=0---\left(1\right)$

$\frac{{\∂\stackrel{\‾}{u}}_i}{\∂t}+\frac{\∂{\stackrel{\‾}{u}}_i{\stackrel{\‾}{u}}_j}{\∂x_j}=-\frac{1}{\mathrm{\ρ}}\frac{\∂\stackrel{\‾}{p}}{{\∂x}_i}+v\frac{{\∂}^2{\stackrel{\‾}{u}}_i}{{\∂x}_j{\∂x}_j}+g_i+\frac{\∂{\stackrel{\‾}{\mathrm{\τ}}}_{\mathrm{ij}}}{{\∂x}_j}---\left(2\right)$

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:

${\stackrel{\‾}{\mathrm{\τ}}}_{\mathrm{ij}}={2v}_t{\stackrel{\‾}{S}}_{\mathrm{ij}}-\frac{2}{3}{\mathrm{k\δ}}_{\mathrm{ij}}---\left(3\right)$

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:

$v_t={\left(C_s\mathrm{\Δ}\right)}^2{\left(2{\stackrel{\‾}{S}}_{\mathrm{ij}}{\stackrel{\‾}{S}}_{\mathrm{ij}}\right)}^{1/2}---\left(4\right)$

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 Δ ₁Filter with subscript "-" expression, large scale Δ ₂Filter with subscript "～" expression.Use Δ ₁, Δ ₂The subparticle stress that double filtration obtains can be than separately with the yardstick Δ ₂Filter 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:

$L_{\mathrm{ij}}={\widetilde{\mathrm{\τ}}}_{\mathrm{ij}}-\stackrel{~}{{\stackrel{\‾}{\mathrm{\τ}}}_{\mathrm{ij}}}---\left(6\right)$

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

$C_D=\frac{M_{\mathrm{ij}}{L}_{\mathrm{ij}}}{M_{\mathrm{ij}}{M}_{\mathrm{ij}}}---\left(8\right)$

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:

$I_{\mathrm{LM}}={\∫}_{t^{\′}=0}^{t^{\′}=t}{M}_{\mathrm{ij}}\left(t^{\′}\right)L_{\mathrm{ij}}\left(t^{\′}\right)W(t-t^{\′}){\mathrm{dt}}^{\′}---\left(9\right)$

$I_{\mathrm{MM}}={\∫}_{t^{\′}=0}^{t^{\′}=t}{M}_{\mathrm{ij}}\left(t^{\′}\right)M_{\mathrm{ij}}\left(t^{\′}\right)W(t-t^{\′}){\mathrm{dt}}^{\′}---\left(10\right)$

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:

$\frac{I_{\mathrm{LM}}^{n+1}-I_{\mathrm{LM}}^n}{\mathrm{\Δt}}=\frac{1}{T}[{\left(L_{\mathrm{ij}}{M}_{\mathrm{ij}}\right)}^{n+1}-I_{\mathrm{LM}}^{n+1}]---\left(11\right)$

$\frac{I_{\mathrm{MM}}^{n+1}-I_{\mathrm{MM}}^n}{\mathrm{\Δt}}=\frac{1}{T}[{\left(M_{\mathrm{ij}}{M}_{\mathrm{ij}}\right)}^{n+1}-I_{\mathrm{MM}}^{n+1}]---\left(12\right)$

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Δ ₂(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:

$C_D=C_s^2=\frac{I_{\mathrm{LM}}}{I_{\mathrm{MM}}}---\left(14\right)$

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:

$\stackrel{\widetilde{\‾}}{\mathrm{\Φ}}\left(x\right)=\frac{1}{\underset{j}{\mathrm{\Σ}}G\left(\right||x_j-x_i||,\mathrm{\Δ})}\underset{j}{\mathrm{\Σ}}{\stackrel{\‾}{\mathrm{\Φ}}}_{j}G\left(\right||x_j-x_i||,\mathrm{\Δ})---\left(15\right)$

Wherein, the expression formula of Gaussian function G (r, Δ) is as follows under two-dimensional case:

$G(r,\mathrm{\Δ})=\frac{9}{\mathrm{\π}{\mathrm{\Δ}}^2}\exp (-\frac{9r^2}{{\mathrm{\Δ}}^2})$

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 ³, kinematic viscosity is 0.00001m ²/ 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

$\mathrm{\γ}=\sqrt{2{\left(\frac{\∂u}{\∂x}\right)}^2+2{\left(\frac{\∂v}{\∂y}\right)}^2+{(\frac{\∂v}{\∂x}+\frac{\∂u}{\∂y})}^2}---\left(17\right)$

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:

$\frac{{\&PartialD;u}_i}{\&PartialD;x}=\underset{j\&NotEqual;i}{\mathrm{\&Sigma;}}\frac{x_j-x_i}{|r_j-r_i|}<{\&dtri;u}_i>=\frac{d}{n_0}\underset{j\&NotEqual;i}{\mathrm{\&Sigma;}}\left[\frac{x_j-x_i}{|r_j-r_i|}(u_j-u_i)\frac{(r_j-r_i)}{{|r_j-r_i|}^2}w\left(\right|r_j-r_i\left|\right)\right]---\left(18\right)$

$\frac{{\&PartialD;u}_i}{\&PartialD;y}=\underset{j\&NotEqual;i}{\mathrm{\&Sigma;}}\frac{y_j-y_i}{|r_j-r_i|}<{\&dtri;u}_i>=\frac{d}{n_0}\underset{j\&NotEqual;i}{\mathrm{\&Sigma;}}\left[\frac{y_j-y_i}{|r_j-r_i|}(u_j-u_i)\frac{(r_j-r_i)}{{|r_j-r_i|}^2}w\left(\right|r_j-r_i\left|\right)\right]---\left(19\right)$

$\frac{{\&PartialD;v}_i}{\&PartialD;x}=\underset{j\&NotEqual;i}{\mathrm{\&Sigma;}}\frac{x_j-x_i}{|r_j-r_i|}<{\&dtri;v}_i>=\frac{d}{n_0}\underset{j\&NotEqual;i}{\mathrm{\&Sigma;}}\left[\frac{x_j-x_i}{|r_j-r_i|}(v_j-v_i)\frac{(r_j-r_i)}{{|r_j-r_i|}^2}w\left(\right|r_j-r_i\left|\right)\right]---\left(20\right)$

$\frac{{\&PartialD;v}_i}{\&PartialD;y}=\underset{j\&NotEqual;i}{\mathrm{\&Sigma;}}\frac{y_j-y_i}{|r_j-r_i|}<{\&dtri;v}_i>=\frac{d}{n_0}\underset{j\&NotEqual;i}{\mathrm{\&Sigma;}}\left[\frac{y_j-y_i}{|r_j-r_i|}(v_j-v_i)\frac{(r_j-r_i)}{{|r_j-r_i|}^2}w\left(\right|r_j-r_i\left|\right)\right]---\left(21\right)$

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 ₀Represent 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:

$\mathrm{\&rho;}\frac{\mathrm{Du}}{\mathrm{Dt}}=-\&dtri;p+{\mathrm{\&mu;}\&dtri;}^{2}u+F---\left(22\right)$

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:

$u_i^{*}=u_i^n+\frac{\mathrm{\&Delta;t}}{\mathrm{\&rho;}}(\mathrm{\&mu;}\left(\mathrm{\&gamma;}\right){\&dtri;}^{2}u+F)---\left(23\right)$

$r_i^{*}=r_i^n+\mathrm{\&Delta;t}\&CenterDot;u_i^{*}---\left(24\right)$

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:

$\frac{\&PartialD;u}{\&PartialD;t}=-\frac{1}{\mathrm{\&rho;}}\&dtri;p+\frac{1}{\mathrm{\&rho;}}\&dtri;\&CenterDot;\mathrm{\&tau;}+f---\left(25\right)$

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:

$w\left(r\right)=\left\{\begin{array}{cc}\frac{r_e}{r}-1& (0<\&le;r_e)\\ 0& (r_e<r)\end{array}\right.---\left(26\right)$

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).

$w\left(r\right)=\left\{\begin{array}{cc}\frac{10}{{7\mathrm{\&pi;h}}^2}(1-\frac{3}{2}{q}^2+\frac{3}{4}{q}^3)& (0<q<1)\\ \frac{10}{28\mathrm{\&pi;}{h}^2}& (1\&le;q\&le;2)\\ 0& (q\&GreaterEqual;2)\end{array}\right.---\left(27\right)$

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)

${<\frac{1}{\mathrm{\&rho;}}\&dtri;\&CenterDot;\mathrm{\&tau;}>}_i=\underset{j}{\overset{N}{\mathrm{\&Sigma;}}}{m}_j[\frac{{\mathrm{\&tau;}}_i}{{\mathrm{\&rho;}}_i^2}+\frac{{\mathrm{\&tau;}}_j}{{\mathrm{\&rho;}}_j^2}]\&CenterDot;{\&dtri;}_i{w}_{\mathrm{ij}}---\left(28\right)$

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

${\mathrm{\&rho;}}_i=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{m}_j{w}_{\mathrm{ij}}---\left(29\right)$

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:

${<\&dtri;\&CenterDot;\mathrm{\&tau;}>}_i=\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}({\mathrm{\&tau;}}_i+{\mathrm{\&tau;}}_j)\&CenterDot;{\&dtri;}_i{w}_{\mathrm{ij}}/\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{w}_{\mathrm{ij}}---\left(30\right)$

Wherein:

${\&dtri;}_i{w}_{\mathrm{ij}}=\mathrm{dw}/{\mathrm{dr}}_{\mathrm{ij}}\&CenterDot;[(x_i-x_j)i+(y_i-y_j)j]/\left|r_{\mathrm{ij}}\right|---\left(31\right)$

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:

$D=\left[\begin{array}{cc}\frac{\&PartialD;u}{\&PartialD;x}& \frac{1}{2}(\frac{\&PartialD;v}{\&PartialD;x}+\frac{\&PartialD;u}{\&PartialD;y})\\ \frac{1}{2}(\frac{\&PartialD;v}{\&PartialD;x}+\frac{\&PartialD;u}{\&PartialD;y})& \frac{\&PartialD;v}{\&PartialD;y}\end{array}\right]---\left(33\right)$

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 ₀Mud-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:

$\mathrm{\&tau;}=({\mathrm{\&mu;}}_{\&infin;}+\frac{{\mathrm{\&mu;}}_0-{\mathrm{\&mu;}}_{\&infin;}}{1+{\left(K\right|\mathrm{\&gamma;}\left|\right)}^m})\mathrm{\&gamma;}---\left(34\right)$

Simplifying method in the list of references [1], formula (34) can be write following form:

$\mathrm{\&mu;}=({1000\mathrm{\&mu;}}_B+\frac{{{1000\mathrm{\&mu;}}_B}^2}{{\mathrm{\&tau;}}_B}|\mathrm{\&gamma;}\left|\right)/(1+\frac{1000{\mathrm{\&mu;}}_B}{{\mathrm{\&tau;}}_B}|\mathrm{\&gamma;}\left|\right)---\left(35\right)$

Wherein, μ _BBe Ben-Hur surrender viscosity, τ _BBe the Ben-Hur yield stress.

Fluid density is 1200kg/m ³, yield stress is taken as 25.0Pa, and Ben-Hur surrender viscosity is 0.07Ns/m ²

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.

Claims (7)

1. the construction method of Free Surface flow model in the improved semi implicit algorithm, the method may further comprise the steps:

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.

2. the construction method of Free Surface flow model in described a kind of improved semi implicit algorithm according to claim 1 is characterized in that: in the described step 3), adopt the interim velocity field of the explicit calculating of particle method large eddy simulation turbulence model of introducing.

3. the construction method of Free Surface flow model in described a kind of improved semi implicit algorithm according to claim 2, it is characterized in that: the Large eddy simulation method of Eulerian mesh method simulation turbulent flow is coupled in the particle method, 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.

4. the construction method of Free Surface flow model in described a kind of improved semi implicit algorithm according to claim 3, it is characterized in that: 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:

${\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{ij}}={2v}_t{\stackrel{\&OverBar;}{S}}_{\mathrm{ij}}-\frac{2}{3}{\mathrm{k\&delta;}}_{\mathrm{ij}}---\left(3\right)$

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:

$v_t={\left(C_s\mathrm{\&Delta;}\right)}^2{\left(2{\stackrel{\&OverBar;}{S}}_{\mathrm{ij}}{\stackrel{\&OverBar;}{S}}_{\mathrm{ij}}\right)}^{1/2}---\left(4\right)$

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:

$L_{\mathrm{ij}}={\widetilde{\mathrm{\&tau;}}}_{\mathrm{ij}}-\stackrel{~}{{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_{\mathrm{ij}}}---\left(6\right)$

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

$C_D=\frac{M_{\mathrm{ij}}{L}_{\mathrm{ij}}}{M_{\mathrm{ij}}{M}_{\mathrm{ij}}}---\left(8\right)$

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:

$C_D=C_s^2=\frac{I_{\mathrm{LM}}}{I_{\mathrm{MM}}}---\left(14\right)$

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.

5. the construction method of Free Surface flow model in described a kind of improved semi implicit algorithm according to claim 1, it is characterized in that: in the interim velocity field of described explicit calculating, adopt becoming viscosity Newtonian fluid model processes constitutive equation and is the non-Newtonian fluid of μ=f (| γ |), μ is kinetic viscosity, and γ is shearing rate.

6. the construction method of Free Surface flow model in described a kind of improved semi implicit algorithm according to claim 5, it is characterized in that: non-Newtonian fluid is considered as becoming the viscosity Newtonian fluid, 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.

7. the construction method of Free Surface flow model in described a kind of improved semi implicit algorithm according to claim 1, it is characterized in that: in the interim velocity field of described explicit calculating, adopt the shear stress of the discrete non-Newtonian fluid of cubic spline kernel function of SPH method.

Images