An analytical simulation method and system for gas-liquid-solid three-phase flow

Abstract

The application relates to a kind of gas-liquid-solid three-phase flow analytical simulation method and system, method includes: on Euler grid, establish gas, liquid, solid three-phase volume fraction field, meet volume constraint requirement;Solid phase volume fraction field is carried out fluid volume method advection transport, gas-liquid interface is terminated at solid boundary;The momentum source of immersed boundary is added to momentum equation and participates in pressure-velocity coupling solution, obtains flow field velocity and pressure;According to the fluid stress of revised and immersed boundary momentum source, the hydrodynamic force and torque suffered by each particle are calculated, as external force is applied to discrete element particle, and the particle motion and contact solution are promoted;Output three-phase flow field, phase distribution, particle motion state and force information.The application improves the prediction accuracy of solid into water, wetting, three-phase contact line and other problems, avoids the phase fraction out-of-bound and mass loss when two-phase VOF and discrete particle are mixedly described, and expands the applicability in liquid drop, bubble and particle interaction scene.

Description

Technical Field

[0001] This invention relates to the field of computational fluid dynamics, specifically to an analytical simulation method and system for gas-liquid-solid three-phase flow. Background Technology

[0002] Gas-liquid-solid three-phase flow is widely present in natural and industrial processes such as river and coastal sediment transport, slurry pipeline transportation, bubble column reactors, and hydraulic fracturing proppant transport. In terms of numerical simulation, existing technologies mainly develop along two routes: one is unanalyzed or semi-analyzed CFD-DEM, which treats particles as Lagrangian points and uses averaged governing equations for fluid flow, closed by empirical drag, lift, and buoyancy correlations; the other is analytical CFD-DEM, which directly analyzes the fluid around each particle on an Eulerian grid, with particle motion tracked by the DEM. In recent years, to capture the gas-liquid interface, researchers have further introduced the VOF (Volume of Fluid) method into CFD-DEM, forming DEM-VOF-like methods.

[0003] Sun & Sakai coupled DEM with VOF, solving fluid motion using CFD, capturing the gas-liquid interface with VOF, and tracking solid particles using DEM; the coupling between solid and fluid employed volume averaging and introduced an empirical drag correlation. Pozzetti & Peters proposed a multi-scale DEM-VOF coupling method, employing a dual-grid strategy with both volumetric and fine-scale fluid measurements to achieve a balance between accuracy and efficiency.

[0004] Saeedipour et al. coupled VOF with the fictitious domain method to study droplet-particle interactions; the liquid-gas flow was described by VOF, and the particles were represented as a volume fraction field on the Eulerian grid, with solid boundary conditions implemented through a penalty term in the momentum equation. Davydzenka et al. developed the open-source solver cfdemSolverIB_VOF, implementing analytical IB / fictitious domain-VOF CFD-DEM within the OpenFOAM / CFDEM framework for multiphase flows where deformable interfaces coexist with solid particles, but did not introduce an independent solid phase fraction field in the Eulerian framework to physically shield the gas-liquid interface. Wagner & Higgs III used the fictitious domain immersion boundary method to achieve fully analytical CFD-DEM-VOF coupling and proposed an extended VOF-IB algorithm to apply contact angle conditions on the three-phase contact line.

[0005] Lai, Zhao, Zhao & Huang introduced the signed distance field (SDF) into a fully analytical CFD-DEM to handle multiphase fluids and irregularly shaped particles, but did not disclose the use of the solid phase as an independent VOF phase within a unified, conserved three-phase framework of gas and liquid phases. Patel et al. and O'Brien & Bussmann respectively proposed methods for coupling VOF with IBM to handle contact line dynamics.

[0006] However, the aforementioned existing technologies still have the following drawbacks:

[0007] (1) The solid phase is not treated as an independent VOF phase. Existing DEM-VOF methods usually treat solid particles as discrete Lagrangian points or local porosity. The volume fraction of the gas-liquid two phases is still calculated as a two-phase VOF. An independent solid phase fraction field parallel to gas and liquid is not introduced, which leads to an incomplete three-phase volume conservation equation. Non-physical phase fraction overruns or mass defects are prone to occur near the interface.

[0008] (2) Multiphase flow interfaces are prone to penetrating into the particle interior. Traditional penetrating boundary methods only force the solid domain velocity through momentum sources and do not provide explicit volume occupancy for the solid phase on the Eulerian grid; the gas-liquid VOF interface lacks effective shielding near the particle surface and is prone to penetrating into the particle interior due to numerical diffusion, insufficient IBM velocity reconstruction or interface advection error, resulting in distortion of solid entry into water, wetting, and three-phase contact line position prediction.

[0009] (3) Fluid-particle coupling is mostly unanalyzed or semi-analyzed, relying on empirical correlations. Traditional CFD-DEM calculates the force of fluid on particles through empirical drag and lift correlations, which has limited predictive ability for irregular particles, strong particle-fluid-particle interactions, and flows near the interface.

