General simulation method for complex particle two-phase flow

Abstract

The invention belongs to the technical field of flow fields, and particularly relates to a universal simulation method for complex particle two-phase flow. Comprising the following steps: S1, constructing a geometric representation model to represent geometric information of complex particles and dividing grids; s2, on the basis of the grids divided in S1, solving the flow field by adopting a finite element method based on characteristic line splitting to obtain the flow field speed at the pseudo (n + 1) moment and the flow field pressure at the (n + 1) moment; s3, calculating the stress and moment of force of the complex particles; and S4, calculating a complex particle attitude through a six-degree-of-freedom model. The method has the advantages that the calculation efficiency is high, the calculation complexity and cost are reduced, and the simulation efficiency is improved. The translation and rotation motion of complex particles can be accurately described, and comprehensive kinetic analysis is provided. And stability is good, stable solution of fluid and particle motion coupling is achieved, and numerical oscillation and instability are avoided. And multi-scale simulation capability: accurate simulation can be carried out on different scales, and the requirements of multi-scale numerical simulation are met.

Description

Technical Field

[0001] The invention belongs to the technical field of flow fields, and in particular relates to a general simulation method for complex particle two-phase flow. Background Art

[0002] The interaction between particles and fluids is a common physical phenomenon in nature and industrial production. In order to promote the development and innovation of science and technology, people have been committed to studying and exploring the problem of particle two-phase flow. This includes theoretical research, numerical simulation and experimental research. In the past few decades, great research progress has been made in the problem of particle two-phase flow. Researchers have proposed various mathematical models and calculation methods, such as the Euler-Lagrange method and multi-scale simulation, to describe the interaction between particles and fluids and the movement laws of particles. With the deepening of experimental research, complex particles have attracted much attention due to their unique morphology and functions. In numerical simulation research, the main difficulty of complex particle two-phase flow lies in establishing mathematical models of complex particles and studying the complex interaction relationship between complex particles and fluids.

[0003] Mathematical models of particles are the key to describing the morphological structure and motion state of particles. Previous studies have regarded particles as point masses or spherical particles, and considered the impact of particle translation on the flow field. Spherical particles have been widely used and studied, but 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, these models are often limited to specific situations and cannot accurately capture the motion of complex particles.

[0004] At present, the numerical simulation of two-phase flow with non-spherical particles is mainly based on the Euler-Lagrange solution framework. In this framework, the Navier-Stokes equations of the fluid phase are solved on the Euler grid, while the particle phase uses the Lagrange tracking method to describe the behavior of each particle. The solution of particles from the Lagrange point of view can accurately describe the translation, orientation and rotation behavior of particles, which can help researchers to more intuitively understand the relationship between the complex motion behavior of non-spherical particles and the flow field. Numerical simulation methods for two-phase flow with particles include body-fitting grid methods based on the Arbitrary Lagrangian Eulerian (ALE) method, non-body-fitting grid methods including the Immersed Boundary Method (IBM) method.

[0005] Arbitrary Lagrangian Euler (ALE) method:

[0006] How it works: The ALE method handles boundary problems by moving the mesh 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 movement of solid particles in fluids.

[0008] Disadvantages: Mesh regeneration and movement bring high computational complexity and instability, especially when simulating complex-shaped particles, the mesh quality is difficult to guarantee. In addition, frequent mesh updates increase the computational cost and reduce simulation efficiency.

[0009] Defects of the prior art:

[0010] 1. Inaccurate geometric description: The ALE method has low geometric description accuracy when dealing with complex-shaped particles. The mesh regeneration process of the ALE method can easily lead to a decrease in mesh quality, especially for irregular-shaped particles.

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

[0012] 3. Inadequate description of particle motion: Existing methods have difficulty accurately describing the motion of complex-shaped particles in fluids, especially rotational motion. The ALE method has limitations when dealing with six-degree-of-freedom motion, which affects the accuracy of the simulation results.

[0013] 4. Poor stability: When dealing with large deformation and high-speed moving particles, the existing methods have poor stability and are prone to numerical errors and computational instability. The grid update process of the ALE method is prone to numerical oscillation. Summary of the invention

