A numerical simulation method and system for gas-solid two-phase flow

Abstract

This invention discloses a numerical simulation method and system for gas-solid two-phase flow, applicable to the field of numerical simulation in fluid mechanics. The method includes: acquiring flow field variables comprising a fluid phase and a particle phase; constructing a preprocessing matrix to compress the eigenvalue differences of the control equations to obtain preprocessed control equations; numerically solving the preprocessed control equations to update the fluid phase flow field variables; calculating the particle motion state based on the updated fluid phase flow field variables; generating gas-solid coupling source terms based on the particle motion state and adding them to the preprocessed control equations; repeating the above iterations until the flow field variables converge to obtain the numerical simulation results. The numerical simulation method for gas-solid two-phase flow provided by this invention eliminates numerical rigidity at low Mach numbers through a preprocessing matrix, achieving stable and efficient solutions, while simultaneously realizing bidirectional coupling between the gas and solid phases.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology in fluid mechanics, and in particular to a numerical simulation method and system for gas-solid two-phase flow. Background Technology

[0002] Gas-solid two-phase flow is widely used in engineering fields such as chemical engineering, energy, metallurgy, and environmental protection. Numerical simulation technology is of great significance for the design and optimization of related equipment. When the particle concentration transitions from a rarefied state to a dense state, the compressible multiphase particle mesh method, by introducing a gas-solid energy coupling mechanism and a phase discontinuity handling strategy, achieves effective simulation of particle two-phase flow under compressible conditions, significantly improving physical fidelity and applicability.

[0003] However, existing technologies, when dealing with low-speed flow scenarios, still rely on the governing equations framework constructed under compressible conditions and directly employ density-based solvers for numerical solutions within this framework. Since the sound velocity term dominates at low Mach numbers, while the flow term is relatively weak, the eigenvalues ​​of the equation set vary significantly, resulting in a substantial increase in numerical rigidity. This makes it difficult to maintain stability in the solution process, leading to a sharp decline in convergence efficiency, and in severe cases, even computational divergence, making it impossible to obtain effective simulation results. Summary of the Invention

[0004] This invention provides a numerical simulation method and system for gas-solid two-phase flow to solve the technical problem of numerical rigidity caused by the significant difference in the characteristic values ​​of the governing equations in the numerical simulation of gas-solid two-phase flow under low Mach number conditions, so as to achieve stable and efficient numerical simulation of gas-solid two-phase flow under low-speed conditions.

[0005] To address the aforementioned technical problems, embodiments of the present invention provide a numerical simulation method for gas-solid two-phase flow, comprising: The flow field variables of gas-solid two-phase flow under low Mach number conditions are obtained, including fluid phase flow field variables and particulate phase flow field variables. Entering the gas-solid coupling iterative mode, a preprocessing matrix is ​​constructed based on the flow field variables at the current moment, and a preprocessing control equation is obtained based on the preprocessing matrix. The preprocessing matrix is ​​used to compress the eigenvalue differences of the preprocessing control equation. The preprocessing control equations are numerically solved, and the fluid phase flow field variables at the current moment are updated based on the solution results. The motion state of the particles corresponding to the particle phase flow field variables is obtained based on the updated fluid phase flow field variables; Based on the motion state of the particles and the fluid phase flow field variables, a gas-solid coupling source term is generated and added to the preprocessing control equation; Repeat the gas-solid coupling iterative mode until the flow field variables meet the preset convergence conditions to obtain the numerical simulation results of the gas-solid two-phase flow.

[0006] As one preferred embodiment, the process of constructing the preprocessing matrix includes: The equivalent sound velocity is calculated based on the flow field variables, and the sound velocity term in the preprocessing control equation is corrected using the equivalent sound velocity. The preprocessing matrix is ​​generated based on the corrected sound speed term.

[0007] As one preferred embodiment, the step of numerically solving the preprocessing control equations and updating the fluid phase flow field variables at the current moment based on the solution results includes... The convection terms in the preprocessed control equations are spatially discretized, and the interface flux is calculated using a flux splitting scheme. The viscous flux is obtained by calculating the viscous term in the preprocessing control equation. Based on the interface flux and the viscous flux, update the fluid phase flow field variables at the current moment.

[0008] As one preferred embodiment, the generation of gas-solid coupling source terms based on the particle motion state and the fluid phase flow field variables includes: The particle motion equations are solved based on the particle motion state and the fluid phase flow field variables, and the positions of each particle are updated based on the solution results. The particle phase flow field variables are calculated based on the updated particle positions. Based on the particle phase flow field variables and the fluid phase flow field variables, the momentum exchange between the particles and the fluid is calculated; The momentum exchange quantity is allocated to the preprocessing control equation to obtain the gas-solid coupling source term.

[0009] As one preferred embodiment, solving the particle motion equations based on the particle's motion state and the fluid phase flow field variables includes: The position and velocity of the particle at the current moment are obtained based on the particle's motion state. The resultant force on the particle is calculated based on the particle's position and velocity and the fluid phase flow field variables; The acceleration of the particle is calculated based on the resultant force, and the acceleration is integrated to obtain the position of the particle at the next moment.

[0010] Another embodiment of the present invention provides a numerical simulation system for gas-solid two-phase flow, comprising: The flow field variable acquisition module is used to acquire the flow field variables of gas-solid two-phase flow under low Mach number conditions. The flow field variables include fluid phase flow field variables and particulate phase flow field variables. The control equation generation module is used to enter the gas-solid coupling iterative mode, construct a preprocessing matrix based on the flow field variables at the current moment, and obtain the preprocessed control equation based on the preprocessing matrix. The preprocessing matrix is ​​used to compress the eigenvalue differences of the preprocessed control equation. The numerical solution module is used to numerically solve the preprocessed control equations and update the fluid phase flow field variables at the current moment based on the solution results. The particle motion solution module is used to obtain the motion state of the particles corresponding to the updated fluid phase flow field variables based on the updated fluid phase flow field variables. The coupling source term generation module is used to generate gas-solid coupling source terms based on the motion state of the particles and the fluid phase flow field variables, and add the gas-solid coupling source terms to the preprocessing control equations; The iterative control module is used to repeatedly execute the gas-solid coupling iterative mode until the flow field variables meet the preset convergence conditions, and obtain the numerical simulation calculation results of the gas-solid two-phase flow.

