Quasi-brittle material fracture simulation method considering shear softening effect

Abstract

The invention provides a quasi-brittle material fracture simulation method considering a shear softening effect. The method comprises the following steps: step 10, constructing an extended bond-based near-field kinetic model of a quasi-brittle material; the extended bond-based near-field kinetic model combines a bond stretching softening effect and a shear softening effect; and step 20, solving the extended bond-based near-field kinetic model by adopting an arc length method. According to the quasi-brittle material fracture simulation method considering the shear softening effect, the shear softening effect is considered, the springback instability behavior of cracks of the quasi-brittle material can be accurately captured, the physical accuracy and numerical stability of crack simulation of the quasi-brittle material are improved, and the method is suitable for popularization and application. The invention discloses the influence of the type II critical energy release rate on a force-displacement curve and a material springback instability behavior.

Description

Technical Field

[0001] This invention belongs to the field of computational mechanics technology, specifically relating to a method for simulating the fracture of quasi-brittle materials that considers shear softening effects. Background Technology

[0002] Quasi-brittle materials (such as rock and concrete) are widely present in practical engineering. Their fracture process often exhibits a typical tension-shear hybrid fracture mode and may experience snap-back instability during failure, posing a significant challenge to structural safety analysis. Traditional finite element methods (FEM) are based on the theory of local continuum media, which makes it difficult to handle the discontinuities in geometric and physical fields during crack initiation and propagation, thus limiting their ability to simulate complex fracture behavior.

[0003] Peridynamics (PD) is an emerging theory of nonlocal continuum mechanics. Through its governing equations based on spatial integrals, it successfully overcomes the singularity problem caused by the absence of spatial derivatives in traditional continuum mechanics methods when dealing with discontinuous problems, thus possessing superior capabilities in simulating crack propagation. PD can be specifically divided into two main categories: bond-based PD and state-based PD. Bond-based PD, with its simple and clear physical concepts and relatively easy numerical implementation, has been widely used in simulating crack propagation. However, traditional bond-based models only consider the tensile and compressive forces between point pairs along the bond direction, failing to reflect the shear deformation mechanism. Furthermore, traditional models generally employ displacement-controlled loading methods, making it difficult to stably capture the springback instability behavior during the unstable crack propagation stage, and lacking consideration of type II energy release rates. G II System modeling and analysis of the coupling relationship between fracture process and fracture. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a method for simulating the fracture of quasi-brittle materials that takes into account the shear softening effect. Considering the shear softening effect, the method can accurately capture the springback instability behavior of quasi-brittle material cracks, thereby improving the physical accuracy and numerical stability of quasi-brittle material crack simulation.

[0005] To address the aforementioned technical problems, this invention provides a method for simulating the fracture of quasi-brittle materials considering shear softening effects, comprising the following steps: Step 10: Construct an extended bond-based peri-field dynamics model for the quasi-brittle material; the extended bond-based peri-field dynamics model combines bond stretching softening effect and shear softening effect; Step 20: Solve the near-field dynamics model of the extended bond base using the arc length method.

[0006] As a further improvement of the present invention, in the extended bond-base near-field dynamics model, the bond force expression is given by equation (1): Equation (1) In the formula, Indicates bond force. Indicates the normal stiffness of the key. Indicates bond elongation. Indicates the normal bond damage parameters. Indicates the tangential stiffness of the bond. Indicates the key angle value. This represents the tangential bond damage parameter.

[0007] As a further improvement of the present invention, the expression for the normal bond damage parameter in equation (1) is equation (2): Equation (2) In the formula, ; This represents the critical tensile value of the normal bond. , Indicates the tensile strength of quasi-brittle materials. Indicates the normal stiffness of the key. Indicates the thickness of the two-dimensional model. Indicates the near-field radius; This represents the elongation corresponding to the inflection point of the softening curve in the trilinear softening model. ; Indicates the maximum elongation of the normal bond. ; This represents the area under the softening curve in the trilinear softening model. , Indicates the type I critical energy release rate; The expression for the tangential bond damage parameter in equation (1) is equation (3): Equation (3) In the formula, This represents the critical rotation angle value of the tangential bond. , Indicates the shear strength of quasi-brittle materials. Indicates the tensile strength of quasi-brittle materials. Indicates the shear modulus of quasi-brittle materials; This represents the maximum rotation angle of the tangential key. , This indicates the Type II critical energy release rate.

[0008] As a further improvement of the present invention, in step 20, the overall stiffness matrix, load factor and displacement vector are obtained by solving.

