An adaptive phase-field method for three-dimensional brittle material fracture simulation

Abstract

The present invention discloses an adaptive phase field method for simulating the fracture of three-dimensional brittle materials, comprising the following steps: obtaining material property parameters by conducting physical experiments, establishing a phase field control equation applicable to the fracture of three-dimensional brittle materials, establishing a trilinear adaptive unit and adaptive criteria applicable to the phase field model, establishing a finite element model with the same size and shape as the simulation object, assigning the measured material parameters and the boundary conditions to be applied to the finite element model, selecting a softening curve suitable for the simulated material, solving the control equation using a high-precision iterative algorithm to obtain the stress field and phase field distribution, judging the degree of damage of the specimen based on the phase field value, thereby ultimately achieving accurate simulation of the crack propagation process, crack propagation direction and crack propagation path of the composite material under load. The present invention ensures that the shape function form always maintains the trilinear characteristic in the unit containing any number of nodes, solving the computational efficiency problem of the traditional phase field model in solving the crack process.

Description

Technical Field

[0001] The present invention belongs to the field of computational mechanics, and in particular relates to an adaptive phase field method for three-dimensional brittle material fracture simulation. Background Art

[0002] As an emerging dispersion method based on finite elements, the phase field model has many advantages in simulating the fracture failure of brittle materials. In recent years, it has been widely used in the analysis and research of various solid materials such as rock, concrete, ceramics, and sea ice. One of the main tasks of the phase field model is to predict the crack propagation path of material failure. In order to achieve simulation accuracy, the mesh needs to be finely divided at the crack propagation location. However, fine meshing of the entire computational domain will greatly increase the amount of calculation and time consumption. Manual area division makes it difficult to find reasonable areas for mesh refinement without knowing the crack propagation path. The quality of the mesh division directly determines the accuracy of the simulation. How to perform reasonable mesh division has always been a difficult problem to solve for the phase field model and even other methods of simulating crack propagation paths.

[0003] The present invention mainly solves the problem of limitations of meshing and calculation efficiency of phase field models when simulating crack propagation of brittle materials.

[0004] Phase field models have been dedicated to accurately predicting and simulating the fracture process of brittle materials. Reasonable meshing is one of the key factors to ensure computational accuracy. When the crack propagation path is unknown, manually refining the mesh for the entire crack propagation area is almost impossible, and even more so for complex cracking paths. Adaptive meshing is the perfect solution to the mesh refinement problem. In the course of its development, there are roughly several categories of adaptive meshing: h-method, p-method, r-method, and methods formed by combining the above methods. However, these traditional methods have more or less the following problems: (1) It is difficult to ensure computational accuracy during the adaptive process; (2) The algorithm is complex and difficult to implement; (3) Every time a node is inserted into a unit during the adaptive process, the unit order increases once, and the computational complexity increases exponentially, making the solution difficult, and only a few nodes can be inserted. Summary of the Invention

[0005] The purpose of the present invention is to provide an adaptive method for simulating the three-dimensional fracture process of brittle materials such as rock and concrete using a phase field model.

[0006] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0007] An adaptive phase field method for three-dimensional brittle material fracture simulation, characterized by comprising the following steps:

[0008] S1. Physical testing using a universal material testing machine to obtain material property parameters, wherein the material property parameters include (1) elastic modulus, (2) Poisson's ratio, and (3) density, and energy release rate and uniaxial tensile strength of the material are obtained by a load-fracture test;

[0009] S2. Establish the phase field governing equation applicable to the fracture of three-dimensional brittle materials. The governing equation is as follows:

[0010]

[0011] Among them: the first equation is called the elastic equation, which is used to solve the displacement field and stress field of the simulated object. The second equation is called the crack evolution equation, which is used to solve the damage field of the simulated object. The third and fourth equations are Newman-type boundary conditions. Ω and Θ represent the area where stress changes and damage occur, respectively. u and s are the displacement field and phase field, respectively. The displacement field has three components and is used to express motion in three-dimensional space. The phase field is a quantity that changes in the interval [0,1]. When s=0, it means that the material is intact. When s=1, it means that the material is completely broken and damaged. is the partial derivative of the phase field with respect to the spatial coordinates, l0 is the crack diffusion width, σ(u,s) is the degenerate stress field, div is the divergence sign, ω′(s) is the derivative of the degenerate function ω(s) with respect to the phase field s, ψ0 is the non-degenerate elastic potential energy, G c is the measured energy release rate, is a scale parameter, is the crack geometry function, is the force boundary condition, n and n Θ are the external boundary normal vectors of Ω and Θ respectively;

