A calculation method for multi-flow regime partition identification and dynamic update of DFNS based on VOF-CFD-DEM coupling

Abstract

This invention relates to the field of numerical simulation technology. Addressing the limitations of existing CFD-DEM coupled methods in simultaneously describing the entire process of sudden water inrush and the gas-liquid interaction and free surface evolution within the tunnel space, this invention proposes a computational method for multi-fluid regime identification and dynamic updating of DFNS based on VOF-CFD-DEM coupling. This method includes constructing a fluid-solid-gas multi-field coupled computational domain and a discretized model; establishing a unified DFNS control equation system; establishing a multi-fluid regime criterion system to generate fluid regime partition labels; formulating parameter mapping logic to adaptively switch the control equations; discretely solving the control equations and iteratively updating them at each time step; updating particle states and unit porosity; reconstructing and updating porosity and permeability; and outputting the optimal parameter combination. This invention can more realistically characterize key processes such as fluid regime evolution and reconstructing dominant permeability channels under sudden water inrush loading, improving the physical consistency and stability of simulations of water and mud inrush processes in fault fracture zones.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation technology, and in particular to a computational method for multi-fluid state partition identification and dynamic updating of DFNS based on VOF–CFD–DEM coupling. Background Technology

[0002] Fault fracture zones, characterized by fractured structures, well-developed pores and fissures, and abundant groundwater, are typical adverse geological formations that can induce water and mud inrush disasters in tunnel engineering. Driven by tunnel excavation disturbances and hydraulic head differences, the fracture zone is prone to particle erosion and migration, pore structure reconstruction, and seepage channel evolution, resulting in significant spatiotemporal unsteady characteristics in seepage morphology, pressure distribution, and flow processes. Especially when water inrushes into the tunnel excavation space, the initial air and inrushing water form a gas-liquid free surface. The coupling effect of interface advancement, phase ratio changes, and pressure release further exacerbates the complexity of the flow state.

[0003] In existing technologies, seepage theory is the fundamental method for studying groundwater movement. The linear Darcy model is suitable for low-velocity seepage dominated by viscosity. However, when pore structure expands, pressure gradients increase, or flow velocity increases, inertial effects intensify, and seepage often exhibits nonlinear characteristics, requiring the use of non-Darcy models such as Forchheimer. When water enters the tunnel excavation space and coexists with air, flow control more closely resembles the Navier-Stokes free flow model and involves the evolution of the gas-liquid interface. Computational Fluid Dynamics (CFD) can solve for fluid velocity fields, pressure fields, and phase interface propagation processes, while the Discrete Element Method (DEM) can characterize the motion, separation, and structural reconstruction of granular systems. Establishing a two-way coupling between CFD and DEM can realize the driving effect of fluid on particles and the inverse influence of particle structure changes on seepage parameters, which has significant application value in the study of nonlinear seepage and water / mud inrush processes in fault fracture zones.

[0004] However, existing CFD-DEM coupling methods for analyzing the evolution of water and mud inrush in fault fracture zones and tunnel systems still have the following shortcomings: First, many methods employ a single seepage model across the entire domain, or use fixed partitions / fixed equation settings. This makes it difficult to identify, partition, and switch the Darcy, Forchheimer, and Navier-Stokes governing equations in real time based on local flow states and medium evolution. This leads to a mismatch between the spatial and temporal configuration of porous resistance terms, inertial terms, and free flow equations, resulting in calculations of flow rate, pressure gradient, and fluid-structure interaction forces that deviate from the actual evolution process. Second, some methods still approximate the water inrush process with single-phase water flow, failing to integrate the gas-liquid two-phase free surface model and the fluid-structure coupling process into the same computational link. This easily leads to problems with the consistency of pressure release, boundary conditions, and phase ratio evolution within the tunnel, making it difficult to accurately describe the interface propagation of water after entering the tunnel air domain and the resulting reconstruction of the pressure and velocity fields. Third, multi-software coupled computation involves grid / cell correspondence, cross-platform data transmission, parallel consistency and solution stability control. Existing technologies often lack an engineering update link that can read external parameters and write them into the computational cell in real time during the iteration / coupling process. This leads to mismatch, asynchrony or numerical instability in parameter write-back, which restricts the practical application of multi-flow dynamic switching and parameter closed-loop update.

[0005] Therefore, it is necessary to establish a method that can simultaneously introduce a gas-liquid two-phase free surface model into CFD–DEM bidirectional coupled calculations, identify and partition local unit flow states based on quantifiable criteria, form a dynamic switching rule for partition-enabling item mapping, realize closed-loop updates of seepage parameters driven by pore structure evolution caused by particle migration, and have an implementable real-time update mechanism to more realistically characterize the seepage field, particle structure, and tunnel water and mud inrush disaster evolution process in fault fracture zones. Summary of the Invention

[0006] To address the limitations of traditional single-phase seepage or single governing equation models in simultaneously describing gas-liquid interaction, free surface evolution, and the coexistence of seepage, non-Darcy, and free flow within the tunnel space during the water inrush stage, this invention proposes a computational method for DFNS multi-flow regime partition identification and dynamic updating based on VOF–CFD–DEM coupling. This method aims to solve the challenges of refined numerical simulation arising from the coexistence of gas-liquid-solid multiphase flows, rapid evolution of free interfaces, and dramatic spatiotemporal transformations of seepage flow regimes during the evolution of water and mud inrush disasters in fault fracture zones.

[0007] The technical concept of this invention is as follows: a two-phase flow free surface model is used to uniformly describe the gas-liquid coexistence and interface propagation process in tunnels; fluid-solid bidirectional coupling is used to drive particle migration by fluid action and reverse update of seepage parameters by particle structure evolution; and a non-Darcy effect criterion is introduced to identify the flow state of local units, dividing the computational domain into Darcy seepage zone, non-Darcy seepage zone and free flow zone. Then, under the same computational framework, the porous resistance term and the governing equation term are dynamically configured to realize the continuous switching and transition of multiple flow states in space and time, thereby realizing the continuous simulation of the entire process of water inrush and mud inrush disasters and the adaptive update of the parameter field.

[0008] To achieve the above-mentioned technical objectives, the present invention adopts the following technical solution: This invention discloses a computational method for multi-flow regime partition identification and dynamic updating based on VOF–CFD–DEM coupling in DFNS, comprising the following steps: Construct a fluid-solid-gas multi-field coupled computational domain and a discretized model; A unified governing equation system for DFNS is established based on the volume averaging method. The medium properties in different regions of the computational domain are characterized by parameterization of porosity, porous medium resistance source term and momentum source term. Based on the medium characteristics of different regions within the computational domain, flow regime partition labels are generated using an established multi-flow regime criterion system. Based on the flow regime partition labels, parameter mapping logic is formulated to enable adaptive switching of the DFNS unified control equation system; The DFNS unified control equation system is discretized and solved using real physical velocities. The discretized equation system is iteratively updated at each time step using a pressure-velocity coupling algorithm until the preset convergence condition is met, thereby obtaining the basic field variable input parameters for fluid-particle interaction force calculation. The drag force of the fluid is calculated using a drag force model that considers local porosity correction, and the translational and rotational states of the particles are updated. The sum of particle volumes in each updated fluid unit is statistically analyzed to update the unit porosity. At the same time, the relationship between porosity and the average permeability of the fluid unit is established. Porosity and permeability are reconstructed and updated at the end of each coupling step until the current simulation time or coupling step number reaches the preset termination condition, and the optimal parameter combination is output.

[0009] A further preferred approach involves constructing a fluid-solid-gas multi-field coupled computational domain and a discretized model, including the following steps: Establish a global computational domain that includes the surrounding rock matrix, fault fracture zone, and tunnel excavation space; Within the global computational domain, the fluid computational domain is discretized using the finite volume method, and a discrete element model is constructed for the solid particle system in the fault fracture zone, forming a coupled computational model of the fluid control volume and the particle assembly. The VOF multiphase flow model is introduced into the fluid computation domain to achieve a unified flow field simulation of the gas phase in the tunnel excavation space, as well as the liquid phase in the surrounding rock and fracture zone.

[0010] A further preferred approach is to introduce a VOF multiphase flow model into the fluid computational domain, including the following steps: A discretized model is established based on the local volume averaging method, and porosity is defined in the discretized model. The phase volume fraction from the VOF multiphase flow model is introduced into the discretized model to weight the calculation of the mixture properties parameters of the gas-liquid two-phase mixture; For the gas-liquid interface during the water inrush process, a continuous surface force model is used to calculate the surface tension.