[0011] As one preferred embodiment, the governing equation generation module includes: A sound velocity correction unit is used to calculate the equivalent sound velocity based on the flow field variables, and to use the equivalent sound velocity to correct the sound velocity term in the preprocessing control equation. A preprocessing matrix construction unit is used to generate the preprocessing matrix based on the corrected sound speed term.

[0012] As one preferred embodiment, the numerical solution module includes: The convection term discretization unit is used to spatially discretize the convection terms in the preprocessing control equations and calculate the interface flux using a flux splitting scheme. The viscosity term calculation unit is used to calculate the viscosity term in the preprocessing control equation to obtain the viscous flux; The flow field update unit is used to update the fluid phase flow field variables at the current moment based on the interface flux and the viscous flux.

[0013] As one preferred embodiment, the coupling source term generation module includes: The particle position update unit is used to solve the particle motion equation based on the particle's motion state and the fluid phase flow field variables, and update the position of each particle based on the solution result; The particle phase field statistics unit is used to calculate the particle phase flow field variables based on the updated position of the particles. The momentum exchange calculation unit is used to calculate the momentum exchange between the particles and the fluid based on the particle phase flow field variables and the fluid phase flow field variables. The source term generation unit is used to allocate the momentum exchange quantity to the preprocessing control equation to obtain the gas-solid coupling source term.

[0014] As one preferred embodiment, the particle position update unit is further configured to: The position and velocity of the particle at the current moment are obtained based on the particle's motion state. The resultant force on the particle is calculated based on the particle's position and velocity and the fluid phase flow field variables; The acceleration of the particle is calculated based on the resultant force, and the acceleration is integrated to obtain the position of the particle at the next moment.

[0015] Compared with the prior art, the beneficial effects of the embodiments of the present invention are at least one of the following: (1) This invention initializes the flow field variables of the gas-solid two-phase flow and then enters a gas-solid coupling iterative mode: a preprocessing matrix is ​​constructed based on the current flow field variables to compress the eigenvalue differences of the control equation, resulting in a preprocessed control equation; the equation is numerically solved to update the fluid phase flow field variables; the particle motion state is calculated based on the updated fluid phase flow field variables, and a gas-solid coupling source term is generated based on this motion state and the fluid phase flow field variables, which is then added to the preprocessed control equation; the above iterative process is repeated until the flow field variables meet the convergence condition. This invention effectively solves the problem of instability caused by numerical rigidity under low Mach number conditions by compressing eigenvalue differences through the preprocessing matrix, and can achieve stable convergence while maintaining computational accuracy.

[0016] (2) This invention integrates preprocessing matrix construction, numerical solution, particle motion calculation, and gas-solid coupling source term feedback into the same iterative framework, realizing bidirectional coupled calculation of fluid and particle phases. Compared with the existing technology that follows the compressible framework and directly uses the density solver, this invention can avoid numerical divergence caused by the dominance of the sound velocity term under low Mach number conditions, significantly improving computational stability and convergence efficiency. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating a numerical simulation method for gas-solid two-phase flow in one embodiment of the present invention. Figure 2 This is a schematic diagram of a numerical simulation system for gas-solid two-phase flow in one embodiment of the present invention.

[0018] Figure label: Among them, 11 is the flow field variable acquisition module, 12 is the control equation generation module, 13 is the numerical solution module, 14 is the particle motion solution module, 15 is the coupled source term generation module, and 16 is the iterative control module. Detailed Implementation

[0019] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The purpose of providing these embodiments is to make the disclosure of the present invention more thorough and comprehensive. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0020] In the description of this application, the terms "first," "second," "third," etc., are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Therefore, a feature defined with "first," "second," "third," etc., may explicitly or implicitly include one or more of that feature. In the description of this application, unless otherwise stated, "a plurality of" means two or more.

[0021] In the description of this application, it should be noted that, unless otherwise expressly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to fixed connections, detachable connections, or integral connections; they can refer to mechanical connections or electrical connections; they can refer to direct connections or indirect connections through an intermediate medium; and they can refer to the internal communication between two components. The terms "vertical," "horizontal," "left," "right," "upper," "lower," and similar expressions used herein are for illustrative purposes only and do not indicate or imply that the device or component referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as limiting the invention. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items. Those skilled in the art can understand the specific meaning of the above terms in this application based on the specific circumstances.

[0022] In the description of this application, it should be noted that, unless otherwise defined, all technical and scientific terms used in this invention have the same meaning as commonly understood by one of ordinary skill in the art. The terminology used in this specification is for the purpose of describing specific embodiments only and is not intended to limit the invention. Those skilled in the art can understand the specific meaning of the above terms in this application based on the specific circumstances.

[0023] One embodiment of the present invention provides a numerical simulation method for gas-solid two-phase flow. For details, please refer to [link to documentation]. Figure 1 , Figure 1 The diagram shown is a flowchart illustrating a numerical simulation method for gas-solid two-phase flow according to one embodiment of the present invention, which includes steps S1 to S6: S1: Obtain the flow field variables of gas-solid two-phase flow under low Mach number conditions. The flow field variables include fluid phase flow field variables and particulate phase flow field variables. Step S1 is the starting point of the entire numerical simulation method. Its purpose is to establish a digital model of the computational domain in the computer, providing an initial state for subsequent gas-solid coupling iterative calculations. Low Mach number conditions refer to flow states where the flow Mach number is much less than 1. Typical scenarios include powder mixing chambers, microreactor mixing chambers, and porous media pore units. In these scenarios, the fluid velocity is much lower than the speed of sound, causing traditional solvers based on the compressibility assumption to struggle to converge stably due to significant differences in eigenvalues. This embodiment is designed specifically for this scenario, thus explicitly defining the flow being processed as a low Mach number condition in the initial stage.