[0010] (4) Immersion boundary momentum sources are mostly explicit penalties, which limit stability. Existing analytical IBM often uses explicit penalty terms to force the solid domain velocity to the particle velocity, and the time step is strictly limited by stability conditions; at the same time, hydrodynamic calculations often do not deduct the fictitious force applied by IBM inside the particle, resulting in particle force error.

[0011] (5) The Eulerian solid interface and the Lagrangian DEM surface are prone to misalignment. VOF transport will naturally diffuse the solid interface, while the DEM particle surface is accurately described by SDF; misalignment between the two can lead to non-physical flow, false surface tension or contact line position deviation at the solid-liquid / solid-gas interface.

[0012] (6) Lack of conservation treatment when particle volume changes. For scenarios where particle volume changes due to chemical expansion, dissolution or hygroscopic expansion, existing methods usually only modify the particle geometry or DEM radius, without simultaneously adding volume source terms to the VOF transport equation and fluid continuity equation, resulting in non-conservation of surrounding fluid volume, abnormal pressure field or interface distortion. Summary of the Invention

[0013] The purpose of this invention is to provide an analytical simulation method and system for gas-liquid-solid three-phase flow, in order to solve the problems in the prior art such as the lack of independent sites for the solid phase, the easy penetration of the multiphase flow interface into the particle interior, reliance on empirical correlations, poor stability of the immersion boundary method, and easy misalignment of the Eulerian-Lagrange interface.

[0014] To achieve the above objectives, the present invention adopts the following technical solution:

[0015] An analytical simulation method for gas-liquid-solid three-phase flow, including

[0016] Obtain computational domain mesh, particle shape sign distance field, fluid and solid material properties, and discrete element contact parameters;

[0017] A three-phase volume fraction field for gas, liquid, and solid phases is established on an Eulerian grid, where the solid phase volume fraction field is... With liquid phase volume fraction field Gas phase volume fraction field satisfy The volume constraints are used to calculate the mixture density ρ and dynamic viscosity μ by weighting the three-phase volume fractions.

[0018] The solid phase volume fraction field of each Euler grid element is calculated based on the sign distance field of the discrete element particles. The solid-occupied areas are explicitly marked on the Euler grid;

[0019] For solid phase volume fraction field Advection transport using the fluid volume method confines the sum of the gas-liquid two-phase volume fraction fields to 1- Within this range, the gas-liquid interface terminates at the solid boundary;

[0020] The density ρ and dynamic viscosity μ of the mixture are updated based on the three-phase volume fraction field, and the density-weighted solid phase fraction is calculated. and target speed Constructing an immersion boundary momentum source ,in The diagonal coefficients of the discrete momentum equation are used to incorporate the immersion boundary momentum source into the momentum equation and participate in the pressure-velocity coupling solution to obtain the flow field velocity and pressure.

[0021] The hydrodynamic forces and torques acting on each particle are calculated based on the corrected fluid stress and the momentum source at the immersion boundary.

[0022] Hydrodynamics and torque are applied as external forces to discrete element particles to drive particle motion and contact solutions.

[0023] Update the particle position, attitude, and sign distance field, and output the three-phase flow field, phase distribution, particle motion state, and force information.

[0024] Furthermore, the solid phase volume fraction field of each Eulerian grid element is calculated based on the sign distance field of the discrete element particles. The steps include:

[0025] The preliminary estimate of the solid phase fraction is obtained by counting the number of grid vertices located inside the particles;

[0026] For mesh elements partially cut by particles, a volume-weighted correction is performed using vertex symbol distance values, calculated as follows: The summation process iterates through the unit vertices. This is the symbolic distance value at the grid vertex after taking the minimum value among all particles.

[0027] Furthermore, in the solid phase volume fraction field Before or during the process of fluid volumetric advection transport, a solid-phase interface flux correction step is also included:

[0028] Solid phase fraction calculated from the symbolic distance field With the current solid phase volume fraction field deviation Construction correction speed ,in Here, ε represents the interface weights, and ε is a small amount to prevent division by zero; projecting the corrected velocity yields the corrected flux that satisfies the divergence-free condition. The modified flux is superimposed on the fluid volume method translation flux to drive the solid phase volume fraction field. Align with the solid surface defined by the symbolic distance field.

[0029] Furthermore, in the step of calculating the hydrodynamic force and torque on each particle based on the corrected fluid stress and the momentum source of the immersion boundary, the hydrodynamic force on particle i is... and torque By integrating over the Euler mesh elements of the covering particles, we obtain:

[0030] ;

[0031] ;

[0032] in For fluid stress tensor, Let i be the set of elements that intersect with particle i. For the unit volume, -(1- ) The term is used to subtract the fictitious force exerted by the immersion boundary momentum source inside the particle.

