Analysis method and system for dynamic fracture of quasi-brittle material

Abstract

The invention provides an analysis method and system for dynamic fracture of a quasi-brittle material. The method is used for solving the technical problem that an existing method is insufficient in dynamic fracture capturing capacity in quasi-brittle materials. The method comprises the following steps: constructing a representative volume element RVE of a quasi-brittle material with two embedded fracture planes; a dual-scale model of the RVE is constructed, the dual-scale model comprises a macroscale sub-model and a mesoscale sub-model, and the crack traction force and the displacement jump increment in the mesoscale sub-model meet the rate-dependent cohesive friction law; identifying the parameters of the dual-scale model; and embedding the dual-scale model into finite element software to simulate the dynamic fracture of the quasi-brittle material. The contained mesoscale crack plane has the key characteristics of the dynamic fracture behavior of the quasi-brittle material, and has remarkable advantages compared with a classical empirical model.

Description

Technical Field

[0001] The present invention relates to the field of material property testing, and particularly to an analysis method and system for dynamic fracture of quasi-brittle materials. Background Art

[0002] The dynamic fracture of quasi-brittle materials caused by factors such as earthquakes, rock explosions, blasting, protective structure design, or meteorite impacts is a common problem in geotechnical engineering and engineering applications. Generally speaking, dynamic fracture is a complex process that is significantly different from quasi-static fracture, including crack initiation, orientation, and propagation due to the effects of high strain rate, inertial effects, and stress wave interactions. At the same time, the dynamic fracture mechanical parameters of quasi-brittle materials, including mode I, II, and III fracture toughness, as well as various strength parameters (i.e., tensile, compressive, shear), exhibit significant but varying degrees of rate dependence. This rate dependence of the mechanical properties of quasi-brittle materials has been extensively studied by experimental methods. When subjected to impact loading, the dynamic fracture process of quasi-brittle materials can be characterized by the propagation of local bands, where obvious discontinuities occur, such as those observed on stressed rocks by high-speed digital image correlation method (DIC) at high strain rates. Given the significant strain rate sensitivity of the fracture process (including initiation and propagation), the rate-dependent cohesive and frictional interactions along the local fracture plane constitute the key control factors for dynamic fracture behavior.

[0003] In recent decades, the constitutive modeling of dynamic fracture of quasi-brittle materials has been widely explored. Traditional phenomenological constitutive models only consider the mechanical response at the macroscopic level while ignoring the underlying local failure mechanisms. As an alternative, the cohesive zone model (CZM) has been more widely applied. This is a continuum-based method that represents the fracture process zone (FPZ) as a local region of damage and softening that evolves along the fracture path. In this region, damage, microcracks, and material degradation occur, and its mechanical response is controlled by the traction-separation law. The traction-separation law provides a basis for the further implementation of rate-dependent characteristics. Although the cohesive zone model (CZM) has significant capabilities and simplicity in simulating dynamic fracture problems in various numerical methods, its application is still limited by various factors: (1) A predefined fracture surface is required (except for the extended finite element method (XFEM)), which is both time-consuming and unrealistic in cases involving complex multiple cracks; (2) The stress locking problem occurs because the fracture direction remains unchanged after crack initiation, resulting in inaccurate stress transfer under non-proportional loading; (3) In the discrete element method (DEM) and phase field model (PFM), the computational efficiency is further limited due to the computational burden of the bonded interface, making its application in large-scale or high-precision simulations less practical. Therefore, due to the inherent limitations of the fixed fracture orientation in the standard cohesive model, the stress locking problem in cohesive elements remains unresolved.

[0004] In addition to the cohesive zone model, significant progress has also been made in other continuum-based rate-dependent constitutive modeling. Existing research has developed a rate-sensitive microplane model that captures the dual rate effects: microcrack propagation kinetics and concrete viscosity. Existing research has proposed a damage plasticity model that combines the effects of confining pressure and strain rate on rock strength. Their model incorporates the DIF formula into the constitutive model through a radial enhancement method to ensure uniform radial enhancement of the yield surface. Although these continuum-based methods can effectively capture the dynamic fracture of brittle materials, they operate under the assumption of uniform material deformation. This inherent limitation reduces their physical significance in dealing with dynamic fracture phenomena and cannot explain significant discontinuities.

[0005] To overcome these limitations, existing research has proposed a two-scale constitutive model that has made significant progress in simulating the fracture process of quasi-brittle materials such as concrete and rock. This model uses a meso-to-macro scale framework to inherently capture the high strain discontinuities on the crack plane through kinematic enrichment. The model initially contains a crack and uses a damage model to describe the inelastic behavior of the material in a local area. Then, it is extended to include two embedded cohesive friction cracks, represented by a damage-plasticity cohesive friction model that naturally provides the initiation and direction of the localization band while effectively alleviating the stress locking problem. The two-scale framework has been proven to be effective in simulating various fracture phenomena of quasi-brittle materials under various mixed-mode loading scenarios. This includes laboratory-scale simulations of the basic load paths in rocks and concretes under quasi-static loads, such as uniaxial tension, triaxial compression, shear, and triaxial tension. However, its ability to capture dynamic fracture in quasi-brittle materials is still insufficient. Summary of the Invention

[0006] The object of the present invention is to provide an analysis method and system for the dynamic fracture of quasi-brittle materials, which is used to solve the technical problem of the insufficient ability of existing methods to capture dynamic fracture in quasi-brittle materials.

[0007] An analysis method for the dynamic fracture of quasi-brittle materials, the specific steps are as follows:

[0008] S1: Construct a representative volume element (RVE) of the quasi-brittle material with two embedded fracture planes;

