Underwater shield tunnel catastrophe simulation method based on CFD-DEM

Abstract

The invention discloses an underwater shield tunnel catastrophe simulation method based on CFD-DEM, and relates to the technical field of tunnel and underground engineering, and the method comprises the steps: constructing a fluid CFD model of an underwater shield tunnel, a BPM contact model concrete lining model, and a soil body DEM model based on soil body particles; carrying out particle mesoscopic calibration and constructing an integral DEM model considering the interaction between the soil body and the structure; and setting a time step length of the fluid CFD model and the overall DEM model and a coupling frequency between the fluid CFD model and the overall DEM model, and carrying out coupling calculation to obtain an underwater shield tunnel catastrophe simulation result. According to the method, by setting the time step length and the coupling gap, the whole process of catastrophe of the underwater shield tunnel can be clearly observed according to the simulation result, fracture and non-fracture distinguishing can be carried out on a Bond key in a final BPM model of the concrete lining, and the damage occurrence area of the lining is observed.

Description

Technical Field

[0001] The present invention relates to the technical field of tunnel and underground engineering, and more specifically, to a method for simulating disasters of underwater shield tunnels based on CFD-DEM. Background Art

[0002] The shield method has become the most main construction method for urban subways and underwater tunnels due to its safety and high efficiency. Compared with land tunnels, underwater tunnel construction faces more complex geological conditions, such as soft upper and hard lower composite strata, fault fracture zones, and waters with infinite recharge. When the shield machine crosses such strata, the shield attitude control is poor, and the jacks are prone to uneven thrust, resulting in a relatively high degree of local stress concentration in the segments, causing problems such as segment cracking and damage. In addition, the high underwater osmotic pressure will cause further development of the initial damage of the segments, thus affecting the construction quality of underwater tunnels. In addition, natural disasters such as earthquakes also pose great safety hazards to the tunnel structure. To explore the "lining-soil" integrity structure failure problem caused by underwater high osmotic pressure and natural disasters such as earthquakes in underwater shield tunnels, an efficient and accurate simulation and analysis method is urgently needed.

[0003] There are existing simulation methods for analyzing the seismic damage of underwater shield tunnels, but there are still a series of problems that need to be solved urgently, including: (1) Most of the existing methods for simulating underwater shield tunnels use the finite element method, mainly focusing on macroscopic phenomena and rarely considering the damage development at the mesoscopic level; (2) The existing finite element models for simulating underwater shield tunnels are difficult to truly restore the seawater flow and overburden conditions, and cannot truly restore the disaster damage under fluid-structure coupling; (3) Most of the existing discrete element methods for simulating tunnel structure disasters do not consider the influence of fluids and only consider the particle model; (4) Most of the existing discrete element methods for simulating tunnel structure disasters rarely consider the influence of the "lining-soil" interaction and mainly focus on the loss of soil particles.

[0004] Therefore, how to propose a method for simulating disasters of underwater shield tunnels based on CFD-DEM, considering the "lining-soil" interaction, seawater-groundwater, and soil-structure coupling, and better simulating the disasters of underwater shield tunnels under complex geological conditions is an urgent problem to be solved by those skilled in the art. Summary of the Invention

[0005] In view of this, the present invention provides a method for simulating disasters of underwater shield tunnels based on CFD-DEM, considering the "lining-soil" interaction, seawater-groundwater, and soil-structure coupling, and better simulating the disasters of underwater shield tunnels under complex geological conditions. To achieve the above object, the present invention adopts the following technical solutions:

[0006] A method for simulating disasters of underwater shield tunnels based on CFD-DEM includes:

[0007] Define the fluid using the VOF method and construct a fluid CFD model for the underwater shield tunnel;

[0008] Establish a BPM contact model for the concrete lining model to simulate the mesoscopic development of the tunnel structure damage under seawater pressure and seismic action;

[0009] Define the concrete lining model as a porous medium model, conduct tests on the porosity and permeability coefficient of the concrete standard specimen, and then calculate the permeability coefficient and porosity of the porous medium model based on the porosity and permeability coefficient obtained from the tests;

[0010] Adopt the cohesive force Hertz contact to establish a DEM model of the soil based on soil particles to simulate the interaction between the soil and the structure and the ground settlement under seawater pressure and seismic action;

