A method for solving damage evolution parameters in a damage constitutive model

Abstract

This application discloses a method for solving damage evolution parameters in a damage constitutive model, belonging to structural mechanics and fracture mechanics. The method includes: establishing a mathematical relationship between damage variables and associated variables based on thermodynamic principles, and obtaining an integral expression for the energy dissipation potential per unit volume by combining damage evolution criteria; simplifying the integral expression for the energy dissipation potential per unit volume based on the principle of equivalent strain, obtaining a scalar analytical expression; numerically integrating the scalar analytical expression using the three-point Newton-Cotes integration algorithm, and iteratively solving for the damage evolution parameters using the secant method. This method transforms complex tensor operations into simple scalar analysis, reducing the computational burden. Furthermore, the use of the Newton-Cotes integration method and the secant method iterative solution strategy ensures accurate calculation of damage parameters, enhances computational stability and reliability, and significantly improves the efficiency and accuracy of the solution.

Description

Technical Field

[0001] This application relates to the fields of structural mechanics and fracture mechanics, and in particular to a method for solving damage evolution parameters in a damage constitutive model. Background Technology

[0002] In structural mechanics and fracture mechanics, continuum damage mechanics is widely used to describe the degradation of mechanical properties of materials under the process of microscopic damage accumulation. Damage constitutive models based on effective stress are widely used in the strength and failure analysis of engineering materials such as metals, composites, and porous media because they can accurately describe the damage evolution process. In these models, the damage variable is usually controlled by one or more evolution parameters. Although this method has good adaptability and generalization in theory, it faces a long-standing challenge in numerical implementation—the mesh dependency problem. The crack band model establishes a direct link between the fracture toughness of the material and the damage evolution law, achieving mesh-independent energy dissipation by controlling the energy consumed in damage propagation. However, to achieve this goal, the key lies in accurately calculating the damage evolution parameter A. m .

[0003] However, in existing technologies, the determination of this parameter relies on the integral expression of the damage energy density corresponding to the force. This expression typically has an integral structure in tensor form, making it difficult to solve explicitly. Especially in the finite element numerical implementation process, without a clear A... m The format of the expression will lead to model parameters needing to be obtained through repeated trial and error or parameter inversion, which is not only computationally inefficient but also highly dependent on user experience. Furthermore, inaccurate parameter assignment may cause damage propagation behavior to be highly dependent on mesh generation, thereby weakening the physical reliability of the model's prediction results.

[0004] Therefore, it is necessary to develop an A function that is simple in form, highly versatile, and easy to implement numerically. m Expression and solution methods have become key links in constructing high-precision damage models, eliminating mesh dependencies, and improving simulation credibility. Summary of the Invention

[0005] In view of the above-mentioned shortcomings in the prior art, this application provides a method for solving the damage evolution parameters in the damage constitutive model, which solves the problems of strong network dependence and complex solution in the prior art.

[0006] To achieve the aforementioned objectives, the technical solution adopted in this application is as follows:

[0007] This application provides a method for solving damage evolution parameters in a damage constitutive model, including:

[0008] S1: Based on the principles of thermodynamics, establish the mathematical relationship between the damage variable and the associated variable, and combine it with the damage evolution criterion to obtain the integral expression of the energy dissipation potential per unit volume.

[0009] S2: Based on the principle of equivalent strain, the integral expression of the energy dissipation potential per unit volume is simplified to obtain the scalar analytical expression;

[0010] S3: The scalar analytical expression is numerically integrated using the three-point Newton-Cotes integration algorithm, and the damage evolution parameters are iteratively solved using the secant method.

[0011] Further, S1 includes:

[0012] S101: Based on thermodynamic principles, the complementary free energy density of the damage constitutive model is defined;

[0013] S102: Based on the complementary free energy density, obtain the associated variable:

[0014]

[0015] in, For complementary free energy densities, As a damage variable, For damage evolution parameters, These are internal variables controlled by the laws governing damage evolution. To be related to damage variables Conjugate covariates, The Young's modulus of the matrix. It is the shear stiffness of the material. , , , , and The shear component of the true stress tensor. Representative material type;