[0009] S2: Construct a two-scale model of the representative volume element (RVE), the two-scale model includes a macroscale and a mesoscale submodel, and the crack traction force and displacement jump increment in the mesoscale submodel satisfy the rate-dependent cohesive friction law;

[0010] S3: Identify the dual-scale model parameters, which include material property parameters, elastoplastic parameters, and rate-dependence parameters;

[0011] S4: Embed the dual-scale model into finite element software to simulate the dynamic fracture of quasi-brittle materials.

[0012] Optionally, the specific steps for constructing the dual-scale model in step S2 are as follows:

[0013] S2.1: Establish the mathematical relationship between the crack strain increment in the macro-scale sub-model and the crack displacement jump increment in the meso-scale sub-model. The crack strain increment is equal to the product of the crack displacement jump increment and the crack normal vector;

[0014] S2.2: Establish the mathematical relationship between the crack displacement jump increment and the crack traction force in the meso-scale sub-model. The crack traction force and the displacement jump increment in the meso-scale sub-model satisfy the rate-dependent cohesive friction law;

[0015] S2.3: Establish the mathematical model of the crack stress in the macro-scale sub-model. The virtual work generated by the crack stress and the crack strain increment in the macro-scale sub-model should be equal to the total work done by the stress and strain increments inside and outside the fracture surface in the meso-scale sub-model.

[0016] Optionally, the rate-dependent cohesive friction law is a dynamic yield function combining hyperbolic yield and failure:

[0017]

[0018] where t n and t s are the normal and tangential traction forces respectively, D is the damage variable, DIF T and DIF S are the tensile and shear dynamic enhancement factors respectively, A = (1 - D)f t + B, B = c(1 - D) / [2(1 - D)+2Du 2 , μ is the internal friction coefficient, f t and c are the static tensile strength and static cohesion of the material respectively.

[0019] Optionally, the damage variable D is:

[0020]

[0021] where is the plastic displacement of the crack in the meso-scale sub-model, and α and β are the coefficients controlling the influence of normal and shear displacements on damage; and They are the plastic displacements in the normal and shear directions respectively, and δ0 is the parameter for normalizing the relative displacement of u p

[0022] Optionally, in the rate-dependent cohesive friction law, the tensile and shear dynamic enhancement factors DIF T and DIF S are respectively:

[0023]

[0024] wherein, is the tensile strain rate at the mesoscopic scale, is the shear strain rate at the mesoscopic scale.

[0025] Optionally, the material property parameters in step S3 include the static tensile strength f t and the static cohesion c;

[0026] The elastoplastic parameters include the elastic tensile stiffness K n and the elastic shear stiffness K s , and the rate-dependent parameters include the parameters α and β respectively related to the fracture energy, and the relative displacement parameter δ0.

[0027] An analysis system for dynamic fracture of quasi-brittle materials, used to implement the above-mentioned analysis method for dynamic fracture of quasi-brittle materials, includes:

[0028] Element construction module: used to construct a representative volume element RVE of the quasi-brittle material with two embedded fracture planes;

[0029] Model construction module: used to construct a two-scale model of the representative volume element RVE, the two-scale model includes macroscopic scale and mesoscopic scale sub-models, and in the mesoscopic scale sub-model, the crack traction force and the displacement jump increment satisfy the rate-dependent cohesive friction law;

[0030] Parameter identification module: used to identify the two-scale model parameters, and the two-scale model parameters include material property parameters, elastoplastic parameters and rate-dependent parameters;

[0031] Performance analysis module: embed the two-scale model into finite element software to simulate the dynamic fracture of the quasi-brittle material.

[0032] Due to the adoption of the above technical solutions, the present invention has the following advantages:

[0033] ​1. This application extends the rate-dependent fracture characteristics of quasi-brittle materials from the macroscopic scale to the mesoscopic scale for description. By adopting the DIF strain rate empirical law at the mesoscopic scale, a cross-scale correlation is established, and a refined representation of the interaction between dynamic cracks and the external matrix is achieved based on the kinematic enhancement method.

[0034] 2. By incorporating the rate-dependent enhancement effects of tensile strength and shear strength into the cohesive friction law respectively, this application can effectively capture the dynamic fracture response under mixed modes. This characteristic is crucial for its wide application in simulating complex and arbitrary stress loading conditions.

[0035] 3. The two-scale model of this application contains two embedded fracture planes, which effectively solves the stress locking problem. At the same time, the constructed representative volume element RVE has an inherent length scale parameter related to the size of the volume element, thus ensuring numerical convergence during discretization refinement.

[0036] Other advantages, objectives, and features of the present invention will, to some extent, be elaborated in the subsequent specification, and to some extent, will be obvious to those skilled in the art based on the study of the following text, or can be taught from the practice of the present invention. The objectives and other advantages of the present invention can be realized and obtained through the following specification. Brief Description of the Drawings

[0037] The brief description of the drawings of the present invention is as follows.

[0038] Figure 1 It is a flowchart of the analysis method for the dynamic fracture of quasi-brittle materials of the present invention.

[0039] Figure 2 It is a schematic structural diagram of the representative volume element RVE with two embedded fracture planes of the present invention.

[0040] Figure 3 It is a schematic structural diagram of the two-scale model of the present invention.

[0041] Figure 4 It is a graph of the tensile and shear DIF of the quasi-brittle materials of the present invention at different strain rates and the corresponding fitting curves.

[0042] Figure 5 It is an evolution diagram of the rate-dependent yield surface of the present invention.

[0043] Figure 6 It is the specimen geometry for the dynamic direct shear test of the present invention and its strain rate DIF S relationship diagram.

[0044] Figure 7 It is a simulation diagram of the shear stress-displacement curve of the present invention.

[0045] Figure 8 It is a simulation diagram of the rate dependence of the mechanical parameters of the present invention on the shear rate.

[0046] Figure 9 It is a simulation diagram of the dynamic spalling test of the present invention.

[0047] Figure 10 It is a comparison diagram of the dynamic tensile strength and strain rate between the experiment and simulation of the present invention.

[0048] Figure 11 It is an FPZ diagram of two typical crack types of the present invention.

[0049] Figure 12 It is an experimental setup diagram of the mixed-mode three-point bending test of the present invention.

[0050] Figure 13 It is a comparison diagram of the experimental results and numerical results of the mixed-mode three-point bending test.

[0051] Figure 14 It is a comparison diagram of the fracture mode and simulation of the three-point bending test at different notch positions under static and dynamic loads. Detailed implementation manners

[0052] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0053] Embodiment 1:

[0054] As Figure 1 shown, an analysis method for dynamic fracture of a quasi-brittle material, the specific steps are as follows:

[0055] S1: Construct a representative volume element RVE of the quasi-brittle material with two embedded fracture planes;

[0056] In this embodiment, as Figure 2 shown, in the two-scale framework, the representative volume element (RVE) Ω with two embedded fracture planes consists of two fracture process zones (FPZ) Ω k (k = 1, 2) and the external volume Ω0, and the fracture plane is characterized by its thickness l k , area A k and direction θ k , and its normal vector is defined by the vector n k . The nominal size of the representative volume element RVE is defined as H k = Ω / A k .

[0057] S2: Construct a two-scale model of the representative volume element (RVE). The two-scale model includes a macro-scale and a meso-scale sub-model. In the meso-scale sub-model, the crack traction force and the displacement jump increment satisfy the rate-dependent cohesive friction law. The specific steps are as follows:

[0058] S2.1: Establish the mathematical relationship between the crack strain increment in the macro-scale sub-model and the crack displacement jump increment in the meso-scale sub-model. The specific steps are as follows:

[0059] In this embodiment, as Figure 3 shown, dε and are the crack strain increment and the crack strain rate respectively. The macro-strain increment dε of the representative volume element (RVE) consists of the strain increment dε ik of the crack and the external volume dε0. Here, i represents the interior, k represents the crack number. Considering its volume fraction, the expression is as follows:

[0060] dε = (1 - η1 - η2)dε0 + η1dε i1 + η2dε i2 (1)

[0061] In the formula, η k = Ω k / Ω = l k / H k represents the volume fraction of the crack. The fracture surface is simplified as an extremely thin region of a quasi-brittle material (l k << H k ). The strain increment within the FPZ can be approximately calculated using the displacement jump increment at the meso-scale level:

[0062]

[0063] In the formula, du k represents the displacement jump increment at the fracture surface, and n k is the crack normal vector at an angle θ to the horizontal axis;

[0064] By substituting Equation (2) into Equation (1) and performing some transformations, the strain increment of the external block can be expressed as:

[0065]

[0066] S2.2: Establish the mathematical relationship between the crack displacement jump increment and the crack traction force in the meso-scale sub-model. In the meso-scale sub-model, the crack traction force and the displacement jump increment satisfy the rate-dependent cohesive friction law;

[0067] In this embodiment, the dynamic response of the fracture surface is described by the rate-dependent cohesive friction law, where the subscript c represents a variable in the local coordinate system at the mesoscopic scale. When these variables are related to the global coordinate system, the transformation matrix R should be used. The total displacement jump increment vector du c can be decomposed into normal and tangential elastic components and plastic components

[0068]