[0011] As a preferred embodiment of the above technical solution, in the DFNS unified governing equation system, based on the difference in values ​​between the porous medium resistance source term and the volume force source term, the unified momentum conservation equation corresponds to different flow characteristics in different regions of the computational domain, and its expression is as follows: ; In the formula: For fluid effective occupancy rate, The density of the mixed fluid within the grid cell. For numerical computation time, For hybrid viscosity, u is the fluid physical velocity vector, p is the fluid pressure, g is the gravitational acceleration vector, and S is the porous medium resistance source term; For vector differential operators, For fluid pressure gradient, Let be the fluid velocity gradient tensor. For volumetric force source terms; Among them, the fluid motion in the surrounding rock matrix zone is driven by the pressure gradient, dominated by viscous resistance, and the unified momentum conservation equation degenerates into the classical Darcy's law. The pore structure in the fault fracture zone is heterogeneous, and the increase in flow velocity leads to a significant local inertial effect. The pressure drop and flow velocity have a nonlinear relationship, and the unified momentum conservation equation degenerates into the Forchheimer nonlinear seepage equation. The actual flow velocity of the fluid in the free flow zone of the tunnel excavation space increases significantly between the particle skeleton due to the reduction of the flow cross section. The unified momentum conservation equation is solved by using the Navier-Stokes equation with porosity term and combined with the standard k-epsilon turbulence model.

[0012] Based on the above technical solution, the further steps to establish a multi-flow criterion system include: The Fochheimer number Fo is selected to characterize the nonlinear intensity of fluid motion in porous media; Based on the Fochheimer number Fo, a non-Darcy effect index E is introduced to normalize the nonlinearity of the flow to the [0,1] interval, quantitatively defining the flow regime transition boundary; the discrimination criteria of the multi-flow regime criterion system are as follows: when When the inertial term is negligible, it is determined to be Darcy laminar flow; when At that time, the corresponding critical Fochheimer number At this point, the inertial resistance is significant, and it is determined to be Forchheimer nonlinear seepage; when At this point, the influence of the viscous term weakens, and the flow is completely dominated by inertia, which is determined to be Navier-Stokes free turbulence.

[0013] As a further preferred technical solution, in the multi-flow criterion system, a cross-scale coupling architecture and mesh mapping of fluid-particle are constructed, and the false gas phase drag in the fluid-solid-gas multi-field coupling is corrected based on the velocity field correction algorithm of phase threshold. When the water phase volume fraction of the mesh cell in the fluid computation domain is determined to be the gas phase dominant cell, the velocity is zeroed and corrected.

[0014] A further preferred technical solution is that the parameter mapping logic based on the flow regime partitioning label is as follows: When the flow is determined to be Darcy laminar flow, a flow regime label "Darcy" is generated, linear drag is activated, nonlinear drag is deactivated, and the porous media drag source term degenerates into... ; When the flow is determined to be non-Darcy flow, a flow regime label (Forchheimer) is generated, and both linear drag and inertial drag are activated. The drag source term in porous media degenerates into... ; When the flow is determined to be free turbulence, a flow regime label Navier-Stokes is generated, the momentum source term is set to zero, and the equations are automatically reduced to Navier-Stokes equations.

[0015] As a further preferred embodiment of the above technical solution, the steps for obtaining the input parameters of the basic field variables for fluid-particle force calculation include: Construct a multiphase flow control equation that includes porosity parameters based on the DFNS unified control equation system; The multiphase flow control equations, which include porosity parameters, are iteratively updated at each time step using a pressure-velocity coupling algorithm. After convergence, the mixed pressure field, velocity field, pressure gradient field, and gas-liquid phase field information of the entire field grid are obtained, which serve as the input parameters of the basic field variables for calculating fluid-particle interaction forces.

[0016] A further technical solution is to consider the following expression for the drag force model, which takes into account local porosity correction: ; ; In the formula: The forces acting on the particles in the gas-liquid-solid three-phase coupled field. For fluid drag force, For the drag force of a single particle, The average porosity of the fluid unit containing the particle. This is an empirical factor for correcting the local average porosity.

[0017] A further technical solution involves the following steps for outputting the optimal parameter combination: A particle-scale similarity and parameter equivalence correction mechanism is introduced. At the end of each coupling step, the porosity and permeability are reconstructed and updated. The updated medium parameters are fed back to the fluid solution end as input parameters for the next coupling step momentum source term construction and flow regime identification calculation. Determine whether the current simulation time or coupling step count has reached the preset termination condition. If it has, output the optimal parameter combination. If it has not, enter the next coupling loop until the preset termination condition is reached, and then output the optimal parameter combination.

[0018] In addition, the present invention also discloses an electronic device, the electronic device including at least one processor and a memory communicatively connected to the processor; wherein the memory stores instructions executable by the processor, the instructions being executed by the processor to enable the processor to perform the steps of the above-described calculation method.

[0019] Compared with the prior art, the beneficial effects of this embodiment are as follows: To address the rapid transformation and highly unsteady evolution of parameters in the "porous medium seepage—non-Darcy seepage—free flow in tunnel excavation space" characteristics of water and mud inrush disasters in fault fracture zones, this invention introduces a gas-liquid two-phase free surface description within a CFD–DEM bidirectional coupling framework. This enables unified calculation of interface propagation and phase ratio evolution during water inrush into the tunnel air domain. Simultaneously, based on quantifiable non-Darcy effect criteria, local flow regimes in CFD grid cells are identified in real time, establishing a "partition-enabling term mapping" rule. This allows for dynamic configuration and continuous switching of the Darcy, Forchheimer, and Navier–Stokes governing equations, as well as porous resistance / source terms, avoiding equation configuration mismatches caused by a single flow regime or fixed partitions across the entire domain. Furthermore, the porosity, permeability, and non-Darcy parameters are updated in a closed loop using DEM particle migration and pore structure evolution results, and fed back to the CFD solution in real time. This ensures consistent evolution of the fluid equations and particle structure during coupled iteration, and, combined with an implementable real-time parameter write-back mechanism, improves data consistency and computational stability under parallel conditions. The method of this invention can more realistically depict key processes such as flow regime evolution, reconstruction of dominant permeable channels, and water-air interface interaction under water inrush loading, thereby improving the physical consistency, numerical stability, and engineering applicability of the simulation of the entire process of water inrush and mud inrush disasters in fault fracture zones. Attached Figure Description

[0020] To more clearly illustrate the technical solutions of the embodiments of this application, the drawings used in the description of the embodiments of this application will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0021] Figure 1 This is a schematic diagram of the geometric model and mesh generation of the CFD computational domain; Figure 2 The diagram shows the CFD flow field model, the DEM particle model, and the region grouping. (a) Schematic diagram of fluid computation domain meshing, (b) Global computation domain, including fluid computation domain and particle computation domain, (c) Particle computation domain model. Figure 3 This is a flowchart of the bidirectional coupling calculation of VOF–CFD–DEM for DFNS multi-flow state identification and dynamic updating in this invention; Figure 4 The initial phase distribution and the volume fraction cloud map of the gas-liquid two-phase (VOF) phase during the tunnel water inrush process; Figure 5 The spatial distribution of the non-Darcy effect index E; Figure 6 This is a typical velocity field distribution diagram; Figure 7This is a typical pressure field distribution diagram; Figure 8 This is a map showing the migration and displacement distribution of particles in the fault fracture zone. Figure 9 A schematic diagram showing the layout of monitoring points in the computing domain; Figure 10 Velocity time history curves for representative monitoring points in different flow regime zones; Figure 11 Pressure time history curves for representative monitoring points in different flow regime zones; Figure 12 The velocity-time history curve of the monitoring section at the working face; Figure 13 The pressure time history curve of the monitoring section at the working face. Detailed Implementation

[0022] Embodiments of the present invention are described in detail below, examples of which are illustrated in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and intended to explain the present invention, and should not be construed as limiting the present invention.

