ADAPTIVE MULTI-LEVEL PHASE-FIELD METHOD FOR BRITTLE FRACTURE OF ELASTIC MATERIALS UNDER THERMAL SHOCK

Abstract

The present invention belongs to the field of computational mechanics, and provides an adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock, thus to provide a new numerical calculation method for efficient fracture study. In the present invention, damage degree of a material is characterized by a phase field, and the multi-level hp-FEM is used as an adaptive strategy to implement dynamic discretization in a process of crack propagation under thermal shock. Compared with previous adaptive phase-field fracture models, the present invention has no hanging nodes, and has simple numerical implementation and high calculation efficiency, which can effectively predict the fracture behavior of a structure. At the same time, three pre-existing crack treatment techniques are developed to initialize the phase field. The temperature, displacement and phase fields are solved by an alternating minimization algorithm, which can stably and effectively solve the multi-field coupling large-scale fracture problems.

Description

TECHNICAL FIELD

The present invention belongs to the field of computational mechanics, and relates to an adaptive phase-field method for brittle fracture of elastic materials under thermal shock.

BACKGROUND

Thermal shock fracture exists widely in engineering practice and daily life, for example thermal cracking caused by rapid temperature changes in machining and manufacturing processes (such as additive manufacturing, quenching and welding), microcracks of a multiphase composite material produced during machining due to different thermal expansion coefficients, and self-explosion of window glass in hot weather. Therefore, the study of the failure mechanism under thermal shock is helpful to the safety assessment of a structural device. If an analysis is carried out by experimental means, it is difficult and expensive to implement the analysis, difficult to consider the influence of multiple external factors on a fracture behavior, and difficult to observe a crack evolution process. As a new high performance simulation computing technology, numerical simulation can help us to further explore the failure process under thermal shock, so as to predict the failure behavior of a material comprehensively and accurately. However, the study of multi-field coupling fracture simulation under thermal shock is still a very challenging topic, mainly because it requires comprehensive and detailed consideration of stable crack propagation simulation, multi-physical field coupling solution and effective integration of efficient adaptive methods.

For the numerical simulation of fracture, two types of methods are mainly used at present: discrete crack characterization methods and smeared crack characterization methods. Among which, the discrete crack characterization methods are to describe cracks by explicitly defining crack collection information, which include mesh type treatment methods such as the element deletion method, the cohesive zone model (CZM) method and the extended finite element method (XFEM). In a simulation process, a crack will propagate, and the geometry thereof will change with the load and time, so that a mesh needs to be redivided constantly in a solution process, which greatly consumes the computing resources. At the same time, due to the limitation of their own algorithm, the discrete crack characterization methods cannot treat complex fracture behaviors such as crack initiation, bifurcation and intersection easily and efficiently. Whereas the smeared crack characterization methods are to characterize cracks by field variables without the need of separate crack modeling or mesh redivision, which improves the computational efficiency and reduces the difficulty of numerical implementation to a certain extent, and has significant advantages in treating complex fractures.

In the smeared crack characterization methods, Bourdin et al. has developed a regularized approximate smeared fracture model, i.e. a phase field model, in “Numerical experiments in revisited brittle fracture” based on the fracture variational principle proposed by Francfort and Marigo. The model describes cracks by a phase field variable to characterize the damage value of a material, and has been widely used in treating crack initiation and propagation. Amor et al. has conducted volume-deviatoric strain decomposition of strain energy in “Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments”, assuming that crack evolution can only be driven by a deviatoric stress and a positive spherical stress. In 2010, Miehe et al. has re-established a phase-field fracture model in “A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits” under the framework of thermodynamics. Strain energy is decomposed into tensile and compressive parts by the strain tensor spectrum, and it is believed that crack evolution can only be driven by tensile deformation, so that unrealistic crack propagation during compression is prevented. In order to simulate a fracture process more accurately, Borden et al. has proposed a fourth-order phase-field fracture equation in “A phase-field description of dynamic brittle fracture” under the framework of the isogeometric analysis method to simulate the fracture of a brittle material. In order to treat the quasi-brittle fracture of a material more accurately, Professor Wu Jianying has introduced the concept of cohesive force into the phase-field fracture model and proposed the phase-field regularized cohesive zone model (PF-CZM), which is combined with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm to further improve the computational efficiency. When the phase field model is used for simulating crack evolution, it is not necessary to update the geometric topology of the model. The phase field model has advantages in treating complex fracture modes such as crack initiation, propagation and intersection, and is widely used in various types of fractures, such as thermal cracking, dynamic fracture, hyperelastic fracture, hydraulic fracture and shock fracture. Among which, the study on thermal cracking is helpful to the safety assessment of manufacturing processes such as additive manufacturing, quenching and casting, and has drawn the attention of numerous scholars. Tangella et al. has adopted the hybrid phase-field method in “Hybrid phase-field modeling of thermo-elastic crack propagation” to simulate the fracture of an elastic-brittle material. Ruan et al. has proposed a thermo-mechanical 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 the fracture mode, and applied the model to the simulation of small crack initiation in additive manufacturing. Mandal et al. has used PF-CZM to simulate the thermal cracking phenomenon of elastic-brittle material in “Fracture of thermo-elastic solids: Phase-field modeling and new results with an efficient monolithic solver”, and further improved the computational efficiency by the BFGS algorithm.

The phase field model has been widely used in the simulation of the fracture process, but the numerical simulation requires extremely fine mesh, which makes the calculation amount large. In order to meet the mesh requirements without significantly increasing the computational cost, an adaptive mesh refinement scheme is widely used in phase-field fracture simulation. Heister et al. has proposed a predictor-corrector mesh adaptivity method with h-refinement in “A primal-dual active set method and predictor-corrector mesh adaptivity for computing fracture propagation using a phase-field approach” to locally refine the mesh in real time during crack propagation. Patil et al. has proposed a phase-field fracture model under an adaptive multiscale finite element method (MsFEM), which can significantly reduce the computational cost and improve the computational efficiency. Recently, Zander et al. has proposed a new adaptive method in “Multi-level hp-adaptivity: high-order mesh adaptivity without the difficulties of constraining hanging nodes”, which is called the multi-level hp-FEM. The method uses superimposed fine meshes to realize local refinement. In most of the adaptive methods in the past, fine meshes are used instead of coarse meshes, which will cause hanging node problems and lead to complex numerical implementation. Whereas the multi-level hp-FEM uses the idea of refinement and the high-order integrated Legendre polynomial for shape function, so that the multi-level hp-FEM has no hanging nodes, and the numerical implementation is easier. At the same time, the advantages of h-refinement and p-refinement are fully utilized to allow high convergence in the process of treating the problem of high gradients. The characteristics are suitable for phase-field fracture simulation, and the cracks are locally covered by multi-level fine meshes, thus to reduce the consumption of the computing resources.

