A Method for Dynamic Response Simulation and Damage Classification of Underground Structures under Explosive Load Based on ALE Fluid-Structure Coupling

Abstract

This invention relates to the field of underground engineering protection and computational mechanics simulation technology, and provides a method for simulating the dynamic response and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction. The method includes the following steps: Step 1, unified domain modeling and optimization of working parameters; Step 2, selection of material model and equation of state; Step 3, application of initial field, boundary, and initiation conditions; Step 4, explicit dynamic ALE remapping solution; Step 5, result extraction and multi-index fusion classification; Step 6, multi-condition batch processing and sensitivity analysis; Step 7, engineering calibration of parameters and thresholds. This scheme establishes an integrated ALE fluid-structure interaction calculation framework encompassing the "explosive gas-surrounding rock-underground structure," constructing a robust, efficient, and quantifiable three-dimensional ALE fluid-structure interaction simulation and damage classification technology system that is engineering-applicable, robust, efficient, and quantifiable.

Description

Technical Field

[0001] This invention relates to the field of underground engineering protection and computational mechanics simulation technology, specifically to a method for simulating the dynamic response of underground structures to explosive loads and classifying damage based on ALE fluid-structure interaction. Background Technology

[0002] Underground structures (such as subway tunnels, integrated utility tunnels, and underground stations) exhibit strong transient, strong nonlinear, and multi-domain coupled dynamic characteristics under near-field blast loads. The blast shock wave penetrates the overburden and surrounding rock in an extremely short time, reflecting, superimposing, and diffracting with the structure, gradually transforming into a stress wave. During this process, significant inertial effects, plastic deformation, and damage accumulation occur between the structure and the surrounding rock. Common failure modes include circumferential cracking of the lining, local crushing, interface slippage, and instability at critical nodes. These failures exhibit significant spatiotemporal nonhomogeneity and suddenness, making accurate assessment difficult using traditional steady-state or quasi-static analysis methods.

[0003] Current blast-resistant design and evaluation in engineering mainly rely on three types of methods: simplified equivalent methods (such as equivalent static load methods) are easy to calculate but rely on experience and are difficult to reflect the propagation and interaction of real blast waves; single-domain numerical methods (pure Lagrangian or pure Eulerian methods) have their own limitations in characterizing large fluid deformation or solid damage details and are difficult to handle fluid-solid interfaces; and while the partitioned coupling method improves efficiency to some extent, the spatiotemporal distribution of the applied blast load is difficult to generalize and cannot truly reproduce the cross-domain transfer and dissipation of energy between gas, surrounding rock, structure and interface.

[0004] From a materials and constitutive perspective, the nonlinear behavior of materials under high strain rates (such as surrounding rock softening, concrete rate effect and damage, and steel reinforcement bond slip) are coupled with each other, making the selection of constitutive parameters and model stability key bottlenecks. In terms of boundaries and solutions, improper handling of complex geostress, groundwater, finite domain boundaries and contact conditions can easily lead to numerical distortion and sensitivity to calculation parameters, affecting the reliability of engineering results.

[0005] From the perspective of engineering assessment needs, practical engineering projects require rapid scheme comparison and selection under multiple parameters (such as charging method, burial depth, surrounding rock conditions, reinforcement, etc.). Existing research is mostly limited to typical single working conditions, lacking systematic multi-working-condition batch processing, sensitivity analysis, and surrogate model capabilities, making it difficult to support efficient decision-making. At the same time, damage assessment usually relies on single indicators such as peak vibration velocity or stress, failing to integrate multi-dimensional information such as cracks, plastic work, and damage volume, resulting in insufficient robustness in identifying weak points. The calibration with field monitoring data also lacks a unified process, limiting the cross-engineering application of models and thresholds.

[0006] Based on the above situation, there is an urgent need to establish a unified three-dimensional fluid-structure interaction computational framework capable of simultaneously describing the dynamic interactions between explosive gases, surrounding rock, and structures, and realistically simulating shock wave propagation, energy attenuation, and damage evolution while ensuring numerical stability and energy conservation. The ALE (Arbitrary Lagrange-Eulerian) method combines the advantages of the Eulerian description for large-deformation fluids with the Lagrange description's ability to characterize solid details. However, in engineering applications, it still requires systematic and standardized solutions in areas such as mesh stability, material parameter determination, non-reflective boundary setting, multi-index fusion evaluation, and automated batch processing to effectively support the blast-resistant design and safety assessment of underground structures. Summary of the Invention

[0007] To address the problems existing in the background technology, this invention proposes a method for dynamic response simulation and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction, and constructs a three-dimensional ALE fluid-structure interaction simulation and damage classification technology system that is engineering-applicable, robust, efficient, and quantifiable.

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

[0009] A method for simulating the dynamic response and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction includes the following steps:

[0010] Step 1: Establish a three-dimensional computational domain that includes the underground structure domain, the surrounding rock domain, and the explosive gas domain, and define the geometric parameters and working parameters and set the mesh.

[0011] Step 2: Define material parameters and select corresponding constitutive models for the underground structural domain and the surrounding rock domain, and define material parameters and select equations of state for the gas domain; at the same time, set the constitutive or contact model of the underground structural-surrounding rock interface and the ALE fluid-structure coupling condition of the gas-solid interface.

