A multiphase flow simulation method based on density-separated SPH model

Abstract

The present invention discloses a multiphase flow simulation method based on a density-separated SPH model, which belongs to the field of fluid mechanics. The method first traverses all particles and finds all particles j that interact with particle i; for fluid particles, only the influence of the same-phase particles on the density of a certain particle is considered, and the density increment of the particle is calculated using the continuity equation after discretization by the SPH method; considering the influence of all particles in the support domain on the velocity of a certain particle, the velocity increment of the particle is calculated using the Navier-Stokes equation after discretization by the SPH method; for boundary particles, Shepard interpolation is used to solve the pressure and velocity of the particle, and the same boundary particle should have multiple different velocities to participate in the calculation of fluid particles of different phases; the density and velocity of the particle are updated, and the solution of the next time step is entered. The present invention can ensure the stability of numerical calculations, thereby greatly improving the calculation efficiency and saving the calculation cost.

Description

Technical Field

[0001] The present invention belongs to the field of computational fluid dynamics, and in particular relates to a multiphase flow simulation method based on a density-separated SPH (Smoothed Particle Hydrodynamics) model. Background Art

[0002] Traditional numerical simulation methods mostly use Euler-type mesh methods, such as the finite difference method (FDM), the finite volume method (FVM), and the finite element method (FEM), to approximate the Navier-Stokes momentum equation based on irregular grid interpolation. However, traditional mesh methods require special methods to deal with free surface flow problems. For example, the finite element method (FEM) cannot effectively solve large fluid deformation problems and requires mesh adjustment. The finite volume method (FDM) needs to track the changes of the free surface, and even if the Lagrangian dynamic mesh is used, the mesh distortion problem may occur.

[0003] Compared with the mesh method, the meshless method has special advantages in simulating incompressible flow with free surface, boundary deformation and moving interface. As a typical representative of the meshless method, SPH has a wide range of applications, including elastic flow, multiphase flow, quasi-incompressible flow, gravity flow, porous media flow, heat conduction, impact simulation, heat transfer and mass flow.

[0004] However, when simulating problems involving propagation and exchange across media, the SPH multiphase model is more complex and challenging than the single-phase model due to the discontinuity of physical fields such as density and velocity at the phase interface. To this end, scholars have proposed and developed many methods for multiphase interface processing, such as density reformulation, surface tension model, interface force model, volume summation method, etc. However, the existing multiphase flow SPH model still has many shortcomings, such as low accuracy, poor stability, and long calculation time. There is an urgent need for a multiphase flow SPH model that can take into account both numerical stability and numerical accuracy. Summary of the invention

[0005] In response to the problem of "insufficient accuracy and long calculation time of SPH multiphase flow model" in the prior art, the present application proposes a multiphase flow simulation method based on the density-separated SPH model. By separating particles of different phases when solving the continuity equation, and using a corrected pressure value instead of the pressure of particles of different phases when solving the momentum equation, the calculation accuracy and stability are improved, allowing a larger CFL number to be used in multiphase flow calculations, thereby improving the calculation efficiency.

[0006] In order to solve the above technical problems, the technical solution adopted by the present invention is:

[0007] A multiphase flow simulation method based on a density-separated SPH model comprises the following steps:

[0008] 1) According to the distribution and specific size of each phase of the fluid at the beginning of the multiphase flow process to be simulated, different types of SPH particles are used to discretize the fluids of different phases and establish the corresponding multiphase flow numerical model; within a certain time step, all particles are traversed to find all particles j that interact with particle i,

[0009] 2) If particle i is a fluid particle, the density increment of particle i is calculated using the continuity equation discretized by the SPH method, considering only the contribution of virtual particles in the support domain or fluid particles in the same phase as particle i to the density increment of particle i.

