A method for evaluating damage evolution of underground structure under sudden explosion load

Abstract

The application discloses a kind of underground structure damage evolution evaluation methods under the action of burst explosion load, it is related to civil engineering and engineering disaster prevention technical field, it is characterized in that, including the following steps: S1: based on explosion load, geometric boundary and material parameters, simulate underground structure, construct unified constitutive model;S2: unified constitutive model introduces the pressure-dependent strength criterion of Drucker-Prager criterion to depict the influence of confining pressure on shear strength.The technical problem to be solved by the application is to provide an underground structure damage evolution evaluation method under the action of burst explosion load, to solve the key pain points such as the difficulty of existing analysis method to simultaneously consider pressure correlation, high strain rate effect, material random discreteness and damage softening energy consistency, to establish a set of damage evolution evaluation method and implementation system that can be implemented from material to structure level, to support quantitative prediction and engineering classification of underground structure from crack initiation, expansion to instability and failure.

Description

Technical Field

[0001] This invention relates to the field of civil engineering and disaster prevention technology, specifically to a method for assessing the damage evolution of underground structures under sudden explosive loads. Background Technology

[0002] With the large-scale development of urban underground space, the number of underground structures such as subway stations and tunnels, underground utility tunnels, municipal and hydraulic engineering caverns, and national defense projects is rapidly increasing. These projects are situated in environments constrained by deep overburden and surrounding rock / soil. Once subjected to sudden explosive loads (such as explosions caused by transportation, storage, construction, or terrorist attacks), they will exhibit dynamic characteristics of high confining pressure, multiaxial stress, and strong velocity effects. The shock wave generated by the detonation will be reflected and superimposed multiple times within enclosed or semi-enclosed spaces. Accompanied by the quasi-static hysteresis of gaseous products, this can easily induce tensile / shear composite cracking and rapid damage propagation at stress concentration points such as the arch crown, arch waist, and corners, leading to load-bearing capacity degradation and connectivity disruption, resulting in secondary disasters and operational interruptions. Therefore, quantitative damage evolution assessment and blast-resistant design methods for explosive impact conditions are urgently needed.

[0003] From the perspective of materials and mechanical mechanisms, commonly used materials in underground structures (concrete, masonry composites, fiber-reinforced cementitious materials, etc.) exhibit significant characteristics under explosive impact: (1) Pressure correlation: Increased confining pressure can enhance shear strength and ductility, and the yield surface changes significantly with the first stress invariant; (2) High strain rate effect: Strength and apparent peak strain increase with increasing loading rate; (3) Quasi-brittleness and damage softening: Microcrack initiation and synergistic propagation lead to stiffness and strength degradation, resulting in softening and localization; (4) Statistical dispersion: Inhomogeneity of internal material defects, aggregate distribution, and interface characteristics leads to significant randomness in crack initiation strain, damage threshold, and failure path. Accordingly, accurately characterizing the coupled behavior of "pressure correlation - rate sensitivity - statistical dispersion - damage softening - energy dissipation" is a key challenge in the analysis of underground structures under explosive conditions.

[0004] Existing technologies mainly include the following types of methods, each with its own obvious shortcomings: 1. Experience and Equivalence Method Based on scaled-down distances, empirical impact response spectra, and the equivalent single degree of freedom (SDOF) method, component displacement or equivalent internal forces can be quickly estimated, making them suitable for early-stage scheme comparison. However, these methods are mostly based on freely propagating shock waves in ground or open environments, making it difficult to reflect multiple reflections, confining pressure effects, and hysteresis in underground environments; moreover, they cannot reveal the damage evolution process of crack initiation, propagation, and penetration, thus providing limited guidance for damage classification and repair decisions.

[0005] 2. Traditional elastoplastic models (such as Mohr–Coulomb, Drucker–Prager, and capped models) These models can consider pressure-dependent yield characteristics and are widely used in soil-structure interaction analysis. However, most implementations do not explicitly introduce damage variables, making it difficult to describe stiffness degradation and unloading path evolution; at high strain rates, fitting only static parameters easily underestimates peak strength and energy dissipation capacity; and they lack characterization of crack localization and penetration paths, failing to provide a damage field consistent with actual failure modes.

[0006] 3. Continuous damage and cracking model Continuous medium models such as isotropic / anisotropic damage, tensile-compressive splitting damage, and concrete damage plasticity (CDP) can describe strength and stiffness degradation, and some models also support tensile-compressive asymmetry. However, common problems include the following steps: (1) The rate effect is simply handled, and the strength is often corrected independently by the calibration factor without being coupled with damage evolution and yield surface; (2) There is a lack of statistical characterization of material discretization and weakening elements, which makes the crack initiation location and failure path sensitive to the mesh and initial defects; (3) The mesh dependence and non-physical energy dissipation problems caused by softening are not systematically solved, which makes the width of the damage zone and energy dissipation unstable with the mesh.