[0024] In this context, fluid phase flow field variables refer to the fluid state parameters described on the Eulerian grid. These specifically include: pressure p, velocity vector (u,v), temperature T, density ρ, and fluid properties such as dynamic viscosity μ, gas constant R, and specific heat ratio γ. These variables are defined on each grid cell, characterizing the instantaneous state of the fluid within that cell. In numerical simulations, these fluid phase flow field variables are the direct solution objects of the governing equations.

[0025] In this scheme, the particulate phase flow field variable specifically refers to the particle volume fraction α. p The particle volume fraction refers to the ratio of the volume occupied by all particles within each grid cell to the volume of that grid cell, ranging from 0 to 1. This variable is also defined on an Eulerian grid, but its source differs from that of the fluid phase—it is not obtained directly by solving equations, but rather through statistical analysis of the positions of individual particles. Specifically, after tracing the position of each particle using the Lagrangian method, the total volume of particles falling within each grid cell is calculated, and this volume is divided by the volume of that grid cell to obtain the particle volume fraction for that cell.

[0026] The distinction between these two types of variables arises because in gas-solid two-phase flow, the fluid and particles exist under different descriptive frameworks. The fluid is described using an Eulerian framework (observing changes in physical quantities on a fixed grid), while particles are described using a Lagrangian framework (tracing the trajectory of each particle). The particle volume fraction serves as a bridge connecting these two frameworks: it transforms discrete particle information into continuous field information, enabling the particle phase to influence the governing equations of the fluid phase as source terms.

[0027] Specifically, the system first generates a computational mesh based on the geometric parameters input by the user. Taking the top cover drive cavity model as an example, the cavity length is L and the height is H. The computer divides it into N... x Multiply by Nᵧ rectangular cells. Each grid cell has a unique number, and the grid node coordinates are determined by linear interpolation. The density of the grid affects the calculation accuracy and efficiency, and users can preset the number of grid cells according to their actual needs.

[0028] After generating the mesh, the system assigns fluid property parameters to each mesh cell. These parameters include: initial pressure p0, initial temperature T0, initial density ρ0, dynamic viscosity μ, gas constant R, and specific heat ratio γ. Simultaneously, the system verifies the consistency of the parameters according to the gas law p=ρRT, ensuring that the initial state satisfies the thermodynamic equilibrium condition. The dynamic viscosity μ is stored as a global constant for subsequent viscous flux calculations.

[0029] Setting boundary conditions is a crucial step in initialization. The system identifies mesh boundary elements and applies a constant velocity U_lid to the upper boundary (top cover), parallel to the long side of the cavity. No-slip boundary conditions, i.e., zero velocity, are applied to the lower, left, and right boundaries. These boundary conditions constitute the external excitation that drives the fluid within the cavity to form a complex vortex structure.

[0030] Through the three processing actions described above, the system outputs the initial flow field variables, including fluid phase flow field variables (pressure distribution, velocity distribution, temperature distribution, density distribution) and particulate phase flow field variables (particle volume fraction distribution), stored in memory as a multidimensional array, with each array index corresponding to a grid number. This data is used to subsequently construct the preprocessing matrix, and the flow field variables can be obtained in ways not limited to the initialization operations described above; stored flow field data can also be read directly from external data files.

[0031] S2: Enter the gas-solid coupling iterative mode, construct a preprocessing matrix based on the flow field variables at the current moment, and obtain the preprocessing control equation based on the preprocessing matrix. The preprocessing matrix is ​​used to compress the eigenvalue differences of the preprocessing control equation. After initialization, the system has acquired the current flow field variables. The next step is to ensure the density-based solver operates stably under low Mach number conditions. The difficulty in low Mach number flow lies in the governing equations themselves. The eigenvalues ​​of the original compressible Navier-Stokes equations are uc, u, and u+c, where c is the speed of sound. Under low-speed conditions, the fluid velocity u is very small, but the speed of sound c is very large, differing by several orders of magnitude. This results in significant differences in the eigenvalues ​​of the equations, leading to extremely strong numerical rigidity. This embodiment addresses this problem by introducing low Mach number preprocessing techniques.

[0032] Preferably, in one embodiment of the present invention, the process of constructing the preprocessing matrix includes: The equivalent sound velocity is calculated based on the flow field variables, and the sound velocity term in the preprocessed control equation is corrected using the equivalent sound velocity. A preprocessing matrix is ​​generated based on the corrected sound speed term.

[0033] The essence of preprocessing is to modify the time derivative terms of the equations. Specifically, a preprocessing matrix is ​​constructed, and the time derivative terms are multiplied by this matrix on the left to transform the original equations into a preprocessed form. The core function of this preprocessing matrix is ​​to compress the eigenvalues ​​of the equation system from {uc, u, u+c} to {u-c̃, u, u+c̃}, where c̃ is the processed equivalent speed of sound.

[0034] Specifically, the equivalent sound velocity is first calculated based on the flow field variables at the current moment. The equivalent sound velocity is defined as c̃ = √min(c², V_ref²), where c is the local sound velocity and V_ref is a user-preset reference velocity. The reference velocity is usually set to be on the same order of magnitude as the characteristic flow velocity; for example, in the case of a top cover driving a square cavity, the velocity of the top cover can be taken. When the flow Mach number is extremely low, c̃ ≈ V_ref, and the equivalent sound velocity is no longer much greater than the flow velocity, thus compressing the difference in eigenvalues.