[0009] As a further improvement of the present invention, in step 20, the displacement vector and load factor of the node are used as variables, and the overall stiffness matrix, load factor and displacement vector of the current loading step are calculated by using an iterative solution method based on the overall stiffness matrix, load factor and displacement vector of the previous loading step. Specifically, the calculation of the (k+1)th iteration in the (n+1)th loading step includes: Step 210: Apply boundary conditions; Step 220: Calculate the internal force matrix and the unbalanced force matrix; Step 230: Calculate the displacement increment per unit load direction at the (k+1)th iteration. Displacement increment caused by unit load increment ; Step 240: Calculate the load factor correction amount for the (k+1)th iteration using a correction algorithm. ; Step 250: Based on the displacement increment in the unit load direction at the (k+1)th iteration. Displacement increment caused by unit load increment and load factor correction The displacement correction amount at the (k+1)th iteration is calculated. ; Step 260, based on the load factor of the nth loading step and the load factor correction at the (k+1)th iteration The load coefficient at the (k+1)th iteration of the (n+1)th loading step is obtained. ; Step 270, based on the displacement vector of the nth loading step and the displacement correction amount at the (k+1)th iteration The displacement vector at the (k+1)th iteration of the (n+1)th loading step is obtained. ; Step 280, calculate bond elongation. and bond angle value ; Step 290, based on the bond elongation Calculated normal bond damage parameters According to the key angle value Tangential key damage parameters were calculated. ; Step 300, based on the normal bond damage parameters The local tensile damage of the mass point is calculated; based on the tangential key damage parameters... The local shear damage of the mass point is calculated; Step 400: Update the overall stiffness matrix; Step 401: If the convergence criterion is met, stop the iteration; otherwise, proceed to the next iteration.

[0010] As a further improvement of the present invention, in step 230, the displacement increment in the unit load direction at the (k+1)th iteration is calculated using equation (4). : Equation (4) In the formula, This represents the tangential stiffness matrix at the k-th iteration of the (n+1)-th loading step. This represents the unbalanced force vector during the k-th iteration of the (n+1)-th loading step; The displacement increment caused by the unit load increment at the (k+1)th iteration is calculated using equation (5). : Equation (5) In the formula, Represents the external force vector; In step 250, the displacement correction amount at the (k+1)th iteration is calculated using equation (6). : Equation (6) In the formula, This represents the load factor correction amount at the (k+1)th iteration.

[0011] As a further improvement of the present invention, the correction algorithm in step 240 specifically includes: if ,but ;otherwise, , ;if and ,but ;if and ,but ;if and If ,but ,otherwise ;otherwise ;in, , , , and for Two different solutions, Indicates the specified arc length.

[0012] As a further improvement of the present invention, in step 260, the load factor at the (k+1)th iteration of the (n+1)th loading step is calculated using equation (7). : Equation (7) In the formula, This represents the load factor for the nth loading step. This represents the load factor increment during the k-th iteration of the (n+1)-th loading step. This represents the load factor correction amount at the (k+1)th iteration; In step 270, the displacement vector at the (k+1)th iteration of the (n+1)th loading step is calculated using equation (8). : Equation (8) In the formula, This represents the displacement vector at the nth loading step. This represents the displacement increment during the k-th iteration of the (n+1)-th loading step. This represents the displacement correction amount at the (k+1)th iteration.

[0013] As a further improvement of the present invention, in step 290, the local tensile damage of mass x is calculated using equation (9): Equation (9) In the formula, This represents the local tensile damage of particle x. This represents the near-field region of particle x. Represents the small volumetric element related to the bond within the near-field region; The local shear damage of particle x is calculated using equation (10): Equation (10) In the formula, This represents the local shear damage of particle x.

[0014] As a further improvement of the present invention, in step 400, the stiffness matrix of the near-field dynamic region is coupled with the stiffness matrix of the finite element region to obtain the overall stiffness matrix.

[0015] Compared with the prior art, the technical solution of the present invention has the following beneficial effects: This invention provides a method for simulating the fracture of quasi-brittle materials considering shear softening effects. First, an extended bond-based near-field dynamics model of the quasi-brittle material is constructed. During model construction, considering the shear softening effect, bond deformation is decomposed into normal and tangential deformations. Normal and tangential bond forces are introduced into the bond force expression, resulting in a bond force expression with directional discrimination capability. The normal force is related to the Type I critical energy release rate, and the tangential force is related to the Type II critical energy release rate, thus achieving physical modeling of the tension-shear coupled fracture behavior of quasi-brittle materials. Then, the arc-length method is used to solve the extended bond-based near-field dynamics model, and a nonlinear coupling path tracking strategy between incremental load and displacement is constructed to accurately capture the springback instability behavior of quasi-brittle material cracks. This improves the physical accuracy and numerical stability of quasi-brittle material crack simulation. Attached Figure Description

