CN111241742A - Multiphase flow calculation method - Google Patents

Multiphase flow calculation method Download PDF

Info

Publication number
CN111241742A
CN111241742A CN201911377984.4A CN201911377984A CN111241742A CN 111241742 A CN111241742 A CN 111241742A CN 201911377984 A CN201911377984 A CN 201911377984A CN 111241742 A CN111241742 A CN 111241742A
Authority
CN
China
Prior art keywords
grid
particle
calculation
follows
calculated
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201911377984.4A
Other languages
Chinese (zh)
Other versions
CN111241742B (en
Inventor
张斌
张盛丰
刘晓星
单建强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN201911377984.4A priority Critical patent/CN111241742B/en
Publication of CN111241742A publication Critical patent/CN111241742A/en
Application granted granted Critical
Publication of CN111241742B publication Critical patent/CN111241742B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Abstract

The invention discloses a multiphase flow calculation method, which comprises the steps of firstly determining a calculation area, and inputting physical property parameters and corresponding boundary conditions; solving a convection equation of an Euler part by using a semi-Lagrange method of constraint interpolation conservation, solving a pressure Poisson equation of a non-convection part of a grid part, updating a numerical value of a next time step of the grid part, transmitting the speed of the grid part to particles for correction to obtain a new particle position, obtaining a grid color function value near an interface according to the particle position, and completing the numerical simulation calculation of multiphase flow if the simulation time is reached. The method has the accuracy of tracking the phase interface by the particle method and the rapidity of the grid method. The method is suitable for multiphase flow motion under the condition of not comprising heat transfer, wherein the advantages can be played most when the phase interface changes violently and the calculation area is large.

Description

Multiphase flow calculation method
Technical Field
The invention belongs to the technical field of fluid mechanics, and particularly relates to a multiphase flow calculation method.
Background
The multiphase flow is a complex flow phenomenon, is widely applied to the industrial fields of petroleum, chemical industry, metallurgy, power, energy and the like, and the research on the multiphase flow has very important significance for the technological progress of related industries. The capture of phase interfaces is a key and difficult point in numerical simulation of multiphase flow. In order to solve the problem of capturing phase interfaces in the conventional euler method, many phase interface Tracking and reconstructing methods such as a VOF method, a Level Set method, a Front-Tracking method, a VOSET method and the like are developed. However, these methods lose the details of the phase interface during the reconstruction or construction of the phase interface, and require complicated meshing operations, which is inconvenient in the calculation process for processing multiphase flow. The lagrangian approach based on the lagrangian framework shows its particular advantages in numerical simulations of multiphase flows. These methods include the SPH method, the MPS method, and the FVP method. These lagrangian methods are used in the calculation of multiphase flow with severe interfacial changes such as bubble rise, dam break, vapor explosion, and solid-liquid coupling motion. While this class of methods can automatically track phase interfaces, no additional effort is required to capture or track the interface. However, these methods have too many calculation particles, which consume calculation time and calculation performance.
Disclosure of Invention
The technical problem to be solved by the present invention is to provide a multiphase flow calculation method, which can not only accurately track the phase interface of multiphase flow, but also reduce the consumption of calculation performance as much as possible and shorten the time required for calculation, in view of the above-mentioned deficiencies in the prior art.
The invention adopts the following technical scheme:
a multiphase flow calculation method comprises the steps of firstly determining a calculation area, and inputting physical property parameters and corresponding boundary conditions; solving a convection part of a control equation under an Euler coordinate system by using a semi-Lagrange method of constrained interpolation conservation, solving a pressure Poisson equation of a non-convection part under the Euler coordinate system, and updating a numerical value of a next time step of a grid part; and transferring the speed of the grid part to a particle calculation part for correction to obtain a new particle position, updating the color function value of the grid near the interface according to the calculation result of the particle part, and completing the simulation calculation of the multiphase flow numerical value after the simulation time is up.
Specifically, the governing equation for viscous fluid in the convection section is:
Figure BDA0002341493140000021
Figure BDA0002341493140000022
Figure BDA0002341493140000023
where p is the density of the particles,
Figure BDA0002341493140000024
is the velocity, p is the pressure, μ is the kinematic viscosity,
Figure BDA0002341493140000025
is the momentum and theta is the mesh volume fraction.
Further, the convection term equation is:
Figure BDA0002341493140000026
Figure BDA0002341493140000027
Figure BDA0002341493140000028
solving for next time step by using semi-Lagrange method of constraint interpolation conservation
Figure BDA0002341493140000029
And
Figure BDA00023414931400000210
the value of (a) is,
Figure BDA00023414931400000211
and
Figure BDA00023414931400000212
is the volume-integrated average and the area-integrated average of the grid density at the next time instant,
Figure BDA00023414931400000213
and
Figure BDA00023414931400000214
is the volume-integrated and area-integrated average of the temporal grid momentum.
Specifically, solving the poisson equation of the pressure of the non-convective part in the euler coordinate system specifically includes:
the semi-discrete non-convective term equation of the grid section viscous fluid control equation is:
Figure BDA0002341493140000031
the non-convective part is calculated as follows:
Figure BDA0002341493140000032
the poisson equation for pressure is calculated as follows:
Figure BDA0002341493140000033
wherein p isn+1Is the pressure value of the next moment, u□Temporary velocities, p, calculated for the convection partn+1Is the grid density of the next time instant, m□Is the temporary momentum calculated for the convection part.
Further, updating the momentum area integral average and the volume integral average of the grid at the next moment specifically includes:
the momentum area integral of the next time-of-day grid averages as:
Figure BDA0002341493140000034
wherein the content of the first and second substances,
Figure BDA0002341493140000035
is the pressure volume integral average of the i +1 grid at the next time instant,
Figure BDA0002341493140000036
is the pressure volume integral average of the i grid at the previous time instant;
the grid volume integral of the next time step of momentum averages:
Figure BDA0002341493140000037
wherein the content of the first and second substances,
Figure BDA0002341493140000038
and
Figure BDA0002341493140000039
the momentum volume integral average of the i-grid at the next and temporal instants, respectively.
Specifically, the particle calculation part is modified specifically as follows:
in the particle calculation of the particle fraction, the mass equation and the momentum conservation equation for a viscous incompressible fluid are as follows:
Figure BDA0002341493140000041
Figure BDA0002341493140000042
where ρ is the fluid density, t is the time,
Figure BDA0002341493140000043
p is the pressure, upsilon is the kinematic viscosity coefficient,
Figure BDA0002341493140000044
in order to be a surface tension force,
Figure BDA0002341493140000045
is the acceleration of gravity;
in the two-dimensional finite volume method, a control volume of moving particles is a circle whose discrete area is calculated as follows:
V＝πR2＝l2
where V is the discrete volume, R is the control volume radius of the particle, and l is the size of the grid.
Further, according to the relationship between the area of the curve and the volume fraction in the gaussian theorem, in the two-dimensional finite volume particle method, the gradient operator and the laplace operator are calculated as follows:
Figure BDA0002341493140000046
Figure BDA0002341493140000047
wherein phi issIs the physical quantity acting between particle i and particle j,
Figure BDA0002341493140000048
is the area of action of particle j on particle i,
Figure BDA0002341493140000049
is the distance between particles i and j, the area of action of j on particle i
Figure BDA00023414931400000410
The calculation is as follows:
Figure BDA00023414931400000411
wherein n is0Is the initial population density of the particles,
Figure BDA00023414931400000412
is the unit vector between particle i and particle j, wijIs the sum function value between particles i and j;
initial kernel function integral value n0The calculation is as follows:
Figure BDA00023414931400000413
unit vector between two particles
Figure BDA00023414931400000414
The calculation is as follows:
Figure BDA0002341493140000051
kernel function wijThe calculation is as follows:
Figure BDA0002341493140000052
wherein r iseTo cut off the radius, among
Figure BDA0002341493140000053
When wij＝0；
The gradient operator and laplacian operator are calculated as follows:
Figure BDA0002341493140000054
Figure BDA0002341493140000055
wherein phi issCan be expressed as a linear function calculation as follows:
Figure BDA0002341493140000056
。
further, the initial particle velocity of the particle fraction is calculated as follows: the temporary velocity and position of the particle are obtained by interpolating the central velocity of the grid calculated from the grid part by using the kernel function, and are calculated as follows:
Figure BDA0002341493140000057
Figure BDA0002341493140000058
wherein the content of the first and second substances,
Figure BDA0002341493140000059
for the initial velocity of the interpolation of the i-particles,
Figure BDA00023414931400000510
is the center velocity, w, of the grid (m, n)i,(m,n)Is the weight of i particle to grid (m, n), the weight of i particle to grid (m, n) is calculated as follows:
Figure BDA00023414931400000511
wherein r iseIs the cutoff radius.
Further, the particle of the particle fraction is calculated by correcting the velocity and position according to the finite volume particle method, and the divergence of the velocity for the incompressible fluid is calculated as follows:
Figure BDA0002341493140000061
the relationship between the correction value of the velocity and the gradient of the pressure is determined according to the momentum equation of the particle motion as follows:
Figure BDA0002341493140000062
pressure poisson equation calculation:
Figure BDA0002341493140000063
the pressure Poisson equation is discretely solved by using the gradient operator and the Laplace operator to obtain the pressure, and the relationship between the corrected particle velocity and the particle pressure is as follows:
Figure BDA0002341493140000064
the position of the particle at the new instant is calculated as follows:
Figure BDA0002341493140000065
。
specifically, the color function value of the grid near the interface in the grid calculation is calculated according to the number density of the particles in the grid, after the particle part is calculated, the color function value of the grid at the interface is updated according to the positions of the particles, and the color function value θ of the grid is:
Figure BDA0002341493140000066
Figure BDA0002341493140000067
Figure BDA0002341493140000068
wherein n is0The particle number density, r, of particles when the grid is filled with particles of 4X 4eThe cutoff radius is a value of 8.1 times the initial distance of the particle.
Compared with the prior art, the invention has at least the following beneficial effects:
the invention provides a mixing method for calculating multiphase flow, which is based on an internal interpolation/multi-moment finite volume method in an Euler method and a finite volume particle method in a Lagrange method, combines the requirement on phase interface tracking in multiphase flow calculation, and arranges Lagrange particles near a phase interface for calculation, thereby being capable of accurately tracking the phase interface, and the particles are only arranged near the phase interface, thereby reducing the number of particles required by calculation, shortening the calculation time and reducing the consumption of performance.
Furthermore, the difficulty of capturing the phase interface is calculated by combining the multiphase flow in the Euler method and the inconvenience of grid adjustment near the phase interface, so that the capturing of the phase interface in the Euler method is avoided.
Furthermore, the particles of the Lagrangian are only arranged near the phase interface, so that the number of the particles required by calculation of the Lagrangian method is reduced, and the calculation time is shortened.
Furthermore, the Lagrange method and the Euler method are combined, the advantages of the Lagrange method and the Euler method are fully combined, the advantage that the phase interface of the Lagrange method is accurately captured is fully utilized, and the advantage that the Euler method is fast in calculation is utilized.
In summary, the invention combines the advantages of the mesh method and the lagrangian method, and tracks the movement of the phase interface by calculating the movement of the particles without a specific interface tracking algorithm. Provides a rapid and accurate multiphase flow numerical simulation calculation method.
The technical solution of the present invention is further described in detail by the accompanying drawings and embodiments.
Drawings
FIG. 1 is a flow chart of the present invention;
FIG. 2 is a definition of grid volume integral averaging and area integral;
FIG. 3 is a discrete region of a finite volume particle method;
FIG. 4 is a schematic diagram of the mesh and particle velocity transfer;
FIG. 5 is a calculated area division of a dam break;
FIG. 6 shows the simulation calculation result of dam break in the hybrid method;
FIG. 7 is a distance time calculation;
FIG. 8 is a height time calculation;
FIG. 9 is a region of settlement in a single bubble ascent;
fig. 10 is a calculation result of the present method for a single bubble rise.
Detailed Description
The invention provides a multiphase flow calculation method, which is based on an internal interpolation/multi-moment finite volume method in an Euler method and a finite volume particle method in a Lagrange method, combines the difficulty of the Euler method in actual numerical simulation calculation for capturing a multiphase flow phase interface and the time consumption of the Lagrange method calculation, obviously reduces the particle number in the Lagrange method by arranging Lagrange particles in a plurality of layers of grids near the phase interface, and avoids the capture of the phase interface in the complex Euler method. The hybrid method combines the advantages of the grid method and the Lagrange method, and the movement of the phase interface is tracked by calculating the movement of the particles by the moving particle method, so the phase interface is automatically determined by the distribution of the particles calculated by the particle method, a specific forward tracking algorithm is not needed, the rest parts are calculated by the grid method, the numerical diffusion cannot be generated, and the calculation time is greatly shortened.
Referring to fig. 1, a multiphase flow calculation method according to the present invention uses a method of mixing a grid method and a moving particle method to calculate multiphase flow. In the mixed method, a grid part is solved by using an internal constraint interpolation/multi-moment finite volume method, and a particle part is solved by using a finite volume particle method. The method comprises the steps of dividing a control equation of viscous fluid into a convection term and a non-convection term by an internal constraint interpolation/multi-moment finite volume method used by a grid part of a mixing method, performing distribution solution, transferring the central speed of a grid to a particle calculation part of the mixing method, then performing calculation of the particle part to correct the position of particles, and obtaining the value of a color function of the grid according to the number density of the particles in the grid near an interface after correction. In the mixing method, the particles are only arranged near the phase interface of the multiphase flow, so that the number of the Lagrangian particles required by calculation can be reduced, the calculation time is shortened, and meanwhile, the advantages of the interface can be tracked by using the Lagrangian method without any special method for tracking the interface. In the numerical simulation calculation, most calculation areas are calculated by using a grid method, only the calculation is carried out by using a moving particle method near a phase interface of the multiphase flow, and after the particle calculation of the Lagrange part is finished, the value of a color function in the grid method is corrected according to the particle number density at the phase interface, so that the quick and accurate calculation of the multiphase flow is realized, and the method comprises the following steps:
s1, determining the calculation region, and the corresponding physical property and boundary condition of the calculation region
Determining a calculation region, and calculating region boundaries of different phases in the region. Giving physical parameters in the calculation regions of the respective grid sections, including: density, viscosity, surface tension coefficient. Lagrangian particles are arranged near the phase interface according to the division of the grid part calculation region. And after the Lagrange particles are arranged, calculating to obtain the value of the initial color function of each grid according to the initial arrangement of the particles.
S2, solving convection part calculation of a control equation under an Euler coordinate system;
in the mixing method, all calculation regions are calculated by using an internal interpolation/multi-moment finite volume method in an Euler method, then a finite volume particle method in a Lagrange method is used for calculating a grid region near a phase interface, the color function values of grids near the phase interface in the grid method are updated according to the particle arrangement obtained by the particle method, and the Lagrange particles clearly describe the boundary of the phase interface.
The internal interpolation/multi-moment finite volume method of Euler part in the mixing method divides the control equation of viscous flow into two parts, one part is convection part, and the other part is non-convection part. The calculation for the convection part is as follows:
the control equation for viscous fluids in the grid method is as follows:
Figure BDA0002341493140000101
Figure BDA0002341493140000102
Figure BDA0002341493140000103
where p is the density of the particles,
Figure BDA0002341493140000104
is the velocity, p is the pressure, μ is the kinematic viscosity,
Figure BDA0002341493140000105
is the momentum and theta is the mesh volume fraction.
Wherein the convection term equation set is as follows:
Figure BDA0002341493140000106
Figure BDA0002341493140000107
Figure BDA0002341493140000108
where p is the density of the particles,
Figure BDA0002341493140000109
is the velocity, p is the pressure, μ is the kinematic viscosity,
Figure BDA00023414931400001010
is the momentum and theta is the mesh volume fraction.
Solving the equation of convection by using a semi-Lagrange method of constraint interpolation conservation to obtain the next time step
Figure BDA00023414931400001011
And
Figure BDA00023414931400001012
the value of (c). Wherein the content of the first and second substances,
Figure BDA00023414931400001013
and
Figure BDA00023414931400001014
is the volume-integrated average and the area-integrated average of the grid density at the next time instant,
Figure BDA00023414931400001015
and
Figure BDA00023414931400001016
is the volume-integrated and area-integrated average of the temporal grid momentum.
Further, the calculation of the convection is followed by the calculation of the viscous force, the surface tension and the effect of gravity on the grid.
The solution to the flow equation uses a cubic constrained interpolation function-the semi-lagrange method. The method is to establish an interpolation function according to the volume integral average and the area integral average of the grid, as shown in fig. 2.
Wherein the grid volume integral mean is defined as:
Figure BDA00023414931400001017
area-integrated average definition of grid:
Figure BDA0002341493140000111
Figure BDA0002341493140000112
wherein, Δ xiAnd Δ yiRespectively the size of the grid in the x-direction and the y-direction, phi (x, y, t) is the value of some physical quantity of the grid (x, y) at the time t, phi (x)i+1/2Y, t) is the value of some physical quantity of the right boundary in the x-direction of the grid (x, y), phi (x, y)j+1/2And t) is a value of some physical quantity of the upper interface in the grid (x, y) y direction.
The method is to establish an upwind interpolation function according to the definition of the grid volume integral average and the area integral average, and in the x direction, under the condition that u is larger than 0, the construction function based on the left side interpolation is as follows:
Figure BDA0002341493140000113
wherein x ∈ [ x ]i-(1/2),xi+(1/2)]。
The four coefficients can be determined according to the continuity of the left and right interfaces of the grid physical quantity, the conservation of the grid volume and the slope of the grid center:
according to the left and right sides of the grid, the following components are continuously arranged:
Figure BDA0002341493140000114
Figure BDA0002341493140000115
wherein
Figure BDA0002341493140000116
Is the area integral average of the left interface in the x direction,
Figure BDA0002341493140000117
is the area integral average of the right interface in the x-direction.
The conservation according to the grid volume integral is:
Figure BDA0002341493140000118
wherein
Figure BDA0002341493140000119
Is the volume-integrated average in the x-direction.
According to the middle phi of the gridi LThe first derivative of (a) is:
Figure BDA0002341493140000121
wherein the slope diAn approximate representation of the fourth order is taken, of the form:
Figure BDA0002341493140000122
further, four coefficients can be calculated according to the four definitions as follows:
Figure BDA0002341493140000123
Figure BDA0002341493140000124
Figure BDA0002341493140000125
Figure BDA0002341493140000126
similarly, the right-side interpolation constructor phi can be obtained when u < 0i R。
Further, cubic interpolation functions on the left side and the right side are obtained by the same method, and a value (SIA) of the physical quantity area integral average of the interface position of the next time step can be obtained by using a half Lagrange method.
The semi-lagrange method used therein is:
Figure BDA0002341493140000127
wherein the content of the first and second substances,
Figure BDA0002341493140000128
the area integral average temporary value of grid (x, y) at the x right interface,
Figure BDA0002341493140000129
is based on the constructor on the right-hand side,
Figure BDA00023414931400001210
is based on the left constructor, where α has the value:
Figure BDA00023414931400001211
wherein the content of the first and second substances,
Figure BDA00023414931400001212
and
Figure BDA00023414931400001213
the values of (d) may be discrete as:
Figure BDA00023414931400001214
Figure BDA0002341493140000131
wherein u isi+1Is the x-direction velocity value of the i +1 grid in the x-direction. u. ofiIs the x-direction velocity value of the i-grid in the x-direction, ui+(1/2)Is the value of the velocity of the x-direction right interface of the i-grid in the x-direction.
Further, the correction value at the next time is:
Figure BDA0002341493140000132
where τ is the motion trajectory of the half-Lagrangian particle in the dt time interval. Solving for the second term in the equation for the next time step
Figure BDA0002341493140000133
The following approximation may be used:
Figure BDA0002341493140000134
wherein the content of the first and second substances,
Figure BDA0002341493140000135
is a temporary value of the speed after the half lagrange solution.
After the area integral average of the grid is obtained, the grid volume integral average value of the next time step is:
Figure BDA0002341493140000136
wherein, gi+(1/2),gi-(1/2)At the boundary position t of the corresponding gridn-1-tnThe flux in the chamber.
Figure BDA0002341493140000137
Figure BDA0002341493140000138
Wherein the content of the first and second substances,
Figure BDA0002341493140000139
is the grid volume integral average of the last instant.
The value of the next grid volume integral average is updated in this way, and the density of the grid volume integral average and the density of the area integral average at the next moment, and the temporary momentum of the grid volume integral average and the momentum of the area integral average can be obtained by the method.
S3, solving a non-convection part under an Euler coordinate system;
after the calculation of the non-convection part, a viscosity term and a source term need to be calculated, and after the calculation of the viscosity term and the source term, a semi-discrete non-convection part equation can be obtained, so that the semi-discrete non-convection part equation is as follows:
Figure BDA0002341493140000141
wherein m isn+1Is the momentum value of the next moment, m□Is a temporary momentum value, pn+1Is the pressure value at the next moment.
Further, obtaining:
Figure BDA0002341493140000142
further, a poisson equation for pressure is obtained:
Figure BDA0002341493140000143
wherein p isn+1Is the pressure value of the next moment, u□Temporary velocities, p, calculated for the convection partn+1Is the grid density of the next time instant, m□Is the temporary momentum calculated for the convection part.
The value of the pressure at the next moment can be obtained by solving the pressure poisson equation. And updating the momentum area integral average and the volume integral average of the grid at the next moment.
The momentum area integral of the next time-of-day grid averages as:
Figure BDA0002341493140000144
wherein the content of the first and second substances,
Figure BDA0002341493140000145
is the pressure volume integral average of the i +1 grid at the next time instant,
Figure BDA0002341493140000146
is the pressure volume integral average of the i-grid at the last instant.
The grid volume integral of the next time step of momentum averages:
Figure BDA0002341493140000147
wherein the content of the first and second substances,
Figure BDA0002341493140000151
and
Figure BDA0002341493140000152
the momentum volume integral average of the i-grid at the next and temporal instants, respectively.
And finishing the calculation of the grid part.
S4, transferring physical quantities of the grid part and the particle calculation part;
in the mixing method, the particles are distributed only on the side of the denser liquid in the phase interface. As shown in fig. 3, a grid is filled with 4 x 4 particles, and the particles are only distributed on one side of the incompressible fluid with a relatively high density in the two phases, and the thickness of the particle layer near the phase interface is generally greater than that of a four-layer grid.
After the grid part is calculated, the particles near the interface are calculated, so that the phase interface is tracked according to the movement of the particles, and the value of the color function near the phase interface is determined according to the density of the kernel function of the particles. As shown in fig. 4, the initial velocity of the particle is interpolated from the center of the grid according to the position of the grid computing part where the particle is located.
Figure BDA0002341493140000153
Wherein
Figure BDA0002341493140000154
For the initial velocity of the interpolation of the i-particles,
Figure BDA0002341493140000155
is the center velocity of the grid (m, n). w is ai,(m,n)Is the weight of the i particle and the grid (m, n).
Figure BDA0002341493140000156
Wherein r iseThe initial distance of the particle whose value is 8.1 times the cutoff radius.
S5, correcting and calculating the particle part;
for a viscous incompressible fluid, the mass equation and conservation of momentum equations are as follows:
Figure BDA0002341493140000157
Figure BDA0002341493140000158
where ρ is the fluid density, t is the time,
Figure BDA0002341493140000161
p is the pressure, upsilon is the kinematic viscosity coefficient,
Figure BDA0002341493140000162
in order to be a surface tension force,
Figure BDA0002341493140000163
is the acceleration of gravity.
In the two-dimensional finite volume method, a control volume of moving particles is a circle whose discrete area is calculated as follows:
V＝πR2＝l2
where V is the discrete volume, R is the control volume radius of the particle, and l is the size of the grid.
Further, according to the relationship between the area of the curve and the volume fraction in the gaussian theorem, in the finite volume particle method, the gradient operator and the laplacian operator are calculated as follows:
Figure BDA0002341493140000164
Figure BDA0002341493140000165
wherein phi issIs the physical quantity acting between particle i and particle j,
Figure BDA0002341493140000166
is the area of action of particle j on particle i,
Figure BDA0002341493140000167
is the distance between particles i and j.
In particular, the area of action of j on the particle i
Figure BDA0002341493140000168
The calculation is as follows:
Figure BDA0002341493140000169
wherein n is0Is the initial population density of the particles,
Figure BDA00023414931400001610
is the unit vector between particle i and particle j, wijIs the sum function value between particles i and j.
In particular, the initial kernel function integral value n0The calculation is as follows:
Figure BDA00023414931400001611
in particular, the unit vector between two particles
Figure BDA00023414931400001612
The calculation is as follows:
Figure BDA00023414931400001613
in particular, the kernel function wijThe calculation is as follows:
Figure BDA0002341493140000171
specifically, reTo cut off the radius, among
Figure BDA0002341493140000172
When wij＝0。
Further, the gradient operator and laplacian operator are calculated as follows:
Figure BDA0002341493140000173
Figure BDA0002341493140000174
where S is the area calculation region of the particle and has a value of 2 π R, V is the discrete volume of the particle, n0Is the initial particle number density, wijIs the kernel function weight value between particles i and j,
Figure BDA0002341493140000175
is the distance between particles i and j, phisAre physical quantities to be dispersed.
φsExpressed as a linear function, is calculated as follows:
Figure BDA0002341493140000176
the initial particle velocity of the particle fraction is calculated as follows:
specifically, the calculation of the temporary velocity and position of the particle by interpolating the central velocity of the grid calculated from the grid part using the kernel function is as follows:
Figure BDA0002341493140000177
Figure BDA0002341493140000178
wherein the content of the first and second substances,
Figure BDA0002341493140000179
for the initial velocity of the interpolation of the i-particles,
Figure BDA00023414931400001710
is the center velocity of the grid (m, n). w is ai,(m,n)Is a particle ofAnd the weight of the grid (m, n).
Wherein the weight of i particle and grid (m, n) is calculated as follows:
Figure BDA0002341493140000181
wherein r iseThe cutoff radius is the initial distance of the particle with a value of 8.1 times.
After interpolation to obtain the initial velocity of the particle, the temporary position of the particle motion can be obtained:
Figure BDA0002341493140000182
wherein the content of the first and second substances,
Figure BDA0002341493140000183
respectively.
After the initial position is obtained, the position of the particles at the interface position is adjusted using the finite volume particle method. For incompressible fluids, it can be known from the continuity equation that the divergence of the velocity needs to be 0.
Figure BDA0002341493140000184
The correction value for the speed is derived from the gradient of the corrected pressure:
Figure BDA0002341493140000185
the following pressure poisson equation can be obtained:
Figure BDA0002341493140000186
specifically, the relationship between the corrected particle velocity and the particle pressure is as follows:
Figure BDA0002341493140000187
the particle position is updated by solving the pressure poisson equation above.
Figure BDA0002341493140000188
And S6, updating the color function value of the grid near the interface according to the calculation result of the particle part, and completing the simulation calculation of the multiphase flow numerical value if the simulation time is reached.
After the particle part is calculated, the value of the color function of the grid at the interface is updated according to the position of the particle.
Figure BDA0002341493140000191
Where θ is the color function value of the grid, where,
Figure BDA0002341493140000192
n0is the particle number density of particles when the grid is full of particles according to 4X 4, where wi,(m,n)The calculation is as follows:
Figure BDA0002341493140000193
wherein r iseThe cutoff radius is the initial distance of the particle with a value of 8.1 times.
The method has the accuracy of tracking the phase interface by the particle method and the rapidity of the grid method. The method is suitable for multiphase flow motion under the condition of not comprising heat transfer, wherein the advantages can be played most when the phase interface changes violently and the calculation area is large.
In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. The components of the embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations. Thus, the following detailed description of the embodiments of the present invention, presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
In order to verify the correctness, the effectiveness and the rapidity of the developed mixing method for calculating the multiphase flow, a numerical simulation mode is adopted for verification, dam break and single bubble rising calculation examples are selected for simulation, and the simulation results are compared with simulation results of experiments or other methods.
Wherein, the grid number in the dam break calculation is 40 multiplied by 40, and the gravity acceleration in the calculation is 9.8m/s2. The size of the grid in the calculation region is 2 × 10-2m, particle size 5X 10-3m, the boundary condition is a no-slip boundary. Table 1 is a detailed parameter list of the dam break condition, and fig. 5 is the division of the calculation area of the dam break.
Table 1 list of parameters of dam break
Figure BDA0002341493140000201
The results of the simulated calculation of the break in the hybrid method are shown in fig. 6. The results of the intermediate moments of the calculation of the grid section of the dam break calculation are presented in fig. 6 (a). Fig. 6(b) shows the result of the intermediate time of the particle portion calculated by the mixing method. It can be seen that the liquid hits the right wall surface at around 0.3 s.
In the dam break calculation process, the moving distance of the liquid level and the height of the liquid level before the liquid collides with the right wall surface are recorded, the moving distance of the liquid level and the height of the liquid level are analyzed in a non-dimensionalization mode, and the non-dimensionalization time is taken as
Figure BDA0002341493140000202
The distance is taken as L/a, wherein g is the acceleration of gravity, a is the initial width of the liquid level, and L is the distance moved by the liquid level. As shown in fig. 7 and 8The moving distance of the liquid surface and the height of the liquid surface before the impact of about 0.3 second were compared with the results of the experiment. From the above figures it can be seen that the results of the calculations herein are close to the experimental results.
However, the grid number required by the calculation is 40 × 40, the method is far more accurate in capturing the phase interface than the VOF method, the method can accurately capture the interfaces of two phases, and the method is not like other Euler methods which need special phase interface capturing technologies, and in addition, because the particles are only arranged near the phase interface, the calculated particle number is saved by 48% to the maximum in the calculation process, so that the fast calculation compared with the Lagrangian method is realized.
Fig. 9 shows the settlement area in the single bubble rising. In the simulation of single bubble rising, the distance between the center of the bubble and the left and right wall surfaces needs at least 2.5 times of the initial radius width of the bubble, and the distance between the center of the bubble and the lower bottom needs at least 1.5 times of the initial radius height of the bubble, so that the influence of the left and right wall surfaces and the bottom wall surface on the movement of the bubble is avoided.
The calculated area herein is a 100 x 100 grid of areas. The size of the grid is 1 × 10-3m, the diameter of the corresponding bubble is 15mm, and the physical properties of the gas phase and the liquid phase are shown in Table 2.
TABLE 2 parameter List for Single bubble rise
Figure BDA0002341493140000211
The calculation results of the method for single bubble rise are shown in fig. 10. Where fig. 10(a) is a result of calculation of the grid part and fig. 10(b) is a result of calculation of the particle part. It can be seen that the final shape of the bubble is finally a spherical crown shape at the corresponding morton and aid numbers. This is in good agreement with the Grace empirical relationship under the corresponding Morton and Aids dimensions, and the results of the calculation are also substantially consistent compared with the results of other methods.
In the simulation of single bubble rising, the method has accurate phase interface capture, is more real than the VOF in the traditional Euler method, does not need a special phase interface capture technology, and if the Lagrangian method is used, 160000 particles are needed by the method, which is a huge particle number for the particle method, and the calculation takes several days, while in the method, the particles are only arranged near the phase interface, so the number of particles needed by the calculation is only 3208 at most. The number of particles required by calculation is greatly reduced, and the calculation of the embodiment can be completed only in a few hours.
In summary, the method of the present invention uses a method of mixing a grid method and a moving particle method to calculate multiphase flow. The lagrangian particles in the mixing method are only arranged near the phase interface of the multiphase flow, so that the number of the lagrangian particles required by calculation can be reduced, the calculation time is shortened, and the advantage of tracking the interface by using the lagrangian method without any special method is utilized to track the interface. Most calculation areas in the numerical calculation are calculated by using a grid method, only the mobile particle method is used for calculating near a phase interface of the multiphase flow, and after the particle calculation of the Lagrange part is finished, the value of a color function in the grid method is corrected according to the particle number density at the phase interface, so that the rapid and accurate calculation of the multiphase flow is realized.
The above-mentioned contents are only for illustrating the technical idea of the present invention, and the protection scope of the present invention is not limited thereby, and any modification made on the basis of the technical idea of the present invention falls within the protection scope of the claims of the present invention.