[0010] 3) If particle i is a fluid particle, the velocity increment of particle i is calculated using the Navier-Stokes equation discretized by the SPH method, taking into account the contribution of fluid particles and virtual particles of all phases in the support domain to the velocity change of particle i. However, when calculating the contribution of fluid particle j of a different phase from particle i to the velocity change of particle i, the Shepard interpolation result of the pressure at the location of particle j is used instead of the pressure of particle j itself to participate in the calculation.

[0011] 4) If i belongs to a boundary particle, use Shepard interpolation to solve the pressure and velocity of particle i; in the pressure interpolation of particle i, all particles in the domain are supported to participate in the interpolation calculation:

[0012]

[0013] In the above formula, f represents the fluid particles, p f ,m f ,ρ f are the pressure, mass and density of the fluid particles, W if is the kernel function value determined by the distance between the boundary particle i and the fluid particle f;

[0014] The same boundary particle i should have multiple different velocities to participate in the calculation of fluid particles of different phases.

[0015]

[0016] Among them, u f is the velocity vector of the fluid particle f, χ i is the set of all particles belonging to this phase;

[0017] For example, when simulating a gas-liquid two-phase flow problem, the boundary particle i should have two different velocities u iair ,u iwater , which are used to participate in the calculation of gas and liquid phases respectively:

[0018]

[0019]

[0020] 5) For all fluid particles, after obtaining the particle velocity increment and density increment according to steps 2) and 3), update the particle density and velocity according to the time advancement format used in the calculation and enter the solution of the next time step.

[0021] According to a preferred embodiment of the present invention, in step 2), the basis for judging whether particles i and j interact is: the distance r between particles i and j ij Less than the radius R of the kernel function multiplied by the smoothing length h of particle i.

[0022] According to a preferred embodiment of the present invention, in step 2), the density increment of particle i is calculated specifically according to the following formula:

[0023]

[0024] Among them, i is the density of particle i, u i ,u j are the velocity vectors of particles i and j, respectively, is the kernel function gradient, m j ,ρ j are the mass and density of particle j, respectively, i is the set of all virtual particles in the support domain of particle i and the fluid particles in the same phase as particle i, and t represents time.

[0025] According to a preferred embodiment of the present invention, in step 2), the density increment of particle i can be solved by directly discretizing the continuity equation, or a density dissipation term can be added to the continuity equation for solution. Specifically, the continuity equation after adding the density dissipation term to the continuity equation is as follows:

[0026]

[0027]

[0028] In the formula, p background is the background pressure value used in multiphase flow calculations, ρ 0,i is the reference density of the particle, j' is a particle in the support domain of particle j, and W jj′ is the kernel function value determined by the distance between particle j and particle j', m j' ,ρ j' are the mass and density of particle j' respectively.

[0029] According to a preferred embodiment of the present invention, in step 3), the velocity increment of particle i is calculated by the following formula:

[0030]

[0031] in,

[0032] π ij =(u j -u i )·(x j -x i ) / |x ij | 2 ,

[0033]

[0034] In the above formula, Ω i ,Ω j are the set formed by all particles in the support domain of particle i and the set formed by all particles in the support domain of particle j, g is the gravitational acceleration vector, c is the speed of sound, ρ 0i is the reference density of particle i, p i ,p j is the pressure of particles i and j, j' is a particle in the support domain of particle j, W jj′ is the kernel function value determined by the distance between particle j and particle j', dim is the dimension of the simulation problem, usually two-dimensional or three-dimensional, μ i and μ j are the dynamic viscosities of particles i and j, respectively, and x i ,x j is the position vector of particle i and particle j, x ij =x i -x j , h is the smooth length of the particle.

[0035] According to a preferred embodiment of the present invention, in step 4), Shepard interpolation can be performed at the location of the boundary particles, or a mapping point method can be used to extend half the particle spacing from the boundary inward along the normal direction of the boundary as a mapping point, and the Shepard interpolation results of the pressure and velocity at the mapping point are used as the velocity and pressure values ​​of the boundary particles.