[0016] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the embodiments of the present invention will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0017] Figure 1 This is a flowchart of a method for simulating the fracture of quasi-brittle materials considering shear softening effects, provided in an embodiment of the present invention. Figure 2 This is a schematic diagram of the near-field dynamics trilinear stretching and softening curve in an embodiment of the present invention; Figure 3 This is a schematic diagram of the near-field dynamics bilinear shear softening curve in an embodiment of the present invention; Figure 4 This is a schematic diagram of a two-dimensional near-field dynamics-finite element coupling model in an embodiment of the present invention; Figure 5 This is a schematic diagram of the geometry and boundary conditions of the four-point bending test of the Gálvez notched beam; Figure 6 This is a schematic diagram of the discrete model of the Gálvez notched beam in the method of this embodiment of the invention; Figure 7 This is a schematic diagram of the failure process of the Gálvez notched beam at four points of bending, obtained by PD simulation. Figure 8 These are comparison diagrams of the crack paths of the Gálvez notched beam at four points of bending, obtained through experimental observation and PD simulation, respectively. Figure 9 This is a comparison chart of the CMOD and CMSD curves of the four-point bend of the Gálvez notched beam; Figure 10The figures are comparisons of the force-displacement curves of the four-point bending of the Gálvez notched beam obtained by simulation using the method of the embodiments of the present invention, displacement method simulation, and experimental observation, respectively. Figure 11 This is a schematic diagram illustrating the effect of the Type II critical energy release rate on the simulation results of the four-point bending test of the Gálvez notched beam in the method of this embodiment of the invention; Figure 12 (a) is a schematic diagram of the geometry and boundary conditions of the Arrea single-notch shear beam. Figure 12 (b) is a schematic diagram of the geometry and boundary conditions of the Arrea single-notch shear beam in the numerical scheme; Figure 13 This is a schematic diagram of the discrete model of the Arrea single-notch shear beam in the method of this embodiment of the invention; Figure 14 This is a schematic diagram of the failure process of the Arrea single-notch shear beam obtained by PD simulation; Figure 15 These are comparison diagrams of crack paths in the Arrea single-notch shear beam obtained through experimental observation and PD simulation, respectively. Figure 16 (a) is a comparison of the force-displacement curves at point A of the Arrea single-notch shear beam obtained by simulation using the method of the present invention and PD simulation, respectively; Figure 12 (b) is a comparison of the force-displacement curves at point B of the Arrea single-notch shear beam obtained by simulation using the method of the present invention and PD simulation, respectively; Figure 17 (a) is a crack path diagram of the Arrea single-notch shear beam simulated using PD. Figure 17 (b) is a diagram of the final crack morphology of the Arrea single-notch shear beam simulated using PD. Figure 18 This is a schematic diagram illustrating the influence of the Type II critical energy release rate on the simulation results of the Arrea single-notch shear beam test in the method of this embodiment of the invention. Detailed Implementation

[0018] The technical solution of the present invention will now be described in detail with reference to the accompanying drawings.

[0019] The failure of quasi-brittle materials typically exhibits a mixed tensile-shear fracture mechanism. However, most existing near-field dynamic models only consider the tensile softening behavior of materials, neglecting shear softening. To more accurately simulate the failure process of quasi-brittle materials, this invention provides a method for simulating the fracture of quasi-brittle materials that considers the shear softening effect, such as... Figure 1 As shown, it includes the following steps: Step 10: Construct an extended bond-based peri-field dynamics model for quasi-brittle materials. The extended bond-based peri-field dynamics model combines bond stretching softening effect and shear softening effect.

[0020] Step 20: Solve the near-field dynamics model of the extended bond base using the arc length method.

[0021] In step 10, preferably, in the extended bond-base near-field dynamics model, the bond force expression is Equation (1): Equation (1) In the formula, Indicates bond force. Indicates the normal bond force. Indicates tangential bond force; , Indicates the normal stiffness of the key. Indicates bond elongation. Indicates the normal bond damage parameters; , Indicates the tangential stiffness of the bond. Indicates the key angle value. This represents the tangential bond damage parameter.