In the present invention, the multi-level hp-FEM is used as an adaptive strategy to conduct dynamic discretization of the thermo-mechanical phase-field fracture model, and three pre-existing crack treatment techniques are developed to ensure the accurate imposition of crack boundaries in the phase field. Moreover, the alternating minimization algorithm is adopted to effectively obtain temperature, displacement and phase fields. The study of an efficient numerical method for fracture analysis of an elastic-brittle material under thermal shock proposed by the present invention can capture the crack evolution process efficiently, has no hanging nodes in traditional adaptive methods, has easy numerical implementation, is used for simulating the elastic-brittle fracture mechanical behavior of a structure under a thermal shock load, and provides an effective way for the accurate analysis of engineering fractures.

SUMMARY

Aiming at the efficient numerical analysis of fracture of an elastic-brittle material under thermal shock, the present invention innovatively proposes an adaptive multi-level phase-field method (AMLPFM). The purposes are that: first, in order to overcome the deficiency of traditional damage model which is difficult to accurately capture a crack propagation path and explore the fracture behavior under thermal shock, the present invention deduces a phase-field fracture model aiming at the thermo-mechanical fracture behavior of an elastic-brittle material based on the energy functional variational principle, and uses the alternating minimization algorithm to solve the coupled governing equations. Second, aiming at the fine mesh requirements of a fracture phase field and the problem of hanging nodes in the traditional adaptive strategy, the present invention proposes a novel adaptive multi-level discrete scheme, which can be more easily implemented numerically. In a simulation process, based on the advance of a phase field contour, a crack tip is tracked to identify a refinement region, and a mesh is dynamically updated in real time to reduce the calculation cost and improve the calculation efficiency. In addition, the present invention provides three pre-existing crack treatment techniques (dividing the geometric model, applying the damage Dirichlet boundary condition weakly and weakening the material properties surrounding the pre-existing cracks) to accurately consider the influence of various pre-existing cracks under the multi-level framework. Finally, the present invention aims to solve the shortcomings of the existing finite element-phase field method that the calculation cost is high due to a fine mesh in analyzing the fracture problem under thermal shock, and the limitation of complex numerical implementation due to the problem of hanging nodes in the previous adaptive methods.

The technical solution of the present invention is as follows:

An adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock, and an adaptive multi-level phase-field method (AMLPFM) for thermal cracking simulation. In order to accurately capture the high gradients of a thermal stress during thermal shock and reduce the calculation cost in fracture simulation analysis, the present invention combines the multi-level hp-method with a thermo-mechanical phase-field fracture model to predict the initiation and propagation of a crack under a thermo-mechanical-displacement load, and provides an adaptive multi-level phase-field method. The specific steps are as follows:

Considering an arbitrary thermo-mechanical solid domain (computational domain) Ω, containing a sharp crack Γ in time t∈[0, t_a], as shown in FIG. 1. Combining with the existing thermo-mechanical phase-field fracture model, the following governing equations can be obtained,

$\begin{array}{cc}\begin{array}{cc}\nabla ·\sigma +b=0,& \mathrm{in}\Omega \times \left[0,t_a\right],\end{array}& \left(1\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}\mathrm{\rho c\theta }+\nabla ·J=\gamma ,& \mathrm{in}\Omega \times \left[0,t_a\right],\end{array}& \left(2\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}-2\left(1-d\right){\psi }_e^{+}+\frac{G_c}{l_0}\left(d-l_0^2\Delta d\right)=0,& \mathrm{in}\Omega \times \left[0,t_a\right],\end{array}& \left(3\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}\sigma ·n=\stackrel{\_}{t},& \mathrm{in}\partial {\Omega }_t\times \left[0,t_a\right],\end{array}& \left(4\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}J·n=\stackrel{\_}{q},& \mathrm{in}\partial {\Omega }_q\times \left[0,t_a\right],\end{array}& \left(5\right)\end{array}$ $\begin{array}{cc}\begin{array}{cc}\nabla d·n=0,& \mathrm{in}\partial \Omega \times \left[0,t_a\right],\end{array}& \left(6\right)\end{array}$

Wherein equations (1), (2) and (3) are the governing equations of a displacement field u=Ω×[0, t_a], a temperature field θ=Ω×[0, t_a] and a phase field d=Ω×[0, t_a]→[0,1], respectively, and equations (4)-(6) are the corresponding boundary conditions; σ is the Cauchy stress and is the strain energy density of a tensile part, ρ is the density, c is the specific heat capacity, {dot over (θ)} is the rate at which temperature changes over time, γ is the internal thermal source, G_c is the critical energy release rate, l₀ is the crack width used for controlling the width of a smeared crack band, b is the physical strength, J is the heat flux, and n is the outward normal vector of the solid domain boundary ∂Ω; displacement, traction, temperature and heat flux boundary conditions are ∂Ω_u, ∂Ω_t, ∂Ω_θ and ∂Ω_q, respectively, and displacement ū, traction t, temperature θ and heat flux q are applied accordingly.

Considering the strain spectrum decomposition format given by Miehe et al., decomposing the strain energy density ψ_e into the tensile part ψ_e⁺ and a compressive part ψ_e⁻, and weakening the tensile part, that is,

$\begin{array}{cc}{\psi }_e\left({\epsilon }_e\left(u,\theta \right)\right)=w\left(d\right){\psi }_e^{+}+{\psi }_e^{-},& \left(7\right)\end{array}$

Wherein ε_e is the elastic strain, w(d)=(1−d)²+ξ is the degradation function, 0<ξ<<1 is added to obtain a better numerical stability.

Then the Cauchy stress σ can be derived as follows,

$\begin{array}{cc}\sigma =w\left(d\right)\frac{\partial {\psi }_e^{+}}{\partial {\epsilon }_e}+\frac{\partial {\psi }_e^{-}}{\partial {\epsilon }_e}& \left(8\right)\end{array}$

The relationship between the elastic strain ε_e and the thermal strain ε_θ under the thermo-mechanical framework is as follows,

$\begin{array}{cc}{\epsilon }_e=\epsilon -{\epsilon }_{\theta },& \left(9\right)\end{array}$

Wherein the total strain under the assumption of small strain is

$\epsilon =\frac{1}{2}\left(\nabla u+{\nabla }^{T}u\right),$

ε_θ=α(θ−θ₀)I is the expression of the thermal strain, α is the thermal expansion coefficient, θ₀ is the initial temperature, and I is the identity tensor.

In order to satisfy the irreversibility of crack growth, a historical maximum strain energy H is introduced to replace ψ_e⁺ in formula (3),

$\begin{array}{cc}H=\underset{\tau \in \left[0,t\right]}{\max }\left\{{\psi }_e^{+}\left({\epsilon }_e,\tau \right)\right\}& \left(10\right)\end{array}$

Thus to prevent non-physical crack healing due to the drop of ψ_e⁺.

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

$\begin{array}{cc}J=-k_d·\nabla \theta ,& \left(11\right)\end{array}$

Wherein, in order to ensure no heat flux across the crack surface, the inherent thermal conductivity k_d is also degraded by the phase field value d, that is,

$\begin{array}{cc}{k}_d=w\left(d\right)k_0,& \left(12\right)\end{array}$