[0036] For any particle i, whether it is a fluid particle or a boundary particle, the pressure p i and density ρ i The weak compressibility equation should be satisfied i =c 2 (ρ i -ρ 0,i), that is: for all boundary particles, in step 4), after obtaining their pressure, their density can be obtained by the weak compressibility equation; for all fluid particles, in step 5), after obtaining their density, their pressure can be obtained by the weak compressibility equation.

[0037] The beneficial effects of the present application are as follows: when solving the interaction between particles of different phases at the interface of multiphase flow, only the influence of particles of the same phase on the density of a certain particle is considered, and it is assumed that particles of different phases have no influence on the density, so that the numerical model is more in line with physical reality and the model accuracy is improved; when considering the influence of particles of different phases on the velocity, the pressure of particles of different phases is replaced by the Shepard interpolation of the pressure at the location of the particles of different phases, thereby improving the numerical stability at the interface. The SPH multiphase flow format in the prior art often requires a very small CFL number (0.05-0.1) and time step to ensure numerical stability due to its low precision and poor numerical stability. The new separated SPH multiphase flow model proposed in the present application can ensure the stability of numerical calculations while ensuring numerical accuracy, and can be calculated with a larger CFL condition number (0.2), thereby greatly improving the calculation efficiency and saving the calculation cost. BRIEF DESCRIPTION OF THE DRAWINGS

[0038] Figure 1 For a central particle near the phase interface, (a) the schematic diagram of the set χ formed by all the particles in the same phase with it, and (b) the schematic diagram of the set Ω formed by all the particles in the support domain;

[0039] Figure 2 Schematic diagram of pressure and velocity interpolation of boundary particles near the phase interface, (a) schematic diagram of pressure interpolation, (b) schematic diagram of velocity interpolation;

[0040] Figure 3 Schematic diagram of the numerical wave tank model layout considering gas-liquid two-phase;

[0041] Figure 4 Diagram of particle arrangement in simulation;

[0042] Figure 5 Time series results of the wave height meter when the relative wave height H0 / d=0.1;

[0043] Figure 6 Time series results of the wave height meter when the relative wave height H0 / d=0.5.

[0044] Figure 7 Schematic diagram of solitary wave breaking simulation.

[0045] Figure 8 When studying the crushing indexes and crushing types at different slopes, the prediction results of this method are compared with the results of the empirical formula.

[0046] Fig. 9When the relative wave height H0 / d0=0.6 and the slope gradient s=1 / 10, the gas-liquid interface curve shape at different times and the gas velocity field distribution. (a) t(g / d0) 1 / 2 =12.6(b)t(g / d0) 1 / 2 =15.4(c)t(g / d0) 1 / 2 =18.2. DETAILED DESCRIPTION

[0047] The present invention is further described and illustrated below in conjunction with specific embodiments. The embodiments are merely exemplary of the present disclosure and do not define the scope of limitation. The technical features of each embodiment of the present invention may be combined accordingly without conflicting with each other.

[0048] Example 1: The following is an example of numerical wave tank wave generation considering gas-liquid two-phase flow to illustrate the technical solution of the present invention. The simulation includes the following steps:

[0049] According to the attached Figure 3 The calculation model shown in the figure is used to establish a numerical water tank and initialize the various parameters required for the simulation. In the numerical simulation, the distance from the wave-making plate to the monitoring position of the wave height is x = 0.4m, the water depth is d = 0.114m, the relative wave height is H0 / d = 0.1,0.5, and the gravity acceleration is g = 9.8m / s 2 Particle spacing Δ = 0.002m, 0.003m, two different particle spacings are used for any calculation condition, and the simulation is performed twice to verify the calculation convergence; the density of gas particles ρ a =1kg / m 3 , artificial viscosity coefficient α a =0.001, pressure value p a =p background , that is, the pressure gradient of the gas is ignored; the reference density of the liquid particles ρ w =1000kg / m 3 , artificial viscosity coefficient α w =0.01, pressure value p w =p background +ρ w gz, z is the distance from the particle to the still water surface, and the density of the liquid particle is given by the weak compressibility equation, that is, The mass of all particles m = ρΔ 2 , smooth length h = 1.33Δ. Figure 4 The particle arrangement in the water tank at the initial moment is given, and the pressure coefficient P = p / ρ is also shown w Distribution of gd in the tank.

