CN102867094B - The construction process of free surface flow model in a kind of improved semi implicit algorithm - Google Patents

Abstract

The present invention provides the construction process of free surface flow model in a kind of improved semi implicit algorithm: the turbulence model of introducing comprises the dynamic Smagorinsky model of static Smagorinsky model and Lagrangian Form; Adopt and become viscosity Newtonian fuid models treated rheological equation shape such as the non-Newtonian fluid of ��=f (| �� |); Introduce cubic spline kernel function, adopt the shear-stress of the discrete non-Newtonian fluid of divergence discrete solution of smoothed particle method method, process non-Newtonian fluid free surface flowing, the present invention can consider the sub-particle stress model of turbulent flow, and may be used for the calculating of non-Newtonian fluid.

Description

Method for constructing free surface flow model in moving particle semi-implicit algorithm

Technical Field

The invention relates to the field of numerical simulation of turbulence and non-Newtonian fluid free surface flow, and relates to construction and optimization of a turbulence and non-Newtonian fluid numerical model.

Background

In ship hydrodynamics, aerospace, hydraulic engineering and petroleum engineering, a large number of free surface flow problems exist. For example: the large carrier-liquid ship must be designed to take into account the influence of the sloshing of the liquid inside on the running stability of the ship, and the sloshing of the liquid in the cabin is a typical free surface flow problem. Also, as a natural disaster of debris flow caused by dam break, the development of debris flow is also a typical free surface flow problem. The common processes of evaporation and condensation, ink-jet printing, chemical reactions, etc. in daily life are accompanied by free surface flow phenomena. Because the free surface flow phenomenon is so common and so important, its research has been one of the hot spots in the field of fluid flow research.

With the development of computer technology, numerical simulation is widely applied in the fields of physics, chemistry, biology, aerospace and the like, and is considered as a third scientific research means in addition to theoretical analysis and experimental observation since the last century. Compared with experimental observation, the numerical simulation has the advantages of high efficiency, low cost, easy reproduction and the like. Since free surface flow mostly involves high-cost, large-scale or micro-scale study objects, experimental studies are difficult to perform, and numerical simulation is also one of the most powerful tools for studying free surface flow.

The most notable traditional grid-like methods used In the field of free surface flow numerical simulation research are the PIC (Particle-In-Cell) method, the mac (markerandcell) method, the vof (volume of fluid) method, and the LevelSet method. Although many advances have been made in dealing with the free surface flow problem, they all suffer from difficulties such as complex algorithms and large computational complexity in terms of irregular interface grid arrangement and accurate tracking or reconstruction of interfaces.

Particle-based methods are algorithms that have been rapidly developed in recent years, and since they use tracking of particle motion under the lagrangian framework, the position of the free surface can be tracked, and thus it is easy to simulate free surface flow. Smooth particle dynamics (SPH) and moving particle semi-implicit algorithms (MPS) are two common particle-based methods, SPH methods explicitly handle pressure terms using equation of state, so the fluid is approximately incompressible, while MPS implicitly handles pressure terms, the fluid is truly incompressible; the discrete scheme of the differential operator of the SPH method directly acts on the kernel function, so that the requirement on the kernel function is high, while the differential operator of the MPS has low requirement on the kernel function through particle weighted average dispersion, and the form is simple and rough.

In addition, the free surface flow common in nature and engineering fields is mostly the turbulent flow of non-newtonian fluid, such as the problem of dam break, the water flow of the dam break carries a large amount of mud, the non-newtonian fluid is required for modeling, and the flow is complex turbulent flow. Therefore, when the particle-based method is applied to the flow, a turbulence model and a non-newtonian fluid model need to be introduced in the particle method. At present, the MPS method is mainly applied to the flow problem of Newtonian fluid with smaller Reynolds number, and has less research on non-Newtonian and turbulent flow models.

Disclosure of Invention

The invention aims to provide a method for constructing a free surface flow model in a moving particle semi-implicit algorithm.

In order to achieve the purpose, the invention adopts the following technical scheme:

1) measuring physical parameters of a flow field and fluid, wherein the physical parameters comprise the size of the flow field, the density and the viscosity of the fluid;

2) setting an initial distribution of particles;