[0033] Furthermore, it also includes a step for conserving particle volume changes:

[0034] When the particle volume changes, calculate the particle volume change rate. The solid volume source is uniformly distributed according to the number of elements in the Euler grid cells intersecting with the particle. ,in The number of units occupied by particle i;

[0035] solid volume source Add a global velocity divergence source to the transport equation based on the volume fraction of the solid phase, and simultaneously add it to the continuity equation or the pressure Poisson equation. ,in To calculate the total volume of the domain, the solid volume changes by displacing or drawing in the surrounding fluid in equal amounts.

[0036] Furthermore, it also includes a three-phase contact line contact angle correction step:

[0037] When calculating the gas-liquid interface normal vector, the mesh surface crossed by the solid-phase interface is identified, and the calculation is based on the solid-phase surface normal vector. With a given contact angle Corrected gas-liquid interface normal vector To satisfy For dynamic contact angle, it depends on the tangential velocity of the contact line relative to the solid surface. At the angle of advance With back angle Interpolation between: ,in This is the static contact angle. The characteristic velocity scale.

[0038] Furthermore, the liquid phase volume fraction is transported by the fluid volume method equation with an interfacial compressibility term: ,in The relative compression velocity; the volume fraction of the solid phase is transported by the pure advection equation: ,in The source term is due to changes in particle volume; the gas phase volume fraction is determined by... Calculate and crop the interval [0,1] to eliminate rounding errors.

[0039] Furthermore, the local solid velocity of the particles Translational velocity of discrete element particles With angular velocity calculate: For deformable particles, The value is obtained by interpolation of the finite element deformation field inside the particle at a local location.

[0040] An analytical simulation system for gas-liquid-solid three-phase flow includes:

[0041] The input layer is used to obtain the computational domain mesh, particle shape sign distance field, fluid and solid material properties, and discrete element contact parameters.

[0042] The computational fluid dynamics solution layer is used to establish the volume fraction fields of the gas, liquid, and solid phases on an Eulerian grid, where the volume fraction field of the solid phase is... With liquid phase volume fraction field Gas phase volume fraction field satisfy The volume constraints are applied, and the mixture density ρ and dynamic viscosity μ are calculated using a three-phase volume fraction weighting. The solid phase volume fraction field for each Eulerian grid element is calculated based on the sign distance field of the discrete element particles. Explicitly identify the solid-occupied region on the Eulerian grid; for the solid phase volume fraction field Adjacent flow transport using the fluid volume method confines the sum of the gas-liquid two-phase volume fraction fields to 1- Within this range, the gas-liquid interface terminates at the solid boundary;

[0043] The coupling layer is used to update the mixture density ρ and dynamic viscosity μ based on the three-phase volume fraction field, and to calculate the density-weighted solid phase fraction. and target speed Constructing an immersion boundary momentum source ,in The diagonal coefficients of the discrete momentum equation are used to incorporate the immersion boundary momentum source into the momentum equation and participate in the pressure-velocity coupling solution to obtain the flow field velocity and pressure. The hydrodynamic force and torque on each particle are calculated based on the corrected fluid stress and the immersion boundary momentum source.

[0044] The discrete element solution layer is used to apply the hydrodynamic forces and torques as external forces to the discrete element particles, thereby promoting particle motion and contact solution.

[0045] The output layer is used to update the particle position, attitude, and sign distance field, and output the three-phase flow field, phase distribution, particle motion state, and force information.

[0046] A storage medium storing a computer program, which, when executed by a processor, implements the analytical simulation method for the gas-liquid-solid three-phase flow described above.

[0047] In summary, the present invention has the following advantages:

[0048] (1) Explicit solid occupancy avoids the multiphase flow interface from penetrating into the particle interior. This is achieved by introducing an independent solid phase fractional field. It provides physical sites for the solid phase on the Eulerian grid; the gas-liquid interface terminates naturally at the solid phase boundary, avoiding the non-physical phenomenon of interface penetration into particles in traditional IBM, thereby improving the prediction accuracy of problems such as solid entry into water, wetting, and three-phase contact lines.

[0049] (2) Three-phase unified conservation. By introducing an independent solid phase fraction field, the gas, liquid and solid phases share a unified definition of mixture properties and volume constraints, avoiding the phase fraction overshoot and mass loss that may occur when describing the mixture of two-phase VOF and discrete particles.

[0050] (3) Fully analytical fluid-particle interaction. The flow field around the particles is directly analyzed on the Eulerian grid. Combined with DEM, the particle contact and irregular shape are accurately handled. It gets rid of the dependence on empirical drag and lift correlation formulas and is suitable for non-spherical, polydisperse particles and strong interaction scenarios.

[0051] (4) Implicit IBM improves stability and accuracy. The momentum source enters the PIMPLE iteration implicitly, allowing for a time step that is compatible with the fluid convection stability; IBM fictitious forces are subtracted in the hydrodynamic calculation, making the force on the particles more accurate.