[0050] 1) After starting the calculation, the calculation steps for each step are:

[0051] 1a) In a certain time step, traverse all particles and find all particles j that interact with particle i. The basis for judging whether particles i and j interact is: the distance r between particles i and j ij Less than the radius R of the kernel function multiplied by the smooth length h of particle i;

[0052] 1b) If i is a fluid particle, the density increment of particle i is calculated using the continuity equation discretized by the SPH method, but only the contribution of virtual particles in the support domain or fluid particles in the same phase as particle i to the density increment of particle i is considered; the velocity increment of particle i is calculated using the Navier-Stokes equation discretized by the SPH method, and the contribution of fluid particles and virtual particles of all phases in the support domain to the velocity change of particle i is considered. However, when calculating the contribution of fluid particle j in a different phase from particle i to the velocity change of particle i, the Shepard interpolation result of the pressure at the location of particle j is used instead of the pressure of particle j itself to participate in the calculation;

[0053] 1c) If i belongs to a boundary particle, use Shepard interpolation to solve the pressure and velocity of particle i. 1d) For all fluid particles, after obtaining the velocity increment and density increment of the particle according to step 2b), update the density and velocity of the particle according to the time advancement format used for calculation (such as prediction-correction time format) and enter the solution of the next time step.

[0054] 2) After the physical time required for calculation is reached, stop the calculation and output the result.

[0055] Attached Figure 5 The time series results of the wave height meter when the relative wave height H0 / d = 0.1 are given. Figure 6 The time series results of the wave height meter when the relative wave height H0 / d=0.5 are given, and the results at different particle spacings Δ are also given. The predicted results from the GN theory and the experimental results are given for comparison. It can be seen from the figure that the separated multiphase flow SPH model proposed in this application has good numerical convergence and numerical accuracy.

[0056] Example 2: When a solitary wave propagates near the shore and encounters a slope, the wave will gradually deform and break in shallow water due to changes in the local water depth during the climbing process. During the breaking process, due to the drastic deformation of the phase interface and the drastic mixing between the gas-liquid two-phase, the existing multiphase flow numerical simulation methods are often difficult to simulate accurately, or require a very small time step to ensure numerical stability during the simulation process. The following takes the solitary wave breaking process of gas-liquid two-phase flow as an example to illustrate the technical solution of the present invention. The simulation includes the following steps:

[0057] According to the attached Figure 7The calculation model shown is used to establish a numerical flume and initialize the parameters required for the simulation. In the numerical simulation, the wave-making plate is 1.0m away from the slope foot, the water depth is d = 0.2m, the relative wave height H0 / d = 0.3-0.6, and different relative wave heights are used for simulation. The slope gradient s = 1 / 8-1 / 20, and different slopes are used for simulation. Gravitational acceleration g = 9.8m / s 2 , particle distance Δ = 0.01m, gas particle density ρ a =1kg / m 3 , artificial viscosity coefficient α a =0.001, pressure value p a =p background , that is, the pressure gradient of the gas is ignored; the reference density of the liquid particles ρ w =1000kg / m 3 , artificial viscosity coefficient α w =0.01, pressure value p w =p background +ρ w gz, z is the distance from the particle to the still water surface, and the density of the liquid particle is given by the weak compressibility equation, that is, The mass of all particles m = ρΔ 2 , smooth length h = 1.33Δ.

[0058] 1) After starting the calculation, the calculation steps for each step are:

[0059] 1a) In a certain time step, traverse all particles and find all particles j that interact with particle i. The basis for judging whether particles i and j interact is: the distance r between particles i and j ij Less than the radius R of the kernel function multiplied by the smooth length h of particle i;