3) after the step 2), explicitly calculating a temporary speed field;

4) after the step 3), estimating the particle position;

5) after the step 4), calculating the particle number density, and then establishing a pressure Poisson equation;

6) after the step 5), solving a pressure gradient item of a pressure Poisson equation according to the pressure field, and then correcting the particle speed and the position.

In the step 3), a particle method large vortex simulation turbulence model (sub-particle stress model) is introduced to explicitly calculate the temporary velocity field, wherein the large vortex simulation turbulence model comprises a static Smagorinsky model or a dynamic Smagorinsky model in a Lagrange form.

A large vortex simulation method for simulating turbulence by using an Euler grid method is coupled to a particle method, a control equation of the large vortex simulation has one sub-particle stress term more than that of a conventional control equation, and the sub-particle stress term of the large vortex simulation turbulence model is constructed by adopting a static Smagorinsky model with a constant model coefficient or a dynamic Smagorinsky model in a Lagrange form.

The particle method large eddy simulation turbulence model explicit calculation temporary velocity field comprises the following steps:

1) when the Smagorinsky model was used, the sub-particle 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})$ is the deformation ratio tensor of the solvable scale,_ijrepresenting the Kronik tensor, v_tAnd k represents the vortex viscosity coefficient and the turbulence kinetic energy, respectively; the turbulence energy term is solved in conjunction with the pressure term, and the vortex viscosity coefficient is calculated by the following equation:

$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 the filtration scale, C_sIs the model coefficient;

2) by C_s=0.15 as model coefficients of a sub-particle stress model, called static Smagorinsky model; or calculating model coefficients at all positions according to a Germano equation by adopting a method of filtering twice continuously, and calling the model coefficients as a dynamic Smagorinsky model; the Germano equation is as follows:

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

whereinRepresenting the sub-particle stress after large-scale filtration,representing the sub-particle stress after two consecutive filtrations of the size scale;

3) for the dynamic Smagorinsky model, assuming that the coefficient of the dynamic Smagorinsky model is independent of the filtering scale, the Smagorinsky model described in the formulas (3) and (4) is substituted into the Germano equation, so that the following relation can be obtained:

definition ofAnd isThe dynamic Smagorinsky model coefficient of each point is

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

Wherein,representing the strain rate tensor, L, after large-scale filtering_kkIndicates the sub-particle stress increment L_ijA diagonal component of (a);

4) dynamic Smagorinsky model coefficients are determined using a time-averaging method of tracking particles, defining a time-averaged M over a particle track_ijL_ijAnd M_ijM_ijAre respectively I_LMAnd I_MM(ii) a The dynamic Smagorinsky model coefficients are as follows:

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

5) after obtaining the model coefficients, the sub-particle stress model is discretized by adopting the combination of the gradient model and the divergence model in the MPS method, after the discretization of the sub-particle stress model is completed, an explicit algorithm is adopted to solve the sub-particle stress term, and the first derivative of the velocity used in the sub-particle stress model is based on the velocity distribution of the previous time layer.

In the explicit calculation temporary velocity field, a variable viscosity Newtonian fluid model is adopted to process the non-Newtonian fluid with the constitutive equation of mu-f (| gamma |), mu is dynamic viscosity, and gamma is shear rate.

Aiming at the problem of free surface flow of non-Newtonian fluid with the constitutive equation of mu-f (| gamma |), the non-Newtonian fluid is regarded as variable-viscosity Newtonian fluid to be approximated, an MPS method suitable for the problem of flow of the non-Newtonian fluid is formed, in each time step, the shear rate gamma of each position of the flow field is calculated by using the speed field of the previous time step, and then the dynamic viscosity mu of each position of the flow field is calculated by using the constitutive equation mu-f (| gamma |), which is equivalent to converting the non-Newtonian fluid into the variable-viscosity Newtonian fluid, and the flow control equation of the Newtonian fluid is still suitable for the non-Newtonian fluid, wherein the difference is that mu is variable.