[0007] 4. Specialized blast-resistant concrete models (such as high-pressure, high-rate constitutive families commonly used in the engineering field). These models comprehensively consider confining pressure, strain rate, and fragmentation evolution, making them well-suited for penetration, intrusion, and near-field impact problems. However, their parameter systems are complex, calibration costs are high, and they often rely on specific test platforms and large-sample impact tests. Furthermore, they primarily focus on material element-level responses, lacking engineering outputs linked to structural damage grading and repairability assessments. Some implementations still suffer from insufficient softening regularization and numerical stability issues.

[0008] 5. Discrete / Lattice and Phase-Field Fracture Methods Discrete element method (DEM), lattice / microbeam network, and phase field fracture method can naturally generate crack networks and penetration paths, and have strong ability to describe fracture topology. However, their parameter and scale mapping relationships are complex and computationally expensive; their application in the overall analysis of large-scale underground structures is limited, and they are not well integrated with soil-structure interaction and engineering classification processes.

[0009] 6. Numerical stability and parameter inversion problems The stress-strain curve caused by explosive impact exhibits a significant softening segment. Without energy regularization constrained by fracture energy or a viscous / overstress stabilization mechanism, non-physical localization and divergent results often occur. On the other hand, data from common engineering tests such as triaxial tests and split Hopkinson bar (SHPB) tests have not established a robust automatic inversion process with the model parameter system, resulting in subjective parameter selection and insufficient repeatability, thus weakening the model's engineering application value.

[0010] Therefore, the engineering community needs a constitutive and assessment framework for underground structure explosion impact scenarios that can uniformly couple pressure-related yielding (such as DP), strain rate amplification, viscous / overstress flow, and Weibull statistical damage. Simultaneously, it needs to be equipped with energy regularization to ensure mesh independence during the softening stage, and to form a parameterizable, calibrable, and graded output engineering process and software system to support the prediction, damage level determination, and reinforcement decisions throughout the entire process from crack initiation to instability failure. This is precisely the urgent direction and practical need for technological development in this field. Summary of the Invention

[0011] The technical problem to be solved by this invention is to provide a method for assessing the damage evolution of underground structures under sudden explosive loads. This method addresses the key pain points of existing analysis methods, such as the inability to simultaneously consider pressure correlation, high strain rate effect, material random discreteness, and consistency of damage softening energy. It establishes a damage evolution assessment method and implementation system that can be implemented from the material to the structural level, supporting the quantitative prediction and engineering classification of the entire process of underground structures from crack initiation to instability and failure.

[0012] The present invention achieves its objective by employing the following technical solution: A method for assessing the damage evolution of underground structures under sudden explosive loads, characterized by the following steps: S1: Based on explosive loads, geometric boundaries, and material parameters, simulate underground structures and construct a unified constitutive model; S2: The unified constitutive model introduces the Drucker-Prager criterion pressure-related strength criterion to characterize the influence of confining pressure on shear strength; S3: The unified constitutive model achieves the coupled evolution of high strain rate-damage-yield through rate effect and viscous flow mechanism; S4: The unified constitutive model introduces the Weiber statistical framework to characterize the strength dispersion and crack initiation probability caused by material defects; S5: The unified constitutive model adopts a softening regularization strategy through fracture energy constraint and characteristic length control to suppress mesh dependence and non-physical localization, ensuring that energy dissipation and crack width in the softening stage are mesh-independent, thereby improving the repeatability and robustness of the simulation. S6: The unified constitutive model after softening regularization is further processed using numerical integration and algorithmic framework; S7: Establish a parameter calibration and inversion process based on S6; S8: Provide damage grading based on S6.