[0012] Step 3: Apply initial and boundary conditions related to the explosion load and geological environment to the computational domain;

[0013] Step 4: Perform explicit dynamic integration on the computational domain, and monitor and record multiple energy parameters in real time;

[0014] Step 5: Construct a multi-damage index system and damage classification by extracting the solution results from Step 4, and formulate targeted engineering design suggestions.

[0015] Step 6: Establish a batch case generation and automated solution process, and build a proxy model;

[0016] Step 7: Based on field monitoring data or experimental data, perform inversion calibration on model parameters and grading thresholds.

[0017] Optionally, step 1 specifically includes:

[0018] Step 1.1: Construct a three-dimensional computational domain, which includes the underground structure domain and surrounding rock domain described by Lagrange, and the explosive gas domain described by Eulerian or ALE.

[0019] Step 1.2: Parametrically set the geometric parameters and working parameters. The geometric parameters include burial depth, overburden thickness and structural joint parameters. The working parameters include load application method, loading sequence, equivalent charge, in-situ stress coefficient K0, pore water pressure and saturation.

[0020] Step 1.3: Implement differentiated grid division, locally densify the key parts of the underground structural domain, and adopt a grid division method of progressively thickening from the cavern to the periphery for the surrounding rock domain, and make the characteristic size of the Eulerian domain unit corresponding to the explosive gas match the wavelength.

[0021] Step 1.4: Set contact and interface conditions. The underground structure and surrounding rock adopt a penalty function-Coulomb friction model with bond-slip degradation function. The gas-solid interface is coupled using ALE contact or volume fraction exchange algorithm.

[0022] Optionally, in step 2, the gas domain adopts the JWL or equivalent EOS equation of state for the explosion products, the surrounding rock domain adopts the constitutive model of the Drucker-Prager criterion or the Mohr-Coulomb criterion, the underground structural domain adopts the concrete damage plasticity or equivalent damage model considering tensile-compressive asymmetry and strain rate effect, and the interface between the lining and the surrounding rock adopts the bond-slip degradation model.

[0023] Optionally, in step 3, an initial triaxial in-situ stress field determined by the burial depth and the in-situ stress coefficient K0 is applied to the surrounding rock domain and / or the underground structure domain, and gravity is superimposed; for the surrounding rock domain, an initial pore water pressure field is applied according to the pore or groundwater conditions, or an equivalent hydrostatic pressure field obtained by solid-fluid coupling calculation is used as the pore pressure loading; a non-reflective boundary is set at the outer boundary of the surrounding rock domain, wherein the non-reflective boundary is a viscous boundary or an infinite element boundary; an initiation zone is defined in the gas domain, and an initial energy and pressure growth history is assigned to the initiation zone.

[0024] Optionally, in step 4, a remapping algorithm is applied to the gas domain, and the time step is controlled based on the Courland-Friedrich-Liouvi CFL stability condition; the solid domain is subjected to mass scaling within a limited range; the energy parameters include external work, kinetic energy, internal energy, viscoplastic energy dissipation, and artificial viscous dissipation.

[0025] Optionally, step 5 specifically includes:

[0026] Step 5.1: Extract the time history response parameters of key measuring points from the numerical calculation results, and output the isochronous cross section or isosurface map of the spatial field.

[0027] Step 5.2: Based on the results obtained in Step 5.1, construct a multi-damage index system including I1 peak vibration velocity PPV, I2 tensile damage volume fraction, I3 plastic work density, I4 crack continuity index and I5 stiffness reduction rate.

[0028] Step 5.3: Preset the grading thresholds and fusion weights for damage indices I1–I5, and output damage levels L0 to L4 based on weighted scoring or logic thresholds, automatically marking weak areas. The grading thresholds are determined based on standard limits, experimental data, or simulation calibration, and the fusion weights are determined based on analytic hierarchy process, entropy weighting, or data-driven methods.

[0029] Step 5.4: Overlay the damage grading results with the weak points to generate engineering design recommendations that include local reinforcement, interface enhancement, allowable equivalent, and minimum safe burial depth suggestions.

[0030] Optionally, the time history response parameters include stress, strain, equivalent plastic strain, damage factor D, particle velocity, and peak velocity PPV; the isochronous cross-sectional diagram or isosurface diagram of the spatial field includes pressure wavefront distribution, energy density distribution, and damage cloud map distribution.

[0031] Optionally, step 6 specifically includes:

[0032] Step 6.1: Generate multiple sets of working condition samples using Latin hypercube design or orthogonal design methods. The parameters covered by the working condition samples include equivalent charge, charge arrangement or detonation arrangement parameters, underground structure burial depth, K0, pore water pressure, and material parameter sets.

[0033] Step 6.2: Perform batch automatic calculation on the multiple sets of working condition samples, and extract the damage index I1–I5 under each working condition after the calculation is completed;

[0034] Step 6.3: Construct a proxy model based on the Kriging model or response surface model.

[0035] Optionally, step 7 specifically includes:

[0036] Step 7.1: Using on-site monitoring data or test data as observations, the material parameters and the interface parameters of the underground structure and surrounding rock are inverted and calibrated using the Bayesian update method or the least squares inversion method.