Claims (10)

1. A multiphase flow calculation method is characterized in that firstly, a calculation area is determined, and physical property parameters and corresponding boundary conditions are input; solving a convection part of a control equation under an Euler coordinate system by using a semi-Lagrange method of constrained interpolation conservation, solving a pressure Poisson equation of a non-convection part under the Euler coordinate system, and updating a numerical value of a next time step of a grid part; and transferring the speed of the grid part to a particle calculation part for correction to obtain a new particle position, updating the color function value of the grid near the interface according to the calculation result of the particle part, and completing the simulation calculation of the multiphase flow numerical value after the simulation time is up.
2. A multiphase flow calculation method as defined in claim 1 wherein the governing equation for the viscous fluid in the convection section is:
Figure FDA0002341493130000011
Figure FDA0002341493130000012
Figure FDA0002341493130000013
where p is the density of the particles,
Figure FDA0002341493130000014
is the velocity, p is the pressure, μ is the kinematic viscosity,
Figure FDA0002341493130000015
is the momentum and theta is the mesh volume fraction.
3. A multiphase flow calculation method as defined in claim 2 wherein the set of convection term equations is:
Figure FDA0002341493130000016
Figure FDA0002341493130000017
Figure FDA0002341493130000018
solving for next time step by using semi-Lagrange method of constraint interpolation conservation
Figure FDA0002341493130000019
And
Figure FDA00023414931300000110
the value of (a) is,
Figure FDA00023414931300000111
and
Figure FDA00023414931300000112
is the volume-integrated average and the area-integrated average of the grid density at the next time instant,
Figure FDA00023414931300000113
and
Figure FDA00023414931300000114
is the volume-integrated and area-integrated average of the temporal grid momentum.
4. The multiphase flow calculation method of claim 1, wherein solving the poisson equation of pressure of the non-convective part in the euler coordinate system is specifically as follows:
the semi-discrete non-convective term equation of the grid section viscous fluid control equation is:
Figure FDA0002341493130000021
the non-convective part is calculated as follows:
Figure FDA0002341493130000022
the poisson equation for pressure is calculated as follows:
Figure FDA0002341493130000023
wherein p isn+1Is as followsPressure value of one moment, u□Temporary velocities, p, calculated for the convection partn+1Is the grid density of the next time instant, m□Is the temporary momentum calculated for the convection part.
5. A multi-phase flow calculation method as claimed in claim 4, wherein updating the momentum area integral average and the volume integral average of the grid at the next time is embodied as:
the momentum area integral of the next time-of-day grid averages as:
Figure FDA0002341493130000024
wherein the content of the first and second substances,
Figure FDA0002341493130000025
is the pressure volume integral average of the i +1 grid at the next time instant,
Figure FDA0002341493130000026
is the pressure volume integral average of the i grid at the previous time instant;
the grid volume integral of the next time step of momentum averages:
Figure FDA0002341493130000027
wherein the content of the first and second substances,
Figure FDA0002341493130000028
and
Figure FDA0002341493130000029
the momentum volume integral average of the i-grid at the next and temporal instants, respectively.
6. The multiphase flow calculation method of claim 1, wherein the particle calculation section modification is specifically:
in the particle calculation of the particle fraction, the mass equation and the momentum conservation equation for a viscous incompressible fluid are as follows:
Figure FDA0002341493130000031
Figure FDA0002341493130000032
where ρ is the fluid density, t is the time,
Figure FDA0002341493130000033
p is the pressure, upsilon is the kinematic viscosity coefficient,
Figure FDA0002341493130000034
in order to be a surface tension force,
Figure FDA0002341493130000035
is the acceleration of gravity;
in the two-dimensional finite volume method, a control volume of moving particles is a circle whose discrete area is calculated as follows:
V＝πR2＝l2
where V is the discrete volume, R is the control volume radius of the particle, and l is the size of the grid.
7. A multiphase flow calculation method according to claim 6, wherein according to the relation of curve area and volume fraction in Gaussian theorem, in the two-dimensional finite volume particle method, the gradient operator and Laplace operator are calculated as follows:
Figure FDA0002341493130000036
Figure FDA0002341493130000037
wherein phi issIs the physical quantity acting between particle i and particle j,
Figure FDA0002341493130000038
is the area of action of particle j on particle i,
Figure FDA0002341493130000039
is the distance between particles i and j, the area of action of j on particle i
Figure FDA00023414931300000310
The calculation is as follows:
Figure FDA00023414931300000311
wherein n is0Is the initial population density of the particles,
Figure FDA00023414931300000312
is the unit vector between particle i and particle j, wijIs the sum function value between particles i and j;
initial kernel function integral value n0The calculation is as follows:
Figure FDA0002341493130000041
unit vector between two particles
Figure FDA0002341493130000042
The calculation is as follows:
Figure FDA0002341493130000043
kernel function wijThe calculation is as follows:
Figure FDA0002341493130000044
wherein r iseTo cut off the radius, among
Figure FDA0002341493130000045
When wij＝0；
The gradient operator and laplacian operator are calculated as follows:
Figure FDA0002341493130000046
Figure FDA0002341493130000047
wherein phi issCan be expressed as a linear function calculation as follows:
Figure FDA0002341493130000048
。
8. a multiphase flow calculation method according to claim 7, wherein the initial particle velocity of the particle fraction is calculated as follows: the temporary velocity and position of the particle are obtained by interpolating the central velocity of the grid calculated from the grid part by using the kernel function, and are calculated as follows:
Figure FDA0002341493130000049
Figure FDA00023414931300000410
wherein the content of the first and second substances,
Figure FDA00023414931300000411
for the initial velocity of the interpolation of the i-particles,
Figure FDA00023414931300000412
is the center velocity, w, of the grid (m, n)i,(m,n)Is the weight of i particle to grid (m, n), the weight of i particle to grid (m, n) is calculated as follows:
Figure FDA0002341493130000051
wherein r iseIs the cutoff radius.
9. A multiphase flow calculation method according to claim 6, wherein the particles of the particle fraction are calculated for velocity and position corrections according to a finite volume particle method, and for incompressible fluids, the divergence of velocity is calculated as follows:
Figure FDA0002341493130000052
the relationship between the correction value of the velocity and the gradient of the pressure is determined according to the momentum equation of the particle motion as follows:
Figure FDA0002341493130000053
pressure poisson equation calculation:
Figure FDA0002341493130000054
the pressure Poisson equation is discretely solved by using the gradient operator and the Laplace operator to obtain the pressure, and the relationship between the corrected particle velocity and the particle pressure is as follows:
Figure FDA0002341493130000055
the position of the particle at the new instant is calculated as follows:
Figure FDA0002341493130000056
。
10. a multiphase flow calculation method according to claim 1, wherein the grid color function value near the interface in the grid calculation is calculated according to the number density of particles in the grid, the particle part is calculated, the value of the grid color function at the interface is updated according to the positions of the particles, and the grid color function value θ is:
Figure FDA0002341493130000061
Figure FDA0002341493130000062
Figure FDA0002341493130000063
wherein n is0The particle number density, r, of particles when the grid is filled with particles of 4X 4eThe cutoff radius is a value of 8.1 times the initial distance of the particle.
CN201911377984.4A 2019-12-27 2019-12-27 Multiphase flow calculation method Active CN111241742B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911377984.4A CN111241742B (en) 2019-12-27 2019-12-27 Multiphase flow calculation method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911377984.4A CN111241742B (en) 2019-12-27 2019-12-27 Multiphase flow calculation method

