Interlaced iteration-based integral transformation cohesion fracture phase field simulation method and device

Abstract

The invention discloses an integral transformation cohesion fracture phase field simulation method and device based on staggered iteration. 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 intrinsic mechanical parameters comprise fracture energy, intrinsic characteristic length, Young modulus, tensile strength, shear modulus and shear strength; 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 staggered iteration; and for each load step, when the residual error calculated by the residual error equation reaches a convergence condition, determining that the iteration of the load step is completed. The method can be implemented in ABAQUS software, and is 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 a method and apparatus for simulating cohesive fracture phase fields based on staggered iteration integral transform. Background Technology

[0002] The fracture phase field method is a fracture numerical simulation technique based on nonlocal theory. It effectively eliminates the need for prior settings of crack paths and propagation criteria, and realizes spontaneous crack tracking. It is widely used in materials science and engineering to describe phase transformation, microstructure evolution and damage processes of materials.

[0003] In the fracture phase-field model, the governing equations consist of two coupled partial differential equations, which control the evolution of the displacement field and the phase field, respectively. When solving within the finite element framework, the nodal degrees of freedom include both displacement and phase field degrees of freedom, thus requiring the use of multi-field elements. While numerical implementation is not difficult, the characteristics of the problem present numerous challenges in numerically solving the phase-field model. When both displacement and phase fields are unknowns, the energy functional inherent in the coupled system equations is non-convex, resulting in poor robustness and difficulty in guaranteeing convergence of the overall solution algorithm.

[0004] With the further development of engineering and scientific problems, the demand for solving complex cohesion laws is becoming increasingly significant. Theoretically, the integral transform phase-field model can solve this problem, but its complex mathematical expressions limit its application. Therefore, developing efficient numerical solution algorithms remains a crucial task for fracture phase-field simulation. Another obstacle limiting the use of phase-field damage models in practical problems is their heavy reliance on internal code for numerical simulation. Currently, the integral transform phase-field model is only used in the finite element software package Fenics. While its automatic differentiation function avoids the complex processes of stiffness matrix splicing in finite element calculations, it is not suitable for large-scale engineering applications. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a method and apparatus for simulating cohesive fracture phase fields based on staggered iteration integral transform. This method can be implemented in the finite element simulation software ABAQUS, improving the feasibility and accuracy of fracture phase field simulation in complex structures and large-scale engineering applications.

[0006] This invention provides a method and apparatus for simulating cohesive fracture phase fields based on staggered iteration integral transform. The technical solution is as follows: On the one hand, a method for simulating the cohesive fracture phase field based on staggered iteration integral transform 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. Load conditions are applied to the outer boundary of the mesh model, and the phase field distribution is obtained by solving the load step by step based on staggered iteration; wherein, for each load step, the iteration of the load step is completed when the residual calculated by the residual equation reaches the convergence condition.

[0007] On the other hand, an integral transform cohesive fracture phase field simulation device based on staggered iteration 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 simulation solution module is used to apply load conditions on the outer boundary of the mesh model and obtain the phase field distribution result by gradually loading load steps based on staggered iteration; wherein, for each load step, when the residual calculated by the residual equation reaches the convergence condition, 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 above-described method for simulating cohesive fracture phase field based on interleaved iteration integral transform.