Aiming at the problem of free surface flow of the non-Newtonian fluid with the constitutive equation of mu-f (| gamma |), the kernel function commonly used by the MPS method can be abandoned, the cubic spline kernel function of the SPH method which can not only reflect the rule of particle interaction but also has correct spatial derivative is adopted, and then the shear stress of the non-Newtonian fluid is dispersed by adopting the dispersion scheme of the SPH, so that the MPS method suitable for the flow problem of the non-Newtonian fluid is formed.

The invention mainly aims at a moving particle semi-implicit algorithm (MPS), introduces a turbulence model of sub-particle stress and two non-Newtonian fluid models into a particle method, so that the particle method can well simulate engineering flow. The invention provides an MPS method based on the large vortex simulation and an MPS method aiming at non-Newtonian fluid with the equation of mu-f (| gamma |), and expands the application range of the MPS method.

Drawings

FIG. 1 is a flow chart of the algorithm of the present invention;

FIG. 2 is a geometry for a dam break problem;

FIG. 3 is a comparison of dam break problem simulation results and experimental results;

FIG. 4 is a schematic view of an exemplary Cross fluid dam break;

FIG. 5 is a graph showing the results of dam break simulation of Cross fluid using the non-Newtonian fluid MPS method approximated by a variable viscosity Newtonian fluid; (a) t =0.1s, (b) t =0.3s, (c) t =0.6 s; t represents the break development time;

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

fig. 7 is a Cross fluid dam free surface distribution at t =0.37s (fig. 7 a) and t =0.6s (fig. 7 b); t represents the break development time.

Detailed Description

The invention will be further explained with reference to the drawings.

The overall calculation flow of the algorithm of the present invention comprises the following steps (see fig. 1):

1. measuring physical parameters of a flow field and fluid, including the size of the flow field, the density and the viscosity of the fluid;

2. setting initial arrangement of particles according to an initial flow field;

3. the method for constructing the turbulence model and the non-Newtonian fluid model based on the moving particle semi-implicit algorithm is adopted to simulate the fluid flow;

4. and outputting the state of the particle model at each time step to a computer screen, and displaying the real-time flow process of the flow field.

The following describes some modifications in the flow. Reference is made to Koshizuka S.and OkaY.moving-particulate semi-iterative method for constructing information on compressing effective particle, NuclearScienceand engineering1996;123: 421-.

Sub-particle stress model for particle method large eddy simulation

Fig. 1, steps 8, 9, 10 are steps for considering turbulence effects, wherein static and dynamic sub-particle stress models are used, and the two models and the discrete mode are described in detail below. For the N-S equation for an incompressible fluid, the LES control equation obtained by filtering in physical space is 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,is the filtered large scale velocity component. The right side of the equation after filtering is added withThis term is referred to as the Sub-particle stress term (SPS) in the particle method.

The filtering operation is performed on a mesh scale in a mesh-like method, and should be performed on a particle scale in a particle-like method. In the invention, a Smagorinsky model widely applied in a grid method is still adopted to model the sub-particle stress model, and two methods, namely a static model and a dynamic model, are adopted to calculate the model coefficient of the Smagorinsky model. The key to introducing a sub-particle stress model into the particle method is: (1) how to adopt the discrete model in the particle method to disperse the Smagorinsky model; (2) how to couple the discrete sub-particle model into the algorithm of MPS.

When the Smagorinsky model was used, the sub-particle 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,is the deformation ratio tensor of the solvable scale,_ijrepresenting the Kronik tensor, v_tAnd k represents the vortex viscosity coefficient and turbulence kinetic energy, respectively, which can both characterize the intensity of the turbulence. The turbulence energy term is coupled with the pressure term to solve for, and the vortex viscosity coefficient can be calculated by the following equation:

$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 the filtration scale, C_sAre the model coefficients.

The invention employs two methods to determine the Smagorinsky model coefficients: first, the widely used empirical constant C in the grid-like method can be used_s=0.15 as coefficients of a sub-particle stress model, such a scheme is called the static Smagorinsky model; another method is to calculate the model coefficients of each position according to the Germano equation by using two filtering methods, and such a scheme is called a dynamic Smagorinsky model. For simplicity, the small scale Δ₁The filtration is indicated by the superscript "-", the large scale Δ₂The filtration is indicated by the superscript "-". By a₁、Δ₂The sub-particle stress obtained by two consecutive filtrations will be larger than that obtained by the scale delta alone₂Filtering the excess stress of the sub-particles; l is used for increasing the sub-particle stress_ijExpressed, its expression is as follows:

the Germano equation is as follows:

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

assuming dynamic Smigorinsky model coefficients C_D(Note: for ease of description, the dynamic Smagorinsky model coefficients and the static Smagorinsky model coefficients are related as followsIndependent of the filtering scale, substituting the Smagorinsky model into the Germano equation can obtain the following relationship:

definition ofThe dynamic Smagorinsky model coefficient of each point is

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

From the above equation, the dynamic Smagorinsky model coefficients at any time and any location can be determined, but C is thus calculated_DThe variation may be large, resulting in unstable calculations, and thus the model coefficients are obtained by averaging. Since the MPS method can easily track the movement of particles, the invention adopts the time-averaging method for tracking particles proposed by Meneveauc, LundTS, Cambot WH, Alagrangen dynamic complex-scalemodelofulency, journal of fluid mechanics1996, 319:353-85 to determine the dynamic Smigorinsky model coefficient, and define the time-averaged M on the particle track_ijL_ijAnd M_ijM_ijThe following were used:

$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 W (t-t') is an exponential weight function, such that I_LMAnd I_MMThe solution can be solved by a discrete method as follows:

$\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)$

where T is the relaxation time at the current time, n represents the current number of iteration steps, and Meveneau proposes the following calculation formula:

T=1.5Δ₂(I_LMI_MM)^-1/8（13）

through the calculation, I at different time and different positions can be obtained_MMAnd I_LMThe dynamic Smagorinsky model coefficients are as follows:

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

after the dynamic Smagorinsky model coefficients are obtained, the coefficients are substituted into the Smagorinsky model (formula (3) and formula (4)), and the sub-particle stress of the dynamic Smagorinsky model can be calculated.

The differentiation terms in the SPS model are all first order derivatives, and may be discretized using a combination of gradient and divergence models in the MPS method. In particular, sub-particle stressThe model of (2) only involves the first derivative terms of the velocity in different coordinate directions; the gradient model can give a first partial derivative of a velocity component in each coordinate direction, so that the gradient discrete model of the MPS can be completely used for calculating the sub-particle stressAfter the sub-particle stress is obtained, the addition of the particle generated by the sub-particle stress term can be calculated according to the divergence operator in the MPS methodSpeed. In the dynamic Smagorinsky model, large-scale filtering on a resolvable scale is also used, and according to the characteristics of the MPS discrete model, the large-scale filtering operation can be dispersed by adopting a weighted average method 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 of the gaussian function G (r, Δ) in two dimensions is as follows:

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

where x represents the position coordinates of the particle and r is the distance between two particles.

After the discretization of the sub-particle stress model is completed, the invention adopts an explicit algorithm to solve the sub-particle stress term, namely, the first derivative of the velocity used in the model is based on the velocity distribution of the previous time layer.

The invention adopts the dam break model shown in figure 2 to calculate the parallel efficiency, and tests the effect of the turbulence model, a liquid column is arranged in a water tank with the height of 0.5m and the length of 1.0m, the height of the liquid column is H = 0.5m, the width of the liquid column is L =0.25m, the liquid is Newtonian fluid, and the density is 1000kg/m³Kinematic viscosity of 0.00001m²The problem of dam break based on Smagorinsky model static and dynamic sub-particle stress was calculated under the condition of 50 × 100 particle arrangement, and for comparison with Martin experimental results, the position Z of the front edge of the dam break was dimensionless with L and the time t was taken asThe results were plotted in a nondimensional manner against experimental results and results considering only molecular viscosity as shown in FIG. 3. As can be seen in fig. 3: only the difference between the molecular viscosity result and the experimental result is considered, and the situation shows that the position of the front edge of the simulated dam break is advanced from the experimental result at the same moment; consider the followingThe simulation result after the particle stress effect is better matched with the experimental result; meanwhile, it can be seen that in the later stage of calculation, the static Smagorinsky model is closer to the experimental result than the dynamic Smagorinsky model. This demonstrates that using the sub-particle stress model does improve the simulation results.

2 non-Newtonian fluid model approximated by variable viscosity Newtonian fluid

FIG. 1, steps 4, 5, and 7 illustrate a first method of treating a non-Newtonian fluid by approximating the non-Newtonian fluid as a variable viscosity Newtonian fluid. And in each time step, calculating the shear rate gamma of each position of the flow field by using the velocity field of the previous time step, and calculating the dynamic viscosity mu of each position of the flow field by using the constitutive equation mu-f (| gamma |). After the above treatment, which is equivalent to converting the non-newtonian fluid into a variable viscosity newtonian fluid, the flow control equation of the newtonian fluid is still applicable to the non-newtonian fluid with the constitutive equation of μ ═ f (| γ |), except that μ is a variable.

For the two-dimensional flow problem, the shear rate γ is expressed as

$\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)$