k₀ is the inherent thermal conductivity for undamaged material; when d=0, the heat transfer capacity of material is not affected, and when d=1, the material no longer has heat transfer capacity.

The above is a thermo-mechanical phase-field fracture model; when a mesh type method is used to solve the problem, the crack need to be divided into very fine meshes, and the calculation cost is high, so it is necessary to combine an adaptive refinement method. In the present invention, a novel adaptive strategy (the multi-level hp-FEM) is adopted, the fine mesh is refined in real time with crack propagation in the simulation process, no hanging node exists, and the numerical implementation is simple, so that the fracture behavior under thermal shock can be efficiently and effectively simulated and analyzed.

Different from the previous adaptive methods, the idea of superposition is adopted to realize adaptive refinement in the present invention, the computational domain Ω is discretized into a base mesh by a coarse mesh, and is locally discretized into a multi-level fine mesh (which is called an overlay mesh) at a crack path, as shown in FIG. 2. According to the idea of superposition, the approximate physical value at arbitrary x₀ can be expressed as follows,

$\begin{array}{cc}u=u_b+u_{o2}+\cdots +u_{{\mathrm{on}}_k}& \left(13\right)\end{array}$

Wherein the subscripts b and o represent the base mesh and the overlay mesh, respectively, and u_b, u_o2 and u_onk represent the approximate physical values on the base mesh and the overlay mesh, respectively; when such refinement strategy is used, the linear independence and compatibility of shape functions have to be considered, so as to eliminate the problem of hanging nodes in the previous adaptive methods.

In the present invention, the adaptive refinement strategy with the idea of superposition is introduced into the governing equations (1)-(6) of the thermo-mechanical phase-field fracture model, and the global computational domain Ω is discretized into the base mesh Ω_b and the multi-level overlay mesh Ω_on, wherein n=2, . . . , n_k. The temperature field θ, the displacement field u, the phase field d and spatial derivatives thereof are expressed as follows,

$\begin{array}{cc}\theta ={\theta }_b+\sum _{n=2}^{n_k}{\theta }_{\mathrm{on}},u=u_b+\sum _{n=2}^{n_k}{u}_{\mathrm{on}},d=d_b+\sum _{n=2}^{n_k}{d}_{\mathrm{on}},& \left(14\right)\end{array}$ $\begin{array}{cc}\nabla \theta =\nabla {\theta }_b+\sum _{n=2}^{n_k}\nabla {\theta }_{\mathrm{on}},\epsilon ={\epsilon }_b+\sum _{n=2}^{n_k}{\epsilon }_{\mathrm{on}},\nabla d=\nabla d_b+\sum _{n=2}^{n_k}\nabla d_{\mathrm{on}},& \left(15\right)\end{array}$

The temperature approximate values Bb and Bon in each level of mesh after discretization can be interpolated as follows,

$\begin{array}{cc}\{\begin{array}{c}{\theta }_b=N_b^N\left({\eta }_b\right){\theta }_b^N+N_b^E\left({\eta }_b\right){\theta }_b^E+N_b^F\left({\eta }_b\right){\theta }_b^F\\ {\theta }_{\mathrm{on}}=N_{\mathrm{on}}^N\left({\eta }_{\mathrm{on}}\right){\theta }_{\mathrm{on}}^N+N_{\mathrm{on}}^E\left({\eta }_{\mathrm{on}}\right){\theta }_{\mathrm{on}}^E+N_{\mathrm{on}}^F\left({\eta }_{\mathrm{on}}\right){\theta }_{\mathrm{on}}^F\end{array}& \left(16\right)\end{array}$

Wherein η. is the local coordinates under each level of mesh, N.^N, N.^E and N.^F are the shape functions attached to points, lines and surfaces, and the corresponding temperature values to be obtained on the shape functions are θ.^N, θ.^E and θ.^F. A consistent interpolation format can be used for the displacement values u_b and u_on and the phase field values d_b and d_on.

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),