[0023] Example 1: Reference Figures 1-13 This invention discloses a computational method for multi-flow regime partition identification and dynamic updating of DFNS based on VOF–CFD–DEM coupling, comprising the following steps: S1: Constructing a fluid-solid-gas multi-field coupled computational domain and discretization model; S2: Establish a unified governing equation system for DFNS based on the momentum conservation equation; S3: Establish a multi-flow criterion system based on the Fochheimer number Fo and the non-Darcy effect index E, and use the multi-flow criterion system to determine the flow morphology partition of the computational domain and generate flow morphology partition labels. S4: Based on the DFNS unified control equation system, the porous medium momentum source term is used to characterize the medium resistance characteristics of different flow regions in the computational domain, and the parameter mapping logic is formulated according to the real-time results of the flow region determination, so that the DFNS unified control equation system can be adaptively switched. S5: The discrete equations are iteratively updated at each time step using a pressure-velocity coupling algorithm until the preset convergence condition is met, thus obtaining the basic field variable input parameters for fluid-particle interaction force calculation. S6: The drag force of the fluid is calculated using a drag force model that considers local porosity correction, and the translational and rotational states of the particles are updated based on the particle dynamics control equation. S7: Based on the updated particle positions, calculate the sum of particle volumes in each fluid unit, update the unit porosity, establish the relationship between porosity and average permeability of the fluid unit based on the KC equation, and reconstruct and update porosity and permeability at the end of each coupling step until the current simulation time or coupling step number reaches the preset termination condition.

[0024] Specifically, in step S1, the evolution of tunnel water and mud inrush disasters is essentially a dynamic coupling process of "continuous flow state transformation—particle migration—tunnel protection structure response." Given the significant spatial distribution differences in water flow states depending on the medium properties, a computational geometric domain is first established, encompassing the surrounding rock matrix, fault fracture zone, and tunnel excavation space. Based on this, the fluid computational domain is discretized using the finite volume method (FVM), and the solid particle system of the fault fracture zone is discretized using the discrete element method (DEM), forming a coupled computational object of fluid and particle aggregation. Gravity and volume force terms are set in the fluid control equations, including vertically downward gravity, with a gravitational acceleration of 9.8 m / s². 2 The fluid-structure interaction momentum exchange source term between the fluid and DEM particles is used to initialize the groundwater level and boundary pressure conditions based on geological survey data.

[0025] To address the phenomena of "groundwater inrush" and "air disturbance within the tunnel" involved in water and mud inrush processes, and to achieve a unified description of the coexistence of gas and liquid phases and the evolution of the free liquid surface, a VOF (Volume of Fluid) multiphase flow model is introduced into the fluid computational domain. This model assumes the existence of a gas phase (air) and a liquid phase (water) within the cell mesh, and defines... and Let be the volume fractions of the gas phase and the liquid phase, respectively, and satisfy . Mixed fluid density within the grid cell viscosity with hybrid power The calculation is as follows: (1)

[0026] In the formula, , These are the densities (kg / m³) of pure water and pure air, respectively. , These are the viscosities (Pa·s) of pure water and pure air, respectively. This is a mixture property parameter. and This will serve as the basic input for solving the momentum equation, thereby automatically simulating gas-phase turbulence within the tunnel cavity and liquid-phase seepage within the surrounding rock.

[0027] To accurately describe the coexistence of fluids, particles, and gas-liquid two phases at the microscopic unit scale, this invention establishes a discretized model based on the local volume averaging method. For any computational grid cell, its total volume... Pore ​​volume divided into fluid-occupied areas Volume occupied by solid particles Define porosity The dynamic evolution relationship of the fluid-flowable space is as follows: (3) (4) In the formula, The interior is further filled with immiscible liquid and gas phases. This definition clarifies that when particles migrate into the grid, the effective flow field space will decrease accordingly, thereby triggering fluid displacement and acceleration effects.

[0028] Introducing the concept of phase volume fraction from the VOF multiphase flow model, weighted calculations are performed on the property parameters of the gas-liquid two-phase mixture, defining... The volume fraction of liquid phase within the fluid pores: (5) It can be seen that when Then the unit pores are filled with water, when The pores of the unit cell are then filled with air. Therefore, the local average density and average dynamic viscosity of the fluid unit cell are calculated by weighting the volume ratios of the two phases.

[0029] In this embodiment, the surface tension is considered in relation to the topological evolution of the gas-liquid interface during the water inrush process. The influence of surface tension on the generation of sliding induced waves is significant and cannot be ignored. This invention employs a continuous surface force (CSF) model to measure surface tension. Calculation: (6) In the formula, Let be the local curvature at the interface, and n be the unit normal vector of the interface. is the surface tension coefficient.

[0030] Based on existing physical experimental data and current research, the surface tension coefficient of the water-air interface can be taken as [value missing]. The introduction of this term helps improve the model's ability to depict local curvature changes and wave surface details of free liquid surfaces.

[0031] In step S2, to overcome the limitations of traditional segmented calculation methods, this invention establishes a unified Darcy-Forchheimer-Navier-Stokes (DFNS) governing equation system based on the volume averaging method, including mass conservation equations and momentum conservation equations. Furthermore, it characterizes the medium properties of different regions parametrically through porosity and momentum source terms, wherein: Unified continuity equation: (7) (8) Unified momentum conservation equation: (9) In the formula, For fluid effective occupancy rate, The density of the mixed fluid within the grid cell. For numerical computation time, Here, is the hybrid viscosity, u is the fluid physical velocity vector, p is the fluid pressure, g is the gravitational acceleration vector, and S is the porous medium resistance source term. For vector differential operators, For fluid pressure gradient, For the fluid velocity gradient tensor, This is a volumetric force source term.

[0032] Among them, the momentum source term for porous media consists of a viscous loss term (Darcy term) and an inertial loss term (Forchheimer term): (10) (11) In the formula, k is the permeability of the medium. The non-Darcy inertial drag coefficient (related to the non-Darcy factor β), u m U is the fluid dynamic viscosity coefficient. x u y u z These are the velocities in the x, y, and z directions, respectively.

[0033] The degradation forms of three typical flow regimes correspond to the physical scenarios of tunnel water and mud inrush disasters, based on porosity. With the momentum source term S of porous media i Due to the differences in the values ​​of the above unified governing equations, the system of flow equations automatically degenerates into the following three types of flow equations in different regions: (a) Darcy laminar flow region; During the water inrush incubation stage, the aquifer tectonic rock system is in its original state or only slightly disturbed by excavation. The hydraulic connection between the surrounding rock and geological defects is relatively stable. The surrounding rock has small pore size and limited connectivity. Fluid movement is mainly driven by pressure gradients, and flow is mainly controlled by viscous resistance. Pressure drop and flow velocity have a linear relationship, which conforms to the Darcy flow characteristics of porous media. In this stage, the seepage calculation adopts Darcy's law: The surrounding rock has small pore size and limited connectivity. Fluid movement is driven by pressure gradients and dominated by viscous drag. The Reynolds number is extremely low, and the inertial drag coefficient is very low. Retaining the viscous term, unifying the momentum conservation equation, neglecting the inertial and convection terms, degenerates into the classical Darcy's law: (12) In the formula, V s The apparent velocity of the fluid.

[0034] (b) Non-Darcy (Forchheimer) high-velocity seepage zone; Fault fracture zones are composed of a mixture of breccia, fragments, and fine particles. Their internal pore structure continuously evolves with particle migration and local scouring, causing changes in permeability, flow resistance coefficient, and local flow regime. Local velocity increases, inertial forces and localized flow around the fault become significant, resulting in a nonlinear flow regime. Traditional CFD-DEM models typically use a fixed permeability, which is insufficient to reflect this dynamic evolution of the flow regime. The Forchheimer equation, by adding an inertial term to the Darcy linear drag, can reflect the nonlinear energy loss caused by pore irregularities and increased flow velocity.

[0035] The fault fracture zone exhibits a heterogeneous pore structure that evolves continuously with particle migration. Increased flow velocity leads to significant local inertial effects, resulting in a nonlinear relationship between pressure drop and flow velocity. Therefore, viscous drag and inertial drag terms are retained in this region, and the Forchheimer nonlinear seepage equation is as follows: (13) In the formula, β is the non-Darcy factor, which is closely related to porosity and particle size distribution and is dynamically updated with the fluid-structure interaction process. The non-Darcy factor is calculated as follows: (14) (c) Free flow zone in tunnel excavation space (Navier-Stokes); Initially, there is no rock mass skeleton obstructing the tunnel cavity, and the flow field is controlled only by the fluid's own viscosity and the gas-liquid interfacial tension. With the occurrence of a water inrush disaster, high-pressure water carries a large number of fractured solid particles into the tunnel, forming a high-concentration liquid-solid two-phase turbulent flow. At this point, the porosity within the tunnel region grid is no longer constant at 1, but exhibits a spatiotemporal dynamic decay as the particle group migrates in. The actual flow velocity of the fluid between the particle skeleton increases significantly due to the reduction in the flow cross-section. Initially, the porosity is close to 1, and the permeability tends to infinity, i.e., the porous medium resistance source term S=0. At each coupled time step, the grid porosity is updated based on the DEM particle position mapping. This part is solved using the Navier-Stokes equations that include the porosity term: (15) (16) In the formula, This represents the gas-liquid surface tension source term calculated based on the VOF multiphase flow model, using the continuous surface force (CSF) model. This indicates the drag force of the particles.

