A general simulation method for complex particle-laden two-phase flow

Abstract

The present application belongs to the technical field of flow field, and particularly relates to a general simulation method of complex particle two-phase flow. The method comprises the following steps: S1: constructing a geometric representation model to represent the geometric information of the complex particle and dividing a grid; S2: on the basis of the grid divided in S1, adopting a characteristic line splitting finite element method to solve the flow field to obtain the flow field velocity at pseudo n+1 time and the flow field pressure at n+1 time; S3: calculating the force and torque of the complex particle; and S4: calculating the attitude of the complex particle through a six-degree-of-freedom model. The present application has the following effects: high calculation efficiency, reduced calculation complexity and cost, and improved simulation efficiency. Strong comprehensiveness, which can accurately describe the translation and rotation motion of the complex particle and provide comprehensive dynamic analysis. Good stability, which realizes stable solution of the coupling of fluid and particle motion and avoids numerical oscillation and instability. Multi-scale simulation capability, which can accurately simulate at different scales and meet the demand of multi-scale numerical simulation.

Description

Technical Field

[0001] This invention belongs to the field of flow field technology, specifically relating to a general simulation method for complex particle two-phase flow. Background Technology

[0002] The interaction between particles and fluids is a common physical phenomenon in nature and industrial production. To promote scientific and technological development and innovation, researchers have been dedicated to studying and exploring the problem of particle two-phase flow. This includes theoretical research, numerical simulation, and experimental research. Over the past few decades, significant progress has been made in the study of particle two-phase flow. Researchers have proposed various mathematical models and computational methods, such as the Euler-Lagrange method and multi-scale simulations, to describe the interaction between particles and fluids and the motion of particles. With the deepening of experimental research, complex particles have attracted considerable attention due to their unique morphology and function. However, in numerical simulation research, the main difficulty in studying complex particle two-phase flow lies in establishing mathematical models of complex particles and studying the complex interaction relationships between complex particles and fluids.

[0003] Particle mathematical models are crucial for describing the morphology and motion of particles. Previous studies have treated particles as point masses or spherical particles, considering the influence of particle translation on the flow field. Spherical particles have been widely used and studied. However, in practical applications, most particles are non-spherical, such as ellipsoidal particles, cylindrical particles, and complex particles. However, due to the high surface area ratio, complex surface morphology, and complex cross-sectional structure of these models, they are often limited to specific situations and cannot accurately capture the motion of complex particles.

[0004] Currently, numerical simulations of two-phase flows involving non-spherical particles are primarily based on the Euler-Lagrange solution framework. Within this framework, the Navier-Stokes equations for the fluid phase are solved on an Euler grid, while the particle phase is described using the Lagrange tracking method to depict the behavior of each individual particle. Solving for particles from the Lagrange perspective accurately describes their translational, orientation, and rotational behaviors, thus helping researchers gain a more intuitive understanding of the relationship between the complex motion of non-spherical particles and the flow field. Numerical simulation methods for particle two-phase flows include body-fitted grid methods based on the Arbitrary Lagrangian Eulerian (ALE) method, and non-body-fitted grid methods including the Immersed Boundary Method (IBM).

[0005] Arbitrary Lagrange-Eulerian (ALE) Method:

[0006] How it works: The ALE method handles boundary problems by moving the mesh along with the solid particles. At each time step, the mesh is regenerated and adapted to the new particle positions, ensuring consistency between the solid boundary and the fluid mesh.

[0007] Advantages: It can accurately capture moving boundaries and is suitable for simulating the motion of solid particles in fluids.

[0008] Disadvantages: Mesh regeneration and movement introduce high computational complexity and instability, especially when simulating complex particle shapes, where mesh quality is difficult to guarantee. Furthermore, frequent mesh updates increase computational costs and reduce simulation efficiency.

[0009] Deficiencies of existing technology:

[0010] 1. Inaccurate geometric description: The ALE method has low accuracy in geometric description when dealing with complex-shaped particles. The mesh regeneration process of the ALE method is prone to degrading mesh quality, especially for irregularly shaped particles.

[0011] 2. High computational complexity: The frequent mesh regeneration and movement in the ALE method leads to a significant increase in computational complexity and cost. This limits the simulation of large-scale complex systems in practical applications.

[0012] 3. Insufficient description of particle motion: Existing methods struggle to accurately describe the motion of complex-shaped particles in fluids, especially rotational motion. The ALE method has limitations in handling six-degree-of-freedom motion, affecting the accuracy of simulation results.

