Adaptive Multilevel Mesh Phase-Field Method for Elastic-Brittle Fracture Analysis under Thermal Shock

Abstract

The present invention belongs to the field of computational mechanics and provides an adaptive multi-level grid phase-field method for analyzing elastic-brittle fracture under thermal shock, offering a brand-new numerical calculation method for efficient fracture failure research. In the present invention, the phase field is used to characterize the damage degree of materials, and the multi-level hp finite element is adopted as an adaptive strategy to perform dynamic discretization during the crack propagation process under thermal shock. Compared with the previous adaptive phase-field fracture models, there is no hanging node problem, the numerical implementation is simple and the calculation efficiency is high, and it can effectively predict the fracture failure behavior of structures. At the same time, three methods for prefabricated crack treatment technology are also developed to initialize the phase field, and the temperature field, displacement field and phase field are solved by staggered iteration, which can stably and efficiently solve the large-scale fracture failure problem of multi-field coupling. In addition, the present invention can be further extended to the research on defect prediction in the additive manufacturing process.

Description

Technical Field

[0001] The present invention belongs to the field of computational mechanics and relates to an adaptive mesh phase-field method for elastobrittle fracture analysis under thermal shock. Background Art

[0002] Problems of thermal shock failure widely exist in engineering practice and daily life, such as thermal cracks caused by rapid temperature changes in manufacturing processes (additive manufacturing, quenching, welding, etc.), microcracks generated due to different expansion coefficients during the processing of multiphase composite materials, and spontaneous explosion of window glass in high-temperature weather. Therefore, studying the failure mechanism under thermal shock is helpful for the safety assessment of structural devices. If analyzed by experimental means, it is not only difficult to implement and costly, but also very difficult to consider the influence of multiple external factors on the fracture behavior simultaneously, and it is difficult to observe the crack evolution process. The numerical simulation method, as a new high-performance simulation calculation technology, can assist us to explore the failure process under thermal shock more deeply, so as to comprehensively and accurately predict the failure behavior of materials. However, for the multi-field coupled fracture simulation study under thermal shock, it is still a very challenging topic, mainly because it needs to comprehensively and carefully consider the effective integration of stable crack propagation simulation, multi-physics field coupled solution and efficient adaptive method.

[0003] For the numerical simulation of fracture problems, there are currently two main types: discrete crack representation method and smeared crack representation method. Among them, the discrete crack representation method describes cracks by explicitly defining crack set information, and there are mesh processing methods such as element deletion method, cohesive zone model (CZM) and extended finite element method (XFEM). During the simulation process, cracks will expand, and their geometric shapes change with the load magnitude and time, which requires continuous mesh redivision during the solution process, greatly consuming computational resources. At the same time, due to the limitations of its own algorithm, it is impossible to simply and efficiently handle complex fracture behaviors such as crack initiation, bifurcation and intersection. The smeared crack representation method characterizes cracks through field variables, does not require separate modeling or mesh redivision of cracks, improves the computational efficiency to a certain extent and reduces the numerical implementation difficulty, and has significant advantages in dealing with complex fracture problems.

[0004] In the method for characterizing diffuse cracks, Bourdin et al. developed a diffuse fracture model with a regularized approximation, namely the phase field model, based on the fracture variational principle proposed by Francfort and Marigo in "Numerical experiments in revisited brittle fracture". This model describes cracks through a phase field scalar to characterize the damage value of materials and has been widely applied in dealing with crack initiation and propagation problems. Amor et al. in "Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments" decomposed the strain energy into volumetric-deviatoric strain components and assumed that only the deviatoric stress and positive spherical stress would drive crack evolution. In 2010, Miehe et al. reformulated the phase field fracture model within a thermodynamic framework in the paper "A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits". By decomposing the strain energy into tensile and compressive parts through the strain tensor spectrum, it was considered that only tensile deformation would drive crack evolution, thus preventing unrealistic crack propagation during compression. Borden et al. in "A phase-field description of dynamic brittle fracture" proposed a fourth-order phase field fracture equation within the framework of isogeometric analysis to more accurately simulate the fracture process of brittle materials. To more accurately handle the quasi-brittle fracture of materials, Professor Wu Jianying introduced the concept of cohesive force into the phase field fracture model, proposed the phase field regularized cohesive crack model (PF-CZM), and further improved the calculation efficiency by combining with the BFGS quasi-Newton iteration algorithm. Simulating crack evolution using the phase field model does not require updating the geometric topology of the model, which has advantages in dealing with complex fracture modes such as crack initiation, propagation, and intersection, and is widely applied to various types of fracture failure problems, such as thermal cracking, dynamic fracture, hyperelastic fracture, hydraulic fracturing, impact failure, etc. Among them, the research on thermal cracking problems helps with the safety assessment of manufacturing processes such as additive manufacturing, quenching, and casting, and has received the attention of many scholars. Tangella et al. in "Hybrid phase-field modeling of thermo-elastic crack propagation" adopted a hybrid phase field model to simulate the fracture of elasto-brittle materials.Ruan et al. proposed a thermo-mechanical coupled phase-field fracture model in "A thermo-mechanical phase-field fracture model: Application to hot cracking simulations in additive manufacturing", discussed the influence of temperature on fracture modes, and applied it to the simulation of microcrack initiation in additive manufacturing. Mandal et al. used the PF-CZM model to simulate the thermal cracking phenomenon of elasto-brittle materials in "Fracture of thermo-elastic solids: Phase-field modeling and new results with an efficient monolithic solver", and further improved the calculation efficiency with the BFGS iterative method.