[0036] Furthermore, at the boundaries of the aforementioned different flow regimes, such as the fault fracture zone-tunnel interface, strict pressure and velocity continuity conditions are defined to ensure the physical consistency of the partitioned coupled calculations: (17) (18) In the formula, The pressure at the Darcy domain boundary. For the pressure at the boundary of the Forchheimer domain, The velocity at the boundary of the Darcy domain. The velocity at the boundary of the Forchheimer domain, The pressure at the boundary of the Navier-Stokes domain. The velocity at the boundary of the Navier-Stokes domain.

[0037] Through the construction of the above unified framework, seamless coupling simulation of the entire domain, from low-speed seepage to high-speed free turbulence and from single-phase flow to gas-liquid two-phase flow, was achieved.

[0038] In step S3, to achieve accurate identification of the flow regime (Darcy / Forchheimer / Navier-Stokes) of each unit throughout the entire water inrush process, this invention selects the Forchheimer number F. oUsing Forchheimernumber and non-Darcy effect index E as core parameters, a multi-flow criterion system is established. Based on the flow velocity, pressure gradient and medium pore structure characteristics of the local grid, the degree of flow deviation from linear Darcy law is quantified.

[0039] Specifically, the Fochheimer number F o As a dimensionless number characterizing the proportion of inertial drag to total drag, it intuitively reflects the nonlinear intensity of fluid motion in porous media, and its expression is as follows: (19) In the formula, β is the non-Darcy factor. For the density of the mixed fluid, v s For flow velocity vectors, Let be the dynamic viscosity of the mixed fluid, and k be the permeability.

[0040] To overcome the limitations of a single index, a non-Darcy effect index E is introduced to normalize the nonlinearity of the flow to the [0,1] interval, and is used to quantitatively delineate the flow transition boundary. The expression for the non-Darcy effect index E is as follows: (20) The flow regime of each unit is determined and classified using the non-Darcy effect index E, and the discrimination criterion is as follows: when At this time, the inertial term can be ignored, and it is determined to be Darcy laminar flow; when At that time, the corresponding critical Fochheimer number At that time, the inertial resistance was significant, and it was determined to be Forchheimer nonlinear seepage; when At this point, the influence of the viscous term weakens, and the flow is completely dominated by inertia, which is determined to be Navier-Stokes free turbulence.

[0041] To achieve efficient bidirectional data interaction between the CFD flow field and the DEM particle field, this invention constructs a C / S (Client-Server) coupled architecture based on the TCP / IP protocol. Utilizing a Python open-source library, the DEM software is configured as the p2plinkServer master server, and the CFD software as the p2plinkClient client. Both systems transmit floating-point numbers, strings, and NumPy arrays at high frequency in real time via local or remote TCP sockets. In this embodiment, the Python open-source libraries used are itasca and pyfluent, the DEM software is PFC3D, and the CFD software is Fluent.

[0042] During the coupling initialization phase, the CFD fluid computational domain mesh information is read in, and a background fluid mesh consistent with the CFD fluid computational domain space is constructed using Python functions. Finally, a mapping index relationship between the fluid elements and the DEM background fluid elements is established.

[0043] To address the "false gas phase drag" problem in gas-liquid-solid three-phase flow coupling calculations, where airborne particles incorrectly calculate the fluid forces based on water medium parameters, this invention uses a velocity field correction algorithm based on phase thresholds to overcome this problem.

[0044] After convergence at each CFD time step, all mesh elements, the fluid velocity vector u, and the water volume fraction within the fluid computational domain are extracted via a user-defined function (UDF) or a data post-processing script. Set a critical threshold for the volume fraction of the aqueous phase. Iterate through all fluid mesh elements, based on the current time. Make a judgment: like The unit is determined to be an effective liquid phase unit, and its original flow field velocity information is retained; like This unit was determined to be a gas-phase dominant unit.

[0045] For regions identified as "gas-dominant units," the coupled velocity field is corrected and transferred. A velocity zeroing correction operation is performed before the data is transmitted to the PFC terminal. In this embodiment, the preferred method is... .

[0046] (twenty one) In the formula, The fluid velocity that is actually written to the swap file and passed to the PFC program.

[0047] This correction step sets the velocity of the gas-dominant region in the flow field data read by the PFC to zero. This means that in the DEM calculation, particles in the air are only controlled by gravity and interparticle contact forces, eliminating the non-dominant aerodynamic drag generated by low-density air. Given that the impact force of water flow in a sudden flood is much greater than the driving force of airflow, this simplification effectively avoids non-physical gas-phase interference, ensuring that the fluid-structure interaction calculation focuses only on the strong interaction between "water flow and particles," significantly improving the physical realism and computational stability of the simulation.

[0048] In step S4, based on the aforementioned unified governing equation system, a generalized porous media resistance source term is adopted within the computational domain. This is used to characterize the resistance features of media in different flow regimes. In the Fluent solver, the source term for porous media resistance... The sink term, appended to the right-hand side of the momentum equation, has its component expression along the i-th direction constructed as follows: (twenty two) In the formula, S is the apparent velocity vector containing three-dimensional directional components; is the equivalent viscous drag coefficient; is the inertial drag coefficient; and is the velocity vector magnitude. The first term in the formula represents the Darcy viscous drag proportional to the velocity, and the second term represents the Forchheimer inertial drag related to the square of the velocity.

[0049] Based on the real-time flow regime discrimination results output above, the adaptive switching of the unified control equation system is achieved through the following parameter mapping logic: When the flow is determined to be Darcy laminar flow, a flow regime label "Darcy" is generated, linear drag is activated, and nonlinear drag is deactivated. Let... ,make The source term degenerates into .

[0050] When the flow is determined to be non-Darcy flow, a flow regime label (Forchheimer) is generated, and both linear drag and inertial drag are activated. Let... ,make The source term degenerates into .

[0051] When the flow is determined to be free turbulence, a flow regime label, Navier-Stokes, is generated, and the resistance source terms for porous media are all zeroed out. Let... , And set the mesh porosity to ,at this time The discrete equations are automatically reduced to Navier-Stokes equations.

[0052] The parameters described above are not statically set, but rather dynamically loaded at the time step level using a Python-Fluent interactive script: at the beginning of each time step, the Python script reads the current flow field velocity and the non-Darcy effect index E through a TCP interface; based on the parameter mapping logic, it generates the corresponding full-field mesh in memory. and Coefficient matrix; using Fluent's UDF macro DEFINE_PROFILE, the calculated coefficient matrix is ​​directly written to the solver memory address UDM, overwriting the medium properties of the previous time step.

[0053] In step S5, based on the flow regime partition labels and their corresponding adaptive momentum source terms formed above, a multiphase flow control equation including porosity parameters is constructed within the framework of a unified control equation system. This multiphase flow control equation including porosity parameters can adaptively solve the flow regime characteristics of the corresponding region by assigning different values ​​of porosity and momentum source terms to different physical regions of the computational domain. In the solution system, in order to eliminate the flow velocity distortion caused by particle size magnification or equivalent medium assumptions, apparent velocity and physical velocity are clearly distinguished, and a momentum equation based on physical velocity is used for discrete solution.

[0054] Define u as the true physical velocity, that is, the actual flow velocity of the fluid within the pore, which is different from the apparent velocity. The conversion relationship is as follows: (twenty three) By explicitly introducing porosity into the transient and convection terms of the governing equations The solver directly calculates u, which is the actual flow velocity within the pores. This process ensures that the inertial force term and momentum convection effect can be correctly captured in non-Darcy flow conditions, avoiding the underestimation of flow velocity and distortion of the physical field caused by the traditional Darcy model's calculation of only the apparent velocity.

[0055] At each time step, the discrete equations are iteratively updated using a pressure-velocity coupling algorithm until a preset convergence criterion is met. After convergence, the hybrid pressure field p, velocity field u, and pressure gradient field of the entire field mesh are obtained. And simultaneously obtain gas-liquid phase field information. This provides the basic field variable inputs for subsequent fluid-particle interaction force calculations.