[0012] S3. Establish a trilinear adaptive element and adaptive criteria suitable for the phase field model, so that the shape function of the element always maintains a trilinear form when any number of nodes are inserted, thus reducing the computational effort;

[0013] S4. Establish a finite element model with the same size and shape as the simulation object, and assign the measured material parameters and the boundary conditions to be applied to the finite element model;

[0014] S5. Select a softening curve suitable for the simulated material to obtain a more realistic fracture process;

[0015] S6. Use high-precision algorithms to solve the governing equations to obtain stress and phase field distributions, and select appropriate solution strategies based on the problem being analyzed;

[0016] S7. In a load step, if the size of an element meets the adaptive criteria, nodes are inserted into the element and adaptive refinement is performed. The element is refined into several smaller sub-elements. After the element refinement, the load step is recalculated until no element in the load step needs to be refined, and the next load step is started.

[0017] S8. Determine the damage level of the specimen based on the phase field value. If the phase field value at a certain position reaches 1, it is determined that fracture has occurred at this displacement. The fracture state within the simulated area during the loading process is then analyzed.

[0018] S9. The calculation ends when all applied boundary conditions have been loaded or the material breaks completely and loses the ability to bear the load.

[0019] Furthermore, the degraded stress field σ(u , s) is obtained by adopting the mixed mode to degenerate, and then degenerate the entire elastic potential energy of the damage zone to simplify the calculation. The degenerate form of the stress field is:

[0020]

[0021] in: is the non-degenerate stress tensor, is the elastic matrix, E0 and ν0 are the elastic modulus and Poisson’s ratio respectively, is the fourth-order unit tensor, 1 is the second-order unit tensor, ε is the strain tensor, is the non-degenerate elastic potential energy, ε kk , ε ii , ε ij , ε ij is the component of ε, λ and μ are Lame constants, and the form of the degradation function is:

[0022]

[0023] R(s)=b1s+b1b2s 2 +b1b2b3s 3 +…=b1s·Q(s) (4)

[0024] Q(s)=1+b2s+b2b3s 2 +… (5)

[0025] Where: χ≥2 is an exponential, R(s) and Q(s) are polynomials, It is an intermediate function used to simplify the expression of the degradation function. b1, b2, b3... are the coefficients of the polynomial. Generally, Q(s) can be taken as a quadratic function to reduce the complexity of the calculation.

[0026] Furthermore, establishing a trilinear adaptive unit suitable for the phase field model specifically includes:

[0027] For a three-dimensional hexahedral element, assuming that there are n nodes inside the element, the shape function needs to meet the following conditions:

[0028] Unit decomposition:

[0029] Linearly complete:

[0030] Kronecker-delta properties:

[0031] The value is 0 on the boundary that does not contain this node: φ j (ξ pe )=0 (9)

[0032] Among them: j, k are node numbers, φ j (ξ) is the shape function corresponding to the jth node, ξ is the local coordinate of any point in the unit, j is the local coordinate of the jth node, ξ pe is the boundary that is not adjacent to node j. The shape function of the element is actually obtained by calculating the basis function group. For each additional node inside the element, an element is added to the basis function group. For an eight-node hexahedral element, the basis function group is in the form of:

[0033] P T =[1 ξ η ζ ξη ηζ ζξ ξηζ] (10)

[0034] Where: P is the basis function group, P T Its transpose matrix, (ξ,η,ζ) are the three spatial components of the coordinate ξ, and each time a node is inserted into the unit, it is necessary to add T When the inserted node κ is on the edge parallel to the ξ direction, the additional basis function is in the form of |ξ-ξ κ |(η+sign(η κ ))(ζ+sign(ζ κ )); When the inserted node κ is on the edge parallel to the η direction, the additional basis function is added in the form of |η-η κ |(ξ+sign(ξ κ ))(ζ+sign(ζ κ )); When the inserted node κ is on the edge parallel to the ζ direction, the additional basis function is added in the form of |ζ-ζ κ |(ξ+sign(ξ κ ))(η+sign(η κ)); When the inserted node κ is on a plane perpendicular to the ξ direction, the additional basis function is in the form of |η-η κ ||ζ-ζ κ |(ξ+sign(ξ κ )); When the inserted node κ is on a plane perpendicular to the η direction, the additional basis function is in the form of |ξ-ξ κ ||ζ-ζ κ |(η+sign(η κ )); When the inserted node κ is on a plane perpendicular to the ζ direction, the additional basis function is added in the form of |ξ-ξ κ ||η-η κ |(ζ+sign(ζ κ )), these basis functions are all trilinear basis, so for a unit with n nodes, the shape function of the jth node can be expressed as

