Seismic liquefaction fluid-solid coupling numerical value method for stratum containing tunnel saturated sand

Abstract

The invention relates to the crossing field of geotechnical engineering and seismic engineering, and provides a seismic liquefaction fluid-solid coupling numerical method for a stratum containing tunnel saturated sand. Comprising the following steps: step 0, establishing particle models of a model box, a stratum and a tunnel in PFC3D; the method comprises the steps that 1, fluid grids are divided in OpenFOAM, and initial calculation is conducted in combination with initial porosity information provided by PFC3D; 2, the initial fluid speed, pressure and pressure gradient obtained through initialization calculation are transmitted to a PFC3D; 3, the PFC3D applies a seismic oscillation acceleration time history to the bottom of the model box, and dynamic calculation is carried out; through speed correction and coordinate system transformation, the problem that computational domains of a particle model and a fluid model cannot coincide under the condition of dynamic change of the computational domains in a vibration table test simulation process is solved, so that application of the CFD-DEM method in seismic research is realized, and fine simulation of the whole saturated sand seismic liquefaction process in the presence of a tunnel can be realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the cross field of geotechnical engineering and earthquake engineering, and in particular to a numerical method for fluid-structure coupling of seismic liquefaction of saturated sand stratum containing a tunnel. BACKGROUND

[0002] Seismic liquefaction refers to the phenomenon that the pore water pressure in saturated sand soil rapidly rises, the effective stress between soil bodies decreases or even becomes zero, resulting in loss of soil strength and stiffness. In the eastern coastal areas, underground structures such as tunnels are usually built in saturated sand stratum, and seismic liquefaction may cause the tunnel to float, deviate or be damaged.

[0003] Traditional liquefaction analysis methods mainly include laboratory tests and numerical simulation methods. The laboratory tests mainly include shaking table tests, which are costly and have single observation means, and are difficult to meet all the requirements of the research. The numerical simulation method is based on the constitutive model and the continuum theory and is carried out by the finite element method, which is limited by the calculation principle and is difficult to fully reveal the nonlinear large deformation process of fluid-structure interaction at the particle scale, and it is even more difficult to reveal the micro interaction between soil particles and fluid.

[0004] In recent years, the coupling calculation of discrete element method (DEM) and finite volume method (FVM) has provided a new idea for the research of fluid-structure coupling problems. This method carries out solid particle calculation through discrete element method (DEM) and fluid calculation through finite volume method (FVM), and realizes coupling calculation through information interaction between the two. The common discrete element method (DEM) software is PFC3D, and the common finite volume method (FVM) software is OpenFOAM. However, due to the fact that the position of the calculation domain in the fluid calculation of OpenFOAM cannot be moved, the fluid calculation domain and the solid particle calculation domain are inconsistent, so it is rarely used for dynamic calculation, which limits its application in the simulation of shaking table test of stratum liquefaction response. The Chinese invention patent application with publication number CN115544911A introduces a karst area water and land surface erosion / underground leakage process simulation method based on CFD-DEM fluid-structure coupling calculation. However, this method requires that the particle model area and the fluid area coincide, and the calculation range does not change, which is not suitable for simulating the fluid-structure coupling of seismic liquefaction of saturated sand stratum containing a tunnel under the condition of a shaking table test.

[0005] Therefore, the previous numerical method using fluid-structure (particle) coupling is not suitable for simulating the working condition under seismic conditions, and the simulation results have large differences compared with the test data. SUMMARY

[0006] To address the problems of existing technologies, the purpose of this invention is to propose an efficient, scalable, and high-precision numerical method that overcomes the limitations of traditional fluid-structure interaction (FSI) calculations, where the fluid computational grid is fixed and cannot move synchronously with solid particles. This enables the simulation of shaking table tests, reduces experimental costs, and provides a method for the macro- and micro-level study of seismic liquefaction processes in tunnel-saturated sandy soil strata. It also reveals the influence mechanism of tunnels on the liquefaction development and pore pressure evolution of the surrounding foundation under seismic loading. This invention improves upon traditional numerical methods through in-depth research into the fluid-structure interaction mechanism, ensuring its application under seismic conditions.