[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 above-described method for simulating cohesive fracture phase field based on interleaved iteration integral transform 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 above-described method for simulating cohesive fracture phase field based on interleaved iteration integral transform.

[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 simplicity and efficiency, utilizing the built-in UEL subroutine in ABAQUS to implement an interleaved iterative solution framework for the integral transform phase-field model. This significantly improves the stability and convergence efficiency of the solution process, while drastically reducing the number of iteration steps. Thus, this implementation framework not only ensures high consistency and accuracy of the method under different working conditions but also lays a solid foundation for the widespread application of the integral transform phase-field model in complex engineering scenarios, demonstrating good scalability and engineering applicability, especially 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 a method for simulating cohesive fracture phase field based on interleaved iteration integral transform according to 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 phase field distribution of concrete at three points during bending, provided by an embodiment of the present invention; Figure 4 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 5 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 6 This is a hardware architecture diagram of a computer device provided in an embodiment of the present invention; Figure 7This is a structural diagram of an integral transform cohesive fracture phase field simulation device based on staggered iteration 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 a method for simulating the cohesive fracture phase field based on staggered iteration integral transform, 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: Apply load conditions to the outer boundary of the mesh model, and obtain the phase field distribution results by gradually loading load steps based on staggered iteration; wherein, for each load step, when the residual calculated by the residual equation reaches the convergence condition, 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 simplicity and efficiency, utilizing the built-in UEL subroutine of ABAQUS to implement an interleaved iterative solution framework for the integral transform phase-field model. This significantly improves the stability and convergence efficiency of the solution process, while drastically reducing the number of iteration steps. Thus, this implementation framework not only ensures high consistency and accuracy of the method under different working conditions but also lays a solid foundation for the widespread application of the integral transform phase-field model in complex engineering scenarios, demonstrating good scalability and engineering applicability, especially suitable for the detailed simulation of large-scale crack propagation and complex failure processes.

[0018] The following description Figure 1 The execution method of 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. , is 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. , is 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, ; Let be 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) 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; 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.

[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, the equivalent effective stress is determined by the following formula: (13) in, , These are the first principal stress and the first stress invariant of the effective stress, respectively; For effective deviatoric stress invariants; For Macauly operators; , Uniaxial compressive strength; For tensile strength, the specific formula for calculating the equivalent effective stress is selected based on the inherent properties of the material being simulated.

[0026] 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.

[0027] 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.

[0028] 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.

[0029] 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).

[0030] 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).

[0031] In step 104, the phase field distribution result is obtained by progressively applying load steps based on staggered iteration, including: For each load step, execute: The residual equations are linearized 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. S0, in the k-th iteration and when k is odd, the phase field value and stress tensor updated after the (k-1)-th iteration are used as the initial phase field value and initial stress tensor of the current iteration, and the residual of the displacement field calculated based on the initial phase field value and initial stress tensor is obtained; where the residual equation is used to calculate the residual of the displacement field and the residual of the phase field. S1, based on the phase field value (i.e. the initial phase field value) and constitutive relation, the first component in the stiffness matrix is ​​calculated (see formula (43)). S2, substitute the residual and first component of the displacement field into the matrix control equation (formula (41)) to calculate the stress increment, and take the sum of the stress increment and the initial stress tensor as the first stress tensor; S3, increment the iteration number by 1, and calculate the residual of the phase field based on the phase field value at the k-th iteration and the first stress tensor at the (k+1)-th iteration; wherein, the phase field value at the k-th iteration is the same as the initial phase field value; S4, based on the first stress tensor, energy degradation function and inherent mechanical parameters, the second component in the stiffness matrix is ​​calculated (see formula (43)). S5, substitute the residual and second component of the phase field into the matrix control equation (formula (41)) to calculate the phase field increment, and take the sum of the phase field increment and the initial phase field value as the first phase field value; S6, Substitute the first stress tensor, the body force and surface traction force determined by the load step into the residual equation (Formula (9)) to obtain the residual of the current displacement field; and Substitute the inherent mechanical parameters, the first phase field value and the first stress tensor into the residual equation (Formula (9)) to obtain the residual of the current phase field; S7, determine whether the convergence condition has been met; if yes, end the iteration of this load step; otherwise, execute step S8. S8, increment the iteration count by 1, update the phase field value and stress tensor of the current iteration, and return to step S0.

[0032] 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 104 include the phase field distribution diagram and stress-strain curve (or load-displacement curve).

[0033] In this embodiment of the invention, an alternating iterative method is used to solve the coupled displacement field and phase field, avoiding the non-convexity of the energy functional during global solution and exhibiting better stability. Specifically, in odd-numbered steps, the phase field value is kept constant, and only the displacement and stress tensors are updated; in even-numbered steps, the stress tensor is kept constant, and only the phase field is updated. For example, in the third iteration, the first stress tensor is obtained using the phase field value updated after the second iteration; in the fourth iteration, the first phase field value is obtained using the phase field value of the first stress tensor updated after the third iteration; this process is repeated alternately until the convergence condition is met.

[0034] In a preferred embodiment, the matrix control equations are 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; The residual of the phase field.

[0035] 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, Newton's iteration method is used to linearize equation (9), where the unknowns 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).

[0036] In one specific implementation, the convergence criteria include: Determine whether the product of the current displacement and the total time-averaged flux of the damage reached in this loading step and the first coefficient is not less than the residual of the current displacement field and the residual of the current phase field. Determine whether the product of the maximum increment of the variable in the current loading step 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.