$\begin{array}{cc}\{\begin{array}{c}{\int }_{\Omega }\left[\left(\mathrm{\rho c}\stackrel{.}{\theta }{\mathrm{\delta \theta }}_b-J·\nabla {\mathrm{\delta \theta }}_b\right)+\sum _{n=2}^{n_k}\left(\mathrm{\rho c}\stackrel{.}{\theta }{\mathrm{\delta \theta }}_{\mathrm{on}}-J·\nabla {\mathrm{\delta \theta }}_{\mathrm{on}}\right)\right]\mathrm{dV}=\delta F^{\theta },\\ {\int }_{\Omega }\left(\sigma :\nabla \delta u_b+\sum _{n=2}^{n_k}\sigma :\nabla \delta u_{\mathrm{on}}\right)\mathrm{dV}=\delta F^u,\\ {\int }_{\Omega }\left[-2\left(1-d\right)H\delta d_b+\frac{G_c}{l_0}\left(l_0^2\nabla d·\nabla \delta d_b\right)+d\delta d_b\right]\mathrm{dV}+\\ \sum _{n=2}^{n_k}{\int }_{\Omega }\left[-2\left(1-d\right)H\delta d_{\mathrm{on}}+\frac{G_c}{l_0}\left(l_0^2\nabla d·\nabla \delta d_{\mathrm{on}}\right)+d\delta d_{\mathrm{on}}\right]\mathrm{dV}=0,\end{array}& \left(17\right)\end{array}$

Wherein the virtual forces of the temperature field and the displacement field are as follows,

$\begin{array}{cc}\delta F^{\theta }={\int }_{\Omega }\left({\mathrm{\gamma \delta \theta }}_b+\sum _{n=2}^{n_k}{\mathrm{\gamma \delta \theta }}_{\mathrm{on}}\right)\mathrm{dV}+{\int }_{\partial \Omega }\left(\stackrel{\_}{q}·{\mathrm{\delta \theta }}_b+\sum _{n=2}^{n_k}\stackrel{\_}{q}·{\mathrm{\delta \theta }}_{\mathrm{on}}\right)\mathrm{dS},& \left(18\right)\end{array}$ $\begin{array}{cc}\delta F^{\theta }={\int }_{\Omega }\left(f·\delta u_b+\sum _{n=2}^{n_k}f·\delta u_{\mathrm{on}}\right)\mathrm{dV}+{\int }_{\partial \Omega }\left(\stackrel{\_}{t}·\delta u_b+\sum _{n=2}^{n_k}\stackrel{\_}{t}·\delta u_{\mathrm{on}}\right)\mathrm{dS}.& \left(19\right)\end{array}$

The basic format of coupled temperature-displacement-phase field in a multi-level hp mesh is derived based on the above theory and is loaded slowly in a quasi-static state, and then the Newton-Raphson alternating solution strategy is used to simulate crack initiation and propagation, so that the study of the adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock proposed by the present invention can be realized.

The present invention is realized by MATLAB software programming, and the visualization of the results is realized by ParaView software. In combination with the schematic diagram of AMLPFM as shown in FIG. 3, the specific implementation steps are as follows:

• • Step 1: establishing a multi-level discrete mesh model and a finite element interpolation format, defining material parameters, considering linear independence and compatibility of shape functions, and judging the activation states of topological components (points, lines and surfaces) in each level of mesh;
  • Step 2: initializing the phase field values by dividing the geometric model, applying the phase-field Dirichlet boundary condition weakly or weakening the material properties surrounding the pre-existing cracks according to the position of the pre-existing cracks;
  • Step 3: using the Newton-Raphson iteration method to implement the alternating solution strategy in the current time step, updating the temperature, displacement and phase fields successively, and obtaining the temperature value, the displacement value and the phase field value in the current time step when iteration converges;
  • Step 4: storing and outputting the information of relevant variables, returning to step 3, and proceeding to a next time step until the calculation is completed.

In step 1, quadrilateral elements are used for discretization and refined into multi-level meshes on the tips or paths of the pre-existing cracks, as shown in FIG. 4. For simplicity, the symbol m is used to represent the base mesh b and the overlay mesh o. Then the temperature field θ_m, the displacement field u_m and the phase field d_m of each level of mesh after discretization can be expressed as follows,

$\begin{array}{cc}\{\begin{array}{c}{\theta }_m=\sum _A{N}_{mA}^{\theta }\left(x\right){\theta }_{mA}=N_m^{\theta }{a}_m,\\ u_m=\sum _A{N}_{mA}^u\left(x\right)u_{mA}=N_m^u{\hat{a}}_m,\\ d_m=\sum _A{N}_{mA}^d\left(x\right)d_{mA}=N_m^d{\stackrel{\_}{a}}_m,\end{array}& \left(20\right)\end{array}$

Wherein N_mA^θ(x), N_mA^u(x) and N_mA^d(x) are the shape functions of the temperature, displacement and phase fields of the corresponding topological components; N_m^θ, N_m^u and N_m^d are the shape function matrices of the three fields, respectively, and the temperature value θ_mA, the displacement value u_mA and the phase field value d_mA attached to the topological components are assembled as vectors a_m, â_m and ā_m, respectively.

Then the spatial derivatives of the three fields can be further expressed as follows,

$\begin{array}{cc}\{\begin{array}{c}\nabla {\theta }_m=\sum _A{B}_{mA}^{\theta }\left(x\right){\theta }_{mA}=B_m^{\theta }{a}_m,\\ {\epsilon }_m=\sum _A{B}_{mA}^u\left(x\right)u_{mA}=B_m^u{â}_m,\\ \nabla d_m=\sum _A{B}_{mA}^d\left(x\right)d_{mA}=B_m^d{\overline{a}}_m,\end{array}& \left(21\right)\end{array}$

Wherein B_mA^θ, B_mA^u and B_mA^d are

$\begin{array}{c}\left(22\right)\end{array}$ $B_{mA}^{\theta }\left(x\right)=\left[\begin{array}{c}{\partial }_x{N}_{mA}^{\theta }\\ {\partial }_y{N}_{mA}^{\theta }\end{array}\right],B_{mA}^u\left(x\right)=\left[\begin{array}{cc}{\partial }_x{N}_{mA}^u& 0\\ 0& {\partial }_y{N}_{mA}^u\\ {\partial }_y{N}_{mA}^u& {\partial }_x{N}_{mA}^u\end{array}\right],B_{mA}^d\left(x\right)=\left[\begin{array}{c}{\partial }_x{N}_{mA}^d\\ {\partial }_y{N}_{mA}^d\end{array}\right].$

Wherein ∂_x and ∂_y represent the derivatives with respect to the x and y axes, respectively.

Considering the above discretization treatment and combining with the weak form (17) to obtain the residual equations of the three fields as follows, respectively,

$\begin{array}{cc}\{\begin{array}{c}{\left({\mathrm{\delta \theta }}_b\right)}^T·r_b^{\theta }+{\left({\mathrm{\delta \theta }}_o\right)}^T·r_o^{\theta }=0,\\ {\left(\delta u_b\right)}^T·r_b^u+{\left(\delta u_o\right)}^T·r_o^u=0,\\ {\left(\delta d_b\right)}^T·r_b^d+{\left(\delta d_o\right)}^T·r_o^d=0,\end{array}& \left(23\right)\end{array}$

Wherein r_b^θ, r_o⁰, r_b^u, r_o^u, r_b^d and r_o^d are the residuals of the temperature, displacement and phase fields in each level of mesh, and b and o are used as subscripts to replace m, then the expressions of the temperature field residual r_m^θ, the displacement field residual r_m^u and phase field residual r_m^d in each level of mesh are as follows,

$\begin{array}{cc}{r}_m^{\theta }={\left(f_m^{\theta }\right)}_{\mathrm{ext}}-{\int }_{{\Omega }_m}{\left(B_m^{\theta }\right)}^T{k}_d\nabla \theta \mathrm{dV}-{\int }_{{\Omega }_m}\rho c{\left(N_m^{\theta }\right)}^T\dot{\theta }\mathrm{dV},& \left(24\right)\end{array}$ $\begin{array}{cc}{r}_m^u={\left(f_m^u\right)}_{\mathrm{ext}}-{\int }_{{\Omega }_m}{\left(B_m^u\right)}^T\sigma \mathrm{dV},& \left(25\right)\end{array}$ $\begin{array}{cc}{r}_m^d=-\left[{\int }_{{\Omega }_m}\frac{G_c}{l_0}{\left(N_m^d\right)}^{T}d\mathrm{dV}+{\int }_{{\Omega }_m}{G}_c{l_0\left(B_m^d\right)}^T\nabla d\mathrm{dV}+{\int }_{{\Omega }_m}{w}^{\prime }\left(d\right){\left(N_m^d\right)}^{T}H\mathrm{dV}\right].& \left(26\right)\end{array}$

At the same time, the external force (f_m^θ)_ext of the temperature field and the external force (f_m^u)_ext of the displacement field are as follows,

$\begin{array}{cc}{\left(f_m^{\theta }\right)}_{ext}={\int }_{{\Omega }_m}{\left(N_m^{\theta }\right)}^T\gamma \mathrm{dV}+{\int }_{\partial {\Omega }_q}{\left(N_m^{\theta }\right)}^T\overline{q}\mathrm{dS},& \left(27\right)\end{array}$ $\begin{array}{cc}{\left(f_m^u\right)}_{ext}={\int }_{{\Omega }_m}{\rho \left(N_m^u\right)}^{T}b\mathrm{dV}+{\int }_{\partial {\Omega }_t}{\left(N_m^u\right)}^T\overline{t}\mathrm{dS}.& \left(28\right)\end{array}$

Under the framework of multi-level hp-FEM, the linear independence and compatibility of the shape functions need to be ensured. The linear independence is the requirement that the shape functions of a coarse scale mesh cannot be linearly combined by the shape functions of a fine scale mesh, and the compatibility is the requirement that the approximate solution after superposition has C⁰-continuity, which can be realized simply by deactivating (i.e. removing) the corresponding mesh topology components.

In step 2, the present invention selects a suitable method to initialize the phase field according to the position of the pre-existing cracks, as shown in FIG. 4. Two types of pre-existing cracks are found, wherein cracks 1 coincide with the edge of an element in each level of mesh, and cracks 2 pass diagonally through the inside of the element. The influence of pre-existing cracks will be considered by three methods below:

Method I: dividing the geometric model; the pre-existing cracks 1 are only refined locally in a crack tip region, and the mesh along the crack path is not treated, thus the pre-existing cracks 1 can be treated simply by dividing the geometric model; it should be noted that method I is only suitable for the first type of pre-existing cracks;

Method II: applying the phase-field Dirichlet boundary condition; the influence of arbitrary cracks can be considered by ensuring the phase field value d=1 at the pre-existing cracks in the simulated process; however, for the cracks 2 passing diagonally through the inside of the element, the Dirichlet boundary condition cannot be applied directly to the corresponding topological components to constrain the phase field value in the element; here, the penalty method is used to apply the phase-field Dirichlet boundary condition weakly to the inside of the element, and the steps are as follows:

First, selecting a number of nodes along the pre-existing cracks and letting the phase field value d=1 at the nodes, then the corresponding constraints are as follows,

$\begin{array}{cc}\left[\begin{array}{cc}{G}_{bb}& G_{bo}\\ G_{ob}& G_{oo}\end{array}\right]\left\{\begin{array}{c}{\stackrel{\_}{a}}_o\\ {\stackrel{\_}{a}}_o\end{array}\right\}=V\to G\overline{a}=V,& \left(29\right)\end{array}$

Wherein G is the shape function matrix of the selected nodes in the base mesh and the overlay mesh, and V is a vector filled by element 1;

Then considering the above constraint equations in a pure phase field diffusion model; the initial phase field distribution can be obtained, and the solution format is as follows,

$\begin{array}{cc}\left(K^d+\alpha G^{T}G\right)\overline{a}=\alpha G^{T}V,& \left(30\right)\end{array}$

Wherein α=1×10⁵ is the penalty factor, and K^d is the total stiffness matrix for phase field;

Finally, using the Newton-Raphson iteration method make V=0 in the subsequent time steps to ensure that the phase field value on the pre-existing cracks alway maintains that d=1;

Method III: weakening the material properties surrounding the pre-existing cracks; for the pre-existing cracks 2, the effect approximate to that of method I can be achieved by weakening the material properties surrounding the pre-existing cracks, which is equivalent to making the material near the crack band particularly soft; here, a minimal scalar β is introduced to weaken the elasticity matrix of the crack band, that is,

$\begin{array}{cc}{D}^{*}=\beta D,& \left(31\right)\end{array}$

Wherein 0<<β<1;

Both method II and method III can be used for treating various types of pre-existing cracks; in a solution process, method II is to introduce a penalty term to the solution format, while method III does not require additional operations.

In step 3, the present invention uses the Newton-Raphson iteration method to implement the alternating solution strategy in the current time step [t_l, t_l+1], and updates the temperature, displacement and phase fields successively; wherein {dot over (θ)} in the residual equation (24) of the temperature field is express as follows by the backward difference method,

$\begin{array}{cc}\dot{\theta }=\frac{{\theta }_{l+1}-{\theta }_l}{\Delta t},& \left(32\right)\end{array}$

Wherein Δt=t_l+1−t_l is the time increment; the coupled nonlinear equations (24)-(26) are solved by the Newton-Raphson iteration method, wherein the solution format of the i^th iteration step is as follows:

First, fixing the displacement field â^(i-1) and the phase field ā^(i-1) of the previous iteration step, and solving the temperature field a⁽ⁱ⁾ of the current iteration step according to equation (24),

$\begin{array}{cc}\left\{\begin{array}{c}{a}_b^{\left(i\right)}\\ a_o^{\left(i\right)}\end{array}\right\}=\left\{\begin{array}{c}{a}_b^{\left(i-1\right)}\\ a_o^{\left(i-1\right)}\end{array}\right\}+{\left[\begin{array}{cc}{K}_{bb}^{\theta }& K_{bo}^{\theta }\\ {\left(K_{bo}^{\theta }\right)}^T& K_{oo}^{\theta }\end{array}\right]}^{-1}\left\{\begin{array}{c}{r}_b^{\theta ,\left(i\right)}\\ r_o^{\theta ,\left(i\right)}\end{array}\right\},& \left(33\right)\end{array}$

Wherein K_bb^θ, K_oo^θ and K_bo^θ are respectively the coupled stiffness matrices for temperature of the base mesh, the overlay mesh and the two-level mesh, and the expressions are as follows,

$\begin{array}{cc}{K}_{bb}^{\theta }=-\frac{\partial r_b^{\theta ,}}{\partial a_b^{\left(i\right)}}={\int }_{{\Omega }_b}{\left(B_b^{\theta }\right)}^T{k}_d{B}_b^{\theta }\mathrm{dV}+{\int }_{{\Omega }_b}\frac{\rho {c\left(N_b^{\theta }\right)}^T{N}_b^{\theta }}{\Delta t}\mathrm{dV},& \left(34\right)\end{array}$ $\begin{array}{cc}{K}_{oo}^{\theta }=-\frac{\partial r_o^{\theta ,}}{\partial a_o^{\left(i\right)}}={\int }_{{\Omega }_o}{\left(B_o^{\theta }\right)}^T{k}_d{B}_o^{\theta }\mathrm{dV}+{\int }_{{\Omega }_o}\frac{\rho {c\left(N_o^{\theta }\right)}^T{N}_o^{\theta }}{\Delta t}\mathrm{dV},& \left(35\right)\end{array}$ $\begin{array}{cc}{K}_{bo}^{\theta }=-\frac{\partial r_b^{\theta ,}}{\partial a_o^{\left(i\right)}}={\int }_{{\Omega }_b}{\left(B_b^{\theta }\right)}^T{k}_d{B}_o^{\theta }\mathrm{dV}+{\int }_{{\Omega }_b}\frac{\rho {c\left(N_b^{\theta }\right)}^T{N}_o^{\theta }}{\Delta t}\mathrm{dV},& \left(36\right)\end{array}$

Wherein the inherent thermal conductivity k_d=k_d(ā^(i-1)) is degraded due to material damage;

Next, fixing the temperature field a⁽ⁱ⁾ of the current iteration step and the phase field ā^(i-1) of the previous iteration step, and solving the displacement field â⁽ⁱ⁾ of the current iteration step according to equation (25),

$$\begin{gathered} \begin{array}{cc}\left\{\begin{array}{c}{\hat{a}}_b^{\left(i\right)}\\ {\hat{a}}_b^{\left(i\right)}\end{array}\right\}=\left\{\begin{array}{c}{\hat{a}}_b^{\left(i-1\right)}\\ {\hat{a}}_o^{\left(i-1\right)}\end{array}\right\}+{\left[\begin{array}{cc}{K}_{bb}^u& K_{bo}^u\\ {\left(K_{bo}^u\right)}^T& K_{oo}^u\end{array}\right]}^{-1}\left\{\begin{array}{c}{r}_b^{u,\left(i\right)}\\ r_o^{u,\left(i\right)}\end{array}\right\},\text{ \\ Kb⁢bu=-∂rbu,(i)∂a^b(i)=∫Ωb(Bbu)T⁢DBbu⁢dV,}& \left(37\right)\end{array} \end{gathered}$$ $\begin{array}{cc}{K}_{bb}^u=-\frac{\partial r_b^{u,\left(i\right)}}{\partial {\hat{a}}_b^{\left(i\right)}}={\int }_{{\Omega }_b}{\left(B_b^u\right)}^T{\mathrm{DB}}_b^u\mathrm{dV},& \left(38\right)\end{array}$

Wherein the elasticity matrix

$D=\frac{\partial \sigma }{\partial \epsilon },$

and the expressions of the stiffness matrices K_oo^u and K_bo^u are similar to that of K_bb^u; it is worth noting that when method III is used to treat the pre-existing cracks, the elasticity matrix near the crack band is weakened, i.e. D*=βD;

Finally, fixing the temperature field a⁽ⁱ⁾ and the displacement field â⁽ⁱ⁾ of the current iteration step, and solving the phase field ā⁽ⁱ⁾ according to equation (26),

$\begin{array}{cc}\left\{\begin{array}{c}{\stackrel{\_}{a}}_b^{\left(i\right)}\\ {\stackrel{\_}{a}}_o^{\left(i\right)}\end{array}\right\}=\left\{\begin{array}{c}{\stackrel{\_}{a}}_b^{\left(i-1\right)}\\ {\stackrel{\_}{a}}_o^{\left(i-1\right)}\end{array}\right\}+{\left[\begin{array}{cc}{K}_{bb}^d& K_{bo}^d\\ {\left(K_{bo}^d\right)}^T& K_{oo}^d\end{array}\right]}^{-1}\left\{\begin{array}{c}{r}_b^{d,\left(i\right)}\\ r_o^{d,\left(i\right)}\end{array}\right\},& \left(39\right)\end{array}$ $\begin{array}{c}\left(40\right)\end{array}$ $K_{bb}^d={\int }_{{\Omega }_b}\left\{w^{″}\left(d\right){\left(N_b^d\right)}^T{N}_b^{d}H+\frac{G_c}{l_0}{\left(N_b^d\right)}^T{N}_b^d+G_c{l_0\left(B_b^d\right)}^T{B}_b^d\right\}\mathrm{dV},$

Wherein the expressions of the stiffness matrices K_bo^d and K_oo^d are similar to that of K_bb^d;

It should be noted that when method II is used to treat the influence of the pre-existing cracks, the constraint is Gā=0, and the constraint equation is considered into the phase field solution format (39) by the penalty method, that is

$\begin{array}{cc}{\stackrel{\_}{a}}^{\left(i\right)}={\stackrel{\_}{a}}^{\left(i-1\right)}+{\left[K^{d,\left(i\right)}+\alpha G^{T}G\right]}^{-1}{r}^{d,\left(i\right)}.& \left(41\right)\end{array}$

In combination with FIG. 3, the specific implementation process of the adaptive multi-level phase-field method (AMLPFM) for brittle fracture of elastic materials under thermal shock proposed by the present invention will be shown in the following pseudo-code form:

• • 1). Defining parameters (density ρ, Young's modulus E, Poisson's ratio ν, inherent thermal conductivity k₀, specific heat capacity c, critical energy release rate G_c), length scale l₀, time increment Δt, penalty factor α, weakening factor β and so on;
  • 2). Discretizing the finite element model of the multi-level hp mesh, selecting an appropriate method (dividing the geometric model, applying the phase-field Dirichlet boundary condition weakly or weakening the material properties) according to the position of the pre-existing cracks, and initializing the phase field value distribution;
  • 3). Judging the activation states of mesh topological components according to the discrete finite element model;
  • 4). Cycling the time steps l=0, 1, 2, . . . , N_a;
    • 4.1). Initializing the iteration variables i=0, with a⁽ⁱ⁾←a_(l), â⁽ⁱ⁾←â_(l) and ā⁽ⁱ⁾←ā_(l);
    • 4.2). Starting iterative solution;
      • 4.2.1). i←i+1;
      • 4.2.2). Fixing the displacement field â⁽ⁱ⁾ and the phase field ā^(i-1) according to equation (33), and updating the temperature field a⁽ⁱ⁾←(a_b⁽ⁱ⁾, a_o⁽ⁱ⁾;
      • 4.2.3). Fixing the temperature field a⁽ⁱ⁾ and the phase field ā^(i-1) according to equation (37), and updating the displacement field â⁽ⁱ⁾←(â_b⁽ⁱ⁾, â_o⁽ⁱ⁾;
      • 4.2.4). Fixing the temperature field a⁽ⁱ⁾ and the displacement field â⁽ⁱ⁾ according to equation (39) (if equation (41) is used to treat the pre-existing cracks, then according to equation (41)), and updating the phase field ā⁽ⁱ⁾←(ā_b⁽ⁱ⁾, ā_o⁽ⁱ⁾);
      • 4.2.5). If ∥â⁽ⁱ⁾−â^(i-1)∥/∥â⁽ⁱ⁾−â⁽⁰⁾∥<ϵ and ∥ā⁽ⁱ⁾−ā^(i-1)∥/∥ā⁽ⁱ⁾−ā⁽⁰⁾∥<ϵ, entering step 4.3); otherwise, returning to step 4.2) to continue iteration;
    • 4.3). Updating the temperature, displacement and phase fields (a_(l+1), â_(l+1), ā_(l+1))←(a⁽ⁱ⁾, â⁽ⁱ⁾, ā⁽ⁱ⁾);
    • 4.4). Outputting a calculation document of the current time step, and conducting post-treatment;
  • 5). Returning to step 4) until the calculation is completed;In the above implementation steps, ϵ is the convergence tolerance, and ϵ=1×10⁻⁵.

The present invention has the following beneficial effects:

(1) The adaptive multi-level phase-field method for brittle fracture of elastic materials 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 range of the multi-level hp-FEM in the field of fracture analysis. By using the phase field model, the shortcoming that it is difficult to accurately capture crack propagation path using a traditional continuum model can be effectively overcome, and the complex crack propagation problems such as crack intersection, bifurcation and free crack propagation in three-dimensional space can be easily treated. In addition, the present invention can further extend the constitutive model of a material to the fracture analysis of other materials, such as massive deformation fracture analysis, plastic material fracture analysis, etc.;

(2) The efficient adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock provided by the present invention uses the multi-level hp-FEM as the adaptive strategy, which can reduce the computing resources and improve the computational efficiency. Compared with the previous adaptive methods, the multi-level hp-FEM realizes local refinement through superimposed meshes, so that no hanging node exists, and the numerical implementation is simple. At the same time, the method fully utilizes the advantages of h-refinement and p-refinement, and can effectively capture high gradients, which is suitable for fracture simulation analysis. In addition, the novel adaptive strategy can also be extended to other large-scale mechanical problems, such as numerical analysis for complex structures, topological optimization, aerodynamic research, etc.;

(3) The efficient adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock provided by the present invention develops three multi-level crack treatment techniques under the multi-level framework. In method II, the pre-existing cracks are used as the boundary conditions, and the solution format is corrected by the penalty method, so as to apply the boundary conditions weakly; In method III, the material properties surrounding the pre-existing cracks are weakened to achieve an effect similar to that of method I, the essence of which is to treat the cracks as a fictitious domain. Method III can be further extended to the fictitious domain embedding technology of complex geometric bodies to conduct efficient numerical simulation for complex structures such as lattices and heat sinks;

(4) The efficient adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock provided by the present invention uses the Newton-Raphson iteration method to alternately solve the temperature, displacement and phase fields, which greatly reduces the difficulty of numerical implementation, improves the calculation efficiency, and provides a feasible scheme for the parallel computing study of large-scale complex fracture problems.

DESCRIPTION OF DRAWINGS

FIG. 1 is a schematic diagram of a thermo-mechanical fracture problem of an arbitrary solid domain Ω containing a sharp crack interface Γ and a smeared crack of the present invention, wherein (a) is a sharp crack, and (b) is a smeared crack;

FIG. 2 is a schematic diagram of a multi-level hp-FEM of the present invention, wherein (a) is a one-dimensional schematic diagram, and (b) is a two-dimensional schematic diagram;

FIG. 3 is an efficient adaptive multi-level phase-field method (AMLPFM) for brittle fracture of elastic materials under thermal shock of the present invention;

FIG. 4 is a multi-level schematic diagram of two types of pre-existing cracks of the present invention;

FIG. 5 shows distribution of normal and weakened integration points in pre-existing crack treatment method III of the present invention;

FIG. 6 is a schematic diagram of geometry and boundary conditions for the double edges notched tension test in embodiment 1 of the present invention;

FIG. 7 is a phase field evolution cloud map of pre-existing cracks treated by method II and method III in embodiment 1 of the present invention, wherein (a) is method II and (b) is method III;

FIG. 8 is a comparison diagram of load-displacement curves of pre-existing cracks treated by method II and method III in embodiment 1 of the present invention;

FIG. 9 shows a phase field evolution cloud map and an experimental result diagram of ceramic plate quenching fracture under thermal shock in embodiment 2 of the present invention, wherein (a) is t=0.5 ms, (b) is t=5 ms, (c) is t=30 ms, (d) is t=75 ms, (e) is t=200 ms, and (f) an experimental result.

DETAILED DESCRIPTION

Specific embodiments of the present invention will be further described below in combination with the drawings and the technical solution.

EMBODIMENTS

The accuracy, reliability and excellent performance of the efficient adaptive multi-level phase-field method (AMLPFM) for brittle fracture of elastic materials under thermal shock proposed by the present invention are further described in detail in combination with FIGS. 6-9.

First, quasi-static tensile brittle fracture simulation of a square plate with inclined initial cracks will be performed to demonstrate the reliability of the developed pre-existing crack treatment techniques and the effectiveness and high efficiency of the adaptive phase-field method for solving fracture problems; then ceramic plate quenching thermal cracking is simulated, the crack path is compared with the experimental results to verify the accuracy of AMLPFM in simulating the fracture behavior under thermal shock. The reference solution in the first embodiment is derived from the work of Kim et al., and the experimental results in the second embodiment are derived from the work of Jiang et al. In all embodiments, quadrilateral elements are used, and the phase field model parameter l₀ is set to be greater than or equal to half of the mesh size h_e/2 to satisfy the rationality of the numerical simulation results of phase-field fracture. All embodiments are considered as plane strain cases, and the gravity effect is not considered.

1) Embodiment 1: Double Edges Notched Tension Test (FIGS. 6-8)

In embodiment 1, the double edges notched tension test as shown in FIG. 6 is simulated, the meshes are refined into 3 levels at the cracks, and the pre-existing cracks are treated by method II (weakly applied phase field boundary conditions) and method III (weakened material properties at the crack), respectively. The lower boundary of the geometric body is fixed, and a vertical displacement load is applied to the upper boundary. The material parameters are shown in Table 1 below.

TABLE 1
Material parameters for double edges notched tension test
	Parameter	E(GPa)	v	l0(mm)	Gc(kN/mm)
	Value	210	0.3	0.015	0.0027

FIG. 7 shows a phase field evolution cloud map of the pre-existing cracks treated by method II and method III. In method II, the phase field value at the pre-existing cracks is ensured to be d=1 before the displacement load is applied, so that the material at the cracks does not have any bearing capacity; in method III, the phase field value is not initialized, and the phase field value at the pre-existing cracks changes gradually from 0 to 1 with the application of the displacement load. Consistent crack propagation paths are obtained by both of the two methods, and crack tips are tracked in the simulation process to dynamically updated the meshes, thus to further reduce the number of degree of freedoms, which makes AMLPFM has a great potential to solve large-scale complex fracture problems. FIG. 8 shows the load-displacement curves of the two methods, and takes the work of Kim et al. as the reference solution. It can be seen that the three curves basically coincide, so that the reliability of the results thereof is verified. In summary, embodiment 1 proves the effectiveness and high efficiency of the adaptive phase-field method proposed by the present invention and the accuracy of the three pre-existing crack treatment techniques.

(2) Embodiment 2: Simulation of Ceramic Plate Quenching Fracture Under Thermal Shock (FIG. 9)

In embodiment 2, ceramic plate quenching fracture under thermal shock is simulated and compared with the experimental results. After a ceramic plate is heated to 300° C. and thrown into a pool of 20° C., the sharp cooling of the ceramic plate leads to crack initiation and propagation. After the ceramic plate is fished out and colored with ink, an internal crack propagation pattern can be observed, as shown in FIG. 9(f). Considering the symmetry of the geometric body, only a quarter model is used here for numerical simulation, and the material parameters thereof are shown in Table 2 below.

TABLE 2
Material parameters of ceramic
plate quenching fracture under thermal shock
			Gc	ρ	k0