[0011] Conduct uniaxial compression tests on concrete, triaxial compression tests on soil, and corresponding DEM numerical simulations according to the BPM contact model concrete lining model, compare the numerical simulation results with the test results, and perform particle mesoscopic calibration;

[0012] Construct an overall DEM model considering the interaction between the soil and the structure according to the particle mesoscopic calibration results;

[0013] Set the time step of the fluid CFD model and the overall DEM model, as well as the coupling frequency between the fluid CFD model and the overall DEM model, and perform coupling calculations to obtain the catastrophe simulation results of the underwater shield tunnel.

[0014] Optionally, the calculating the permeability coefficient and porosity of the porous medium model according to the porosity and permeability coefficient obtained from the tests includes: calculating the Darcy seepage coefficient and viscous resistance coefficient of the porous medium model according to the existing empirical formula and the formula in OpenFOAM, specifically as follows:

[0015]

[0016] Among them, S1 represents the method for calculating the Darcy seepage coefficient and viscous resistance coefficient in the empirical formula, S2 represents the method for calculating the Darcy seepage coefficient and viscous resistance coefficient when defining the porous medium model in OpenFOAM, μ represents the dynamic viscosity of the fluid, φ is the porosity, K is the permeability coefficient, ρ is the fluid density, F is the Forchheimer coefficient, u is the flow velocity, d is the Darcy coefficient, and f is the Forchheimer coefficient.

[0017] Optionally, it further includes: performing CFD mesh division on the fluid CFD model of the underwater shield tunnel, where the fluid CFD model of the underwater shield tunnel is divided into CFD meshes according to the tunnel structure area, soil area, and the air area inside the tunnel.

[0018] Optionally, defining the fluid by the VOF method and constructing the fluid CFD model of the underwater shield tunnel includes:

[0019] Defining multiple fluids by the VOF method, simulating the overlying seawater pressure by setting the pressure boundary of the pressure fluid, and simulating the infinite replenishment of seawater by continuously injecting the fluid;

[0020] Setting the pressure difference between the seawater inlet and the air outlet under the initial region and initial boundary conditions to meet the pressure control, and constructing the fluid CFD model of the underwater shield tunnel.

[0021] Optionally, the calculation method of the seawater pressure includes:

[0022] p = ρgh;

[0023] where ρ is the density of seawater; h is the depth of the overlying seawater; g is the acceleration due to gravity.

[0024] Optionally, solving the motion of the fluid phase in the fluid CFD model:

[0025]

[0026] where φ f is the volume fraction of the fluid; is the fluid velocity; ρ f is the fluid density; is the fluid pressure gradient; is the stress tensor gradient; g is the acceleration due to gravity in the research area; is the particle-fluid interaction force, is the gradient operator.

[0027] Optionally, the mesoscopic parameters of the BPM contact model concrete lining model include:

[0028] F i = F n n i + F s t i ;

[0029] In the formula, F i is the interaction force between particles, F n and F s respectively represent the normal and tangential force components, n i and t i respectively represent the plane unit vectors, and the calculation of the normal force is:

[0030] F n = K n U n ;

[0031] In the formula, U n is the overlapping amount between particles, and K n is the contact normal stiffness, which can be calculated by the following formula:

[0032]

[0033] In the formula, and are the normal stiffnesses of the particles.

[0034] Optionally, the motion control equations of the particles in the soil body DEM model include:

[0035]

[0036] where m pi , I pi , are respectively the mass, moment of inertia, linear velocity, and angular velocity of the i-th particle; n c is the number of particles in contact with the i-th particle; is the contact force between the i-th particle and the j-th particle; is the gravity force received by the i-th particle; is the fluid-particle interaction force received by the i-th particle; is the torque generated on the i-th particle by the interaction between the i-th particle and the j-th particle.

[0037] Optionally, the steps of performing uniaxial compression tests on concrete, triaxial compression tests on soil, and corresponding DEM numerical simulations, and comparing the numerical simulation results with the test results for particle mesoscopic calibration include:

[0038] Determine the particle size distribution of the original soil sample, conduct triaxial compression tests on the soil, and obtain the stress-strain relationship curve, Mohr stress circle, and strength envelope of the original soil sample;