[0013] 4. Poor stability: Existing methods exhibit poor stability when dealing with large deformations and high-speed moving particles, easily leading to numerical errors and computational instability. The mesh update process of the ALE method is prone to numerical oscillations. Summary of the Invention

[0014] To address the aforementioned problems, this invention proposes a novel numerical simulation method aimed at solving the following technical issues:

[0015] 1. Improve the accuracy of geometric description: Introduce a geometric representation model (GRM) to solve the problem of inaccurate geometric description by accurately representing the boundaries of complex particles of various shapes.

[0016] 2. Reduced computational complexity and cost: By introducing the immersed boundary method (IBM), frequent mesh reconstruction and movement are avoided, significantly reducing computational complexity and improving simulation efficiency. This enables the invention to be applied to the numerical simulation of large-scale complex systems.

[0017] 3. Improve the description of particle motion: By using the six-degree-of-freedom (6DOF) model, the translational and rotational motions of complex particles can be accurately described, providing a more comprehensive dynamic analysis and solving the problem of insufficient description of particle motion in existing methods.

[0018] 4. Improve computational stability: Through the application of the characteristic based splitting (CBS) finite element method, stable coupled solutions of fluid and particle motions are achieved, avoiding numerical oscillations and instability problems, and improving the stability and reliability of calculations.

[0019] 5. Achieve multi-scale simulations: By combining the GRM, IBM, 6DOF, and CBS methods, accurate simulations can be carried out at different scales, meeting the requirements of multi-scale numerical simulations, and providing strong technical support for the comprehensive analysis of complex particle-fluid two-phase flow systems.

[0020] Technical solution:

[0021] The present invention provides a general simulation method for complex particle two-phase flow, including: [[ID=,]]

[0022] S1: Divide the flow field region of the particle two-phase flow into grids and construct a geometric characterization model to characterize the geometric information and motion information of the particles;

[0023] S2: Based on the grids divided in S1, use the characteristic line splitting finite element method to solve the flow field to obtain the flow field velocity at the pseudo n+1 moment and the flow field pressure at the n+1 moment;

[0024] S3: Determine the interaction region between the particles and the flow field according to the geometric information of the particles. Combine the flow field velocity and flow field pressure calculated in step S2, and calculate the forces on the particle boundary, the forces on the particle surface, and the torque acting on the particles through the IBM method; and construct a flow field velocity update formula to update the flow field velocity at the next moment according to the flow field velocity update formula;

[0025] S4: Use the forces on the particle surface and the torque acting on the particles calculated in S3 to calculate the particle attitude through the six-degree-of-freedom model; Given time T and time step dt, calculate t = k*dt, where k is the number of loops. When t < T, assign the updated flow field velocity to the flow field velocity at the n moment, and return to S2 to continue the loop; when t ≥ T, exit the loop and the calculation ends.

[0026] Preferably, S1 includes:

[0027] S1.1: Divide the flow field region into grids, denoted as D

[0024] , ,

[0025] , ,

[0022] , mesh , , ,

[0027] ,

[0023] , ,

[0028] , , ,

[0026] ;

[0028] S1.2: Use Xr Representing point X bk At the particle boundary Γ b Sets on:

[0029] X r =A k X r +(IA k )X bk k = 1, 2, ..., N;

[0030] Among them, A k Here is the rotation matrix used to fit the boundaries of all particles. Particle boundary From N boundary segments Γ i Composed of (i = 1, 2, ..., N);

[0031] S1.3: Based on S1.2, solve for the centroid of particles of arbitrary shape in the two-dimensional case:

[0032]

[0033] Among them, X r V represents the coordinates of the particle boundary; V is the area of ​​the particle in two dimensions and the volume of the particle in three dimensions.

[0034] S1.4: Determine the particle's centroid X at time n c Grid D mesh A cell is defined as the unit containing the remaining cells within a square or cube with a side length of 2dh centered on this cell, and the grid node coordinates (x) of these cells are stored. j .

[0035] Preferably, the grid D is divided in S1. mesh Define node velocity and pressure Then the flow field velocity of S2 and flow field pressure p n+1 The following equations are solved sequentially to obtain the following:

[0036]

[0037] in:

[0038]

[0039] K u =∫ Ω (LN u ) T (LN u )dΩ

[0040]

[0041] N u and N p These are the finite element basis functions defined in velocity space and pressure space, respectively; dt is the time step; M u The mass matrix represents the velocity. For the predicted flow field velocity at the node; The matrix represents the convection terms; η is the dynamic viscosity; K u B is the stiffness matrix related to the flow field velocity; u n K represents the second derivative of the flow field velocity. p Here, represents the stiffness matrix of the pressure; G represents the term related to the flow field pressure; H represents the coupling term between the flow field velocity and the flow field pressure. This is the second derivative of the pressure in the flow field; The corrected flow field velocity; and p n+1 These are the corrected flow field velocity and flow field pressure for the next time step, respectively.