[0022] The normal bond force is modeled using a trilinear softening model consistent with the macroscopic response of quasi-brittle materials, corresponding to the three stages of quasi-brittle material failure. For example... Figure 2 As shown, when the bond elongation Less than the critical tensile value of the normal bond At this time, the material is in the elastic response stage. When the bond elongation... greater than the critical tensile value of the normal bond At this point, the bond begins to exhibit a softening effect. The elongation corresponding to the inflection point of the softening curve in the trilinear softening model is... Normal bond force The final elongation corresponding to zero, that is, the maximum elongation of the normal bond, is .

[0023] Therefore, the normal bond damage parameter in equation (1) Expressed using equation (2): Equation (2) In the formula, This represents the parameters of the trilinear softening model. ; This represents the critical tensile value of the normal bond. , Indicates the tensile strength of quasi-brittle materials. Indicates the normal stiffness of the key. Indicates the thickness of the two-dimensional model. Indicates the near-field radius; This represents the elongation corresponding to the inflection point of the softening curve in the trilinear softening model. ; Indicates the maximum elongation of the normal bond. ; This represents the area under the softening curve in the trilinear softening model. , This represents the Type I critical energy release rate.

[0024] Preferably, the modeling of tangential bond forces adopts the following method: Figure 3 The bilinear shear softening model is shown. When the bond angle value... The critical rotation angle greater than or equal to the tangential bond At this point, the tangential bond force begins to soften. Once the bond rotation angle value... Exceeding the maximum rotation angle of the tangential key The tangential bond force rapidly decreases to zero.

[0025] Therefore, the tangential bond damage parameter Expressed using equation (3): Equation (3) In the formula, This represents the critical rotation angle value of the tangential bond. , Indicates the shear strength of quasi-brittle materials. Indicates the tensile strength of quasi-brittle materials. Indicates the shear modulus of quasi-brittle materials; This represents the maximum rotation angle of the tangential key. , This indicates the Type II critical energy release rate.

[0026] The arc length method used in step 20 combines force and displacement increments, and tracks the nonlinear path in the form of arc length to obtain the overall stiffness matrix, load factor and displacement vector, so as to solve the instability behavior caused by material damage or fracture.

[0027] In the arc-length method, the displacement vector and load factor It is considered an unknown quantity. (Point mass) The static equilibrium equation can be reformulated as equation (13): Equation (13) In the formula, Indicates internal force. Indicates the action on the particle. External force vector on, This represents the load factor for the nth loading step.

[0028] Assuming the (n+1)th loading step satisfies the static equilibrium condition, the load factor is: Therefore, the new static equilibrium equation can be expressed as equation (14): Equation (14) In the formula, Represents a point mass The unbalanced force vector is used to correct the displacement vector so that it satisfies the equilibrium condition of the (n+1)th loading step.

[0029] Equation (14) can be expressed in matrix form as Equation (15): Equation (15) In the formula, , , , where nf represents the total degrees of freedom.

[0030] The above nonlinear equation set (15) needs to be solved using an iterative method. Assume that the displacement vector of the (k+1)th iteration step of the (n+1)th loading step is given by equation (16): Equation (16) In the formula, , This represents the displacement correction amount at the (k+1)th iteration. This represents the displacement vector at the nth loading step.

[0031] In order to calculate The system equations are approximated by the first-order Taylor expansion as equation (17): Equation (17) In the formula, Let represent the unbalanced force vector at the k-th iteration of the (n+1)-th loading step. The first term is the tangential stiffness matrix of the system at the k-th iteration, which can be expressed as equation (18): Equation (18) During the iteration process, a convergence criterion needs to be specified to determine whether the equation converges. For example, the force convergence criterion can be expressed as equation (19): Equation (19) In the formula, This indicates the preset convergence tolerance.

[0032] In the arc-length method, the load factor of the (n+1)th loading step It is not a fixed value; it must be calculated in each iteration step. Equation (20) Equation (21) In the formula, This represents the load factor correction amount at the (k+1)th iteration.

[0033] Since the arc-length method treats not only nodal displacements but also load factors as variables, an additional constraint equation is required when solving the nonlinear equation system. The constraint equation can be expressed as equation (22): Equation (22) In the formula, This indicates a specified arc length, which determines the incremental changes in the nodal displacement vector or force vector.