[0014] In view of the above problems, the present invention proposes a new numerical simulation method, aiming to solve the following technical problems:

[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. Reduce computational complexity and cost: By introducing the immersed boundary method (IBM), frequent grid reconstruction and movement are avoided, computational complexity is significantly reduced, and simulation efficiency is improved, which enables the present invention to be applied to numerical simulation of large-scale complex systems.

[0017] 3. Improve the description of particle motion: By adopting the six-degree-of-freedom (6DOF) model, it can accurately describe the translational and rotational motions of complex particles, provide a more comprehensive dynamic analysis, and solve 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 solutions:

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

[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: On the basis of 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. Combining the flow field velocity and flow field pressure calculated in step S2, 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 mesh ;

[0028] S1.2: Use Xr Indicates point X bk At the particle boundary Γ b The collection on:

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

[0030] Among them, A k is the rotation matrix used to fit the boundaries of all particles Boundaries of particles By N segment boundaries Γ i , (i=1,2,…,N);

[0031] S1.3: Based on S1.2, solve the center of gravity of particles of arbitrary shape in two dimensions:

[0032]

[0033] Among them, X r are 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 center of gravity X at time n c Grid D mesh Unit, with this unit as the center, the rest of the units contained in the square or cube with a side length of 2dh, and save the grid node coordinates x of these units j .

[0035] Preferably, the grid D divided in S1 mesh Define the node speed and pressure Then the flow field velocity of S2 is and flow field pressure p n+1 Solve the following equations in turn:

[0036]

[0037] in:

[0038]

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

[0040]

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

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

[0043]

[0044] Among them, x j is the grid node coordinate, X r is the coordinate of the particle boundary point, F lr (X r ,t) is the particle boundary point X r The force at 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;

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

[0046]

[0047] Preferably, S4 includes:

[0048] S4.1: Construct the six-degree-of-freedom model of the particle as:

[0049]

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

[0051] S4.2: Solve 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 is the coordinate of the particle boundary point, F lr (X r ,t) is the particle boundary point X r The force at n+1 (x, t) represents the flow field velocity at time n+1, represents the flow field velocity at pseudo time n+1, is the force transmitted by the particle to the flow field.

[0058] Compared with the prior art, the present invention has the following beneficial effects:

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

[0060] 2. Efficient and stable computing: IBM can avoid frequent grid reconstruction and movement, significantly reduce computing complexity and cost, and improve simulation efficiency and stability.

[0061] 3. Comprehensive motion simulation: Through 6DOF, the translational and rotational motion of complex particles can be accurately described, providing a more comprehensive dynamic analysis.

[0062] 4. Good computational stability: Through the CBS method, a stable coupling solution of fluid and particle motion is achieved, numerical oscillation and instability problems are avoided, and the stability and reliability of the calculation are improved.

[0063] 5. Multi-scale simulation: The method of the present invention can perform accurate simulation at different scales, meet the needs of multi-scale numerical simulation, and provide strong technical support for the comprehensive analysis of complex particle-fluid two-phase flow systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0064] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0065] In the attached picture:

[0066] Figure 1 is a flow chart of the method of the present invention;

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

[0068] Figure 3 Schematic diagrams of particle shapes of the present invention: (a) schematic diagram of central symmetrical particles; (b) schematic diagram of irregular particle shapes; (c) irregular particles generated by GRM; d) symmetrical particles generated by GRM;

[0069] Figure 4 (a) is a schematic diagram of the five-pointed star generation. Figure 4 (b) is a schematic diagram of curved quadrilateral generation;

[0070] Figure 5 Schematic diagram of the immersed boundary method;

[0071] Fig. 6(a) is a velocity cloud diagram of flake particles of different shapes: (a1): pentagon, (a2): hexagon, (a3): enneagon, (a4): tetradecagon, (a5): icosagon, (a6): icosagon, (a7): triangle, (a8): pentagon, (a9): quadrilateral, Fig. 6(b) is a schematic diagram of the streamline diagram;