[0005] The phase-field model has been widely used in simulating the fracture process. However, in the numerical simulation process, an extremely fine mesh is required, resulting in a large computational amount. To meet the mesh requirements without significantly increasing the computational cost, the adaptive mesh method has been widely applied to the phase-field fracture simulation. Heister et al. proposed a predictor-corrector mesh adaptivity method in "A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach". This method is an h-adaptivity strategy that locally subdivides the mesh in real time during crack propagation. Patil et al. proposed a phase-field fracture model under adaptive multi-scale finite element (MsFEM), which can significantly reduce the computational cost and improve the computational efficiency. Recently, Zander et al. proposed a new adaptive method called the multi-level hp finite element method (The multi-level hp-FEM) in "Multi-level hp-adaptivity: high-order mesh adaptivity without the difficulties of constraining hanging nodes". This method uses the superposition of fine meshes to achieve local subdivision. In most previous adaptive methods, fine meshes are used to replace coarse meshes, which will cause the problem of hanging nodes and lead to complex numerical implementation. However, the multi-level hp finite element method constructs the shape function using the superposition idea and the integral Legendre basis function, so there is no hanging node problem, and the numerical implementation is more convenient. At the same time, it makes full use of the advantages of h-subdivision and p-subdivision, allowing high convergence in dealing with problems with high gradient changes. These characteristics are applicable to the phase-field fracture simulation, locally covering multiple levels of fine meshes at the crack, thereby reducing the consumption of computing resources.

[0006] In the present invention, the multi-level hp finite element method is used as an adaptive strategy to perform dynamic discretization on the phase-field fracture model under thermo-mechanical coupling, and three prefabricated crack treatment techniques are developed to accurately apply the phase-field boundary conditions. In addition, the staggered iteration method is used to solve the temperature field, displacement field, and phase-field in sequence. The research on an efficient numerical method for analyzing the fracture failure problem of elastic-brittle materials under thermal shock proposed in the present invention can efficiently capture the crack evolution process, without the hanging node problem in traditional adaptive methods, and the numerical implementation is simple. It is used to simulate the elastic-brittle fracture mechanics behavior of structures under thermal shock loads, providing an effective way for the accurate analysis of engineering fracture problems. Summary of the Invention

[0007] For the efficient numerical analysis of the fracture failure of elastic-brittle materials under thermal shock, the present invention innovatively proposes an adaptive multi-level phase-field fracture simulation method (English name: Adaptive multi-level phase-field model, English abbreviation: AMLPFM). The purposes are as follows: First, in order to overcome the deficiency that traditional damage models are difficult to accurately capture the crack propagation path and explore the fracture behavior under thermal shock, the present invention derives a phase-field fracture model for the fracture failure behavior of elastic-brittle materials under thermo-mechanical coupling based on the variational principle of energy functional, and uses an alternating iteration algorithm to solve the coupled control equations. Second, aiming at the fine grid requirements of the fracture phase field and the problem of hanging nodes existing in traditional adaptive strategies, the present invention proposes a novel adaptive multi-level grid discretization scheme, which can be more simply implemented numerically. During the simulation process, based on the advancement of the phase-field contour, the crack tip is tracked to identify the refinement region, the grid is dynamically updated in real time, the calculation cost is reduced, and the calculation efficiency is improved. In addition, the present invention proposes three prefabricated crack treatment techniques (split geometric model, weakly applied Dirichlet boundary condition, and weakening the material properties at the crack) to accurately consider the influence of various prefabricated cracks within the framework of multi-level grids. Finally, the present invention aims to solve the disadvantages of high calculation cost caused by fine grids when the existing finite element-phase field method analyzes the fracture failure problem under thermal shock, and the limitation of complex numerical implementation caused by the problem of hanging nodes existing in previous adaptive methods.

[0008] The technical solution of the present invention:

[0009] An adaptive multi-level grid phase-field method for elastic-brittle fracture analysis under thermal shock, an adaptive multi-level grid phase-field thermal fracture simulation method (AMLPFM). In order to accurately capture the high-gradient change of thermal stress during thermal shock and reduce the calculation cost in fracture simulation analysis, the present invention combines a new adaptive method based on multi-level grids (The multi-level hp-method) with a phase-field fracture model under thermo-mechanical coupling to predict the initiation and propagation of cracks under thermal-mechanical-displacement loads, and provides an adaptive multi-level grid phase-field thermal fracture simulation method. The specific steps are as follows:

[0010] Consider an arbitrary solid domain (computational domain) Ω under thermo-mechanical coupling, which contains a sharp crack Γ at time t ∈ [0, t a , as Figure 1 shown. Combining the existing phase-field fracture model under thermo-mechanical coupling, the following control equations can be obtained,

[0011]

[0012] Among them, equations (1), (2), and (3) are the displacement fields u = Ω × [0, t a, temperature field $\theta=\Omega\times[0,t a and phase field $d=\Omega\times[0,t a \to [0,1]$ of the governing equations, equations (4)-(6) are the corresponding boundary conditions; $\sigma$ is the Cauchy stress, $\psi$ is the strain energy density of the tensile part, $\rho$ is the density, $c$ is the specific heat capacity, $\dot{\theta}$ is the rate of change of temperature with time, $\gamma$ is the internal heat source, $G c is the critical fracture energy density, $l_0$ is the crack width used to control the width of the smeared crack band, $b$ is the body force, $J$ is the heat flux, $n$ is the unit outer normal vector of the solid domain boundary The boundary conditions for displacement, surface force, temperature, and heat flux are respectively and The displacement boundary surface force temperature and heat flux

[0013] Considering the strain spectral decomposition format given by Miehe et al., the strain energy density $\psi e is decomposed into a tensile part and a compressive part and the tensile part is weakened, that is,

[0014]

[0015] where $\varepsilon e is the elastic strain, $w(d)=(1 - d) 2 +\xi$ is the degradation function, and $0\lt\xi\ll1$ is added to obtain better numerical stability.

[0016] Then, the Cauchy stress $\sigma$ can be derived as

[0017]

[0018] In the thermo-mechanical coupling framework, the relationship between the elastic strain $\varepsilon e and the thermal strain $\varepsilon θ is as follows

[0019] $\varepsilon e =\varepsilon-\varepsilon θ , (9)

[0020] where, under the small strain assumption, the total strain $\varepsilon θ =\alpha(\theta-\theta_0)I$ is the thermal strain expression, $\alpha$ is the coefficient of thermal expansion, $\theta_0$ is the initial temperature, and $I$ is the second-order unit tensor.

[0021] To satisfy the irreversibility of crack growth, a historical maximum strain energy $H$ is introduced to replace

[0022]

[0023] Thus, it prevents non-physical crack healing caused by the decline.

[0024] According to Fourier's law, the heat flux J is

[0025]

[0026] where, to ensure no heat flux passes through the crack surface, the thermal conductivity k d is also degraded by the phase field value d

[0027] k d = w(d)k0, (12)

[0028] k0 is the thermal conductivity of the undamaged material. When d = 0, the heat transfer ability of the material is not affected, and when d = 1, the material no longer has the heat transfer ability.

[0029] The above is the phase field fracture model under thermo-mechanical coupling. When using the grid-like method for solution, extremely fine grids need to be divided at the crack, resulting in a large computational cost. It is necessary to combine the adaptive subdivision method. In the present invention, a novel adaptive strategy (the multi-level hp - finite element method, The multi-level hp-FEM) is adopted. During the simulation, the fine grids are subdivided in real time as the crack propagates, and there is no hanging node problem, and the numerical implementation is simple, so as to carry out efficient and effective numerical simulation analysis on the fracture failure behavior under thermal shock.

[0030] Different from the previous adaptive methods, in the present invention, the idea of superposition is adopted to achieve adaptive subdivision. Here, the computational domain Ω is discretized into a base grid using a coarse grid, and locally discretized into multiple layers of fine grids at the crack path, which is called the covering grid, as Figure 2 shown. According to the superposition idea, the approximate physical value at any x0 can be expressed as

[0031]

[0032] where the subscripts b and o represent the base grid and the covering grid respectively, u b 、u o2 and represent the approximate physical values on the base grid and the covering grid respectively. When adopting this subdivision strategy, the linear independence and compatibility between the shape functions need to be considered to eliminate the hanging node problem in the previous adaptive methods.

[0033] In the present invention, the above adaptive subdivision strategy with the idea of superposition is introduced into the control equations (1)-(6) of the phase field fracture model under thermo-mechanical coupling. The overall computational domain Ω is discretized into the base grid Ωb with the multi-layer covering grid Ω on , where n = 2, …, n k . The temperature field θ, displacement field u, phase field d, and their spatial derivatives are expressed as,

[0034]

[0035] The approximate temperature values θ b 、θ on in each layer of the grid after discretization can be interpolated as,

[0036]

[0037] where η · is the local coordinate under each layer of the grid, and N · N 、N · E and N · F are the shape functions attached to points, lines, and surfaces. The temperature values to be solved on these shape functions are θ · N 、θ · E and θ · F . For the displacement values u b 、u on and the phase field values d b 、d on a consistent interpolation format can be adopted.

[0038] Considering arbitrary virtual variables δθ b , δu b , δd b , δθ on , δu on and δd on , the following weak forms can be obtained from the governing equations (1)-(6),

[0039]

[0040] where the virtual forces of the temperature field and displacement field are,

[0041]

[0042] Based on the above theoretical derivation, the basic format of the coupled temperature-displacement-phase field under the multi-level hp grid is obtained. By applying quasi-static slow loading and then adopting the Newton-Raphson staggered solution strategy to simulate the crack initiation and propagation process, the adaptive multi-level grid phase field method for thermo-elasto-brittle fracture analysis proposed in the present invention can be realized.

[0043] The present invention is implemented through MATLAB software programming, and the result visualization is realized by using ParaView software. Combining with Figure 3 the calculation cycle schematic diagram of the AMLPFM method shown in

[0044] Step 1: Establish a multi-level discrete grid model and a finite element interpolation format, define material parameters, consider the linear independence and compatibility of the shape functions, and judge the activation states of the topological components (points, lines, and surfaces) in each layer of the grid;

[0045] Step 2: According to the position of the prefabricated crack, initialize the phase field value by using a split geometry model, weakly applying a phase field Dirichlet boundary condition, or weakening the material properties at the crack;

[0046] Step 3: Implement an alternating solution strategy using the Newton-Raphson iteration method in the current time step, and update the temperature field, displacement field, and phase field in sequence. When the iteration converges, obtain the temperature value, displacement value, and phase field value at the current time step;

[0047] Step 4: Store and output the relevant variable information, return to Step 3, enter the next time step, and continue until the calculation is completed.

[0048] In Step 1, the discretization is performed using quadrilateral elements, and the area near the tip or path of the prefabricated crack is subdivided into multiple layers of grids, as shown in Figure 4 . For simplicity, the symbols m are used to represent the base grid b and the overlay grid o. Then, the temperature field θ m , displacement field u m , and phase field d m of each layer of the grid after discretization can be expressed as

[0049]

[0050] where and are the shape functions of the temperature, displacement, and phase field of the corresponding topological components respectively, and are the shape function matrices of the three fields respectively. The temperature value θ mA , displacement value u mA , and phase field value d mA attached to the topological components are respectively assembled into vectors a m , and

[0051] Then, the spatial derivatives of the three fields can be further expressed as

[0052]

[0053] where and be

[0054]

[0055] wherein, and respectively represent taking partial derivatives with respect to the x-axis and y-axis.

[0056] Considering the above discretization and combining with the weak form (17), the residual equations for the three fields are obtained as follows.

[0057]

[0058] wherein, and are the residuals of temperature, displacement, and phase field under each layer of mesh. The subscripts are replaced by b and o with m. Then, the temperature field residual displacement field residual and phase field residual are expressed as,

[0059]

[0060] Meanwhile, the external force of the temperature field and the external force of the displacement field are,

[0061]

[0062] In the framework of the multi-level hp finite element, it is necessary to ensure the linear independence and compatibility among the shape functions. Linear independence requires that the shape functions of the coarse-scale mesh cannot be linearly combined by the shape functions of the fine-scale mesh. Compatibility requires that the superimposed approximate solution has C 0 continuity, which can be simply achieved by deactivating (i.e., removing) the corresponding mesh topology components in the numerical implementation.

[0063] In step 2, the present invention selects a suitable method to initialize the phase field according to the position of the prefabricated crack. As Figure 4 shown, there are two types of prefabricated cracks. Crack 1 coincides with the element edges of each level of mesh, and crack 2 obliquely passes through the interior of the element. Next, the influence of the prefabricated crack will be considered through three methods:

[0064] Method I: Split the geometric model. For prefabricated crack 1, only the area near the crack tip is locally refined, and the mesh on the crack path is not otherwise processed. Then, prefabricated crack 1 can be simply processed by splitting the geometric model. It should be noted that Method I is only applicable to the first type of prefabricated crack.

[0065] Method II: Apply the phase field Dirichlet boundary condition. By ensuring that the phase field value d = 1 at the prefabricated crack during the simulation process, the influence of any crack can be considered. However, for the crack 2 that obliquely passes through the interior of the element, the Dirichlet boundary condition cannot be directly applied to the corresponding topological component to constrain the phase field value inside the element. Here, the penalty function method is used to weakly apply the phase field Dirichlet boundary condition inside the element, and the steps are as follows:

[0066] First, select several points along the prefabricated crack and set the phase field value d = 1 at these points. Then the corresponding constraint conditions are as follows,

[0067]

[0068] where G is the shape function matrix of the base grid and the covering grid at the selected points, and V is the vector filled with element 1;

[0069] Then, considering the above constraint equations into the pure phase field diffusion model, the initial phase field distribution can be obtained, and its solution format is,

[0070]

[0071] where α = 1×10 5 is the penalty factor, and K d is the overall phase field stiffness matrix;

[0072] Finally, when using the Newton - Raphson iteration method in subsequent time steps, make V = 0 to ensure that the phase field value on the prefabricated crack always remains d = 1;

[0073] Method III: Weaken the material properties. For the prefabricated crack 2, the material properties at the prefabricated crack can also be weakened, which is equivalent to making the material near the crack band extremely soft, and the effect similar to Method I can be achieved. Here, a very small scalar β is introduced to weaken the elastic matrix of the crack band, that is,

[0074] D * = βD, (31)

[0075] where 0 << β < 1;

[0076] Both Method II and Method III can handle various prefabricated cracks. During the solution process, Method II introduces a penalty term to the solution format, while Method III does not require additional operations.

[0077] In step 3, the present invention adopts the Newton - Raphson iteration method to implement the staggered solution strategy in the current time step [t l , t l+1 , and updates the temperature field, displacement field, and phase field in turn. Among them, for the temperature field residual equation (24), It is expressed by the backward difference method as

[0078]

[0079] where Δt = t l+1 - t l is the time increment; the coupled nonlinear equations (24)-(26) are solved by the Newton-Raphson iteration method, and the solution format at the i-th iteration step is as follows:

[0080] First, fix the displacement field and the phase field at the previous iteration step, and solve the temperature field a at the current iteration step according to Equation (24) (i) ,

[0081]

[0082] where and are the temperature stiffness matrices coupled by the base grid, the covering grid and the two-level grid respectively, and the expressions are

[0083]

[0084] where the thermal conductivity degrades due to material damage;

[0085] Next, fix the temperature field a (i) at the current iteration step and the phase field at the previous iteration step, and solve the displacement field

[0086]

[0087] where the elastic matrix the stiffness matrix and have expressions similar to ; it should be noted that when Method III is used to handle the prefabricated crack, the elastic matrix near the crack band is weakened, i.e., D * = βD;

[0088] Finally, fix the temperature field a (i) at the current iteration step and the displacement field and solve the phase field

[0089]

[0090] where the phase field stiffness matrix and The expression of is similar to that of

[0091] It should be noted that when dealing with the influence of prefabricated cracks by Method II, there are constraint conditions This constraint equation is considered into the phase-field solution format (39) through the penalty function method, that is

[0092]

[0093] Combined with Figure 3 , the specific implementation process of the adaptive multi-level grid phase-field method (AMLPFM) for elastic-brittle fracture analysis under thermal shock proposed by the present invention will be shown in the following form of pseudocode:

[0094] 1), Define parameters (density ρ, elastic modulus E, Poisson's ratio ν, thermal conductivity k0, specific heat capacity c, critical fracture energy density G c ), crack characteristic width l0, time increment Δt, penalty factor α, weakening coefficient β, etc.;

[0095] 2), Discretize the finite element model of the multi-level hp grid, and select a suitable method (split geometric model, weakly apply phase-field Dirichlet boundary conditions or weaken material properties) according to the position of the prefabricated crack, and initialize the distribution of phase-field values;

[0096] 3), According to the discretized finite element model, judge the activation state of the grid topology components;

[0097] 4) Time step loop l = 0, 1, 2,..., N a ;

[0098] 4.1), Initialize the iteration variable i = 0, and let a (i) ←a (l) , and

[0099] 4.2), Start the iterative solution;