[0056] After the CFD solver completes the calculation and converges at the current time step, this invention utilizes an embedded Python script interface to traverse and extract the center velocity vector, pressure gradient, and water volume fraction of all grid cells within the fluid computational domain. Subsequently, based on a pre-constructed spatial mapping index of "CFD fluid cell - DEM background fluid grid," the extracted field variable data is rearranged and serialized using Python's NumPy library to generate a standard field variable array that strictly corresponds one-to-one with the background grid ID of the discrete element computation end. Finally, the serialized flow field data packet is sent to the DEM particle calculation module in real time via inter-process communication (IPC). In this embodiment, TCP Socket communication is preferred, but shared memory or file exchange methods can also be used to ensure efficient and low-latency transmission of massive coupled data between heterogeneous software.

[0057] In step S6, the particles are subjected to forces in the gas-liquid-solid three-phase coupled field. It is usually composed of fluid drag force, fluid pressure gradient force and buoyancy. Among them, drag force is the key force term characterizing fluid-driven particle migration and channel evolution.

[0058] To address the variable porosity characteristics of fault fracture zones, this invention employs a drag force model that considers local porosity correction to assess fluid drag force. Calculations are performed. As a non-limiting preferred embodiment, the expression for the drag force model considering local porosity correction is as follows: (twenty four) Fluid drag force Defined as the resistance exerted by a fluid on a particle, it always acts on the particle's center of mass; the particle does not exert a rotational torque; the fluid drag force. Calculate using the following formula: (25) In the formula: For the drag force of a single particle, The average porosity of the fluid unit containing the particle. This is an empirical factor for correcting the local average porosity. This correction term will be applicable to both high and low porosity systems as well as a wide range of Reynolds numbers.

[0059] (26) In the formula: The drag coefficient, Let r be the fluid density, r be the particle radius, and u be the fluid velocity. Particle velocity. Drag coefficient. Calculate using the following formula: (27) In the formula: is the particle Reynolds number.

[0060] In addition, adjust empirical factors Defined as: (28) Reynolds number Calculate using the following formula: (29) In the formula: is the fluid dynamic viscosity coefficient.

[0061] In step S7, the DEM software acquires the flow velocity, pressure, and pressure gradient of the fluid element via TCP socket transmission. Based on Newton's equations of motion, the translational and rotational states of the particles are updated. The particle dynamics control equations are the translational equation and the rotational equation, respectively, and their calculation expressions are as follows: (30) (31) In the formula: Indicates the velocity of particle movement; Indicates the time required for numerical computation; It represents the particle contact force between particles, including particle-particle contact force and particle-wall contact force, and reflects the interaction between particles through a contact model; This represents non-contact forces between particles, mainly including van der Waals forces, electromagnetic forces, and capillary forces. Indicates an externally applied force; This represents the total fluid-particle interaction force exerted by the fluid on the particles; Indicates particle mass; Represents gravitational acceleration; Indicates the angular velocity of the particle; This represents the contact torque experienced by the particle; This represents the moment of inertia of a particle.

[0062] After completing the dynamic solution of the coupling step, the discrete element computing terminal outputs updated information such as particle position, velocity and contact network state, which is used for porosity reconstruction and permeability update, and then proceeds to the next coupling time step until the preset total coupling time or coupling step number termination condition is met.

[0063] After the particle positions are updated, the solid volume fraction within the fluid grid cell changes, causing porosity and permeability to evolve over time. Therefore, at the end of each coupling step, this invention calculates the sum of particle volumes within each fluid cell based on the particle positions output from the discrete element, according to a pre-established particle-cell mapping relationship, and updates the cell porosity. The calculation formula is as follows: (32) In the formula: Porosity of the fluid unit; This is the sum of the volumes of all particles within the measurement area; To measure the volume of the fluid grid within the measurement area.

[0064] In the CFD-DEM coupled analysis system, the Kozeny-Carman (KC) equation is used to establish the relationship between porosity and the average permeability of fluid elements, characterizing the permeability changes caused by particle motion and realizing the influence of DEM particle motion on the seepage field distribution. (33) In the formula: The permeability of the fluid unit; denoted as ρ, where ρ is the porosity of the fluid element; r is the radius of the particle.

[0065] Permeability coefficient The conversion formula with permeability k is: (34) To avoid non-physical permeability under extreme porosity conditions, this invention sets the initial porosity to 1 for free-flow regions such as tunnel excavation spaces, to ensure that the model is consistent with the physical scenario.

[0066] To address the issue of "mismatch between model permeability and engineering prototype" caused by scaling up particle size due to computational efficiency limitations in discrete element numerical simulations, this invention introduces a particle size similarity and parameter equivalence correction mechanism. At the end of each coupling step, porosity and permeability are reconstructed and updated, and the updated medium parameters are fed back to the fluid solution end, realizing a closed-loop coupling of "particle migration - pore structure change - permeability parameter update - flow field redistribution".

[0067] In the discretized model, the model particle radius r is defined as the same as the actual particle radius of the engineering prototype. Scale ratio between for: (35) Since enlarging the particle geometry alters the characteristic scale of the pore structure, the geometric permeability k naturally formed by the enlarged particle geometry packing in the model differs from the original permeability. There is a theoretical conversion relationship: (36) That is, the geometric permeability of the model is theoretically magnified. The particle drag force is also affected by the amplification of the particle radius, according to the conversion relationship.

[0068] (37) (38) To ensure that the seepage resistance in the fluid solver matches the preset Reynolds number similarity or Darcy flow similarity criterion, this invention uses script control to perform consistency correction on the viscous resistance coefficient fed back to the fluid end of the DEM. The corrected equivalent input parameters... The resistance source term S of porous media must satisfy: (39) (40) in, To simulate fluid viscosity, This is the permeability correction factor determined based on the fluid viscosity scaling strategy. When fluid viscosity is only amplified... To maintain Reynolds number similarity, we need to take... .

[0069] It should be noted that porosity It is a dimensionless quantity and generally does not change with scale ratio. Transformation; while the permeability k and its corresponding viscous resistance coefficient 1 / k are directly affected by the particle characteristic scale. This invention calculates the above formulas (36)-(38) through an embedded script, realizes the automatic correction of permeability and resistance parameters, and provides a consistent parameter basis for subsequent flow regime partitioning and source term invocation.

[0070] After completing the porosity and permeability update, the updated values ​​will be... , The viscous drag coefficient, non-Darcy coefficient, and inertial drag parameters derived therefrom are sorted and encapsulated according to the global element number to form a parameter array that corresponds one-to-one with the fluid mesh. This array is then fed back to the fluid solver as input for the construction of the next coupled step momentum source term and the flow regime identification calculation.

[0071] Then it is determined whether the current simulation time or coupling step count has reached the set coupling physical time; if not, the next coupling loop is entered.

[0072] Example 2: This embodiment takes a water and mud inrush disaster occurring in a tunnel traversing a fault fracture zone as its background. It addresses the technical problems encountered during the disaster's evolution, such as the rapid transformation of seepage patterns between porous media seepage and free flow in the tunnel excavation space, and the continuous evolution of pore structure and permeability due to particle migration, making it difficult to fix and set flow field parameters. It combines... Figures 1-13 The technical solution of the present invention will be further described below.

[0073] The overall calculation process is as follows: Figure 2 As shown, the specific implementation steps are as follows: (1) Engineering background and geometric model construction; This embodiment selects a tunnel traversing a fault fracture zone as the engineering application object. The average burial depth of this tunnel section is approximately 179m, the groundwater level is approximately 50m below the surface, and the tunnel cross-section is approximately 13m high and 15.5m in diameter. During construction, a fault fracture zone was exposed, the surrounding rock integrity was poor, and groundwater was abundant, presenting typical engineering conditions for water and mud inrush disasters, making it suitable for verifying the technical solution of this invention.

[0074] To reduce the impact of boundary conditions on tunnel seepage and stress response, based on the empirical criterion for the influence range of underground cavern boundaries, when the boundary of the computational domain is more than three times the tunnel diameter from the center, its influence on the surrounding rock stress and seepage field can be ignored. Therefore, in the horizontal direction, the model boundary is set at 60m from the center of the tunnel section; in the vertical direction, the lower boundary is set 48m below the tunnel bottom, and the upper boundary is set to the ground surface; in the tunnel axial direction, the computational domain extends 70m before and after the fault. Considering the symmetrical distribution of the engineering structure along the tunnel axis, a symmetrical boundary is set on the plane of the tunnel axis, and only half of the overall engineering model is used for calculation, ultimately establishing... Figure 1 The three-dimensional computational model shown has an overall size of 60m×190m×190m.