[0016] S103: Based on the associated variables and complementary free energy densities, the energy dissipation potential per unit volume is obtained:

[0017]

[0018] in, The energy dissipation potential per unit volume, superscript Represents the derivative with respect to time;

[0019] S104: Based on the characteristic length and fracture toughness of the finite element method, the energy dissipation potential per unit volume is regularized to obtain the integral expression for the energy dissipation potential per unit volume:

[0020]

[0021] in, Energy dissipated per unit volume It is the energy release rate of the matrix. It is the length of the feature element. For time, Represents the element in finite element calculation. It represents a break.

[0022] Furthermore, based on the principle of equivalent strain, the integral expression of the energy dissipation potential per unit volume is simplified to obtain a scalar analytical expression, including:

[0023] S201: Considering the special case of applying a uniaxial tensile load to the material, a damage activation function is defined for the uniaxial tensile condition. The components of the effective stress tensor are obtained by finding the zero point of the damage activation function. The expressions for the damage activation function and the components of the effective stress tensor are as follows:

[0024]

[0025]

[0026] in, For damage activation function, and These represent the tensile and compressive ultimate strengths of the matrix material, respectively. These are the components of the effective stress tensor;

[0027] S202: Based on the principle of equivalent strain, the true stress tensor is defined, and the relationship between the components of the effective stress tensor and the true stress tensor is established:

[0028]

[0029]

[0030] in, For the effective stress tensor, Poisson's ratio of the matrix For transpose;

[0031] S203: Based on the relationship between the components of the effective stress tensor and the true stress tensor, To simplify the conditions, the associated variable per unit volume under uniaxial tension is obtained:

[0032]

[0033] in, The complementary free energy density under uniaxial tension;

[0034] S204: Based on the damage evolution law, the derivatives of the damage variable and the internal variables are obtained:

[0035]

[0036] S205: Based on the associated variables per unit volume, the effective strain tensor, and the derivatives of the damage variable and internal variables under uniaxial tension, the scalar analytical expression is obtained:

[0037] .

[0038] Furthermore, the step of numerically integrating the scalar analytical expression using the three-point Newton-Cotes integration algorithm and iteratively solving the damage evolution parameters using the secant method includes:

[0039] S401: Constructing the solution function based on scalar analytical expressions:

[0040] S402: Define the first two damage evolution parameters for the initial iteration;

[0041] S403: Based on the first two damage evolution parameters of the initial iteration, the three-point Newton-Cotes integral algorithm is used to solve the value of the improper integral in the solution function during the iteration process, and the function value of the solution function of the initial iteration is obtained;

[0042] S404: Based on the function value of the solution function in the initial iteration, the damage evolution parameters for the next iteration are obtained by iteratively solving the logarithmic secant method.

[0043] S405: Based on the damage evolution parameters of the next iteration, the three-point Newton-Cotes integral algorithm is used to solve the value of the improper integral in the solution function during the iteration process, and the function value of the solution function in the next iteration is obtained.

[0044] S406: Substitute the function value of the solution function for the next iteration into the solution function. If the absolute value of the solution function is less than the given integral tolerance, the iteration ends and the final damage evolution parameters are obtained; otherwise, update the damage evolution parameters of the iteration and proceed to S404 to continue the iterative solution.

[0045] Furthermore, the solution function is:

[0046]

[0047] in, For improper integrals, For the independent variable The relevant energy release rate, To and Relevant damage evolution parameters.

[0048] Furthermore, the first two damage evolution parameters of the initial iteration are:

[0049]

[0050] in, and These represent the material's elastic modulus and tensile strength, respectively. and These are the two damage evolution parameters for the initial iteration.

[0051] Furthermore, the three-point Newton-Cotes integration algorithm is as follows:

[0052]

[0053] in, The integration step size is... The lower limit of integration, This is the maximum number of points. Let be the integrand.

[0054] Furthermore, the increment of the integration step size in the three-point Newton-Cotes integration algorithm is a function of the number of integration steps:

[0055]

[0056] in, The integration step size is... These are parameter values ​​set by the user. For tensile or compressive strength, Represents stretching. Represents compression. This is the integration step.