[0100] 4.2.1), i←i + 1;

[0101] 4.2.2), According to Equation (33), fix the displacement field and the phase field and update the temperature field

[0102] 4.2.3), According to Equation (37), fix the temperature field a (i) and the phase field and update the displacement field

[0103] 4.2.4), according to Equation (39) (if Equation (41) is used to process the prefabricated crack, then according to Equation (41)), fix the temperature field a (i) and the displacement field Update the phase field

[0104] 4.2.5), if and Go to step 4.3), otherwise, return to step 4.2) to continue the iteration;

[0105] 4.3), Update the temperature, displacement and phase field

[0106] 4.4), Output the calculation file of the current time step and perform post-processing;

[0107] 5), Return to step 4) until the calculation ends;

[0108] In the above implementation steps, ∈ is the convergence tolerance, and take ∈ = 1×10 -5 .

[0109] Advantages of the present invention:

[0110] (1) The adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock provided by the present invention provides a simple and efficient numerical simulation method for thermal shock fracture analysis, and broadens the application scope of multi-level hp finite elements in the field of fracture analysis. Through the phase field model, the deficiency that it is difficult to accurately capture the crack propagation path by using the traditional continuous medium model can be effectively overcome, and complex crack propagation problems can be relatively simply processed, such as crack crossing, bifurcation, and free crack propagation in three-dimensional space, etc. Moreover, the present invention can further extend the constitutive model of materials to the fracture failure analysis of other materials, such as large deformation fracture analysis, plastic material fracture analysis, etc.;

[0111] (2) The efficient adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock provided by the present invention. In this method, multi-level hp finite elements are used as an adaptive strategy, which can reduce computing resources and improve computing efficiency. Compared with the previous adaptive methods, multi-level hp finite elements achieve local refinement by superimposing grids, and there is no hanging node problem, making the numerical implementation simple. At the same time, this method makes full use of the advantages of h refinement and p refinement, can effectively capture high-gradient changes, and is suitable for fracture simulation analysis. In addition, this new adaptive strategy can also be extended to other large-scale mechanical problems, such as numerical analysis of complex structures, topology optimization, and aerodynamics, etc.;

[0112] (3) The efficient adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock provided by the present invention develops three multi-level crack treatment techniques under the multi-level grid framework. In Method II, the prefabricated crack is used as a boundary condition, and the solution format is corrected by the penalty function method to weakly apply the boundary condition; Method III weakens the material properties at the crack to achieve an effect similar to Method I. Its essence is to regard the crack as a virtual domain, which can be further extended to the virtual domain embedding technology for complex geometric bodies, and efficient numerical simulation analysis is carried out for complex structures such as lattices and heat sinks.

[0113] (4) The efficient adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock provided by the present invention uses the Newton-Raphson iteration method to solve the temperature field-displacement field-phase field alternately, greatly reducing the difficulty of numerical implementation, improving the calculation efficiency, and providing a feasible solution for the parallel computing research of large-scale complex fracture problems. Description of the Drawings

[0114] Figure 1 It is a schematic diagram of the fracture problem under thermo-mechanical coupling of an arbitrary solid domain Ω containing a sharp crack interface Γ and a diffuse crack of the present invention. Among them, (a) is a sharp crack, and (b) is a diffuse crack;

[0115] Figure 2 It is a schematic diagram of the multi-level hp finite element method of the present invention. Among them, (a) is a one-dimensional schematic diagram, and (b) is a two-dimensional schematic diagram;

[0116] Figure 3 It is an operation flow chart of an efficient adaptive multi-level grid phase field method (AMLPFM) for elastic-brittle fracture analysis under thermal shock of the present invention;

[0117] Figure 4 It is a schematic diagram of the multi-level grid at two prefabricated cracks of the present invention;

[0118] Figure 5 It is the Gaussian integration point distribution in Method III for prefabricated crack treatment of the present invention;

[0119] Figure 6 It is a schematic diagram of the structure and boundary conditions of Example 1 of quasi-static tensile fracture of a prefabricated double-inclined crack rectangular plate of the present invention;

[0120] Figure 7 It is a cloud diagram of the crack phase field evolution of Example 1 of the present invention using Method II and Method III to treat prefabricated cracks. Among them, (a) is Method II, and (b) is Method III;

[0121] Figure 8 It is a comparison diagram of the force-displacement curves of Example 1 of the present invention using Method II and Method III to treat prefabricated cracks;