[0072] Figure 6(b) is the streamline diagram of flake particles of different shapes: (b1): pentagon, (b2): hexagon, (b3): ​​enneagon, (b4): tetradecagon, (b5): icosagon, (b6): icosagon, (b7): triangle, (b8): pentagon, (b9): quadrilateral:

[0073] Figure 7 Schematic diagram of the coordinate system of the six-degree-of-freedom model;

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

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

[0076] Figure 9 (1) shows the vorticity cloud diagram 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 diagram 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: 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 resistance coefficient of two particles moving towards each other over time. DETAILED DESCRIPTION

[0081] The following is combined with Figure 1 - Figure 9 illustrates a preferred embodiment of the present invention. It should be understood that the preferred embodiment described herein is only used to illustrate and explain the present invention, and is not used to limit the present invention.

[0082] A general simulation method for complex particle two-phase flows, including:

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

[0084] S1.1: Divide the flow field area into grids, denoted as D mesh ;

[0085] S1.2: Assume that the particle boundaries N boundary segments Γ i ,(i=1,2,…,N). Boundary Γ i The length is l i , boundary Γ i and Γ i+1 The angle between them is denoted as θ i We choose a length of The boundary of Γ b To approximate the boundary

[0086] Boundary Γ b is divided into N boundary segments Γ bi , and assume that the boundary segment Γ b1 The starting coordinate is 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 As the Lagrange point of the particle boundary. Therefore, we have At the boundary Γ b The coordinates on .

[0087] S1.2: Next, we use the rotation matrix A k Fitting the Boundaries of Particles Specifically, we can use the rotation matrix A k Transform each point x bk From Γ b The local coordinate system of the particle boundary The global coordinate system of . The rotation matrix A k can be calculated based on the desired 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 r Indicates point X bk At the particle boundary Γ b A k is the rotation matrix used to fit the boundaries of all particles Boundaries of particles By N segment boundaries Γ i , (i=1,2,…,N);

[0091] (like Figure 2 ) When dealing with the boundary of a curve, you need to know the curvature of the boundary of the curve. The curvature of a curve is an indicator that describes the degree of curvature of the curve. The curvature of a curve can be calculated using the curvature formula or numerical approximation method. If the boundary of the curve can be represented by the equation y = f(x), the curvature calculation formula is:

[0092]

[0093] But usually, the curved boundary is Figure 2 It cannot always be expressed by an equation, so a numerical approximation method is needed to solve the curvature of the curve. By solving the radius R of the circumscribed circle of the triangle formed by the three points E, A, and B on the curve boundary l, the curvature of the curve boundary can be obtained, that is:

[0094]

[0095] Where S represents the area of ​​triangle ΔEAB. In this way, the fitting boundary of the curve boundary can be obtained through parameter coordinates. Figure 3 , is a schematic diagram of the particle shape of the present invention, and the geometric characterization process can refer to the S1 geometric characterization model. Figure 4 As shown in the figure, it can be seen that GRM can well characterize complex particle morphologies with different shapes. Therefore, this method can be extended to high-dimensional situations and is applicable to various particle boundaries.

[0096] In previous studies on particle two-phase flow, the particles studied have central symmetry, such as Figure 3 (a), the geometric center of the rotating particle is usually used as the rotation axis. This way of choosing the rotation axis simplifies the description and calculation of the problem and is consistent with the actual situation. It should be noted that in practical applications, if the shape or mass distribution of the particles is uneven, such as Figure 3 (b) or is subject to external torques, further considerations may be required to determine the position and orientation of the axis of rotation.

[0097] S1.3: Here is a method to solve the centroid of complex shaped particles. The formula for calculating the centroid of particles of arbitrary shape in two dimensions is:

[0098]

[0099] Among them, X r are 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 solving the centroid of a two-dimensional polygon can be to decompose the polygon into several triangles, then use the area of ​​each triangle as the weight, calculate the centroid of each sub-triangle, and finally perform weighted average of these centroids to obtain the centroid of the polygon. This method is based on segmenting the polygon and then using the characteristics of each segmented triangle to calculate the centroid position of the entire polygon. It is an effective and commonly used mathematical modeling technique. Through such a decomposition and weighted average process, the centroid position of the polygon can be determined more accurately.