[0042] Preferably, the surface force of the S3 particles is:

[0043]

[0044] Where, x j X represents the coordinates of the grid nodes. r Let F be the coordinates of the particle boundary point. lr (X r ,t) represents the particle boundary point X r The force at point U n+1 (X r F(t) represents the particle surface velocity at time n+1. p Let d be the total net external force acting on the particle, and let dh,dg and drg be the Euler mesh step size, Lagrange point spacing and particle boundary thickness, respectively.

[0045] Based on the center of gravity of particle S1, the torque T acting on particle p for:

[0046]

[0047] Preferably, S4 includes:

[0048] S4.1: The six-degree-of-freedom model of the particle is constructed as follows:

[0049]

[0050] Where F is F in S3 p M is the particle mass, U c,ω,Q,X c These are particle velocity, particle angular velocity, particle Euler angle, and particle displacement, respectively; I p This represents the moment of inertia of the particle. and Represents the coordinate transformation matrix between coordinate systems;

[0051] S4.2: Solving the discrete form of the 6DOF model using the fourth-order Runge-Kutta method:

[0052]

[0053] in:

[0054]

[0055] Preferably, the flow field velocity update formula in S4 is:

[0056]

[0057] Where x is the grid node coordinate, X r Let F be the coordinates of the particle boundary point. lr (X r ,t) represents the particle boundary point X r The force at the point, u n+1 (x,t) represents the flow velocity at time n+1. This represents the flow velocity at the pseudo-n+1 time. This represents the force transmitted from the particles into the flow field.

[0058] Compared with the prior art, the beneficial effects of the present invention are:

[0059] 1. High-precision geometric description: GRM can accurately describe the boundaries of complex particles of various shapes, improving the accuracy of geometric description and making it suitable for a variety of practical application scenarios.

[0060] 2. Highly efficient and stable computing: By using IBM, frequent grid reconstruction and movement are avoided, significantly reducing computational complexity and cost, and improving simulation efficiency and stability.

[0061] 3. Comprehensive motion simulation: 6DOF enables accurate description of the translational and rotational motion of complex particles, providing a more comprehensive dynamic analysis.

[0062] 4. Good computational stability: The CBS method achieves stable coupled solution of fluid and particle motion, avoiding numerical oscillation and instability problems, and improving the stability and reliability of the calculation.

[0063] 5. Multi-scale simulation: The method of this invention can perform accurate simulations at different scales, meeting the needs of multi-scale numerical simulation and providing strong technical support for the comprehensive analysis of complex particulate fluid two-phase flow systems. Attached Figure Description

[0064] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used together with the embodiments of the invention to explain the invention and do not constitute a limitation thereof.

[0065] In the attached diagram:

[0066] Figure 1 This is a flowchart of the method of the present invention;

[0067] Figure 2 This is the curved boundary of the present invention;

[0068] Figure 3 These are schematic diagrams of the particle shapes of the present invention: (a) schematic diagram of centrally symmetric particles; (b) schematic diagram of irregular particle shapes; (c) irregular particles generated by GRM; (d) symmetric particles generated by GRM;

[0069] Figure 4 (a) is a schematic diagram of the formation of a pentagram. Figure 4 (b) is a schematic diagram of the generation of the curved quadrilateral;

[0070] Figure 5 This is a schematic diagram of the submerged boundary method;

[0071] Figure 6(a) shows velocity cloud diagrams of sheet-like particles of different shapes: (a1): pentagon, (a2): hexagon, (a3): nonagon, (a4): fourteen-sided, (a5): icosahedron, (a6): twenty-nine-sided, (a7): triangle, (a8): pentagram, (a9): curved quadrilateral. Figure 6(b) is a schematic diagram of streamlines.

[0072] Figure 6(b) shows the streamlines of sheet-like particles of different shapes: (b1): pentagon, (b2): hexagon, (b3): ​​nonagon, (b4): fourteen-sided, (b5): icosagon, (b6): twenty-nine-sided, (b7): triangle, (b8): pentagram, (b9): curvilinear quadrilateral.

[0073] Figure 7 This is a schematic diagram of the coordinate system for a six-degree-of-freedom model.

[0074] Figure 8(a) shows the sedimentation process of square particles with different initial angles: (a1) angle is 0, (a2) angle is π / 4, (a3) ​​angle is π / 2, and (a4) angle is 3π / 4.

[0075] Figure 8(b) shows the evolution of the sedimentation process of the circular particles at times t = 0.2, 0.5, 0.7, and 0.8, respectively.

[0076] Figure 9(1) shows the vorticity cloud diagrams of two circular particles moving towards each other at different times: (a) t = 24, (b) t = 32;

[0077] Figure 9(2) shows the vorticity cloud diagrams of square and round particles moving towards each other at different times: (a) t = 24, (b) t = 32;

[0078] Figure 9(3) shows the vorticity cloud diagrams of two square particles moving towards each other at different times: (a) t = 24, (b) t = 32;

[0079] Figure 9(4) shows the evolution of the lift coefficient of two particles moving towards each other over time: where CCd represents the force on the lower circular particle in the case of two circular particles, CCu represents the force on the upper circular particle in the case of two circular particles, RCd represents the force on the lower circular particle in the case of square and circular particles, and RCu represents the force on the upper square particle in the case of square and circular particles.

[0080] Figure 9(5) shows the evolution of the drag coefficient of two particles moving towards each other over time. Detailed Implementation

[0081] The following is in conjunction with the appendix Figure 1 - Figure 9 illustrates a preferred embodiment of the present invention. It should be understood that the preferred embodiments described herein are for illustrative and explanatory purposes only and are not intended to limit the scope of the invention.

[0082] A general simulation method for complex particulate two-phase flow includes:

[0083] S1: Divide the flow field region of the particle two-phase flow into a grid and construct a geometric representation model to represent the geometric and motion information of the particles.

[0084] S1.1: Mesh the flow field region, denoted as D. mesh ;

[0085] S1.2: Assuming particle boundaries Composed of N boundary segments Γ i The composition of (i = 1, 2, ..., N). Boundary Γ i The length is l i , boundary Γ i and Γ i+1 The angle between them is represented by θ. i We choose a length of Boundary Γ b To approximate the boundary

[0086] Boundary Γ b Divided into N boundary segments Γ bi And assume the boundary segment Γ b1 The starting coordinates are x b1 Then the boundary segment Γ bi The starting coordinates are In this way, we can set m k coordinate point x ij ,(j=1,2,…,m k ) in each boundary segment Γ bi The Lagrange point inside serves as the particle boundary. Therefore, we have At the boundary Γ b The coordinates on the [space].

[0087] S1.2: Next, we use rotation matrix A k Boundary of fitted particles Specifically, we can use rotation matrix A k Transform each point x bk From Γ b Local coordinate system to particle boundary The global coordinate system. Rotation matrix A k It can be calculated based on the required rotation angle and can be applied to all points x. bk .

[0088]

[0089] X r =A k X r +(IA k )x bk k = 1, 2, ..., N. (1.2)

[0090] Where X is used r Representing point X bk At the particle boundary Γ b A is a set on the set. Where A is... k Here is the rotation matrix used to fit the boundaries of all particles. Particle boundary From N boundary segments Γ i Composed of (i = 1, 2, ..., N);

[0091] (like Figure 2 When dealing with curve boundaries, it is necessary to know the curvature of the curve boundary. The curvature of a curve is an indicator describing the degree of bending of the curve. The curvature of a curve can be calculated using curvature formulas or numerical approximation methods. If the curve boundary can be represented by the equation y = f(x), the formula for calculating the curvature is:

[0092]

[0093] However, in normal circumstances, the curve boundary is like... Figure 2 Curve curvature cannot always be represented by equations; therefore, a numerical approximation method is needed to solve for it. The curvature of the curve boundary can be obtained by solving for the radius R of the circumcircle of the triangle formed by three points E, A, and B on the curve boundary l. That is:

[0094]

[0095] Where S represents the area of ​​triangle ΔEAB. This allows us to obtain the fitted boundary of the curve using the parametric coordinates. For example... Figure 3 This is a schematic diagram of the particle shape of the present invention. The geometric characterization process can refer to the S1 geometric characterization model. Figure 4 As shown, GRM can effectively characterize complex particle morphologies under different shapes. Therefore, this method can be extended to high-dimensional cases and is applicable to various particle boundaries.

[0096] In previous studies of particulate two-phase flow, the particles under study had centrosymmetry, such as... Figure 3 (a) The geometric center of the rotating particle is usually used as the axis of rotation. This choice of axis of rotation simplifies the description and calculation of the problem and is consistent with reality. It should be noted that in practical applications, if the shape or mass distribution of the particles is not uniform, such as... Figure 3 (b) or, if subjected to an external torque, other factors may need to be considered to determine the position and orientation of the rotation axis.

[0097] S1.3: This section presents a method for calculating the centroid of particles with complex shapes. The formula for calculating the centroid of particles of arbitrary shapes in two dimensions is as follows:

[0098]

[0099] Among them, X r V represents the coordinates of the particle boundary; V is the area of ​​the particle in two dimensions and the volume of the particle in three dimensions.

[0100] S1.5: The method for determining the centroid of a 2D polygon involves decomposing the polygon into several triangles, then calculating the centroid of each sub-triangle using its area as a weight, and finally taking a weighted average of these centroids to obtain the polygon's centroid. This method, based on segmenting the polygon and utilizing the characteristics of each segmented triangle to calculate the centroid position of the entire polygon, is an effective and commonly used mathematical modeling technique. Through this decomposition and weighted averaging process, the centroid position of the polygon can be determined more accurately.

[0101] Assuming that polygonal particles are known, such as Figure 3 The coordinates of each vertex on the right (x) A ,x B ,x C ,x D ,x E ,x F ,x G Divide it into There are 5 triangles in total. The formulas for the area and centroid of a triangle are given below:

[0102]

[0103] Where, x i ,x j ,x k Let be the coordinates of the three vertices of the triangle.

[0104] S1.5: The centroid coordinates of the polygonal particle are:

[0105]

[0106] in |V i | These represent the centroid and area of ​​the triangle, respectively.

[0107] S1.6: Determine the centroid X of the particle at time n. c Located in mesh partition D mesh Within which cell, and within a square (cube) centered on this cell with a side length of 2dh, identify the remaining cells and store their node information x. j This provides flow field information for IBM interpolation calculations, such as... Figure 5 As shown.

[0108] S2: Based on the mesh generated in S1, the flow field is solved using the finite element method based on the characteristic line splitting to obtain the flow field velocity and flow field pressure at pseudo time n+1.

[0109] The flow field velocity at pseudo-n+1 time is obtained by solving the characteristic line splitting finite element method (CBS). The CBS method achieves stable coupling of fluid and particle motion, avoiding numerical oscillations and instabilities. The CBS method improves the stability and reliability of the calculation. Specific methods include:

[0110] S2.1: The Navier-Stokes equations are widely used to describe incompressible viscous flow problems and can be combined with the direct force method of the submerged boundary method to derive:

[0111]

[0112] Where, ρ f u, p, and η represent the fluid's density, velocity, pressure, and viscosity, respectively (μ = η / ρ). f Indicates the kinematic viscosity of the fluid. b U and F represent the force density near the immersion boundary. l These are the Lagrange point velocity and force of the particle, respectively. The function δ(·) is an interpolation function used to transfer physical information between the flow field and the particle, u D and u N These represent the Dirichlet boundary conditions and the Neuman boundary conditions, respectively.

[0113] The CBS method effectively overcomes the problems of numerical oscillation and instability, providing a stable and accurate solution for velocity and pressure.

[0114] We choose a linear finite element shape function to approximate the velocity u and pressure p:

[0115]

[0116] in and This represents the quality of the node.

[0117] S2.2: By discretizing the momentum equation (2.1) along the characteristic line in time, a semi-discrete form is obtained:

[0118]

[0119] The underlined part on the right side of equation (2.2) represents the stabilizing term, which effectively eliminates numerical oscillations and instabilities caused by the dominant convection.

[0120] S2.3: Solving for the prediction step, we obtain an approximate intermediate velocity u by simplifying equation (2.2). * .

[0121]

[0122] By defining an auxiliary matrix:

[0123]

[0124] The matrix form of equation (2.3) can be obtained:

[0125]

[0126] in:

[0127]

[0128] K u =∫Ω (LN u ) T (LN u )dΩ.

[0129] S2.4: Next, use the intermediate speed u * To calculate the pressure, we need to solve the Poisson equation for pressure derived from the divergence of the continuity and momentum equations. The pressure field is then updated based on intermediate velocities to ensure that the incompressibility condition is satisfied.

[0130]

[0131] because but:

[0132]

[0133] S2.5: The matrix form is as follows:

[0134]

[0135] in,

[0136] S2.6: Equation (2.8) is used to calculate the pseudo-velocity at the next time step.

[0137]

[0138] S2.7: Define nodal velocities on the mesh defined in S1. and pressure S2 flow field velocity and pressure p n +1 The following equations are obtained by solving them sequentially:

[0139]

[0140] in:

[0141]

[0142] K u =∫ Ω (LN u ) T (LN u )dΩ

[0143]

[0144] N u and N p These are the finite element basis functions defined in velocity space and pressure space, respectively, where dt is the time step and M is the pressure space.u The mass matrix is ​​the velocity matrix, derived from the fluid density ρ. f Finite element basis functions N for flow field velocity u Obtained through integration; The predicted flow field velocity at the node is the node flow field velocity at the current time step. and the nodal flow velocity at the next time step An intermediate estimate between; The matrix represents the convection term, and is related to the nodal velocities of the flow field. Correlation reflects the convection effect of the fluid; η is the dynamic viscosity, a measure of the internal friction of the fluid; K u B is the stiffness matrix related to the flow field velocity; u n The second derivative of the flow field velocity is given by the nodal velocity of the flow field. Related; K p Let N be the stiffness matrix of the pressure, and N be the finite element basis function of the flow field pressure. p G is obtained through integration; G is a term related to flow field pressure, obtained through integration, reflecting the influence of flow field velocity on flow field pressure; H is the coupling term between flow field velocity and flow field pressure. The second derivative of the flow field pressure is used to predict the flow field velocity. Related; The corrected flow field velocity is used to calculate the flow field velocity in the next time step; and p n+1 These are the corrected flow field velocity and flow field pressure for the next time step, respectively.

[0145] S3: Determine the interaction region between the particle and the flow field based on the particle's geometric information. Combine the flow field velocity and pressure calculated in step S2, calculate the particle boundary force, particle surface force, and torque acting on the particle using the IBM method. Construct a flow field velocity update formula and update the flow field velocity at the next moment based on the flow field velocity update formula.

[0146] Calculating forces and moments in complex particles (IBM): Utilizing IBM technology avoids frequent mesh reconstruction and movement, thereby reducing computational complexity and cost. By simulating fluid motion on Euler meshes and combining IBM with particle boundary treatment, efficient and stable numerical simulations are achieved. Specifically, this includes:

[0147] IBM's implementation typically involves two steps:

[0148] 1. Interpolation: The flow field information is interpolated from the Euler grid to the Lagrange points using a regularization function. Then, the IB force at the Lagrange points is calculated to update the state of the particle boundaries.

[0149] 2. Propagation: The flow field information is updated using a regularization function, and the IB force obtained at the Lagrange point is propagated back to the Euler mesh.

[0150] These two steps facilitate information exchange between the fluid and solid domains, making fluid-solid coupling easier to achieve. Interpolation and propagation in Figure 5 Performed at the Lagrange point and the adjacent Euler grid.

[0151] S3.1: Assume the number of local Euler grid points around the particle boundary is N. E The number of Lagrange points is N. L .

[0152] Calculate the velocity at the Lagrange point and the force density at the Euler grid point using the following equations.

[0153]

[0154] Here, f b F represents the force density at a point on the Euler grid, while F li This represents the force acting at the i-th Lagrange point. The Lagrange point is denoted by X, and its velocity is denoted by U. The parameters dh, dg, and drg represent the Euler grid step size, the Lagrange point spacing, and the particle boundary thickness, respectively.

[0155] S3.2: The equations corresponding to the implicit velocity correction method are shown below:

[0156]

[0157] S3.3: Substituting equation (3.1) into equation (3.2), we get:

[0158]

[0159] S3.4: The interpolation and propagation operations follow the unit partition condition, where the force at the Lagrange point is first propagated to the Euler grid point and then interpolated back to obtain the value at the Lagrange point. This condition is expressed by equation (3.4):

[0160]

[0161] S3.5: Combining equations (3.1) and (3.4), and then substituting them into equation (3.3), we obtain the expression for the Lagrange point force as follows:

[0162]

[0163] S3.6: The regularization function plays a crucial role in establishing communication between the particles and the fluid. It should be a continuous and differentiable function with second-order accuracy. We provide the following options for the regularization function in IBM to meet these requirements:

[0164]

[0165] where,

[0166] S3.7: The surface force on the particle can be calculated using Equation (3.7):

[0167]

[0168] where, x j is the grid node coordinate, X r is the particle boundary point coordinate, F lr (X r ,t) is the force on the particle boundary point X r , U n+1 (X r ,t) represents the particle surface velocity at time n+1, F p is the total external force on the particle, and the parameters dh, dg, and drg are the Euler grid step size, Lagrange point spacing, and particle boundary thickness, respectively;

[0169] S3.9: Combining the center of gravity of the particle in S1, the torque T p acting on the particle can be obtained from Equation (3.8):

[0170]

[0171] S4: Using the surface force on the particle and the torque acting on the particle calculated in S3, calculate the particle attitude through the six-degree-of-freedom model; given time T and time step dt, calculate t = k*dt, where k is the loop count. When t < T, assign the updated flow field velocity to the flow field velocity at time n, and return to S2 to continue the loop; when t ≥ T, exit the loop and the calculation ends.

[0172] Calculate the complex particle attitude through the six-degree-of-freedom model; use the 6DOF model to accurately describe the motion of complex particles in the fluid, including translational and rotational motions. Through the 6DOF model, the dynamic behavior of the particles can be comprehensively captured, improving the accuracy of the simulation.

[0173] To accurately describe the absolute position of the complex particle in space, the relative position with respect to its center of mass, and the orientation (or attitude) of the complex particle. The following Figure 7 reference frames will be used to explain the kinematic behavior of the complex particle, which are respectively:

[0174] 1) Fixed spatial coordinate system S g =(Ox g y g z g ): The origin of the coordinate system is located at the particle's center of mass, and the z-axis is... g Vertically downwards, x-axis g and y g Located in the horizontal plane, axis x g Parallel to the particle axis pointing forward, axis y g Perpendicular to the particle plane, pointing to the right;

[0175] 2) Particle body coordinate system S b =(Ox b y b z b ): The origin of the coordinate system is located at the particle's center of mass, and the x-axis is... b Parallel to the particle axis and pointing in the direction of fluid flow, axis y b Perpendicular to the particle plane pointing to the right, axis z b Perpendicular to plane Ox b y b ;