[0037] Specifically, the convergence condition is determined by the following formula: (44) in, It is the maximum value between the residual of the current displacement field and the residual of the current phase field; 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; (45) 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.

[0038] This invention applies, for the first time, the staggered iteration algorithm to the complex mechanical model of the phase field in integral transformation, and verifies the feasibility and accuracy of the method. Although the staggered iteration algorithm is not provided in the general-purpose finite element simulation software ABAQUS, the user-defined element subroutine UEL provided by ABAQUS can implement the algorithm. Therefore, this invention develops a multi-field custom element coupled with the phase field and displacement field. Taking two dimensions as an example, degrees of freedom 1 and 2 of the element are defined as displacement degrees of freedom, 6 as an additional degree of freedom used to control the staggered iteration, and 11 as the phase field degree of freedom. In the UEL subroutine, the stiffness matrix and residual vector of the element (including the residuals of the displacement field and the phase field) are given by formulas (9) and (43). However, it should be noted that the above formulas need to be summed at each integration point in the UEL subroutine. The specific details are not elaborated here. Among them, the extra degrees of freedom control the convergence of the iteration. When it is an odd number of iterations (fixed phase field), the residual of the displacement field is calculated by formula (9), while the residual of the phase field is set to 0. Usually, the displacement field basically converges after one iteration, while the residual of the extra degrees of freedom is set to be very large to ensure the next iteration (even number of iterations). When it is an even number of iterations (fixed displacement field), the residual of the extra degrees of freedom is set to 0. The residual of the displacement field has been calculated to be basically 0 in the previous step. Therefore, the overall solution will only converge when the calculated residual of the phase field is also small enough.

[0039] It should be noted that since the phase field of the integral transform is nonlinear, the solution usually requires a large number of iterations when using the staggered iterative algorithm. Therefore, when using ABAQUS software, it is necessary to modify the default control parameters and use the fixed increment step calculation.

[0040] 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.

[0041] In one specific implementation, a crack propagation test of a notched concrete beam under three-point bending was used as a standard example, and the method of this invention was verified. First, the geometric dimensions and boundary conditions of the notched concrete beam were defined as follows: length L = 450 mm, width 100 mm, thickness 100 mm, notch width 5 mm, and notch depth 50 mm. A displacement load of 1 mm was applied directly above the notch. During the test, as the load increased, a crack initiated and evolved at the initial notch tip due to stress concentration. The material parameters used in the simulation were: fracture energy 0.113 N / mm², tensile strength 2.4 MPa, Young's modulus 20 GPa, and Poisson's ratio 0.2; the mesh size near the notch was 0.5 mm, and the remaining mesh size was 10 mm to reduce the computational scale. The final experimental and simulation results are as follows: Figures 3 to 6 As shown, Figure 3 The phase field distribution diagram after the simulation calculation is shown, in which the crack extends along the direction of the pre-fabricated notch, consistent with the experimental observation. Figure 4 A comparison is given between the model of this invention (i.e., the integral transform phase-field model) and the cohesive phase-field model proposed in the prior art. Both types of models adopt the UEL interleaved iteration method, such as... Figure 4 As can be seen, the calculation results of the two models are almost identical overall, with only minor differences near the peak value, thus demonstrating the correctness of the method proposed in this invention. Furthermore, Figure 5 The computational efficiency of the two models was compared, by Figure 5 It can be seen that the integral transform phase field model of this invention is significantly more efficient than the cohesive phase field model overall, especially near certain critical steps. The cohesive phase field model has more than 200 iterations, while the integral transform phase field model only has a high number of iterations at key points. In summary, the integral transform cohesive fracture phase field simulation method based on staggered iteration proposed in this invention exhibits excellent performance in both accuracy and efficiency of phase field distribution results, demonstrating significant advantages.

[0042] This invention utilizes the user-defined element subroutine (UEL) in ABAQUS software to implement an interleaved iterative solution framework for the integral transform phase-field model. Compared to previous simplified implementations based on UMAT, the proposed method exhibits significant advantages in both numerical stability and computational efficiency. It is particularly suitable for the detailed simulation of large-scale crack propagation and complex failure processes, thus providing a solid numerical foundation and practical value for efficient and reliable failure analysis of engineering structures.