[0052] (5) The solid-phase interface is geometrically consistent with the DEM surface. The flux correction algorithm actively suppresses VOF diffusion, keeps the Eulerian solid boundary sharp, and is consistent with the Lagrangian DEM surface, thereby improving the accuracy of the flow field and contact line position near the interface.

[0053] (6) Supports particle volume change conservation. The volume change of expanding / contracting particles is processed synchronously through solid VOF source terms and fluid divergence source terms to ensure the conservation of the total volume of the three phases and the physical consistency of the pressure field and interface evolution.

[0054] (7) It can handle three-phase contact angles and wetting behavior. By modifying the gas-liquid interface normal vector, static or dynamic contact angle conditions can be applied to the surface of moving particles, which expands the applicability of the method in scenarios involving the interaction of droplets, bubbles and particles. Attached Figure Description

[0055] Figure 1 The diagram shows a comparison between the traditional immersion boundary method and the solid-phase field occupancy method of the present invention, where (a) is the traditional immersion boundary method and (b) is the present invention.

[0056] Figure 2 This is a block diagram of the module composition of the system of the present invention.

[0057] Figure 3 This is a flowchart illustrating the method of the present invention.

[0058] Figure 4 This is a schematic diagram of the calculation of solid phase fraction based on SDF.

[0059] Figure 5 This is a schematic diagram of flux correction at the solid-phase interface.

[0060] Figure 6 This is a schematic diagram illustrating the conservation of particle volume change.

[0061] Figure 7 This is a schematic diagram of three-phase contact angle correction. Detailed Implementation

[0062] The present invention will now be described in further detail.

[0063] This invention addresses the numerical simulation of three-phase flows involving discrete solid particles, continuous liquids, and continuous gases. Typical applications include: particle entry / exit at the water-gas interface, bubble movement through a particle bed, movement of expanding / contracting particles in a liquid-containing environment, and transport of irregular particles in slurry or bubble columns.

[0064] Traditional immersion boundary methods for handling such problems typically apply only virtual forces to the momentum equation to force the solid domain velocity to equal the particle velocity, without providing independent volume occupancy for the solid phase on the Eulerian grid. As a result, the VOF field describing the gas-liquid interface may "immerse" into the particle interior during interfacial advection and IBM velocity reconstruction: the gas or liquid phase volume fraction is non-zero in the region where the particle should occupy, leading to incorrect descriptions of key physical phenomena such as solid ingress into water, particle wetting, and the location of the three-phase contact line.

[0065] This invention introduces an independent solid phase volume fraction field. Explicitly mark the regions occupied by solids on the Eulerian grid; the volume fractions of the gas and liquid phases satisfy... Under constraints, the flow is confined outside the solid domain, thus naturally preventing the multiphase flow interface from entering the particle interior. Simultaneously, the three-phase VOF method, the immersion boundary method, and the SDF-based DEM are organically coupled to form a method and system for analytically simulating gas-liquid-solid three-phase flow on a fixed Eulerian grid. Among these, the solid phase volume fraction field... It is not a simple label or mask, but an independent physical variable that participates in the conservation of momentum, conservation of mass, and interface evolution. Its introduction changes the mathematical structure of the entire set of governing equations, rather than just a post-processing modification.

[0066] An analytical simulation system for gas-liquid-solid three-phase flow includes:

[0067] The input layer is used to obtain the computational domain mesh, particle shape sign distance field, fluid and solid material properties, and discrete element contact parameters.

[0068] The computational fluid dynamics solution layer is used to establish the volume fraction fields of the gas, liquid, and solid phases on an Eulerian grid, where the volume fraction field of the solid phase is... With liquid phase volume fraction field Gas phase volume fraction field satisfy The volume constraints are applied, and the mixture density ρ and dynamic viscosity μ are calculated using a three-phase volume fraction weighting. The solid phase volume fraction field for each Eulerian grid element is calculated based on the sign distance field of the discrete element particles. Explicitly identify the solid-occupied region on the Eulerian grid; for the solid phase volume fraction field Adjacent flow transport using the fluid volume method confines the sum of the gas-liquid two-phase volume fraction fields to 1- Within this range, the gas-liquid interface terminates at the solid boundary;

[0069] The coupling layer is used to update the mixture density ρ and dynamic viscosity μ based on the three-phase volume fraction field, and to calculate the density-weighted solid phase fraction. and target speed Constructing an immersion boundary momentum source ,in The diagonal coefficients of the discrete momentum equation are used to incorporate the immersion boundary momentum source into the momentum equation and participate in the pressure-velocity coupling solution to obtain the flow field velocity and pressure. The hydrodynamic force and torque on each particle are calculated based on the corrected fluid stress and the immersion boundary momentum source.