[0176] 3) Inertial coordinate system S = (Oxyz).

[0177] Particle body coordinate system S b With spatial fixed coordinate system S g Euler angles exist

[0178] 1) Yaw angle ψ: x b When projected onto the horizontal plane, it is perpendicular to the x-axis. g The angle between x and x. b The projection of the positive semi-axis is located at x g When on the right, ψ > 0, and the range of ψ is ψ∈[-π,π].

[0179] 2) Pitch angle θ: x b The angle with the horizontal plane. When x b The positive semi-axis lies in the plane Ox g y g When above, θ > 0, and the range of θ is θ∈[-π / 2,π / 2];

[0180] 3) Roll angle z b With through axis x g The angle between z and the vertical plane. b When the positive semi-axis is located on the left side of the vertical plane, The range of values ​​is

[0181] The motion of rigid, complex particles in fluids, studied in this paper, is highly complex. Due to the shape and size of the particles, they undergo both translation and rotation during their motion. Translation causes changes in the particle's position, while rotation results in changes in its angle. In the three-dimensional case, describing these two motions requires velocities in three directions and three Euler angles; therefore, we employ a six-degree-of-freedom equation to describe the particle's motion in the fluid. By solving this six-degree-of-freedom equation, we can obtain the particle's trajectory and attitude changes in the fluid. All motions studied in this paper obey Newton's second law.