[0035] The sound speed term in the governing equation is corrected using an equivalent sound speed, and this correction is reflected in the piecewise function of the reference velocity. The reference velocity U_{r,i} is defined as follows: when the local flow velocity is less than ε times the sound speed, U_{r,i} takes the value of ε·a_i; when the local flow velocity is between ε times the sound speed and the sound speed, U_{r,i} takes the local flow velocity; when the local flow velocity is greater than the sound speed, U_{r,i} takes the sound speed. Here, ε is taken as 10⁻. 5 where a_i is the speed of sound. This piecewise function ensures that preprocessing plays a dominant role in the low-speed region, gradually degenerating into the original compressible format in the high-speed region.

[0036] A preprocessing matrix is ​​generated based on the corrected sound velocity term, and the original governing equations are transformed into a preprocessed form. The transformed equation is J·∂P / ∂t +∇·F_c -∇·F_v = D_p, where J is the preprocessing matrix (i.e., the Jacobian matrix from the conserved variables to the original variables), P is the original variable vector (p,u,v,T), F_c and F_v are the convection and viscous fluxes, respectively, and D_p is the source term.

[0037] After the above processing, the preprocessing matrix compresses the eigenvalue differences of the equation system from the order of O(1 / Ma) to the order of O(1), thus eliminating numerical rigidity.

[0038] It should be noted that the preprocessing matrix constructed in step S2 is not static. During subsequent iterations, this matrix is ​​dynamically updated based on the current flow field variable state to ensure that the preprocessing effect always matches the flow field. This iterative coupling mechanism is key to achieving stability in extremely low Mach numbers or high shear regions in this embodiment. Secondly, the above preprocessing operation only changes the time derivative term of the equation and does not affect the steady-state solution. That is, once the calculation converges, the effect of the preprocessing matrix automatically disappears, and the obtained solution is consistent with the steady-state solution of the original equation.

[0039] S3: Numerically solve the preprocessing control equations and update the fluid phase flow field variables at the current moment based on the solution results; Step S3 essentially involves solving the governing equations to calculate the pressure, velocity, and temperature distributions at the next moment, given the current flow field state. After step S2, the numerical rigidity in the original governing equations has been eliminated, resulting in preprocessed governing equations. However, these equations are still in continuous partial differential form, which the computer cannot directly process. The next step is to discretize these continuous equations into algebraic equations and solve them to obtain the updated flow field variables for each grid cell.

[0040] Preferably, in one embodiment of the present invention, the preprocessing control equations are numerically solved, and the fluid phase flow field variables at the current moment are updated based on the solution results, including... The convection terms in the preprocessed control equations are spatially discretized, and the interface flux is calculated using a flux splitting scheme. The viscous flux is obtained by calculating the viscous term in the preprocessing governing equations. Update the fluid phase flow field variables at the current moment based on interfacial flux and viscous flux.

[0041] Specifically, the convection terms in the preprocessing control equations are spatially discretized, and the interface flux is calculated using a flux splitting scheme. The purpose of spatial discretization is to transform continuous partial differential equations into algebraic equations that can be processed by a computer. Specifically, the system divides the computational domain into several grid cells and calculates the flux at each grid interface. In this embodiment, a flux splitting scheme is used, which decomposes the convection flux into two parts: one part propagating upstream and the other downstream, and then processes them separately according to the flow direction.

[0042] In practice, the fluxes on the interface are calculated using the flow field variables on both sides of the interface as input. To improve calculation accuracy, this embodiment also uses the MUSCL method to extrapolate the higher-order reconstructed values ​​on both sides of the interface from the cell center values. The interface Mach number is defined as a combination function of the Mach numbers on both sides. The flux splitting scheme divides the numerical flux into convective flux and pressure flux, and its expression is: in, For interface quality traffic, For the quantity of flow, This refers to pressure flux. The formula for calculating pressure flux introduces a pressure correction term related to the Mach number, which is expressed as: in, It is an adjustable parameter. and Let the pressure splitting function be... , These are the Mach numbers for the left and right sides of the interface. The average sound velocity at the interface. This pressure correction term is a core feature of the AUSM⁺-UP scheme, specifically designed for low Mach number flows, and can effectively suppress non-physical numerical dissipation and pressure oscillations generated by traditional compressible flux schemes under low-speed conditions.

[0043] The viscous term in the preprocessing governing equations is calculated to obtain the viscous flux, which originates from the shear stress within the fluid and is driven by the velocity gradient. This embodiment extends the original inviscid Euler equations to the Navier-Stokes equations, which include viscous terms. The calculation of the viscous flux depends on the fluid dynamic viscosity μ and the gradient of the velocity field. The system traverses all mesh interfaces and calculates the viscous flux at each interface based on the current velocity distribution of the flow field. This enables the invention to handle two-phase flow problems involving viscous particles, significantly expanding its applicability.

[0044] The fluid phase flow field variables at the current moment are updated based on interfacial flux and viscous flux. The system substitutes the calculated convective interfacial flux and viscous flux into the discretized control equations, and obtains the flow field variables at the next moment through time progression. Specifically, the pressure, velocity, temperature, and density of each grid cell are updated using the finite difference method. For the non-conservative term -p∇α in the gas-solid two-phase flow control equations... p In this embodiment, the non-conservation term is corrected based on the particle volume fraction, and the corrected non-conservation term is added as an independent source term to the right-hand side of the momentum equation, thereby achieving stable calculation of gas-solid two-phase coupling. The necessity of this correction lies in the fact that the particle phase does not participate in pressure propagation but affects momentum exchange, and directly applying the treatment method for gas-liquid two-phase flow would lead to inaccurate calculation of the non-conservation term.

[0045] Through the above processing, the system outputs updated fluid phase flow variables, including pressure distribution, velocity distribution, temperature distribution, and density distribution. These variables will serve as inputs for subsequent steps, including time progression and particle phase calculations.