[0035]

[0036] Where: sign(x) is the sign function, when x is greater than 0, sign(x) = 1, when x is equal to 0, sign(x) = 0, when x is less than 0, sign(x) = -1, P j (ξ) are the components of the basis function group P, β j For P j The coefficient corresponding to (ξ) changes with the change of the node distribution inside the unit and can be calculated according to the conditions of formula (6)-(9). Therefore, for units with arbitrary nodes, their shape functions remain in trilinear form.

[0037] Furthermore, the established grid adaptation criteria are as follows:

[0038] When the unit size meets

[0039]

[0040] When , the unit needs to be refined into smaller units, where d is the size of the unit at the end of the previous load step calculation, and its value is taken as the length of the longest side of the unit. The element size required to meet the calculation accuracy of this load step can be obtained by the following formula:

[0041]

[0042] Where: d0 is the initial size of the unit generated by the finite element mesh division before calculation, s e is the unit phase field value calculated in this load step, s minIt is the minimum phase field threshold for starting mesh refinement, which is generally set to 0.01-0.1. c is the ratio of l0 set to meet the calculation accuracy to the minimum mesh size at the complete fracture of the model.

[0043] Furthermore, a softening curve suitable for the simulated material is selected to obtain a more realistic fracture process. The softening curve can be applied by changing the coefficients in formula (5). The values of each coefficient can be calculated using the following formula:

[0044]

[0045]

[0046]

[0047] Where: E0 is the elastic modulus, f t is the tensile strength of the material, k0<0 is the initial slope of the selected softening curve, w c is the crack opening when the stress at the crack tip disappears, k0 and w c Through physical test measurements, for different softening curves, the value can be

[0048] Linear softening curve:

[0049] Exponential softening curve:

[0050] Cornelissen softening curve:

[0051] Among them: the linear softening curve and the exponential softening curve are more suitable for brittle materials, and the Cornelissen softening curve is more suitable for quasi-brittle materials such as concrete.

[0052] Furthermore, when a high-precision iterative algorithm is used to solve the control equation to obtain the stress field and phase field distribution, a coupled iterative solution algorithm, a staggered solution algorithm, or a partial iterative solution algorithm is used.

[0053] Compared with the existing technology, the application of trilinear multi-node elements in this method avoids the disadvantage of traditional multi-node elements that cannot insert a large number of nodes due to the limitation of shape function calculation. It allows any number of nodes to be inserted into a unit without significantly increasing the calculation amount, thereby expanding the application of the unit in non-matching mesh problems, adaptive problems, and multi-scale problems. At the same time, the application of adaptive methods avoids the problem of manual mesh refinement around cracks in fracture solution problems. In most cases, the fracture location cannot be known in advance, making effective meshing difficult. This method can automatically perform meshing based on the calculation results of each load step during the calculation process, which not only ensures the accuracy of the calculation but also eliminates the meshing problem. At the same time, this method also minimizes the number of finite element meshes, greatly improving the calculation efficiency. It has been verified that this method can reduce the calculation time by more than 90% compared with non-adaptive methods. The present invention has a simple implementation process, strong practicality, a wide range of applications, and can even be extended to other fields. It can be built in commercial finite element software such as ABAQUS, COMSOL, and ANSYS, and can also be implemented in other open source finite element software or self-programming. BRIEF DESCRIPTION OF THE DRAWINGS

[0054] Attachment Figure 1 It is a schematic diagram of the technical route for implementing the three-dimensional adaptive phase field method for simulating material failure;

[0055] Attachment Figure 2 Schematic diagram of the shape and boundary conditions of the organic glass single-notch beam;

[0056] Attachment Figure 3 FIG. 5 is a diagram of the INSTRON 5969 universal materials testing machine used to measure the material parameters of this example;