[0075] To standardize the geometric description, the following coordinate system is established: the tunnel axis is the y-axis, the vertical axis is the z-axis, and the horizontal axis is the x-axis; the tunnel centerline is located along the y-axis, and a fault fracture zone is exposed at y=92m during tunnel excavation; the groundwater level is located approximately 50m below the surface (corresponding elevation is as follows). Figure 2 (Given). In this embodiment, the tunnel geometry and fault zone spatial locations in the CFD computational domain and the DEM particle domain are kept consistent to ensure the consistency of the coupling object between the fluid control volume and the particle system. The model as a whole is considered as a fixed boundary, and geometric updates caused by the overall displacement of the surrounding rock are not considered.

[0076] The finite volume method was used for mesh discretization of the fluid computational domain. To improve the resolution at the fault-tunnel connection and the free liquid surface propagation process within the tunnel, the mesh was refined in the tunnel region with a mesh size of approximately 3m; the mesh in the fault fracture zone region was refined with a mesh size of approximately 7m; and the mesh in the surrounding rock region was approximately 12m, which could be adjusted appropriately within this order of magnitude according to computational resources. The total number of meshes was 3312. The mesh generation and regional refinement distribution are as follows: Figure 1 As shown.

[0077] Meanwhile, to balance computational efficiency and local response accuracy, the DEM particle system uses relatively smaller particle sizes around tunnels and faults, while appropriately increasing particle sizes in areas far from tunnels and faults.

[0078] (2) Constitutive parameters and boundary conditions settings; To ensure the reproducibility of the calculation process, this embodiment provides the fluid computational domain model and parameter settings, the constitutive and initial equilibrium settings of the DEM particle system, as well as the boundary conditions, initial conditions, and solution control parameters.

[0079] 1) Fluid properties and multiphase flow model settings: The fluid computational domain employs a water-air two-phase system. To describe the free surface advancement and gas-liquid coexistence state within the tunnel excavation space during water inrush, the VOF (Volume of Fluid) multiphase flow model is activated in Fluent 2024 R2, defining the fluid computational domain as an immiscible gas-liquid two-phase system. The VOF multiphase flow model defaults to 2 phases, with air as the primary phase, used to fill the cavity region after tunnel excavation, and liquid water as the secondary phase, serving as the source fluid for fault-induced water inrush. Physical property parameters, including density and dynamic viscosity, are assigned to both phases to ensure that the mixing density and viscosity vary with the phase volume fraction and are calculated based on actual material parameters.

[0080] To reflect the hydrostatic pressure gradient and head-driven effect of groundwater, a gravity term is enabled in the solver, and the gravitational acceleration is taken as... The direction is set to the −z direction of the coordinate system. To be consistent with the hydrostatic pressure reference, this embodiment sets the operating pressure reference or operating density of the solver accordingly. In this embodiment, the pressure reference is set to atmospheric pressure (101325 Pa), and the operating density is set to 0. The hydrostatic pressure distribution is calculated based on the elevation and gravity term of the VOF multiphase flow model in the fluid computational domain to ensure that the formation of the hydrostatic pressure gradient is consistent with the definition of subsequent boundary pressures. These settings collectively form the basis for defining the multiphase flow system and its physical properties, providing a unified physical framework for subsequent calculations of the free surface evolution during water intrusion processes.

[0081] Table 1. Fluid phase properties (Fluent 2024R2)

[0082] 2) Porous media model and regional parameter assignment: The surrounding rock matrix zone, fault fracture zone, and tunnel excavation space are treated in the fluid computation domain according to a unified framework of "porous media - free flow": Figure 2 In this region, the surrounding rock matrix zone is approximated using Darcy flow, while the fault fracture zone dynamically updates porosity parameters based on coupled calculations, allowing for non-Darcy effects. The tunnel excavation space zone is solved using the Navier-Stokes governing equations when free flow occurs. In this embodiment, the parameters of the surrounding rock matrix zone are fixed values, while the parameters of the fault fracture zone and tunnel excavation space zone are dynamically updated with the coupling cycle. Parameter selection is based on the interval corresponding to the non-Darcy effect E value, and the updated parameters are then used. , , The drag force terms Fx, Fy, and Fz are attached to the Fluent porous medium and source terms via UDF.

[0083] Table 2 Parameters of Fluid Phase Porous Media

[0084] To balance computational efficiency with the microscopic response resolution of fault-tunnel, the PFC3D particle system adopts a zoned and graded particle size setting: establishing particles for the surrounding rock, fault, tunnel, and surrounding rock-fault transition zone, and setting small / medium / large particle size combinations in each particle zone. Relatively smaller particle sizes are used around the tunnel and fault, while the particle size is appropriately increased in areas far from the tunnel and fault to improve the simulation accuracy of local particle migration and pore evolution.

[0085] Considering the mud bonding effect in the soil-rock mixture, a linear parallel bonding model is used for the contact between particles inside the fault and normal rock particles; a linear contact model is used between particles and the wall. Gravity is applied in the model to achieve initial in-situ stress equilibrium; for particles located below the groundwater level, buoyancy correction is introduced to reflect the effective unit weight under saturation. After achieving in-situ stress equilibrium, displacement is removed to avoid confusion caused by the superposition of displacements during subsequent excavation.

[0086] The tunnel excavation process was treated as a one-time excavation: the total axial excavation length was 92m, including 70m of normal rock and 22m of fault rock; the fault fracture zone was exposed at y=92m. To reflect the lack of local support at the tunnel face, the excavation face within a range of 0.5 times the tunnel cross-section diameter, approximately 7.75m, was set as an unsupported section; the remaining tunnel wall length of approximately 84.25m was set as a one-time formed lining support, implemented as a wall in the lining model, and its geometric and mechanical parameters are given in Table 3. The total number of particles in the DEM model in this embodiment is 56,234. After the excavation and support settings were completed, a period of balancing was performed to allow the DEM particle model to reach a static equilibrium state, and no significant deformation occurred in the fault surrounding rock in front of the tunnel face; then the displacement was cleared again as the initial moment for coupled calculation.

[0087] Table 3. Basic parameters and region grouping settings of the DEM particle system (PFC3D 6.30.0)

[0088] Note: Particle density is based on... calculate, The density of the rock (the density of the rock is 2000 kg / m³). 3 The density of the fault is 1850 kg / m³. 3 ), Porosity.

[0089] Table 4 DEM Contact Model and Parameter Settings

[0090] Note: The coefficient of contact friction is calculated according to... calculate, The friction angle of the rock and soil mass is denoted as .

[0091] Based on the engineering water level conditions, the groundwater level is initialized in the calculation domain, approximately 50m above the ground surface. The specific elevation is determined according to... Figure 2 Given the following initial fluid states: the tunnel region is initially filled with air; the water-bearing region outside the tunnel and the fault-supplying region are initially water-phase. Specifically, after initialization, a local initialization / patch is applied to the region outside the tunnel, setting the water phase volume fraction to 1 and the air phase volume fraction to 0, while setting the air phase volume fraction to 1 for the tunnel region. The upper boundary of the computational domain is set as a free water surface boundary to reflect the connectivity between the surface water level and the atmosphere; the remaining boundaries are set as impermeable boundaries to ensure that water flow is mainly controlled by gravity head and fault seepage channels entering the tunnel. The boundary types and numerical values ​​corresponding to each boundary in Fluent are given in Table 5. A pressure outlet condition is used at the tunnel exit to represent the pressure release condition between the excavated space and the external atmosphere.

[0092] Table 5 Fluent Boundary Conditions and Initial Conditions Settings

[0093] The fluid computational domain employs transient solution. The pressure-velocity PISO coupling algorithm, VOF interface reconstruction format, time step Δt, maximum number of iterations per step, and residual convergence criterion, among other solution control parameters, are set according to the implementation example, as detailed in Table 6. To achieve bidirectional CFD-DEM coupling, this embodiment defines the coupling period as every... The fluid time step performs a data exchange and parameter update: PFC3D outputs microscopic field quantities such as porosity, permeability, and drag force, while Python calculates the drag coefficient and source term and completes the flow regime partitioning. Then, it updates to Fluent via UDF and enters the next coupling cycle.