[0046] S4: Obtain the motion state of the particles corresponding to the particle phase flow field variables based on the updated fluid phase flow field variables; After updating the fluid phase flow variables in step S3, obtaining only the evolution information of the fluid phase is insufficient for numerical simulation of gas-solid two-phase flow. The particulate phase moves under the influence of fluid drag and pressure gradient forces in the flow field, and the presence of particles also affects the flow pattern of the fluid phase through momentum exchange. Therefore, after completing the fluid phase calculations, it is necessary to further solve for the motion state of the particulate phase.

[0047] This embodiment uses the Lagrange method to track individual particles. Specifically, each particle is treated as an independent computational object, and its position and velocity are recorded over time. This can accurately describe the particle's trajectory and its interaction with the fluid.

[0048] Specifically, the system obtains the particle motion state and fluid phase flow field variables at the current moment. The particle motion state includes the position coordinates and velocity components of each particle, and the fluid phase flow field variables are obtained from the calculation results in step S3. Since particles are usually located inside grid cells rather than exactly on grid nodes, the system needs to interpolate the fluid velocity, pressure, and other physical quantities on the surrounding grid nodes based on the particle's current location to obtain the fluid parameter values ​​at the particle's location. The interpolation method can also use bilinear interpolation or higher-order interpolation formats.

[0049] The net force acting on a particle is calculated as the particle moves through a fluid and is subjected to various forces. This embodiment considers drag, gravity, Magnus force, Saffman force, and virtual mass force. Drag originates from friction and pressure difference between the fluid and particle surfaces; its direction is opposite to the particle's motion relative to the fluid, and its magnitude depends on relative velocity, particle diameter, and fluid viscosity. Gravity is the volume force exerted by the Earth's gravitational field, and its direction is vertically downward. Magnus force originates from the velocity difference on both sides caused by particle rotation. Saffman force originates from the shear lift generated by the velocity gradient in the flow field. Virtual mass force originates from the additional inertial effect generated when the particle accelerates, causing the surrounding fluid to accelerate as well. The calculation formulas for these forces adopt the standard form for particle two-phase flow, and the specific coefficients are determined based on the particle Reynolds number and flow field characteristics. After calculating each force sequentially, the system vector-superimposes them to obtain the net force acting on the particle.

[0050] The system solves the particle's motion equations based on the resultant force, updating the particle's position and velocity. The essence of these equations is Newton's second law: the mass of a particle multiplied by its acceleration equals the net force acting on it. Starting from the particle's current position and velocity, the system calculates the acceleration based on the resultant force, then integrates the acceleration over time to obtain the velocity at the next time step, and finally integrates the velocity over time to obtain the position at the next time step. This embodiment uses a first-order explicit Euler method for time integration. This method requires only one calculation per time step, is simple to implement, has low computational cost, and meets the accuracy requirements of engineering applications.

[0051] Through the above steps, the system outputs updated particle positions and velocities. This data will be used as input to S5 to statistically analyze particle phase flow field variables and generate gas-solid coupling source terms.

[0052] S5: Generate gas-solid coupling source terms based on particle motion state and fluid phase flow field variables, and add the gas-solid coupling source terms to the preprocessing control equations; Step S4 completes the calculation of the particle motion state, obtaining the updated position and velocity of each particle. However, at this point, the information of the particle phase remains at the discrete individual level, while the governing equations of the fluid phase are solved on an Eulerian grid; the two are not under the same descriptive framework. In order for the particle phase to influence the fluid phase, the discrete particle information must be converted into continuous field information on the grid and added to the fluid phase governing equations as source terms. Step S5's task is to complete this conversion.

[0053] Preferably, in one embodiment of the present invention, generating gas-solid coupling source terms based on the particle motion state and fluid phase flow field variables includes: The particle motion equations are solved based on the particle motion state and fluid phase flow field variables, and the positions of each particle are updated based on the solution results. Calculate the particle phase flow field variables based on the updated particle positions; Based on the particulate phase flow field variables and the fluid phase flow field variables, the momentum exchange between particles and fluid is calculated. The momentum exchange quantity is allocated to the preprocessing control equations to obtain the gas-solid coupling source term.

[0054] Specifically, the particle motion equations are solved based on the particle motion state and fluid phase flow field variables, and the positions of each particle are updated based on the solution results. The specific implementation of this step has been detailed in step S4. The second step is to calculate the particle phase flow field variables based on the updated particle positions. In this embodiment, the particle phase flow field variables specifically refer to the particle volume fraction α. p The volume fraction is the ratio of the volume occupied by a particle within each grid cell to the volume of that grid cell, ranging from 0 to 1. The system iterates through all particles, determining which grid cell each particle falls into based on its position, and adding the particle's volume to the corresponding grid cell. After the iteration is complete, the total volume of particles in each grid cell is divided by the volume of that grid cell to obtain the particle volume fraction for that cell. This step only uses the particle's position information and does not require its velocity. Particles with the same velocity and particles with different velocities contribute the same amount to the volume fraction.

[0055] Based on the particulate and fluid phase flow field variables, the momentum exchange between particles and fluid is calculated. This momentum exchange originates from the drag force, gravity, Magnus force, Suffman force, and virtual mass force acting on the particles. According to Newton's third law, the reaction force exerted by the fluid on the particles is the force exerted by the particles on the fluid. The system distributes the resultant force on each particle calculated in S4 to its corresponding grid cell according to the particle's location, obtaining the momentum exchange on each grid cell. The relative position of the particles to the grid cells must be considered during the distribution; a distance-weighted method is typically used.