X and y respectively represent two coordinate directions of a two-dimensional coordinate, and u and v respectively represent the speeds of a speed vector in the two coordinate directions;

in the MPS method, the partial derivative of velocity can be found using the 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 is_i、x_j、y_i、y_jRespectively, the position of the particle in the coordinate direction, the index i representing the central particle, the index j representing the surrounding particles, r_i、r_jCoordinate vector representing particle, w represents kernel function, n₀Which represents the initial population density of the particles,representing a hamiltonian;

after the shear rate γ is obtained, the dynamic viscosity μ of the flow field at the current time level can be calculated by using the constitutive equation μ ═ f (| γ |). According to the newton fluid momentum equation:

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

wherein,representing the random derivative of velocity, ρ representing the particle density, p representing the pressure, u representing the velocity vector, and F representing the volumetric force;

in this case, equation (22) is applied as long as the non-newtonian fluid is regarded as a variable viscosity newtonian fluid, except that μ is changed from a constant to a variable. Thus, the explicit partial calculation of the MPS method for calculating non-Newtonian fluid flow 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,which represents the temporary velocity of the particle i,representing the provisional velocity of the nth time step of the particle i, at represents the time step,which represents the temporary position of the particle i,a temporary position representing the nth time step of particle i;

3 non-Newtonian model using cubic spline kernel function

Fig. 1, steps 4, 5, and 6 show a second processing method for a non-newtonian fluid, because the kernel function form of the MPS method is simple and rough, which results in a narrow application range of a discrete format and a poor expansibility, and for the problem of free surface flow of a non-newtonian fluid whose constitutive equation is μ ═ f (| γ |), the present invention may also adopt a cubic spline kernel function which is commonly used in the SPH method and is capable of not using a kernel function commonly used in the MPS method, but reflecting the rule of particle interaction and having a correct spatial derivative, and further adopt a dispersion scheme of the dispersion of the SPH to disperse the dispersion of the shear stress τ of the non-newtonian fluid, thereby simulating the free surface flow of the non-newtonian fluid.

For non-newtonian fluids, the momentum equation is written as:

$\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 a shear stress tensor, and F = F/ρ;

the interactions between particles in both MPS and SPH methods are performed by kernel functions, but the kernel function requirements of the two methods are different. In the MPS method, the kernel function is only a weight function, and therefore the kernel function only needs to reflect the rule of the strength of the interaction between the particles. The most commonly used MPS method kernel function is shown by the following formula:

$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 is_eFor particles to effectively influence the radius, r is the inter-particle distance.