[0039] Conduct DEM numerical simulations of triaxial compression tests on soil by the discrete element method;

[0040] Use C55 cement and multi-graded aggregates to cast a cube with a side length of 150 mm for uniaxial compression tests, and obtain the failure mode and stress-strain curve of the concrete specimen;

[0041] Use the discrete element method to establish a corresponding numerical model to carry out uniaxial compression simulations and compare with the test results;

[0042] Establish an overall DEM model, and calibrate the particle parameters required for the overall DEM model, the mesoscopic parameters of the concrete lining BPM model, and the fracture parameters of the corresponding bond bonds.

[0043] Optionally, constructing an overall DEM model considering the interaction between soil and structure according to the particle mesoscopic calibration results further includes:

[0044] The grid of the concrete lining model area of the BPM contact model is encrypted, and the fluid grid size is set to 2-3 times that of the particles;

[0045] In the non-seepage area, particles with a diameter 2-3 times the average diameter of the particles in the area near the tunnel are selected.

[0046] Optionally, the coupled calculation includes:

[0047] The fluid CFD model performs independent calculations according to the time step of the fluid CFD model;

[0048] The overall DEM model performs independent calculations according to the time step of the overall DEM model;

[0049] The fluid CFD model and the overall DEM model perform coupled calculations according to the coupling frequency between the fluid CFD model and the overall DEM model. After updating the particle-fluid interaction force and particle positions, the independent calculations of the fluid CFD model and the overall DEM model are performed respectively.

[0050] It can be seen from the above technical solutions that compared with the prior art, the present invention discloses a method for simulating underwater shield tunnel disasters based on CFD-DEM, which has the following beneficial effects:

[0051] The present invention proposes a method for simulating disasters of underwater shield tunnels based on CFD-DEM, including: defining fluids using the VOF method and constructing a fluid CFD model of the underwater shield tunnel; establishing a concrete lining model with a BPM contact model to simulate the mesoscopic development of damage to the tunnel structure under seawater pressure and seismic action; defining the concrete lining model as a porous medium model, conducting tests on the porosity and permeability coefficient of concrete standard specimens, and then calculating the permeability coefficient and porosity of the porous medium model based on the porosity and permeability coefficient obtained from the tests; establishing a soil DEM model based on soil particles using cohesive Hertz contact to simulate the interaction between the soil and the structure and ground settlement under seawater pressure and seismic action; conducting uniaxial compression tests on concrete, triaxial compression tests on soil, and corresponding DEM numerical simulations according to the concrete lining model with the BPM contact model, comparing the numerical simulation results with the test results, and performing particle mesoscopic calibration; constructing an overall DEM model considering the interaction between the soil and the structure according to the particle mesoscopic calibration results; setting the time step of the fluid CFD model and the overall DEM model and the coupling frequency between the fluid CFD model and the overall DEM model, and performing coupling calculations to obtain the simulation results of disasters of the underwater shield tunnel. By setting the time step and coupling gap, the present invention enables the entire process of disasters of the underwater shield tunnel to be clearly observed in the simulation results, and can distinguish between broken and unbroken Bond bonds in the BPM model of the final concrete lining, observe the damage occurrence area of the lining, and at the same time, other data of the soil or the lining can be extracted. Description of the Drawings

[0052] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for use in the description of the embodiments or the prior art. Obviously, the drawings described below are only the embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on the provided drawings without creative efforts.

[0053] Figure 1 It is a schematic flow chart of a method for simulating disasters of an underwater shield tunnel based on CFD-DEM provided by the present invention.

[0054] Figure 2 It is a schematic plan view of the geometric model of the CFD-DEM underwater shield tunnel provided by the present invention.

[0055] Figure 3 It is a schematic diagram of the DEM numerical simulation of the uniaxial compression test provided by the present invention.

[0056] Figure 4 It is a schematic diagram of the DEM numerical simulation of the triaxial compression test provided by the present invention.

[0057] Figure 5 Schematic diagram of the triaxial compression simulation specimen of the soil mass in the underwater shield tunnel provided by the present invention.

[0058] Figure 6 Schematic diagram of the DEM numerical simulation of the soil-structure interaction in the underwater shield tunnel provided by the present invention. Specific embodiments