[0056] The momentum exchange is allocated to the preprocessed control equations to obtain the gas-solid coupling source term. This embodiment modifies the non-conservation terms in the control equations to address the specific characteristics of gas-solid two-phase flow. The original control equations contain a term -p∇α. p This term originates from the influence of the particle volume fraction gradient on fluid momentum. In gas-solid two-phase flow, the particle phase does not participate in pressure propagation but does affect momentum exchange. Directly applying the treatment methods for gas-liquid two-phase flow would lead to inaccurate calculations of non-conservation terms. Therefore, this embodiment uses the particle volume fraction α at the current moment... p The gradient of the non-conservative term is used to correct it. The corrected non-conservative term, together with the momentum exchange, constitutes the gas-solid coupling source term, which is stored in the source term array of the corresponding grid cell and then added to the right-hand side of the preprocessed governing equation. In subsequent iterative calculations, this source term will participate in the time advancement of the momentum equation, thereby achieving reverse coupling between the particles and the fluid phase.

[0057] Preferably, in one embodiment of the present invention, solving the particle motion equations based on the particle's motion state and the fluid phase flow field variables includes: The position and velocity of the particle at the current moment are obtained based on the particle's motion state; The resultant force on the particle is calculated based on the particle's position and velocity, as well as the fluid phase flow field variables. The particle's acceleration is calculated based on the resultant force, and the acceleration is integrated to obtain the particle's position at the next moment.

[0058] S6: Repeat the gas-solid coupling iterative mode until the flow field variables meet the preset convergence conditions to obtain the numerical simulation results of the gas-solid two-phase flow.

[0059] Steps S2 to S5 constitute a complete iterative calculation process for gas-solid coupling. A single execution of this process completes the calculation within one time step: starting from the flow field variables at the current moment, it sequentially undergoes preprocessing, fluid phase solution, particle motion calculation, source term generation and feedback, resulting in updated fluid phase flow field variables and particle states.

[0060] However, the calculation results of a single time step are far from meeting the requirements of engineering analysis. For unsteady flow problems, it is necessary to start from the initial moment and advance sequentially until the preset simulation duration is reached. For the solution within each time step, due to the nonlinearity of the equations and the strong coupling between the two phases, multiple sub-iterations are usually required to obtain a convergent solution that meets the accuracy requirements. Step S6 controls the entire time advancement process to ensure that the calculation is completed under stable and convergent conditions, and finally outputs the calculation results.

[0061] In this embodiment, the loop execution and convergence determination are implemented using a dual-time-step advancement strategy.

[0062] Dual-timestep simulation discretizes time into several physical time steps, and then executes several virtual time sub-iterations within each physical time step. The physical time steps correspond to the actual flow's time progression, and their step size is set by the user according to the simulation requirements. The virtual time sub-iterations are internal loops introduced to obtain a convergent solution within each physical time step; their step size only serves the purpose of iterative convergence and has no physical meaning. The advantage of this dual-timestep structure is that it allows for larger physical time steps in unsteady simulations, while ensuring computational accuracy within each time step through virtual time sub-iterations, eliminating the impact of preprocessing operations on time accuracy.

[0063] In each virtual time sub-iteration of the physical time step, the preprocessing matrix is ​​dynamically updated and the interface flux is recalculated. One of the core innovations of this embodiment lies in the iterative coupling mechanism between the preprocessing matrix and the interface flux. In traditional methods, the preprocessing matrix is ​​calculated once at the beginning of the time step based on the initial flow field and remains unchanged throughout the entire time step. When the flow field changes rapidly, the initially calculated preprocessing matrix gradually deviates from the actual flow field state, leading to a mismatch between eigenvalue correction and flux calculation, continuous accumulation of numerical errors, and decreased computational stability.

[0064] To address this issue, this embodiment performs the following operations in each virtual time sub-iteration: First, the preprocessing matrix is ​​recalculated based on the flow field state of the current iteration step, including recalculating the equivalent sound velocity and preprocessing parameters; then, based on the newly constructed preprocessing matrix, the convection and pressure fluxes on all mesh interfaces are recalculated. Through the iterative cycle of "flow field update - preprocessing matrix update - flux recalculation - flow field update," the preprocessing matrix is ​​ensured to remain consistent with the current flow field state.

[0065] A third-order Runge-Kutta method is employed for virtual time advancement. Within each virtual time sub-iteration, this scheme uses the third-order Runge-Kutta method for time advancement. This method offers third-order time accuracy, good stability, and allows for a larger virtual time step size while maintaining computational accuracy. Specifically, three sub-steps are executed within each sub-iteration step. Each sub-step updates the flow field variables based on the current residual, and the updated result for that sub-iteration step is obtained after all three sub-steps are completed.

[0066] The convergence criteria are determined to decide whether to continue the loop. In this embodiment, the convergence determination includes two levels of conditions: The first level is the convergence of virtual time sub-iterations. When the change in the flow field variables is less than a preset threshold, or when the preset maximum number of sub-iterations is reached, the calculation of the current physical time step is considered to have converged, the sub-iteration within that physical time step ends, and the next physical time step begins. The second level is the termination of the overall time progression. When the physical time step length progresses to the preset simulation duration, or when the user actively terminates the calculation, the entire loop process ends.

[0067] The dual-time-step strategy, combined with dynamic updates to the preprocessing matrix, ensures both the temporal accuracy of unsteady flow simulations and good numerical stability in extremely low Mach numbers or high shear regions. Through this cyclic control mechanism, this embodiment can stably advance calculations under low Mach number conditions, ultimately yielding computational results suitable for engineering analysis.

[0068] When the convergence condition is met, the computer outputs the numerical simulation results of the gas-solid two-phase flow. These results include: fluid phase flow field variables (pressure field, velocity field, temperature field, density field) and particle phase states (particle position, velocity, particle volume fraction distribution) at each time step. These data can be further used for post-processing analyses such as flow field visualization, particle enrichment region identification, and mixing uniformity assessment.