Publications (2)

Publication Number Publication Date
CN111241742A true CN111241742A (en) 2020-06-05
CN111241742B CN111241742B (en) 2021-11-19

Family

ID=70877582

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911377984.4A Active CN111241742B (en) 2019-12-27 2019-12-27 Multiphase flow calculation method

Country Status (1)

Country Link
CN (1) CN111241742B (en)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111967149A (en) * 2020-08-03 2020-11-20 电子科技大学 Particle motion semi-interpolation solving method for particle simulation algorithm
CN113158520A (en) * 2021-04-09 2021-07-23 西安交通大学 Fuel ice layer interface tracking simulation method for freezing target system
CN113761812A (en) * 2021-09-09 2021-12-07 西安交通大学 Method and system for solving single-layer wall surface complex geometric watershed based on gridless method
CN113792496A (en) * 2021-08-20 2021-12-14 华南理工大学 Free surface tension modeling method based on combination of particles and grids
CN114357907A (en) * 2022-01-07 2022-04-15 中国空气动力研究与发展中心计算空气动力研究所 Parallel method suitable for Lagrange type particle numerical simulation
CN114757082A (en) * 2022-03-07 2022-07-15 同济大学 Fluid-solid coupling calculation method based on Lagrange-Euler stable distribution point
CN115655648A (en) * 2022-10-13 2023-01-31 哈尔滨工程大学 Bubble information acquisition system and method for measuring and inverting pulsation type bubble motion characteristics through flow field pressure

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010096204A (en) * 2008-10-14 2010-04-30 Mori Sachiko Transport pipe
CN102867094A (en) * 2012-09-19 2013-01-09 西安交通大学 Establishment method for free surface flow model in moving particle semi-implicit algorithm
CN103714575A (en) * 2013-12-30 2014-04-09 北京大学 Fluid simulation method for combining SPH and dynamic surface mesh
CN104268943A (en) * 2014-09-28 2015-01-07 北京航空航天大学 Fluid simulation method based on Eulerian-Lagrangian coupling method
KR20160129671A (en) * 2015-04-30 2016-11-09 주식회사 엘쏠텍 Prediction and assessment apparatus for pipe-thinning based on three-dimensional computational fluid dynamics and method thereof
CN107657131A (en) * 2017-10-18 2018-02-02 北方工业大学 Fluid interactive simulation method and system based on GPUs (general purpose computing) clusters
CN108491619A (en) * 2018-03-19 2018-09-04 浙江大学 Based on physics and the non-physical complex scene fluid structurecoupling efficient analogy method mixed
CN110147575A (en) * 2019-04-18 2019-08-20 河海大学 A kind of calculation method that the two-phase stream interface based on single layer particle levels collection captures
CN110287590A (en) * 2019-06-24 2019-09-27 天津大学 Based on operator splitting and the method for improving semi-Lagrange solution dispersion of pollutants

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010096204A (en) * 2008-10-14 2010-04-30 Mori Sachiko Transport pipe
CN102867094A (en) * 2012-09-19 2013-01-09 西安交通大学 Establishment method for free surface flow model in moving particle semi-implicit algorithm
CN103714575A (en) * 2013-12-30 2014-04-09 北京大学 Fluid simulation method for combining SPH and dynamic surface mesh
CN104268943A (en) * 2014-09-28 2015-01-07 北京航空航天大学 Fluid simulation method based on Eulerian-Lagrangian coupling method
KR20160129671A (en) * 2015-04-30 2016-11-09 주식회사 엘쏠텍 Prediction and assessment apparatus for pipe-thinning based on three-dimensional computational fluid dynamics and method thereof
CN107657131A (en) * 2017-10-18 2018-02-02 北方工业大学 Fluid interactive simulation method and system based on GPUs (general purpose computing) clusters
CN108491619A (en) * 2018-03-19 2018-09-04 浙江大学 Based on physics and the non-physical complex scene fluid structurecoupling efficient analogy method mixed
CN110147575A (en) * 2019-04-18 2019-08-20 河海大学 A kind of calculation method that the two-phase stream interface based on single layer particle levels collection captures
CN110287590A (en) * 2019-06-24 2019-09-27 天津大学 Based on operator splitting and the method for improving semi-Lagrange solution dispersion of pollutants

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JIE LIU 等: "A hybrid particle-mesh method for viscous,incompressible, multiphase flows", 《JOURNAL OF COMPUTATIONAL PHYSICS》 *
李勇霖 等: "对流传热问题的粒子-网格混合方法数值模拟", 《原子能科学技术》 *
王艺萍 等: "粒子法与网格法的耦合及在固液两相流中的应用", 《中国优秀博硕士学位论文全文数据库(硕士)工程科技Ⅱ辑》 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111967149A (en) * 2020-08-03 2020-11-20 电子科技大学 Particle motion semi-interpolation solving method for particle simulation algorithm
CN111967149B (en) * 2020-08-03 2022-11-04 电子科技大学 Particle motion semi-interpolation solving method for particle simulation algorithm
CN113158520A (en) * 2021-04-09 2021-07-23 西安交通大学 Fuel ice layer interface tracking simulation method for freezing target system
CN113158520B (en) * 2021-04-09 2022-10-28 西安交通大学 Fuel ice layer interface tracking simulation method for freezing target system
CN113792496A (en) * 2021-08-20 2021-12-14 华南理工大学 Free surface tension modeling method based on combination of particles and grids
CN113761812A (en) * 2021-09-09 2021-12-07 西安交通大学 Method and system for solving single-layer wall surface complex geometric watershed based on gridless method
CN113761812B (en) * 2021-09-09 2024-03-29 西安交通大学 Solving method and system for single-layer wall complex geometric drainage basin based on gridless method
CN114357907A (en) * 2022-01-07 2022-04-15 中国空气动力研究与发展中心计算空气动力研究所 Parallel method suitable for Lagrange type particle numerical simulation
CN114357907B (en) * 2022-01-07 2023-03-21 中国空气动力研究与发展中心计算空气动力研究所 Parallel method suitable for Lagrange type particle numerical simulation
CN114757082A (en) * 2022-03-07 2022-07-15 同济大学 Fluid-solid coupling calculation method based on Lagrange-Euler stable distribution point
CN114757082B (en) * 2022-03-07 2024-04-12 同济大学 Fluid-solid coupling calculation method based on Lagrangian-Euler stable distribution point
CN115655648A (en) * 2022-10-13 2023-01-31 哈尔滨工程大学 Bubble information acquisition system and method for measuring and inverting pulsation type bubble motion characteristics through flow field pressure