[0059] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0060] The embodiment of the present invention discloses a method for simulating disasters in an underwater shield tunnel based on CFD-DEM, including: defining fluids by the VOF method and constructing a fluid CFD model of the underwater shield tunnel; establishing a BPM contact model concrete lining model to simulate the mesoscopic development of damage to the tunnel structure under seawater pressure and seismic action; defining the concrete lining model as a porous medium model, conducting tests on the porosity and permeability coefficient of concrete standard specimens, and then calculating the permeability coefficient and porosity of the porous medium model according to the porosity and permeability coefficient obtained from the tests; establishing a soil DEM model based on soil particles by using cohesive force Hertz contact to simulate the interaction between the soil and the structure and the ground settlement under seawater pressure and seismic action; conducting uniaxial compression tests on concrete, triaxial compression tests on soil, and corresponding DEM numerical simulations according to the BPM contact model concrete lining model, comparing the numerical simulation results with the test results, and performing particle mesoscopic calibration; constructing an overall DEM model considering the interaction between the soil and the structure according to the particle mesoscopic calibration results; setting the time step of the fluid CFD model and the overall DEM model and the coupling frequency between the fluid CFD model and the overall DEM model, and performing coupling calculations to obtain the simulation results of disasters in the underwater shield tunnel.

[0061] In the specific embodiments, a method for simulating disasters in an underwater shield tunnel based on CFD-DEM includes the following specific steps:

[0062] For the CFD geometric model of the underwater shield tunnel, the VOF method is used to define various fluids, including: seawater-air, seawater-air-groundwater, and the overlying seawater pressure on the ground surface is simulated through fluid pressure, and the infinite replenishment of seawater is simulated by continuously injecting seawater into the fluid domain, thereby constructing the fluid CFD model of the underwater shield tunnel. The plane structure of the CFD-DEM underwater shield tunnel model geometric model is as Figure 2 shown;

[0063] Concrete lining BPM model, a concrete lining model with bond bonds is established based on the BPM contact model, which can simulate the mesoscopic development of damage to the tunnel structure under seawater pressure and seismic action;

[0064] Define the porous medium model of the concrete lining and calculate the Darcy seepage coefficient and viscous resistance coefficient;

[0065] Soil DEM model, a soil particle model is established using cohesive Hertz contact, which can simulate the interaction between the soil and the structure and the ground settlement under seawater pressure and seismic action;

[0066] Calibration of particle mesoscopic parameters, including uniaxial compression tests of concrete, triaxial compression tests of soil, and corresponding DEM numerical simulations, and comparing the numerical simulation results with the test results to achieve the calibration of particle mesoscopic parameters;

[0067] According to the calibration results of the particle mesoscopic parameters, construct the overall DEM model of the lining-soil considering soil-structure interaction. The DEM numerical simulation of the "soil-structure" interaction of the underwater shield tunnel is as Figure 6 shown;

[0068] Set the time step of the fluid CFD model, the time step of the overall lining-soil DEM model, and the coupling frequency between the fluid CFD model and the overall lining-soil DEM model, and perform coupling calculations until the set end time is reached to obtain the catastrophe simulation results of the underwater shield tunnel.

[0069] Furthermore, calculate the Darcy seepage coefficient and viscous resistance coefficient of the porous medium model according to the existing empirical formula and the formula in OpenFOAM, as follows:

[0070]

[0071] Among them, S1 represents the method for calculating the Darcy seepage coefficient and viscous resistance coefficient in the empirical formula, S2 represents the method for calculating the Darcy seepage coefficient and viscous resistance coefficient when defining the porous medium model in OpenFOAM, μ represents the dynamic viscosity of the fluid, φ is the porosity, K is the permeability coefficient, ρ is the fluid density, F is the Forchheimer coefficient, u is the flow velocity, d is the Darcy coefficient, and f is the Forchheimer coefficient.

[0072] Furthermore, perform CFD mesh division on the CFD geometric model of the underwater shield tunnel, specifically:

[0073] The CFD geometric model of the underwater shield tunnel is used to set and analyze the influence of the fluid;

[0074] The CFD geometry model of the underwater shield tunnel is meshed by CFD according to the tunnel structure area, soil area, and the air area inside the tunnel;