Since the SPH method directly acts the dispersion of various operators on the kernel function, the kernel function not only reflects the rule of particle interaction, but also ensures that the spatial derivative is correct. In order to adopt the divergence operator dispersion scheme of the SPH method, the invention adopts the SPH common cubic spline kernel function shown as the formula (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 is the smoothing length. Since the fluid in the MPS method is a truly incompressible fluid, h is a constant value, h = r_e/2。

The gradient operator and the Laplace operator still adopt the traditional MPS method discrete format, and the dispersion of the SPH common dispersion operator discrete format discrete shear stress tensor tau is adopted as follows: (τ represents the shear stress tensor, which may be referred to simply as shear stress when the tensor property is not emphasized.)

${<\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 the mass of the particles, w_ijDenotes the kernel function of particle j to particle i, and N denotes the total number of particles around particle i. The expression for the particle density p is

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

Similarly, since the MPS method can be used to characterize a truly incompressible fluid, m, ρ in the above equation are constants, and the particle density is equal to the macroscopic density of the fluid. Substituting equation (29) for equation (28) yields the divergence operator discretization format herein:

${<\&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 is_ijIndicates the distance between the particle i and the particle j.

Note that τ is a second order tensor, which is related to the one-dimensional vectorThe result of the dot product of (d) is a one-dimensional vector.

Thus, the second MPS algorithm for calculating the free surface flow of a non-Newtonian fluid is organized as follows: in each time step, the shear rate gamma of each position of the flow field is calculated by using the velocity field of the previous time step, the dynamic viscosity mu of each position of the flow field is calculated by using the constitutive equation mu-f (| gamma |), then the shear stress tensor tau is calculated according to the formula (32), and the divergence of the shear stress tensor tau is dispersed by using the formula (30).

τ=2μ(γ)D（32）

Mu is viscosity and D is strain rate tensor. In the two-dimensional problem, the strain rate tensor D is represented 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)$

the two non-Newtonian fluid models proposed by the present invention were verified using the same mud-water mixture dam as the documents [1] (Shao, S.D. and E.Lo, Incompression PHMethod for modeling Newtonian and not-Newtonianflowwith the surface water resources, advanced water resources2003.26(7): 787-.

Example model dimensions as shown in FIG. 4, at a length of 2m and a slope of S₀A 0.1% water tank is filled with a mud-water mixture of depth 0.1 m. And document [1]]Keeping the same, consider it as Cross fluid, whose 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)$

the simplification method in reference [1], formula (34) can be written in the form of:

$\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, mu_BIs Bingham yield viscosity, τ_BIs the bingham yield stress.

The fluid density was 1200kg/m³The yield stress is 25.0Pa, and the Bingham yield viscosity is 0.07Ns/m²。

The invention adopts the initial particle arrangement consistent with the document [1], the initial particle distance is 0.01m, 2000 fluid particles and 1256 wall surface particles are uniformly arranged, the calculation results when t is 0.1s, 0.3s and 0.6s are shown in fig. 5 and fig. 6, the calculation results are well consistent with the document [1], and the correctness of two non-Newtonian fluid algorithms provided by the invention is proved. For further comparison, the particle distributions at t-0.37 s and 0.6s were compared, as shown in fig. 7. It can be seen that the simulation results of the two non-Newtonian fluid algorithms provided by the invention are well consistent with the experimental results and the simulation results of the SPH method. In contrast, the free surface of the front end of the break dam obtained by the non-Newtonian fluid algorithm adopting the cubic spline kernel function is best in accordance with the experimental result.

Claims (6)

1. A method for constructing a free surface flow model in a moving particle semi-implicit algorithm comprises the following steps:

1) measuring physical parameters of a flow field and fluid, wherein the physical parameters comprise the size of the flow field, the density and the viscosity of the fluid;

2) setting an initial distribution of particles;

3) after the step 2), explicitly calculating a temporary speed field;

in the explicit calculation temporary velocity field, the shear stress of the SPH cubic spline kernel function discrete non-Newtonian fluid is adopted, and in order to adopt the divergence operator discrete scheme of the SPH method, the SPH common cubic spline kernel function shown in the formula (27) is adopted

$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}{(2-q)}^3& (1\&le;q<2)\\ 0& (q\&GreaterEqual;2)\end{array}\right.---\left(27\right)$

Where q is r/h and h is the smooth length, h is a constant value since the fluid in the MPS method is a truly incompressible fluid, h is r_e/2，r_eEffectively influencing the radius for the particle;

the gradient operator and the laplacian operator still adopt the traditional MPS method discrete format, and the dispersion of the SPH common dispersion operator discrete format discrete shear stress tensor tau is adopted as follows:

wherein m represents the mass of the particles, w_ijDenotes the kernel function of particle j to particle i, N denotes the total number of particles surrounding particle i, and the particle density ρ_iIs expressed as

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

Similarly, since the MPS method can represent the characteristics of a true incompressible fluid, where m, ρ are constants and the particle density is equal to the fluid density, equation (29) is substituted for equation (28) to obtain the divergence operator discrete format:

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