[0070] The discrete element solution layer is used to apply the hydrodynamic forces and torques as external forces to the discrete element particles, thereby promoting particle motion and contact solution.

[0071] The output layer is used to update the particle position, attitude, and sign distance field, and output the three-phase flow field, phase distribution, particle motion state, and force information.

[0072] Specifically, the analytical simulation method for gas-liquid-solid three-phase flow includes the following steps:

[0073] Step 1: At the beginning of each fluid time step, map the current position and orientation of the DEM particles to the CFD background mesh to establish the affiliation relationship between the mesh nodes and the particles.

[0074] Step 2: Calculate the required solid phase fraction for each CFD mesh element based on the particle SDF. The phase fraction explicitly identifies the solid-occupied region on the Eulerian mesh, providing a physically shielded boundary for the multiphase flow interface. For cells partially cut by particles, a volume-weighted formula based on vertex SDF values ​​is used to improve accuracy.

[0075] Step 3: Apply fractional field to the solid phase Perform VOF advection; if flux correction is enabled, then perform advection according to... and The bias constructs a corrected flux.

[0076] Step 4: Drive the process using a flux correction algorithm. Towards Alignment was achieved to suppress interfacial diffusion caused by VOF transport; the volume fractions of the gas and liquid phases were then resynchronized after correction to meet the requirements. .because By occupying the solid domain, the gas-liquid interface is confined outside the solid boundary, preventing the interface from penetrating into the particle's interior.

[0077] Step 5: Update the mixture density based on the three-phase volume fractions With dynamic viscosity ; Calculate density-weighted solid phase fraction and target speed .

[0078] Step 6: Construct an implicit IBM momentum source ,in The diagonal coefficients of the discrete momentum equation; Incorporate the momentum equation and participate in the PIMPLE pressure-velocity coupling solution.

[0079] Step 7: Solve for the velocity and pressure fields in the PIMPLE iteration to obtain the flow field that satisfies the incompressibility condition (including the source term of particle volume change).

[0080] Step 8: Based on the corrected fluid stress and IBM source term, calculate the hydrodynamic force and torque acting on each particle: , .

[0081] Step 9: Apply hydrodynamic force as an external force to the DEM particles, and use a DEM sub-cycle with a time step much smaller than the fluid time step to drive particle motion and contact solution.

[0082] Step 10: If the particles undergo volume changes (expansion or contraction), calculate the rate of volume change and... Adding a volume source to the transport equation Simultaneously, a velocity divergence source is added to the continuity equation / pressure Poisson equation. This causes the solid volume to change in an equal amount, displacing or drawing in the surrounding fluid.

[0083] Step 11: Update particle position, attitude and SDF, and output three-phase flow field, phase distribution, particle motion state and force information, and proceed to the next fluid time step.

[0084] This invention has the following innovations:

[0085] (1) Three-phase VOF unified framework: explicit occupancy of solid-state field and interface shielding

[0086] In the continuous medium region, liquid phase volume fraction is defined. Gas phase volume fraction With solid volume fraction ,satisfy ;

[0087] The density and dynamic viscosity of a mixture are defined as follows: , ;

[0088] in , , Densities for liquid, gas, and solid, respectively. , , This corresponds to the dynamic viscosity. Solid phase viscosity. Choose a sufficiently large value to suppress spurious flow inside the particles.

[0089] Unlike existing IBM-VOF methods, this invention will... As with , Parallel independent phase fractional fields participate in transport and conservation. The volume occupied by solid particles was explicitly marked on the Eulerian grid, thus limiting the sum of the volume fractions of the gas and liquid phases to 1- Within. Thus, when the gas-liquid interface adsorbs to the vicinity of the solid boundary, due to... The area has been occupied. and The natural tendency is towards 0, and the interface will not allow non-physical penetration inside the particles.

[0090] Liquid volume fraction is transported by the VOF equation with interfacial compressibility terms:

[0091] ;

[0092] in The relative compression velocity is used to maintain a sharp gas-liquid interface. The solid volume fraction is transported by the pure advection equation:

[0093] ;

[0094] in This is the source term caused by changes in particle volume. The gas phase volume fraction is determined by... Calculate and crop the interval [0,1] to eliminate rounding errors.

[0095] (2) Calculation of solid phase fraction based on SDF

[0096] Each solid particle is defined by a signed distance field. describe, >0 indicates the interior of the particle. <0 indicates external, =0 represents the surface. For each CFD mesh element, the number of vertices located inside the particles is first counted to obtain a preliminary estimate of the solid phase fraction; for elements partially cut by particles, volume-weighted correction is further performed using the vertex SDF values.

[0097] ;

[0098] The summation process iterates through the unit vertices. This is the minimum signed distance over all particles at that vertex. This formula provides a more accurate solid volume fraction than simple node counting, even when the particle size is only slightly larger than the mesh size, thus ensuring the accurate placement of the solid field in the particle geometry.