[0043] like Figure 6 , Figure 7 As shown, this embodiment of the invention provides a device for simulating the cohesive fracture phase field based on interleaved iteration integral transform. The device embodiment can be implemented in software, hardware, or a combination of both. From a hardware perspective, as... Figure 7 The diagram shown is a hardware architecture diagram of a computing device for simulating cohesive fracture phase field based on interleaved iteration integral transform, as provided in an embodiment of the present invention. (Except for...) Figure 7 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 7 As shown, a logical device is formed by the CPU of its computing device reading the corresponding computer program from non-volatile memory into memory and running it. This embodiment provides a device for simulating the cohesive fracture phase field based on interleaved iteration integral transform, comprising: The first construction module 700 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 702 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 simulation solution module 704 is used to apply load conditions on the outer boundary of the mesh model and obtain the phase field distribution results by gradually loading load steps based on staggered iteration; wherein, for each load step, the iteration of the load step is completed when the residual calculated by the residual equation reaches the convergence condition.

[0044] In some specific implementations, the first construction module 700 can be used to execute the above step 100, the second construction module 702 can be used to execute the above step 102, and the simulation solution module 704 can be used to execute the above step 104.

[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 a cohesive fracture phase-field simulation device based on interleaved iteration integral transform. In other embodiments of the present invention, a cohesive fracture phase-field simulation device based on interleaved iteration integral transform 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 a method for simulating cohesive fracture phase field based on interleaved iteration integral transform according to any embodiment of this invention.

[0049] This invention also provides a computer-readable storage medium storing a computer program. When executed by a processor, the computer program causes the processor to perform a method for simulating cohesive fracture phase field based on interleaved iteration integral transform 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 any of the methods for simulating cohesive fracture phase fields based on interleaved iteration integral transform as described in 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 method for simulating the cohesive fracture phase field based on staggered iteration integral transform, 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. Load conditions are applied to the outer boundary of the mesh model, and the phase field distribution is obtained by solving the load step by step based on staggered iteration; wherein, for each load step, the iteration of the load step is completed when the residual calculated by the residual equation reaches the convergence condition.

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; These are eigenfunctions of energy density; 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 variation of the dissipative function functional; The fracture energy is mentioned above; 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 problem through staggered iteration by gradually applying load steps includes: For each of the aforementioned load steps, the following is performed: 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 variational values ​​of the displacement degrees of freedom, the variational values ​​of the phase field degrees of freedom, and the residual through the stiffness matrix; S0, in the k-th iteration and when k is odd, the phase field value and stress tensor updated after the (k-1)-th iteration are used as the initial phase field value and initial stress tensor of the current iteration, and the residual of the displacement field calculated based on the initial phase field value and the initial stress tensor is obtained; wherein, the residual equation is used to calculate the residual of the displacement field and the residual of the phase field. S1, the first component of the stiffness matrix is ​​calculated based on the phase field value and the constitutive relation; S2, Substitute the residual of the displacement field and the first component into the matrix control equation to calculate the stress increment, and use the sum of the stress increment and the initial stress tensor as the first stress tensor; S3, increment the iteration number by 1, and at the (k+1)th iteration, calculate the residual of the phase field based on the phase field value at the kth iteration and the first stress tensor; wherein, the phase field value at the kth iteration is the same as the initial phase field value; S4. The second component of the stiffness matrix is ​​calculated based on the first stress tensor, the energy degradation function, and the inherent mechanical parameters. S5, Substitute the residual of the phase field and the second component into the matrix control equation to calculate the phase field increment, and take the sum of the phase field increment and the initial phase field value as the first phase field value; S6, Substitute the first stress tensor, the body force and surface traction force determined by the load step into the residual equation to obtain the residual of the current displacement field; and substitute the inherent mechanical parameters, the first phase field value and the first stress tensor into the residual equation to obtain the residual of the current phase field; S7. Determine whether the convergence condition has been met; if yes, end the iteration of this load step; otherwise, proceed to step S8. S8, increment the iteration count by 1, update the phase field value and stress tensor of 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.

8. A device for simulating the cohesive fracture phase field based on staggered iteration integral transform, 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 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 gradually loading load steps based on staggered iteration; wherein, for each load step, when the residual calculated by the residual equation reaches the convergence condition, 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