Underground structure damage evolution evaluation method under sudden explosion load effect

Abstract

The invention discloses a method for evaluating damage evolution of an underground structure under the action of a sudden explosion load, and relates to the technical field of civil engineering and engineering disaster prevention, and the method is characterized by comprising the following steps: S1, simulating the underground structure based on an explosion load, a geometric boundary and material parameters, and constructing a unified constitutive model; and S2, introducing a pressure-related strength criterion of a Dwark-Prasugrel criterion into the unified constitutive model so as to describe the influence of confining pressure on the shear strength. The technical problem to be solved by the invention is to provide the method for evaluating the damage evolution of the underground structure under the action of the burst explosion load, so that the key pain points such as pressure correlation, high strain rate effect, material random discreteness, damage softening energy consistency and the like are difficult to consider at the same time in the existing analysis method. A set of damage evolution evaluation method and implementation system which can be implemented from the material level to the structure level is established, and quantitative prediction and engineering grading judgment of the whole process of crack initiation, expansion to instability failure of the underground structure are supported.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of civil engineering and engineering disaster prevention technology, in particular to a method for evaluating damage evolution of underground structures under sudden explosion load. BACKGROUND

[0002] With the large-scale development of urban underground space, the number of underground structures such as subway stations and tunnels, underground comprehensive pipe galleries, municipal and hydraulic caverns, and national defense projects has rapidly increased. These projects are in a deep overburden and surrounding rock / soil constraint environment. Once subjected to sudden explosion load (such as explosions caused by transportation, storage, construction, or terrorist attacks), they will exhibit dynamic characteristics of high confining pressure, multi-axial stress, and strong rate effect: the shock wave generated by the explosion will be reflected and superimposed multiple times in a closed or semi-closed space, accompanied by quasi-static lag pressure of gas products, which is easy to induce tensile / shear composite cracking and rapid damage propagation at stress concentration sites such as arch crown, arch waist, and corners, thereby causing load bearing capacity degradation and connectivity damage, leading to secondary disasters and operational interruptions. Therefore, there is an urgent need for quantitative damage evolution evaluation and anti-explosion design methods for explosion impact conditions.

[0003] From the perspective of material and mechanical mechanism, the commonly used materials for underground structures (concrete, mortar composite, fiber reinforced cement-based materials, etc.) exhibit significant characteristics under explosion impact: (1) pressure dependence: the increase of confining pressure can enhance the shear strength and ductility, and the yield surface changes significantly with the first stress invariant; (2) high strain rate effect: the strength and apparent peak strain increase with the increase of loading rate; (3) quasi-brittle and damage softening: micro-crack initiation and cooperative propagation lead to stiffness and strength degradation, resulting in softening and localization; (4) statistical dispersion: the internal defects, aggregate distribution, and interface characteristics of the material are not uniform, resulting in significant randomness in the crack initiation strain, damage threshold, and failure path. Accordingly, accurately characterizing the coupling behavior of "pressure dependence-rate sensitivity-statistical dispersion-damage softening-energy dissipation" is a key challenge for the analysis of underground structures under explosion conditions.

[0004] The existing technologies mainly include the following methods, each of which has obvious shortcomings: 1. Empirical and equivalent methods Relying on scaled distance and empirical shock response spectrum, equivalent single degree of freedom (SDOF) method, etc., the displacement or equivalent internal force of the component can be quickly estimated, which is suitable for early scheme comparison and selection. However, this method is based on free propagation of shock waves on the ground or in open space, and cannot reflect the multiple reflections, confining pressure effect, and lag pressure action in underground environment; it also cannot reveal the damage evolution process of crack initiation, propagation, and penetration, and has limited guidance for damage classification and repair decision-making.

[0005] 2. Traditional elastic-plastic models (such as Mohr-Coulomb, Drucker-Prager, and capped models) Such 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; only static parameters are fitted at high strain rates, which may underestimate peak strength and energy dissipation capacity; and the localization and penetration of cracks are not well described, making it impossible to give a consistent damage field with the actual failure pattern.

[0006] 3. Continuous damage and cracking model Isotropic / anisotropic damage, tension-compression split damage, and concrete damage plasticity (CDP) can describe strength and stiffness degradation, and some models also support tension-compression asymmetry. However, common problems include the following steps: (1) the rate effect is simply handled by a calibration factor that independently modifies the strength without coupling with damage evolution and yield surface; (2) there is a lack of statistical description of material discreteness and weakened elements, resulting in sensitivity of crack initiation location and failure path to mesh and initial defects; (3) the mesh dependence and non-physical energy dissipation caused by softening are not systematically solved, resulting in unstable damage zone width and energy dissipation with mesh changes.