[0057] Furthermore, the damage evolution parameters for the next iteration are obtained by iteratively solving using the logarithmic secant method:

[0058]

[0059] in, Representing the The damage evolution parameters obtained in the second iteration are also used as the first... i Initial values ​​of damage evolution parameters for +1 iteration. Representing the The first initial value of the damage evolution parameter in the next iteration. Representing the The second initial value of the damage evolution parameter in the next iteration. and Representing the first The integral values ​​corresponding to the two initial values ​​in the next iteration.

[0060] The beneficial effects of this application are:

[0061] This application provides a method for solving damage evolution parameters in a damage constitutive model. By transforming complex tensor operations into simple scalar analysis, the computational burden is significantly reduced. Furthermore, the use of the Newton-Cotes integration method and the secant method iterative solution strategy not only ensures the accurate calculation of damage parameters but also enhances the stability and reliability of the numerical algorithm, avoiding convergence problems common in traditional methods and significantly improving the efficiency and accuracy of solving damage evolution parameters. Attached Figure Description

[0062] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other embodiments can be obtained based on these drawings.

[0063] Figure 1 This is a flowchart illustrating a method for solving damage evolution parameters in a damage constitutive model provided in an embodiment of this application.

[0064] Figure 2 An iterative solution for A using the secant method provided in this application embodiment m A schematic diagram of the algorithm flow.

[0065] Figure 3 This is a flowchart illustrating a three-point Newton-Cotes integration algorithm provided in an embodiment of this application.

[0066] Figure 4 Damage evolution parameter A provided in the embodiments of this application m Diagram illustrating the impact of grid dependency.

[0067] Figure 5 The diagram shows the out-of-plane shear calculation results of the unidirectional fiber composite material RVE provided in the embodiments of this application. Detailed Implementation

[0068] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art based on this application are within the scope of protection of this application.

[0069] This application provides a method for solving damage evolution parameters in a damage constitutive model. This method can be found in [reference needed]. Figure 1 , Figure 1 The diagram shown is a flowchart illustrating a method for solving damage evolution parameters in a damage constitutive model according to an embodiment of this application, including:

[0070] S1: Based on the principles of thermodynamics, establish the mathematical relationship between the damage variable and the associated variable, and combine it with the damage evolution criterion to obtain the integral expression of the energy dissipation potential per unit volume.

[0071] In one embodiment of this application, to introduce the method of this application, the basic principles of continuous damage mechanics are briefly described, taking an isotropic elastic damage constitutive model as an example. Following a thermodynamic approach, the complementary free energy density in the material is first defined. This is a positive definite scalar function and must be zero at the origin of stress. The complementary free energy density of the isotropic damage model is defined as follows:

[0072]

[0073] in, For complementary free energy densities, and These are Young's modulus and Poisson's ratio of the matrix, respectively. , , , , and The shear component of the true stress tensor. This is the damage variable. We will only consider one damage variable here. It only affects the Young's modulus of the material. Therefore, it should be assumed that the Poisson's ratio of the material is not affected by the damage variable.

[0074] To ensure the irreversibility of the damage process, the rate of change of the complementary free energy must be greater than the rate of change of the residual energy of the applied stress.

[0075]

[0076] in, For stress tensor, For strain tensor, superscript It represents the derivative with respect to time.

[0077] After applying the chain derivation rule, the equations are expanded while ensuring the symmetry of the stress and strain tensors:

[0078]

[0079] The strain tensor is given by the derivative of the complementary free energy density with respect to the stress tensor:

[0080]

[0081] in, It is the first invariant of stress. It is a second-order unit tensor. This is the shear stiffness of the material. For ease of implementation, the strain tensor will be defined using engineering shear strain. The fourth-order compliance tensor and fourth-order stiffness tensor of the material can be obtained from the above equation, and for convenience, are given in Vogit form:

[0082]

[0083]

[0084] in, It is a flexibility matrix. It is the stiffness matrix. It is the shear modulus. and intermediate variables, and It is given by the following formula:

[0085]

[0086] The aforementioned stiffness tensor will allow for stress updates based on the increment of elastic strain and damage evolution. Prior to this, a damage initiation criterion needs to be defined, which will be introduced using the parabolic criterion given by Melro et al.

[0087] The damage activation function is defined by the following equation:

[0088]

[0089] in, For damage activation function, It is an internal variable controlled by the law of damage evolution. The load function is defined by the following formula:

[0090]

[0091] in, and These represent the tensile and compressive ultimate strengths of the matrix material, respectively, two invariants. and It is determined using the concept of effective stress. Effective stress tensor The calculation is as follows:

[0092]

[0093] in, It is a compliance matrix that does not consider damage.

[0094] Once damage begins, its evolution can be measured by the energy dissipation potential per unit volume:

[0095]

[0096] in, The energy dissipation potential per unit volume, To be related to damage variables The conjugate covariate, also known as the generalized thermodynamic force, This is an internal variable controlled by the damage evolution law. Based on the definition of complementary free energy density in the isotropic damage model, it is ensured that the thermodynamic force is always positive:

[0097]

[0098] The Kuhn-Tucker conditions for the damage model are determined using the Lagrange multiplier method:

[0099]

[0100] in, Represents a specific time The specific value of the load function.

[0101] When the material is completely damaged, the damage variable will take the value of 1, while the internal variable... The value will tend towards infinity. When the tangent stiffness tensor is negative definite, the damage will be confined to a narrow band with a thickness equal to that of the element activating the damage. Therefore, the structural response depends on the mesh size—the smaller the element is in the local damage band, the smaller the calculated dissipation potential. To circumvent this phenomenon, by combining the crack band model and utilizing the characteristic length and fracture toughness of the finite element, the calculated dissipation potential can be regularized:

[0102]

[0103] in, Energy dissipated per unit volume It is the energy release rate of the matrix. It is the length of the feature element. For time, Represents the element in finite element calculation. This represents fracture. For the completeness of the formula, a damage evolution criterion needs to be added; here, an exponential model is used for introduction:

[0104]

[0105] Among them, damage evolution parameters This requires solving the regularized dissipation potential equation. This formula is a function of the element characteristic length. Therefore, the parameters... Each finite element in the mesh is unique, and calculations only need to be performed in the initial stage of the finite element computation. Although the solution has been provided... The expression for the integral is given, but the first term in the integral is a function of the actual stress, and the second term is a function of the effective stress. If the integral expression is given in tensor form, it is very difficult to solve. Therefore, a concise method for calculation is required. .

[0106] S2: Based on the principle of equivalent strain, the integral expression of the energy dissipation potential per unit volume is simplified to obtain the scalar analytical expression.

[0107] In one embodiment of this application, the relationship between effective stress and true stress under uniaxial conditions is established based on the principle of equivalent strain:

[0108]

[0109] in, It is the undamaged stiffness tensor. For the true stress tensor, For strain tensor, For the effective stress tensor, It is the stiffness tensor. yes For isotropic materials in the analysis, the relationship between the stress tensor and the effective stress tensor is as follows:

[0110]

[0111] Among them, the fourth-order tensor The Voigt form is defined as:

[0112]

[0113] Define a special case where a uniaxial tensile load is applied to a material, where the effective stress tensor is defined as:

[0114]

[0115] in, These are the components of the effective stress tensor. This is a transpose.

[0116] When the three shear components are zero Combining the formula The associated variables per unit volume under uniaxial tension are obtained:

[0117]

[0118] in, This represents the complementary free energy density under uniaxial tension.

[0119] The three normal components of the true stress tensor are given by the following equation:

[0120]

[0121] Damage evolution law on internal variables The derivative is given by the following formula:

[0122]

[0123] The damage activation function defined for uniaxial tension is given by the following equation:

[0124]

[0125] Solving the zeros of the equation yields the expression for the effective stress:

[0126]

[0127] Based on this, according to the associated variables per unit volume, the effective strain tensor, and the derivatives of the damage variable and the internal variables under uniaxial tension, the scalar analytical expression is obtained:

[0128]

[0129] This is the final analytical expression derived from the previous one. Compared with the previous one, this expression simplifies the relationship between the actual stress and the effective stress. Although the formula is more complex overall, all of them are scalars, making it easier to solve.

[0130] S3: The scalar analytical expression is numerically integrated using the three-point Newton-Cotes integration algorithm, and the damage evolution parameters are iteratively solved using the secant method.

[0131] Furthermore, the step of numerically integrating the scalar analytical expression using the three-point Newton-Cotes integration algorithm and iteratively solving the damage evolution parameters using the secant method includes:

[0132] S401: Constructing the solution function based on scalar analytical expressions:

[0133] S402: Define the first two damage evolution parameters for the initial iteration;

[0134] S403: Based on the first two damage evolution parameters of the initial iteration, the three-point Newton-Cotes integral algorithm is used to solve the value of the improper integral in the solution function during the iteration process, and the function value of the solution function of the initial iteration is obtained;

[0135] S404: Based on the function value of the solution function in the initial iteration, the damage evolution parameters for the next iteration are obtained by iteratively solving the logarithmic secant method.

[0136] S405: Based on the damage evolution parameters of the next iteration, the three-point Newton-Cotes integral algorithm is used to solve the value of the improper integral in the solution function during the iteration process, and the function value of the solution function in the next iteration is obtained.

[0137] S406: Substitute the function value of the solution function for the next iteration into the solution function. If the absolute value of the solution function is less than the given integral tolerance, the iteration ends and the final damage evolution parameters are obtained; otherwise, update the damage evolution parameters of the iteration and proceed to S404 to continue the iterative solution.

[0138] In one embodiment of this application, the analytical expression derived by S2 is based on r Since the integral is an improper integral of the independent variable, it cannot be solved analytically directly. Therefore, numerical integration is used to solve it.

[0139] Without loss of generality, the solution function is written in the following form:

[0140]

[0141] in, Represents the anomalous integral, For the independent variable The relevant energy release rate, To and The relevant damage evolution parameters are determined by solving the above equation. For the above nonlinear equations, since their derivatives cannot be directly calculated, the secant method is used for iterative solution, such as... Figure 2 As shown, Figure 2 An iterative solution for A using the secant method provided in this application embodiment m A schematic diagram of the algorithm flow. The first two parameters to start the iteration process:

[0142]

[0143] in, and These represent the material's elastic modulus and tensile strength, respectively.

[0144] Damage parameters It only has physical meaning when it is greater than zero, therefore a logarithmic secant solution method is used:

[0145]

[0146] in, and Representing the Two initial values ​​for the damage evolution parameters in the next iteration. Representing the The damage evolution parameters obtained in the next iteration are also used as the initial values ​​of the damage evolution parameters in the (i+1)th iteration. Representing the The first initial value of the damage evolution parameter in the next iteration. Representing the The second initial value of the damage evolution parameter in the next iteration. and Representing the first The integral values ​​corresponding to the two initial values ​​in the next iteration.