wherein:

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

wherein r is_ijRepresents the distance between the particle i and the particle j; x is the number of_i、x_j、y_i、y_jRespectively indicate the positions of the particles in the coordinate direction, the subscript i indicates a central particle, and the subscript j indicates a peripheral particle;

4) after the step 3), estimating the particle position;

5) after the step 4), calculating the particle number density, and then establishing a pressure Poisson equation;

6) after the step 5), solving a pressure gradient item of a pressure Poisson equation according to the pressure field, and then correcting the particle speed and the position.

2. The method for constructing the free surface flow model in the moving particle semi-implicit algorithm according to claim 1, wherein the method comprises the following steps: in the step 3), a temporary velocity field is explicitly calculated by adopting an introduced particle method large eddy simulation turbulence model.

3. The method for constructing the free surface flow model in the moving particle semi-implicit algorithm according to claim 2, wherein the method comprises the following steps: the large vortex simulation method for simulating turbulence by using the Euler grid method is coupled to the particle method, and a sub-particle stress term of the large vortex simulation turbulence model is constructed by adopting a static Smagorinsky model with a constant model coefficient or a dynamic Smagorinsky model in a Lagrange form.

4. The method for constructing the free surface flow model in the moving particle semi-implicit algorithm according to claim 3, wherein the method comprises the following steps: the particle method large eddy simulation turbulence model explicit calculation temporary velocity field comprises the following steps:

1) when the Smagorinsky model was used, the sub-particle stress model is as follows:

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

wherein,is the deformation ratio tensor of the solvable scale,_ijrepresenting the Kronik tensor, v_tAnd k represents the vortex viscosity coefficient and the turbulence kinetic energy, respectively; the turbulence energy term is solved in conjunction with the pressure term, and the vortex viscosity coefficient is calculated by the following equation:

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

wherein Δ represents the filtration scale, C_sIs the model coefficient;

2) by C_s0.15 as a model coefficient of the sub-particle stress model, called static Smagorinsky model; or calculating model coefficients at all positions according to a Germano equation by adopting a method of filtering twice continuously, and calling the model coefficients as a dynamic Smagorinsky model; the Germano equation is as follows:

whereinRepresenting the sub-particle stress after large-scale filtration,ruler for indicating sizeSub-particle stress after twice continuous filtration;

3) for the dynamic Smagorinsky model, assuming that the coefficient of the dynamic Smagorinsky model is independent of the filtering scale, the Smagorinsky model described in the formulas (3) and (4) is substituted into the Germano equation, so that the following relation can be obtained:

definition ofAnd isThe dynamic Smagorinsky model coefficient of each point is

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

Wherein,representing the strain rate tensor, L, after large-scale filtering_kkIndicates the sub-particle stress increment L_ijDiagonal line ofA component;

4) dynamic Smagorinsky model coefficients are determined using a time-averaging method of tracking particles, defining a time-averaged M over a particle track_ijL_ijAnd M_ijM_ijAre respectively I_LMAnd I_MM(ii) a The dynamic Smagorinsky model coefficients are as follows:

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

5) after obtaining the model coefficients, the sub-particle stress model is discretized by adopting the combination of the gradient model and the divergence model in the MPS method, after the discretization of the sub-particle stress model is completed, an explicit algorithm is adopted to solve the sub-particle stress term, and the first derivative of the velocity used in the sub-particle stress model is based on the velocity distribution of the previous time layer.

5. The method for constructing the free surface flow model in the moving particle semi-implicit algorithm according to claim 1, wherein the method comprises the following steps: in the explicit calculation temporary velocity field, a variable viscosity Newtonian fluid model is adopted to process the non-Newtonian fluid with the constitutive equation of mu-f (| gamma |), mu is dynamic viscosity, and gamma is shear rate.

6. The method for constructing the free surface flow model in the moving particle semi-implicit algorithm according to claim 5, wherein the method comprises the following steps: and (3) regarding the non-Newtonian fluid as variable-viscosity Newtonian fluid, calculating the shear rate gamma of each position of the flow field by using the velocity field of the previous time step in each time step, and calculating the dynamic viscosity mu of each position of the flow field by using the constitutive equation mu-f (| gamma |).