[0094] Table 6 Solving for control and coupling parameters

[0095] (3) Coupling implementation; Numerical analysis of water and mud inrush processes in fault fracture zones is achieved based on CFD-DEM bidirectional coupling. By establishing a linkage calculation mechanism between the fluid field and the particle field, the seepage state, particle migration, and pore structure evolution can be updated collaboratively within the same framework. This method can distinguish different flow regimes based on local flow characteristics and implement corresponding parameter updates and equation configurations, thereby achieving a continuous transition between porous media seepage and free flow in tunnel excavation space, improving the applicability and stability of the whole disaster simulation.

[0096] In this embodiment, a Python environment was built using Fluent 2024R2, PFC3D 6.00.30 and Python 3.12, and a C language environment was built using Visual Studio 2022 to perform CFD-DEM fluid-structure interaction calculations.

[0097] 1) Unit-level consistency association: To ensure spatial consistency of data write-back between the fluid computational domain and the particle domain, a cell association method based on geometric location features is adopted to establish the correspondence between computational cells on both sides using cell spatial location identifiers.

[0098] 2) Two-way data exchange and coupling cycle advancement: This embodiment uses the itasca module provided by the PFC3D embedded Python environment to implement TCP socket communication. It uses itasca.util.p2pLinkServer and itasca.util.p2pLinkClient to establish a connection between the server and client Python environments. It uses send_data() and read_data() to transmit numerical values, strings and NumPy arrays between the two Python instance objects to achieve bidirectional exchange of the field quantities required for coupling.

[0099] The coupling calculation is executed cyclically according to the preset coupling cycle. The fluid side outputs field quantities such as velocity, pressure and gradient for flow regime discrimination and fluid-solid interaction calculation. The particle side outputs pore structure characterization quantities for seepage parameter update. The updated parameter field is then fed back to the fluid side to participate in the solution of the next cycle. The coupling cycle and update frequency are executed according to the parameter settings given in the embodiment.

[0100] 3) Multi-flow-state partitioning identification based on the non-Darcy effect index E: The Darcy flow, non-Darcy flow, and free flow that may occur during the water inrush evolution process are regarded as a multi-flow process that can be continuously transformed. In each coupling cycle, the non-Darcy effect index E is calculated and the flow state is partitioned according to this index. This drives the dynamic configuration of the control equation terms and porous media parameters, realizing the continuous transition between porous media seepage and free flow in the tunnel excavation space. For the free flow region in the tunnel excavation space, auxiliary criteria such as regional attributes or proportion can be used to constrain it, so as to improve the physical consistency and stability of the partition boundary.

[0101] 4) Dynamic update of parameter field and partition mapping: Dynamic updates involve two aspects: first, the evolution of the pore structure of the granular system drives the updating of seepage parameters; second, based on the flow regime partitioning results, the governing equations and porous media parameter terms are dynamically configured, enabling continuous transition and switching between different flow regimes within a unified framework. In the coupled calculation, the pore structure characterization quantity (such as porosity) is output from the particle side, from which the permeability parameters and corresponding drag coefficients are calculated, and necessary scale correlations are established with engineering scales. Subsequently, the parameters of each region are standardized according to the flow regime partitioning results, forming an updated parameter field, which is then written back to the fluid side for subsequent solutions. Through the above closed-loop mechanism of "structural evolution—parameter update—partition configuration—write-back solution," the adaptive updating of porous media seepage parameters with particle migration is achieved, ensuring that the equation configuration of different flow regime regions matches their local flow states.

[0102] 5) Fluid-side runtime update mechanism: To achieve real-time updates of coupling parameters during the fluid solution process, a user-defined function (UDF) provided by the solver enables runtime write-back of external parameters. A user-defined memory (UDM) is set up within the solver to store updated variables such as porosity, drag coefficient, and momentum source term. A parameter update operation is triggered in each coupling cycle, writing the updated parameter field to the corresponding computational unit. Subsequently, the solver calls this stored variable in the calculation of porous medium parameters and source terms, dynamically assigning values ​​to the drag and source terms. For parallel computing scenarios, a consistent distribution and write strategy is adopted to ensure that each computational partition obtains consistent parameter update results, thereby reducing the risks of write-back mismatches, asynchrony, and the resulting numerical instability. This runtime update mechanism supports dynamic parameter configuration and continuous coupled computation driven by multiple flow regime partitions.

[0103] (4) Simulation process and effect verification; To verify the feasibility and technical effectiveness of the method of the present invention in simulating the entire process of water and mud inrush in tunnels, this embodiment, based on the engineering model, parameter settings and coupling implementation described in sections (1) to (3), provides a simulation process, comparative examples and quantitative verification indicators. The coupling simulation process of this embodiment is executed according to the following steps, and the overall process is shown in [link to full description]. Figure 3.

[0104] 1) Simulation process and operation steps: (I) Model initialization and static equilibrium. Gravity is applied in PFC3D and the buoyancy correction of particles below the water level is considered to perform ground stress equilibrium calculation; after equilibrium is completed, the displacement is cleared and saved as the initial state of the coupling calculation. Then, the excavation and support settings are set as described in (2) above to form the tunnel excavation space and lining boundary; after excavation, the equilibrium is continued until the system reaches static stability, and the displacement is cleared again as the starting time of the coupling calculation.

[0105] (II) Fluid field initialization and initial phase distribution setting. Enable the VOF multiphase flow model (VOF two-phase flow model in this embodiment) in Fluent and complete the standard initialization; set the tunnel area to an air phase volume fraction of 1, and set the water-bearing area outside the tunnel to a water phase volume fraction of 1 through local initialization, so as to form the initial phase distribution of the gas phase inside the tunnel, the fault and the water phase outside the surrounding rock.

[0106] (III) Initial Coupling Parameter Writing. The initial period is based on the porosity output by PFC3D. viscous resistance formed by permeability k and its reciprocal Inertial drag The drag force terms Fx, Fy, Fz, etc., are written into Fluent, and the porous medium properties and source terms are loaded. During the initial coupling stage, since there is no stable velocity field output, the three-flow-state discrimination based on the non-Darcy effect index E is not performed. The partitions are initialized according to the regional attributes, with the surrounding rock area being the Darcy area and the fault area and tunnel neighborhood being the areas to be updated.

[0107] (IV) Coupled Cycle Propulsion. Entering the coupled cycle: Fluent solves the propulsion flow field transiently and outputs velocity, pressure, pressure gradient, and water phase volume; PFC3D solves the propulsion particle system using the DEM and updates the pore structure; the Python side calculates viscous drag, inertial drag, drag force, and non-Darcy effect index E in each coupled cycle, completes the DFNS (Darcy–Forchheimer–Navier–Stokes) three-flow-state identification and dynamic partitioning configuration, and generates a data file; then, UDF is called to complete parameter write-back, entering the next coupled cycle. This cycle continues until the total simulation time T=300s.

[0108] (V) Output and Monitoring. Time history information for key sections and monitoring points was recorded during the simulation. Results and analysis are shown below. Figure 4-13 Including: water phase volume fraction in the tunnel The evolution, velocity and pressure changes within the tunnel excavation space, flow rate at the fault passage, and spatial distribution and temporal variation of the non-Darcy effect index E in the fault zone are used for subsequent effect verification and comparative analysis.

[0109] 2) Definition of verification indicators: The effectiveness of this embodiment is verified from two dimensions: physical consistency and the characterizability of disaster evolution. First, regarding physical consistency, by identifying and dynamically configuring the non-Darcy effect index E in different regions, the porous resistance terms in different regions are matched with the free-flow control equations, thereby suppressing local non-physical velocity anomalies and making the spatial distribution of the velocity and pressure fields more reasonable; the relevant results are illustrated by velocity field distribution maps and pressure field distribution maps, respectively. Second, regarding disaster evolution characterization, this embodiment utilizes the water phase volume fraction... The velocity and pressure time history curves of the cloud map and the monitoring section at the tunnel face characterize the propulsion characteristics and flow response of the water inrush process. Furthermore, by combining the velocity and pressure time history curves of representative monitoring points in different flow regime zones, the dynamic response evolution law of key areas during the water inrush propulsion process is verified. In summary, the results show that the method of this invention can achieve a two-phase interface propulsion description under the condition of water inrush entering the tunnel air domain, and realize multi-flow regime partitioning and dynamic switching in coupled calculations, thereby improving the physical consistency and process characterization of the entire disaster simulation.