Also Published As

Publication number Publication date
CN111241742B (en) 2021-11-19

Similar Documents

Publication Publication Date Title
CN111241742B (en) Multiphase flow calculation method
Yanke et al. Simulation of slag-skin formation in electroslag remelting using a volume-of-fluid method
CN107220399A (en) Weight the whole flow field analogy method of non-oscillatory scheme substantially based on Hermite interpolation
CN105760602A (en) Total flow field numerical simulation method for finite volume weighted essentially non-oscillatory scheme
CN107423511A (en) Meet to immerse border implicit iterative solving method without sliding boundary condition and the condition of continuity
CN115587551B (en) Multi-scale prediction method for ablation behavior of heat-proof structure of hypersonic aircraft
CN112613243B (en) Method, device and computer readable storage medium for hydrodynamic simulation
CN115238611A (en) Multiphase flow simulation efficiency optimization method based on multiphase lattice Boltzmann flux method
CN110717269A (en) Fluid surface detail protection method based on grid and particle coupling
Essadki et al. Adaptive mesh refinement and high order geometrical moment method for the simulation of polydisperse evaporating sprays
WO2020244597A1 (en) Boltzmann-based method for simulating cvi densification process of composite material
CN105022928A (en) Digitized real-time determination method for center-of-gravity position of fuel system of aircraft
Sato et al. A comparative study of the cumulant lattice Boltzmann method in a single-phase free-surface model of violent flows
CN102063551B (en) Method for simulating segregation numerical value of blast furnace distributing charge level
CN108509741B (en) Finite element numerical solution method of debris flow equation
CN104749625A (en) Regularization technology based seismic data dig angle estimation method and device
Piña et al. Air injection in vertical water column: Experimental test and numerical simulation using volume of fluid and two-fluid methods
Maliska et al. An unstructured finite volume procedure for simulating flows with moving fronts
CN105183965B (en) For predicting the Large eddy simulation method of atomization process
Cao et al. Smoothed particle hydrodynamics modeling and simulation of foundry filling process
Wang et al. An improved free surface capturing method based on Cartesian cut cell mesh for water-entry and-exit problems
CN114611423A (en) Three-dimensional multiphase compressible fluid-solid coupling rapid calculation method
Sijoy et al. An Eulerian multi-material scheme for elastic–plastic impact and penetration problems involving large material deformations
CN104517299A (en) Method for restoring and resimulating physical video fluid driving model
Razavi et al. Numerical model of slug development on horizontal two-phase flow

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant