Efficient integral transformation cohesion fracture phase field simulation method and device

Abstract

The invention discloses an efficient integral transformation cohesion fracture phase field simulation method and device. The method is applied to finite element simulation software and comprises the following steps: establishing a geometric model of a material to be simulated, and performing mesh generation on the geometric model to obtain a mesh model; determining a constitutive relation between strain and stress based on the energy degradation function, the elastic stiffness tensor of the to-be-simulated material and inherent mechanical parameters; determining a control equation based on the dissipative function functional, the energy storage function and the thermodynamic law; obtaining a residual equation according to the control equation and the constitutive relation; the residual equation is linearized, and the control equation is converted into a matrix control equation; applying a load condition on the outer boundary of the grid model, and solving to obtain a phase field distribution result by gradually loading a load step based on a matrix control equation; for each load step, it is determined to complete the iteration of the load step when convergence is reached. The method can be implemented in ABAQUS software, and is efficient and suitable for large-scale engineering application.

Description

Technical Field

[0001] This invention relates to the field of computational mechanics, particularly to the field of fracture mechanics simulation, and specifically to an efficient integral transform cohesive fracture phase field simulation method and apparatus. Background Technology

[0002] The fracture phase-field method has become a promising variational framework for simulating complex crack problems. By solving for the auxiliary variable—the phase-field parameter—the fracture problem can be re-considered as an energy minimization problem. Complex fracture features, such as crack branching, crack initiation at arbitrary locations, or the merging of multiple cracks, can be naturally captured in the original finite mesh. With the continuous development of the fracture phase-field method theory, researchers have also devoted considerable effort to developing effective schemes for solving displacement-phase-field coupled problems. Mathematically, it can be proven that the total potential energy of such coupled problems is not globally convex. It is precisely because of this convexity that the uniform tangent stiffness matrix in Newton's iteration method becomes indeterminate, the overall solution algorithm becomes unstable and has poor convergence, and the energy minimization problem remains extremely challenging.

[0003] In fact, for this kind of non-global convex problem, the subspace staggered iteration scheme is more popular. That is, by fixing the original displacement field (i.e., phase field), the total potential energy becomes a convex function with respect to the phase field variable (i.e., displacement field variable). This method has been verified as a robust method, but the problem it faces is that the computational cost is huge, especially in some key steps, where the number of iterations may reach thousands.

[0004] With the further development of engineering and scientific problems, the demand for solving complex cohesion laws is becoming increasingly significant. Theoretically, cohesion curves of arbitrary shapes can be realized through integral transform phase-field models. However, due to the complexity of the theoretical formulas, many challenges are faced in practical numerical implementation. Therefore, developing efficient numerical solution algorithms remains an urgent task for fracture phase-field simulation. Another obstacle limiting the use of phase-field damage models in practical problems is that their numerical simulation methods are highly dependent on internal code. Currently, integral transform phase-field models are only used in the finite element software package Fenics. The automatic differentiation function of this package avoids the complex process of splicing stiffness matrices in finite element calculations, and users only need to provide the weak form of the governing equations. However, this results in its inapplicability to large-scale engineering applications. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides an efficient integral transform cohesive fracture phase field simulation method and apparatus, which can be implemented in the general-purpose finite element simulation software ABAQUS, thereby improving the computational efficiency and accuracy of fracture phase field simulation in complex structures and large-scale engineering applications.

[0006] This invention provides an efficient method and apparatus for simulating the cohesive fracture phase field using integral transform. The technical solution is as follows: On the one hand, an efficient integral transform cohesive fracture phase field simulation method is provided and applied to finite element simulation software. The method includes: A geometric model of the material to be simulated is established, and the geometric model is meshed to obtain a mesh model; The constitutive relationship between strain and stress is determined based on the energy degradation function, the elastic stiffness tensor of the material to be simulated, and the inherent mechanical parameters; and the governing equations are determined based on the dissipation function functional, the storage function, and the laws of thermodynamics; and the residual equations are obtained according to the governing equations and the constitutive relationship; wherein the inherent mechanical parameters include fracture energy, intrinsic characteristic length, Young's modulus, tensile strength, shear modulus, and shear strength. The residual equation is linearized to transform the governing equation into a matrix governing equation; the matrix governing equation is used to describe the relationship between the variation of the displacement degree of freedom, the variation of the phase field degree of freedom, and the residual through the stiffness matrix; Load conditions are applied to the outer boundary of the mesh model, and the phase field distribution is obtained by solving the matrix control equation step by step by loading the load. For each load step, the iteration of the load step is completed when convergence is reached.

[0007] On the other hand, an efficient integral transform cohesive fracture phase field simulation device is provided, the device comprising: The first construction module is used to establish a geometric model of the material to be simulated, and to mesh the geometric model to obtain a mesh model; The second construction module is used to determine the constitutive relationship between strain and stress based on the energy degradation function, the elastic stiffness tensor of the material to be simulated, and the inherent mechanical parameters; and to determine the governing equations based on the dissipation function functional, the storage function, and the thermodynamic laws; and to obtain the residual equations according to the governing equations and the constitutive relationship; wherein the inherent mechanical parameters include fracture energy, intrinsic characteristic length, Young's modulus, tensile strength, shear modulus, and shear strength; The transformation module is used to linearize the residual equation and convert the control equation into a matrix control equation; the matrix control equation is used to describe the relationship between the variation of the displacement degree of freedom, the variation of the phase field degree of freedom, and the residual through the stiffness matrix; The simulation solution module is used to apply load conditions on the outer boundary of the mesh model and obtain the phase field distribution results by progressively loading load steps based on the matrix control equations; wherein, for each load step, when convergence is reached, it is determined that the iteration of the load step is completed.

[0008] On the other hand, a computer device is provided, the computer device including a memory and a processor, the memory for storing computer programs, and the processor for executing the computer programs stored in the memory to implement the steps of the efficient integral transform cohesive fracture phase field simulation method described above.

[0009] On the other hand, a computer-readable storage medium is provided, wherein a computer program is stored therein, and when the computer program is executed by a processor, the steps of the efficient integral transform cohesive fracture phase field simulation method described above are implemented.

[0010] On the other hand, a computer program product is provided, including a computer program that, when executed by a processor, implements the steps of the efficient integral transform cohesive fracture phase field simulation method described above.

[0011] The technical solution provided by this invention can bring at least the following beneficial effects: This invention, based on the laws of thermodynamics and combined with the phase-field characterization of cracks, systematically derives the displacement field equilibrium equation and the integral transform phase-field evolution equation (i.e., the governing equations) within a unified thermodynamic framework. This allows for the gradual application of load steps to obtain the phase-field distribution results. In terms of numerical implementation, this invention prioritizes efficiency, utilizing the built-in UEL subroutine of ABAQUS to numerically solve the integral transform phase-field model. This significantly improves both solution stability and convergence efficiency while drastically reducing the number of iterations, further enhancing computational efficiency. Thus, this implementation framework not only ensures accuracy but also offers higher computational efficiency, laying a solid foundation for the widespread application of the integral transform phase-field model in complex engineering scenarios. It demonstrates good scalability and engineering applicability, particularly suitable for the detailed simulation of large-scale crack propagation and complex failure processes. Attached Figure Description

[0012] To more clearly illustrate the technical solutions in the embodiments of the present invention 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 some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0013] Figure 1 This is a flowchart of an efficient integral transform cohesive fracture phase field simulation method provided by an embodiment of the present invention; Figure 2 This is a schematic diagram of a solid control body provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the geometric model and boundary conditions of an L-shaped concrete slab provided in an embodiment of the present invention; Figure 4 This is a schematic diagram of the phase field distribution after fracture of an L-shaped concrete slab according to an embodiment of the present invention; Figure 5 This is a comparison chart of load-displacement curves using an integral transform phase-field model and a cohesive phase-field model, provided by an embodiment of the present invention. Figure 6 This is a comparison chart of the computational efficiency results of the integral transform phase field model and the cohesive phase field model provided in an embodiment of the present invention; Figure 7 This is a hardware architecture diagram of a computer device provided in an embodiment of the present invention; Figure 8 This is a structural diagram of an efficient integral transform cohesive fracture phase field simulation device provided by an embodiment of the present invention. Detailed Implementation

[0014] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0015] The following is the specific concept of the present invention.

[0016] Please refer to Figure 1 This invention provides an efficient integral transform cohesive fracture phase field simulation method, applied to finite element simulation software. The method includes: Step 100: Establish a geometric model of the material to be simulated, and mesh the geometric model to obtain a mesh model; Step 102: Determine the constitutive relationship between strain and stress based on the energy degradation function, the elastic stiffness tensor of the material to be simulated, and the inherent mechanical parameters; determine the governing equations based on the dissipation function functional, the storage function, and the laws of thermodynamics; and obtain the residual equations based on the governing equations and the constitutive relationship; wherein, the inherent mechanical parameters include fracture energy, intrinsic characteristic length, Young's modulus, tensile strength, shear modulus, and shear strength; Step 104: Linearize the residual equations to transform the governing equations into matrix governing equations; the matrix governing equations are used to describe the relationship between the variational values ​​of the displacement degrees of freedom, the variational values ​​of the phase field degrees of freedom, and the residuals through the stiffness matrix. Step 106: Apply load conditions to the outer boundary of the mesh model, and obtain the phase field distribution results by solving the matrix control equations step by step by loading loads; wherein, for each load step, when convergence is reached, the iteration of that load step is completed.

[0017] In this embodiment of the invention, based on the laws of thermodynamics and combined with the phase-field characterization of cracks, the displacement field equilibrium equation and the integral transform phase-field evolution equation (i.e., the governing equations) are systematically derived within a unified thermodynamic framework. This allows for the gradual application of load steps to obtain the phase-field distribution results. In terms of numerical implementation, this invention prioritizes efficiency, utilizing the built-in UEL subroutine of ABAQUS to numerically solve the integral transform phase-field model. This significantly improves both solution stability and convergence efficiency while drastically reducing the number of iterations, further enhancing computational efficiency. Thus, this implementation framework not only ensures accuracy but also offers higher computational efficiency, laying a solid foundation for the widespread application of the integral transform phase-field model in complex engineering scenarios. It demonstrates good scalability and engineering applicability, particularly suitable for the detailed simulation of large-scale crack propagation and complex failure processes.

[0018] The following description Figure 1 The execution method for each step is shown.

[0019] First, for step 100, a geometric model of the material to be simulated is established, and the geometric model is meshed to obtain a mesh model. Specifically, the geometric model is a two-dimensional or three-dimensional geometric model of the material to be simulated. A continuous domain is divided using multi-field elements (taking two-dimensional as an example, two degrees of freedom are used for the displacement field, and one degree of freedom is used for the phase field) to form a finite element computational domain, which contains a number of nodes and a number of elements.

[0020] In a preferred embodiment, in step 102, the constitutive relation is determined by the following formula: (1) in, For stress tensor; It is the energy degradation function; For effective stress tensor; Let be the elastic stiffness tensor of the material to be simulated; For strain tensor; The energy degradation function is determined by the following formula: (2) (3) (4) in, These are intermediate material parameters; These are eigenfunctions of energy density; This is the fracture energy; The intrinsic feature length; Young's modulus; Tensile strength; Shear modulus; Shear strength; d This represents the phase field value.

[0021] In a preferred embodiment, in step 102, the governing equations are determined based on the dissipation function functional, the storage function, and the laws of thermodynamics, including: Based on the variational relationships of fracture energy and crack density function, determine the variational relationship of dissipative function functional. , represented as: (5) in, This is the fracture energy; For the variation of the crack density function; , This is the cracked area; The material to be simulated is a solid system; The energy storage function of the solid system of the material to be simulated is transformed into a function of the strain tensor and phase field variables, and the variation of the energy storage density functional is determined based on the stress tensor and the damage energy release rate. , represented as: (6) in, For stress tensor; It is a symmetric gradient operator; For the variation of the displacement field; The rate of energy release from damage; The variation of the crack phase field; Substituting the variational functions of the dissipation function functional, the storage density functional, and the virtual work of external forces expressed by volume forces and surface traction forces into the laws of thermodynamics, we obtain the first equation. Using the divergence theorem, the first equation is transformed into the governing equation, expressed as: (7) (8) in, This is the damage flux vector; Net damage source per unit volume; For stress tensor; , These are the body force and surface traction force, determined by the load step, respectively. This is the cracked area; The material to be simulated is a solid system; For the force boundary of a solid system, ; The normal vector of the solid system; This is the normal vector of the crack region; For Hamiltonian operators.

[0022] In a preferred embodiment, the mesh model includes several multi-field elements and several nodes; in step 102, the residual equation is obtained based on the governing equations and constitutive relations, including: Based on the mesh type and number of nodes of the multi-field element, the interpolation function matrix of the displacement field, the interpolation function matrix of the phase field, the strain-displacement compatibility matrix, and the phase field gradient compatibility matrix are determined respectively. Substituting the interpolation function matrix of the displacement field, the interpolation function matrix of the phase field, the strain-displacement compatibility matrix, and the phase field gradient compatibility matrix into the governing equations yields the residual equations; the residual equations are determined by the following formula: (9) in, For the residuals of the displacement field; The residual of the phase field; and These are the interpolation function matrices for the displacement field and the phase field, respectively. , These are the body force and surface traction force, determined by the load step, respectively. , These are the strain-displacement compatibility matrix and the phase-field gradient compatibility matrix, respectively. For stress tensor; The finite element computational domain of the mesh model. The boundary of the finite element computational domain, The finite element computational domain is defined for the crack region. Let be the derivative of the energy degradation function with respect to space; This represents the historical maximum effective damage energy release rate. This is the fracture energy; The intrinsic feature length; The derivative of the phase field profile function with respect to space; d This refers to the phase field value; For Hamiltonian operators.

[0023] In a more preferred embodiment, each load step further includes: The equivalent effective stress is calculated using the stress tensor from the previous iteration. The effective energy damage release rate of the previous iteration step is calculated based on the equivalent effective stress and Young's modulus. The critical effective energy damage release rate is calculated based on tensile strength and Young's modulus. The maximum value between the effective energy damage release rate and the critical effective energy damage release rate of the previous iteration step is selected as the historical maximum value of the effective damage energy release rate of the current iteration step. The historical maximum effective damage energy release rate was determined by the following formula: (10) (11) (12) (13) in, This represents the historical maximum effective damage energy release rate. The critical effective energy damage release rate; The effective energy damage release rate of the previous iteration step; This is the equivalent effective stress; Young's modulus; Tensile strength; , Uniaxial compressive strength; This is the first principal stress; It is the second invariant; For the Macauly operator; the effective stress tensor in formula (1) The first principal stress and the second invariant are derived from the stress tensor of the previous iteration step. Sure; Specifically, for the equivalent effective stress, the first principal stress and the second invariant of the effective stress tensor are first determined by the stress tensor of the previous iteration step; the first principal stress, the second invariant, and the uniaxial compressive strength and tensile strength of the material to be simulated are substituted into formula (13) to calculate the equivalent effective stress; among them, the inherent mechanical parameters also include uniaxial compressive strength and tensile strength.

[0024] It should be noted that when the current iteration number is 1, the stress tensor of the previous iteration step is the final stress tensor output by the previous load step.

[0025] Specifically, based on the mesh model, a solid control volume of the material to be simulated is selected, such as... Figure 2 The image shows a diffuse crack. solid systems In physical strength Mechanical description under action and boundary constraints. The outer boundary and normal vector of the solid system are respectively represented as: and The outer boundary can be represented as a non-intersecting displacement boundary. Force Boundary Displacement loads are applied to the displacement boundaries. Surface traction force Apply to the force boundary. Represents the coordinates of various locations in space, displacement field Characterizing the motion and deformation of a solid system. Under the assumption of small deformation, the strain field of the solid system can be expressed as: ,in It is a symmetric gradient operator, representing the gradient with respect to the spatial system; The displacement field and its variational space are determined by the following formula: (14) For the space of the displacement field; Let be the variational space of the displacement field; For the variation of the displacement field; Within the framework of phase-field theory, cracks are also considered as a phase distinct from intact materials, i.e., a continuously changing quantity, namely the phase-field variable, is introduced. Describe its state (the phase field value will be used in the following text) d (This indicates that) the sharp crack shown is diffused into a crack band with a certain width. , These represent the undamaged and fully damaged states of the solid, respectively. Therefore, the space of the crack phase field and its variation is determined by the following formula: (15) in, The space of the crack phase field; For the variational space of the crack phase field; It is the variation of the crack phase field.

[0026] In this invention, the phase field evolution process is described using thermodynamics, and the first and second laws of thermodynamics are as follows: (16) External force virtual work Represented by the sum of volume force and surface traction force: (17) The crack regularization process is determined based on the crack density function through a geometric description of the crack region, where the crack density function is: : (18) in, The intrinsic feature length; These are the energy density eigenfunctions, which are the core of the integral transform phase field; This is the Hamiltonian operator. Since crack propagation is essentially an energy dissipation process, meaning it is irreversible, that is: (19) Thus, the crack density function is a monotonically non-decreasing function, which fully characterizes the irreversibility of the crack; Cracked area The actual crack area. In this invention, rate-dependent constitutive behavior is not considered, and the dissipative function functional... It is given by the following formula: (20) Therefore, the variation of the dissipative function functional is determined by formula (5); the variation of the crack density function Determined by the following formula: (twenty one) in, This is called the phase-field profile function. Without loss of generality, the energy function for a solid-state system is... Determined by the following formula: (twenty two) in, That is, displacement field Local energy density functional The following formula is used to determine it: (twenty three) Among them, in formula (23) The energy degradation function characterizes the reduction of the initial strain energy as the phase field evolves. Under the framework of the integral transformation phase field theory, it is determined by the following formula (2), where the intermediate material parameters are calculated by substituting the inherent mechanical parameters into formula (3), and the energy density eigenfunction is obtained by formula (4). In formula (23) The initial free energy of a material in its undamaged state is represented as follows: For an elastic material, it is expressed as: (twenty four) in, For effective stress tensor; These are the elastic stiffness tensor and flexibility tensor of the material to be simulated, respectively.

[0027] In this embodiment of the invention, the variation of the energy storage density functional is... It is given by the following formula: (6) And it has the following constitutive relations, (25) in, For stress tensor; This is the first derivative of the energy degradation function with respect to the phase field variables; The thermodynamic driving force, which is the dual to the phase field variable, is generally referred to as the damage energy release rate. To effectively reduce the energy release rate of damage; The constitutive relation (i.e., equation (1)) and the thermodynamic driving force related to phase field evolution are defined as follows: (26) Among them, effective damage energy release rate Elastic strain energy of representative non-destructive materials: (27) The constitutive relations described above do not distinguish between asymmetric mechanical behaviors under tensile and compressive stress states, resulting in crack initiation and evolution that do not conform to physical reality. To address this issue, formulas (12) and (13) are used to consider the different mechanical behaviors under tensile and compressive stress states.

[0028] In this embodiment of the invention, the constitutive relation can still be used for brittle and quasi-brittle materials, that is, formula (26) can still be used, but the effective energy damage release rate needs to be redefined, as shown in formula (12). Thus, by substituting the variational formula (5) of the dissipative function functional and the variational formula (6) of the storage density function functional into the thermodynamic law (16), and combining it with the expression (17) of the virtual work of external forces, we obtain: (28) (29) Formula (28) is derived from formula (17). Substituting formula (28), formula (5) and formula (6) into formula (16) yields formula (29). Formulas (28) and (29) are the first equation. Using the divergence theorem, the above weak form can be written as: (30) (31) Therefore, based on formulas (30) and (31), the control equation of the system can be written as: (32) (33) The variation of the energy density function can be written as: , If the operator is Laplace, then the natural boundary conditions are: (34) in, Indicates the boundary of the crack region The unit outward normal vector. For the convenience of subsequent numerical implementation, the strong form governing equations of the phase field and displacement field are written in a unified form, as shown in equations (7) and (8); similar to the form of the heat conduction equation, the damage flux vector Net damage source per unit volume This can be expressed by the following formula: (35) (36) It can be seen from formulas (35), (36) and (2) that the energy density eigenfunction establishes the coupling relationship between the displacement field equilibrium equation and the phase field evolution equation; In this embodiment of the invention, to ensure the irreversibility of the crack evolution process, a phenomenological maximum historical variable is constructed as the phase field evolution threshold, theoretically guaranteeing the physical rationality of damage evolution. Therefore, the historical maximum value of the effective damage energy release rate is introduced. Replace the effective damage energy release rate in the above formula (36) As shown in formula (10), the critical effective energy damage release rate is In order to meet the conditions It can be seen that the material's breaking strength is the key parameter controlling crack initiation. At this point, substituting formula (10) into formula (35), the net damage source per unit volume is equivalent to: In the governing equations, the only unknown function to be determined is the derivative of the phase field profile function. It is related to the energy density eigenfunction. The cohesive eigenfunction and the energy density eigenfunction in the formula can be linked by the Feng-Li integral transform, but the present invention adopts the linear cohesive law, so the energy density eigenfunction is as shown in formula (4); According to the weighted residual method, under the above spatial constraints, the weak form of the governing equations can be written as: (37).

[0029] In this embodiment of the invention, specifically, a feature size of [missing information] is used. Multi-field elements (taking two dimensions as an example, with two degrees of freedom for the displacement field and one degree of freedom for the phase field) divide the continuous domain To form the finite element computational domain , which contains Each node and Each element. Based on the mesh type and number of nodes of the multi-field element, the displacement vectors of all element nodes within the finite element computational domain are calculated. and phase field variables The vectors are respectively (38) Based on this, the displacement and phase field distributions of the nodes are interpolated to the element integration points using interpolation functions. Since the functions are continuous, they can be solved using Gaussian integration at the integration points. Therefore, the displacement and phase fields in the finite element computational domain are obtained through nodal variable interpolation, i.e. (39) and These are the interpolation function matrices for the displacement field and the phase field, respectively. and It is the interpolation function related to node I. It is a unit shape function. It is a unit vector; It is the nodal displacement vector of node I. These are the nodal phase field values, and the global displacement vector of the finite element computational domain can be obtained through global assembly. and the overall phase field vector The displacement field in the computational domain can be obtained through the above interpolation. and phase field ; Correspondingly, strain field and phase field gradient field It can be represented as (40) In the formula, The strain-displacement compatibility matrix, The phase field gradient coordination matrix; Then, the above finite element discretization schemes (38) to (40) are substituted into the system's control equations and constitutive relations to obtain the residual equations as shown in formula (9).

[0030] In a preferred embodiment, for step 104, the matrix control equation is determined by the following formula: (41) in, This is the first component; This is the second component; For the variation of displacement degrees of freedom; For the variation of the phase field degrees of freedom; For the residuals of the displacement field; This represents the residual of the phase field.

[0031] In this invention, the aforementioned coupled nonlinear equations can be solved using the finite element simulation software ABAQUS. This method is simple and easy to implement, and does not require writing multiple subroutines. In this method, the global quasi-Newton iteration method is used to linearize equation (9), where the unknowns (i.e., variables, including displacement and phase field values) can be uniformly written as... The governing equations can be written as: (42) The components of the stiffness matrix can be written as: (43) Since the integrands in formula (43) are all continuous, these matrices can be calculated using the standard Gaussian quadrature formula. In this invention, in order to linearize the equations and transform the above strongly coupled problem into a weakly coupled problem, ignoring the inter-field coupling can eliminate the off-diagonal matrix, thus obtaining formula (41).

[0032] In a preferred embodiment, the overall control equations are determined by the following formula: (44) in, Let $\mathbf{k}$ be the global stiffness matrix for the $k$-th iteration. Let this be the global variable for the k-th iteration; For the (k-1)th iteration, this is the global variable; The global residual for the k-th iteration; This represents the overall residual from the (k-1)th iteration. It should be noted that the overall variable increment... The actual variables include displacement increments and phase field increments, but since stress increments can be calculated from displacement increments, the overall variable increments here include both stress increments and phase field increments. This overall governing equation is used to describe the relationship between the overall variables and the overall residuals through the overall stiffness matrix.

[0033] Specifically, in the quasi-Newton iterative framework, the residual equations are represented by g. When the solution to the nonlinear residual equations approaches 0, the original stiffness matrix is ​​replaced by the global stiffness matrix, resulting in formula (44), to reduce computational overhead and improve convergence efficiency. Then, the global stiffness matrix is ​​updated by a rank-2 correction matrix, resulting in the corrected equation: (45) Wherein, the superscript -1 is used to indicate finding the inverse matrix, the superscript T is used to indicate taking the transpose matrix, and k is the number of iterations; This represents the overall change in variables, specifically the increment of the overall variables in the previous iteration (i.e., the (k-1)th iteration). The change in the overall residual is the change between the overall residual of the current iteration (i.e., the kth iteration) and the overall residual of the previous iteration (i.e., the (k-1th iteration)).

[0034] In step 106, the phase field distribution results are obtained by solving the matrix control equations through stepwise loading of loads, including: The matrix control equation (Equation (41)) is transformed by the overall stiffness matrix to obtain the overall control equation (Equation (44)); the overall control equation is modified to obtain the modified equation (Equation (45)). For each load step, execute: S0, the overall residual is calculated based on the initial phase field value and the initial stress tensor; wherein, the overall residual includes the residual of the displacement field and the residual of the phase field, and is calculated by the residual equation (formula (9)); S1, substitute the overall variable increment, overall residual change and overall stiffness matrix of the previous iteration into the correction equation (formula (45)) to obtain the overall stiffness matrix of the current iteration; S2, substitute the overall residual obtained in step S0, the overall residual of the previous iteration, and the overall stiffness matrix of the current iteration obtained in step S1 into the overall control equation (formula (44)) to calculate the overall variable increment of the current iteration; where the overall variable increment includes stress increment and phase field increment; S3, based on the stress increment and the initial stress tensor, obtain the first stress tensor; and based on the phase field increment and the initial phase field value, determine the first phase field value; S4, Substitute the first stress tensor and the first phase field value into the residual equation (formula (9)) to obtain the current overall residual; S5. Determine whether convergence has been achieved based on the current overall residual; if yes, end the iteration of this load step; otherwise, proceed to step S6. S6, increment the iteration count by 1, and based on the first phase field value and the first stress tensor, update the initial phase field value and the initial stress tensor for the current iteration, and return to step S0.

[0035] It should be noted that when the current iteration number is 2, the first and second components of the stiffness matrix are calculated using formula (43), and this stiffness matrix is ​​the global stiffness matrix of the previous iteration. In step S3, the updated solution is calculated using a scaling factor, i.e. ; , These are the initial variables for the k-th and (k-1)-th iterations, namely the initial stress tensor and the initial phase field value; 's' represents the overall variable change in the current iteration, i.e., the overall variable increment obtained in step S2; 's' is the scaling factor, which can be automatically provided by the ABAQUS software. For example, the first stress tensor is the sum of the product of the stress increment and the scaling factor and the initial stress tensor; the first phase field value is the sum of the product of the phase field increment and the scaling factor and the initial phase field value. In step S4, the first stress tensor, the body force and surface traction force determined by this load step are substituted into the residual equation to obtain the residual of the displacement field; and the inherent mechanical parameters, the first phase field value and the first stress tensor are substituted into the residual equation to obtain the residual of the phase field.

[0036] It should be noted that when k=1, the phase field value and stress tensor output from the previous load step are used as the initial phase field value and initial stress tensor for the first iteration. The phase field distribution results output in step 106 include the phase field distribution diagram and stress-strain curve (or load-displacement curve).

[0037] In this invention, the line search method in step S3 further avoids the divergence of equilibrium iteration caused by the inaccuracy of the stiffness matrix, so that the value of the residual equation in the search direction is zero within a certain tolerance range.

[0038] In one specific implementation, determining whether the convergence condition has been met includes: Determine whether the product of the total time-averaged flux of the displacement and damage currently achieved in this loading step and the first coefficient is not less than the current overall residual (i.e., the maximum value of the residual in the current displacement field or the field). Determine whether the product of the maximum increment change of the variable in the current loading step (i.e., the maximum value of the overall variable increment) and the second coefficient is not less than the maximum correction value of the variable in the current iteration step; If all the judgment results are yes, then the convergence condition has been met.

[0039] Specifically, the convergence condition is determined by the following formula: (46) in, It is the maximum value among the residuals of the current displacement field and the residuals of the current phase field in the current overall residual; The total time-averaged flux of displacement and damage achieved so far in this loading step; This is the first coefficient, used to represent the set tolerance, for example, 0.005; (47) in, This represents the maximum correction value of the variable in the current iteration step. This represents the maximum incremental change of the variable in the current loading step. The second coefficient represents the set tolerance, for example, 0.01; the variable here refers to the stress increment or phase field increment.

[0040] In this invention, the quasi-Newton iterative algorithm is applied for the first time to the complex mechanical model of the integral transform phase field, and the feasibility and accuracy of the method are verified. Furthermore, in the UEL subroutine of the general-purpose finite element simulation software ABAQUS, the stiffness matrix and residual vector of the element (given by formulas (9), (43), and (45)) are given. However, it should be noted that the above formulas need to be summed at each integration point in the UEL subroutine; specific details are not elaborated here. It should be noted that since the integral transform phase field is strongly nonlinear, exceeding the number of iterations required for conventional finite element calculations, it is necessary to modify the default control parameters built into the ABAQUS software.

[0041] It should be noted that the same characters in all the above formulas have the same meaning; the " " above the character ", located above and to the right of the characters Both ":" and ":=" are used to indicate differentiation of a function; the ":" on the right side of the equals sign is used to indicate the double dot product of a tensor; ":=" is used to indicate a user-defined expression, which is equivalent to the equals sign.

[0042] In one specific embodiment, a hybrid failure test of an L-shaped concrete slab is used as an example. The feasibility and accuracy of the method provided by this invention are verified through tensile numerical simulation calculations of the specimen. The geometric dimensions and boundary conditions of the specimen are as follows: Figure 3 As shown, F and u represent the applied load conditions and loading point, respectively. F represents the applied external force and u represents the applied displacement, corresponding to the force control and displacement control methods in numerical simulation. All values ​​are in mm. Although the L-shaped specimen is subjected to simple forces, the numerical simulation of its failure mode remains a highly challenging problem. First, the material parameters of the specimen are: fracture energy of 0.09 N / mm², tensile strength of 2.7 MPa, Young's modulus of 25.85 GPa, and Poisson's ratio of 0.18. Figure 4 The numerically predicted failure modes of the L-shaped specimen are shown. Figure 5 The model of this invention (i.e., the integral transform phase field model) is compared with the cohesive phase field model proposed in the prior art. The force-displacement curves at the loading point of the specimen obtained by numerical simulation show that the calculation results of the two models are almost identical overall, with only slight differences near the peak value, thus verifying the correctness of the method. Figure 6 The computational efficiency of these two phase-field models is further illustrated. It can be seen that the integral transform phase-field model is significantly more efficient overall than the cohesive phase-field model, and the number of iterations for the integral transform phase-field model does not exceed 20, demonstrating superior performance. It should be noted that... Figure 3In the figure, F and u represent the applied load conditions, where F represents the applied external force and u represents the applied displacement, corresponding to the force control and displacement control forms in numerical simulation, respectively.

[0043] like Figure 7 , Figure 8 As shown, this embodiment of the invention provides a highly efficient integral transform cohesive fracture phase field simulation device. The device embodiment can be implemented through software, hardware, or a combination of both. From a hardware perspective, as... Figure 8 The diagram shown is a hardware architecture diagram of a computing device housing a high-efficiency integral transform cohesive fracture phase field simulation device provided in an embodiment of the present invention. (Except for...) Figure 8 In addition to the processor, memory, network interface, and non-volatile memory shown, the computing device in the embodiment may also include other hardware, such as a forwarding chip responsible for processing packets. Taking software implementation as an example, such as... Figure 8 As shown, as a logical device, it is formed by the CPU of its computing device reading the corresponding computer program from the non-volatile memory into memory and running it. This embodiment provides a high-efficiency integral transform cohesive fracture phase field simulation device, comprising: The first construction module 800 is used to establish the geometric model of the material to be simulated and to mesh the geometric model to obtain a mesh model; The second construction module 802 is used to determine the constitutive relationship between strain and stress based on the energy degradation function, the elastic stiffness tensor of the material to be simulated, and the inherent mechanical parameters; and to determine the governing equations based on the dissipation function functional, the storage function, and the laws of thermodynamics; and to obtain the residual equations according to the governing equations and the constitutive relationship; wherein the inherent mechanical parameters include fracture energy, intrinsic characteristic length, Young's modulus, tensile strength, shear modulus, and shear strength; The transformation module 804 is used to linearize the residual equations and convert the governing equations into matrix governing equations. The matrix governing equations are used to describe the relationship between the variational values ​​of the displacement degrees of freedom, the variational values ​​of the phase field degrees of freedom, and the residuals through the stiffness matrix. The simulation solution module 806 is used to apply load conditions on the outer boundary of the mesh model and obtain the phase field distribution results by solving the matrix control equations step by step by loading loads. Specifically, for each load step, the iteration of the load step is completed when convergence is reached.

[0044] In some specific implementations, the first construction module 800 can be used to execute the above step 100, the second construction module 802 can be used to execute the above step 102, the conversion module 804 can be used to execute the above step 104, and the simulation solution module 806 can be used to execute the above step 106.

[0045] Since the contents of the above-described apparatus are based on the same concept as the method embodiments of the present invention, the specific contents can be found in the descriptions in the method embodiments of the present invention, and will not be repeated here.

[0046] It is understood that the structures illustrated in the embodiments of the present invention do not constitute a specific limitation on an efficient integral transform cohesive fracture phase-field simulation device. In other embodiments of the present invention, an efficient integral transform cohesive fracture phase-field simulation device may include more or fewer components than illustrated, or combine some components, or split some components, or have different component arrangements. The illustrated components may be implemented in hardware, software, or a combination of software and hardware.

[0047] The information interaction and execution process between the modules in the above-mentioned device are based on the same concept as the method embodiment of the present invention, and the specific details can be found in the description of the method embodiment of the present invention, and will not be repeated here.

[0048] This invention also provides a computing device, including a memory and a processor. The memory stores a computer program, and when the processor executes the computer program, it implements an efficient integral transform cohesive fracture phase field simulation method according to any embodiment of this invention.

[0049] This invention also provides a computer-readable storage medium storing a computer program, which, when executed by a processor, causes the processor to perform an efficient integral transform cohesive fracture phase field simulation method according to any embodiment of this invention.

[0050] Embodiments of this application also provide a computer program product, which includes a computer program. A processor of a computer device reads the computer program from a computer-readable storage medium and executes the computer program, causing the computer device to perform an efficient integral transform cohesive fracture phase field simulation method as described in any of the above embodiments.

[0051] Specifically, a system or apparatus equipped with a storage medium may be provided, on which software program code implementing the functions of any of the embodiments described above is stored, and the computer (or CPU or MPU) of the system or apparatus may read and execute the program code stored in the storage medium.

[0052] In this case, the program code read from the storage medium can itself implement the function of any of the above embodiments, and therefore the program code and the storage medium storing the program code constitute part of the present invention.

[0053] Storage media embodiments for providing program code include floppy disks, hard disks, magneto-optical disks, optical disks (such as CD-ROM, CD-R, CD-RW, DVD-ROM, DVD-RAM, DVD-RW, DVD+RW), magnetic tapes, non-volatile memory cards, and ROMs. Alternatively, program code can be downloaded from a server computer via a communication network.

[0054] Computer-readable signal media may include data signals propagated in baseband or as part of a carrier wave, carrying computer-readable program code. Such propagated data signals may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. Computer-readable signal media may also be any computer-readable medium other than computer-readable storage media, which can send, propagate, or transmit programs for use by or in conjunction with an instruction execution system, system, or device.

[0055] The program code contained on a computer-readable medium may be transmitted using any suitable medium, including, but not limited to, wireless, wire, optical fiber, RF, etc., or any suitable combination thereof.

[0056] Computer program code for performing the operations of this invention can be written in one or more programming languages ​​or a combination thereof, including object-oriented programming languages ​​such as Java, Smalltalk, and C++, as well as conventional procedural programming languages ​​such as "C" or similar programming languages. The program code can be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computer (e.g., via the Internet using an Internet service provider).

[0057] Furthermore, it should be clear that not only can the program code read by the computer be executed, but also the operating system or other components operating on the computer can be instructed based on the program code to perform some or all of the actual operations, thereby realizing the function of any of the embodiments described above.

[0058] Furthermore, it is understood that the program code read from the storage medium is written to the memory set in the expansion board inserted into the computer or to the memory set in the expansion module connected to the computer. Then, based on the instructions of the program code, the CPU or other components installed on the expansion board or expansion module execute some and all of the actual operations, thereby realizing the function of any of the above embodiments.

[0059] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitations, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element.

[0060] Those skilled in the art will understand that all or part of the steps of the above method embodiments can be implemented by hardware related to program instructions. The aforementioned program can be stored in a computer-readable storage medium. When the program is executed, it performs the steps of the above method embodiments. The aforementioned storage medium includes various media that can store program code, such as ROM, RAM, magnetic disk, or optical disk.

[0061] 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 foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A highly efficient integral transform method for simulating cohesive fracture phase fields, characterized in that, Applied to finite element simulation software, including: A geometric model of the material to be simulated is established, and the geometric model is meshed to obtain a mesh model; The constitutive relationship between strain and stress is determined based on the energy degradation function, the elastic stiffness tensor of the material to be simulated, and the inherent mechanical parameters; and the governing equations are determined based on the dissipation function functional, the storage function, and the laws of thermodynamics; and the residual equations are obtained according to the governing equations and the constitutive relationship; wherein the inherent mechanical parameters include fracture energy, intrinsic characteristic length, Young's modulus, tensile strength, shear modulus, and shear strength. The residual equation is linearized to transform the governing equation into a matrix governing equation; the matrix governing equation is used to describe the relationship between the variation of the displacement degree of freedom, the variation of the phase field degree of freedom, and the residual through the stiffness matrix; Load conditions are applied to the outer boundary of the mesh model, and the phase field distribution is obtained by solving the matrix control equation step by step by loading the load. For each load step, the iteration of the load step is completed when convergence is reached.

2. The method according to claim 1, characterized in that, The constitutive relation is determined by the following formula: in, For stress tensor; It is the energy degradation function; For effective stress tensor; Let be the elastic stiffness tensor of the material to be simulated; For strain tensor; The energy degradation function is determined by the following formula: in, These are intermediate material parameters; Let be the energy density eigenfunction; The fracture energy is mentioned above; The length of the intrinsic feature; The Young's modulus; The tensile strength is mentioned above; The shear modulus; The shear strength is mentioned. d This refers to the phase field value; And / or, The determination of governing equations based on dissipation function functionals, storage functions, and thermodynamic laws includes: Based on the variation of the fracture energy and crack density function, determine the variation of the dissipation function functional; The energy storage function of the solid system of the material to be simulated is converted into a function of the strain tensor and the phase field variables, and the variation of the energy storage density functional is determined based on the stress tensor and the damage energy release rate. Substituting the variation of the dissipation function functional, the variation of the energy storage density functional, and the virtual work of external forces represented by volume forces and surface traction forces into the laws of thermodynamics, we obtain the first equation; using the divergence theorem, the first equation is transformed into the governing equation.

3. The method according to claim 2, characterized in that, The variational factor of the dissipative function functional is determined by the following formula: in, For the variational of the dissipative function functional; The fracture energy is mentioned above; This is the variation of the crack density function; , This is the cracked area; The material to be simulated is a solid system; The variational function of the energy storage density functional is determined by the following formula: in, For the variation of the energy storage density functional; For stress tensor; It is a symmetric gradient operator; For the variation of the displacement field; The damage energy release rate; The variation of the crack phase field; The governing equations are determined by the following formula: in, This is the damage flux vector; Net damage source per unit volume; For stress tensor; , These are the body force and surface traction force determined by the load step, respectively; This is the cracked area; The material to be simulated is a solid system; This is the force boundary of the solid system. ; Let be the normal vector of the solid system; This is the normal vector of the crack region; For Hamiltonian operators.

4. The method according to claim 1, characterized in that, The mesh model includes several multi-field elements and several nodes; Based on the governing equations and the constitutive relation, the residual equations are obtained, including: Based on the mesh type of the multi-field element and the number of nodes, the interpolation function matrix of the displacement field, the interpolation function matrix of the phase field, the strain-displacement compatibility matrix, and the phase field gradient compatibility matrix are determined respectively. Substituting the interpolation function matrix of the displacement field, the interpolation function matrix of the phase field, the strain-displacement compatibility matrix, and the phase field gradient compatibility matrix into the governing equation, the residual equation is obtained; the residual equation is determined by the following formula: in, For the residuals of the displacement field; The residual of the phase field; and These are the interpolation function matrices for the displacement field and the phase field, respectively. , These are the body force and surface traction force determined by the load step, respectively; , These are the strain-displacement compatibility matrix and the phase-field gradient compatibility matrix, respectively. For stress tensor; The finite element computational domain of the mesh model is... The boundary of the finite element computational domain, The finite element computational domain is defined for the crack region. Let be the derivative of the energy degradation function with respect to space; This represents the historical maximum effective damage energy release rate. The fracture energy is mentioned above; The length of the intrinsic feature; The derivative of the phase field profile function with respect to space; d This refers to the phase field value; For Hamiltonian operators.

5. The method according to claim 4, characterized in that, Each load step also includes: The equivalent effective stress is calculated using the stress tensor from the previous iteration. The effective energy damage release rate of the previous iteration step is calculated based on the equivalent effective stress and the Young's modulus. The critical effective energy damage release rate is calculated based on the tensile strength and the Young's modulus. The maximum value between the effective energy damage release rate of the previous iteration step and the critical effective energy damage release rate is selected as the historical maximum value of the effective damage energy release rate of the current iteration step.

6. The method according to any one of claims 1 to 5, characterized in that, The phase field distribution result obtained by solving the matrix control equation step by step by applying load steps includes: The matrix control equations are transformed using the overall stiffness matrix to obtain the overall control equations; the overall control equations are then modified to obtain the modified equations. For each of the aforementioned load steps, the following is performed: S0, the overall residual is calculated based on the initial phase field value and the initial stress tensor; wherein, the overall residual includes the residual of the displacement field and the residual of the phase field, and is calculated by the residual equation; S1, Substitute the overall variable increment, overall residual change and overall stiffness matrix of the previous iteration into the correction equation to obtain the overall stiffness matrix of the current iteration; S2, Substitute the overall residual, the overall residual of the previous iteration, and the overall stiffness matrix of the current iteration into the overall control equation to calculate the overall variable increment of the current iteration; wherein, the overall variable increment includes stress increment and phase field increment; S3, based on the stress increment and the initial stress tensor, obtain the first stress tensor; and based on the phase field increment and the initial phase field value, determine the first phase field value; S4, Substitute the first stress tensor and the first phase field value into the residual equation to obtain the current overall residual; S5. Determine whether convergence has been achieved based on the current overall residual; if yes, end the iteration of this load step; otherwise, proceed to step S6. S6, increment the iteration count by 1, and based on the first phase field value and the first stress tensor, update the initial phase field value and the initial stress tensor for the current iteration, and return to step S0.

7. The method according to claim 6, characterized in that, The matrix control equation is determined by the following formula: in, This is the first component; This is the second component; For the variation of displacement degrees of freedom; For the variation of the phase field degrees of freedom; For the residuals of the displacement field; The residual of the phase field; The overall control equation is determined by the following formula: in, Let $\mathbf{k}$ be the global stiffness matrix for the $k$-th iteration. Let this be the global variable for the k-th iteration; For the (k-1)th iteration, this is the global variable; The global residual for the k-th iteration; Let be the global residual of the (k-1)th iteration.

8. A highly efficient integral transform cohesive fracture phase field simulation device, characterized in that, include: The first construction module is used to establish a geometric model of the material to be simulated, and to mesh the geometric model to obtain a mesh model; The second construction module is used to determine the constitutive relationship between strain and stress based on the energy degradation function, the elastic stiffness tensor of the material to be simulated, and the inherent mechanical parameters; and to determine the governing equations based on the dissipation function functional, the storage function, and the thermodynamic laws; and to obtain the residual equations according to the governing equations and the constitutive relationship; wherein the inherent mechanical parameters include fracture energy, intrinsic characteristic length, Young's modulus, tensile strength, shear modulus, and shear strength; The transformation module is used to linearize the residual equation and convert the control equation into a matrix control equation; the matrix control equation is used to describe the relationship between the variation of the displacement degree of freedom, the variation of the phase field degree of freedom, and the residual through the stiffness matrix; The simulation solution module is used to apply load conditions on the outer boundary of the mesh model and obtain the phase field distribution results by progressively loading load steps based on the matrix control equations; wherein, for each load step, when convergence is reached, it is determined that the iteration of the load step is completed.

9. A computer device, characterized in that, The computer device includes a memory and a processor. The memory is used to store computer programs, and the processor is used to execute the computer programs stored in the memory to implement the steps of the method according to any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The storage medium stores a computer program, which, when executed by a processor, implements the steps of the method described in any one of claims 1-7.

Images