[0069] The final calculation results can be used for: (1) Evaluate particle mixing uniformity. Calculate the mixing uniformity index by statistically analyzing the particle volume fraction within each grid cell to determine whether the particle distribution within the cavity is uniform. When the mixing uniformity index is lower than a set threshold, it indicates a mixing inhomogeneity problem, requiring adjustment of operating parameters.

[0070] (2) Identify particle deposition risk areas. Based on the spatiotemporal distribution of particle volume fraction, identify particle enrichment areas and deposition areas. These areas are usually flow dead zones or vortex core locations, and long-term deposition may lead to equipment blockage or decreased reaction efficiency. The calculation results can guide structural improvements or adjustments to operating strategies.

[0071] (3) Predicting energy loss. Based on the velocity field of the fluid phase and the motion state of the particle phase, the energy loss during stirring or conveying is calculated. The energy loss curve reflects the energy consumption level under different operating conditions, providing a basis for optimizing process parameters.

[0072] The above engineering applications demonstrate that the numerical simulation method provided by this invention can not only obtain calculation results of gas-solid two-phase flow under low Mach number conditions, but also transform these results into specific engineering decision-making basis, providing technical support for the design and optimization of equipment such as powder stirring chambers, microreactor mixing chambers, and porous media pore units.

[0073] Another embodiment of the present invention provides a numerical simulation system for gas-solid two-phase flow. For details, please refer to [link to documentation]. Figure 2 , Figure 2 The diagram shown illustrates a numerical simulation system for gas-solid two-phase flow according to one embodiment of the present invention, comprising: The flow field variable acquisition module 11 is used to acquire the flow field variables of gas-solid two-phase flow under low Mach number conditions. The flow field variables include fluid phase flow field variables and particulate phase flow field variables. The control equation generation module 12 is used to enter the gas-solid coupling iterative mode, construct a preprocessing matrix based on the flow field variables at the current moment, and obtain the preprocessed control equation based on the preprocessing matrix. The preprocessing matrix is ​​used to compress the eigenvalue differences of the preprocessed control equation. The numerical solution module 13 is used to numerically solve the preprocessed control equations and update the fluid phase flow field variables at the current moment based on the solution results. The particle motion solution module 14 is used to obtain the motion state of the particles corresponding to the updated fluid phase flow field variables based on the updated fluid phase flow field variables. The coupling source term generation module 15 is used to generate gas-solid coupling source terms based on the motion state of particles and fluid phase flow field variables, and add the gas-solid coupling source terms to the preprocessing control equations. The iterative control module 16 is used to repeatedly execute the gas-solid coupling iterative mode until the flow field variables meet the preset convergence conditions, and obtain the numerical simulation calculation results of the gas-solid two-phase flow.

[0074] Preferably, in one embodiment of the present invention, the governing equation generation module includes: The sound velocity correction unit is used to calculate the equivalent sound velocity based on the flow field variables and to correct the sound velocity term in the preprocessed control equation using the equivalent sound velocity. The preprocessing matrix construction unit is used to generate a preprocessing matrix based on the corrected sound velocity term.

[0075] Preferably, in one embodiment of the present invention, the numerical solution module includes: The convection term discrete element is used to spatially discretize the convection terms in the preprocessed control equations and calculate the interface flux using a flux splitting scheme. The viscosity term calculation unit is used to calculate the viscosity term in the preprocessing governing equations to obtain the viscous flux; The flow field update unit is used to update the fluid phase flow field variables at the current moment based on the interface flux and viscous flux.

[0076] Preferably, in one embodiment of the present invention, the coupling source term generation module includes: The particle position update unit is used to solve the particle motion equations based on the particle's motion state and fluid phase flow field variables, and update the position of each particle based on the solution results. The particle phase field statistics unit is used to calculate the particle phase flow field variables based on the updated particle positions. The momentum exchange calculation unit is used to calculate the momentum exchange between particles and fluid based on the particle phase flow field variables and the fluid phase flow field variables. The source term generation unit is used to allocate momentum exchange quantities to the preprocessing control equations to obtain gas-solid coupling source terms.

[0077] Preferably, in one embodiment of the present invention, the particle position update unit is further configured to: The position and velocity of the particle at the current moment are obtained based on the particle's motion state; The resultant force on the particle is calculated based on the particle's position and velocity, as well as the fluid phase flow field variables. The particle's acceleration is calculated based on the resultant force, and the acceleration is integrated to obtain the particle's position at the next moment.

[0078] Compared with the prior art, the beneficial effects of the embodiments of the present invention are at least one of the following: (1) In the gas-solid coupling iteration process, the present invention obtains the particle position by solving the particle motion equation, calculates the particle phase flow field variables based on the particle position, and then generates a gas-solid coupling source term to be added to the control equation, thereby realizing momentum feedback from the particle phase to the fluid phase. The present invention fully considers the momentum exchange mechanism between the gas and solid phases. Compared with the existing technology that ignores the particle feedback effect or only considers unidirectional coupling, it can more realistically reflect the interphase interaction in dense gas-solid two-phase flow and effectively improve the physical fidelity of numerical simulation.

[0079] (2) The embodiments of the present invention employ eigenvalue compression technology based on preprocessing matrices, enabling density-based solvers to operate stably under low Mach number conditions without the need to switch to pressure-based solvers or other low Mach number approximation methods. Compared to the limitations of existing technologies that require switching solution strategies across different Mach number ranges, the present invention broadens the applicability of a single solution framework and reduces the implementation complexity of numerical simulations.

[0080] The embodiments described above are merely illustrative of several implementations of the present invention, and while the descriptions are specific and detailed, they should not be construed as limiting the scope of the present invention. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of the present invention, and these modifications and improvements all fall within the scope of protection of the present invention. Therefore, the scope of protection of this patent should be determined by the appended claims.

Claims