[0037] Step 7.2: Based on the results of the material parameters and interface parameters inversion calibration in Step 7.1, update the grading thresholds of L0–L4 corresponding to the I1–I5 damage indices in Step 5.4 and the multi-index fusion weights; and propagate the parameter uncertainty to the damage indices and grading judgment results, outputting the damage grading results with uncertainty intervals, and update and optimize the engineering design suggestions in Step 5.5 accordingly.

[0038] Optionally, the method further includes: automatically generating visual charts and textual conclusions related to numerical calculations, and realizing model version traceability and cross-site reuse by saving metadata.

[0039] The beneficial effects of this invention are as follows: This solution aims to address the engineering verification and optimization of underground structures under explosive loads, and constructs a robust, efficient, and quantifiable three-dimensional ALE fluid-structure interaction simulation and damage classification technology system that is engineering-applicable, robust, and quantifiable. Specifically, it includes:

[0040] (1) An ALE fluid-structure interaction calculation framework integrating the three domains of “explosive gas-surrounding rock-underground structure” was established. Under conditions such as single-point / multi-point detonation, different burial depths, charge equivalent and in-situ stress, the shock wave propagation, reflection / diffraction and energy transfer process were stably reproduced.

[0041] (2) It integrates the equation of state of explosive gas (EOS), the Mohr-Coulomb / DP constitutive model of surrounding rock (including softening and shear dilatation) and the damage plasticity of concrete (including tensile-compressive asymmetry and strain rate effect), and couples the bond-slip degradation model of lining-surrounding rock, which realistically depicts key mechanisms such as cracking-slip-reclosure and local crushing.

[0042] (3) A grid densification strategy, non-reflective boundary configuration, ALE remapping and time step control principles were formed, and an energy monitoring criterion of external work-internal energy-kinetic energy-viscoplastic energy consumption was established to ensure computational stability and repeatability under strong impact conditions.

[0043] (4) A comprehensive classification system (L0–L4) consisting of peak vibration velocity (PPV), tensile damage volume fraction, plastic work density, crack continuity index (CI) and stiffness reduction rate was constructed, and the weak parts such as arch waist, arch foot and joint were automatically identified by damage hot spot connectivity and energy density threshold.

[0044] (5) A batch case generation and automated solution process based on Latin hypercube / orthogonal design was established. A proxy model was constructed for sensitivity ranking, rapid pre-evaluation and scheme selection to support decision-making at the "scheme space" level.

[0045] (6) A closed-loop calibration method with field monitoring / comparison test data is provided, and Bayesian updates are performed on material parameters and classification thresholds to improve the transferability and result confidence across sites and working conditions.

[0046] Meanwhile, this method can generate reportable results such as damage level maps, reinforcement recommendations for key components, minimum safe burial depth and allowable equivalent range required for blast resistance verification. It achieves efficient transformation from complex numerical simulation to intuitive engineering conclusions, which can directly serve design optimization and risk management, and meet the actual needs of blast resistance verification and scheme optimization for underground structures. Attached Figure Description

[0047] Figure 1 This is a flowchart of the method of the present invention;

[0048] Figure 2 This is a schematic diagram of the three-dimensional computational domain and mesh division in an embodiment of the present invention;

[0049] Figure 3 This is a schematic diagram showing the initial field, boundary conditions, and interface contact / gas-solid coupling conditions in an embodiment of the present invention.

[0050] Figure 4 This is a schematic diagram illustrating the evolution of the pressure wavefront distribution and energy density distribution isosurface under typical working conditions in an embodiment of the present invention.

[0051] Figure 5 This is a schematic diagram of the spatial distribution of lining damage factor D and PPV in an embodiment of the present invention;

[0052] Figure 6 This is a schematic diagram of the damage grading criteria and thresholds in an embodiment of the present invention. Detailed Implementation

[0053] To make the present invention clearer and more understandable, the present invention will be described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the given embodiments are only one implementation method and do not represent all embodiments.

[0054] Example

[0055] Combination Figure 1 This invention provides a method for simulating the dynamic response and damage classification of underground structures under explosion loads based on ALE fluid-structure interaction. By establishing an integrated ALE fluid-structure interaction calculation framework encompassing the three domains of "explosive gas-surrounding rock-underground structure", a robust, efficient, and quantifiable three-dimensional ALE fluid-structure interaction simulation and damage classification technology system is constructed, which is engineering-applicable, robust, efficient, and quantifiable.

[0056] Specifically, the steps include the following:

[0057] Step 1, Unified Domain Modeling and Operating Parameter Optimization: Establish a 3D computational domain including the underground structure domain, surrounding rock domain, and explosive gas domain, and define and mesh the geometric and operating parameters accordingly. Figure 2 This step overcomes the limitations of traditional partitioned coupling or single-domain simulation in accurately characterizing the cross-medium propagation, reflection, and energy exchange of shock waves, laying the foundation for dynamic response analysis. Specifically, it includes:

[0058] Step 1.1: Construct a three-dimensional computational domain, which includes an underground structure domain and a surrounding rock domain described by Lagrange, and an explosive gas domain described by Eulerian or ALE. The underground structure includes at least one of annular lining, box lining, and irregular lining.