[0013] As a further limitation of this technical solution, step S1 includes the following steps: S11: Explosion Load History: Input Near / Mid / Far Field Equivalent Pressure – Time History p(t), or import the wall pressure distribution obtained from computational fluid dynamics / fluid-structure interaction calculations. p (x, t), allowing a given rise time t r Peak p max hysteresis p.s. With decay constant k d ; S12: Geometry and Boundaries: Establish a 3D underground structure—surrounding rock / soil model, and set the contact (normal stiffness K). n tangential friction coefficient m ) and constraint / radiation boundary; S13: Materials include cement-based, rock-based, or composite materials: Elasticity: Elasticity model E∈[15,60]GPa, Poisson's ratio n ∈[0.15,0.25]; Strength (DP equivalent): The effect of friction angle Φ on yield criterion α∈[0.05,0.35], a material constant related to cohesion c. k ∈[1,20] MPa, or converted from (c, Φ); Statistical damage (Weibull): shape parameter m∈[2,20], scale parameter ∈[5×10⁻⁴, 5×10⁻³]; Viscosity / Overstress: Coefficient of Viscosity or ∈[10 -3 10 -1 ]s, exponent n∈[1,6], threshold σ0∈[0.1,5]MPa; Strain rate amplification: Strain rate sensitivity coefficient C ∈[0.5,50]s⁻¹, exponential parameter / curve shape parameter p ∈[0.2,1.2]; Fracture energy: tensile crack G f t ∈[50,250]N / m, shear G f s ∈[100,500]N / m; Feature length (cell-dependent): Feature length l ch Values ​​are automatically retrieved based on the grid scale.

[0014] As a further limitation of this technical solution, step S2 includes the following steps: S21: Damage equivalent stress; (1) in: For the damage equivalent stress tensor, D As a damage variable, For total strain tensor and It is the plastic strain tensor; S22: Drucker-Prager criterion for yielding and hardening degradation; (2) in: I 1 is the first invariant. J 2 represents the second deviatoric stress invariant; For intensity rate correction; F This represents a Drucker–Prager type yield function.

[0015] As a further limitation of this technical solution, step S3 includes the following steps: S31: Intensity rate amplification adopts an example power law or exponential type; The example power law is: Exponential type is: (3) Wherein: parameters (a,b) are fitted by experiments; denoted as equivalent strain rate; a is the strain rate sensitivity coefficient; b is the strain rate exponent parameter. Equivalent plastic strain; S32: Overstressed viscoplasticity; (4) in: For the plastic strain rate tensor, The derivative of the yield function, These are components of the stress tensor; As a further limitation of this technical solution, step S4 includes the following steps: S41: Probability of crack initiation; (5) S42: Damage Growth Law; (6) in: l and β For calibration coefficients, F 0 is the normalization constant; D For damage variables; Optional tensile / compressive splitting damage D t and D c Driven by the principal tensile / compressive strain respectively, satisfying: D = xD t +(1- x ) D c (7) Where: χ is the tension / compression weighting factor.

[0016] As a further limitation of this technical solution, step S5 includes the following steps: To eliminate mesh dependence, a fracture energy-controlled softening modulus correction is employed:

[0017] in: s pk Peak intensity; according to l ch Automatically adjust the softening slope to ensure conservation of unit fracture energy; G f Values G f t or G f s .

[0018] As a further limitation of this technical solution, step S6 includes the following steps: Time integration: using explicit central difference or implicit incremental Newton method; Return mapping: Tangentially consistent plastic correction is applied to the yield surface of the Drucker-Prager criterion; rate amplification and viscous flow are coupled and updated in the sub-iteration layer.

[0019] Substep size adaptive: based on F , With energy error control time step; Contact / interface: Combining the normal penalty function with Coulomb friction, and introducing an interfacial bond-slip-softening segment where necessary to describe the degradation of the lining-soil interface.

[0020] As a further limitation of this technical solution, step S7 includes the following steps: S71: Data source; Static / quasi-static: Single-axis and triaxial (E, ν, c, Φ) and peak / residual values ​​are obtained. High strain rate: Rate sensitivity curves and fracture energy were obtained from split Hopkinson bar tension / compression and shear bar tests.

[0021] S72: Inversion objective function:

[0022] Where: parameter vector ; J(p) is the objective function / cost function used for parameter inversion (identification); w i In order to be with the first i Weighting coefficients (dimensionless) corresponding to the group data. For numerically predicted stress vectors / sequences; For the first i Measured stress vectors / sequences from the group of tests; w g The weight (dimensionless) of the fracture energy (or energy dissipation) constraint term. The fracture energy / energy dissipation per unit area is calculated by the numerical model under parameter p. The fracture energy was measured experimentally. In response to the situation; t For time; S73 Algorithm: Global and local hybrid optimization, cross-validation to prevent overfitting; output parameter confidence intervals.