[0122] Figure 9 This is the phase-field evolution nephogram and experimental result diagram of the quenching fracture of the ceramic plate under thermal shock in Embodiment 2 of the present invention. Among them, (a) is at t = 0.5 ms, (b) is at t = 5 ms, (c) is at t = 30 ms, (d) is at t = 75 ms, (e) is at t = 200 ms, and (f) is the experimental result. Detailed implementation manners

[0123] The following further illustrates the detailed implementation manners of the present invention in combination with the attached drawings and technical solutions.

[0124] Embodiment

[0125] Combined with Figures 6 to 9 The accuracy, reliability, and excellent performance of the efficient adaptive multi-level grid phase-field method (AMLPFM) for the analysis of elastic-brittle fracture under thermal shock proposed by the present invention are further described in detail.

[0126] First, a quasi-static tensile brittle fracture simulation of a square plate with an inclined initial crack will be implemented to illustrate the reliability of the developed pre-crack treatment technology, as well as the effectiveness and efficiency of the adaptive grid-phase field method proposed by the present invention for solving fracture problems; then, the thermal cracking of the ceramic plate will be simulated, and the crack path will be compared with the experimental results to verify the accuracy of the AMLPFM for simulating the fracture failure behavior under thermal shock. The reference solution in the first embodiment comes from the work of Kim et al., and the experimental results in the second embodiment come from the work of Jiang et al. In all embodiments, quadrilateral elements are used, and the phase-field model parameter l0 is set to be greater than or equal to half of the grid size h e / 2 to meet the rationality of the numerical simulation results of phase-field fracture. All embodiments are considered under the plane strain assumption, and the gravity effect is not considered.

[0127] 1) Embodiment 1: Quasi-static tensile fracture of a rectangular plate with prefabricated double inclined cracks (attached Figures 6 - 8 )

[0128] Embodiment 1 simulates the tensile fracture of a rectangular plate with double inclined cracks as Figure 6 shown. The grid is subdivided into 3 levels at the crack, and Method II (weakly applying the phase-field boundary condition) and Method III (weakening the material properties at the crack) are respectively used to process the prefabricated crack. The lower boundary of the geometry is fixed, and a vertical displacement load is applied to the upper boundary. The material parameters are shown in Table 1 below.

[0129] Table 1. Material parameter table for quasi-static tensile brittle fracture of a rectangular plate with prefabricated double inclined cracks

[0130]

[0131] Figure 7 The phase-field evolution contour maps of the prefabricated cracks processed by Method II and Method III are given. Before the displacement loading in Method II, the phase-field value d = 1 at the prefabricated crack is ensured, so that the material at the crack does not have the bearing capacity; in Method III, the phase-field value is not initialized, and the phase-field value at the prefabricated crack gradually changes from 0 to 1 with the displacement loading. Both methods obtain the same crack propagation path, and the mesh is dynamically updated by tracking the crack tip during the simulation, further reducing the total degrees of freedom, making the AMLPFM method have great potential to solve large-scale complex fracture problems. In Figure 8 the force-displacement curves of the two methods are given, and the work of Kim et al. is taken as the reference solution. It can be seen that the three curves basically coincide, verifying the reliability of the results. Through the above analysis, Example 1 proves the effectiveness and efficiency of the adaptive mesh-phase field method proposed in the present invention, as well as the accuracy of the three prefabricated crack treatment techniques.

[0132] (2) Example 2: Simulation of quenching fracture of ceramic plates under thermal shock (attached Figure 9 )

[0133] Example 3 simulates the quenching fracture phenomenon of ceramic plates and compares it with the experimental results. After the ceramic plate is heated to 300 °C, it is thrown into a pool of water at 20 °C. The rapid cooling of the ceramic plate causes crack initiation and propagation. After being fished out and dyed with ink, the internal crack propagation mode can be observed, and the results are as shown in Figure 9 (f). Considering the symmetry of the geometry, only a quarter model is used for numerical simulation here, and its material parameters are shown in Table 2 below.

[0134] Table 2. Material parameter table for quenching fracture of ceramic plates under thermal shock

[0135]

[0136] Figure 9 The evolution contour map of the quenching fracture of the ceramic plate is given and compared with the experimental results. It can be seen that the mesh will be dynamically and adaptively refined as the crack propagates, and numerous cracks are generated inside the ceramic plate, which is close to the experimental results. This verifies that the AMLPFM proposed in the present invention can efficiently and effectively conduct numerical simulation research on the fracture failure behavior of elastic-brittle materials under thermal shock. As shown in Figure 9 (a), at the beginning of the simulation, a consistent damage zone will be generated on the outer boundary of the ceramic plate, and no crack initiation occurs; then in 9(b), tiny cracks ① with equal distance and equal length rapidly initiate on the boundary; then as shown in Figure 9 (c), about half of the tiny cracks ① continue to expand inward to crack ②, while the other half of the cracks ① stop expanding due to insufficient strain energy supply; this phenomenon is shown in Figure 9(d) occurs again until the simulation stops after the crack ③ is extended. The final crack path is as shown in Figure 9 (e), and the result is close to the experimental result. This fracture process is extremely short. It is difficult to observe the crack evolution process through experiments, while numerical simulation can be used to more easily and deeply explore the crack propagation process under thermal shock. This further demonstrates the practicality of the AMLPFM method of the present invention and has certain scientific and engineering practical significance.