[0069] In the formula, du c =[du n du s1 du s2 T du n and du s1 、du s2 respectively represent the normal and tangential displacement jump increments on the fracture surface.

[0070] The relationship between the crack traction force t c and the crack displacement can be calculated as:

[0071]

[0072] where, t n 、t s1 、t s2 respectively represent the normal and tangential traction forces on the fracture surface, K n and K s respectively represent the normal and tangential cohesive elastic stiffnesses of the fracture surface; is the secant stiffness matrix; D is the damage variable, and H(t n ) is the Heaviside function used to consider the crack closure effect; and respectively represent the normal and tangential plastic displacement jumps.

[0073] In this embodiment, the ratio of the dynamic strength to the static strength, i.e., the dynamic enhancement factor DIF, is introduced to capture the strength enhancement under dynamic loading. The dynamic enhancement factor DIF includes the tensile dynamic enhancement factor DIF T and the shear dynamic enhancement factor DIF S . The DIF within a wide range of dynamic strain rates introduces the cohesive-friction mechanical law to comprehensively simulate the dynamic crack behavior of the mixed mode. The tensile dynamic enhancement factor DIF T is defined as the ratio of the dynamic tensile strength f t D to the static tensile strength f t , while the shear dynamic enhancement factor DIF S ​is the dynamic shear strength and the ratio of the static shear strength f s . A logarithmic function is used to describe the relationship between DIF and the strain rate:

[0074]

[0075] As Figure 4 shown, it shows the variations of DIF T and DIF S with the tensile and shear strain rates. The rapid increase in the DIF value indicates that both the tensile strength and the shear strength exhibit significant sensitivity to the dynamic strain rate.

[0076] In this embodiment, the rate-dependent cohesion-friction law is a dynamic yield function of a hyperbolic form combination of yield and failure, and damage is introduced as an evolution parameter to represent the accumulation of plastic displacement:

[0077]

[0078] where μ is the coefficient of internal friction, t s 2 = t s1 2 + t s2 2 , A = (1 - D)f t + B, B = c(1 - D) / [2(1 - D)+2Du 2 , f t and c are respectively the static tensile strength and the static cohesion of the material, m is a parameter controlling the shape of the initial yield surface, DIF T and DIF S serve as the proportional parameters in the yield function, and by changing the intersections of the yield surface with the normal traction axis (X-axis) and the shear traction axis (Y-axis), the expansion and contraction of the yield surface are respectively controlled.

[0079] As Figure 5 shown, the hyperbolic form of the yield function is beneficial for mixed-mode fracture because it integrates the normal stress and shear stress components into a unified and continuous formula, thus allowing smooth coupling and transition between pure tension, pure shear, and mixed-mode conditions. When DIF = 1, the yield function represents the static yield surface. As the strain rate on the fracture plane enters the dynamic range, DIF T and DIF S increase logarithmically, resulting in a proportional expansion of the yield surface to represent the dynamic strength enhancement. Specifically, DIF T acts on the scaling of the dynamic tensile strength, while DIF SScaling for dynamic shear strength, with each parameter adjusted based on the tensile and shear strain rates on the fracture plane. This arrangement ensures that the model can accurately and independently capture the dynamic mechanical responses in the normal and shear directions. Note that due to the time-varying impact load, before fracture, DIF T and DIF S will vary dynamically according to the real-time strain rate until the material begins to yield. After yielding, the two parameters remain unchanged during the post-failure evolution.

[0080] In this embodiment, when the material is subjected to a dynamic load, the normal and shear strain rates are calculated based on the most critical potential fracture plane at the mesoscopic scale. At this time, the static yield surface (corresponding to DIF T = DIF S = 1) will expand to the dynamic yield surface (stage 1 shown in Figure 5 ), characterizing the dynamic strength enhancement effect. Once the traction state reaches the dynamic yield surface, damage begins to accumulate, driving the yield surface to gradually contract until it reaches the final failure surface represented by the dashed line (stage 2 shown in Figure 5 ), and this failure surface corresponds to D = 1 and follows the classical frictional Mohr-Coulomb criterion. The traction-displacement evolution is jointly controlled by the dynamic yield surface, the flow rule, and the damage evolution. In stage 1, the expansion of the dynamic yield surface is controlled by DIF T and DIF S ; while in stage 2, the contraction of the yield surface is dominated by the damage variable D. The plastic potential function and the non-associated flow rule are defined as:

[0081]

[0082] where γ is the parameter controlling non-associativity; dλ is the plastic multiplier.

[0083] To complete the cohesive force model, the traction-separation law is formulated based on the coupling of damage and plasticity. Given that irreversible deformation and material degradation occur simultaneously in quasi-brittle geotechnical materials, a damage variable D is defined, which is the result of the propagation of plastic displacement jumps in the normal and shear directions when the traction stress state reaches the yield surface. In addition, since the mechanical response of rock or concrete after failure usually shows a rapid decrease in stress after yielding and then gradually decreases under continuous loading, an exponential-form damage evolution function is introduced to capture this behavior:

[0084]

[0085] where is the cumulative plastic displacement in dimensionless form, and α and β are the coefficients controlling the contributions of normal displacement and shear displacement to damage, and They are the plastic displacements in the normal and shear directions respectively, and δ0 is the relative displacement used to normalize u p whose magnitude is determined by the material-specific plastic displacement characteristics. δ0 provides the necessary flexibility in simulating various material behaviors from brittle to ductile responses.

[0086] S2.3: Construct a mathematical model of crack stress in the macro-scale sub-model, and the specific steps are as follows:

[0087] In this embodiment, in order to relate the mechanical responses inside and outside the fracture surface to the overall behavior of the representative volume element (RVE), the Hill-Mandel condition is adopted, that is, the virtual work generated by the macroscopic stress σ and the strain increment dε should be equal to the total work done by the stress and strain increment inside and outside the fracture surface at the mesoscopic scale:

[0088]

[0089] where σ ik and σ0 are the stress tensors inside and outside the fracture surface respectively.

[0090] By substituting equations (1) and (2) into equation (11) and considering the following equation can be obtained:

[0091]

[0092] where t ik is the traction force on the fracture surface. In order to make the crack strain and crack displacement increment satisfy the above Hill-Mandel condition at any value, the following condition should be met:

[0093]

[0094] where a0 is the elastic stiffness matrix of the material. Assuming that the thickness of the fracture surface is extremely thin (η k →0), equation (3) can be expressed as:

[0095]

[0096] The traction-displacement relationship of the fracture surface in the global coordinate system is in the form of:

[0097]

[0098] where is the tangent stiffness matrix of the fracture surface. Substituting equation (15) into the incremental form of equation (13), we get:

[0099]

[0100] Substituting Equation (14) into Equation (16), the strain-displacement relationship can be obtained as follows:

[0101]

[0102] where M i (i = 1, 4) is a 2×2 matrix in two dimensions and a 3×3 matrix in three dimensions. The displacement of the fracture surface is obtained through the given macroscopic strain increment:

[0103]

[0104] where V i (i = 1, 2) is a 4×4 matrix in two dimensions and a 6×6 matrix in three dimensions. Specifically, in two dimensions:

[0105]

[0106] By substituting Equation (14) and Equation (18) into the incremental form of Equation (13), the constitutive relationship between the macroscopic stress and strain can be obtained:

[0107]

[0108] Equation (21) provides an explicit form of the macroscopic stress-strain relationship for the two-scale constitutive model embedding two fracture surfaces. It can be seen that the macroscopic behavior is jointly controlled by the responses of the fracture surface and the external block. In addition, by introducing the characteristic length H k of the fracture surface, the model intrinsically considers the size effect at the constitutive level, so as to be able to capture the size-dependent behavior of the material.

[0109] In this embodiment, the two-scale constitutive model code can be written into the user subroutine VUMAT of the ABAQUS finite element software to accurately simulate the dynamic fracture of quasi-brittle materials.

[0110] S3: Identify the two-scale model parameters, where the two-scale model parameters include material property parameters, elastoplastic parameters, and rate-dependent parameters;

[0111] The material property parameters include the static tensile strength f t , the static compressive strength f c , the Young's modulus E, the Poisson's ratio ν, the friction coefficient and the parameter m for controlling the shape of the static yield surface.

[0112] In this embodiment, the material property parameters are obtained through basic mechanical tests (such as pure tension tests and uniaxial compression tests). The parameter m ensures that its intersection with the Y-axis represents the initial static cohesion. By obtaining the experimental data sets (such as yield points and failure directions) of triaxial compression tests or direct shear tests, m can be better calibrated. Existing research shows that Young's modulus, Poisson's ratio, and friction angle are rate-independent parameters.

[0113] The elastoplastic parameters include the elastic tensile stiffness K n and the elastic shear stiffness K s , and the elastic tensile stiffness K n and the elastic shear stiffness K s are calibrated through the linear segments of the stress-displacement curves obtained from pure tension tests and direct shear tests respectively.

[0114] The rate-dependent parameters include the parameters α and β related to the mode-I and mode-II fracture energies (G Ι and G ΙΙ ) respectively, and the displacement parameter δ0 used to normalize the plastic displacement u p and regulate the damage evolution, as well as the dilatancy parameter γ used to control the ratio between the plastic normal displacement and the shear displacement.

[0115] In this embodiment, the parameter α is calibrated by ensuring that the projected area of the stress-displacement curve under pure tension loading on the displacement axis is equal to G Ι . Similarly, β is calibrated by matching it with G ΙΙ through shear tests. Larger values of α or β may lead to more rapid softening behaviors in tension and shear, corresponding to smaller G Ι and G ΙΙ respectively. The dilatancy parameter γ is calibrated based on the dilatancy behavior. In this embodiment, the calibration of the rate-dependent parameters focuses on capturing the rate-dependence of the dynamic tensile strength and shear strength, and these characteristics are represented by the empirical formulas DIF T and DIF S , which vary with the tensile strain rate and shear strain rate respectively. For experiments with sufficient dynamic strength-strain rate data, by utilizing their linear relationship on the logarithmic scale and fitting the curve to the experimental data points, the empirical coefficients of DIF can be easily calibrated. For experiments with limited data, the coefficients are calibrated through simulation iterations to make the strain rate and dynamic structural strength consistent with the experimental results. Overall, the effectiveness of the calibrated coefficients is verified by comparing the simulated predicted dynamic crack paths and structural mechanical responses with the experimental observations.

[0116] S4: Used to embed the two-scale model into finite element software to simulate the dynamic fracture of quasi-brittle materials.

[0117] In this embodiment, the finite element software selects ABAQUS, and embeds the two-scale model into the user-defined explicit subroutine VUMAT of the finite element software ABAQUS. Specifically, for any dynamic impact problem in an actual situation, a model is built in ABAQUS. This numerical model contains multiple representative volume elements (RVEs), and each RVE responds according to the two-scale model (the two-scale constitutive model written into VUMAT is selected during operation). Finally, the dynamic fracture of the structure under the impact is reflected. Specifically, in the built model, each RVE will be subjected to different stress and strain effects, and each RVE will calculate the cracking angle according to the dynamic yield equation.

[0118] In this example, before cracking, the material is regarded as homogeneous, and its behavior can be described as elastic behavior. Once the traction force on the potential fracture surface satisfies the yield surface condition related to the strain rate, the material will be divided into a fracture process zone (FPZ) and an external main part. For any potential fracture surface, the normal traction force and shear traction force at the mesoscopic scale can be obtained by calculating from the macroscopic stress through the continuum mechanics equation:

[0119]

[0120] In the formula, n represents the normal vector at an angle θ with the horizontal axis, and n1 and n2 respectively represent the corresponding unit vector components in the global coordinate system; the dynamic yield function is:

[0121]

[0122] For each direction angle θ, different yield function values y will be generated for the traction force and DIF corresponding to the potential fracture plane d , with y d The direction angle with the maximum value should be regarded as the critical plane, and the direction angle of the critical plane corresponds to the maximum value of y d . When the derivative of y d with respect to θ is equal to zero, this direction angle can be determined analytically:

[0123]

[0124] Taking the derivative of both sides of the above formula with respect to θ, the critical crack angle θ can be solved and analyzed cr . The analytical critical crack angle θ cr is the angle corresponding to the maximum yield function value . If the material is still in the elastic state, the solution at this time represents the angle of the most critical potential fracture plane. If then the crack yields, and the cracking direction is judged according to the critical crack angle. DIF T and DIF Sis determined based on the normal strain rate and shear strain rate on the plane:

[0125]

[0126] wherein, and are the normal strain rate and shear strain rate on the potential fracture plane respectively; is the strain rate at the macroscopic scale.

[0127] S5: Experiment and Simulation:

[0128] S5.1 Experimental Setup: In this embodiment, three laboratory tests were simulated to demonstrate the performance of the proposed model in predicting rate-dependent problems, including dynamic direct shear tests, dynamic spalling tests, and dynamic mixed-mode three-point bending tests under constant normal load (CNL). The ability of the proposed model to predict dynamic fracture modes and mechanical behaviors was examined by comparing the numerical results with the experimental results. The model was integrated into the commercial finite element analysis method (FEM) software ABAQUS through a user-defined explicit subroutine VUMAT. The mechanical properties and constitutive parameters of the material are listed in Table 2.

[0129] The above simulations used three-node triangular elements and four-node quadrilateral elements under two-dimensional plane stress conditions. The characteristic length H k is calculated according to the area of each element through H k = Ω k / A k In two-dimensional problems, H k is calculated through H k = A k / l k where A k is the element area and l k is the length of the embedded crack in the element. To improve the computational efficiency during implementation, the calculation method is adopted, and this simple approximation method helps to achieve the independence of the meshless method and the mesh-based method in discretization.

[0130] Table 2 Mechanical Properties and Constitutive Parameters of the Material

[0131]

[0132] S5.2 Dynamic Direct Shear Experiment:

[0133] Simulate the shear rate-dependent behavior of coarse-grained granite under dynamic direct shear tests with a conventional normal load σ n = 2 Mpa. The test setup is as Figure 6(a). A vertical normal load is applied to the sample, and the bottom surface can slide horizontally freely (free slip), while the upper half of the sample is fixed horizontally. The experimental shear rate ranges from 1 μm / s to 10 mm / s, corresponding to a strain rate range from 1.67×10 -7 s -1 to 0.167 s -1 . This range represents a low dynamic shear strain rate spectrum, which is usually classified as "sub-seismic shear rate" in geomechanics. In this example, the mechanical properties and constitutive parameters are calibrated with data at a shear rate of 1 μm / s (strain rate of 1.67×10 -7 s -1 ). The strain rate - DIF S empirical equation is fitted according to the results of dynamic direct shear tests at five different shear rate levels, as shown in Figure 6 (b). Considering the symmetric macroscopic loading conditions, the test is simulated as a representative volume element (RVE) at a known horizontal fracture plane and characteristic length H = Ω / A k = 100 mm, where the contact surface is not included. Five specific experimental shear strain rate levels are replicated using different shear rates by applying displacement-controlled motion on the horizontal shear fracture plane.

[0134] As shown in Figure 7 (a), the predicted shear stress - displacement results show a high degree of agreement with the experimental curve. It should be noted that, in order to improve the comparability between the simulation results and the experimental data, the initial loading stage of the experimental curve has been replaced by a linear curve. This adjustment is necessary because the original curve exhibits inelastic behavior due to the compressive displacement between the loading system and the specimen during the experiment. As shown in Figure 7 (b), the shear stress - displacement simulation results at five different shear rate levels are shown, indicating three distinct stages: First, the shear stress increases within the elastic range until it reaches a peak. At this stage, the traction state of the shear band (i.e., shear stress t s and normal stress t n ) is still below the dynamic yield surface. The second stage is the peak shear stress stage, where shear cracks begin to form and penetrate the rock. At this stage, the traction state reaches the dynamic yield surface, and damage and dilation start to accumulate. The last stage is the post-peak stage, where the cracks exhibit softening behavior until the residual strength is reached. At this stage, as the dynamic yield surface shrinks, the damage and dilation values continue to increase.

[0135] As shown in Figure 8 (a), the trend of increasing shear strength with increasing shear strain rate in the simulation is consistent with the experimental results. As shown in Figure 8(b), the simulation results of the breakdown stress drop (defined as the difference between the peak shear strength and the residual shear strength) at five different shear rate levels are also completely within the range of the experimental results.

[0136] S5.3 Dynamic spalling experiment:

[0137] In this embodiment, the ability of the proposed model to capture the rate-dependent behavior of concrete under dynamic tensile loading is demonstrated. A series of dynamic spalling tests using a Split Hopkinson Pressure Bar (SHPB) system are conducted to achieve this. As Figure 9 (a) shows, the experimental geometry and loading configuration are presented. The setup includes an air launcher with a cylindrical projectile, an incident Hopkinson bar, and a concrete specimen (120 mm in length × 40 mm in diameter) in direct contact with the Hopkinson bar. The specimen is constrained in the vertical direction and allowed to move freely in the horizontal direction. The incident compressive stress wave generated by the projectile is transmitted through the incident bar to the specimen and reflected as a tensile wave at the free end, causing spalling of the specimen. The experiments cover a range of tensile strain rates over an order of magnitude, from 10 s-1 to 120 s-1, achieved by applying different impact velocities to the projectile.

[0138] To evaluate the effect of element size on the dynamic tensile strength of concrete at a strain rate of approximately 100 s-1, a convergence study is carried out. As Figure 9 (c) shows, the results of the study indicate that the numerical results tend to converge within the element size range of 0.1 mm to 1 mm, suggesting that the simulation results are robust within this range. Therefore, to optimize computational resources, an element size of 1 mm is selected for the specimen, while a coarser mesh is used for the incident bar and the projectile. As Figure 13 (b) shows, under simplified one-dimensional stress conditions and in the absence of mixed-mode fracture components, a single-crack model is adopted to improve computational efficiency while maintaining the accuracy of the solution.

[0139] In the simulation, the impact velocity of the projectile ranges from 6.0 m / s to 14.0 m / s. A typical stress-time curve recorded in the middle of the incident bar during the simulation is shown in Figure 9 (d). The static material properties and constitutive parameters are listed in Table 2. The Young's modulus of the aluminum alloy bar and the projectile is 69.5 GPa, and the Poisson's ratio is 0.32.

[0140] The dynamic tensile strength in the simulation is calculated based on the free surface velocity using the formula where ρ is the material density, is the stress wave propagation velocity, △v is the pull-back velocity at the free end of the specimen, L is the length of the specimen, and △t is the time required for the stress wave to pass through the specimen. The pull-back velocity is determined from the measured velocity at the free end (as shown in Figure 13(as shown in (d)). The tensile stress wave is the reflected wave of the transmitted compression wave at the free end of the specimen, which will cause the specimen to fracture when the critical tensile strength is reached. The strain rate is calculated based on the rise time of the tensile stress wave, and the formula is where σ F is the critical stress, E is the Young's modulus, and t c is the critical time from the start of specimen loading to fracture. For the higher strain rate range, the empirical formula of DIF T - strain rate is obtained by fitting the experimental results, while for the smaller strain rate range, such as Figure 9 (as shown in (e)).

[0141] such as Figure 10 shown, the rate dependence of dynamic tensile strength on strain rate among experimental data, current simulation results, and simulation results of existing studies is compared. The results show that the model proposed in this application effectively captures the trend of dynamic tensile strength increasing with the increase of strain rate in the experiment. In addition, some studies have shown that existing data processing methods may overestimate the dynamic tensile strength, which may lead to the overestimation of empirical DIF T , thus resulting in the difference between our simulation results and the reported results. As Figure 11 shown, it shows the fracture process zone (FPZ) of two typical fracture types observed in the simulation. Consistent with the experimental observations, only one crack is observed when the impact velocity is lower than 8.0 m / s, while two cracks are observed at higher impact velocities. The distance from the first crack to the free end ranges from 30 mm to 35 mm, and the distance of the second crack ranges from 77 mm to 86 mm.

[0142] S5.4 Mixed - mode three - point bending test:

[0143] In the foregoing embodiments, this application has systematically studied the dynamic responses of brittle materials under pure shear and pure tensile stress states. To verify the ability of the proposed model to capture the dynamic mixed - mode fracture process more commonly existing in practical engineering, a numerical simulation of the mixed - mode three - point bending test under dynamic loading is carried out.

[0144] The tested structural geometry and boundary conditions are as Figure 12 (as shown in (a)). The specimen length is 203.2 mm, the height is 76.2 mm, and the thickness is 25.4 mm. A vertical notch with a length of 19 mm is cut at the bottom of the specimen. The distance between the notch and the specimen center line is γL / 2, where γ is the coefficient determining the notch offset position. Specimens with five different notch positions are tested under impact loads, and the impact load is applied by imposing a half - sine acceleration - time relationship on the middle top element (as Figure 12(as shown in (b)). The bottom boundary allows free horizontal movement, while the vertical displacement is fixed. To evaluate the numerical convergence during mesh refinement, two different meshes were used along the crack path: Mesh 1 (element size 0.66 mm, 11,177 elements) and Mesh 2 (element size 1.52 mm, 4,679 elements), as Figure 12 (c) shown.

[0145] The material properties and constitutive parameters of the concrete are shown in Table 2. The elastic-plastic parameters were calibrated based on the force-crack opening displacement curve of the center-notched specimen under static loading conditions. Based on the peak force of the specimen with a notch offset ratio γ = 0.5 under dynamic loading, the DIF T / S -strain rate relationship (as Figure 12 (d) shown) was calibrated.

[0146] In the experiment, the existing literature estimated the strain rate based on the initial slope of the strain-time curve at the notch tip and assumed it to be constant along the beam span due to instrument limitations. In contrast, the numerical simulation captures the strain rate at the notch tip at each time step, revealing a trend of gradual increase throughout the loading process. To facilitate direct comparison, the simulated strain rate corresponding to approximately 15% of the peak force of the notch specimen with γ = 0.5 was selected to match the experimental measurement. In this study, a dynamic notch tip strain rate of 0.5 s-1 was used for all notched beams to study the interaction between the peak force and γ under static and dynamic conditions. It should be noted that when the strain rate is lower than 0.5 s-1, the stress wave propagation effect in concrete can be ignored. Therefore, dynamic analyses considering stress wave propagation, dispersion, and reflection effects were not included in this study.

[0147] As Figure 13 (a) shown, it compares the experimental and simulated load-crack opening displacement (COD) response curves and the final fracture morphology of the specimens with a center notch under static loading. The simulated load-COD curve shows a high degree of agreement with the experimental curve. The difference between the simulated and experimental values of the fracture energy (calculated by normalizing the area under the load-COD curve by the crack surface area) is approximately 3.86%, which is within an acceptable range. This difference may be due to the discreteness of the material parameters used in the simulation and the fact that the experimental curve essentially represents the envelope of the cyclic loading-unloading test at a static rate - the unloading process may have an impact on the mechanical response. In addition, Figure 13 (a) shows that the numerical results tend to converge with mesh refinement, which is attributed to the introduction of the characteristic length H. As Figure 13(b), which compares the interaction relationship between the γ parameter and the peak load under static and dynamic loads. In the legend, "dynamic, notch" indicates that the beam fractures at the notch under impact load, and "static, middle" represents that the fracture occurs in the mid-span region under static load, and so on. The analysis shows that the proposed model can effectively capture the evolution trend of the gradually increasing peak load under both static and dynamic loads when the notch position deviates from the mid-span.

[0148] As Figure 14 shown, Figure 14 (a1)-(a5) show the comparison of the simulated and experimental fracture modes of three-point bending tests with different notch offsets under static loading conditions. Figure 14 (b1)-(b5) show the comparison of the simulated and experimental fracture modes of three-point bending tests with different notch offsets under dynamic impact loading conditions. Generally, as the offset ratio increases, the fracture position moves from the notch tip to the mid-span. However, the transition threshold between static loading and dynamic loading is different. The simulation of this application accurately captures the influence of the loading rate on the failure mode transition point: under static loading, the transition from mixed-mode failure to tensile-mode failure occurs at the offset ratio γ = 0.7 (as shown in Figure 14 (a2)), while under dynamic loading, it occurs at γ = 0.766 (as shown in Figure 14 (b4)). Under dynamic loading with γ = 0.766, the fracture position of the beam in the experiment is either at the mid-span or at the notch tip, and both show the same peak load. In the simulation, fractures occur at both positions, but the mid-span fracture dominates the overall failure mode. This difference may be attributed to the inhomogeneity of the concrete specimens in the experiment, which results in different fractures for specimens with the same γ offset.

[0149] Example 2:

[0150] An analysis system for dynamic fracture of quasi-brittle materials, used to implement the analysis method for dynamic fracture of quasi-brittle materials described in Example 1, includes:

[0151] Unit construction module: used to construct a representative volume element (RVE) of quasi-brittle materials with two embedded fracture planes;

[0152] Model construction module: used to construct a two-scale model of the representative volume element (RVE), the two-scale model includes a macro-scale and a meso-scale sub-model, and the crack traction force and displacement jump increment in the meso-scale sub-model satisfy the rate-dependent cohesive friction law; the specific steps are:

[0153] Parameter identification module: used to identify the two-scale model parameters, the two-scale model parameters include material property parameters, elastoplastic parameters, and rate-dependent parameters;

[0154] Performance analysis module: Embed the dual-scale model into the finite element software to simulate the dynamic fracture of quasi-brittle materials.

[0155] In summary, the mesoscale crack plane accommodated by this application has the key characteristics of the dynamic fracture behavior of quasi-brittle materials. The cracking process in these local regions is described by DIF (Dynamic Fracture Factor) and damage plasticity mechanics, providing physically meaningful insights and having significant advantages over classical empirical models.

[0156] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to the above embodiments, those of ordinary skill in the art should understand that: the specific implementation manners of the present invention can still be modified or equivalently replaced, and any modification or equivalent replacement without departing from the spirit and scope of the present invention shall be covered by the protection scope of the claims of the present invention.

Claims

1. A method for analyzing dynamic fracture of quasi-brittle materials, characterized in that: The specific steps are: S1: Construct a representative volume element RVE of a quasi-brittle material with two embedded fracture planes; S2: constructing a dual-scale model of a representative volume unit RVE, wherein the dual-scale model includes a macroscopic scale and a mesoscopic scale sub-model, and the crack traction force and displacement jump increment in the mesoscopic scale sub-model satisfy the rate-dependent cohesive friction law; S3: Identifying the dual-scale model parameters, wherein the dual-scale model parameters include material property parameters, elastic-plastic parameters, and rate-dependent parameters; S4: The dual-scale model is embedded in the finite element software to simulate the dynamic fracture of quasi-brittle materials.

2. The method for analyzing dynamic fracture of quasi-brittle materials according to claim 1, characterized in that: The specific steps of constructing the dual-scale model in step S2 are: S2.1: Construct a mathematical relationship between the crack strain increment in the macroscopic sub-model and the crack displacement jump increment in the mesoscopic sub-model, wherein the crack strain increment is equal to the product of the crack displacement jump increment and the crack normal vector; S2.2: constructing a mathematical relationship between the crack displacement jump increment and the crack traction force in the mesoscopic sub-model, wherein the crack traction force and the displacement jump increment in the mesoscopic sub-model satisfy the rate-dependent cohesive friction law; S2.3: Construct a mathematical model of crack stress in the macroscopic sub-model. The virtual work generated by the crack stress and crack strain increments in the macroscopic sub-model should be equal to the total work completed by the stress and strain increments inside and outside the fracture surface in the mesoscopic sub-model.

3. The method for analyzing dynamic fracture of quasi-brittle materials according to claim 1 or 2, characterized in that: The rate-dependent cohesive friction law is a dynamic yield function of a combination of yield and failure in the form of a hyperbolic curve: Where, t n and t s are normal and tangential traction respectively, D is the damage variable, DIF T and DIF S are the tensile and shear dynamic enhancement factors, respectively, A=(1-D)f t + B, B = c (1-D) / [2 (1-D) + 2Du 2 ], μ is the internal friction coefficient, f t and c are the static tensile strength and static cohesion of the material, respectively.

4. The method for analyzing dynamic fracture of quasi-brittle materials according to claim 3, characterized in that: The damage variable D is: In the formula, is the plastic displacement of the crack in the mesoscopic sub-model, α and β are the coefficients that control the influence of normal displacement and shear displacement on damage; u n p and u s p are the plastic displacements in the normal and shear directions, respectively, and δ0 is used to normalize u p Parameters of relative displacement.

5. The method for analyzing dynamic fracture of quasi-brittle materials according to claim 4, characterized in that: The rate-dependent cohesive friction law describes the tensile and shear dynamic enhancement factors DIF T and DIF S They are: In the formula, is the tensile strain rate at the mesoscopic scale, is the shear strain rate at the mesoscopic scale.

6. The method for analyzing dynamic fracture of quasi-brittle materials according to claim 4, characterized in that: The material property parameters in step S3 include static tensile strength f t and static cohesion c; The elastic-plastic parameters include elastic tensile stiffness K n and elastic shear stiffness K s , the rate-dependent parameters include parameters α and β, which are respectively related to the fracture energy, and a relative displacement parameter δ0.

7. A system for analyzing dynamic fracture of quasi-brittle materials, characterized in that: The method for analyzing the dynamic fracture of a quasi-brittle material according to any one of claims 1 to 6 comprises: Unit building module: used to construct a representative volume unit RVE with two embedded fracture planes for quasi-brittle materials; Model construction module: used to construct a dual-scale model of a representative volume unit RVE, wherein the dual-scale model includes a macroscopic scale and a mesoscopic scale sub-model, and the crack traction force and displacement jump increment in the mesoscopic scale sub-model satisfy the rate-dependent cohesive friction law; Parameter identification module: used to identify the dual-scale model parameters, the dual-scale model parameters include material property parameters, elastic-plastic parameters and rate-dependent parameters; Performance Analysis Module: Used to embed the dual-scale model into the finite element software to simulate the dynamic fracture of quasi-brittle materials.

Images