[0099] (3) Implicit momentum source immersion boundary method

[0100] Define density-weighted solid phase fraction

[0101] ;

[0102] and target speed

[0103] ;

[0104] IBM's momentum source is

[0105] ;

[0106] in These are the diagonal coefficients of the discrete momentum equation, reflecting the implicit feedback strength of the local momentum equation to velocity. Choosing these coefficients ensures that the IBM source term and the momentum equation converge consistently within the same iteration step, without needing to further reduce the time step. Because... The source term is implicitly introduced into the velocity-pressure coupling by iteratively solving within PIMPLE, avoiding the strict time step restrictions imposed by explicit penalty methods.

[0107] Local solid velocity of particles Calculated from the translational and angular velocities of DEM particles:

[0108] ;

[0109] For deformable particles It can be obtained by interpolation at local locations based on the finite element deformation field inside the particle.

[0110] (4) Modified hydrodynamic transmission

[0111] The hydrodynamic forces and torques acting on the particles are obtained by integrating over the CFD mesh elements covering the particles:

[0112] ;

[0113] ;

[0114] in For fluid stress tensor, Let i be the set of elements that intersect with particle i. Let be the unit volume. Where -(1- ) This term is used to subtract the fictitious forces applied by IBM inside the particles, improving the accuracy of hydrodynamic calculations.

[0115] (5) Solid-phase interface flux correction

[0116] VOF advection diffuses at the solid interface. Define the solid fraction deviation.

[0117] ;

[0118] Constructing modified velocities near the solid interface

[0119] ;

[0120] in For interface weight, To minimize the amount of flux and avoid division by zero, the correction rate is limited to an upper limit determined by the local mesh size and the maximum allowable Courant number. The corrected flux satisfying the divergence-free condition is obtained by projecting u* onto the solution of the pressure Poisson equation. And add it to the VOF horizontal flux, thereby driving Towards Alignment.

[0121] (6) Source terms of particle volume change

[0122] For particle i whose volume changes with time, let its volume change rate be... First, identify all CFD mesh elements that intersect with the particle, and then uniformly distribute the solid volume source among these elements:

[0123] ;

[0124] in Let be the number of units occupied by particle i. To ensure the overall volume is conserved, a uniform velocity divergence source is added to the continuity equation:

[0125] ;

[0126] in This is for the total volume of the computational domain. The source term is treated as an extra term in the pressure Poisson equation, causing expanding particles to displace an equal amount of surrounding fluid and contracting particles to draw in an equal amount of fluid.

[0127] (7) Correction of contact angle of three-phase contact line

[0128] Calculating the normal vector of the gas-liquid interface At that time, identify the mesh surfaces traversed by the solid-state interface. For such surfaces, based on the solid-state surface normal vector... With a given static or dynamic contact angle , correct Make it satisfy

[0129] ;

[0130] The dynamic contact angle can be modeled using a hysteresis model: based on the tangential velocity of the contact line relative to the solid surface. At the forward angle With back angle Interpolation between:

[0131] ;

[0132] in The characteristic velocity scale. Corrected. Used for CSF surface tension calculations to achieve correct wetting behavior at the three-phase contact line.

[0133] The key technical parameters are shown in the table below:

[0134]

[0135] The present invention also provides a storage medium storing a computer program, which, when executed by a processor, implements the above-described analytical simulation method for gas-liquid-solid three-phase flow.

[0136] The storage medium described in this embodiment can be a disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), USB flash drive, portable hard drive, etc.

[0137] Example 1: Rigid sphere passes through the gas-liquid interface

[0138] The computational domain is a rectangular tank, filled with liquid at the bottom and gas at the top, with an initial horizontal gas-liquid interface. A tank with a diameter of... A rigid sphere is released vertically downwards with a given initial velocity, allowing it to pass through the gas-liquid interface and enter the liquid. The particles are tracked by a DEM, with their SDF being spherical; the CFD background mesh is fixed, employing the three-phase VOF framework, implicit IBM, and solid-phase interface flux correction of this invention.

[0139] because Explicit occupancy on the Eulerian grid ensures that the gas-liquid interface terminates naturally at the surface of the sphere, without any interface penetration into the interior of the sphere. Verification results show that the sphere's trajectory, interface deformation morphology, and the position of the three-phase contact line are consistent with experimental observations and literature reference solutions. The solid-phase interface remains sharp during traversal, with no significant numerical diffusion. Quantifiable verification metrics include: the relative error between the particle centroid trajectory and the reference solution, the interface penetration depth (which should be close to 0), and the solid-phase interface thickness (which should be maintained within 2–4 grid cells).

[0140] Example 2: Verification of volume conservation of expanding particles in a gas-liquid environment