[0059] Step 1.2 involves setting the geometric parameters and working parameters. The geometric parameters include parameters such as burial depth, overburden thickness, and structural joints. The working parameters include load application method (single-point or multi-point application), loading sequence (synchronous or asynchronous loading), equivalent charge, in-situ stress coefficient K0, pore water pressure, and saturation, thereby providing a basis for the stable reproduction of shock wave propagation, reflection / diffraction, and energy transdomain transfer processes.

[0060] Step 1.3: Implement differentiated meshing, locally densify the key parts of the underground structural domain (i.e., stress concentration areas such as arch waist / arch foot, opening corner, and interface), and adopt a meshing method of progressively thickening the mesh from the cavern to the outside for the surrounding rock domain, and make the characteristic size of the Eulerian domain element corresponding to the explosive gas match the wavelength to ensure the accuracy of wave propagation calculation.

[0061] Step 1.4: Set contact and interface conditions. The underground structure and the surrounding rock adopt a penalty function-Coulomb friction model with bond-slip degradation function, where the friction coefficient μ of the model can be parameterized and adjusted. The gas-solid interface is coupled using ALE contact or volume fraction exchange algorithm.

[0062] Step 2, Material Model and Equation of State Selection: Define material parameters and select corresponding constitutive models for the underground structural domain and the surrounding rock domain, and define material parameters and select equation of state for the gas domain; at the same time, set the constitutive or contact model of the underground structural-surrounding rock interface and the ALE fluid-structure coupling condition of the gas-solid interface.

[0063] Specifically, the gas domain adopts the JWL or equivalent EOS equation of state for the explosion products, wherein the JWL equation is:

[0064] ,

[0065] In the formula, p For pressure, VLet A and B be the relative specific volume, E be the unit initial internal energy, and A and B be the values ​​of A, B, and E. R 1, R 2, ω The material constant is used only to determine the initial energy level and the selection of the JWL constant, and is not related to the manufacture of the explosive device.

[0066] The surrounding rock domain adopts a constitutive model based on the Drucker-Prager criterion or the Mohr-Coulomb criterion, and introduces a strain softening effect, specifically manifested as the gradual degradation of cohesion and friction angle with the increase of equivalent plastic strain. For soft soil or fractured rock conditions, volume compressibility-related parameters and shear modulus degradation characteristics can be further introduced to adapt to the mechanical response simulation under complex geological conditions.

[0067] The underground structural domain adopts a concrete damage plastic or equivalent damage model that considers tension-compression asymmetry and strain rate effect. The strain rate effect is used to correct the strength and fracture energy of concrete through a dynamic gain coefficient. For reinforcing bars and connectors, embedded reinforcement or "solid-reinforcement coupling" modeling method is adopted.

[0068] In addition, the interface between the lining and the surrounding rock adopts a bond-slip degradation model, which includes core parameters such as peak bond strength, critical slip amount and residual bond strength range, to accurately simulate the mechanical behavior of the interface from bond bearing to slip degradation.

[0069] This step integrates the Explosive Gas Equation of State (EOS), the Mohr-Coulomb / DP constitutive model of the surrounding rock (including softening and dilatation), and the damage plasticity of concrete (including tensile-compressive asymmetry and strain rate effects), and couples it with the lining-surrounding rock bond-slip degradation model, thereby realistically depicting key mechanisms such as cracking-slip-reclosure and local crushing, significantly improving the simulation realism of material-level response.

[0070] Step 3, Application of Initial Field, Boundary and Initiation Conditions: Initial and boundary conditions related to the explosion load and geological environment are applied to the computational domain. Specifically, by setting a non-reflective boundary at the outer boundary of the surrounding rock and applying an initial geostress / pore pressure field, the stress wave propagation process is ensured to be unaffected by boundary reflection. At the same time, by combining the lining-surrounding rock contact / interface model and the gas-solid ALE coupling condition, an energy transfer and interaction path is formed from the explosion gas to the surrounding rock to the lining, thereby providing a consistent physical basis for subsequent dynamic response and damage assessment.

[0071] like Figure 3As shown, to ensure the accuracy and stability of the numerical solution for the dynamic response of underground structures under explosive loading, this embodiment collaboratively sets the initial field, boundary conditions, and interface contact / coupling conditions within the same computational domain. These three elements together constitute the solution framework for the explosive dynamic response. Firstly, the initial field (such as initial ground stress / pore pressure) characterizes the in-situ state of the geological environment, avoiding non-physical abrupt changes at the beginning of the calculation. Secondly, boundary conditions (such as non-reflective boundaries and necessary constraint conditions) suppress artificial reflections of stress waves at the boundaries of the computational domain and eliminate the influence of rigid body motion, thereby ensuring the authenticity of the wave propagation process. Thirdly, the interface contact / coupling conditions described in step 2 describe the interaction and force transmission mechanism (including contact, separation, friction, and gas-solid coupling) between the explosive gas domain and the surrounding rock / structure domain, ensuring that the explosive pressure and energy can be continuously transmitted along the "gas-surrounding rock-structure" path. Therefore, Figure 3 This reflects the way the present invention unifies loading, constraint and interaction into the same framework, providing a foundation for the reliable output of subsequent pressure wave propagation, energy distribution and damage grading results.