[0075] Furthermore, the VOF method is used to define multiple fluids, including: seawater-air, seawater-air-groundwater. The overlying seawater pressure is simulated by setting the pressure boundary of the pressure fluid, and the continuous injection of fluid is used to simulate the infinite replenishment of seawater. Specifically:

[0076] Set the initial area, including: the initial saturated subsoil area, overlying seawater area, and air area;

[0077] Set the pressure difference between the seawater inlet and air outlet under the initial boundary conditions to meet the pressure control.

[0078] Furthermore, the calculation method of the seawater pressure is as follows:

[0079] p = ρgh;

[0080] Where ρ is the density of seawater; h is the depth of the overlying seawater; g is the acceleration due to gravity.

[0081] Furthermore, through the comparison of the results of the triaxial compression test of the original soil sample, the uniaxial compression test of the concrete lining, and the DEM numerical simulation of the corresponding compression test, the particle parameters required for calibrating the soil DEM model are determined. The schematic diagram of the DEM numerical simulation of the uniaxial compression test is as Figure 3 shown, the schematic diagram of the DEM numerical simulation of the triaxial compression test is as Figure 4 shown, and the schematic diagram of the triaxial compression simulation specimen of the underwater shield tunnel soil is as Figure 5 shown. Specifically:

[0082] Determine the particle size distribution of the original soil sample, conduct the triaxial compression test, and obtain the stress-strain relationship curve, Mohr stress diagram, and strength envelope of the original soil sample;

[0083] Conduct the DEM numerical simulation of the triaxial compression test by the discrete element method;

[0084] Use C55 cement and multi-graded aggregates to cast a cube with a side length of 150 mm for the uniaxial compression test, and obtain the failure mode and stress-strain curve of the concrete specimen;

[0085] Use the discrete element method to establish a corresponding numerical model to carry out the uniaxial compression simulation and compare it with the test results;

[0086] Establish the overall DEM model of the lining-soil, and calibrate the particle parameters required for the soil DEM model, the mesoscopic parameters of the BPM model, and the fracture parameters of the corresponding bond.

[0087] Further, after generating the particle model within the required research scope based on the particle calibration results, it further includes:

[0088] The mesh of the concrete lining model area is encrypted, and the fluid mesh size is set to 2 - 3 times that of the particles;

[0089] In the non - seepage area, particles with a diameter 2 - 3 times the average diameter of the particles in the area near the tunnel are selected.

[0090] Further, the coupled calculation is specifically as follows:

[0091] The fluid CFD model performs independent calculations according to the time step of the fluid CFD model;

[0092] The lining - soil integrated DEM model performs independent calculations according to the time step of the lining - soil integrated DEM model;

[0093] The fluid CFD model and the lining - soil integrated DEM model perform coupled calculations according to the coupling frequency between the fluid CFD model and the lining - soil integrated DEM model. After updating the particle - fluid interaction force and particle positions, the independent calculations of the fluid CFD model and the lining - soil integrated DEM model are carried out respectively.

[0094] Further, the solution of the motion of the fluid phase in the fluid CFD model is as follows:

[0095]

[0096] where φ f is the volume fraction of the fluid; is the fluid velocity; ρ f is the fluid density; is the fluid pressure gradient; is the stress tensor gradient; g is the gravitational acceleration of the research area; is the particle - fluid interaction force.

[0097] Further, the motion control equation of the particles in the soil DEM model is as follows:

[0098]

[0099] where m pi , I pi , are the mass, moment of inertia, linear velocity, and angular velocity of the i - th particle respectively; n c is the number of particles in contact with the i - th particle; is the contact force between the i - th particle and the j - th particle; is the gravity force acting on the \(i\)-th particle; is the fluid-particle interaction force acting on the \(i\)-th particle; is the torque generated on the \(i\)-th particle due to the interaction between the \(i\)-th particle and the \(j\)-th particle.

[0100] Furthermore, the mesoscopic parameters of the concrete lining BPM model are set as follows:

[0101] F i = F n n i + F s t i ;

[0102] In the formula, F n and F s represent the normal and tangential force components respectively, n i and t i represent the plane unit vectors respectively. The calculation of the normal force is:

[0103] F n = K n U n ;