[0057] Attachment Figure 4 is the initial finite element mesh division diagram of the sample;

[0058] Attachment Figure 5 is the final crack path and unit distribution diagram after the calculation;

[0059] Attachment Figure 6 is the three-dimensional crack shape state diagram under different loading stages;

[0060] Attachment Figure 7 It is a schematic diagram of the unit node insertion and refinement process. DETAILED DESCRIPTION

[0061] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0062] The preferred embodiments of the present invention will be described in detail below in conjunction with the accompanying drawings, wherein the accompanying drawings constitute a part of this application and are used together with the embodiments of the present invention to illustrate the principles of the present invention, and are not used to limit the scope of the present invention.

[0063] <Example>

[0064] In order to have a clearer understanding of the technical features, purposes and effects of the present invention, specific embodiments of the present invention are now described in detail with reference to the accompanying drawings.

[0065] See Figure 1 The present invention provides an adaptive phase field method for three-dimensional brittle material fracture simulation. Taking the three-point bending failure process simulation of a single-notched organic glass beam specimen as an example, the specimen is a cubic block with a size of 260mm×60mm×10mm. An initial crack with a 45-degree inclination is located at the bottom center of the beam. A vertical downward displacement load is applied to the top center of the beam. Figure 2 As shown, the total displacement applied is u* = 0.75 mm and is divided into one hundred sub-load steps.

[0066] 1. Physical tests are performed on a universal material testing machine to obtain material property parameters, wherein the material property parameters include (1) elastic modulus, (2) Poisson's ratio, and (3) density. Loading fracture tests are performed to obtain the energy release rate and tensile strength of the material, such as Figure 3 As shown;

[0067] 2. Establish the phase field governing equation applicable to the fracture of three-dimensional brittle materials. The governing equation is as follows:

[0068]

[0069] Among them: the first equation is called the elastic equation, which is used to solve the displacement field and stress field of the simulated object. The second equation is called the crack evolution equation, which is used to solve the damage field of the simulated object. The third and fourth equations are Newman-type boundary conditions. Ω and Θ represent the area where stress changes and damage occur, respectively. u and s are the displacement field and phase field, respectively. The displacement field has three components and is used to express motion in three-dimensional space. The phase field is a quantity that changes in the interval [0,1]. When s=0, it means that the material is intact. When s=1, it means that the material is completely broken and damaged. is the partial derivative of the phase field with respect to the spatial coordinates, l0 is the crack diffusion width, σ(u,s) is the degenerate stress field, div is the divergence sign, ω′(s) is the derivative of the degenerate function ω(s) with respect to the phase field s, ψ0 is the non-degenerate elastic potential energy, G c is the measured energy release rate, is a scale parameter, is the crack geometry function, is the force boundary condition, n and n Θ are the external boundary normal vectors of Ω and Θ respectively.

[0070] 3. Establish a trilinear adaptive element and adaptive criteria suitable for phase field models, so that the shape function of the element always maintains a trilinear form when any number of nodes are inserted, reducing the amount of calculation;

[0071] 4. Establish a finite element model with the same size and shape as the simulation object and perform mesh division, such as Figure 4 As shown in FIG, the finite element model is given the measured material parameters and the boundary conditions to be applied. The material parameters of each component of the composite material required by this method are: (1) elastic modulus, (2) Poisson's ratio, (3) uniaxial tensile strength, and (4) energy release rate. The uniaxial tensile strength and energy release rate are obtained through physical tests on a material testing machine. The final material parameters are determined as follows: elastic modulus E0 = 28 GPa, Poisson's ratio ν0 = 0.38, tensile strength f t =0.02GPa, energy release rate G c =5×10 -4 kN / mm.

[0072] 5. Select a softening curve suitable for the simulated material to obtain a more realistic fracture process;

[0073] 6. Use a high-precision iterative algorithm to solve the governing equations to obtain the stress field and phase field distribution, and select an appropriate iterative solution strategy based on the problem being analyzed. This example adopts the strategy of solving the elasticity equation and the crack evolution equation simultaneously and coupling them and using the Newton iteration method;

[0074] 7. In a load step, if the size of a unit meets the adaptive criteria, nodes are inserted into the unit and adaptive refinement is performed to refine the unit into several smaller sub-units. After the unit refinement, the load step is recalculated until there are no units in the load step that need to be refined and the next load step is started.