[0137] In summary, the above two embodiments have respectively verified the accuracy and effectiveness of the AMLPFM method proposed by the present invention in detail from different levels, and demonstrated the significant advantages of this method in terms of computational efficiency and computational scale, proving the necessity of its invention. At the same time, using the multi-level hp finite element as an adaptive strategy, its mesh is dynamically refined as the crack propagates. Compared with other adaptive strategies, there is no hanging node problem, and the numerical implementation is simple. It can also be used in other numerical simulations with high gradient changes. Therefore, the efficient adaptive multi-level grid phase field method (AMLPFM) for elastic-brittle fracture analysis under thermal shock proposed by the present invention is a high-performance numerical method with great development prospects.

[0138] The embodiments of the present invention are given for purposes of illustration and description, and are not exhaustive or limit the invention to the disclosed form. Many modifications and variations are obvious to those of ordinary skill in the art. The embodiments are chosen and described in order to better illustrate the principles and practical applications of the present invention, and to enable those of ordinary skill in the art to understand the present invention and design various embodiments with various modifications suitable for specific purposes.

Claims

1. An adaptive multi-level grid phase field method for analyzing elastic-brittle fracture under thermal shock, which is an efficient numerical method for analyzing the fracture failure problem of elastic-brittle materials under thermal shock. It is characterized in that, The steps are as follows: An arbitrary solid domain under thermo-mechanical coupling, taking the solid domain as the computational domain Ω, contains a sharp crack Γ within the time t ∈ [0, t a a ; combining the existing phase-field fracture model under thermo-mechanical coupling, the following governing equations are obtained, wherein, equations (1), (2) and (3) are the governing equations of the displacement field u = Ω×[0, t a , the temperature field θ = Ω×[0, t a , and the phase field d = Ω×[0, t a → [0, 1], and equations (4)-(6) are the corresponding boundary conditions; σ is the Cauchy stress, is the strain energy density of the tensile part; ρ is the density, c is the specific heat capacity, is the rate of change of temperature with time, γ is the internal heat source, G c is the critical fracture energy density, l0 is the crack width used to control the width of the smeared crack band, b is the body force, J is the heat flux, n is the unit outer normal vector of the boundary θΩ of the solid domain; the boundary conditions for displacement, surface force, temperature, and heat flux are respectively and The displacement boundary surface force temperature and heat flux According to the given format of strain spectrum decomposition, decompose the strain energy density ψ e into a tensile part and a compressive part and weaken the tensile part, that is, where ε e is the elastic strain, w(d) = (1 - d) 2 + ξ is the degradation function, and 0 < ξ << 1 is added to obtain better numerical stability; Then, the Cauchy stress σ can be derived as Under the thermo-mechanical coupling framework, the elastic strain ε e and the thermal strain ε θ have the following relationship: ε e = ε - ε θ (9) Among them, under the small strain assumption, the total strain ε θ = α(θ - θ0)I is the expression of thermal strain, where α is the coefficient of thermal expansion, θ0 is the initial temperature, and I is the second-order unit tensor; To satisfy the irreversibility of crack growth, a historical maximum strain energy H is introduced to replace the Thereby preventing non-physical crack healing caused by a decrease; According to Fourier's law, the heat flux J is Among them, to ensure that no heat flux passes through the crack surface, the thermal conductivity k d is also degraded by the phase field value d k d = w(d)k0 (12) where k0 is the thermal conductivity of the undamaged material. When d = 0, the heat transfer ability of the material is not affected, while when d = 1, the material no longer has the heat transfer ability; An adaptive strategy is adopted to solve the above phase-field fracture model under thermo-mechanical coupling, which is specifically as follows: The superposition idea is used to achieve adaptive refinement. The computational domain Ω is discretized into a base grid with a coarse grid, and locally discretized into multiple layers of fine grids at the crack path, which is called the covering grid; According to the superposition idea, the approximate physical value at any x0 is expressed as Among them, the subscripts b and o represent the base grid and the overlay grid respectively, and u b , u o2 and represent the approximate physical values on the base grid and the overlay grid respectively; Introduce the above adaptive subdivision strategy with the superposition idea into the governing equations (1)-(6) of the phase-field fracture model under thermo-mechanical coupling. The overall computational domain Ω is discretized into the base mesh Ω b and the multi-layer covering mesh Ω on , where n = 2, …, n k ; the temperature field θ, displacement field u, phase field d, and their spatial derivatives are expressed as The approximate values of temperature θ in each layer of the discretized grid b and θ on can be interpolated as Among them, η. is the local coordinate under each layer of grid, N. N , N. E and N. F are the shape functions attached to points, lines and surfaces, and the temperature values to be solved on these shape functions are θ. N , θ. E and θ. F ; for the displacement values u b , u on and the phase field values d b , d on adopt a consistent interpolation format; Consider arbitrary virtual variables δθ b , δu b , δd b , δθ on , δu on and δd on , the following weak forms are obtained from the governing equations (1)-(6): where the virtual forces of the temperature field and displacement field are Based on the above theoretical derivation, the basic format of coupling temperature-displacement-phase field under multi-level hp grids is obtained. Through quasi-static slow loading, and then adopting the Newton-Raphson staggered solution strategy, the crack initiation and propagation processes are simulated to realize the proposed adaptive multi-level grid phase-field method for elastic-brittle fracture analysis under thermal shock.

2. The adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock according to claim 1, wherein The adaptive multi-level grid phase-field for elastic-brittle fracture analysis under thermal shock is realized by programming with MATLAB software, and the result visualization is realized by using ParaView software. The specific steps are as follows: Step 1: Establish a multi-level discrete grid model and a finite element interpolation format, define material parameters, consider the linear independence and compatibility of the shape functions, and judge the activation states of topological components in each layer of the grid, where the topological components include points, lines, and faces; Step 2: According to the position of the prefabricated crack, initialize the phase field value by using the split geometry model, weakly applying the phase-field Dirichlet boundary condition, or weakening the material properties at the crack; Step 3: Implement the staggered solution strategy by using the Newton-Raphson iteration method in the current time step, and update the temperature field, displacement field, and phase field in turn. When the iteration converges, obtain the temperature value, displacement value, and phase field value at the current time step; Step 4: Store and output relevant variable information, return to Step 3, enter the next time step, and continue until the calculation is completed.

3. The adaptive multi-level grid phase-field method for elastic-brittle fracture analysis under thermal shock according to claim 2, wherein In Step 1, quadrilateral elements are used for discretization, and are subdivided into multiple layers of grids at the tip or path of the pre-crack; let the symbol m represent the base grid b and the covering grid o; then, the temperature field θ m , displacement field u m and phase field d m are expressed as Among them, and are the shape functions of the temperature, displacement, and phase field of the corresponding topological components respectively, and are the shape function matrices of the three fields respectively. The temperature value θ mA , displacement value u mA and phase field value d mA are respectively assembled into vectors a m , and Then the spatial derivatives of the three fields are further expressed as Among them, and are Among them, and respectively represent taking partial derivatives with respect to the x-axis and the y-axis; Considering the above discretization treatment and combining with the weak form (17), the residual equations of the three fields are obtained respectively as Among them, and are the residuals of temperature, displacement and phase field under each layer of grid. The subscripts b and o are replaced by m. Then the temperature field residual in each layer of grid displacement field residual and phase field residual The expressions are Meanwhile, the external force of the temperature field and the external force of the displacement field are 4. The adaptive multi-level grid phase field method for analyzing elastic-brittle fracture under thermal shock according to claim 2, characterized in that In Step 2, according to the position of the prefabricated crack, the phase field is initialized. There are two types of prefabricated cracks. The prefabricated crack 1 coincides with the element edge of each level of the grid, and the prefabricated crack 2 obliquely passes through the interior of the element; The influence of the prefabricated crack will be considered by three methods: Method I: Split geometry model. For the prefabricated crack 1, only local refinement is performed in the crack tip region, and the grids on the crack path are not otherwise processed. The prefabricated crack 1 is processed by the split geometry model; Method I is only applicable to the first type of prefabricated crack; Method II: Apply the phase field Dirichlet boundary condition. By ensuring that the phase field value d = 1 at the prefabricated crack during the simulation, the influence of any crack can be considered. However, for the prefabricated crack 2 that obliquely passes through the interior of the element, the Dirichlet boundary condition cannot be directly applied to the corresponding topological member to constrain the phase field value inside the element. The penalty function method is used to weakly apply the phase field Dirichlet boundary condition inside the element, and the steps are as follows: First, select several points along the prefabricated crack and set the phase field value d = 1 at these points. Then the corresponding constraint conditions are as follows, where G is the shape function matrix of the base grid and the covering grid at the selected points, and V is the vector filled with element 1; Then, consider the above constraint equation into the pure phase field diffusion model to obtain the initial phase field distribution, and its solution format is, Among them, α = 1×10 5 is the penalty factor, and K d is the overall stiffness matrix of the phase field; Finally, when using the Newton - Raphson iteration method in subsequent time steps, make V = 0 to ensure that the phase field value on the prefabricated crack always remains d = 1; Method III: Weaken the material properties. For the prefabricated crack 2, the effect of Method I can also be achieved by weakening the material properties at the prefabricated crack, making the material near the crack band particularly soft. Here, a very small scalar β is introduced to weaken the elastic matrix of the crack band, that is, D * = βD (31) where 0 < β < 1; Both Method II and Method III can handle various prefabricated cracks. During the solution process, Method II introduces a penalty term to the solution format, while Method III does not require additional operations.

5. The adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock according to claim 2, characterized in that In step 3, in the current time step [t l , t l+1 , the staggered solution strategy is implemented using the Newton-Raphson iteration method to update the temperature field, displacement field, and phase field in sequence; among them, for in the temperature field residual equation (24) is expressed using the backward difference method as, where, Δt = t l+1 - t l is the time increment; the coupled non - linear equations (24)-(26) are solved by the Newton - Raphson iteration method, and the solution format at the i - th iteration step is as follows: First, fix the displacement field of the previous iteration step and the phase field Solve the temperature field a of the current iteration step according to Equation (24) (i) , Among them, and are the temperature stiffness matrices of the base grid, the covering grid, and the coupling of the two-level grid respectively, and the expression is Among them, the thermal conductivity coefficient degrades due to material damage; Next, fix the temperature field a at the current iteration step (i) and the phase field at the previous iteration step Solve for the displacement field at the current iteration step according to Equation (25) Among them, the elastic matrix stiffness matrix and have expressions similar to those of ; when the prefabricated crack is treated by Method III, the elastic matrix near the crack band is weakened, i.e., D * = βD; Finally, fix the temperature field a at the current iteration step (i) and the displacement field Solve the phase field according to Equation (26) Among them, the phase field stiffness matrix and have expressions similar to those of ; When dealing with the influence of prefabricated cracks by Method II, there are constraint conditions The constraint equation is considered in the phase-field solution format (39) by the penalty function method, that is

Images