[0034] In order to solve equations (17) and (22). The form can be assumed to be equation (6): Equation (6) The expression is given by equation (4): Equation (4) The expression is given by equation (5): Equation (5) In equation (22) Depend on Replacement yields information about The quadratic equation is given by equation (23): Equation (23) In the formula, , , .

[0035] Equation (23) generally has two different solutions. To obtain the correct solution, the displacement increment at the (k+1)th iteration can be required. Displacement increment at the kth iteration They form an acute angle. The angle between them is: Equation (24) This indicates the existence of a correct solution. If... and This indicates the existence of two acute angles, in which case a linear solution closer to equation (23) needs to be selected. If equation (23) has two imaginary roots, the specified arc length can be reduced and the calculation recalculated.

[0036] This invention employs a coupled near-field dynamics and finite element method for implicit solution. For example... Figure 4As shown, the square plate is discretized at its center with a uniform near-field dynamics (PD) mesh (orange nodes) with a spacing of Δx, while the rest of the plate is discretized with a finite element mesh. The finite element portion near the PD region uses a uniform mesh with the same spacing as the PD nodes, and its thickness (green rhombus) must be greater than or equal to the radius of the near-field region. The external finite element part can be discretized using a coarser non-uniform mesh, depending on the specific boundary shape.

[0037] The global equilibrium equation for the PD region is: Equation (11) In the formula, The stiffness matrix represents the near-field dynamic region. This represents the displacement matrix of the near-field dynamic nodes. This represents the internal force matrix of the near-field dynamic nodes.

[0038] The equilibrium equations of the finite element region are in the form of: Equation (12) In the formula, Represents the stiffness matrix of the finite element region. This represents the displacement matrix of the finite element nodes. This represents the internal force matrix of a finite element node.

[0039] Finally, the stiffness matrices of the FEM and PD parts are coupled to obtain the overall stiffness matrix.

[0040] Finally, the stiffness matrices of the FEM and PD parts are coupled to obtain the overall stiffness matrix.

[0041] In step 20, the displacement vector and load factor of the node are used as variables. Based on the overall stiffness matrix, load factor and displacement vector of the previous loading step, the overall stiffness matrix, load factor and displacement vector of the current loading step are calculated by an iterative solution method.

[0042] Specifically, the calculation of the (k+1)th iteration in the (n+1)th loading step includes: Step 210: Apply boundary conditions. Specifically, apply fixed displacement constraints to some of the mass points.

[0043] Step 220: Calculate the internal force matrix based on the overall stiffness matrix and displacement matrix. Calculate the unbalanced force matrix based on the internal force matrix and external force matrix.

[0044] Step 230: Calculate the displacement increment per unit load direction at the (k+1)th iteration using equation (4). : Equation (4) In the formula, This represents the tangential stiffness matrix at the k-th iteration of the (n+1)-th loading step. This represents the unbalanced force vector during the k-th iteration of the (n+1)-th loading step.

[0045] The displacement increment caused by the unit load increment at the (k+1)th iteration is calculated using equation (5). : Equation (5) In the formula, This represents the external force vector.

[0046] Step 240: Calculate the load factor correction amount for the (k+1)th iteration using a correction algorithm. .

[0047] Specifically, if ,but .otherwise, There are two solutions, respectively , . , , , ;if and ,but ;if and ,but ;if and If ,but ,otherwise ;otherwise .

[0048] Step 250: Based on the displacement increment in the unit load direction at the (k+1)th iteration. Displacement increment caused by unit load increment and load factor correction The displacement correction amount at the (k+1)th iteration is calculated using equation (6). : Equation (6) In the formula, This represents the displacement correction amount at the (k+1)th iteration. This represents the displacement increment per unit load direction at the (k+1)th iteration.

[0049] Step 260, based on the load factor of the nth loading step and the load factor correction at the (k+1)th iteration The load factor at the (k+1)th iteration of the (n+1)th loading step is calculated using equation (7): Equation (7) In the formula, This represents the load factor in the (k+1)th iteration of the (n+1)th loading step. This represents the load factor for the nth loading step. This represents the load factor increment during the k-th iteration of the (n+1)-th loading step. This represents the load factor correction amount at the (k+1)th iteration.

[0050] Step 270, based on the displacement vector of the nth loading step and the displacement correction amount at the (k+1)th iteration The displacement vector at the (k+1)th iteration of the (n+1)th loading step is calculated using equation (8): Equation (8) In the formula, This represents the displacement vector at the (k+1)th iteration of the (n+1)th loading step. This represents the displacement vector at the nth loading step. This represents the displacement increment during the k-th iteration of the (n+1)-th loading step. This represents the displacement correction amount at the (k+1)th iteration.

