Simulation method and system for concrete gravity dam damage assessment under blast load effect

Abstract

The invention discloses a simulation method for evaluating damage of a concrete gravity dam under the action of an explosion load. The simulation method comprises the following steps: establishing geometric models of the concrete gravity dam, water and explosives, performing grid division on the geometric models, and dispersing units into particles to obtain a discretization model; giving material parameters, boundary conditions and contact to the discretization model, and setting a detonation point position and a detonation mode; coupling effects between the explosive and water, between the water and the concrete gravity dam and between the explosive and the concrete gravity dam are set; key information is initialized; a fluid-structure interaction force is calculated by adopting an asynchronous solution IHPD-SPH coupling method; an SPH solver is adopted to update the speed and density information of the fluid; an IHPD solver is adopted to update internal force and acceleration information of the solid; checking whether the calculation duration is met or not, and if yes, ending calculation; if not, recalculating the fluid-solid coupling force according to the currently obtained particle state; and obtaining a dam damage result, and performing damage grade evaluation.

Description

Technical Field

[0001] The invention belongs to the field of hydraulic structure protection, and in particular relates to a simulation method and system for damage assessment of a concrete gravity dam under explosion load. Background Art

[0002] In the field of water conservancy engineering, concrete gravity dams, as important infrastructure, often face a variety of environmental threats, among which underwater explosion loads are one of the most complex and destructive effects. The pressure waves caused by underwater explosions will have a strong impact on the concrete gravity dam, causing local damage or even overall structural instability, seriously threatening the safety and functionality of the dam. Therefore, the damage assessment of concrete gravity dams under underwater explosion loads has important engineering significance and application value.

[0003] At present, the damage assessment of dams under underwater explosion loads is mainly carried out through three methods: experiments, theoretical analysis and numerical simulation. Due to the high cost of experiments and the limitations of existing technologies, some key data under underwater explosions are difficult to obtain, so experiments are not an ideal assessment method. Theoretical analysis has oversimplified the actual situation and can only solve simple working conditions, and often deviates greatly from the actual results. Existing mesh numerical simulation methods have limitations in describing problems involving large deformations, discontinuities, etc. Additional numerical techniques are needed to solve the distortion of results caused by mesh distortion and it is difficult to solve discontinuities such as crack propagation. Therefore, we urgently need a new simulation method to evaluate the damage of concrete gravity dams under explosion loads in order to predict the protection capacity of concrete gravity dams. Summary of the invention

[0004] In view of the shortcomings of the prior art, the present invention proposes a simulation method and system for damage assessment of concrete gravity dams under explosive loads. The asynchronous solution technology is used to solve the solid domain and the fluid domain respectively, which greatly improves the calculation efficiency and accuracy, so as to predict the safety protection capability of the concrete gravity dam faster and more accurately.

[0005] To achieve the above object, the present invention provides the following solutions:

[0006] A simulation method for damage assessment of a concrete gravity dam under blast loads comprises the following steps:

[0007] S1: Establish the geometric model of concrete gravity dam, water and explosives, mesh the geometric model, discretize the units into particles, and obtain the discretized model;

[0008] S2: Assign material parameters, boundary conditions and contacts to the discretized model, and set the detonation point and detonation method;

[0009] S3: setting the coupling between explosives and water, water and concrete gravity dam, and explosives and concrete gravity dam;

[0010] S4: Initialize key information;

[0011] S5: The fluid-solid coupling force is calculated using the asynchronous solution IHPD-SPH coupling method;

[0012] S6: Use SPH solver to update the velocity and density information of the fluid;

[0013] S7: Use the IHPD solver to update the internal force and acceleration information of the solid;

[0014] S8: Check whether the calculation time is met. If so, the calculation is terminated; if not, S5 is recalculated based on the current particle state.

[0015] S9: Obtain the dam damage results and conduct damage level assessment.

[0016] Preferably, in S2, the method of assigning boundary conditions and contacts to the discretized model includes:

[0017] The bottom of the dam is set as a fixed boundary, the two sides of the dam are set as fixed boundary conditions, and a solid wall boundary is added to the water supply. Different dam sections are independent of each other, and there is friction contact between each dam section.

[0018] Preferably, in said S3, the method of setting the coupling effect between explosives and water, water and concrete gravity dam, and explosives and concrete gravity dam comprises:

[0019] The coupling method between explosive and water is SPH-SPH coupling;

[0020] The coupling method between water and concrete gravity dam, and between explosives and concrete gravity dam is IHPD-SPH coupling.

[0021] Preferably, in S4, the method for initializing key information includes:

[0022] Considering the three-phase components of concrete, there are six types of "bonds", including three "component bonds", i.e., two material points are in the same phase region, and three "interface bonds", i.e., two material points are in different phase regions. The six "bonds" are defined as aggregate-aggregate "bond", mortar-mortar "bond", interface-interface "bond", aggregate-mortar "bond", aggregate-interface "bond", and mortar-interface "bond".

[0023] Preferably, in S5, the method of calculating the fluid-solid coupling force using the asynchronous solution IHPD-SPH coupling method includes:

[0024] While performing one fluid time integration, j solid time integrations are performed. Assume that the time steps of the fluid and solid are Δt1 and Δt2 respectively. Δt1 and Δt2 should satisfy the following conditions:

[0025]

[0026] Where h is the smooth length, c0 is the speed of sound in water, |v| max is the maximum flow rate of the fluid, is the maximum acceleration of the fluid, δ is the radius of the near field, c′ is the speed of sound of the concrete dam, and the time step ratio j is defined as:

[0027]

[0028] The explicit Newmark method is used for time integration, and the field variable information at time n+1 is calculated based on the field variable information at time n, including:

[0029] Predict the speed at time n+1 based on the speed and acceleration at time n And displacement increment Δx:

[0030]

[0031] In the formula, is the speed at time n, is the speed at time n+1, is the acceleration at time n, is the acceleration at time n+1, β and α are variable parameters, and the time step Δt must satisfy the Courant-Friedrichs-Levy condition:

[0032]

[0033] Where c is the speed of sound, C max Generally, it is taken as 0.1;

[0034] Calculate the SPH quantity of virtual particles simulated by IHPD particles, including: the calculation of virtual particle velocity needs to consider the influence of IHPD particle velocity, smooth velocity for:

[0035]

[0036] Speed b* for:

[0037]

[0038] In the formula, v a is the velocity of the neighboring particle a of the virtual particle b*, W is the kernel function, V ais the volume of particle a, is the velocity of IHPD particle b;

[0039] When the SPH particle approaches the interface, a force F is exerted from the IHPD particle b to the particle a. btoa for:

[0040]

[0041] In the formula, m a is the mass of particle a, m b* is the mass of particle b*, σ a is the stress of particle a, σ b* is the stress of particle b*, ρ a is the density of particle a, ρ b* is the density of particle b*, is the gradient of the kernel function, Ω a is the support domain of a;

[0042] Repulsive force f c It should also be applied to particle a, and the expression is:

[0043]

[0044] In the formula, x a is the initial coordinate of particle a, K c and n c is a custom parameter, K c Similar to the contact stiffness in the FEM contact algorithm, Δp avg is the average value of the smooth length between particles; x ab* =x a -x b* is the relative position vector of the particle in the current configuration;

[0045] The force of IHPD particle b on SPH particle a is:

[0046] F PDtoSPHa =Σ b [F btoa -f c (x a )]

[0047] The force of SPH particle a on IHPD particle b is:

[0048]

[0049] Preferably, in S6, the method of updating the velocity and density information of the fluid using the SPH solver includes:

[0050] Update the density of the explosive and water particles using the continuity equation for fluid motion:

[0051]

[0052] Where D / Dt represents the material derivative, ρ a is the density of particle a, m b is the mass of particle b, v ab is the relative velocity of particles a and b, v ab =v a -v b , v a , v b are the speeds of particle a and particle b respectively, is the gradient of the kernel function;

[0053] Update the velocities of the explosive and water particles using the momentum equation for fluid motion:

[0054]

[0055] In the formula, ρ b is the density of particle b, σ a , σ b are the stress of particle a and particle b, respectively, a For physical strength.

[0056] Preferably, in S7, the method of using the IHPD solver to update the internal force and acceleration information of the solid includes:

[0057] The internal force and acceleration of the dam particles are updated through the equations in IHPD, and the equation of motion is:

[0058]

[0059] Ω a ={x a ∈Ω:||x a -x b ||≤δ}

[0060] Where ρ is the density, is the acceleration of particle a, x a is the initial coordinate of particle a, x b is the initial coordinate of particle b, t is the time, V b is the volume of particle b within the near field of particle a, f a is the volume force at particle a, δ is the size of the near field, Ω a is the near field range of particle a, T [x a ,t] <x b -x a> is the force vector state of particle b acting on particle a at time t, T [x b ,t] <x a -x b > is the force vector state of particle a acting on particle b at time t. The angle brackets indicate the bonds on which the state acts. a -x b || is the distance between particle a and particle b;

[0061] Preset force vector state T and force scalar state t Satisfy between:

[0062]

[0063] In the formula, M is a unit vector pointing from the deformed position of particle a to the deformed position of particle b, t is a scalar state, and the calculation formula is:

[0064]

[0065] In the formula, ω is the influence function; p is the peridynamic pressure scalar function; k is the bulk modulus; x is the reference position scalar state; m is the weighted scalar state, expressed as m = ( ωx ) x ; · represents the inner product symbol, α is the equivalent material parameter, expressed as μ is the shear modulus of the material; e d is the deviatoric state of the stretched scalar state, expressed as θ is the weighted expansion scalar state, expressed as e is the stretched scalar state, expressed as e = y - x , y is the deformed scalar state.

[0066] Preferably, in S9, the method of obtaining the dam damage result and conducting damage level assessment includes:

[0067] According to the peridynamic damage model, the critical elongation criterion is used to calculate the damage of concrete dams, including:

[0068] The elongation s of the "key" is:

[0069]

[0070] In the formula, || Y <xb -x a >|| is the distance between particle a and particle b after deformation, || X <x b -x a >|| is the distance between particle a and particle b before deformation;

[0071] The breaking of the “bond” is characterized by introducing a history-dependent Boolean function:

[0072]

[0073] In the formula, ξ is the "bond" between the particle pairs, s0 is the critical elongation of the "bond" between the particle pairs, and when s ≥ s0, it means that the bond is broken;

[0074] The damage at the particle a is defined by the ratio of the number of broken "bonds" to the total number of "bonds" in the near field of the particle a, which is expressed as:

[0075]

[0076] In the formula, It means that there is no damage at particle a. Indicates that particle a is completely destroyed.

[0077] Compared with the prior art, the present invention has the following beneficial effects:

[0078] The present invention mainly adopts the asynchronous solution of subhomogeneous peridynamics-smoothed particle hydrodynamics method (hereinafter referred to as subhomogeneous peridynamics for short IHPD, smoothed particle hydrodynamics for short SPH) coupling method to carry out damage assessment of concrete gravity dam under underwater explosion load. Among them, the IHPD method has significant advantages in simulating the dynamic crack extension of gravity dam and the heterogeneity of concrete, and the SPH method can avoid the grid distortion phenomenon, and has potential advantages in simulating the propagation of explosion gas and explosion waves in water and air. The IHPD-SPH coupling method has the advantages of both the IHPD method and the SPH method, and the underwater explosion is a fully coupled process. Therefore, the IHPD-SPH coupling method is a good means to evaluate the damage of concrete dams under explosion loads. In the simulation process, since the time steps of concrete dam particles and water particles are inconsistent, if a unified time step is used, the total calculation time will be too long, which will lead to low calculation efficiency. Therefore, the present invention adopts asynchronous solution technology to solve the solid domain and the fluid domain respectively, which greatly improves the calculation efficiency and accuracy, so as to predict the safety protection capability of the concrete gravity dam faster and more accurately. BRIEF DESCRIPTION OF THE DRAWINGS