1. A numerical simulation method for gas-solid two-phase flow, characterized in that, include: The flow field variables of gas-solid two-phase flow under low Mach number conditions are obtained, including fluid phase flow field variables and particulate phase flow field variables. Entering the gas-solid coupling iterative mode, a preprocessing matrix is ​​constructed based on the flow field variables at the current moment, and a preprocessing control equation is obtained based on the preprocessing matrix. The preprocessing matrix is ​​used to compress the eigenvalue differences of the preprocessing control equation. The preprocessing control equations are numerically solved, and the fluid phase flow field variables at the current moment are updated based on the solution results. The motion state of the particles corresponding to the particle phase flow field variables is obtained based on the updated fluid phase flow field variables; Based on the motion state of the particles and the fluid phase flow field variables, a gas-solid coupling source term is generated and added to the preprocessing control equation; Repeat the gas-solid coupling iterative mode until the flow field variables meet the preset convergence conditions to obtain the numerical simulation results of the gas-solid two-phase flow.

2. The numerical simulation method for gas-solid two-phase flow as described in claim 1, characterized in that, The process of constructing the preprocessing matrix includes: The equivalent sound velocity is calculated based on the flow field variables, and the sound velocity term in the preprocessing control equation is corrected using the equivalent sound velocity. The preprocessing matrix is ​​generated based on the corrected sound speed term.

3. The numerical simulation method for gas-solid two-phase flow as described in claim 1, characterized in that, The process involves numerically solving the preprocessing control equations and updating the fluid phase flow field variables at the current moment based on the solution results, including... The convection terms in the preprocessed control equations are spatially discretized, and the interface flux is calculated using a flux splitting scheme. The viscous flux is obtained by calculating the viscous term in the preprocessing control equation. Based on the interface flux and the viscous flux, update the fluid phase flow field variables at the current moment.

4. The numerical simulation method for gas-solid two-phase flow as described in claim 1, characterized in that, The generation of gas-solid coupling source terms based on the particle motion state and the fluid phase flow field variables includes: The particle motion equations are solved based on the particle motion state and the fluid phase flow field variables, and the positions of each particle are updated based on the solution results. The particle phase flow field variables are calculated based on the updated particle positions. Based on the particle phase flow field variables and the fluid phase flow field variables, the momentum exchange between the particles and the fluid is calculated; The momentum exchange quantity is allocated to the preprocessing control equation to obtain the gas-solid coupling source term.

5. The numerical simulation method for gas-solid two-phase flow as described in claim 4, characterized in that, Solving the particle motion equations based on the particle's motion state and the fluid phase flow field variables includes: The position and velocity of the particle at the current moment are obtained based on the particle's motion state. The resultant force on the particle is calculated based on the particle's position and velocity and the fluid phase flow field variables; The acceleration of the particle is calculated based on the resultant force, and the acceleration is integrated to obtain the position of the particle at the next moment.

6. A numerical simulation system for gas-solid two-phase flow, characterized in that, include: The flow field variable acquisition module is used to acquire the flow field variables of gas-solid two-phase flow under low Mach number conditions. The flow field variables include fluid phase flow field variables and particulate phase flow field variables. The control equation generation module is used to enter the gas-solid coupling iterative mode, construct a preprocessing matrix based on the flow field variables at the current moment, and obtain the preprocessed control equation based on the preprocessing matrix. The preprocessing matrix is ​​used to compress the eigenvalue differences of the preprocessed control equation. The numerical solution module is used to numerically solve the preprocessed control equations and update the fluid phase flow field variables at the current moment based on the solution results. The particle motion solution module is used to obtain the motion state of the particles corresponding to the updated fluid phase flow field variables based on the updated fluid phase flow field variables. The coupling source term generation module is used to generate gas-solid coupling source terms based on the motion state of the particles and the fluid phase flow field variables, and add the gas-solid coupling source terms to the preprocessing control equations; The iterative control module is used to repeatedly execute the gas-solid coupling iterative mode until the flow field variables meet the preset convergence conditions, and obtain the numerical simulation calculation results of the gas-solid two-phase flow.

7. The numerical simulation system for gas-solid two-phase flow as described in claim 6, characterized in that, The governing equation generation module includes: A sound velocity correction unit is used to calculate the equivalent sound velocity based on the flow field variables, and to use the equivalent sound velocity to correct the sound velocity term in the preprocessing control equation. A preprocessing matrix construction unit is used to generate the preprocessing matrix based on the corrected sound speed term.

8. The numerical simulation system for gas-solid two-phase flow as described in claim 6, characterized in that, The numerical solution module includes: The convection term discretization unit is used to spatially discretize the convection terms in the preprocessing control equations and calculate the interface flux using a flux splitting scheme. The viscosity term calculation unit is used to calculate the viscosity term in the preprocessing control equation to obtain the viscous flux; The flow field update unit is used to update the fluid phase flow field variables at the current moment based on the interface flux and the viscous flux.

9. The numerical simulation system for gas-solid two-phase flow as described in claim 6, characterized in that, The coupling source term generation module includes: The particle position update unit is used to solve the particle motion equation based on the particle's motion state and the fluid phase flow field variables, and update the position of each particle based on the solution result; The particle phase field statistics unit is used to calculate the particle phase flow field variables based on the updated position of the particles. The momentum exchange calculation unit is used to calculate the momentum exchange between the particles and the fluid based on the particle phase flow field variables and the fluid phase flow field variables. The source term generation unit is used to allocate the momentum exchange quantity to the preprocessing control equation to obtain the gas-solid coupling source term.

10. The numerical simulation system for gas-solid two-phase flow as described in claim 9, characterized in that, The particle position update unit is also used for: The position and velocity of the particle at the current moment are obtained based on the particle's motion state. The resultant force on the particle is calculated based on the particle's position and velocity and the fluid phase flow field variables; The acceleration of the particle is calculated based on the resultant force, and the acceleration is integrated to obtain the position of the particle at the next moment.

Images