[0075] 8. Determine the degree of damage to the specimen based on the phase field value. If the phase field value at a certain position reaches 1, it is determined that fracture damage has occurred at this displacement, and then analyze the fracture state in the simulation area during the loading process.

[0076] 9. The calculation ends when all the applied boundary conditions have been loaded or the material breaks and loses the ability to bear the load. The grid distribution and crack path are finally obtained. Figure 5 As shown in Figure 2, the three-dimensional crack shapes at different loading stages are as follows: Figure 6 shown.

[0077] For step 2, a phase field governing equation suitable for three-dimensional brittle material fracture is established. Unlike other phase field models that use energy decomposition and partial energy degradation, the stress field σ(u,s) in this method is obtained by hybrid mode degradation, and then the entire elastic potential energy of the damage zone is degraded to simplify the calculation. The degraded form of the stress field is:

[0078]

[0079] in: is the non-degenerate stress tensor, is the elastic matrix, E0 and ν0 are the elastic modulus and Poisson’s ratio respectively, is the fourth-order unit tensor, 1 is the second-order unit tensor, ε is the strain tensor, is the non-degenerate elastic potential energy, λ and μ are Lame constants, and the degradation function is in the form of:

[0080]

[0081] R(s)=b1s+b1b2s 2 +b1b2b3s 3 +…=b1s·Q(s) (4)

[0082] Q(s)=1+b2s+b2b3s 2 +… (5)

[0083] Where: χ≥2 is an exponential, R(s) and Q(s) are polynomials, It is an intermediate function used to simplify the expression of the degradation function. b1, b2, b3... are the coefficients of the polynomial. Generally, Q(s) can be taken as a quadratic function to reduce the complexity of the calculation.

[0084] For step 3, a trilinear adaptive element suitable for the phase field model is established so that the shape function of the element always maintains a trilinear form, thereby avoiding the problem that the polynomial degree of the shape function of the traditional multi-node element increases with the number of nodes. For a three-dimensional hexahedral element, assuming that there are n nodes inside the element, the shape function needs to meet the following conditions:

[0085] Unit decomposition:

[0086] Linearly complete:

[0087] Kronecker-delta properties:

[0088] The value is 0 on the boundary that does not contain this node: φ j (ξ pe )=0 (9)

[0089] Among them: j, k are node numbers, φ j (ξ) is the shape function corresponding to the jth node, ξ is the local coordinate of any point in the unit, j is the local coordinate of the jth node, ξ pe is the boundary that is not adjacent to node j. The shape function of the element is actually obtained by calculating the basis function group. For each additional node inside the element, an element is added to the basis function group. For an eight-node hexahedral element, the basis function group is in the form of:

[0090] P T =[1 ξ η ζ ξη ηζ ζξ ξηζ] (10)

[0091] Where: P is the basis function group, P T Its transpose matrix, (ξ,η,ζ) are the three spatial components of the coordinate ξ, and each time a node is inserted into the unit, it is necessary to add T When the inserted node κ is on the edge parallel to the ξ direction, the additional basis function is in the form of |ξ-ξ κ |(η+sign(η κ ))(ζ+sign(ζ κ )); When the inserted node κ is on the edge parallel to the η direction, the additional basis function is added in the form of |η-η κ |(ξ+sign(ξ κ ))(ζ+sign(ζ κ )); When the inserted node κ is on the edge parallel to the ζ direction, the additional basis function is added in the form of |ζ-ζ κ |(ξ+sign(ξ κ ))(η+sign(η κ )); When the inserted node κ is on a plane perpendicular to the ξ direction, the additional basis function is in the form of |η-η κ ||ζ-ζ κ |(ξ+sign(ξ κ )); When the inserted node κ is on a plane perpendicular to the η direction, the additional basis function is in the form of |ξ-ξ κ ||ζ-ζ κ |(η+sign(η κ)); When the inserted node κ is on a plane perpendicular to the ζ direction, the additional basis function is added in the form of |ξ-ξ κ ||η-η κ |(ζ+sign(ζ κ )), these basis functions are all trilinear basis, so for a unit with n nodes, the shape function of the jth node can be expressed as

[0092]