[0051] Step 280: Calculate the bond elongation using equation (25). : Equation (25) In the formula, Represents a relative position vector. This represents the relative displacement vector.

[0052] The key rotation angle value is calculated using equation (26). : Equation (26) In the formula, Indicates the current direction of the key. .

[0053] Step 290, based on the bond elongation The normal bond damage parameters are calculated using equation (2). Based on the angle value The tangential bond damage parameters are calculated using equation (3). .

[0054] Step 300, based on the normal bond damage parameters The local tensile damage of particle x is calculated using equation (9): Equation (9) In the formula, This represents the local tensile damage of particle x. This represents the near-field region of particle x. This represents a tiny volume element related to bonds within the near-field region.

[0055] Based on tangential key damage parameters The local shear damage of particle x is calculated using equation (10): Equation (10) In the formula, This represents the local shear damage of particle x.

[0056] Step 400: Couple the stiffness matrix of the near-field dynamic region with the stiffness matrix of the finite element region to obtain the overall stiffness matrix.

[0057] Step 401: If the convergence criterion is met, stop the iteration; otherwise, proceed to the next iteration.

[0058] After the iterative calculation is completed, the overall stiffness matrix, load factor and displacement vector of the (n+1)th loading step are output.

[0059] Two specific examples are provided below.

[0060] Example 1: Simulation of a four-point bending test on a Gálvez notched beam The geometry and boundary conditions of the specimen are as follows: Figure 5 As shown, the specimen is 675 mm long, 150 mm high, and 50 mm thick, with a 75 mm high vertical notch at the center of the bottom. The lower left end of the specimen... y The direction is constrained, and the lower right end is subject to... x and y The upper left end is subject to directional constraints. y Directional constraints.

[0061] The material parameters of the Gálvez notched beam are as follows: Young's modulus E = 38 GPa, Poisson's ratio v = 0.2, tensile strength... =3MPa, Type I critical energy release rate =2700 N / m. Type II critical energy release rate of materials such as rock and concrete. Much greater than the Type I critical energy release rate In this example, we select Considering the balance between computational cost and simulation accuracy, the discretization parameter chosen in this example is: near-field radius. The particle spacing is 3 mm. It is 1mm. Figure 6 The coupled discrete case of PD-FEM is shown. The current model contains a total of 15,795 nodes, including 10,287 PD nodes and 5,508 FEM nodes, with a total degree of freedom of 31,590.

[0062] Figure 7 The crack propagation process in a four-point bend test is illustrated. It can be seen that the crack first initiates at the notch, then propagates at an inclined angle upwards and to the right until it reaches the boundary of the model. From... Figure 8 It can be seen that the crack path simulated by PD matches the experimental observation results well. Figure 9 The force-crack mouth opening displacement (CMOD) curve and the force-crack mouth sliding displacement (CMSD) curve were further shown. It can be seen that the CMSD at the peak strength of the four-point bend test is slightly smaller than the CMOD, but the final CMSD is larger than the CMOD. This indicates that in the four-point bend test, the specimen fractures in a tension-shear mixed mode, and the effect of shear is not negligible.

[0063] The force-displacement curve of the four-point bend test shows obvious springback. Figure 10 The simulation results of the displacement control method and the method of the present invention were compared. It can be seen that the prediction results of the displacement control method cannot capture the springback instability behavior of the material, but the method of the present invention can.

[0064] Figure 11 Further demonstrating the Type II critical energy release rate The impact on the simulation results. It can be seen that... The peak load of the curve has a significant impact, showing that the peak load increases with... The trend is increasing monotonically. It is noteworthy that the Type II critical energy release rate... It also has a significant impact on the springback phenomenon of the curve. It can be seen that, with... As the value increases, the rebound of the curve gradually weakens. When At that time, the force-displacement curve showed almost no rebound phenomenon.

[0065] Example 2: Simulation of the Arrea single-notch shear beam test The geometry and boundary conditions of the pattern are as follows: Figure 12As shown in (a). The experiment was conducted under antisymmetric four-point shear loading conditions. By applying a load to point C of the AB-type steel beam and using crack slip displacement as the control parameter, the fracture process exhibited significant springback behavior. In PD and other numerical methods, two steel blocks, 40 mm long and 20 mm wide, are typically placed at the lower left support and the right loading point to mitigate stress concentration effects. Figure 12 As shown in (b), the loading conditions are asymmetric, meaning that cracks initiating and propagating from the pre-existing defect end will exhibit a mixed pattern of opening and sliding. Therefore, the tensile softening and shear softening behavior of the material requires special attention.