[0079] In order to more clearly illustrate the technical solution of the present invention, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0080] Figure 1 It is a schematic flow chart of a simulation method for damage assessment of a concrete gravity dam under explosion load according to an embodiment of the present invention;

[0081] Figure 2 Schematic diagram of different "key" types of the concrete sub-homogeneous model according to an embodiment of the present invention;

[0082] Figure 3 is a flow chart of a preprocessing algorithm of a subhomogeneous peridynamic model according to an embodiment of the present invention;

[0083] Figure 4 is a schematic diagram of a numerical model of an embodiment of the present invention;

[0084] Figure 5 This is a total damage diagram of the dam according to an embodiment of the present invention. DETAILED DESCRIPTION

[0085] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0086] In order to make the above-mentioned objects, features and advantages of the present invention more obvious and easy to understand, the present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0087] Embodiment 1

[0088] like Figure 1 As shown, the present invention provides a simulation method for damage assessment of concrete gravity dam under explosion load, and stipulates the assessment standard of damage level, including the following steps:

[0089] S1: Establish the geometric model of concrete gravity dam, water and explosives, mesh the geometric model, discretize the units into particles, and obtain the discretized model;

[0090] S2: Assign material parameters, boundary conditions and contacts to the discretized model, and set the detonation point and detonation method;

[0091] S3: setting the coupling between explosives and water, water and concrete gravity dam, and explosives and concrete gravity dam;

[0092] S4: Initialize key information;

[0093] S5: The fluid-solid coupling force is calculated using the asynchronous solution IHPD-SPH coupling method;

[0094] S6: Use SPH solver to update the velocity and density information of the fluid;

[0095] S7: Use the IHPD solver to update the internal force and acceleration information of the solid;

[0096] S8: Check whether the calculation time is met. If so, the calculation is terminated; if not, S5 is recalculated based on the current particle state.

[0097] S9: Obtain the dam damage results and conduct damage level assessment.

[0098] The specific implementation plan is as follows:

[0099] In this embodiment, according to the given size conditions of concrete gravity dam, water, and explosives, their original-size geometric models are respectively established, and the geometric models are meshed, and the units are discretized into particles to obtain discretized concrete gravity dam models, water models, and explosives models.

[0100] In this embodiment, the geometric model is assigned material parameters, boundary conditions and contacts, which involve:

[0101] (2.1) Assigning material parameters to the geometric model;

[0102] (2.2) Assign boundary conditions and contacts, set the bottom of the dam as a fixed boundary, set the two sides of the dam as fixed boundary conditions, and add a solid wall boundary to the water supply. Different dam sections are independent of each other, and there is friction contact between each dam section.

[0103] (2.3) Set the location of the explosive detonation point and the detonation method.

[0104] In this embodiment, coupling effects between explosives and water, water and concrete gravity dam, and explosives and concrete gravity dam are set, and the following are involved:

[0105] (3.1) The coupling method between explosive and water is SPH-SPH coupling;

[0106] (3.2) The coupling method between water and concrete gravity dam, and between explosives and concrete gravity dam is IHPD-SPH coupling.

[0107] In this embodiment, the key information is initialized, and the following is involved:

[0108] Considering the heterogeneity of concrete, a subhomogeneous peridynamic model is used to study the dynamic fracture of concrete. This model considers the aggregate, mortar and interface transition zone inside the concrete, and links the phase components of concrete at the mesoscopic scale with the fracture behavior at the macroscopic scale at the "bond" level without increasing the computational effort in the calculation process.

[0109] The sub-homogeneous peridynamic model works at the "bond" level, that is, by defining different types of "bonds" to characterize the various phase components of concrete at the microscopic scale, but different types of "bonds" are still calculated isotropically during the calculation process, and anisotropy is not considered at the "bond" level. Figure 2 As shown in Figure 1, concrete is considered to be composed of three phase components (aggregate, mortar matrix, and interface transition zone). Any "bond" connecting two material points can be located in the aggregate zone, mortar zone, or interface transition zone. The two material points can be in the same phase region or in different phase regions, so the type of "bond" is defined according to the location of the material points. The probability of different "bonds" occurring depends only on the volume fraction of each phase component.

[0110] Considering the three-phase components of concrete, there are six types of "bonds", including three "component bonds" (two material points are in the same phase region) and three "interface bonds" (two material points are in different phase regions). The six "bonds" are defined as aggregate-aggregate "bond", mortar-mortar "bond", interface-interface "bond", aggregate-mortar "bond", aggregate-interface "bond", and mortar-interface "bond".

[0111] It is assumed that concrete material is homogeneous at the macroscopic scale, so the volume fraction of each phase component is fixed and the proportion of each type of "bond" depends only on the volume fraction of each phase component. Once the type of "bond" is determined, the material parameters of the corresponding component (such as density, elastic modulus, fracture energy, etc.) can be assigned to the corresponding "bond". For the three types of "component bonds", aggregate-aggregate "bond", mortar-mortar "bond" and interface-interface "bond" correspond to the material properties of aggregate, mortar and interface transition zone respectively. Considering that there are three types of "interface bonds", it is assumed that the mechanical properties of aggregate-mortar "bond" are the same as those of the homogeneous concrete model at the macroscopic scale. Aggregate-interface "bond" and mortar-interface "bond" are connected by aggregate, mortar and interface transition zone respectively, so their mechanical properties are assumed to be the same as those of interface-interface "bond".

[0112] After different types of "bonds" are determined, it is necessary to further determine the proportion of each type of "bond". This chapter adopts a simplified method, assuming that the proportion of each "bond" is linearly related to the volume fraction of each phase. If the volume fraction of each phase remains unchanged, the proportion of each "bond" will also be fixed. For example, the volume fraction of aggregate, the volume fraction of mortar, and the volume fraction of the interface transition zone are a, b, and c, respectively. The proportions of the three "component bonds" are aa (aggregate-aggregate "bond"), bb (mortar-mortar "bond"), and cc (interface-interface "bond"), and the proportions of the three "interface bonds" are 2*ab (aggregate-mortar "bond"), 2*ac (aggregate-interface "bond"), and 2*bc (mortar-interface "bond"). As mentioned above, it is assumed that the mechanical properties of aggregate-interface "bonds" and mortar-interface "bonds" are the same as those of interface-interface "bonds", so there are four types of "bonds" in the end, with proportions of aa, bb, cc+2*(ac+bc), and 2*ab. The determination of different types of "bonds" in the subhomogeneous model is implemented as a preprocessing process. The flowchart of the preprocessing process is shown in Figure 3 .

[0113] In this embodiment, the fluid-solid coupling force is calculated using the asynchronous solution IHPD-SPH coupling method, which involves the following:

[0114] During the coupling process, the type of coupled virtual particles corresponding to each IHPD particle is first determined according to the initial configuration of the solid. If the IHPD particle is located on the surface of the structure, it will be disguised as a boundary virtual particle, and other IHPD particles will be disguised as internal virtual particles. Then, at the initial moment of each time step, the position information of the current IHPD particle is passed to the coupled virtual particle. When the distance between the SPH fluid particle and the coupled virtual particle is less than twice the smooth length (2h), the particles come into contact. At this time, it is considered that the SPH fluid particle and the coupled virtual particle are both located in the fluid-solid coupling area. The coupled virtual particle will participate in the field variable information required for SPH calculation, such as density and velocity.

[0115] (5.1) Usually, the time step of solid calculation is smaller than that of fluid calculation. If a unified time step is used, the total calculation time will be too long, which will lead to low calculation efficiency. Therefore, asynchronous solution technology is used. In order to ensure the consistency of the total calculation time of fluid and solid, j solid time integrations are performed while performing one fluid time integration. Assume that the time steps of fluid and solid are Δt1 and Δt2 respectively. Δt1 and Δt2 should satisfy:

[0116]

[0117] In the above formula, h is the smooth length, c0 is the speed of sound in water, |v| max is the maximum flow rate of the fluid, is the maximum acceleration of the fluid, δ is the radius of the near field, c′ is the speed of sound of the concrete dam, and the time step ratio j is defined as:

[0118]

[0119] The explicit Newmark method is used for time integration. The basic idea is to calculate the field variable information at time n+1 based on the field variable information at time n, including:

[0120] Predict the speed at time n+1 based on the speed and acceleration at time n And displacement increment Δx:

[0121]

[0122]

[0123] In the above formula is the speed at time n, is the speed at time n+1, is the acceleration at time n, is the acceleration at time n+1, β and α are variable parameters. To ensure the stability of the calculation process, the time step Δt must satisfy the Courant-Friedrichs-Levy condition:

[0124]

[0125] In the above formula, c is the speed of sound, C max Generally, it is taken as 0.1.

[0126] (5.2) In the IHPD-SPH coupling model, if the SPH particle is located at the interface between the fluid and the solid, the IHPD particle in its support domain is disguised as an SPH virtual particle to achieve SPH boundary defect correction and apply the velocity boundary to the SPH particle. For example, the IHPD particle b is used as a virtual particle b* to participate in the approximation of the SPH particle a, and a repulsive force is applied to the SPH particle when it is very close to the IHPD particle to prevent the occurrence of particle penetration. First, the SPH quantities of the virtual particles simulated by the IHPD particles, such as density and velocity, are calculated. Taking the no-slip boundary condition as an example, the calculation of the virtual particle velocity needs to consider the influence of the IHPD particle velocity, the smooth velocity for:

[0127]

[0128] Speed b* for:

[0129]

[0130] In the above formula, va is the velocity of the neighboring particle a of the virtual particle b*, W is the kernel function, V a is the volume of particle a, is the velocity of IHPD particle b.

[0131] When the SPH particle approaches the interface, a force F is exerted from the IHPD particle b to the particle a. btoa for:

[0132]

[0133] In the above formula, m a is the mass of particle a, m b* is the mass of particle b*, σ a is the stress of particle a, σ b* is the stress of particle b*, ρ a is the density of particle a, ρ b* is the density of particle b*, is the gradient of the kernel function, Ω a is the support domain of a.

[0134] At this time, the repulsive force f c It should also be applied to particle a, and its expression is:

[0135]

[0136] In the above formula, K c and n c is a custom parameter, K c Similar to the contact stiffness in the FEM contact algorithm, Δp avg is the average smooth length between particles, x ab* =x a -x b* is the relative position vector of the particle in the current configuration. Therefore, the force of IHPD particle b on SPH particle a is:

[0137] F PDtoSPHa =∑ b [F btoa -f c (x a )]

[0138] The force of SPH particle a on IHPD particle b is:

[0139]

[0140] In this embodiment, the SPH solver is used to update the velocity, density and other information of the explosive particles and water particles, which involves:

[0141] (6.1) Update the density of explosive particles and water particles through the continuity equation of fluid motion:

[0142]

[0143] In the above formula, D / Dt represents the material derivative, ρ a is the density of particle a, m b is the mass of particle b, v ab is the relative velocity of particles a and b, v ab =v a -v b , v a , v b are the speeds of particle a and particle b respectively, is the gradient of the kernel function.

[0144] (6.2) Update the velocity of explosive particles and water particles through the momentum equation of fluid motion:

[0145]

[0146] In the above formula, ρ b is the density of particle b, σ a , σ b are the stress of particle a and particle b, respectively, a For physical strength.

[0147] In this embodiment, the IHPD solver is used to update the internal force, acceleration and other information of the dam particles, which involves the following:

[0148] The internal force and acceleration of the dam particles are updated through the equations in IHPD. The equation of motion is:

[0149]

[0150] Ω a ={x a ∈Ω:||x a -x b ||≤δ}

[0151] In the above formula, ρ is the density, is the acceleration of particle a, x a is the initial coordinate of particle a, x b is the initial coordinate of particle b, t is the time, V b is the volume of particle b within the near field of particle a, f a is the volume force at particle a, δ is the size of the near field, Ω a is the near field range of particle a, T [x a ,t] <x b -xa > is the force vector state of particle b acting on particle a at time t, T [x b ,t] <x a -x b > is the force vector state of particle a acting on particle b at time t. The angle brackets indicate the bonds on which the state acts. a -x b || is the distance between particle a and particle b.

[0152] Defining the force vector state T and force scalar state t The following conditions are met:

[0153]

[0154] In the above formula M is a unit vector pointing from the deformed position of particle a to the deformed position of particle b, t is a scalar state and can be calculated by the following formula:

[0155]

[0156] In the above formula ω is the influence function; p is the peridynamic pressure scalar function; k is the bulk modulus; x is the reference position scalar state; m is the weighted scalar state, expressed as m = ( ωx ) x ; · represents the inner product symbol, α is the equivalent material parameter, expressed as μ is the shear modulus of the material; e d is the deviatoric state of the stretched scalar state, expressed as θ is the weighted expansion scalar state, expressed as e is the stretched scalar state, expressed as e = y - x , y is the deformed scalar state.

[0157] It is worth noting that there are six types of "bonds" in concrete, namely aggregate-aggregate "bond", mortar-mortar "bond", interface-interface "bond", aggregate-mortar "bond", aggregate-interface "bond", and mortar-interface "bond". These six types of bonds need to be calculated using the material parameters of the corresponding components (such as density, elastic modulus, fracture energy, etc.).

[0158] In this embodiment, the present invention uses the critical elongation criterion to calculate the damage of the concrete dam based on the peridynamic damage model. The elongation s of the "key" is:

[0159]

[0160] In the above formula, || Y <x b -x a >|| is the distance between particle a and particle b after deformation, || X <x b -x a >|| is the distance between particle a and particle b before deformation.

[0161] The breaking of the “bond” is characterized by introducing a history-dependent Boolean function:

[0162]

[0163] In the above formula, ξ is the "bond" between the particle pairs, and s0 is the critical elongation of the "bond" between the particle pairs. When s ≥ s0, it means that the bond is broken.

[0164] The damage at the particle a is defined by the ratio of the number of broken "bonds" to the total number of "bonds" in the near field of the particle a, which is expressed as:

[0165]

[0166] In the above formula, It means that there is no damage at particle a. Indicates that particle a is completely destroyed.

[0167] The total damage range is defined as 0 to 1, where 0 to 0.25 is minor damage, 0.25 to 0.5 is mild damage, 0.5 to 0.75 is moderate damage, and 0.75 to 1 is severe damage. Different colors represent different damage levels so that the damage distribution of the dam can be observed more directly and its weak points can be understood. When the damage reaches the intermediate level or above, necessary protective measures should be taken to improve the overall explosion resistance of the concrete dam.

[0168] Embodiment 2

[0169] In order to verify the effect of the method proposed in the present invention, this case provides a test example.

[0170] 1. Establish the geometric model of concrete gravity dam, water and explosives, mesh the geometric model, discretize the units into particles, and obtain the discretized model. The following involves:

[0171] (1.1) Set as Figure 4The model shown in the figure is composed of three identical dam sections. The model is placed in an aluminum model box of 1200mm×720mm×850mm (length×width×height). Each dam section is 235mm wide and the total width is 705mm. The model dam is 550mm high, 50mm wide at the top and 375mm wide at the bottom. The upstream face of the dam is vertical, and the downstream slope is 1:0.75. The dam foundation is a rectangular parallelepiped, 100mm high and 600mm long.

[0172] (1.2) The model is meshed, water, explosives, and gravity dam are discretized into particles, and the particle size is set to 6 mm.

[0173] 2. Assign material parameters, boundary conditions and contacts to the geometric model, including the following steps:

[0174] (2.1) Determine the material parameters of water, TNT explosives, and concrete dam. The water depth is 600 mm, and the Mie-Gruneisen state equation is used for water:

[0175]

[0176] In the above formula, p is pressure, ρ0 is the initial density of water, C is the speed of sound, γ0 is the Gruneisen parameter, a is the volume correction parameter, S1, S2, S3 are fixed parameters, e is the initial specific internal energy, μ is the volume strain and μ = ρ / ρ0-1, ρ is the current density of water. For water, the specific parameter value is: ρ0 = 1000 kg / m 3 , C=1480m / s, γ0=0.5, a=0, S1=2.56, S2=1.986, S3=1.2268.

[0177] TNT explosive equivalent W = 2.2g, density ρ TNT =1630kg / m 3 , the explosive state equation adopts the JWL state equation:

[0178]

[0179] In the above formula, p is the explosion pressure, η = ρ / ρ0, ρ0 is the initial density of the explosive, A, B, R1, R2, ω are fixed parameters, and e is the detonation specific internal energy. For TNT, the specific parameter value is: ρ0 = 1630 kg / m 3 ,

[0180] A=3.712×10 11 Pa, B = 3.21 × 10 9 Pa, R1=4.15, R2=0.95, ω=0.3, e=4.29×10 6 J / kg.

[0181] (2.2) Assign boundary conditions and contacts, set the bottom of the dam as a fixed boundary, and add a solid wall boundary to the water. Different dam sections are independent of each other, and there is friction contact between each dam section.

[0182] (2.3) Set the detonation mode to point detonation, the blasting distance L = 100mm, and the blasting depth R = 20mm.

[0183] 3. Setting the coupling between explosives and water, water and concrete gravity dam, and explosives and concrete gravity dam involves the following:

[0184] (3.1) The coupling method between explosive and water is SPH-SPH coupling;

[0185] (3.2) The coupling method between water and concrete gravity dam, and between explosives and concrete gravity dam is IHPD-SPH coupling.

[0186] 4. Press as follows Figure 3 Flowchart initialization key information:

[0187] 5. The fluid-solid coupling force is calculated using the asynchronous IHPD-SPH coupling method, which involves the following:

[0188] (5.1) During the simulation, since the time steps of the concrete dam particles and the water particles are inconsistent, if a unified time step is used, the total calculation time will be too long, which will lead to low calculation efficiency. Therefore, the asynchronous solution technology is used. In order to ensure the consistency of the total calculation time of the fluid and solid, j solid time integrations are performed while performing one fluid time integration. Assume that the time steps of the fluid and solid are Δt1 and Δt2 respectively. Δt1 and Δt2 should satisfy:

[0189]

[0190] In the above formula, h is the smooth length, c0 is the reference sound speed, |v| max is the maximum flow rate of the fluid, is the maximum acceleration of the fluid, δ is the radius of the near field, c′ is the sound velocity of the solid material, and the time step ratio j is defined as:

[0191]

[0192] The explicit Newmark method is used for time integration. The basic idea is to calculate the field variable information at time n+1 based on the field variable information at time n, including:

[0193] Predict the speed at time n+1 based on the speed and acceleration at time n And displacement increment Δx:

[0194]

[0195]

[0196] In the above formula is the speed at time n, is the speed at time n+1, is the acceleration at time n, is the acceleration at time n+1, β and α are variable parameters. To ensure the stability of the calculation process, the time step Δt must satisfy the Courant-Friedrichs-Levy condition:

[0197]

[0198] In the above formula, c is the speed of sound, and C max The value generally does not exceed 1.0.

[0199] (5.2) In the IHPD-SPH coupling model, if the SPH particle is located at the interface between the fluid and the solid, the IHPD particle in its support domain is disguised as an SPH virtual particle to achieve SPH boundary defect correction and apply the velocity boundary to the SPH particle. For example, the IHPD particle b is used as a virtual particle b* to participate in the approximation of the SPH particle a, and a repulsive force is applied to the SPH particle when it is very close to the IHPD particle to prevent the occurrence of particle penetration. First, the SPH quantities of the virtual particles simulated by the IHPD particles, such as density and velocity, are calculated. Taking the no-slip boundary condition as an example, the calculation of the virtual particle velocity needs to consider the influence of the IHPD particle velocity, the smooth velocity for:

[0200]

[0201] Speed b* for:

[0202]

[0203] In the above formula, v a is the velocity of the neighboring particle a of the virtual particle b*, W is the kernel function, V a is the volume of particle a, is the velocity of IHPD particle b.

[0204] When the SPH particle approaches the interface, a force F is exerted from the IHPD particle b to the particle a. btoa for:

[0205]

[0206] In the above formula, m a is the mass of particle a, m b*is the mass of particle b*, σ a is the stress of particle a, σ b* is the stress of particle b*, ρ a is the density of particle a, ρ b* is the density of particle b*, is the gradient of the kernel function, Ω a is the support domain of a.

[0207] At this time, the repulsive force f c It should also be applied to particle a, and its expression is:

[0208]

[0209] In the above formula, K c and n c is a custom parameter, K c Similar to the contact stiffness in the FEM contact algorithm, Δp avg is the average smooth length between particles, x ab* =x a -x b* is the relative position vector of the particle in the current configuration. Therefore, the force of IHPD particle b on SPH particle a is:

[0210] F PDtoSPHa =∑ b [F btoa -f c (x a )]

[0211] The force of SPH particle a on IHPD particle b is:

[0212]

[0213] 6. Use the SPH solver to update the velocity, density and other information of explosive particles and water particles, which involves:

[0214] (6.1) Update the density of explosive particles and water particles through the continuity equation of fluid motion:

[0215]

[0216] In the above formula, D / Dt represents the material derivative, ρ a is the density of particle a, m b is the mass of particle b, v ab is the relative velocity of particles a and b, v ab =v a -v b , v a , v b are the speeds of particle a and particle b respectively, is the gradient of the kernel function.

[0217] (6.2) Update the velocity of explosive particles and water particles through the momentum equation of fluid motion:

[0218]

[0219] In the above formula, ρ b is the density of particle b, σ a , σ b are the stress of particle a and particle b, respectively, a For physical strength.

[0220] 7. Use the IHPD solver to update the internal force, acceleration and other information of the dam particles, which involves the following:

[0221] The internal force and acceleration of the dam particles are updated through the equations in IHPD. The equation of motion is:

[0222]

[0223] Ω a ={x a ∈Ω:||x a -x b ||≤δ}

[0224] In the above formula, ρ is the density, is the acceleration of particle a, x a is the initial coordinate of particle a, x b is the initial coordinate of particle b, t is the time, V b is the volume of particle b within the near field of particle a, f a is the volume force at particle a, δ is the size of the near field, Ω a is the near field range of particle a, T [x a ,t] <x b -x a > is the force vector state of particle b acting on particle a at time t, T [x b ,t] <x a -x b > is the force vector state of particle a acting on particle b at time t. The angle brackets indicate the bonds on which the state acts. a -x b || is the distance between particle a and particle b.

[0225] Defining the force vector state T and force scalar state t The following conditions are met:

[0226] T = t M

[0227] In the above formula M Right now M ( Y ) <x a -x b > is a unit vector pointing from the deformed position of particle a to the deformed position of particle b, t is a scalar state and can be calculated by the following formula:

[0228]

[0229] In the above formula ω is the influence function; p is the peridynamic pressure scalar function; k is the bulk modulus; x is the reference position scalar state; m is the weighted scalar state, expressed as m = ( ωx ) x ; · represents the inner product symbol, α is the equivalent material parameter, expressed as μ is the shear modulus of the material; e d is the deviatoric state of the stretched scalar state, expressed as θ is the weighted expansion scalar state, expressed as e is the stretched scalar state, expressed as e = y - x , y is the deformed scalar state.

[0230] It is worth noting that there are six types of "bonds" in concrete, namely aggregate-aggregate "bond", mortar-mortar "bond", interface-interface "bond", aggregate-mortar "bond", aggregate-interface "bond", and mortar-interface "bond". These six bonds need to be calculated using the material parameters of the corresponding components, as shown in Table 1.

[0231] Table 1 Material parameters of each phase component of concrete

[0232]

[0233] 8. Check whether the calculation time is met. If it is met, the calculation is terminated. If not, recalculate step 5 based on the current particle state.

[0234] 9. The present invention uses the critical elongation criterion to calculate the damage of concrete dams based on the peridynamic damage model. The elongation s of the "key" is:

[0235]

[0236] In the above formula, || Y <x b -x a >|| is the distance between particle a and particle b after deformation, || X <x b -x a >|| is the distance between particle a and particle b before deformation.

[0237] The breaking of the “bond” is characterized by introducing a history-dependent Boolean function:

[0238]

[0239] In the above formula, ξ is the "bond" between the particle pairs, and s0 is the critical elongation of the "bond" between the particle pairs. When s ≥ s0, it means that the bond is broken.

[0240] The damage at the particle a is defined by the ratio of the number of broken "bonds" to the total number of "bonds" in the near field of the particle a, which is expressed as:

[0241]

[0242] In the above formula, It means that there is no damage at particle a. Indicates that particle a is completely destroyed.

[0243] The total damage range is defined as 0 to 1, where 0 to 0.25 is minor damage, 0.25 to 0.5 is mild damage, 0.5 to 0.75 is moderate damage, and 0.75 to 1 is severe damage. Different colors are used to represent different damage levels, so that the damage distribution of the dam can be observed more directly, its weak points can be understood, and more effective protective measures can be taken.

[0244] Depend on Figure 5 It can be seen that the damage level at the dam head and the variable section is high and the damage is serious. Therefore, these are the weak parts of the concrete dam. It is necessary to focus on strengthening protection to improve the overall explosion resistance of the concrete dam.

[0245] The simulation method proposed in the present invention fully demonstrates the process of detonation, blast wave propagation, and bubble pulsation under underwater explosion loads, solves the grid distortion problem existing in existing numerical simulation methods when fluid-solid coupling and large deformation are involved, and improves the solution accuracy and efficiency. At the same time, the present invention also formulates a set of standards for judging the degree of damage of various parts of concrete, so as to more intuitively understand its weak parts and take effective protective measures.

[0246] The embodiments described above are only descriptions of the preferred embodiments of the present invention and are not intended to limit the scope of the present invention. Without departing from the design spirit of the present invention, various modifications and improvements made to the technical solutions of the present invention by ordinary technicians in this field should all fall within the protection scope determined by the claims of the present invention.

Claims

1. A simulation method for damage assessment of concrete gravity dams under blast loads, characterized in that: The following steps are involved: S1: Establish the geometric model of concrete gravity dam, water and explosives, mesh the geometric model, discretize the units into particles, and obtain the discretized model; S2: Assign material parameters, boundary conditions and contacts to the discretized model, and set the detonation point and detonation method; S3: setting the coupling between explosives and water, water and concrete gravity dam, and explosives and concrete gravity dam; S4: Initialize key information; S5: The fluid-solid coupling force is calculated using the asynchronous solution IHPD-SPH coupling method; S6: Use SPH solver to update the velocity and density information of the fluid; S7: Use the IHPD solver to update the internal force and acceleration information of the solid; S8: Check whether the calculation time is met. If so, the calculation is terminated; if not, S5 is recalculated based on the current particle state. S9: Obtain the dam damage results and conduct damage level assessment.

2. The simulation method for damage assessment of concrete gravity dam under explosion load according to claim 1 is characterized in that: In S2, the method of assigning boundary conditions and contacts to the discretized model includes: The bottom of the dam is set as a fixed boundary, the two sides of the dam are set as fixed boundary conditions, and a solid wall boundary is added to the water supply. Different dam sections are independent of each other, and there is friction contact between each dam section.

3. The simulation method for damage assessment of concrete gravity dam under explosion load according to claim 1 is characterized in that: In S3, the method of setting the coupling effect between explosives and water, water and concrete gravity dam, and explosives and concrete gravity dam includes: The coupling method between explosive and water is SPH-SPH coupling; The coupling method between water and concrete gravity dam, and between explosives and concrete gravity dam is IHPD-SPH coupling.

4. The simulation method for damage assessment of concrete gravity dam under explosion load according to claim 1, characterized in that: In S4, the method for initializing key information includes: Considering the three-phase components of concrete, there are six types of "bonds", including three "component bonds", i.e., two material points are in the same phase region, and three "interface bonds", i.e., two material points are in different phase regions. The six "bonds" are defined as aggregate-aggregate "bond", mortar-mortar "bond", interface-interface "bond", aggregate-mortar "bond", aggregate-interface "bond", and mortar-interface "bond".

5. The simulation method for damage assessment of concrete gravity dam under explosion load according to claim 1, characterized in that: In S5, the method for calculating the fluid-solid coupling force using the asynchronously solved IHPD-SPH coupling method includes: While performing one fluid time integration, j solid time integrations are performed. Assume that the time steps of the fluid and solid are Δt1 and Δt2 respectively. Δt1 and Δt2 should satisfy the following conditions: Where h is the smooth length, c0 is the speed of sound in water, |v| max is the maximum flow rate of the fluid, is the maximum acceleration of the fluid, δ is the radius of the near field, c′ is the speed of sound of the concrete dam, and the time step ratio j is defined as: The explicit Newmark method is used for time integration, and the field variable information at time n+1 is calculated based on the field variable information at time n, including: Predict the speed at time n+1 based on the speed and acceleration at time n And displacement increment Δx: In the formula, is the speed at time n, is the speed at time n+1, is the acceleration at time n, is the acceleration at time n+1, β and α are variable parameters, and the time step Δt must satisfy the Courant-Friedrichs-Levy condition: Where c is the speed of sound, C max Generally, it is taken as 0.1; Calculate the SPH quantity of virtual particles simulated by IHPD particles, including: the calculation of virtual particle velocity needs to consider the influence of IHPD particle velocity, smooth velocity for: speed for: In the formula, v a is the velocity of the neighboring particle a of the virtual particle b*, W is the kernel function, V a is the volume of particle a, is the velocity of IHPD particle b; When the SPH particle approaches the interface, a force F is exerted from the IHPD particle b to the particle a. btoa for: In the formula, m a is the mass of particle a, is the mass of particle b*, σ a is the stress of particle a, is the stress of particle b*, ρ a is the density of particle a, is the density of particle b*, is the gradient of the kernel function, Ω a is the support domain of a; Repulsive force f c It should also be applied to particle a, and the expression is: In the formula, x a is the initial coordinate of particle a, K c and n c is a custom parameter, K c Similar to the contact stiffness in the FEM contact algorithm, Δp avg is the average value of the smooth length between particles; is the relative position vector of the particle in the current configuration; The force of IHPD particle b on SPH particle a is: F PDtoSPHa =∑ b [F btoa -f c (x a )] The force of SPH particle a on IHPD particle b is:

6. The simulation method for damage assessment of concrete gravity dam under explosion load according to claim 1, characterized in that: In S6, the method of updating the velocity and density information of the fluid using the SPH solver includes: Update the density of the explosive and water particles using the continuity equation for fluid motion: Where D / Dt represents the material derivative, ρ a is the density of particle a, m b is the mass of particle b, v ab is the relative velocity of particles a and b, v ab =v a -v b , v a , v b are the speeds of particle a and particle b respectively, is the gradient of the kernel function; Update the velocities of the explosive and water particles using the momentum equation for fluid motion: In the formula, ρ b is the density of particle b, σ a , σ b are the stress of particle a and particle b, respectively, a For physical strength.

7. The simulation method for damage assessment of concrete gravity dam under explosion load according to claim 1, characterized in that: In S7, the method of using the IHPD solver to update the internal force and acceleration information of the solid includes: The internal force and acceleration of the dam particles are updated through the equations in IHPD, and the equation of motion is: Ω a ={x a ∈Ω:||x a -x b ||≤δ} Where ρ is the density, is the acceleration of particle a, x a is the initial coordinate of particle a, x b is the initial coordinate of particle b, t is the time, V b is the volume of particle b within the near field of particle a, f a is the volume force at particle a, δ is the size of the near field, Ω a is the near field range of particle a, T [x a ,t] <x b -x a > is the force vector state of particle b acting on particle a at time t, T [x b ,t] <x a -x b > is the force vector state of particle a acting on particle b at time t. The angle brackets indicate the bonds on which the state acts. ||x a -x b || is the distance between particle a and particle b; Preset force vector state T and force scalar state t Satisfy between: T ＝ T In the formula, M is a unit vector pointing from the deformed position of particle a to the deformed position of particle b, t is a scalar state, and the calculation formula is: In the formula, ω is the influence function; p is the peridynamic pressure scalar function; k is the bulk modulus; x is the reference position scalar state; m is the weighted scalar state, expressed as m = ( ωx ) x ; · represents the inner product symbol, α is the equivalent material parameter, expressed as μ is the shear modulus of the material; e d is the deviatoric state of the stretched scalar state, expressed as θ is the weighted expansion scalar state, expressed as e is the stretched scalar state, expressed as e = y - x , y is the deformed scalar state.

8. The simulation method for damage assessment of concrete gravity dam under explosion load according to claim 7, characterized in that: In S9, the method of obtaining the dam damage result and conducting damage level assessment includes: According to the peridynamic damage model, the critical elongation criterion is used to calculate the damage of concrete dams, including: The elongation s of the "key" is: In the formula, || Y <x b -x a >|| is the distance between particle a and particle b after deformation, || X <x b -x a >|| is the distance between particle a and particle b before deformation; The breaking of the "bond" is characterized by introducing a history-dependent Boolean function: In the formula, ξ is the "bond" between the particle pairs, s0 is the critical elongation of the "bond" between the particle pairs, and when s≥s0, it means that the bond is broken; The damage at the particle a is defined by the ratio of the number of broken "bonds" to the total number of "bonds" in the near field of the particle a, expressed as: In the formula, It means that there is no damage at particle a. Indicates that particle a is completely destroyed.

Images