[0101] Assume that the polygonal particles are known as Figure 3 The coordinates x of the right vertices A ,x B ,x C ,x D ,x E ,x F ,x G , which is divided into There are 5 triangles in total. The formulas for the area and centroid of the triangles are given below:

[0102]

[0103] Among them, x i ,x j ,x k are the coordinates of the three vertices of the triangle.

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

[0105]

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

[0107] S1.6: Determine the particle center of gravity X at time n c Located in the mesh D mesh Which unit is in which, and with this unit as the center, the rest of the units contained in the square (cube) with a side length of 2dh, and save the node information x of these units j , providing flow field information for IBM interpolation calculation, such as Figure 5 shown.

[0108] S2: Based on the grid divided by S1, the flow field is solved by the characteristic line splitting finite element method to obtain the flow field velocity at the pseudo n+1 moment and the flow field pressure at the n+1 moment;

[0109] The characteristic line splitting finite element method (CBS) is used to solve the flow field and obtain the flow field velocity at the pseudo n+1 moment. The CBS method is used to achieve a stable coupling solution of fluid and particle motion to avoid numerical oscillation and instability. The CBS method can improve the stability and reliability of the calculation. The 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 immersed boundary method to obtain:

[0111]

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

[0113] The CBS method effectively overcomes the problems of numerical oscillation and instability and provides stable and accurate solutions for velocity and pressure.

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

[0115]

[0116] in and Represents the node quality.

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

[0118]

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

[0120] S2.3: To solve the prediction step, we simplify equation (2.2) to produce an approximate intermediate speed u * .

[0121]

[0122] By defining the 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 compute the pressure. This involves solving Poisson's equation for pressure derived from the divergence of the continuity and momentum equations. The pressure field is updated based on the intermediate velocities to ensure that the incompressibility condition is met:

[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 for the next time step

[0137]

[0138] S2.7: Define nodal velocities on the grid divided by S1 and pressure Flow velocity of S2 and pressure p n +1 By solving the following equations in sequence:

[0139]

[0140] in:

[0141]

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

[0143]

[0144] N u and N p are the finite element basis functions defined in velocity space and pressure space respectively, dt is the time step, Mu is the mass matrix of velocity, given by the fluid density ρ f and the finite element basis function N of the flow field velocity u Obtained through points; is the predicted flow field velocity at the node, which is the node flow field velocity at the current time step and the node flow velocity at the next time step An intermediate estimate between is the matrix of the convection term, and is related to the flow field node velocity is related, reflecting the convection effect of the fluid; η is the dynamic viscosity, a measure of the internal friction of the fluid; K u is the stiffness matrix related to the flow field velocity; B u n is the second-order derivative of the flow field velocity, which is related to the flow field node velocity Related; K p is the stiffness matrix of pressure, and the finite element basis function N of flow field pressure p Obtained by integration; G is a term related to flow field pressure, obtained by integration, reflecting the influence of flow field velocity on flow field pressure; H is a coupling term between flow field velocity and flow field pressure; is the second-order derivative of the flow field pressure, which is related to the predicted flow field velocity Related; is the corrected flow field velocity, which is used to calculate the flow field velocity of the next time step; and p n+1 are the corrected flow field velocity and flow field pressure for the next time step respectively.

[0145] S3: Determine the interaction area between the particle and the flow field based on the geometric information of the particle. Combined with the flow field velocity and flow field pressure calculated in step S2, the particle boundary force, particle surface force and torque acting on the particle are calculated by the IBM method. Construct a flow field velocity update formula, and update the flow field velocity at the next moment according to the flow field velocity update formula.

[0146] Calculation of complex particle forces and moments (IBM): IBM is used to avoid frequent grid reconstruction and movement, thereby reducing computational complexity and cost. By simulating fluid motion on the Euler grid and combining IBM to process particle boundaries, efficient and stable numerical simulation is achieved. Specifically, it 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 point using a regularization function. Then, the IB force at the Lagrange point is calculated to update the state of the particle boundary.

[0149] 2. Propagation: Use the regularization function to update the flow field information and propagate the IB force obtained at the Lagrange point back to the Euler grid.

[0150] These two steps facilitate the exchange of information between the fluid and solid domains, which makes it easier to implement fluid-solid coupling. Figure 5 It is performed at the Lagrange points in and adjacent Euler grids.

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

[0152] The velocity at the Lagrange point and the force density at the Euler grid point are calculated using the following equations

[0153]

[0154] Here, f b represents the force density at the Euler grid point, and F li represents the force on the ith Lagrange point. The Lagrange point is represented by X and its velocity is represented by U. The parameters dh, dg and drg are the Euler mesh step size, the Lagrange point spacing and the particle boundary thickness.

[0155] S3.2: The equation corresponding to the implicit velocity correction method is as follows:

[0156]

[0157] S3.3: Substituting equation (3.1) into equation (3.2) yields:

[0158]

[0159] S3.4: The interpolation and propagation operations follow the unit partitioning condition, where the Lagrange point force is first propagated to the Euler grid points 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 equation (3.1) and equation (3.4) and 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 resultant 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 particle surface force 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 number of loops. 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 complex particles in space, the relative position relative to their center of mass, and the orientation (or attitude) of complex particles. The following Figure 7 reference frames will be used to explain the kinematic behavior of complex particles, which are respectively:

[0174] 1) Space fixed coordinate system S g =(Ox g y g z g ): The origin of the coordinate is located at the particle mass center, and the axis z g Vertically downward, axis x g and g Located in the horizontal plane, axis x g Parallel to the particle axis pointing forward, axis y g perpendicular to the plane of the particle pointing to the right;