[0066] The material parameters of the Arrea single-notch shear beam are: Young's modulus E = 24.8 GPa, Poisson's ratio v = 0.18, and tensile strength... =4 MPa, Type I critical energy release rate =150 N / m, Assume the steel block is a linearly elastic material with Young's modulus E = 200 GPa and Poisson's ratio v = 0.18. Discrete parameters for this example: near-field radius. The particle spacing is 10 mm. It is 2 mm. Figure 13 The coupled discrete array contains a total of 19,087 nodes, including 10,670 PD nodes and 8,395 FEM nodes, with a total degree of freedom of 38,174.

[0067] Figure 14 The crack propagation process under different loading steps is illustrated. The crack initially initiates at the end of the pre-existing defect and then propagates towards the right side of the steel block in a hybrid fracture mode. A comparison between the crack path predicted by PD and the experimental observation results is shown below. Figure 15 The two showed good consistency.

[0068] To describe the springback phenomenon in this example, the force-displacement curves at points A and B are as follows: Figure 16 As shown in the figure, it can be seen that the method of this embodiment can successfully capture the springback instability phenomenon of quasi-brittle materials. From Figure 16 It can be seen that in the early stages of loading, point A is... y The displacement is downwards; with further loading, point A... y The displacement gradually increases in the direction of displacement. The same phenomenon was also observed in the finite element numerical calculation results, see [link to details]. Figure 17 .

[0069] Figure 18 Demonstrated Type II critical energy release rate The effect on the force-displacement curve. It can be seen that under different Type II critical energy release rates, the force-displacement curves obtained by the method of this embodiment are basically the same before the peak load, indicating that the pre-peak stage is mainly affected by the Type I critical energy release rate. Control. It mainly has a significant impact on the springback phenomenon of the force-displacement curve, and the springback instability of the curve is similar to... They are negatively correlated, that is The larger the value, the less pronounced the rebound phenomenon, indicating that the post-peak state is mainly determined by the type II critical energy release rate. control.

[0070] The present invention constructs a novel extended bond-based peri-field dynamics model that can consider bond stretching and shear softening effects; introduces the arc length method into peri-field dynamics, constructs an implicit solution scheme based on the arc length method, successfully captures the springback instability phenomenon caused by strain softening in quasi-brittle materials; and reveals the influence of the type II critical energy release rate on the springback instability behavior.

[0071] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for simulating the fracture of quasi-brittle materials considering shear softening effects, characterized in that, Includes the following steps: Step 10: Construct an extended bond-based peri-field dynamics model for the quasi-brittle material; the extended bond-based peri-field dynamics model combines bond stretching softening effect and shear softening effect; Step 20: Solve the near-field dynamics model of the extended bond base using the arc length method.

2. The method for simulating quasi-brittle material fracture considering shear softening effect according to claim 1, characterized in that, In the extended bond base near-field dynamics model, the bond force expression is given by equation (1): Equation (1) In the formula, Indicates bond force. Indicates the normal stiffness of the key. Indicates bond elongation. Indicates the normal bond damage parameters. Indicates the tangential stiffness of the bond. Indicates the key angle value. This represents the tangential bond damage parameter.

3. The method for simulating quasi-brittle material fracture considering shear softening effect according to claim 2, characterized in that, The expression for the normal bond damage parameter in equation (1) is given by equation (2): Equation (2) In the formula, ; This represents the critical tensile value of the normal bond. , Indicates the tensile strength of quasi-brittle materials. Indicates the normal stiffness of the key. Indicates the thickness of the two-dimensional model. Indicates the near-field radius; This represents the elongation corresponding to the inflection point of the softening curve in the trilinear softening model. ; Indicates the maximum elongation of the normal bond. ; This represents the area under the softening curve in the trilinear softening model. , Indicates the type I critical energy release rate; The expression for the tangential bond damage parameter in equation (1) is equation (3): Equation (3) In the formula, This represents the critical rotation angle value of the tangential bond. , Indicates the shear strength of quasi-brittle materials. Indicates the tensile strength of quasi-brittle materials. Indicates the shear modulus of quasi-brittle materials; This represents the maximum rotation angle of the tangential key. , This indicates the Type II critical energy release rate.