[0093] Where: sign(x) is the sign function, when x is greater than 0, sign(x) = 1, when x is equal to 0, sign(x) = 0, when x is less than 0, sign(x) = -1, P j (ξ) are the components of the basis function group P, β j For P j The coefficient corresponding to (ξ) changes with the change of the node distribution inside the unit and can be calculated according to the conditions of formula (6)-(9). Therefore, for units with arbitrary nodes, their shape functions remain in trilinear form.

[0094] For step 3, it is necessary to establish an efficient grid adaptation criterion. When the grid size meets the criterion, the original unit is refined into smaller units. In the criterion established by this method, when the unit size meets

[0095]

[0096] When , the unit needs to be refined into smaller units, where d is the size of the unit at the end of the previous load step calculation, and its value is taken as the length of the longest side of the unit. The element size required to meet the calculation accuracy of this load step can be obtained by the following formula:

[0097]

[0098] Where: d0 is the initial size of the unit generated by the finite element mesh division before calculation, s e is the unit phase field value calculated in this load step, s min The minimum phase field threshold for mesh refinement is set, which is generally set to 0.01-0.1. c is the ratio of l0 to the minimum mesh size at the point where the model is completely broken, which is set to meet the calculation accuracy. It is recommended to set it to c≥4.0. During the refinement process, the parent unit is generally divided into 8 sub-units of 2×2×2 to maximize the calculation accuracy. The unit refinement process is as follows: Figure 7 As shown, in this example, l0=1mm, c=4.0.

[0099] For step 5, a softening curve suitable for the simulated material is selected to obtain a more realistic fracture process. The softening curve can be applied by changing the coefficients in formula (5). The values of each coefficient can be calculated using the following formula:

[0100]

[0101]

[0102]

[0103] Where: E0 is the elastic modulus, f t is the tensile strength of the material, k0<0 is the initial slope of the selected softening curve, w c is the crack opening when the stress at the crack tip disappears, k0 and w c It can be measured through physical tests. For different softening curves, the value can be

[0104] Linear softening curve:

[0105] Exponential softening curve:

[0106] Cornelissen softening curve:

[0107] Among them, the linear softening curve and exponential softening curve are more suitable for brittle materials, and the Cornelissen softening curve is more suitable for quasi-brittle materials such as concrete. In this example, the linear softening curve is used.

[0108] For step seven, a high-precision iterative algorithm is used to solve the control equation to obtain the stress field and phase field distribution, and an appropriate iterative solution strategy is selected according to the problem being analyzed. On the one hand, the solution of the displacement field and the phase field can be regarded as a multi-field coupling problem, and the elastic equation and the crack evolution equation in formula (1) are solved simultaneously, and the Newton iteration method is used to iteratively solve them, and the displacement field and the phase field are obtained at the same time. This is called a coupled iterative solution algorithm. On the other hand, the elastic equation and the crack evolution equation can also be solved alternately. The displacement field is obtained by solving the elastic equation, and the obtained displacement field is substituted into the crack evolution equation to obtain the phase field, and then the phase field is substituted into the elastic equation in the next step without iteration. This is called an interleaved solution algorithm. In order to balance the computational efficiency and computational accuracy, only the elastic equation or only the crack evolution equation can be iterated in the interleaved solution algorithm, which is called a partial iterative solution algorithm. Among them, the coupled iterative solution algorithm has the highest accuracy and the interleaved solution algorithm has the highest efficiency.

[0109] Compared with the existing technology, the application of trilinear multi-node elements in this method avoids the disadvantage of traditional multi-node elements that cannot insert a large number of nodes due to the limitation of shape function calculation. It allows any number of nodes to be inserted into a unit without significantly increasing the calculation amount, thereby expanding the application of the unit in non-matching mesh problems, adaptive problems, and multi-scale problems. At the same time, the application of adaptive methods avoids the problem of manual mesh refinement around cracks in fracture solution problems. In most cases, the fracture location cannot be known in advance, making effective meshing difficult. This method can automatically perform meshing based on the calculation results of each load step during the calculation process, which not only ensures the accuracy of the calculation but also eliminates the meshing problem. At the same time, this method also minimizes the number of finite element meshes, greatly improving the calculation efficiency. It has been verified that this method can reduce the calculation time by more than 90% compared with non-adaptive methods. The present invention has a simple implementation process, strong practicality, a wide range of applications, and can even be extended to other fields. It can be built in commercial finite element software such as ABAQUS, COMSOL, and ANSYS, and can also be implemented in other open source finite element software or self-programming.