[0182] The translational motion of particles can be described by Newton's second law:

[0183]

[0184] Where m,g,F=(F x ,F y ,F z ) represent the mass, acceleration, and net force acting on the particle, respectively.

[0185] Simplifying equation (4.1) yields:

[0186]

[0187] This yields the following translation equation:

[0188]

[0189] in:

[0190]

[0191] The rotation of the particles can be described by Newton's second law of motion for rigid bodies:

[0192]

[0193] Among them, I xx ,I xy ,I xz ,I yx ,I yy ,I yz ,I zx ,I zy ,I zz For the moment of inertia, T = (T x ,T y ,T z Let be the torque acting on the particle; simplifying equation (4.4) yields:

[0194]

[0195] This yields the following rotational equation:

[0196]

[0197] in

[0198]

[0199] I xx =∫(y 2 +z 2 )dm,I yy =∫(x 2 +z 2 )dm,I zz =∫(x 2 +y 2 )dm,

[0200] I xy =∫xydm,I yz =∫yzdm,I zx =∫zxdm.

[0201] Equation (4.6) uses angular acceleration as the variable, but the six-degree-of-freedom equations solve for Euler angles. Therefore, the relationship between angular acceleration and Euler angles is given:

[0202]

[0203] Simplifying equation (4.7) yields:

[0204]

[0205] Therefore, we can conclude that:

[0206]

[0207] in

[0208]

[0209] S4.1: Therefore, the six-degree-of-freedom model (6DOF model) of the particle is constructed as follows:

[0210]

[0211] Where F is F in S3 p M is the particle mass, U c ,ω,Q,X c These are particle velocity, particle angular velocity, particle Euler angle, and particle displacement, respectively; I p This represents the moment of inertia of the particle. and Represents the coordinate transformation matrix between coordinate systems;

[0212] S4.8: Solve the discrete form of the 6DOF model using the fourth-order Runge-Kutta method:

[0213]

[0214] in:

[0215]

[0216] The speed update formula described in S4 is:

[0217]

[0218] Where x is the grid node coordinate, X r Let F be the coordinates of the particle boundary point. lr (X r ,t) represents the force at the particle boundary, u n+1 (x,t) represents the flow velocity at time n+1. This represents the flow velocity at the pseudo-n+1 time. This represents the force transmitted from the particles into the flow field.

[0219] Figures 6, 8, and 9 show the results obtained using the present invention. It can be seen that no matter how the shape of the particles changes, GRM can accurately capture their geometric features and achieve good simulation results for moving particles.

[0220] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of this invention is defined by the appended claims and their equivalents.

Claims

1. A general simulation method for complex granular two-phase flow, characterized in that: Comprise: S1: the flow field region of the particle two-phase flow is divided into a grid, and a geometric characterization model is constructed to characterize the geometric information and motion information of the particles; S2: on the basis of the grid divided in S1, the flow field velocity at pseudo n+1 time and the flow field pressure at n+1 time are obtained by solving the flow field by using the characteristic line splitting finite element method; S3: the particle interaction region with the flow field is determined according to the geometric information of the particles, and the particle boundary force, the particle surface force and the torque acting on the particles are calculated by the IBM method in combination with the flow field velocity and the flow field pressure calculated in the step S2; and a flow field velocity updating formula is constructed, and the flow field velocity at the next time is updated according to the flow field velocity updating formula; S4: using the particle surface stress and torque acting on the particle calculated by S3, the particle attitude is calculated by the six-degree-of-freedom model; given time T, time step dt, calculate , k is the number of cycles, when t < T, the updated flow field velocity is assigned to the flow field velocity at time n, and the cycle continues to S2; when t ≥ T, the cycle is exited and the calculation is completed.

2. The method of claim 1, wherein: S1 includes: S1.1: meshing the flow field region, denoted as ; S1.2: with representative points at the grain boundaries on the set: ； wherein, is a rotation matrix used to fit the boundary of all particles , the boundary of a particle is composed of segments boundaries where i takes values 1, 2,..., N; S1.3: on the basis of S1.2, the center of gravity of the particle of any shape in the two-dimensional case is solved: ； wherein is the coordinate of the particle boundary; is the area of the particle in the two-dimensional case, or the volume of the particle in the three-dimensional case; S1.4: determining the center of gravity of the particle at time n grid in which the cell is located unit, and the remaining cells contained in a square or cubic body with the unit as the center and with the side length of the unit grid node coordinates of the remaining cells .

3. A general simulation method of a complex granular two-phase flow according to claim 2, characterized in that: Grid D at S1 division mesh Node velocities are defined above And pressures Then the flow field velocities at S2 are And flow field pressures Solved in turn by the following equations: ； ； ； ； Wherein: ； ； ； and are the finite element basis functions defined in the velocity space and the pressure space, respectively; is the time step; is the mass matrix of the velocity; is the predicted flow field velocity at the nodes; is the matrix of the convection term; is the dynamic viscosity; is the fluid density; is the stiffness matrix related to the flow field velocity; is the second derivative term of the flow field velocity; is the stiffness matrix of the pressure; is the term related to the flow field pressure; is the coupling term between the flow field velocity and the flow field pressure; is the second derivative term of the flow field pressure; is the corrected flow field velocity; and are the corrected flow field velocity and the flow field pressure of the next time step, respectively.

4. The method of claim 3, wherein: The surface force of the particle in S3 is: ； ； where, is the grid node coordinate, is the particle boundary point coordinate, is the particle boundary point force at, represents the particle surface velocity at time n+1, is the total external force on the particle, parameter , and are the Euler grid step, the Lagrange point spacing and the particle boundary thickness, respectively. The center of gravity of the particle in combination with S1 acts a torque on the particle Is: 。 5. The method of claim 4, wherein: S4 includes: S4.1: a six-degree-of-freedom model of the particle is constructed as: ； where F is the force in S3 , M is the particle mass, are the particle velocity, angular velocity, Euler angles and displacement, respectively; denotes the moment of inertia of the particle, and denotes the coordinate transformation matrix between coordinate systems; S4.2: the discrete form of the 6DOF model is solved by using the fourth-order Runge-Kutta method: ； Wherein: 。 6. The method of claim 5, wherein: The flow field velocity updating formula in S4 is: ； ； wherein is the grid node coordinate, is the particle boundary point coordinate, is the force on the particle boundary point at time n+1, is the flow field velocity at time n+1, is the flow field velocity at time n+1, is the flow field velocity at pseudo-time n+1.

Images