[0104] In the formula, U n is the overlap between particles, and K n is the contact normal stiffness, which can be calculated by the following formula:

[0105]

[0106] In the formula, and are the normal stiffnesses of the particles.

[0107] Furthermore, the incremental method is used to calculate the shear force. When contact is formed between particles, the initial shear force F s is 0. As each relative shear-displacement increment increases, the shear force can be obtained as follows:

[0108]

[0109] In the formula, is the initial shear force, and ΔF s is the shear action increment, and its calculation method is as follows:

[0110] ΔF s = -k s ΔU s ;

[0111] In the formula, k s is the tangential contact stiffness, and ΔU s is the tangential displacement increment, which can be calculated by the following formula:

[0112]

[0113] In the formula, and are the tangential stiffness of the particle, and the relative displacement increment ΔU i is calculated according to the time increment step and the contact velocity:

[0114] ΔU i =V i Δt;

[0115] In the formula, Δt is the time step increment, and V i is the contact velocity, and its calculation method is as follows:

[0116]

[0117] In the formula, x i and ω j are the translational velocity and rotational velocity of the particle respectively, and e ijk is an operation symbol, and the relative shear-displacement increment vector can be obtained:

[0118]

[0119] where n j , V j is the normal unit vector.

[0120] Furthermore, when U n ≤0 (there is a gap between the particles), at this time, the normal and tangential forces of the particles are both 0, otherwise, the slip between the particles will be adjusted by calculating the friction coefficient. The value of the friction coefficient between the particles refers to the following formula:

[0121] μ = min(μ A , μ B );

[0122] In the formula, μ A , μ B are the friction coefficients of particle A and particle B respectively.

[0123] Furthermore, the BPM model, which is used to simulate the crack initiation and propagation of concrete lining, defines the following parameters for the bonds in the BPM model: tensile strength σ t , shear strength τ s , normal stiffness per unit area and tangential stiffness and the truncation distance multiplier of the particle radius to particle radius

[0124] In the formula, R A,R B They are the radii of particle A and particle B respectively, is the minimum distance of the bond. The parallel bond between particles is similar to bonding with a brittle material between two particles. It can transfer the force and moment between particles and can be imagined as a bridge connecting two particles with a length of L.

[0125] The total force and bending moment of the parallel bond are represented by and , indicating the action of the bond on the particle, and decomposed into normal and tangential components as:

[0126]

[0127] In the formula, are the force and bending moment in the normal and tangential directions respectively, and n i and t i are the unit vectors of the contact surface respectively. When the parallel bond is formed, the initial and are 0 and increase with the increase of displacement and rotation increments. Their calculation methods are as follows:

[0128]

[0129] In the formula, ΔU n , ΔU s , Δθ n , Δθ s represent the normal and tangential displacement increments and the rotation increment respectively. A, I, and J are the cross-sectional area, moment of inertia, and polar moment of inertia of the parallel bond cross-section respectively. Their calculations are as follows:

[0130]

[0131] Therefore, the maximum tensile stress and shear stress acting on the parallel bond can be obtained:

[0132]

[0133] When the maximum tensile stress exceeds the tensile strength or the maximum shear stress is greater than the shear strength , bond fracture will occur, and at the same time, the force and moment acting on the particle by the bond will disappear, and thus the crack propagation during the damage process of concrete can be simulated.

[0134] In a specific embodiment, a method for simulating the catastrophe of an underwater shield tunnel based on CFD-DEM includes:

[0135] S1: Use 3D modeling software to establish the geometric model of the underwater shield tunnel, and delimit the studied soil range, tunnel burial depth, lining, seawater, groundwater, and air regions by setting the size parameters of the geometric model. Name each boundary of the geometric model for subsequent establishment of CFD and DEM models.

[0136] S2: Use CFD mesh generation software to mesh the geometric model of the underwater shield tunnel. The overall model considers using the non - analytic (Unresolved - resolved) CFD - DEM model, and densify the mesh in the area near the lining within the preset distance threshold, focusing on the damage development behavior of the concrete lining. Use sparse meshes for areas outside the densification process.