[0023] As a further limitation of this technical solution, step S8 includes the following steps: Let the threshold be 0 < D 1< D 2< D 3≤1: D ∈[0, D 1): Micro-damage; [ D 1, D 2: Moderate injury; [ D 2, D 3): Severe injury; ≥ D 3: Penetration / Instability; Combined with residual bearing capacity or R = R res / R 0 ( or R Residual bearing capacity R res For residual bearing capacity, R The initial bearing capacity is 0, and the crack connectivity and functional indicators are comprehensively evaluated.

[0024] Recommended treatment methods: thickening lining, circumferential / longitudinal fiber composite reinforcement, spraying backfill, interface reinforcement, pressure relief groove, etc.

[0025] As a further limitation of this technical solution, it also includes S9: result output and report, outputting the results and reports of S7 and S8.

[0026] Compared with the prior art, the advantages and positive effects of the present invention are: This invention constructs a unified constitutive model: it introduces the Drucker–Prager (DP) pressure-related strength criterion to characterize the influence of confining pressure on shear strength, introduces the Weibull statistical framework to characterize the strength dispersion and crack initiation probability caused by material defects, and realizes the coupled evolution of high strain rate-damage-yield through rate amplification and viscous / overstress flow mechanism, thereby obtaining a description of the strength and stiffness degradation process that matches the characteristics of explosive impact.

[0027] This invention achieves numerical stability and energy consistency: by using a softening regularization strategy that constrains fracture energy and controls characteristic length, it suppresses mesh dependence and non-physical localization, ensuring that energy dissipation and crack width during the softening stage are mesh-independent, thereby improving the repeatability and robustness of the simulation.

[0028] This invention establishes a parameter calibration and inversion process: it forms a parameter identification method based on conventional triaxial tests and high strain rate (such as SHPB) tests as data sources, automatically inverts statistical damage parameters, rate-sensitive parameters and viscosity parameters, and realizes a traceable and reusable engineering closed loop from test data to model parameters.

[0029] This invention provides engineering classification and decision output: setting classification thresholds that match damage variables, residual load and crack continuity, outputting engineering criteria such as "repairable / reinforcement required / replacement required", and providing corresponding structural reinforcement and material optimization suggestions, which can be directly applied in design, evaluation and emergency response.

[0030] This invention provides a practical software and system implementation: it offers standalone applications or user subroutines (UMAT / VUMAT) and supporting modules (load generation, material constitutive modeling, solution stabilization, parameter identification, evaluation classification, and report visualization) for general-purpose finite element platforms, to meet the needs of various scenarios such as scientific research analysis, engineering design, and rapid evaluation.

[0031] In summary, this invention, through an integrated technical approach of "DP pressure-related yielding—Weibull statistical damage—rate and viscosity coupling—energy regularization—engineering graded output," strives to achieve a balance between mechanistic integrity, numerical robustness, and engineering usability, providing reliable technical support for blast-resistant design and post-damage assessment of underground structures. Attached Figure Description

[0032] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a schematic diagram of the unified constitutive model and yield surface of the present invention; Figure 3 The Weibull probability density and damage growth curves of this invention are shown below; Figure 4 This is a schematic diagram of crack width control under energy regularization according to the present invention; Figure 5 This is a typical underground lining structure grid and boundary of the present invention; Figure 6 This is a schematic diagram of the damage grading of the present invention.

[0033] Figure 7 This is a schematic diagram of the simulated underground structure of the present invention. Detailed Implementation

[0034] The following detailed description of a specific embodiment of the present invention is provided in conjunction with the accompanying drawings. However, it should be understood that the scope of protection of the present invention is not limited to the specific embodiment.

[0035] This invention includes the following steps: S1: Based on explosive loads, geometric boundaries, and material parameters, simulate underground structures and construct a unified constitutive model; S2: The unified constitutive model introduces the Drucker-Prager (DP) pressure-related strength criterion to characterize the influence of confining pressure on shear strength. S3: The unified constitutive model achieves the coupled evolution of high strain rate-damage-yield through rate effect and viscous flow mechanism; S4: The unified constitutive model introduces the Weibull statistical framework to characterize the strength dispersion and crack initiation probability caused by material defects; S5: The unified constitutive model adopts a softening regularization strategy through fracture energy constraint and characteristic length control to suppress mesh dependence and non-physical localization, ensuring that energy dissipation and crack width in the softening stage are mesh-independent, thereby improving the repeatability and robustness of the simulation. S6: The unified constitutive model after softening regularization is further processed using numerical integration and algorithmic framework; S7: Establish a parameter calibration and inversion process based on S6; S8: Provide damage grading based on S6.

[0036] S1 includes the following steps: S11: Explosion Load History: Input Near / Mid / Far Field Equivalent Pressure – Time History p (t), or import the wall pressure distribution obtained from computational fluid dynamics / fluid-structure interaction (CFD / FSI) calculations. p (x, t), allowing a given rise time t r Peak p max hysteresis p.s. With decay constant kd ; S12: Geometry and Boundaries: Establish a 3D underground structure—surrounding rock / soil model, and set the contact (normal stiffness K). n tangential friction coefficient m ) and constraint / radiation boundary; S13: Materials include cement-based, rock-based, or composite materials: Elasticity: Elasticity model E∈[15,60]GPa, Poisson's ratio n ∈[0.15,0.25]; Strength (DP equivalent): The effect of friction angle Φ on yield criterion α∈[0.05,0.35], a material constant related to cohesion c. k ∈[1,20] MPa, or converted from (c, Φ); Statistical damage (Weibull): shape parameter m∈[2,20], scale parameter ∈[5×10⁻⁴, 5×10⁻³]; Viscosity / Overstress: Coefficient of Viscosity or ∈[10 -3 10 -1 ]s, exponent n∈[1,6], threshold σ0∈[0.1,5]MPa; Strain rate amplification: Strain rate sensitivity coefficient C ∈[0.5,50]s⁻¹, exponential parameter / curve shape parameter p ∈[0.2,1.2]; Fracture energy: tensile crack G f t ∈[50,250]N / m, shear G f s ∈[100,500]N / m; Feature length (cell-dependent): Feature length l ch Values ​​are automatically retrieved based on the grid scale.

[0037] S2 includes the following steps: S21: Damage equivalent stress; (1) in: For the damage equivalent stress tensor, D As a damage variable, For total strain tensor and It is the plastic strain tensor; S22: Drucker-Prager criterion for yielding and hardening degradation; (2) in: I 1 is the first invariant. J 2 represents the second deviatoric stress invariant; For intensity rate correction; F This represents a Drucker–Prager type yield function.

[0038] S3 includes the following steps: S31: Intensity rate amplification adopts an example power law or exponential type; The example power law is: Exponential type is: (3) Wherein: parameters (a,b) are fitted by experiments; denoted as equivalent strain rate; a is the strain rate sensitivity coefficient; b is the strain rate exponent parameter. Equivalent plastic strain; S32: Overstressed viscoplasticity; (4) in: For the plastic strain rate tensor, The derivative of the yield function, These are components of the stress tensor; pass or The control of rate sensitivity and numerical stability, n, and σ0 avoids non-physical jumps.

[0039] S4 includes the following steps: S41: Probability of crack initiation; (5) S42: Damage growth law (probability-overstress coupling); (6) in: l and β For calibration coefficients, F 0 is the normalization constant; D For damage variables; Optional tensile / compressive splitting damage D t and D c Driven by the principal tensile / compressive strain respectively, satisfying: D = xD t +(1- x ) D c (7) Where: χ is the tension / compression weighting factor.

[0040] S5 includes the following steps: To eliminate mesh dependence, a fracture energy-controlled softening modulus correction is employed:

[0041] in: s pk Peak intensity; according to l ch Automatically adjust the softening slope to ensure conservation of unit fracture energy; G f Values G f t or G f s .

[0042] S6 includes the following steps: Time integration: Explicit central difference (suitable for strong transient impacts) or implicit incremental Newton method (suitable for coupled hysteresis / contact). Return mapping: Tangentially consistent plastic correction is applied to the yield surface of the Drucker-Prager criterion; rate amplification and viscous flow are coupled and updated in the sub-iteration layer.

[0043] Substep size adaptive: based on F , With energy error control time step; Contact / interface: Combining the normal penalty function with Coulomb friction, and introducing an interfacial bond-slip-softening segment where necessary to describe the degradation of the lining-soil (rock) interface.

[0044] S7 includes the following steps: S71: Data source; Static / quasi-static: Uniaxial and triaxial (different confining pressures) yield (E, ν, c, Φ) and peak / residual values; High strain rate: Rate sensitivity curves and fracture energy were obtained from the split Hopkinson bar (SHPB) tension / compression and shear bar tests.

[0045] S72: Inversion objective function:

[0046] Where: parameter vector ; J (p) is the objective function / cost function used for parameter inversion (identification); w i In order to be with the first i Weighting coefficients (dimensionless) corresponding to the group data. For numerically predicted stress vectors / sequences; For the first i Measured stress vectors / sequences from the group of tests; w g The weight (dimensionless) of the fracture energy (or energy dissipation) constraint term. The fracture energy / energy dissipation per unit area is calculated by the numerical model under parameter p. The fracture energy was measured experimentally. In response to the situation; t For time; α is the yield function parameter, which is related to the friction angle Φ; k The strength constant is related to the cohesion. c Related; m For Weibull shape parameters; For Weibull scale parameters; or The viscosity coefficient; n It is the overstress index; C The strain rate sensitivity coefficient; p The strain rate amplification index; G f t It is the tensile fracture energy; G f s This represents the shear fracture energy.

[0047] S73 Algorithm: Hybrid optimization of global (genetic / particle swarm optimization) and local (quasi-Newton / LM) approaches, cross-validation to prevent overfitting; output parameter confidence intervals.

[0048] S8 includes the following steps: Let the threshold be 0 < D 1< D 2< D 3≤1: D ∈[0, D 1): Minor damage (repairable); [ D 1, D 2: Moderate injury (requires local reinforcement); [ D 2, D 3): Severe damage (significant degradation of structural load-bearing capacity); ≥ D 3: Through-through / Instability (requires replacement or major reinforcement); Combined with residual bearing capacity or R = R res / R 0 ( or R Residual bearing capacityR res For residual bearing capacity, R The initial bearing capacity (0), crack connectivity (based on damage isosurface connectivity analysis), and functional indicators (displacement / crack width limit) are used to comprehensively determine the grade.

[0049] Recommended treatment methods: thickening lining, circumferential / longitudinal fiber composite reinforcement, spraying backfill, interface reinforcement, pressure relief groove, etc.

[0050] It also includes S9: Result Output and Report, which outputs the results and reports of S7 and S8.

[0051] S9 includes the following steps: Time history curves: response spectrum, displacement / acceleration, energy dissipation decomposition (elastic / plastic / viscous / damage).

[0052] Spatial field: equivalent plastic strain, damage isosurface, main crack path, and ranking of critical components.

[0053] Reports: Automatically generate portable document format / Hypertext Markup Language (PDF / HTML), including parameter tables, grading maps and hardening recommendation lists, and export neutral data (Visualization Toolkit / comma-separated values, VTK / CSV).

[0054] Example: Project Background and Objectives A subway section in a certain city uses a ring-shaped reinforced concrete lining of equal thickness, with an inner radius of... R =3.30m, lining thickness t =0.35m, with a soil cover thickness of approximately 12-15m. The design requires verification of blast resistance under the equivalent impact of 8-15kg TNT near-field blast, and obtaining the damage classification.

[0055] Geometry, Mesh and Boundary A three-dimensional 60° periodic sector model (considering circumferential symmetry) with a length of 6.0m was established, featuring lining-soil coupling; the lining utilizes 8-node solid elements with a minimum mesh size h. min =40mm (locally reinforced at the arch waist / arch foot), finest surrounding rock h min =80mm; Lining-soil normal penalty function, tangential Coulomb friction coefficient m =0.65, interface normal stiffness K n =1.0×10 9 N / m 3 The outer domain is set with a viscoelastic / radiative boundary, and the damping is assumed to be ξ=4%.

[0056] Materials and model parameters (according to the technical solution of this invention) Concrete E=32GPa, ν=0.20; Equivalent elasticity of surrounding rock Es=120MPa, ν s =0.30. Drucker–Prager (DP) strength (equivalent to (c,ϕ): α=0.12, k =6.5MPa (corresponding to an equivalent conversion of c≈1.5MPa, Φ≈45°). Weibull statistical damage, shape parameter m=7, dimensional strain. =0.0025. Viscosity / Overstress, or =0.020s, n=3, σ0=0.5MPa. Strain rate amplification (power law). C =5s -1 , p =0.6. Fracture energy and regularization, tensile fracture. G f t =120N / m, shear crack G f s =220N / m, characteristic length l ch =1.2 h min Interfacial bonding-slip-softening (optional), peak bond strength 1.2 MPa, critical slip 0.08 mm, residual 20%.

[0057] Explosive load input Using equivalent wall pressure time history: Operating condition definition: Operating Condition A (8kg TNT): p max =0.80MPa, rise time t r =1.6ms, decay constant k d =1.1×10 3 s -1 hysteresis p s =0.12MPa.

[0058] Operating Condition B (12kg TNT): p max =1.05MPa, rise time t r =1.8ms, decay constant k d =1.0×10 3 s -1 hysteresis p s =0.15MPa.

[0059] Operating Condition C (15kg TNT): pmax =1.20MPa, rise time t r =2.0ms, decay constant k d =9.0×10 2 s -1 hysteresis p s =0.18MPa.

[0060] Spatial distribution: The inner surface of the vault is the peak area, which decreases in a cosine manner along the circumference to ±30°.

[0061] Numerical Implementation and Solution Control: Time integration: explicit central difference, initial step size Δt0 = 0.5 μs Based on F, the energy error is adaptively reduced to 0.1. μs ; Return mapping: DP yield surface tangent consistency algorithm; Stabilization: Energy error threshold <3%, total kinetic energy drops back to below 10% of the peak value within 20ms after peak value.

[0062] Calculation results and damage evolution: Structural response (conditions A / B / C); Peak acceleration (inner edge of the dome): A / B / C approximately 92 g / 118 g / 136 g respectively; Equivalent plastic strain ε p,eq Peak values ​​(arched waist range of 25°-35°): A=0.35%, B=0.62%, C=0.91%; The direction of the principal tensile stress and the crack zone: initially along the circumferential direction, then connecting longitudinally, forming a channelization trend of "arch waist - arch foot" (significant in B and C).

[0063] The damage variable distribution is D(x,t); Condition A: Local formation of arched waist D A sprouting zone of approximately 0.18 was not penetrated. Condition B: The arch waist and sidewalls form a continuous structure. D ∈[0.25,0.45] extended band, local D max ≈0.58; Condition C: The arched waist shows a tendency to penetrate, peak value D max ≈0.74, the circumferential crack zone and the longitudinal connection zone form an "L"-shaped damage core.

[0064] Rate and viscosity contribution; Compared to "no rate amplification / no viscosity": Peak load capacity increased by approximately 12-18%, and the energy consumption (viscous + damage) ratio increased from 0.31 to 0.44. The softening initiation is delayed, the damage bandwidth increases by about 20-35%, and the numerical localization is suppressed, resulting in a more consistent outcome with the measured damage morphology.

[0065] Grading assessment and engineering judgment By threshold D 1 = 0.20 D 2 = 0.50, D 3 = 0.70 (with residual bearing capacity) or R (Joint judgment) Operating Condition A: Most areas D < D 1. Some hot topics D ≈0.18, or R ≈0.90; Determined "repairable", monitoring + surface sealing recommended.

[0066] Working Condition B: Arch Waist - Side Wall D ∈[ D 1, D 2 and locally Dmax≈0.58 (> D 2), or R ≈0.72; determined "reinforcement is required".

[0067] Condition C: Formation of a near-continuous zone. D max ≈0.74(> D 3), or R ≈0.63; judged as "major reinforcement or partial replacement".

[0068] The above description is merely an embodiment of the present invention and does not limit the patent scope of the present invention. Any equivalent structural or procedural transformations made based on the content of the present invention specification and drawings, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of the present invention.

Claims

1. A method for assessing the damage evolution of underground structures under sudden explosive loads, characterized in that, Includes the following steps: S1: Based on explosive loads, geometric boundaries, and material parameters, simulate underground structures and construct a unified constitutive model; S2: The unified constitutive model introduces the Drucker-Prager criterion pressure-related strength criterion to characterize the influence of confining pressure on shear strength; S3: The unified constitutive model achieves the coupled evolution of high strain rate-damage-yield through rate effect and viscous flow mechanism; S4: The unified constitutive model introduces the Weiber statistical framework to characterize the strength dispersion and crack initiation probability caused by material defects; S5: The unified constitutive model adopts a softening regularization strategy through fracture energy constraint and characteristic length control to suppress mesh dependence and non-physical localization, ensuring that energy dissipation and crack width in the softening stage are mesh-independent, thereby improving the repeatability and robustness of the simulation. S6: The unified constitutive model after softening regularization is further processed using numerical integration and algorithmic framework; S7: Establish a parameter calibration and inversion process based on S6; S8: Provide damage grading based on S6; S2 includes the following steps: S21: Damage equivalent stress; （1） in: For the damage equivalent stress tensor, D As a damage variable, For total strain tensor and It is the plastic strain tensor; S22: Drucker-Prager criterion for yielding and hardening degradation; （2） in: I 1 is the first invariant. J 2 represents the second deviatoric stress invariant; For intensity rate correction; F This represents a Drucker–Prager type yield function; S3 includes the following steps: S31: Intensity rate amplification adopts an example power law or exponential type; The example power law is: Exponential type is: （3） Wherein: parameters (a,b) are fitted by experiments; denoted as equivalent strain rate; a is the strain rate sensitivity coefficient; b is the strain rate exponent parameter. Equivalent plastic strain; S32: Overstressed viscoplasticity; （4） in: For the plastic strain rate tensor, The derivative of the yield function, These are components of the stress tensor; S4 includes the following steps: S41: Probability of crack initiation; （5） S42: Damage Growth Law; （6） in: λ and β For calibration coefficients, F 0 is the normalization constant; D For damage variables; Optional tensile / compressive splitting damage D t and D c Driven by the principal tensile / compressive strain respectively, satisfying: D = χD t +(1 χ ) D c （7） Where: χ is the tension / compression weighting factor; S5 includes the following steps: To eliminate mesh dependence, a fracture energy-controlled softening modulus correction is employed: in: σ pk Peak intensity; according to l ch Automatically adjust the softening slope to ensure conservation of unit fracture energy; G f Values G f t or G f s ; S7 includes the following steps: S71: Data source; Static / quasi-static: Single-axis and triaxial (E, ν, c, Φ) and peak / residual values ​​are obtained. High strain rate: Rate sensitivity curves and fracture energy were obtained from split Hopkinson bar tension / compression and shear bar tests; S72: Inversion objective function: Where: parameter vector ; J (p) is the objective function / cost function used for parameter inversion; w i In order to be with the first i Weighting coefficients corresponding to groups of data; For numerically predicted stress vectors / sequences; For the first i Measured stress vectors / sequences from the group of tests; w g The weights of the fracture energy constraint terms; The fracture energy / energy dissipation per unit area is calculated by the numerical model under parameter p. The fracture energy was measured experimentally. In response to the situation; t For time; S73 Algorithm: Global and local hybrid optimization, cross-validation to prevent overfitting; output parameter confidence intervals.

2. The method for assessing the damage evolution of underground structures under sudden explosive loads as described in claim 1, characterized in that: S1 includes the following steps: S11: Explosion Load History: Input Near / Mid / Far Field Equivalent Pressure – Time History p (t), or import the wall pressure distribution obtained from computational fluid dynamics / fluid-structure interaction calculations. p (x, t), allowing a given rise time t r Peak p max hysteresis ps With decay constant k d ; S12: Geometry and Boundaries: Establish a 3D underground structure—surrounding rock / soil model, and set contact and constraint / radial boundaries; S13: Materials include cement-based, rock-based, or composite materials: Elasticity: Elasticity model E∈[15,60]GPa, Poisson's ratio ν ∈[0.15,0.25]; Strength: The effect of friction angle Φ on yield criterion α∈[0.05,0.35], a material constant related to cohesion c. k ∈[1,20] MPa, or converted from (c, Φ); Statistical damage: shape parameter m∈[2,20], scale parameter ∈[5×10⁻⁴, 5×10⁻³]; Viscosity / Overstress: Coefficient of Viscosity η ∈[10 -3 10 -1 ]s, exponent n∈[1,6], threshold σ0∈[0.1,5]MPa; Strain rate amplification: Strain rate sensitivity coefficient C ∈[0.5,50]s⁻¹, exponential parameter / curve shape parameter p ∈[0.2,1.2]; Fracture energy: tensile crack G f t ∈[50,250]N / m, shear G f s ∈[100,500]N / m; Feature length: Feature length l ch Values ​​are automatically retrieved based on the grid scale.

3. The method for assessing the damage evolution of underground structures under sudden explosive loads as described in claim 1, characterized in that: S6 includes the following steps: Time integration: using explicit central difference or implicit incremental Newton method; Return mapping: Tangentially consistent plastic correction is applied to the yield surface of the Drucker-Prag criterion; rate amplification and viscous flow are coupled and updated in the sub-iteration layer; Substep size adaptive: based on F , With energy error control time step; Contact / Interface: Combining the normal penalty function with Coulomb friction, an interfacial bond-slip-softening segment is introduced to describe the degradation of the lining-soil interface.

4. The method for assessing the damage evolution of underground structures under sudden explosive loads according to claim 1, characterized in that: S8 includes the following steps: Let the threshold be 0 < D 1< D 2< D 3≤1: D ∈[0, D 1): Micro-damage; [ D 1, D 2: Moderate injury; [ D 2, D 3) Severe injury; ≥ D 3: Penetration / Instability; Combined with residual bearing rate η R = R res / R 0. Comprehensive rating based on crack connectivity and functional indicators; Recommended treatment methods: thickening lining, circumferential / longitudinal fiber composite reinforcement, spraying backfill, interface reinforcement, and pressure relief groove.

5. The method for assessing the damage evolution of underground structures under sudden explosive loads according to claim 1, characterized in that: It also includes S9: Result Output and Report, which outputs the results and reports of S7 and S8.

Images