Para-	E		(kN/	(kg/	(W/	c	α
meter	(GPa)	ν	mm)	m3)	m · K)	(J/(kg · K)	(—)
Value	370	0.3	42.47	3980	31	880	7.5 × 10−6

FIG. 9 shows an evolution cloud map of ceramic plate quenching fracture and the comparison with the experimental results. It can be seen that the meshes are refined dynamically and adaptively with crack propagation, and numerous cracks are produced in the ceramic plate, which are similar to the experimental results, thus it is verified that AMLPFM proposed in the present invention can carry out an efficient and effective numerical simulation study on the fracture behavior of an elastic-brittle material under thermal shock. As shown in FIG. 9(a), at the beginning of the simulation, a consistent damage band is produced on the outer boundary of the ceramic plate, and no crack initiates. Then in FIG. 9(b), small cracks {circle around (1)} of equal length initiate rapidly at equal spacing on the boundary. Subsequently, as shown in FIG. 9(c), about half of the small cracks {circle around (1)} continue to propagate inward to cracks {circle around (2)}, and the other half of the cracks {circle around (1)} stop propagating due to insufficient supply of strain energy. The phenomenon occurs again in FIG. 9(d), and the simulation stops after the small cracks propagate to the cracks {circle around (3)}. The final crack paths are shown in FIG. 9(e), and the results are similar to the experimental results. The fracture process is extremely short, so that it is difficult to observe the crack evolution process through experiments, but the crack propagation process under thermal shock can be more easily and deeply explored by the numerical simulation, which further demonstrates the practicability of AMLPFM of the present invention and has certain scientific and engineering practical significance.

In summary, the above two embodiments verify in detail the accuracy and effectiveness of AMLPFM proposed by the present invention from different levels, respectively, and show the significant advantages of the method in terms of computational efficiency and computational scale, proving the necessity of the present invention. At the same time, the multi-level hp-FEM is used as an adaptive strategy, and the meshes are dynamically refined with crack propagation. Compared with other adaptive strategies, the present invention has no hanging nodes, but has simple numerical implementation, and can also be used in other numerical simulations with high gradients. Therefore, the efficient adaptive multi-level phase-field method (AMLPFM) for brittle fracture of elastic materials under thermal shock proposed by the present invention is a high performance numerical method with a great development prospect.

Embodiments of the present invention are given for example and description purposes, but are not exhaustive or used to limit the present invention to the disclosed forms. Many modifications and changes are apparent to those skilled in the art. The purpose of selecting and describing the embodiments is to preferably illustrate the principles and practical applications of the present invention, so that those skilled in the art can understand the present invention, thereby designing various modified embodiments applied to specific uses.

Claims

1. An adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock, comprising the following steps: considering an arbitrary thermo-mechanical solid domain (computational domain) Ω, containing a sharp crack Γ in time t∈[0, ta]; combining with the existing thermo-mechanical phase-field fracture model, the following governing equations is obtained,

2. The adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock according to claim 1, wherein the adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock is realized by MATLAB software programming, and the visualization of the results is realized by ParaView software; the specific steps are as follows: step 1: establishing a multi-level discrete mesh model and a finite element interpolation format, defining material parameters, considering linear independence and compatibility of shape functions, and judging the activation states of topological components (i.e. points, lines and surfaces) in each level of mesh;step 2: initializing the phase field values by dividing the geometric model, applying the phase-field Dirichlet boundary condition weakly or weakening the material properties surrounding the pre-existing cracks according to the position of the pre-existing cracks;step 3: using the Newton-Raphson iteration method to implement the alternating solution strategy in the current time step, updating the temperature, displacement and phase fields successively, and obtaining the temperature value, the displacement value and the phase field value in the current time step when iteration converges;step 4: storing and outputting the information of relevant variables, returning to step 3, and proceeding to a next time step until the calculation is completed.

3. The adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock according to claim 2, wherein, in step 1, quadrilateral elements are used for discretization and refined into multi-level meshes on the tips or paths of the pre-existing cracks; the symbol m is used to represent the base mesh b and the overlay mesh o; then the temperature field θm, the displacement field um and the phase field dm of each level of mesh after discretization are expressed as follows,

4. The adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock according to claim 2, wherein in step 2, the phase field is initialized according to the position of the pre-existing cracks; two types of pre-existing cracks are found, wherein cracks 1 coincide with the edge of an element in each level of mesh, and cracks 2 pass diagonally through the inside of the element; the influence of pre-existing cracks will be considered by three methods: method I: dividing the geometric model; the pre-existing cracks 1 are only refined locally in a crack tip region, and the mesh along the crack path is not treated, thus the pre-existing cracks 1 is treated by dividing the geometric model; method I is only suitable for the first type of pre-existing cracks;method II: applying the phase-field Dirichlet boundary condition; the influence of arbitrary cracks can be considered by ensuring the phase field value d=1 at the pre-existing cracks in the simulated process; however, for the pre-existing cracks 2 passing diagonally through the inside of the element, the Dirichlet boundary condition cannot be applied directly to the corresponding topological components to constrain the phase field value in the element; the penalty method is used to apply the phase-field Dirichlet boundary condition weakly to the inside of the element, and the steps are as follows:first, selecting a number of nodes along the pre-existing cracks and letting the phase field value d=1 at the nodes, then the corresponding constraints are as follows,

5. The adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock according to claim 2, wherein in step 3, the Newton-Raphson iteration method is used to implement the alternating solution strategy in the current time step [tl, tl+1], and the temperature, displacement and phase fields are updated successively; wherein {dot over (θ)} in the residual equation (24) of the temperature field is express as follows by the backward difference method,

6. The adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock according to claim 2, wherein the specific implementation process of the adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock will be shown in the following pseudo-code form: 1). defining parameters (density ρ, Young's modulus E, Poisson's ratio ν, inherent thermal conductivity k0, specific heat capacity c, critical energy release rate Gc), length scale l0, time increment Δt, penalty factor α, weakening factor β and so on;2). discretizing the finite element model of the multi-level hp mesh, selecting the method of dividing the geometric model, applying the phase-field Dirichlet boundary condition weakly or weakening the material properties according to the position of the pre-existing cracks, and initializing the phase field value distribution;3). judging the activation states of mesh topological components according to the discrete finite element model;4). cycling the time steps l=0, 1, 2, . . . , Na; 4.1). initializing the iteration variables i=0, with a(i)←a(l), â(i)←â(l) and ā(i)←ā(l);4.2). starting iterative solution; 4.2.1). i←i+1;4.2.2). fixing the displacement field â(i) and the phase field ā(i-1) according to equation (33), and updating the temperature field a(i)←(ab(i), ao(i);4.2.3). fixing the temperature field a(i) and the phase field ā(i-1) according to equation (37), and updating the displacement field â(i)←(âb(i), âo(i);4.2.4). fixing the temperature field a(i) and the displacement field â(i) according to equation (39) (if equation (41) is used to treat the pre-existing cracks, then according to equation (41)), and updating the phase field ā(i)←(āb(i), āo(i));4.2.5). if ∥â(i)−â(i-1)∥/∥â(i)−â(0)∥<ϵ and ∥ā(i)−ā(i-1)∥/∥ā(i)−ā(0)∥<ϵ, entering step 4.3); otherwise, returning to step 4.2) to continue iteration;4.3). updating the temperature, displacement and phase fields (a(l+1), â(l+1), ā(l+1))←(a(i), â(i), ā(i));4.4). outputting a calculation document of the current time step, and conducting post-treatment;5). returning to step 4) until the calculation is completed;

7. The adaptive multi-level phase-field method for brittle fracture of elastic materials under thermal shock according to claim 6, wherein, ϵ is the convergence tolerance, and ϵ=1×10−5.