4. The method for simulating quasi-brittle material fracture considering shear softening effect according to claim 1, characterized in that, In step 20, the overall stiffness matrix, load factor, and displacement vector are obtained by solving.

5. The method for simulating quasi-brittle material fracture considering shear softening effect according to claim 4, characterized in that, In step 20, the displacement vector and load factor of the node are used as variables. Based on the overall stiffness matrix, load factor and displacement vector of the previous loading step, the overall stiffness matrix, load factor and displacement vector of the current loading step are calculated by an iterative solution method. Specifically, the calculation of the (k+1)th iteration in the (n+1)th loading step includes: Step 210: Apply boundary conditions; Step 220: Calculate the internal force matrix and the unbalanced force matrix; Step 230: Calculate the displacement increment per unit load direction at the (k+1)th iteration. Displacement increment caused by unit load increment ; Step 240: Calculate the load factor correction amount for the (k+1)th iteration using a correction algorithm. ; Step 250: Based on the displacement increment in the unit load direction at the (k+1)th iteration. Displacement increment caused by unit load increment and load factor correction The displacement correction amount at the (k+1)th iteration is calculated. ; Step 260, based on the load factor of the nth loading step and the load factor correction at the (k+1)th iteration The load coefficient at the (k+1)th iteration of the (n+1)th loading step is obtained. ; Step 270, based on the displacement vector of the nth loading step and the displacement correction amount at the (k+1)th iteration The displacement vector at the (k+1)th iteration of the (n+1)th loading step is obtained. ; Step 280, calculate bond elongation. and bond angle value ; Step 290, based on the bond elongation Calculated normal bond damage parameters According to the key angle value Tangential key damage parameters were calculated. ; Step 300, based on the normal bond damage parameters The local tensile damage of the mass point is calculated; based on the tangential key damage parameters... The local shear damage of the mass point is calculated; Step 400: Update the overall stiffness matrix; Step 401: If the convergence criterion is met, stop the iteration; otherwise, proceed to the next iteration.

6. The method for simulating quasi-brittle material fracture considering shear softening effect according to claim 5, characterized in that, In step 230, the displacement increment in the unit load direction at the (k+1)th iteration is calculated using equation (4). : Equation (4) In the formula, This represents the tangential stiffness matrix at the k-th iteration of the (n+1)-th loading step. This represents the unbalanced force vector during the k-th iteration of the (n+1)-th loading step; The displacement increment caused by the unit load increment at the (k+1)th iteration is calculated using equation (5). : Equation (5) In the formula, Represents the external force vector; In step 250, the displacement correction amount at the (k+1)th iteration is calculated using equation (6). : Equation (6) In the formula, This represents the load factor correction amount at the (k+1)th iteration.

7. The method for simulating quasi-brittle material fracture considering shear softening effect according to claim 5, characterized in that, In step 240, the correction algorithm specifically includes: if ,but ;otherwise, , ;if and ,but ;if and ,but ;if and If ,but ,otherwise ;otherwise ;in, , , , and for Two different solutions, Indicates the specified arc length.

8. The method for simulating quasi-brittle material fracture considering shear softening effect according to claim 5, characterized in that, In step 260, the load factor at the (k+1)th iteration of the (n+1)th loading step is calculated using equation (7). : Equation (7) In the formula, This represents the load factor for the nth loading step. This represents the load factor increment during the k-th iteration of the (n+1)-th loading step. This represents the load factor correction amount at the (k+1)th iteration; In step 270, the displacement vector at the (k+1)th iteration of the (n+1)th loading step is calculated using equation (8). : Equation (8) In the formula, This represents the displacement vector at the nth loading step. This represents the displacement increment during the k-th iteration of the (n+1)-th loading step. This represents the displacement correction amount at the (k+1)th iteration.

9. The method for simulating quasi-brittle material fracture considering shear softening effect according to claim 5, characterized in that, In step 290, the local tensile damage of mass x is calculated using equation (9): Equation (9) In the formula, This represents the local tensile damage of particle x. This represents the near-field region of particle x. Represents the small volumetric element related to the bond within the near-field region; The local shear damage of particle x is calculated using equation (10): Equation (10) In the formula, This represents the local shear damage of particle x.

10. The method for simulating quasi-brittle material fracture considering shear softening effect according to claim 5, characterized in that, In step 400, the stiffness matrix of the near-field dynamic region is coupled with the stiffness matrix of the finite element region to obtain the overall stiffness matrix.

Images