[0147] Understandably, each iteration requires two initial values ​​for two points to form a secant line. The first iteration ( It was used and Create a secant line; the intersection of the secant line and the horizontal axis is the value obtained in the first iteration. Second iteration ( What is used is and And so on.

[0148] For one-dimensional numerical integrals like the one above, the three-point Newton-Cotes integration algorithm is used for solving. Its advantages include third-order algebraic accuracy and high computational efficiency. The solution process is as follows: Figure 3 As shown, this algorithm is used to calculate improper integrals. First, initialize the integration step. , The parameter value is set manually; in this embodiment, it is taken as 100. `count` is a loop counter. It is an internal variable. Let be the integration step size. The expression for the three-point Newton-Cotes integral is:

[0149]

[0150] in, Let be the integration step size, and for improper integrals, the lower limit of integration. =1, upper limit Assuming the stress becomes positive infinity, and that the stress becomes less than the tensile / compressive strength. KWhen the step size is multiplied by 1, the remaining energy is negligible. Therefore, the increment of the integration step size... h It can be defined as a function of the number of integration steps:

[0151]

[0152] in, The integration step size is... These are parameter values ​​set by the user. For tensile or compressive strength, Represents stretching. Represents compression. This is the integration step.

[0153] In one embodiment of this application, the above method is numerically solved using the UMAT and UEXTERNALDB subroutine modules of the finite element software ABAQUS. The UMAT subroutine calculates the constitutive response of the material, while the UEXTERNALDB subroutine calculates the damage parameters. As mentioned earlier, for each finite element element, The value is fixed, meaning it needs to be calculated before the incremental step begins. Therefore, in the UEXTERNALDB subroutine, LOP=0 is set to ensure that the entire finite element calculation process is performed only once. The damage evolution parameter A... m For the impact on grid dependency, please refer to Figure 4 The calculation results of out-of-plane shear for unidirectional fiber composite materials using RVE can be found in [reference]. Figure 5 .

[0154] This application provides a method for solving damage evolution parameters in a damage constitutive model. By transforming complex tensor operations into simple scalar analysis, it overcomes the difficulty of tensor integral solving in existing models, significantly reducing the computational burden. Furthermore, by employing the Newton-Cotes integration method and the secant method iterative solution strategy, it not only ensures the accurate calculation of damage parameters but also enhances the stability and reliability of the numerical algorithm, avoiding convergence problems common in traditional methods and significantly improving the efficiency and accuracy of solving damage evolution parameters. In addition, this method can accurately calculate the damage evolution parameter A based on the input fracture toughness parameters. m This method effectively solves the mesh dependency problem and provides correct input parameters for subsequent finite element strength calculations. Through this method, researchers can more accurately assess the mechanical response of materials and structures, providing theoretical and technical support for predicting the reliability of structures under service conditions.

[0155] It should be noted that those skilled in the art will recognize that the embodiments described herein are for the purpose of helping readers understand the principles of this application, and should be understood as not limiting the scope of protection of this application to such specific statements and embodiments. Those skilled in the art can make various other specific modifications and combinations based on the technical teachings disclosed in this application without departing from the essence of this application, and these modifications and combinations are still within the scope of protection of this application.

Claims

1. A method for solving damage evolution parameters in a damage constitutive model, characterized in that, include: S1: Based on the principles of thermodynamics, establish the mathematical relationship between the damage variable and the associated variable, and combine it with the damage evolution criterion to obtain the integral expression of the energy dissipation potential per unit volume. S2: Based on the principle of equivalent strain, the integral expression of the energy dissipation potential per unit volume is simplified to obtain the scalar analytical expression; Based on the principle of equivalent strain, the integral expression of the energy dissipation potential per unit volume is simplified to obtain a scalar analytical expression, including: S201: Considering the special case of applying a uniaxial tensile load to the material, a damage activation function is defined for the uniaxial tensile condition. The components of the effective stress tensor are obtained by finding the zero point of the damage activation function. The expressions for the damage activation function and the components of the effective stress tensor are as follows: in, For damage activation function, and These represent the tensile and compressive ultimate strengths of the matrix material, respectively. These are the components of the effective stress tensor. These are internal variables controlled by the laws governing damage evolution. Representative material type; S202: Based on the principle of equivalent strain, the true stress tensor is defined, and the relationship between the components of the effective stress tensor and the true stress tensor is established: in, For the effective stress tensor, Poisson's ratio of the matrix For transpose, , and The shear component of the true stress tensor. As a damage variable, For damage evolution parameters, These are internal variables controlled by the laws governing damage evolution; S203: Based on the relationship between the components of the effective stress tensor and the true stress tensor, To simplify the conditions, the associated variable per unit volume under uniaxial tension is obtained: in, To be related to damage variables Conjugate covariates, The complementary free energy density under uniaxial tension. The Young's modulus of the matrix. , and The shear component of the true stress tensor; S204: Based on the damage evolution law, the derivatives of the damage variable and the internal variables are obtained: S205: Based on the associated variables per unit volume, the effective strain tensor, and the derivatives of the damage variable and internal variables under uniaxial tension, the scalar analytical expression is obtained: In the formula, It is the energy release rate of the matrix. It is the length of the feature element. Represents the element in finite element calculation. Represents a break; S3: The scalar analytical expression is numerically integrated using the three-point Newton-Cotes integration algorithm, and the damage evolution parameters are iteratively solved using the secant method.

2. The method for solving the damage evolution parameters in the damage constitutive model according to claim 1, characterized in that, S1 includes: S101: Based on thermodynamic principles, the complementary free energy density of the damage constitutive model is defined; S102: Based on the complementary free energy density, obtain the associated variable: in, For complementary free energy densities, As a damage variable, For damage evolution parameters, These are internal variables controlled by the laws governing damage evolution. To be related to damage variables Conjugate covariates, The Young's modulus of the matrix. It is the shear stiffness of the material. , , , , and The shear component of the true stress tensor. Representative material type; S103: Based on the associated variables and complementary free energy densities, the energy dissipation potential per unit volume is obtained: in, The energy dissipation potential per unit volume, superscript Represents the derivative with respect to time; S104: Based on the characteristic length and fracture toughness of the finite element method, the energy dissipation potential per unit volume is regularized to obtain the integral expression for the energy dissipation potential per unit volume: in, Energy dissipated per unit volume It is the energy release rate of the matrix. It is the length of the feature element. For time, Represents the element in finite element calculation. It represents a break.

3. The method for solving the damage evolution parameters in the damage constitutive model according to claim 2, characterized in that, The method of numerically integrating the scalar analytical expression using the three-point Newton-Cotes integration algorithm and iteratively solving for the damage evolution parameters using the secant method includes: S401: Constructing the solution function based on scalar analytical expressions: S402: Define the first two damage evolution parameters for the initial iteration; S403: Based on the first two damage evolution parameters of the initial iteration, the three-point Newton-Cotes integral algorithm is used to solve the value of the improper integral in the solution function during the iteration process, and the function value of the solution function of the initial iteration is obtained; S404: Based on the function value of the solution function in the initial iteration, the damage evolution parameters for the next iteration are obtained by iteratively solving the logarithmic secant method. S405: Based on the damage evolution parameters of the next iteration, the three-point Newton-Cotes integral algorithm is used to solve the value of the improper integral in the solution function during the iteration process, and the function value of the solution function in the next iteration is obtained. S406: Substitute the function value of the solution function for the next iteration into the solution function. If the absolute value of the solution function is less than the given integral tolerance, the iteration ends and the final damage evolution parameters are obtained; otherwise, update the damage evolution parameters of the iteration and proceed to S404 to continue the iterative solution.

4. The method for solving the damage evolution parameters in the damage constitutive model according to claim 3, characterized in that, The solution function is: in, For improper integrals, For the independent variable The relevant energy release rate, To and Relevant damage evolution parameters.

5. The method for solving the damage evolution parameters in the damage constitutive model according to claim 4, characterized in that, The first two damage evolution parameters of the initial iteration are: in, and These represent the material's elastic modulus and tensile strength, respectively. and These are the two damage evolution parameters for the initial iteration.

6. The method for solving the damage evolution parameters in the damage constitutive model according to claim 5, characterized in that, The three-point Newton-Cotes integral algorithm is as follows: in, The integration step size is... The lower limit of integration, This is the maximum number of points. Let be the integrand.

7. The method for solving the damage evolution parameters in the damage constitutive model according to claim 6, characterized in that, The step size increment of the three-point Newton-Cotes integration algorithm is a function of the number of integration steps: in, The integration step size is... These are parameter values ​​set by the user. For tensile or compressive strength, Represents stretching. Represents compression. This is the integration step.

8. The method for solving the damage evolution parameters in the damage constitutive model according to claim 7, characterized in that, The damage evolution parameters for the next iteration are obtained by iteratively solving using the logarithmic secant method: in, Representing the The damage evolution parameters obtained in the second iteration are also used as the first... i Initial values ​​of damage evolution parameters for +1 iteration. Representing the The first initial value of the damage evolution parameter in the next iteration. Representing the The second initial value of the damage evolution parameter in the next iteration. and Representing the first The integral values ​​corresponding to the two initial values ​​in the next iteration.

Images