[0110] The above embodiments are merely illustrative of the technical solutions of the present invention. The method and apparatus for numerically simulating concealed pipe drainage and salt removal under varying inlet resistance conditions, as described herein, are not limited solely to those described in the above embodiments but are subject to the scope defined by the claims. Any modifications, supplements, or equivalent substitutions made by persons skilled in the art based on these embodiments are within the scope of protection claimed by the claims.

Claims

1. An adaptive phase field method for three-dimensional brittle material fracture simulation, characterized in that: The following steps are involved: S1. Physical testing using a universal material testing machine to obtain material property parameters, wherein the material property parameters include (1) elastic modulus, (2) Poisson's ratio, and (3) density, and energy release rate and uniaxial tensile strength of the material are obtained by a load-fracture test; S2. Establish the phase field governing equation applicable to the fracture of three-dimensional brittle materials. The governing equation is as follows: Among them: the first equation is called the elastic equation, which is used to solve the displacement field and stress field of the simulated object. The second equation is called the crack evolution equation, which is used to solve the damage field of the simulated object. The third and fourth equations are Newman-type boundary conditions. Ω and Θ represent the area where stress changes and damage occur, respectively. u and s are the displacement field and phase field, respectively. The displacement field has three components and is used to express motion in three-dimensional space. The phase field is a quantity that changes in the interval [0,1]. When s=0, it means that the material is intact. When s=1, it means that the material is completely broken and damaged. is the partial derivative of the phase field with respect to the spatial coordinates, l0 is the crack diffusion width, σ(u,s) is the degenerate stress field, div is the divergence sign, ω′(s) is the derivative of the degenerate function ω(s) with respect to the phase field s, ψ0 is the non-degenerate elastic potential energy, G c is the measured energy release rate, is a scale parameter, is the crack geometry function, is the force boundary condition, n and n Θ are the external boundary normal vectors of Ω and Θ respectively; S3. Establish a trilinear adaptive element and adaptive criteria suitable for the phase field model, so that the shape function of the element always maintains a trilinear form when any number of nodes are inserted, thus reducing the computational effort; S4. Establish a finite element model with the same size and shape as the simulation object, and assign the measured material parameters and the boundary conditions to be applied to the finite element model; S5. Select a softening curve suitable for the simulated material to obtain a more realistic fracture process; S6. Use high-precision algorithms to solve the governing equations to obtain stress and phase field distributions, and select appropriate solution strategies based on the problem being analyzed; S7. In a load step, if the size of an element meets the adaptive criteria, nodes are inserted into the element and adaptive refinement is performed. The element is refined into several smaller sub-elements. After the element refinement, the load step is recalculated until no element in the load step needs to be refined, and the next load step is started. S8. Determine the damage level of the specimen based on the phase field value. If the phase field value at a certain position reaches 1, it is determined that fracture has occurred at this displacement. The fracture state within the simulated area during the loading process is then analyzed. S9. The calculation ends when all the applied boundary conditions have been loaded or the material breaks completely and loses the ability to bear the load.

2. The adaptive phase field method for three-dimensional brittle material fracture simulation according to claim 1 is characterized in that: The degraded stress field σ(u,s) in step S2 is obtained by degrading in a mixed mode, and then degrading the entire elastic potential energy of the damaged area to simplify the calculation. The degraded form of the stress field is: in: is the non-degenerate stress tensor, is the elastic matrix, E0 and ν0 are the elastic modulus and Poisson’s ratio respectively, is the fourth-order unit tensor, 1 is the second-order unit tensor, ε is the strain tensor, is the non-degenerate elastic potential energy, ε kk , ε ii , ε ij is the component of ε, λ and μ are Lame constants, and the form of the degradation function is: R(s)=b1s+b1b2s 2 +b1b2b3s 3 +…=b1s·Q(s) (4) Q(s)=1+b2s+b2b3s 2 +…(5) Where: χ≥2 is an exponential, R(s) and Q(s) are polynomials, is an intermediate function used to simplify the expression of the degradation function, b1, b2, b3... are the coefficients of the polynomial, and Q(s) is taken as a quadratic function to reduce the complexity of the calculation.