[0175] 2) Particle body coordinate system S b =(Ox b y b z b ): The origin of the coordinate is located at the particle mass center, and the axis x b The y axis is parallel to the particle axis and points in the direction of fluid flow. 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 the space fixed coordinate system S g There are Euler angles

[0178] 1) Yaw angle ψ: x b When projected onto the horizontal plane, it is aligned with the axis x g When x b The projection of the positive semi-axis is located at x g On the right side, ψ>0, and the value range of ψ is ψ∈[-π,π];

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

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

[0181] The motion of the rigid complex particles studied in this paper in the fluid is very complex. Due to the shape and size of the particles, the particles will produce both translation and rotation during the motion. Translation causes the position of the particles to change, while rotation causes the angle of the particles to change. In the three-dimensional case, describing these two motions requires velocities in three directions and three Euler angles, so we use the six-degree-of-freedom equation to describe the motion of the particles in the fluid. By solving the six-degree-of-freedom equation, we can obtain the motion trajectory and posture changes of the particles in the fluid. The motions studied in this paper all obey Newton's second law.

[0182] The translation process of particles can be described according to Newton's second law:

[0183]

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

[0185] Simplifying equation (4.1) yields:

[0186]

[0187] This results in the following translation equation:

[0188]

[0189] in:

[0190]

[0191] The rotation process of the particle can be described according to Newton's second law of rotation of rigid bodies:

[0192]

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

[0194]

[0195] This results in the following rotation 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 a variable, but the Euler angle is solved in the six-degree-of-freedom equation. So the relationship between angular acceleration and Euler angle is given as:

[0202]

[0203] Simplifying equation (4.7) yields:

[0204]

[0205] From this we get:

[0206]

[0207] in

[0208]

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

[0210]

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

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

[0213]

[0214] in:

[0215]

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

[0217]

[0218] Where x is the grid node coordinate, X r is the coordinate of the particle boundary point, F lr (X r ,t) is the particle boundary force, u n+1 (x, t) represents the flow field velocity at time n+1, represents the flow field velocity at pseudo time n+1, is the force transmitted by the particle to the flow field.