[0137] S3: Set the basic parameters of the CFD model, including initial boundary conditions, initial fields, basic fluid properties, solution control parameters, etc. Use the VOF method to define multiple fluids, including seawater - air, seawater - air - groundwater. Use the setfields function to set the initial regions, including the initial soil region of groundwater, initial seawater region and depth, and initial air region. Set the inlet of overlying seawater and the outlet of air, and set the pressure difference between the inlet and outlet in the initial boundary conditions to achieve the effect of overlying seawater pressure, thus constructing the fluid CFD model of the underwater shield tunnel.

[0138] S4: Set the calculation parameters of the fluid CFD model and the soil DEM model respectively, including: the time step of the fluid CFD model, the time step of the soil DEM model, and the coupling frequency between the fluid CFD model and the soil DEM model. The fluid CFD model performs independent calculations according to the time step of the fluid CFD model, the soil DEM model performs independent calculations according to the time step of the soil DEM model, and the fluid CFD model and the soil DEM model perform coupling calculations according to the coupling frequency between the fluid CFD model and the soil DEM model. After updating the particle - fluid interaction force and particle positions, then perform independent calculations of the fluid CFD model and the soil DEM model respectively until the set end time is reached to obtain the simulation results of the underwater shield tunnel catastrophe. The time step and coupling gap have a great influence on the simulation results. Through the simulation results, the entire process of the underwater shield tunnel catastrophe can be clearly observed, and the distinction between broken and unbroken Bond bonds in the BPM model of the final concrete lining can be made, the damage occurrence area of the lining can be observed, and at the same time, other data of the soil or lining can be extracted.

[0139] In the present specification, the various embodiments are described in a progressive manner. Each embodiment focuses on the differences from other embodiments. For the same or similar parts among the various embodiments, reference can be made to each other. For the apparatuses disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple. For the relevant parts, reference can be made to the description in the method section.

[0140] The foregoing description of the disclosed embodiments enables those skilled in the art to practice or use the present invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for simulating disasters of underwater shield tunnels based on CFD-DEM, characterized in that Including: Define the fluid using the VOF method and construct a fluid CFD model of the underwater shield tunnel; Establish a BPM contact model for the concrete lining model to simulate the mesoscopic development of the tunnel structure damage under seawater pressure and seismic action; Define the concrete lining model as a porous medium model, conduct tests on the porosity and permeability coefficient of the concrete standard specimens, and then calculate the permeability coefficient and porosity of the porous medium model based on the porosity and permeability coefficient obtained from the tests; Adopt the cohesive force Hertz contact to establish a soil DEM model based on soil particles to simulate the interaction between the soil and the structure and the ground settlement under seawater pressure and seismic action; Conduct uniaxial compression tests on concrete, triaxial compression tests on soil, and corresponding DEM numerical simulations according to the BPM contact model for the concrete lining model, compare the numerical simulation results with the test results, and conduct particle mesoscopic calibration; Construct an overall DEM model considering the interaction between the soil and the structure according to the particle mesoscopic calibration results; Set the time step of the fluid CFD model and the overall DEM model and the coupling frequency between the fluid CFD model and the overall DEM model, conduct coupling calculations, and obtain the catastrophe simulation results of the underwater shield tunnel.

2. The underwater shield tunnel disaster simulation method based on CFD-DEM according to claim 1, wherein, The calculating the permeability coefficient and porosity of the porous medium model based on the porosity and permeability coefficient obtained from the tests includes: Calculate the Darcy seepage coefficient and viscous resistance coefficient of the porous medium model according to the existing empirical formula and the formula in OpenFOAM: Where, S1 represents the method for calculating the Darcy seepage coefficient and viscous resistance coefficient in the empirical formula, S2 represents the method for calculating the Darcy seepage coefficient and viscous resistance coefficient when defining the porous medium model in OpenFOAM, μ represents the dynamic viscosity of the fluid, φ is the porosity, K is the permeability coefficient, ρ is the fluid density, F is the Forchheimer coefficient, u is the flow velocity, d is the Darcy coefficient, and f is the Forchheimer coefficient.

3. A CFD-DEM-based underwater shield tunnel disaster simulation method according to claim 1, characterized in that Also including: Conduct CFD mesh division on the fluid CFD model of the underwater shield tunnel, where the fluid CFD model of the underwater shield tunnel is divided into CFD meshes according to the tunnel structure area, soil area, and the air area inside the tunnel.