[0007] 4. Specialized anti-blast concrete models (such as the high-pressure high-rate constitutive family commonly used in engineering) Such models consider confining pressure, strain rate, and fragmentation evolution, and have good applicability to penetration, penetration, and near-field impact problems. However, their parameter system is complex and the calibration cost is high, often relying on specific test platforms and large-sample impact tests; and most focus on material unit-level responses, lacking engineering outputs that link to damage classification and repairability assessment at the structure level. Some implementations still have problems with regularization of softening and numerical stability.

[0008] 5. Discrete / lattice and phase-field fracture methods Discrete elements, lattice / micro-beam networks, and phase-field fracture methods can naturally generate crack networks and penetration paths, and have strong description capabilities for fracture topology. However, the parameter and scale mapping relationship is complex and the computational cost is high; their application is limited in large-scale underground structure analysis, and they lack sufficient connection to soil-structure interaction and engineering classification processes.

[0009] 6. Numerical stability and parameter inversion problems The stress-strain curve under explosive impact shows a significant softening segment, and if there is no energy regularization or viscosity / overstress stabilization mechanism constrained by fracture energy, non-physical localization and result divergence often occur. On the other hand, the commonly used triaxial test, split Hopkinson pressure bar (SHPB) test, and other data in engineering cannot establish a high-robustness automatic inversion process with the model parameter system, resulting in subjective parameter selection, poor repeatability, and reduced model engineering promotion value.

[0010] Therefore, the engineering field needs a constitutive and evaluation framework for underground structure explosion impact scenarios, which can uniformly couple pressure-related yield (such as DP), strain rate amplification, viscosity / overstress flow and Weibull statistical damage; at the same time, energy regularization is matched to ensure mesh independence in the softening stage, and a parameterized, calibrated and hierarchical output engineering process and software system is formed to support the whole process prediction, damage level judgment and reinforcement decision from crack initiation, propagation to instability failure. This is the urgent direction and real demand of the development of the field. SUMMARY

[0011] The technical problem to be solved by the present application is to provide an underground structure damage evolution evaluation method under sudden explosion load, which solves the key pain points such as pressure correlation, high strain rate effect, material random discreteness and damage softening energy consistency that cannot be considered simultaneously by existing analysis methods, establishes a damage evolution evaluation method and implementation system that can be implemented from the material to the structure level, and supports quantitative prediction and engineering hierarchical judgment of the whole process from crack initiation, propagation to instability failure of underground structure.

[0012] The present application realizes the purpose of the application by adopting the following technical scheme: An underground structure damage evolution evaluation method under sudden explosion load, characterized in that it comprises the following steps: S1: simulating an underground structure based on explosion load, geometric boundary and material parameters, and constructing a unified constitutive model; S2: introducing the Drucker-Prager criterion pressure-dependent strength criterion into the unified constitutive model to depict the influence of confining pressure on shear strength; S3: realizing the coupling evolution of high strain rate-damage-yield through the rate effect and viscous flow mechanism of the unified constitutive model; S4: introducing the Weibull statistical framework into the unified constitutive model to represent the strength discreteness and crack initiation probability caused by material defects; S5: adopting a softening regularization strategy that is constrained by fracture energy and controlled by characteristic length in the unified constitutive model to suppress grid dependence and non-physical localization, ensure the grid independence of energy dissipation and crack band width in the softening stage, and improve the repeatability and robustness of simulation; S6: further processing the unified constitutive model after softening regularization by using numerical integration and algorithm framework; S7: establishing a parameter calibration and inversion process based on S6; S8: providing damage grading based on S6.

[0013] As a further limitation of the technical scheme, S1 comprises the following steps: S11: explosion load history: input near / middle / far field equivalent pressure-time history p(t), or importing the wall pressure distribution from computational fluid dynamics / fluid-structure interaction calculation p (x, t), allowing a given rise time t r , peak p max , hysteresis ps and decay constant k d ; S12: Geometry and boundary: Establish a three-dimensional underground structure-surrounding rock / soil model, set contact (normal stiffness K n , tangential friction coefficient μ ) and constraint / radiation boundary; S13: Material includes cement-based, rock or composite material: Elasticity: Elastic model E ∈ [15, 60] GPa, Poisson's ratio ν ∈ [0.15, 0.25]; Strength (DP equivalent): Influence of friction angle Φ on yield criterion α ∈ [0.05, 0.35], 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-4, 5 × 10-3]; Viscosity / overstress: Viscosity coefficient η ∈ [10 -3 ,10 -1 ] s, index n ∈ [1, 6], threshold σ0 ∈ [0.1, 5] MPa; Strain rate amplification: Strain rate sensitivity coefficient C ∈ [0.5, 50] s-1, index parameter / curve shape parameter p ∈ [0.2, 1.2]; Fracture energy: Tensile fracture G f t ∈ [50, 250] N / m, shear fracture G f s ∈ [100, 500] N / m; Characteristic length (element-dependent): Characteristic length l ch Automatically selected according to grid size.

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

[0015] As a further limitation of the technical solution, the S3 includes the following steps: S31: Strength rate amplification employs an example power law or exponential type; where the example power law is The exponential type is: (3) where: parameters (a, b) are fitted by experiments; is the equivalent strain rate; a is the strain rate sensitivity coefficient, b is the strain rate index parameter, is the equivalent plastic strain; S32: Over-stress viscoplasticity; (4) where: is the plastic strain rate tensor, is the differential of the yield function, is the component of the stress tensor; As a further limitation of the technical solution, the S4 includes the following steps: S41: Crack initiation probability; (5) S42: Damage growth law; (6) where: λ and β are calibration coefficients, F 0 is a normalization constant; D is the damage variable; Possible tensile / compressive splitting damage D t and D c are driven by the principal tensile / compressive strain, respectively, and satisfy: D =χD t + (1- χ ) D c (7) where χ is the tension-compression weight factor.

[0016] As a further limitation of the technical solution, the S5 comprises the following steps: To eliminate the grid dependence, the softening modulus is modified by the fracture energy control:

[0017] where: σ pk peak strength; according to l ch The softening slope is automatically adjusted to ensure the conservation of unit fracture energy; G f The value is G f t or G f s .

[0018] As a further limitation of the technical solution, the S6 comprises the following steps: Time integration: explicit central difference or implicit incremental Newton method; Return mapping: tangent consistent plastic correction on the Drucker-Prager yield surface; rate amplification and viscous flow are coupled and updated at the sub-iteration level.

[0019] Sub-step self-adaptation: according to F , and energy error control time step; Contact / interface: combined with normal penalty function and Coulomb friction, interface bonding-slippage-softening segment is introduced as necessary to describe the degradation of lining-soil interface.

[0020] As a further limitation of the technical solution, the S7 comprises the following steps: S71: data source; Static / quasi-static: uniaxial and triaxial tests to obtain (E, v, c, Φ) and peak / residual; High strain rate: split Hopkinson pressure bar tension / compression, shear bar test to obtain rate sensitive curve and fracture energy.

[0021] S72: inversion objective function:

[0022] where the parameter vector ; J(p) is the objective function / cost function for parameter inversion (identification); w i is the weight coefficient (dimensionless) corresponding to the i-th i group of data; is the stress vector / sequence of numerical prediction; is the stress vector / sequence of the i-th i group of test; w g is the weight (dimensionless) of the fracture energy (or energy dissipation) constraint term; is the fracture energy / energy dissipation per unit area calculated by the numerical model under the parameter p; is the fracture energy measured by the test; is the strain; t is the time; S73 Algorithm: global and local mixed optimization, cross-validation to prevent overfitting; output parameter confidence interval.

[0023] As a further limitation of the technical solution, the S8 comprises the following steps: Set the threshold value 0 D 1 D 2 D 3≤1: D ∈[0, D 1) : micro-damage; [ D 1, D 2: moderate damage; [ D 2, D 3) : severe damage; ≥ D 3: through / instability; Combined with the residual bearing ratio η R = R res / R 0( η R is the residual bearing ratio, R res is the residual bearing capacity, R 0 is the initial bearing capacity), crack connectivity and functional index comprehensive judgment level.

[0024] Output treatment suggestion: thickening lining, ring / longitudinal fiber composite reinforcement, jet backfill, interface reinforcement, pressure relief groove, etc.

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

[0026] Compared with the prior art, the application has the advantages and positive effects that: The application constructs a unified constitutive model: the Drucker-Prager (DP) pressure-dependent strength criterion is introduced to depict the influence of confining pressure on the shear strength, the Weibull statistical framework is introduced to represent the strength dispersion and cracking probability caused by material defects, and the coupling evolution of high strain rate-damage-yield is realized through the rate amplification and viscous / overstress flow mechanism, so that the description of the strength and stiffness degradation process matched with the explosion impact characteristics is obtained.

[0027] The application realizes numerical stability and energy consistency: through the softening regularization strategy of fracture energy constraint and characteristic length control, the grid dependence and non-physical localization are inhibited, the energy dissipation and crack band width in the softening stage are ensured to be grid-independent, and the repeatability and robustness of the simulation are improved.

[0028] The application establishes a parameter calibration and inversion process: a parameter identification method is formed with conventional triaxial test and high strain rate (such as SHPB) test as data source, statistical damage parameters, rate sensitive parameters and viscous parameters are automatically inverted, and traceable and reusable engineering closed loop from test data to model parameters is realized.

[0029] The application provides engineering classification and decision output: the grading threshold value matched with damage variable, residual bearing and crack penetration connectivity is set, the engineering criteria such as “repairable / need reinforcement / need replacement” are output, and the corresponding structure reinforcement and material optimization suggestions are given, which is convenient for direct application in design, evaluation and emergency disposal.

[0030] The application forms a software and system implementation that can be landed: independent application programs or user subprograms (UMAT / VUMAT) and supporting modules (load generation, material constitutive, solution stabilization, parameter identification, evaluation grading and report visualization) are provided for general finite element platforms to meet the multi-scene needs of scientific research analysis, engineering design and rapid evaluation.

[0031] In summary, through the integrated technical route of “DP pressure-dependent yield-Weibull statistical damage-rate and viscosity coupling-energy regularization-engineering grading output”, the application strives to balance between mechanism integrity, numerical robustness and engineering usability, and provides reliable technical support for underground structure anti-blast design and post-damage evaluation. BRIEF DESCRIPTION OF DRAWINGS

[0032] Figure 1 The method flowchart of the application is shown in the figure; Figure 2 The unified constitutive and yield surface diagram of the application is shown in the figure; Figure 3 The Weibull probability density and damage growth curve of the application is shown in the figure; Figure 4 A schematic diagram of crack band width control under energy regularization of the present application; Figure 5 A typical underground lining structure grid and boundary of the present application; Figure 6 A schematic diagram of damage grading of the present application.

[0033] Figure 7 A schematic diagram of a simulated underground structure of the present application. DETAILED DESCRIPTION

[0034] One specific embodiment of the present application will be described in detail below with reference to the accompanying drawings, but it should be understood that the scope of protection of the present application is not limited by the specific embodiment.

[0035] The present application comprises the following steps: S1: Simulate an underground structure based on an explosion load, a geometric boundary and material parameters, and construct a unified constitutive model; S2: Introduce a Drucker-Prager (DP) pressure-dependent strength criterion to the unified constitutive model to depict the influence of confining pressure on shear strength; S3: The unified constitutive model realizes the coupled evolution of high strain rate-damage-yield through the rate effect and viscous flow mechanism; S4: The unified constitutive model introduces a Weibull statistical framework to represent the strength dispersion and crack initiation probability caused by material defects; S5: The unified constitutive model adopts a softening regularization strategy that is constrained by fracture energy and controlled by characteristic length, suppresses grid dependence and non-physical localization, ensures the grid independence of energy dissipation and crack band width in the softening stage, and improves the repeatability and robustness of the simulation; S6: Further process the unified constitutive model after softening regularization using numerical integration and algorithm framework; S7: Establish a parameter calibration and inversion process based on S6; S8: Provide damage grading based on S6.

[0036] The S1 comprises the following steps: S11: Explosion load history: input near / middle / far-field equivalent pressure-time history p (t), or import wall pressure distribution calculated by computational fluid dynamics / fluid-structure interaction (CFD / FSI) p (x, t) with a given rise time t r , peak p max , lag pressure ps and decay constant kd ; S12: Geometry and boundary: Establish a 3D model of underground structure - surrounding rock / soil, set contact (normal stiffness K n , tangential friction coefficient μ ) and constraint / radiation boundary; S13: Material including cement-based, rock mass or composite material: Elasticity: Elastic model E ∈ [15, 60] GPa, Poisson's ratio ν ∈ [0.15, 0.25]; Strength (DP equivalent): Influence of friction angle Φ on yield criterion α ∈ [0.05, 0.35], 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-4, 5 × 10-3]; Viscosity / overstress: Viscosity coefficient η ∈ [10 -3 , 10 -1 ] s, index n ∈ [1, 6], threshold σ0 ∈ [0.1, 5] MPa; Strain rate amplification: Strain rate sensitive coefficient C ∈ [0.5, 50] s-1, index parameter / curve shape parameter p ∈ [0.2, 1.2]; Fracture energy: Tensile fracture G f t ∈ [50, 250] N / m, shear fracture G f s ∈ [100, 500] N / m; Characteristic length (element-dependent): Characteristic length l ch Automatically taken according to grid scale.

[0037] The S2 comprises the following steps: S21: Damage equivalent stress; (1) Wherein: is the damage equivalent stress tensor, D is the damage variable, is the total strain tensor and is the plastic strain tensor; S22: Drucker-Prager criterion yield and hardening degradation; (2) wherein: I 1 is a first invariant, J 2 is a second deviatoric stress invariant; is a strength rate correction; F denotes a Drucker-Prager type yield function.

[0038] The S3 comprises the following steps: S31: Strength rate amplification employs an example power law or exponential type; wherein the example power law is The exponential type is: (3) wherein: parameters (a, b) are fitted by experiments; is an equivalent strain rate; a is a strain rate sensitivity coefficient, b is a strain rate index parameter, is an equivalent plastic strain; S32: Overstress viscoplasticity; (4) wherein: is a plastic strain rate tensor, is a differential of a yield function, is a component of a stress tensor; is controlled by η , n, σ0, the rate sensitivity and numerical stability are avoided, and non-physical jumps are avoided.

[0039] The S4 comprises the following steps: S41: Crack initiation probability; (5) S42: Damage growth law (probability-overstress coupling); (6) wherein: λ and β are calibration coefficients, F 0 is a normalization constant; D is a damage variable; The tensile / compressive split damage D t and D c are driven by the principal tensile / compressive strain, respectively, and satisfy: D = χD t +(1- χ ) D c (7) wherein: χ is a tensile-compressive weight factor.

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

[0041] 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 .

[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 i measured stress vector / sequence of the group test; w g For the weight (dimensionless) of the fracture energy (or energy dissipation) constraint term; For the fracture energy / energy dissipation per unit area calculated by the numerical model at parameter p; For the fracture energy measured by the test; For the strain; t For the time; α is the yield function parameter, which is related to the friction angle Φ; k For the strength constant, which is related to the cohesion c ; m For the Weibull shape parameter; For the Weibull scale parameter; η For the viscosity coefficient; n For the over-stress exponent; C For the strain rate sensitivity coefficient; p For the strain rate amplification exponent; G f t For the tensile fracture energy; G f s For the shear fracture energy.

[0047] S73 Algorithm: Global (genetic / particle swarm) and local (quasi-Newton / LM) hybrid optimization, cross-validation to prevent overfitting; output parameter confidence interval.

[0048] The S8 comprises the following steps: Set the threshold value 0 D 1 D 2 D 3≤1: D ∈[0, D 1) Micro-damage (repairable); [ D 1, D 2: Moderate damage (local reinforcement is required); [ D 2, D 3) Severe damage (structure bearing is greatly degraded); ≥ D 3: Through / instability (replacement or major reinforcement is required); Combined with the residual bearing rate η R = R res / R 0( η R Residual bearing rate,R res For residual bearing capacity, R 0 is the initial bearing capacity, crack connectivity (based on damage equivalent surface connectivity analysis) and functional index (displacement / crack width limit) comprehensive level.

[0049] Output treatment suggestion: thickening lining, ring / longitudinal fiber composite reinforcement, jet backfill, interface reinforcement, pressure relief groove, etc.

[0050] Also includes S9: result output and report, output the results and report of S7 and S8.

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

[0052] Spatial field: equivalent plastic strain, damage equivalent surface, main crack path, dangerous component sorting.

[0053] Report: automatically generate portable document format / hypertext markup language (PDF / HTML), including parameter table, level map and reinforcement suggestion list, can export neutral data (visualization toolkit / comma separated value, VTK / CSV).

[0054] Embodiment: Engineering background and goal A certain city subway section adopts equal-thickness ring reinforced concrete lining, inner radius R =3.30m, lining thickness t =0.35m, overburden thickness about 12-15m. The design requires to perform anti-blast performance verification under the equivalent impact of "ground near-field explosion equivalent to 8-15kg TNT", and obtain damage classification.

[0055] Geometry, mesh and boundary A three-dimensional 60° periodic sector model (considering ring symmetry) is established, with a length of 6.0m, lining-soil coupling; the lining adopts 8-node solid element, the finest mesh size h min =40mm (local encryption at haunch / springing), the finest h min =80mm of surrounding rock; lining-soil normal penalty function, tangential Coulomb friction coefficient μ =0.65, interface normal stiffness K n =1.0×10 9 N / m 3 ; external domain is set to viscoelastic / radiation boundary, with assumed damping ξ=4%.

[0056] Material and model parameters (according to the technical scheme of the application) 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, η =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.20 MPa, rise time t r = 2.0 ms, decay constant k d = 9.0 x 10 2 s -1 , pressure p s = 0.18 MPa.

[0060] Spatial distribution: peak zone at the inner surface of the vault, cosine-type decay along the hoop direction to ±30°.

[0061] Numerical implementation and control: Time integration: explicit central difference, initial step size At0= 0.5 μs , adaptive reduction to 0.1 based on F, energy error μs ; Return mapping: consistent algorithm for the tangent of the DP yield surface; Stabilization: energy error threshold < 3%, total kinetic energy peak back to below 10% within 20 ms after the peak.

[0062] Computational results and damage evolution: Structural response (case A / B / C); Peak acceleration (inner edge of the vault): A / B / C about 92 g / 118 g / 136 g; Equivalent plastic strain e p,eq peak (range of 25°-35° at the haunch): A = 0.35%, B = 0.62%, C = 0.91%; Principal tensile stress direction and crack band: initially along the hoop direction, then connected to the longitudinal direction, forming a "haunch-keel" channelization trend (B, C significantly).

[0063] Damage variable distribution D(x, t); Case A: local formation of D ≈0.18 initiation zone at the haunch, not through; Case B: coherent D ∈[0.25, 0.45] expansion zone formed at the haunch-wall, locally D max ≈0.58; Case C: through trend at the haunch, peak D max ≈0.74, hoop crack band and longitudinal connection zone forming an "L" shaped damage core.

[0064] Rate and viscous contribution; Comparison with "no rate amplification / no viscosity": The peak load is increased by about 12-18%, and the energy consumption (viscosity + damage) is increased from 0.31 to 0.44; The softening initiation is delayed, the damage bandwidth is increased by about 20-35%, the numerical localization is inhibited, and the results are more consistent with the actual failure mode.

[0065] Grading evaluation and engineering judgment According to the threshold D 1=0.20, D 2=0.50, D 3=0.70 (with residual bearing rate η R Joint judgment): Working condition A: most areas D D 1, individual hot spots D ≈0.18, η R ≈0.90; judgment "repairable", recommend monitoring + surface sealing.

[0066] Working condition B: arch waist - side wall D ∈[ D 1, D 2 and locally reaches Dmax≈0.58(> D 2), η R ≈0.72; judgment "need reinforcement".

[0067] Working condition C: forming a nearly through zone, D max ≈0.74(> D 3), η R ≈0.63; judgment "major reinforcement or local replacement".

[0068] The above is only an embodiment of the present application, and does not limit the patent range of the present application, and any equivalent structure or equivalent process transformation using the content of the specification and drawings of the present application, or directly or indirectly applied in other related technical fields, are also included in the patent protection range of the present application.​

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.

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 according to claim 2, characterized in that: 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 For 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.

4. The method for assessing the damage evolution of underground structures under sudden explosive loads according to claim 3, characterized in that: 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 the components of the stress tensor.

5. The method for assessing the damage evolution of underground structures under sudden explosive loads according to claim 4, characterized in that: 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.

6. The method for assessing the damage evolution of underground structures under sudden explosive loads according to claim 5, characterized in that: 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 .

7. The method for assessing the damage evolution of underground structures under sudden explosive loads according to claim 6, 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, and introducing an interfacial bond-slip-softening segment where necessary to describe the degradation of the lining-soil interface.

8. The method for assessing the damage evolution of underground structures under sudden explosive loads according to claim 7, characterized in that: 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.

9. The method for assessing the damage evolution of underground structures under sudden explosive loads according to claim 7, 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 capacity η 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.

10. 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