[0219] FIG6 , FIG8 , and FIG9 are respectively the results obtained by 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 can also achieve good simulation results for moving particles.

[0220] The above shows and describes the basic principles, main features and advantages of the present invention. It should be understood by those skilled in the art that the present invention is not limited to the above embodiments, and the above embodiments and descriptions are only for explaining the principles of the present invention. Without departing from the spirit and scope of the present invention, the present invention may have various changes and improvements, which fall within the scope of the present invention to be protected. The scope of protection of the present invention is defined by the attached claims and their equivalents.

Claims

1. A general simulation method for complex particle two-phase flow, characterized by: Including: 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; S2: Based on the grids divided in S1, use the method of splitting finite elements based on characteristic lines 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; S3: Determine the interaction region between the particles and the flow field according to the geometric information of the particles. Combining the flow field velocity and flow field pressure calculated in step S2, calculate the force on the particle boundary, the force on the particle surface, and the torque acting on the particle through the IBM method; And construct a flow field velocity update formula, and update the flow field velocity at the next moment according to the flow field velocity update formula; S4: Use the force on the particle surface and the torque acting on the particle 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.

2. A general simulation method for particle two-phase flow according to claim 1, characterized in that: S1 includes: S1.1: Divide the flow field area into grids, denoted as D mesh ; S1.2: Use X r Indicates point X bk At the particle boundary Γ b The collection on: X r =A k X r +(I-A k )X bk ,k=1,2,…,N; Among them, A k is the rotation matrix used to fit the boundaries of all particles Boundaries of particles By N segment boundaries Γ i , (i=1,2,…,N); S1.3: Based on S1.2, solve the centroid of particles with arbitrary shapes in two-dimensional cases: Among them, X r are 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; S1.4: Determine the particle center of gravity X at time n c Grid D mesh Unit, with this unit as the center, the rest of the units contained in the square or cube with a side length of 2dh, and save the grid node coordinates x of these units j .

3. A general simulation method for complex particle two-phase flow according to claim 2, characterized in that: The grid D divided in S1 mesh Define node velocity on and pressure Then the flow field velocity of S2 is and flow field pressure p n+1 Solve the following equations in turn: Where: N u and N p are the finite element basis functions defined in velocity space and pressure space respectively; dt is the time step; M u is the mass matrix of velocity; is the predicted flow field velocity at the node; is the matrix of convection terms; η is the dynamic viscosity; K u is the stiffness matrix related to the flow field velocity; B u n is the second-order derivative of the flow field velocity; K p is the stiffness matrix of pressure; G is the term related to flow field pressure; H is the coupling term between flow field velocity and flow field pressure; is the second-order derivative of the flow field pressure; is the corrected flow field velocity; and p n +1 are the corrected flow field velocity and flow field pressure for the next time step respectively.

4. A general simulation method for complex particle two-phase flow according to claim 3, characterized in that: The force on the surface of the particle in S3 is: Among them, x j is the grid node coordinate, X r is the coordinate of the particle boundary point, F lr (X r ,t) is the particle boundary point X r The force at 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; Combined with the center of gravity of the particle S1, the torque T acting on the particle p for:

5. A general simulation method for complex particle two-phase flow according to claim 4, characterized in that: S4 includes: S4.1: Construct the six-degree-of-freedom model of the particle as: Where F is F in S3 p , M is the particle mass, U c ,ω,Q,X c They are particle velocity, particle angular velocity, particle Euler angle and particle displacement; I p represents the moment of inertia of the particle, and Represents the coordinate transformation matrix between coordinate systems; S4.2: Use the fourth-order Runge-Kutta method to solve the discrete form of the 6DOF model: Where:

6. A general simulation method for complex particle two-phase flow according to claim 5, characterized in that: The flow field velocity update formula in S4 is: Where x is the grid node coordinate, X r is the coordinate of the particle boundary point, F lr (X r ,t) is the particle boundary point X r The force at n+1 (x, t) represents the flow field velocity at time n+1, represents the flow field velocity at pseudo time n+1, is the force transmitted by the particle to the flow field.

Images