[0110] The above-disclosed embodiments are merely preferred embodiments of the present invention. These embodiments are only used to illustrate the technical solutions of the present invention and should not be construed as limiting the scope of the present invention. Those skilled in the art can understand that implementing all or part of the processes of the above embodiments and making equivalent changes in accordance with the claims of the present invention still fall within the scope of the invention.

Claims

1. A computational method for multi-flow regime partition identification and dynamic updating of DFNS based on VOF–CFD–DEM coupling, characterized in that, Includes the following steps: Construct a fluid-solid-gas multi-field coupled computational domain and a discretized model; A unified governing equation system for DFNS is established based on the volume average method. The medium properties in different regions of the computational domain are characterized by parameterization of porosity and porous medium resistance source terms and volume force source terms. Based on the medium characteristics of different regions within the computational domain, flow regime partition labels are generated using an established multi-flow regime criterion system. Based on the flow regime partition labels, parameter mapping logic is formulated to enable adaptive switching of the DFNS unified control equation system; The DFNS unified control equation system is discretized and solved using real physical velocities. The discretized equation system is iteratively updated at each time step using a pressure-velocity coupling algorithm until the preset convergence condition is met, thereby obtaining the basic field variable input parameters for fluid-particle interaction force calculation. The drag force of the fluid is calculated using a drag force model that considers local porosity correction, and the translational and rotational states of the particles are updated. The sum of particle volumes in each updated fluid unit is statistically analyzed to update the unit porosity. At the same time, the relationship between porosity and the average permeability of the fluid unit is established. Porosity and permeability are reconstructed and updated at the end of each coupling step until the current simulation time or coupling step number reaches the preset termination condition, and the optimal parameter combination is output.

2. The calculation method for DFNS multi-flow regime partition identification and dynamic update based on VOF–CFD–DEM coupling according to claim 1, characterized in that, The steps for constructing the fluid-solid-gas multi-field coupled computational domain and discretization model include: Establish a global computational domain that includes the surrounding rock matrix, fault fracture zone, and tunnel excavation space; Within the global computational domain, the fluid computational domain is discretized using the finite volume method, and a discrete element model is constructed for the solid particle system in the fault fracture zone, forming a coupled computational model of the fluid control volume and the particle assembly. The VOF multiphase flow model is introduced into the fluid computation domain to achieve a unified flow field simulation of the gas phase in the tunnel excavation space, as well as the liquid phase in the surrounding rock and fracture zone.

3. The calculation method for DFNS multi-flow regime partition identification and dynamic update based on VOF–CFD–DEM coupling according to claim 2, characterized in that, The steps for introducing a VOF multiphase flow model into the fluid computational domain include: A discretized model is established based on the local volume averaging method, and porosity is defined in the discretized model. The phase volume fraction from the VOF multiphase flow model is introduced into the discretized model to weight the calculation of the mixture properties parameters of the gas-liquid two-phase mixture; For the gas-liquid interface during the water inrush process, a continuous surface force model is used to calculate the surface tension.

4. The calculation method for DFNS multi-flow regime partition identification and dynamic update based on VOF–CFD–DEM coupling according to claim 3, characterized in that, In the unified governing equation system of DFNS, based on the differences in the values ​​of porosity, porous media resistance source terms, and momentum source terms, the unified momentum conservation equation corresponds to different flow regime characteristics in different regions of the computational domain, and its expression is as follows: In the formula: For fluid effective occupancy rate, The density of the mixed fluid within the grid cell. For numerical computation time, For hybrid viscosity, u is the fluid physical velocity vector, p is the fluid pressure, g is the gravitational acceleration vector, and S is the porous medium resistance source term; For vector differential operators, For fluid pressure gradient, Let be the fluid velocity gradient tensor. For volumetric force source terms; Among them, the fluid motion in the surrounding rock matrix zone is driven by the pressure gradient, dominated by viscous resistance, and the unified momentum conservation equation degenerates into the classical Darcy's law. The pore structure in the fault fracture zone is heterogeneous, and the increase in flow velocity leads to a significant local inertial effect. The pressure drop and flow velocity have a nonlinear relationship, and the unified momentum conservation equation degenerates into the Forchheimer nonlinear seepage equation. The actual flow velocity of the fluid in the free flow zone of the tunnel excavation space increases significantly between the particle skeleton due to the reduction of the flow cross section. The unified momentum conservation equation is solved by using the Navier-Stokes equation with porosity term and combined with the standard k-epsilon turbulence model.

5. The calculation method for DFNS multi-flow regime partition identification and dynamic update based on VOF–CFD–DEM coupling according to claim 3, characterized in that, The steps to establish a multi-flow criterion system include: The Fochheimer number Fo is selected to characterize the nonlinear intensity of fluid motion in porous media; Based on the Fochheimer number Fo, a non-Darcy effect index E is introduced to normalize the nonlinearity of the flow to the [0,1] interval, quantitatively defining the flow regime transition boundary; the discrimination criteria of the multi-flow regime criterion system are as follows: when At this time, the inertial term can be ignored, and it is determined to be Darcy laminar flow; when At that time, the corresponding critical Fochheimer number At that time, the inertial resistance was significant, and it was determined to be Forchheimer nonlinear seepage; when At this point, the influence of the viscous term weakens, and the flow is completely dominated by inertia, which is determined to be Navier-Stokes free turbulence.

6. The calculation method for DFNS multi-flow regime partition identification and dynamic update based on VOF–CFD–DEM coupling according to claim 5, characterized in that, In the multi-flow criterion system, a cross-scale coupling architecture and mesh mapping of fluid-particle are constructed, and a velocity field correction algorithm based on phase threshold is used to correct the false gas phase drag in the fluid-solid-gas multi-field coupling. When the water phase volume fraction of the mesh cell in the fluid computation domain is determined to be a gas phase dominant cell, a velocity zeroing correction operation is performed on the fluid velocity.

7. The calculation method for DFNS multi-flow regime partition identification and dynamic update based on VOF–CFD–DEM coupling according to claim 5, characterized in that, The parameter mapping logic based on the flow regime partition labels is as follows: When the flow is determined to be Darcy laminar flow, a flow regime label "Darcy" is generated, linear drag is activated, nonlinear drag is deactivated, and the porous media drag source term degenerates into... ; When the flow is determined to be non-Darcy flow, a flow regime label (Forchheimer) is generated, and both linear drag and inertial drag are activated. The drag source term in porous media degenerates into... ; When the flow is determined to be free turbulence, a flow regime label Navier-Stokes is generated, the momentum source term is set to zero, and the equations are automatically reduced to Navier-Stokes equations.

8. The calculation method for DFNS multi-flow regime partition identification and dynamic update based on VOF–CFD–DEM coupling according to claim 5, characterized in that, The steps to obtain the basic field variable input parameters for fluid-particle force calculations include: Construct a multiphase flow control equation that includes porosity parameters based on the DFNS unified control equation system; The multiphase flow control equations, which include porosity parameters, are iteratively updated at each time step using a pressure-velocity coupling algorithm. After convergence, the mixed pressure field, velocity field, pressure gradient field, and gas-liquid phase field information of the entire field grid are obtained, which serve as the input parameters of the basic field variables for calculating fluid-particle interaction forces.

9. The calculation method for DFNS multi-flow regime partition identification and dynamic updating based on VOF–CFD–DEM coupling according to claim 8, characterized in that, The expression for the drag force model considering local porosity correction is as follows: ； ； In the formula: The forces acting on the particles in the gas-liquid-solid three-phase coupled field. For fluid drag force, For the drag force of a single particle, The average porosity of the fluid unit containing the particle. This is an empirical factor for correcting the local average porosity, where r is the particle radius. For fluid pressure gradient, This is the gravitational acceleration vector.

10. The calculation method for DFNS multi-flow regime partition identification and dynamic updating based on VOF–CFD–DEM coupling according to claim 9, characterized in that, The steps to output the optimal parameter combination include: A particle-scale similarity and parameter equivalence correction mechanism is introduced. At the end of each coupling step, the porosity and permeability are reconstructed and updated. The updated medium parameters are fed back to the fluid solution end as input parameters for the next coupling step momentum source term construction and flow regime identification calculation. Determine whether the current simulation time or coupling step count has reached the preset termination condition. If it has, output the optimal parameter combination. If it has not, enter the next coupling loop until the preset termination condition is reached, and then output the optimal parameter combination.

Images