[0060] 1b) If i is a fluid particle, the density increment of particle i is calculated using the continuity equation discretized by the SPH method, but only the contribution of virtual particles in the support domain or fluid particles in the same phase as particle i to the density increment of particle i is considered; the velocity increment of particle i is calculated using the Navier-Stokes equation discretized by the SPH method, and the contribution of fluid particles and virtual particles of all phases in the support domain to the velocity change of particle i is considered. However, when calculating the contribution of fluid particle j in a different phase from particle i to the velocity change of particle i, the Shepard interpolation result of the pressure at the location of particle j is used instead of the pressure of particle j itself to participate in the calculation;

[0061] 1c) If i belongs to a boundary particle, use Shepard interpolation to solve the pressure and velocity of particle i. 1d) For all fluid particles, after obtaining the velocity increment and density increment of the particle according to step 2b), update the density and velocity of the particle according to the time advancement format used for calculation (such as prediction-correction time format) and enter the solution of the next time step.

[0062] 2) After the physical time required for calculation is reached, stop the calculation and output the result.

[0063] Generally speaking, if the slope is gentler or the height of the solitary wave is higher, that is, when the slope gradient s decreases and the relative wave height H0 / d= increases, the breaking morphology of the solitary wave will change from no breaking (hereinafter referred to as NB) to surging (hereinafter referred to as SU), then to plunging (hereinafter referred to as PL), and then to spilling (hereinafter referred to as SP), until it becomes no breaking (NB) again when there is almost no slope. In order to comprehensively consider the influence of the slope gradient and the relative wave height of the solitary wave on the breaking morphology when the solitary wave climbs, some scholars have proposed a dimensionless constant To replace the original slope, S0 and crushing index H are fitted. b / d b The empirical formula for the time interval is given, and the specific breaking forms of waves at different S0 are pointed out, as shown in the attached figure. Figure 8 At the same time, Figure 8 The prediction results obtained by the multiphase flow simulation method based on the density separation SPH model proposed in this application are given. It can be seen that the prediction results obtained by this method are basically consistent with the empirical formula and have good accuracy. Fig. 9 The wave shape at different times during propagation is given when the slope is s = 1 / 10 and the relative wave height is H0 / d = 0.6. The curve of the gas-liquid interface at different times is drawn with a dotted line, and the vortex formed by the gas at the wave crest is shown at the same time, and the exact coordinates of the vortex center are given. It can be seen that even in multiphase flow simulations involving severe deformation of the phase interface and mixing of the two-phase fluid, this method can still predict the multiphase flow process well.

[0064] The above-mentioned embodiments only express several implementation methods of the present invention, and the description is relatively specific and detailed, but it cannot be understood as limiting the scope of the present invention. For ordinary technicians in this field, several modifications and improvements can be made without departing from the concept of the present invention, which all belong to the protection scope of the present invention.

Claims

1. A multiphase flow simulation method based on a density-separated SPH model, characterized in that: The steps include: 1) According to the distribution and specific size of each phase of the fluid at the beginning of the multiphase flow process to be simulated, different types of SPH particles are used to discretize the fluids of different phases and establish the corresponding multiphase flow numerical model; within a certain time step, all particles are traversed to find all particles j that interact with particle i, 2) If particle i is a fluid particle, the density increment of particle i is calculated using the continuity equation discretized by the SPH method, considering only the contribution of virtual particles in the support domain or fluid particles in the same phase as particle i to the density increment of particle i. 3) If particle i is a fluid particle, the velocity increment of particle i is calculated using the Navier-Stokes equation discretized by the SPH method, taking into account the contribution of fluid particles and virtual particles of all phases in the support domain to the velocity change of particle i. However, when calculating the contribution of fluid particle j of a different phase from particle i to the velocity change of particle i, the Shepard interpolation result of the pressure at the location of particle j is used instead of the pressure of particle j itself to participate in the calculation. 4) If i belongs to a boundary particle, use Shepard interpolation to solve the pressure and velocity of particle i; in the pressure interpolation of particle i, all particles in the domain are supported to participate in the interpolation calculation: In the above formula, f represents the fluid particles, p f ,m f ,ρ f are the pressure, mass and density of the fluid particles, W if is the kernel function value determined by the distance between the boundary particle i and the fluid particle f; The same boundary particle i should have multiple different velocities to participate in the calculation of fluid particles of different phases. Among them, u f is the velocity vector of the fluid particle f, χ i is the set of all particles belonging to this phase; 5) For all fluid particles, after obtaining the particle velocity increment and density increment according to steps 2) and 3), update the particle density and velocity according to the time advancement format used in the calculation and enter the solution of the next time step.