[0007] Technical solution A numerical method for seismic liquefaction fluid-structure interaction in saturated sandy soil strata containing tunnels, wherein the fluid dynamics part is implemented based on the open-source software OpenFOAM, and the particle stress and calculation are implemented based on the commercial software PFC3D, including the following steps: Step 0: In PFC3D, create particle models of the model box, strata, and tunnel, set boundary conditions, assign corresponding strata and tunnel parameters, set the calculation time step, and perform initialization calculations; after initialization, provide the initial porosity information to OpenFOAM. Step 1: Generate a fluid mesh in OpenFOAM and assign fluid density. and viscosity Set the calculation step size and perform initial calculations based on the initial porosity information provided by PFC3D; Step 2: Initialize the calculated initial fluid velocity. ,pressure and pressure gradient ( x (For distance) transmitted to PFC3D; Step 3: PFC3D receives fluid information, starts the PFC3D fluid dynamics calculation module to apply it to the particle model, and applies the seismic acceleration time history to the bottom of the model box to perform dynamic calculation; Step 4: Dynamically update the position of the particle model in PFC3D so that it completely overlaps with the computational domain of OpenFOAM; Step 5: Update particle and flow field velocity information in PFC3D; Step 6: Update the porosity and fluid-structure interaction of PFC3D and transmit them to OpenFOAM for flow field information update; Step 7: OpenFOAM performs fluid calculations to obtain the corrected fluid velocity, pressure, and pressure gradient information; Step 8: The OpenFOAM field calculation and processing software transmits fluid velocity, pressure, and pressure gradient information to PFC3D; Step 9: update the flow field information in the PFC3D particle flow program module, jump to step 3 to repeat the fluid-structure coupling calculation until the iteration reaches the end of the input time history curve, that is, stop the calculation.

[0008] Advantages The application realizes the fine simulation of the whole process of saturated sand liquefaction under the condition of the tunnel by introducing the tunnel structure under the CFD-DEM framework, and has the following remarkable technical effects: (1) The application solves the problem that the particle model and the fluid model calculation domain cannot coincide in the calculation domain dynamic change in the simulation process of the shaking table test through velocity correction and coordinate system transformation, thereby realizing the application of the CFD-DEM method in seismic research.

[0009] (2) The method can realize the two-way coupling dynamic response process between the tunnel and the saturated foundation, and can accurately capture the whole process of the rise of pore water pressure, the decay of effective stress and the occurrence of liquefaction under the action of seismic loading; (3) The influence law of the tunnel structure on the liquefaction mode of the surrounding stratum is clear, and the spatial difference of the liquefaction development degree in the soil above and on both sides can be identified; (4) A numerical platform consistent with the shaking table model test results is provided, and theoretical support is provided for complex underground structure seismic analysis; (5) The method has good expansibility and repeatability, and can be used for liquefaction response simulation under different seismic motion, stratum structure and structure size conditions.

[0010] In summary, the application provides a numerical analysis method with innovation and practical value for tunnel seismic resistance and foundation liquefaction prevention and control, which can extend the seismic liquefaction research from macro soil changes to micro fluid-solid interaction, and reduce the dependence on expensive shaking table tests. BRIEF DESCRIPTION OF DRAWINGS

[0011] Fig. 1 is a schematic diagram of the processing steps of the method of the application; Fig. 2 is a schematic diagram of a discrete element model containing a tunnel according to an embodiment of the application; Fig. 3 is a schematic diagram of a flow field model containing a tunnel according to an embodiment of the application; Fig. 4 is a schematic diagram of a fluid-structure coupling model according to an embodiment of the application; Fig. 5 is a graph of the input wave time history curve of the numerical simulation according to an embodiment of the application; Fig. 6 is a schematic diagram of an effective stress monitoring point according to an embodiment of the application; Fig. 7 is a comparison diagram of effective stress in simulation and test according to an embodiment of the application ((a) stratum height 0.2m, near field (b) stratum height 0.6m, near field). DETAILED DESCRIPTION

[0012] The technical solution provided in this application will be further described below with reference to specific embodiments and accompanying drawings. The advantages and features of this application will become clearer from the following description.

[0013] A numerical method for seismic liquefaction fluid-structure interaction in saturated sandy soil strata containing tunnels, with the fluid dynamics part implemented based on the open-source software OpenFOAM, and particle stress and calculation implemented based on the commercial software PFC3D, includes the following steps: (e.g.) Figure 1 ) Step 0: In PFC3D, create particle models of the model box, strata, and tunnel, set boundary conditions, assign corresponding strata and tunnel parameters, set the calculation time step, and perform initialization calculations; after initialization, provide the initial porosity information to OpenFOAM. Step 1: Generate a fluid mesh in OpenFOAM and assign fluid density. and viscosity Set the calculation step size and perform initial calculations based on the initial porosity information provided by PFC3D; Step 2: Initialize the calculated initial fluid velocity. ,pressure and pressure gradient ( x (For distance) transmitted to PFC3D; Step 3: PFC3D receives fluid information, starts the PFC3D fluid dynamics calculation module (CFD module) to apply it to the particle model, and applies the seismic acceleration time history to the bottom of the model box to perform dynamic calculation. Specifically, the earthquake ground acceleration time history is a record of the seismic wave acceleration over time obtained by the earthquake monitoring station at the time of the earthquake. It is either actually measured or selected from the standard. The model box is loaded by applying the acceleration at the corresponding time. The loading time is obtained by accumulating the calculated time step.

[0014] The particle model satisfies Newton's second law, and its translational velocity... and rotational speed The governing equations are: (1) (2) In the formula, It is a particle i quality It is its moment of inertia. It is a particle i and granules j Contact force between them It is a particle i and granules kNon-contact forces between particles (determined according to the specific research object whether to add) and are fluid-structure interaction forces and body forces (such as gravity), respectively, and are the tangential force and rolling friction torque generated between the particles i and the particles j .

[0015] The fluid-structure interaction force mentioned above is specifically: (1) Drag force: (3) In the formula, is the fluid density, is the particle diameter, and are the velocities of the fluid and the particle, respectively, is the porosity. is the drag coefficient, which is determined according to formula 4.

[0016] (4) is the Reynolds number, which is determined according to formula 5.

[0017] (5) In the formula, is the fluid viscosity.

[0018] (2) Buoyancy force: (6) In the formula, V p is the particle volume.

[0019] (3) Pressure gradient force: (7) Step 4: Dynamically update the position of the particle model in PFC3D to completely coincide with the calculation domain of OpenFOAM. Under the action of seismic force, the PFC3D model box will displace, and since the grid coordinates of the OpenFOAM fluid model do not change, it is necessary to dynamically adjust the position of the particle model in PFC3D to correct the particle model coordinates to the initial position, so that it completely coincides with the calculation domain of OpenFOAM.

[0020] Specifically, first record the 8 corner point coordinates of the model box, under dynamic loading, after each time step calculation is completed, calculate the displacement difference between the new corner point coordinates and the original coordinates of the model box, use the displacement difference value to correct the coordinates of each particle, tunnel and model box in the entire PFC3D particle model, so that it is completely coincided with the fluid model.

[0021] Step 5: Update the particle and flow field velocity information in PFC3D (8) (9) Where, v p-b is the particle velocity relative to the model box, v f-b is the fluid velocity relative to the model box; v p is the particle velocity relative to the earth coordinate system; v f is the fluid velocity relative to the earth coordinate system; v b is the velocity of the model box relative to the earth coordinate system.

[0022] It should be noted that in PFC3D software, only the particle velocity needs to be corrected, and the fluid velocity is automatically calculated after the particle velocity is corrected.

[0023] As can be seen from formulas (3)~(7), the interaction force between fluid and solid particles is affected by the velocity difference between the two, so the velocities of the two in the same reference system need to be unified. Therefore, in PFC, the particle velocity is corrected before the particle velocity information is transmitted to OpenFOAM, that is, the particle velocity transmitted to the fluid should be v p-b , and at this time the velocity calculated by the fluid grid should be v f-b . The reason is as follows: in actual test, the particle velocity relative to the earth coordinate system v p can be regarded as the sum of the model box velocity v b and the particle velocity relative to the model box v p-b (formula 8), and similarly, the fluid velocity relative to the earth coordinate system v f can also be regarded as the sum of the model box velocity v b and the fluid velocity relative to the model box v f-b (formula 9).

[0024] Step 6: Update the porosity of PFC3D and the fluid-solid interaction force, and transmit it to OpenFOAM for flow field information update; Step 7: OpenFOAM performs fluid calculation and solving to obtain the corrected fluid velocity, pressure and pressure gradient information.

[0025] Specifically, the OpenFOAM fluid control equation is the average volume Navier-Stokes equation, and the fluid velocity of the Navier-Stokes equation is corrected by formula (9), and formula (10)~formula (11) can be obtained: (10) (11) Among them: (12) (13) In formula (10)~formula (13), And p The corrected fluid velocity and pressure are respectively, The fluid shear stress stress tensor of the calculation unit, The porosity of the calculation unit, g is the acceleration of gravity, F A The particle-fluid interaction force, V cell The volume of the fluid unit, The total fluid-particle interaction force acting on the first i Particle, The sum of other fluid-particle interaction forces.

[0026] Step 8: The OpenFOAM field operation and processing software transmits the fluid velocity, pressure and pressure gradient information to PFC3D; Step 9: Update the flow field information in the PFC3D particle flow program module, jump to step 3 and repeat the fluid-structure coupling calculation until the iteration reaches the end of the input time curve, that is, stop the calculation.

[0027] In summary, in the present application: Through steps 0 and 1: through OpenFOAM and PFC3D coupling modeling, three-dimensional fluid and particle domain synchronous initialization is realized, and the fluid-structure coupling calculation basis is constructed.

[0028] Through step 3: apply actual seismic motion to ensure effective excitation, introduce key forces such as drag force, buoyancy and gradient force in DEM, drive particle motion and respond to fluid pressure change, which is the direct dynamic mechanism of liquefaction.

[0029] Through steps 4~5: update the model calculation domain and use reference frame transformation to ensure that the fluid-structure motion is consistent, and the particle feedback information is transmitted to the fluid domain.

[0030] Through steps 6-9: the fluid response is calculated through the Navier-Stokes equation, and bidirectional data exchange is realized with PFC3D, the particle behavior is updated in real time, and the spatial influence of the structure on the development of liquefaction is reflected.

[0031] Through the technical solution, the application can reproduce the whole process of the rapid growth of pore pressure, the disappearance of effective stress, and the non-uniform development of liquefaction under the influence of the structure in the saturated sand soil around the tunnel in the shaking table test under the seismic excitation, and has high accuracy and applicability. Embodiment Taking the saturated sand soil layer containing the tunnel as the object, the liquefaction response process of the saturated sand soil layer under the action of the Kobe seismic wave is simulated, and the rationality and accuracy of the numerical model are verified. Specifically as follows: (1) Discrete element model construction (corresponding to step 0) A PFC3D particle model with a size of 0.87 m (high) x 0.8 m (long) x 0.3 m (wide) is established (as shown in FIG. 2, wherein different colors represent different sizes of particles), wherein the sand soil adopts spherical particles with a radius distribution of 0.3-0.8 mm, a total of about 250,000 particles, and forms a dense sand soil state with an initial porosity of 0.38. The tunnel adopts a hollow cylindrical structure with a diameter of 0.1 m, which is constructed by an impermeable Wall unit. The buried depth is 0.47 m, and the tunnel axis is parallel to the vibration loading direction. The stiffness value of the tunnel material is 5 x 10 9 Pa, which is regarded as a rigid structure participating in the fluid-structure coupling process. The time step of the PFC3D software is related to the particle parameters and the fluid parameters, and is specifically determined by trial calculation. In this embodiment, the DEM time step is set to: 1 x 10 -6 s.

[0032] (2) Fluid domain establishment and meshing (steps 1-2) A structured fluid mesh is generated in OpenFOAM, which is one-to-one mapped with the particle domain. The mesh size is 80 x 30 x 90, and the total number of meshes is about 210,000. The area near the tunnel is treated by mesh densification, and the minimum unit edge length is about 1.2 mm. The CFD step is: 5 x 10 -5 s; the coupling step frequency is: 50 steps of DEM correspond to 1 step of CFD. The initial porosity field is generated by mapping the particle volume fraction of PFC and is transmitted to the initial field of OpenFOAM. The obtained flow field is shown in FIG. 3, and the flow field information is transmitted to PFC3D after initialization.

[0033] (3) Discrete element model calculation-formation settlement (step 3) After receiving the initial flow field information, the particle model is first subjected to settlement calculation, and the fluid-structure coupling model after settlement is shown in FIG. 4.

[0034] (4) Discrete element model calculation - dynamic loading (step 3) After obtaining the settlement model, a typical component EW (east-west direction) in the Kobe earthquake wave is input, with a maximum peak acceleration of 0.21 g, a sampling interval of 0.005 s, and a total duration of 10 s. The loading method is to assign an acceleration to the model box to force displacement, simulate the physical shaking table surface movement, and the input wave time curve is as shown in Figure 5 .

[0035] (5) Correcting position and velocity information (step 4-5) After each time step loading is completed, the particle model is corrected for position and velocity.

[0036] (6) Fluid-structure interaction information calculation and transmission (step 6) Update the porosity and fluid-structure interaction force and transmit to OpenFOAM.

[0037] (7) Fluid model calculation and information interaction (step 7-9) After receiving the information transmitted by the PFC3D model, the flow field information is updated and calculated. In this embodiment, the fluid viscosity resistance model uses the Ergun correction formula, and the pressure-velocity coupling algorithm uses the PISO method. A velocity smoothing mechanism is used to filter the particle instantaneous velocity by 3 orders to avoid instantaneous fluctuations affecting the stability of the flow field. After the calculation is completed, the new fluid velocity, pressure, pressure gradient, etc. information is transmitted to PFC3D. After receiving, PFC3D updates the relevant information and performs the next time step calculation until the required seismic wave loading is completed.

[0038] (8) Monitoring and data acquisition (data post-processing) Output the global particle and inter-particle contact information, and set 2 groups of monitoring points (0.2 m and 0.6 m from the bottom of the model, respectively, S1 and S2, Figure 6 as shown in FIG. 6) for subsequent calculation of effective stress. Data is output every 0.01 s, with a total of 1000 groups of data points.

[0039] (9) Simulation results The simulation and test data are normalized to obtain the stress time curve at the specified detection position as shown in FIG. 7.

[0040] The results show that, as Figure 7 (a), at a stratum height of 0.2 m, the vertical effective stress at the near-field position (i.e., close to the tunnel) changes little during the entire seismic loading process, and only a slight decrease occurs in the late loading period. The numerical simulation results are basically consistent with the overall trend of the test data.

[0041] As Figure 7(b), at the stratum height of 0.6 m, the vertical effective stress changes significantly at the end of loading, all showing the evolution characteristics of first decreasing, then increasing, and then rapidly reducing, the simulation results are consistent with the test results, reflecting the typical characteristics of the sudden drop of effective stress in the process of liquefaction.

[0042] The above description is only a description of the preferred embodiments of the present application, and is not any limitation on the scope of the present application. Any modification or modification made by any ordinary skilled person in the art according to the above disclosed technical content should be regarded as an equivalent effective embodiment, and belongs to the protection scope of the technical scheme of the present application.

[0043] Appendix: Main Noun Explanation Table 1 Noun Explanation Table

Claims

1. A numerical method for fluid-structure interaction of seismic liquefaction in saturated sand layer with tunnels, the fluid dynamics part is based on open source software OpenFOAM, and the particle stress and calculation are based on commercial software PFC3D, characterized in that, Comprising the following steps: Step 0: Establish the particle model of the model box, stratum and tunnel in PFC3D, set the boundary conditions, assign the corresponding stratum and tunnel parameters, set the calculation time step and perform the initialization calculation; after the initialization is completed, provide the initial porosity information to OpenFOAM; Step 1: Divide the fluid grid in OpenFOAM, give the fluid density and viscosity , set the calculation step, and initialize the calculation combined with the initial porosity information provided by PFC3D; Step 2: The initial fluid velocities calculated from the initialization are transferred to PFC3D, where pressure and pressure gradients are calculated. x is the distance; Step 3: PFC3D receives the fluid information, starts the PFC3D fluid mechanics calculation module to act on the particle model, and applies the seismic acceleration time history at the bottom of the model box to perform dynamic calculation; Step 4: Dynamically update the position of the particle model in PFC3D to completely coincide with the calculation domain of OpenFOAM; Step 5: Update the particle and flow field velocity information in PFC3D; Step 6: Update the porosity of PFC3D and the fluid-structure interaction force, and transmit it to OpenFOAM for flow field information update; Step 7: OpenFOAM performs fluid calculation and solving to obtain the corrected fluid velocity, pressure and pressure gradient information; Step 8: The OpenFOAM field operation and processing software transmits the fluid velocity, pressure and pressure gradient information to PFC3D; Step 9: Update the flow field information in the PFC3D particle flow program module, jump to step 3 to repeat the fluid-structure coupling calculation until the iteration reaches the end of the input time history curve, that is, stop the calculation.

2. The numerical method for fluid-structure interaction of seismic liquefaction in the tunneling saturated sand stratum according to claim 1, wherein, In step 3, the seismic acceleration time history is the record of the seismic wave acceleration about time obtained by the seismic monitoring station when the earthquake occurs, which is actually measured or selected from the specification, and is loaded by applying the acceleration at the corresponding time to the model box. The loading time is obtained by accumulating the calculation time step.

3. The numerical method for fluid-structure interaction of seismic liquefaction in the tunneling saturated sand stratum according to claim 1, wherein, In Step 3, the particle model satisfies Newton's second law, and the control equation of its translational velocity and rotational velocity is （1） （2） wherein is the mass of the particle i , is the moment of inertia of the particle , i is the contact force between the particle j and the particle , i is the non-contact force between the particle k and the particle , and are the fluid-solid interaction force and the body force, respectively, and i are the tangential moment and the rolling friction moment generated between the particle j and the particle The above fluid-structure interaction force is specifically: (1) Drag force: （3） wherein is the fluid density, is the particle diameter, and are the velocities of the fluid and the particle, respectively, is the porosity; is the drag coefficient, determined according to equation 4; （4） Reynolds number, determined according to equation 5; （5） wherein is the fluid viscosity; (2) Buoyancy: （6） wherein V p Vp is the volume of the particles; (3) Gradient force: （7）。 4. The numerical method for fluid-structure interaction of seismic liquefaction in a saturated sand layer with a tunnel according to claim 1, wherein, In step 4, the particle model position adjustment method is as follows: First, record the coordinates of the 8 corner points of the model box. After each time step calculation under dynamic loading, calculate the displacement difference between the new corner point coordinates of the model box and the original coordinates, and use the displacement difference value to correct the coordinates of each particle, tunnel and model box in the entire particle model of PFC3D to completely coincide with the fluid model.

5. The numerical method for fluid-structure interaction of seismic liquefaction in a saturated sand layer with a tunnel according to claim 1, wherein, In step 5, the particle and flow field velocity information updating method is as follows: （8） （9） where v p-b is the velocity of the fluid relative to the model tank; v f-b is the velocity of the fluid relative to the model tank; v p is the velocity of the particle relative to the earth coordinate system; v f is the velocity of the fluid relative to the earth coordinate system; v b is the velocity of the model tank relative to the earth coordinate system.

6. The numerical method for fluid-structure interaction of seismic liquefaction in a saturated sand layer with a tunnel according to claim 1, wherein, In step 7, specifically: The OpenFOAM fluid control equation is the average volume Navier-Stokes equation, which is corrected by formula (9) to obtain formula (10)~formula (11): （10） （11） Wherein: （12） （13） in Equations (10) to (13), and p are the modified fluid velocity and pressure, respectively, is the fluid shear stress stress tensor of the calculation cell, is the porosity of the calculation cell, g is the gravitational acceleration, F A is the inter-particle fluid interaction force, V cell is the volume of the fluid cell, is the total fluid-particle interaction force acting on the i th particle, is the sum of other fluid-particle interaction forces.

Images