[0072] Specifically, an initial triaxial in-situ stress field determined by the burial depth and the in-situ stress coefficient K0 is applied to the surrounding rock domain and / or underground structure domain, and gravity is superimposed to ensure that the initial stress state of the three-dimensional computational domain is consistent with the actual engineering geological conditions. For the surrounding rock domain, an initial pore water pressure field is applied according to the pore or groundwater conditions, or an equivalent hydrostatic pressure field obtained from solid-fluid coupling calculation is used as the pore pressure loading to accurately simulate the influence of groundwater on the mechanical properties of the surrounding rock and soil and the stress state of the structure. A non-reflective boundary is set at the outer boundary of the surrounding rock domain. The non-reflective boundary is a viscous boundary or an infinite element boundary. Furthermore, the top boundary can be a free surface or a road surface-soil composite boundary to adapt to different surface conditions and improve the matching degree between the boundary conditions and the actual engineering conditions.

[0073] In addition, an initiation zone is defined in the gas domain, and an initial energy and pressure growth process is assigned to the initiation zone without involving the manufacture of a specific device, so as to provide a controllable, reproducible and numerically stable equivalent input for the explosion source term.

[0074] Step 4, Explicit Dynamics ALE Remapping Solution: Explicit dynamics integration is performed on the computational domain, with real-time monitoring and recording of various energy parameters. This step implements numerical solution and robustness control measures for the aforementioned three-dimensional computational domain numerical analysis process, and establishes an energy monitoring criterion for external work, internal energy, kinetic energy, and viscoplastic energy dissipation to ensure computational stability and repeatability under strong impact conditions.

[0075] Specifically, a remapping algorithm is employed for the gas domain, and the time step is controlled based on the Courant-Friedrich-Liouvi CFL stability condition to effectively avoid mesh penetration and numerical oscillation problems, ensuring the stability of the computation process. For the solid domain, mass scaling is performed within a limited range to increase the time step, thereby optimizing overall computational efficiency. Furthermore, multiple energy parameters, including external work, kinetic energy, internal energy, viscoplastic energy dissipation, and artificial viscous dissipation, are monitored and recorded in real time. These recorded energy parameters are used as convergence and reliability criteria for the numerical calculation results, verifying their validity.

[0076] Step 5, Result Extraction and Multi-Indicator Fusion and Grading: By extracting the solution results from Step 4, a multi-damage index system and damage grading determination are constructed, and targeted engineering design suggestions are generated. This step constructs a comprehensive grading system (L0–L4) composed of peak vibration velocity (PPV), tensile damage volume fraction, plastic work density, crack continuity index (CI), and stiffness reduction rate, and automatically identifies weak parts such as arch waist, arch foot, and joints by using damage hotspot connectivity and energy density threshold.

[0077] Specifically, it includes the following steps:

[0078] Step 5.1: Extract the time history response parameters of key measurement points from the numerical calculation results, specifically including stress, strain, equivalent plastic strain, damage factor D, particle velocity, and peak velocity (PPV). Figure 5 As shown, the spatial distribution diagrams of the lining damage factor D and PPV are presented. PPV reflects the peak level of particle velocity under explosive dynamic loading and is used to characterize the vibration response intensity. The damage factor D is used to characterize the degree of damage accumulation in the lining material (values ​​range from 0 to 1, with larger values ​​indicating more severe damage). The spatial distribution shows that high PPV areas typically exhibit a significant spatial correspondence and overlap with high D areas in the lining, indicating that the stronger the vibration response and the more concentrated the energy input, the more prone the lining is to damage evolution and accumulation. However, in local locations, there may be situations where PPV is high and D is low, or vice versa. This is because differences in material constitutive parameters, constraint conditions, stress triaxiality, and interface contact states can affect the damage threshold and damage growth rate. Based on the above correspondence, this embodiment uses PPV as a dynamic response intensity index and D as a damage degree index to identify weak points in the lining in step 5.4 through their spatial overlap, providing a basis for subsequent damage classification and risk zoning.

[0079] In addition, this embodiment also outputs isochronous cross-sectional diagrams or isosurface diagrams of the spatial field, specifically covering pressure wavefront distribution, energy density distribution, and damage cloud map, to intuitively present the evolution characteristics of physical quantities within the three-dimensional computational domain. For example... Figure 4This diagram displays the evolution of isosurfaces of pressure wavefront and energy density distribution under typical operating conditions. The pressure wavefront isosurface characterizes the propagation location, arrival time, and attenuation process of the explosion shock wave within the computational domain. The energy density isosurface characterizes the spatial distribution, transmission path, and dissipation accumulation characteristics of the explosion input energy in the gas domain, surrounding rock domain, and lining underground structure domain. As time progresses, the pressure wavefront expands outward from the initiation zone and gradually attenuates. High energy density areas are typically distributed near the pressure wavefront and in its subsequent accompanying region, migrating and spreading along the propagation direction, reflecting the dynamic process of the shock wave's work on the medium and energy input / dissipation. The correspondence between the two indicates that the pressure wavefront determines the "range and timing of the impact," while high energy density areas reveal the "concentrated location of energy input and dissipation." When the wavefront sweeps past the vicinity of the structure, the energy density at the structural surface and interfaces is locally amplified (affected by geometric discontinuities, material differences, and interface interactions). This high-energy area often exhibits spatial consistency with damage hotspots in subsequent damage contour maps, thus providing an intuitive basis for identifying weak points and determining damage grading.