4. A method for simulating underwater shield tunnel disasters based on CFD-DEM according to claim 1, characterized in that, The defining the fluid using the VOF method and constructing a fluid CFD model of the underwater shield tunnel includes: Define multiple fluids using the VOF method, simulate the overlying seawater pressure by setting the pressure fluid pressure boundary, and use the continuous injection of the fluid to simulate the infinite replenishment of seawater; The calculation method of the seawater pressure includes: p = ρgh, where ρ is the density of seawater, h is the depth of the overlying seawater, and g is the acceleration of gravity; Set the pressure difference between the seawater inlet and the air outlet under the initial area and initial boundary conditions to meet the pressure control, and construct a fluid CFD model of the underwater shield tunnel.

5. A method for simulating the catastrophe of an underwater shield tunnel based on CFD-DEM according to claim 1, characterized in that, Solve the motion of the fluid phase in the fluid CFD model: where, φ f is the volume fraction of the fluid; is the fluid velocity; ρ f is the fluid density; is the fluid pressure gradient; is the stress tensor gradient; g is the gravitational acceleration of the research area; is the particle-fluid interaction force, is the gradient operator.

6. A method for simulating disasters of an underwater shield tunnel based on CFD-DEM according to claim 1, characterized in that The mesoscopic parameters of the BPM contact model for the concrete lining model include: F i = F n n i + F s t i ; where, F i is the interaction force between particles, F n and F s represent the normal and tangential force components respectively, n i and t i represent the plane unit vectors respectively, and the calculation of the normal force is: F n = K n U n ; where U n is the overlap amount between particles, and K n is the contact normal stiffness, which can be calculated by the following formula: In the formula, and are the normal stiffnesses of the particles.

7. A CFD-DEM-based underwater shield tunnel disaster simulation method according to claim 1, characterized in that The motion control equations of the particles in the soil DEM model include: where m pi , I pi , are the mass, moment of inertia, linear velocity, and angular velocity of the i-th particle, respectively; n c is the number of particles in contact with the i-th particle; is the contact force between the i-th particle and the j-th particle; is the gravitational force acting on the i-th particle; is the fluid-particle interaction force acting on the i-th particle; is the torque generated on the i-th particle by the interaction between the i-th particle and the j-th particle.

8. A CFD-DEM-based underwater shield tunnel disaster simulation method according to claim 1, characterized in that Perform uniaxial compression tests on concrete, triaxial compression tests on soil, and corresponding DEM numerical simulations, and compare the numerical simulation results with the test results. The particle mesoscopic calibration includes: Determine the particle size distribution of the original soil sample, conduct triaxial compression tests on the soil, and obtain the stress-strain relationship curve, Mohr stress circle, and strength envelope of the original soil sample; Conduct DEM numerical simulations of triaxial compression tests on soil by the discrete element method; Use C55 cement and multi-graded aggregates to cast a cube with a side length of 150 mm for uniaxial compression tests, and obtain the failure mode and stress-strain curve of the concrete specimen; Use the discrete element method to establish a corresponding numerical model to carry out uniaxial compression simulation and compare it with the test results; Establish an overall DEM model and calibrate the particle parameters required for the overall DEM model, the mesoscopic parameters of the concrete lining BPM model, and the fracture parameters of the corresponding bond; 9. A method for simulating underwater shield tunnel disasters based on CFD-DEM according to claim 1, characterized in that The construction of the overall DEM model considering the interaction between soil and structure according to the particle mesoscopic calibration results further includes: Densify the grid in the concrete lining model area of the BPM contact model, and set the fluid grid size to 2-3 times that of the particles; In the non-seepage area, select particles with a diameter 2-3 times the average diameter of the particles in the area near the tunnel; 10. A method for simulating disasters of underwater shield tunnels based on CFD-DEM according to claim 1, characterized in that The coupled calculation includes: The fluid CFD model performs independent calculations according to the time step of the fluid CFD model; The overall DEM model performs independent calculations according to the time step of the overall DEM model; The fluid CFD model and the overall DEM model perform coupled calculations according to the coupling frequency between the fluid CFD model and the overall DEM model. After updating the particle-fluid interaction force and particle position, the independent calculations of the fluid CFD model and the overall DEM model are carried out respectively.

Images