3. The adaptive phase field method for three-dimensional brittle material fracture simulation according to claim 1 is characterized in that: The establishment of a trilinear adaptive element suitable for the phase field model specifically includes: For a three-dimensional hexahedral element, assuming that there are n nodes inside the element, the shape function needs to meet the following conditions: Unit decomposition: Linearly complete: Kronecker-delta properties: The value is 0 on the boundary that does not contain this node: φ j (ξ pe )=0 (9) Among them: j, k are node numbers, φ j (ξ) is the shape function corresponding to the jth node, ξ is the local coordinate of any point in the unit, j is the local coordinate of the jth node, ξ pe is the boundary that is not adjacent to node j. The shape function of the element is actually obtained by calculating the basis function group. For each additional node inside the element, an element is added to the basis function group. For an eight-node hexahedral element, the basis function group is in the form of: P T =[1 ξ η ζ ξη ηζ ζξ ξηζ] (10) Among them: P is the basis function group, P T Its transpose matrix, (ξ,η,ζ) are the three spatial components of the coordinate ξ, and each time a node is inserted into the unit, it is necessary to add T When the inserted node κ is on the edge parallel to the ξ direction, the additional basis function is in the form of |ξ-ξ κ |(η+sign(η κ ))(ζ+sign(ζ κ )); When the inserted node κ is on the edge parallel to the η direction, the additional basis function is added in the form of |η-η κ |(ξ+sign(ξ κ ))(ζ+sign(ζ κ )); When the inserted node κ is on the edge parallel to the ζ direction, the additional basis function is added in the form of |ζ-ζ κ |(ξ+sign(ξ κ ))(η+sign(η κ )); When the inserted node κ is on a plane perpendicular to the ξ direction, the additional basis function is in the form of |η-η κ ||ζ-ζ κ |(ξ+sign(ξ κ )); When the inserted node κ is on a plane perpendicular to the η direction, the additional basis function is in the form of |ξ-ξ κ ||ζ-ζ κ |(η+sign(η κ )); When the inserted node κ is on a plane perpendicular to the ζ direction, the additional basis function is added in the form of |ξ-ξ κ ||η-η κ |(ζ+sign(ζ κ )), these basis functions are all trilinear basis, so for a unit with n nodes, the shape function of the jth node can be expressed as Where: sign(x) is the sign function, when x is greater than 0, sign(x) = 1, when x is equal to 0, sign(x) = 0, when x is less than 0, sign(x) = -1, P j (ξ) are the components of the basis function group P, β j For P j The coefficient corresponding to (ξ) changes with the change of the node distribution inside the unit and can be calculated according to the conditions of formula (6)-(9). Therefore, for units with arbitrary nodes, their shape functions remain in trilinear form.

4. The adaptive phase field method for three-dimensional brittle material fracture simulation according to claim 1 is characterized in that: The established grid adaptation criteria are as follows: When the unit size meets When , the unit needs to be refined into smaller units, where d is the size of the unit at the end of the previous load step calculation, and its value is taken as the length of the longest side of the unit. The element size required to meet the calculation accuracy of this load step can be obtained by the following formula: Where: d0 is the initial size of the unit generated by the finite element mesh division before calculation, s e is the unit phase field value calculated in this load step, s min is the minimum phase field threshold for starting mesh refinement, set to 0.01-0.1, and c is the ratio of l0 set to meet the calculation accuracy to the minimum mesh size at the complete fracture of the model.

5. The adaptive phase field method for three-dimensional brittle material fracture simulation according to claim 1 is characterized in that: Select a softening curve suitable for the simulated material to obtain a more realistic fracture process. The softening curve can be applied by changing the coefficients in formula (5). The values of each coefficient can be calculated using the following formula: Where: E0 is the elastic modulus, f t is the tensile strength of the material, k0<0 is the initial slope of the selected softening curve, w c is the crack opening when the stress at the crack tip disappears, k0 and w c Through physical test measurements, for different softening curves, the value can be Linear softening curve: Exponential softening curve: Cornelissen softening curve: Among them: the linear softening curve and the exponential softening curve are more suitable for brittle materials, and the Cornelissen softening curve is more suitable for quasi-brittle concrete materials.

6. The adaptive phase field method for three-dimensional brittle material fracture simulation according to claim 1 is characterized in that: When a high-precision iterative algorithm is used to solve the control equation to obtain the stress field and phase field distribution, a coupled iterative solution algorithm, a staggered solution algorithm, or a partial iterative solution algorithm is used.

Images