[0080] Step 5.2: Based on the results obtained in Step 5.1, construct a multi-damage index system, where each index is defined as follows:

[0081] I1: Peak velocity (PPV), used as an evaluation index for vibration response intensity; I2: Tensile damage volume fraction, specifically defined as the volume fraction of the lining domain that simultaneously satisfies (i) the maximum principal stress. 1>0 (under tensile stress) and (ii) damage factor > 0 represents the proportion of the unit volume to the total lining volume, used to quantify the tensile damage range; I3: plastic work density, i.e., the cumulative plastic work W per unit volume. p I4: Crack Continuity Index (CI), which is used to characterize the cumulative degree of structural plastic deformation and energy dissipation, and is measured by parameters such as the connectivity of damage hotspots, maximum penetration rate, and path ratio to quantify the development and penetration degree of cracks; I5: Stiffness Reduction Rate, defined as the equivalent circumferential stiffness of the key section of the lining. ( Relative to initial stiffness The relative decrease of 0, i.e. ;in, ( The circumferential stiffness (or equivalent circumferential modulus) can be calculated from the secant stiffness of the circumferential stress-circumferential strain curve of the cross section and is used to evaluate the degree of structural stiffness degradation.

[0082] Step 5.3: Preset the grading thresholds and fusion weights for damage indices I1–I5, and output damage levels from L0 (no obvious damage) to L4 (penetrating damage) based on weighted scoring or logical thresholds, automatically marking weak points (arch / arch foot, interface, seam / hole, etc.). Figure 6 One embodiment of a damage grading criterion is given. The grading threshold is determined based on standard limits, experimental data, or simulation calibration, and the fusion weight is determined based on analytic hierarchy process (AHP), entropy weighting, or data-driven methods.

[0083] Step 5.4: Overlay the damage grading results with the weak points to generate engineering design recommendations that include local reinforcement, interface enhancement, allowable equivalent, and minimum safe burial depth suggestions.

[0084] Step 6, Multi-condition Batch Processing and Sensitivity Analysis: Establish a batch case generation and automated solution process, build a proxy model for sensitivity ranking, rapid pre-assessment and solution selection, and support decision-making at the "solution space" level.

[0085] Specifically, it includes the following steps:

[0086] Step 6.1: Generate multiple sets of working condition samples using Latin hypercube design or orthogonal design methods. The parameters covered by the working condition samples include equivalent charge, charge arrangement or detonation arrangement parameters, underground structure burial depth, K0, pore water pressure, and material zoning combination / material parameter group. Step 6.2: Perform batch automatic calculation on the multiple sets of working condition samples. After the calculation is completed, extract the damage index I1–I5 under each working condition. Step 6.3: Construct a surrogate model based on the Kriging model or response surface model.

[0087] Step 7, Engineering calibration of parameters and thresholds: Based on field monitoring data or experimental data, the model parameters and graded thresholds are inverted and calibrated, which effectively reduces the uncertainty of the model and improves the transferability and confidence of results across sites and working conditions.

[0088] Specifically, this includes: Step 7.1, using on-site monitoring data or test data as observations, employing Bayesian update method or least squares inversion method to invert and calibrate material parameters and interface parameters of underground structure-surrounding rock.

[0089] Step 7.2: Based on the results of the material parameters and interface parameters inversion calibration in Step 7.1, update the grading thresholds of L0–L4 corresponding to the I1–I5 damage indices in Step 5.4 and the multi-index fusion weights; and propagate the parameter uncertainty to the damage indices and grading judgment results, outputting damage grading results with uncertainty intervals (confidence / credibility intervals), and update and optimize the engineering design suggestions in Step 5.5 accordingly.

[0090] Furthermore, this method can automatically generate visual charts and textual conclusions related to numerical calculations, such as damage level diagrams required for blast resistance verification, reinforcement recommendations for key components, and recommended ranges for minimum safe burial depth and allowable equivalent yield. Moreover, by saving metadata (geometric model, mesh, working conditions, parameters), it enables model version tracking and cross-site reuse, directly serving design optimization and risk management.

[0091] Simulation Example:

[0092] The present invention will be further illustrated below with a typical engineering example, but the scope of protection of the present invention is not limited thereto.

[0093] Project background and design objectives: A subway section in a certain city adopts a uniform-thickness annular reinforced concrete lining structure with an inner radius of R=3.30 m, a lining thickness of t=0.35 m, and a soil cover thickness of approximately 12-15 m. Due to its location in the city's core area, to address the potential risk of near-field ground explosions, it is necessary to assess the stress response, vibration velocity level, and damage level of the tunnel structure within the equivalent TNT weight range of 8-15 kg, identify weak components, and propose safety recommendations.

[0094] To meet the above design requirements, the "Dynamic Response Simulation and Damage Classification Method for Underground Structures Based on ALE Fluid-Structure Coupling" proposed in this invention was applied. The specific implementation process is as follows:

[0095] Step A1: Building the Model and Setting Boundaries

[0096] Step A1.1, Geometric model construction and mesh generation: Based on circumferential symmetry, a three-dimensional 60° periodic sector model is constructed with a length of 6.0m; Lagrange solid elements are used for both the underground structure (lining) and the surrounding rock; the explosive gas domain is described using ALE; the lining mesh is refined to 40mm, and the surrounding rock mesh is locally refined to 80mm; the explosive gas domain mesh size is 50mm, matching the main impact wavelength.

[0097] Step A1.2, Initial field and boundary conditions application: A viscous absorbing boundary is set on the outside of the surrounding rock model to suppress reflected waves; the initial in-situ stress is set according to the burial depth (in-situ stress coefficient K0=0.5, the letters γ=18 kN / m). 3 The lining-surrounding rock interface adopts a penalty function contact model, with a friction coefficient μ=0.65 and a normal stiffness set to Kn=1×10. 9 N / m 3 And the adhesion-slip degradation model was selected.

[0098] Step A1.3, Application of Explosion Conditions: The detonation point is located at the center of the ground surface, using a spherical detonation zone, equivalent to...

[0099] The TNT equivalents are 8 kg, 12 kg, and 15 kg, respectively. The initial energy of the explosion gas domain is applied using the JWL equation of state.

[0100] Step A2, configure the material model, as shown in Table 1 below:

[0101]

[0102] Step A3: Perform simulation and energy analysis.

[0103] The explicit dynamics + ALE remapping algorithm was used for the solution, with the simulation time set to 30ms, and the energy conservation error was controlled within ±4%. The output includes: peak particle velocity (PPV) time history; damage factor D inside the lining, equivalent plastic strain distribution; and changes in energy terms such as external work, kinetic energy, internal energy, viscoplastic dissipation, and artificial viscous dissipation over time.

[0104] Step A4: Perform results analysis and multi-indicator fusion and grading.

[0105] The response results under typical operating conditions are shown in Table 2 below (taking 12 kg as an example):

[0106]

[0107] Based on preset grading standards, damage levels from L0 (no obvious damage) to L4 (through failure) are automatically output. In this case, the crown / waist area is classified as L2–L3, indicating obvious tensile cracking and interface slippage; the arch foot area is next (L1–L2). The preset grading standards are shown in Table 3 below:

[0108]

[0109] Step A5: Identify weak points and provide engineering design recommendations.

[0110] The automatically identified weak points include high-damage concentration areas located in the arch waist-arch foot transition area, the lining-surrounding rock contact surface, and construction joint locations; areas where energy density hotspots and damage volume fractions overlap are key areas of concern. Furthermore, the following blast-resistant engineering design recommendations are given: double-layer reinforced concrete lining is recommended for the arch waist / arch foot area; high-strength interface mortar / sprayed binder is recommended to enhance force transmission at the lining-surrounding rock interface; the minimum safe burial depth is recommended to be ≥15.5 m (8-12 kg TNT) and ≥18 m (15 kg TNT); if the burial depth is insufficient or there is surface hard cover, a ground energy buffer layer or blast barrier should be installed.

[0111] Therefore, this simulation embodiment fully demonstrates the entire process of this invention under the same typical working condition, from 3D modeling, material and boundary / contact settings, robust ALE fluid-structure interaction solution, key time history parameter extraction, multi-index system construction and fusion classification, automatic identification of weak points, and design suggestion output. Simulation results show that the pressure wave advance and the migration patterns of high energy density areas have spatial consistency with damage hotspots, which can support the automatic location of weak points; the multi-index system can quantitatively characterize "response intensity - damage range - plastic energy dissipation - crack penetration - stiffness degradation" and output classification results; the reinforcement / optimization suggestions generated accordingly have clear positioning basis and traceable evaluation indicators.

[0112] Thus, it is verified that the method of the present invention can simultaneously output a closed-loop result chain of "impact propagation - energy distribution - damage evolution - classification determination - suggestion generation" in a single calculation case, and has the ability to be implemented, computationally stable and quantitatively evaluated in engineering applications.

[0113] The specific embodiments of the present invention have been described in detail above with reference to the figures, but the present invention is not limited to the described embodiments. For those skilled in the art, various changes, modifications, substitutions, and variations can be made to these embodiments without departing from the principles and spirit of the present invention, and these variations still fall within the protection scope of the present invention.

Claims

1. A method for simulating the dynamic response and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction, characterized in that, Includes the following steps: Step 1: Establish a three-dimensional computational domain that includes the underground structure domain, the surrounding rock domain, and the explosive gas domain, and define the geometric parameters and working parameters and set the mesh. Step 2: Define material parameters and select corresponding constitutive models for the underground structural domain and the surrounding rock domain, and define material parameters and select equations of state for the gas domain; at the same time, set the constitutive or contact model of the underground structural-surrounding rock interface and the ALE fluid-structure coupling condition of the gas-solid interface. Step 3: Apply initial and boundary conditions related to the explosion load and geological environment to the computational domain; Step 4: Perform explicit dynamic integration on the computational domain, and monitor and record multiple energy parameters in real time; Step 5: Construct a multi-damage index system and damage classification by extracting the solution results from Step 4, and formulate targeted engineering design suggestions. Specifically, it includes: Step 5.1: Extract the time history response parameters of key measuring points from the numerical calculation results, and output the isochronous cross section or isosurface map of the spatial field. Step 5.2: Based on the results obtained in Step 5.1, construct a multi-damage index system including I1 peak vibration velocity PPV, I2 tensile damage volume fraction, I3 plastic work density, I4 crack continuity index and I5 stiffness reduction rate. Step 5.3: Preset the grading thresholds and fusion weights for damage indices I1–I5, and output damage levels L0 to L4 based on weighted scoring or logic thresholds, automatically marking weak areas. The grading thresholds are determined based on standard limits, experimental data, or simulation calibration, and the fusion weights are determined based on analytic hierarchy process, entropy weighting, or data-driven methods. Step 5.4: Overlay the damage grading results with the weak points to generate engineering design recommendations including local reinforcement, interface enhancement, allowable equivalent, and minimum safe burial depth recommendations. The time-history response parameters include stress, strain, equivalent plastic strain, damage factor D, particle velocity, and peak velocity PPV; the isochronous cross-sectional diagram or isosurface diagram of the spatial field includes pressure wavefront distribution, energy density distribution, and damage cloud map distribution. Step 6: Establish a batch calculation case generation and automated solution process, and construct a proxy model; specifically, it includes: Step 6.1: Use Latin hypercube design or orthogonal design methods to generate multiple sets of working condition samples. The parameters covered by the working condition samples include equivalent charge, charge arrangement method or detonation arrangement parameters, underground structure burial depth, K0, pore water pressure and material parameter set. Step 6.2: Perform batch automatic calculation on the multiple sets of working condition samples, and extract the damage index I1–I5 under each working condition after the calculation is completed; Step 6.3: Construct a proxy model based on the Kriging model or response surface model; Step 7: Based on field monitoring data or experimental data, perform inversion calibration on model parameters and grading thresholds.

2. The method for simulating dynamic response and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction as described in claim 1, characterized in that: Step 1 specifically includes: Step 1.1: Construct a three-dimensional computational domain, which includes the underground structure domain and surrounding rock domain described by Lagrange, and the explosive gas domain described by Eulerian or ALE. Step 1.2: Parametrically set the geometric parameters and working parameters. The geometric parameters include burial depth, overburden thickness and structural joint parameters. The working parameters include load application method, loading sequence, equivalent charge, in-situ stress coefficient K0, pore water pressure and saturation. Step 1.3: Implement differentiated grid division, locally densify the key parts of the underground structural domain, and adopt a grid division method of progressively thickening from the cavern to the periphery for the surrounding rock domain, and make the characteristic size of the Eulerian domain unit corresponding to the explosive gas match the wavelength. Step 1.4: Set contact and interface conditions. The underground structure and surrounding rock adopt a penalty function-Coulomb friction model with bond-slip degradation function. The gas-solid interface is coupled using ALE contact or volume fraction exchange algorithm.

3. The method for simulating dynamic response and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction as described in claim 1, characterized in that: In step 2, the gas domain adopts the JWL or equivalent EOS equation of state for the explosion products, the surrounding rock domain adopts the constitutive model of the Drucker-Prager criterion or the Mohr-Coulomb criterion, the underground structural domain adopts the concrete damage plasticity or equivalent damage model considering tensile-compressive asymmetry and strain rate effect, and the interface between the lining and the surrounding rock adopts the bond-slip degradation model.

4. The method for simulating dynamic response and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction as described in claim 1, characterized in that: In step 3, an initial triaxial in-situ stress field determined by the burial depth and the in-situ stress coefficient K0 is applied to the surrounding rock domain and / or the underground structure domain, and gravity is superimposed on it; for the surrounding rock domain, an initial pore water pressure field is applied according to the pore or groundwater conditions, or an equivalent hydrostatic pressure field obtained by solid-fluid coupling calculation is used as the pore pressure loading; a non-reflective boundary is set at the outer boundary of the surrounding rock domain, wherein the non-reflective boundary is a viscous boundary or an infinite element boundary; an initiation zone is defined in the gas domain, and the initiation zone is given an initial energy and pressure growth history.

5. The method for simulating dynamic response and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction as described in claim 1, characterized in that: In step 4, a remapping algorithm is applied to the gas domain, and the time step is controlled based on the Courland-Friedrich-Liouvi CFL stability condition; the solid domain is subjected to mass scaling within a limited range; the energy parameters include external work, kinetic energy, internal energy, viscoplastic energy dissipation, and artificial viscous dissipation.

6. The method for simulating dynamic response and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction as described in claim 1, characterized in that: Step 7 specifically includes: Step 7.1: Using on-site monitoring data or test data as observations, the material parameters and the interface parameters of the underground structure and surrounding rock are inverted and calibrated using the Bayesian update method or the least squares inversion method. Step 7.2: Based on the results of the material parameters and interface parameters inversion calibration in Step 7.1, update the grading thresholds of L0–L4 corresponding to the I1–I5 damage indices in Step 5.4 and the multi-index fusion weights; and propagate the parameter uncertainty to the damage indices and grading judgment results, outputting the damage grading results with uncertainty intervals, and update and optimize the engineering design suggestions in Step 5.5 accordingly.

7. The method for simulating dynamic response and damage classification of underground structures under explosive loads based on ALE fluid-structure interaction as described in claim 1, characterized in that: The method also includes: automatically generating visual charts and textual conclusions related to numerical calculations, and realizing model version traceability and cross-site reuse by saving metadata.

Images