[0141] Let an initial radius be Spherical particles are placed in a domain containing a horizontal gas-liquid interface, and the particle radius follows a sinusoidal law. Expansion and contraction. Using the particle volume source term processing method of the present invention, in... Add to the equation Add to the Poisson equation for pressure The sum of the volume fractions of the gas, liquid, and solid phases, as well as the total fluid volume, were monitored throughout the computational domain. The results showed that both remained constant at 1 throughout the simulation (within the limits of machine accuracy), verifying volume conservation. Quantifiable verification metrics included: the maximum deviation of the sum of the three-phase volume fractions from 1, and the integral value of the particle volume change rate. The integral values ​​are consistent, and there are no abnormal oscillations in the pressure field.

[0142] Example 3: Underwater entry of irregular particles

[0143] A hyperellipsoidal SDF was used to describe non-spherical particles, which were injected into a still water body with a certain initial horizontal velocity. The interparticle and particle-wall contacts were processed using a DEM, and the flow field around the particles and the evolution of the free surface were analyzed using CFD. Due to the solid-phase field... By shielding the particle interior, the gas-liquid interface does not penetrate into the particle interior, enabling accurate capture of splashing, cavity formation, and particle trajectory during water entry. This invention eliminates the need to generate a body-fitted mesh for each particle, naturally handling irregular shapes and polydisperse particle groups. It is suitable for simulating the water entry, wetting, and cavity formation processes of non-spherical objects such as soil particles and catalyst particles. Quantifiable verification metrics include: comparison of water entry trajectory with experimental / reference solutions, cavity depth and closure time, and the degree of interface adhesion on the particle surface.

[0144] Example of parameter settings: The following parameters are for implementation reference only and do not constitute a limitation on the invention.

[0145] fluid density kg / m³, kg / m³;

[0146] Hydrodynamic viscosity = Pa·s, = Pa·s;

[0147] solid density kg / m³;

[0148] Solid effective viscosity = Pa·s;

[0149] Surface tension coefficient N / m;

[0150] static contact angle ;

[0151] Flux correction relaxation factor relax = 0.5;

[0152] The maximum correction Courant number is maxAlignCo = 0.5;

[0153] The particle volume source tolerance is volumeRelTol = 0.05.

[0154] The above embodiments are preferred embodiments of the present invention, but the embodiments of the present invention are not limited to the above embodiments. Any changes, modifications, substitutions, combinations, or simplifications made without departing from the spirit and principle of the present invention shall be considered equivalent substitutions and shall be included within the protection scope of the present invention.

Claims

1. An analytical simulation method for gas-liquid-solid three-phase flow, characterized in that, include Obtain computational domain mesh, particle shape sign distance field, fluid and solid material properties, and discrete element contact parameters; A three-phase volume fraction field for gas, liquid, and solid phases is established on an Eulerian grid, where the solid phase volume fraction field is... With liquid phase volume fraction field Gas phase volume fraction field satisfy The volume constraints are used to calculate the mixture density ρ and dynamic viscosity μ by weighting the three-phase volume fractions. The solid phase volume fraction field of each Euler grid element is calculated based on the sign distance field of the discrete element particles. The solid-occupied areas are explicitly marked on the Euler grid; For solid phase volume fraction field Adjacent flow transport using the fluid volume method confines the sum of the gas-liquid two-phase volume fraction fields to 1- Within this range, the gas-liquid interface terminates at the solid boundary; The density ρ and dynamic viscosity μ of the mixture are updated based on the three-phase volume fraction field, and the density-weighted solid phase fraction is calculated. and target speed Constructing an immersion boundary momentum source ,in The diagonal coefficients of the discrete momentum equation are used to incorporate the immersion boundary momentum source into the momentum equation and participate in the pressure-velocity coupling solution to obtain the flow field velocity and pressure. The hydrodynamic forces and torques acting on each particle are calculated based on the corrected fluid stress and the momentum source at the immersion boundary. Hydrodynamics and torque are applied as external forces to discrete element particles to drive particle motion and contact solutions. Update the particle position, attitude, and sign distance field, and output the three-phase flow field, phase distribution, particle motion state, and force information.

2. The analytical simulation method for gas-liquid-solid three-phase flow according to claim 1, characterized in that, The solid phase volume fraction field of each Euler grid element is calculated based on the sign distance field of the discrete element particles. The steps include: The preliminary estimate of the solid phase fraction is obtained by counting the number of grid vertices located inside the particles; For mesh elements partially cut by particles, a volume-weighted correction is performed using vertex symbol distance values, calculated as follows: The summation process iterates through the unit vertices. This is the symbolic distance value at the grid vertex after taking the minimum value among all particles.