2. The multiphase flow simulation method based on the density separation SPH model according to claim 1 is characterized in that: In step 2), the basis for judging whether particles i and j interact is: the distance r between particles i and j ij Less than the radius R of the kernel function multiplied by the smoothing length h of particle i.

3. The multiphase flow simulation method based on the density separation SPH model according to claim 1 is characterized in that: In the step 2), the density increment of particle i is calculated specifically according to the following formula: Among them, i is the density of particle i, u i ,u j are the velocity vectors of particles i and j, respectively, is the kernel function gradient, m j ,ρ j are the mass and density of particle j, respectively, i is the set of all virtual particles in the support domain of particle i and the fluid particles in the same phase as particle i, and t represents time.

4. The multiphase flow simulation method based on the density separation SPH model according to claim 3 is characterized in that: In step 2), the continuity equation can be directly discretized to solve the density increment of particle i, or a density dissipation term can be added to the continuity equation to solve it. Specifically, the continuity equation after adding the density dissipation term to the continuity equation is as follows: In the formula, p background is the background pressure value used in multiphase flow calculations, ρ 0,i is the reference density of the particle, j' is a particle in the support domain of particle j, and W jj′ is the kernel function value determined by the distance between particle j and particle j', m j' ,ρ j' are the mass and density of particle j' respectively.

5. The multiphase flow simulation method based on the density separation SPH model according to claim 1 is characterized in that: In step 3), the velocity increment of particle i is calculated specifically by the following formula: in, π ij =(u j -u i )·(x j -x i ) / |x ij | 2 , In the above formula, Ω i ,Ω j are the set formed by all particles in the support domain of particle i and the set formed by all particles in the support domain of particle j, g is the gravitational acceleration vector, c is the speed of sound, ρ 0,i is the reference density of particle i, p i ,p j is the pressure of particles i and j, j' is a particle in the support domain of particle j, W jj′ is the kernel function value determined by the distance between particle j and particle j', dim is the dimension of the simulation problem, usually two-dimensional or three-dimensional, μ i and μ j are the dynamic viscosities of particles i and j, respectively, and x i ,x j is the position vector of particle i and particle j, x ij =x i -x j , h is the smooth length of the particle.

6. The multiphase flow simulation method based on the density separation SPH model according to claim 1 is characterized in that: In step 4), Shepard interpolation can be performed at the location of the boundary particles, or the mapping point method can be used to extend half the particle spacing from the boundary inward along the normal direction of the boundary as the mapping point, and the Shepard interpolation results of the pressure and velocity at the mapping point are used as the velocity and pressure values ​​of the boundary particles.

7. The multiphase flow simulation method based on the density separation SPH model according to claim 1 is characterized in that: For any particle i, whether it is a fluid particle or a boundary particle, the pressure p i and density ρ i The weak compressibility equation should be satisfied i =c 2 (ρ i -ρ 0,i ), that is: for all boundary particles, in step 4), after obtaining their pressure, their density can be obtained by the weak compressibility equation; for all fluid particles, in step 5), after obtaining their density, their pressure can be obtained by the weak compressibility equation.

Images