3. The analytical simulation method for gas-liquid-solid three-phase flow according to claim 1, characterized in that, In the solid phase volume fraction field Before or during the process of fluid volumetric advection transport, a solid-phase interface flux correction step is also included: Solid phase fraction calculated from the symbolic distance field With the current solid phase volume fraction field deviation Construction correction speed ,in Here, ε represents the interface weights, and ε is a small amount to prevent division by zero; projecting the corrected velocity yields the corrected flux that satisfies the divergence-free condition. The modified flux is superimposed on the fluid volume method translation flux to drive the solid phase volume fraction field. Align with the solid surface defined by the symbolic distance field.

4. The analytical simulation method for gas-liquid-solid three-phase flow according to claim 1, characterized in that, In the step of calculating the hydrodynamic force and torque on each particle based on the corrected fluid stress and the momentum source of the immersion boundary, the hydrodynamic force on particle i is... and torque By integrating over the Euler mesh elements of the covering particles, we obtain: ； ； in For fluid stress tensor, Let i be the set of elements that intersect with particle i. For the unit volume, -(1- ) The term is used to subtract the fictitious force exerted by the immersion boundary momentum source inside the particle.

5. The analytical simulation method for gas-liquid-solid three-phase flow according to claim 1, characterized in that, It also includes a step for conserving particle volume changes: When the particle volume changes, calculate the particle volume change rate. The solid volume source is uniformly distributed according to the number of elements in the Euler grid cells intersecting with the particle. ,in The number of units occupied by particle i; solid volume source Add a global velocity divergence source to the transport equation based on the volume fraction of the solid phase, and simultaneously add it to the continuity equation or the pressure Poisson equation. ,in To calculate the total volume of the domain, the solid volume changes by displacing or drawing in the surrounding fluid in equal amounts.

6. The analytical simulation method for gas-liquid-solid three-phase flow according to claim 1, characterized in that, It also includes the three-phase contact wire contact angle correction steps: When calculating the gas-liquid interface normal vector, the mesh surface crossed by the solid-phase interface is identified, and the calculation is based on the solid-phase surface normal vector. With a given contact angle Corrected gas-liquid interface normal vector To satisfy ; For dynamic contact angle, it depends on the tangential velocity of the contact line relative to the solid surface. At the angle of advance With back angle Interpolation between: ,in This is the static contact angle. The characteristic velocity scale.

7. The analytical simulation method for gas-liquid-solid three-phase flow according to claim 1, characterized in that, The liquid phase volume fraction is transported by the fluid volume method equation with interfacial compressibility terms: ,in This refers to the relative compression speed; The volume fraction of the solid phase is transported by the pure advection equation: ,in The source term is due to changes in particle volume; the gas phase volume fraction is determined by... Calculate and crop the interval [0,1] to eliminate rounding errors.

8. The analytical simulation method for gas-liquid-solid three-phase flow according to claim 1, characterized in that, The local solid velocity of the particles Translational velocity of discrete element particles With angular velocity calculate: ; For deformable particles The value is obtained by interpolation of the finite element deformation field inside the particle at a local location.

9. An analytical simulation system for gas-liquid-solid three-phase flow, used to implement the analytical simulation method for gas-liquid-solid three-phase flow as described in any one of claims 1-8, characterized in that, include: The input layer is used to obtain the computational domain mesh, particle shape sign distance field, fluid and solid material properties, and discrete element contact parameters. The computational fluid dynamics solution layer is used to establish the volume fraction fields of the gas, liquid, and solid phases on an Eulerian grid, where the volume fraction field of the solid phase is... With liquid phase volume fraction field Gas phase volume fraction field satisfy The volume constraints are applied, and the mixture density ρ and dynamic viscosity μ are calculated using a three-phase volume fraction weighting. The solid phase volume fraction field for each Eulerian grid element is calculated based on the sign distance field of the discrete element particles. Explicitly identify the solid-occupied region on the Eulerian grid; for the solid phase volume fraction field Adjacent flow transport using the fluid volume method confines the sum of the gas-liquid two-phase volume fraction fields to 1- Within this range, the gas-liquid interface terminates at the solid boundary; The coupling layer is used to update the mixture density ρ and dynamic viscosity μ based on the three-phase volume fraction field, and to calculate the density-weighted solid phase fraction. and target speed Constructing an immersion boundary momentum source ,in The diagonal coefficients of the discrete momentum equation are used to incorporate the immersion boundary momentum source into the momentum equation and participate in the pressure-velocity coupling solution to obtain the flow field velocity and pressure. The hydrodynamic force and torque on each particle are calculated based on the corrected fluid stress and the immersion boundary momentum source. The discrete element solution layer is used to apply the hydrodynamic forces and torques as external forces to the discrete element particles, thereby promoting particle motion and contact solution. The output layer is used to update the particle position, attitude, and sign distance field, and output the three-phase flow field, phase distribution, particle motion state, and force information.

10. A storage medium, characterized in that, The storage medium stores a computer program, which, when executed by a processor, implements the analytical simulation method for gas-liquid-solid three-phase flow as described in any one of claims 1-8.

Images