Self-adaptive multilevel grid phase field method for analyzing elastic-brittle fracture under thermal shock

Abstract

The invention belongs to the field of computational mechanics, provides a self-adaptive multi-stage grid phase field method for elastic-brittle fracture analysis under thermal shock, and provides a brand-new numerical calculation method for efficient fracture failure research. According to the method, the damage degree of the material is represented by adopting the phase field, the multi-stage hp finite element serves as the self-adaptive strategy, dynamic discretization processing is carried out in the crack expansion process under thermal shock, and compared with a previous self-adaptive phase field fracture model, the method has the advantages that the suspension point problem does not exist, numerical value implementation is simple, the calculation efficiency is high, and the method is suitable for large-scale popularization and application. And the fracture failure behavior of the structure can be effectively predicted. Meanwhile, three methods of pre-crack processing technologies are developed to initialize a phase field, a temperature field, a displacement field and the phase field are subjected to staggered iteration solution, and the multi-field coupling large-scale fracture damage problem can be stably and efficiently solved. In addition, the method can be further expanded to defect prediction research in the additive manufacturing process.

Description

Technical Field

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

[0002] Thermal shock damage is widely present in engineering practice and daily life, such as thermal cracking caused by rapid temperature changes in manufacturing processes (additive manufacturing, quenching and welding, etc.), microcracks caused by different expansion coefficients in multiphase composite materials during processing, and self-explosion of window glass in hot 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 and costly to implement, but also difficult to consider the influence of multiple external factors on fracture failure behavior at the same time, and it is difficult to observe the crack evolution process. As an emerging high-performance simulation computing technology, numerical simulation methods can help us to explore the failure process under thermal shock in more depth, so as to comprehensively and accurately predict the failure behavior of materials. However, the study of multi-field coupled 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-physics field coupling solution and effective integration of efficient adaptive methods.

[0003] For the numerical simulation of fracture problems, there are currently two main types of methods: discrete crack characterization methods and diffuse crack characterization methods. Among them, the discrete crack characterization method describes cracks by explicitly defining crack set information, and there are mesh processing methods such as unit deletion method, cohesive zone model (CZM) and extended finite element method (XFEM). During the simulation process, cracks will expand, and their geometric shape changes with load size and time, which requires continuous mesh re-division during the solution process, which greatly consumes computing 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 diffuse crack characterization method characterizes cracks through field variables, and does not require separate modeling or mesh re-division of cracks. To a certain extent, it improves the computational efficiency and reduces the difficulty of numerical implementation, and has significant advantages in dealing with complex fracture problems.

[0004] In the characterization method of diffuse cracks, Bourdin et al. developed a regularized approximate diffuse fracture model, namely the phase field model, based on the fracture variation 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 the material, and has been widely used in dealing with crack initiation and extension problems. Amor et al. decomposed the strain energy into volume-deviatoric strain in "Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments", assuming that only deviatoric stress and positive spherical stress will drive crack evolution. In 2010, Miehe et al. reformulated the phase field fracture model in the framework of thermodynamics in "A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits", decomposed the strain energy into tension and compression through the strain tensor spectrum, and believed that only tension deformation would drive crack evolution, thereby preventing unrealistic crack extension during compression. In order to simulate the fracture process more accurately, Borden et al. proposed a fourth-order phase field fracture equation under the framework of isogeometric analysis method to simulate the fracture of brittle materials in "A phase-field description of dynamic brittle fracture". In order to more accurately deal with the quasi-brittle fracture of materials, Professor Wu Jianying introduced the concept of cohesion into the phase field fracture model, proposed the phase field regularized cohesive crack model (PF-CZM), and combined the BFGS quasi-Newton iteration algorithm to further improve the computational efficiency. The phase field model does not need to update the geometric topology of the model to simulate crack evolution. It has advantages in dealing with complex fracture modes such as crack initiation, extension and intersection. It is widely used in various types of fracture failure problems, such as thermal cracking, dynamic fracture, hyperelastic fracture, hydraulic fracturing, impact failure, etc. Among them, the study of thermal cracking is helpful for the safety assessment of manufacturing processes such as additive manufacturing, quenching, and casting, and has attracted the attention of many scholars. Tangella et al. used a hybrid phase field model in "Hybrid phase-field modeling of thermo-elastic crack propagation" to simulate the fracture of elastic-brittle materials.Ruan et al. proposed a phase-field fracture model under thermal-mechanical coupling in "A thermo-mechanical phase-field fracture model: Application to hot cracking simulations in additive manufacturing", discussed the effect of temperature on fracture mode, and applied it to the simulation of micro-crack initiation in additive manufacturing. Mandal et al. used the PF-CZM model to simulate the hot cracking phenomenon of elastic-brittle materials in "Fracture of thermo-elastic solids: Phase-fieldmodeling and new results with an efficient monolithic solver", and used the BFGS iteration method to further improve the calculation efficiency.

[0005] Phase field models have been widely used in simulating fracture processes, but the numerical simulation process requires extremely fine meshes, which makes the calculation process large. In order to meet the mesh requirements without significantly increasing the computational cost, adaptive mesh methods are widely used in phase field fracture simulation. Heister et al. proposed a predictor-corrector mesh adaptive 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-adaptive strategy that locally subdivides the mesh in real time during crack propagation. Patil et al. proposed a phase field fracture model under adaptive multiscale finite element (MsFEM), which can significantly reduce the computational cost and improve the computational efficiency. Recently, Zander et al. proposed a new adaptive method in "Multi-level hp-adaptivity: high-order mesh adaptivity without the difficulties of constraining hanging nodes", called the multi-level hp finite element method (The multi-level hp-FEM), which uses superimposed fine meshes to achieve local subdivision. In most previous adaptive methods, fine meshes are used to replace coarse meshes, which will cause hanging point problems and make numerical implementation complicated. The multi-level hp finite element method uses the superposition idea and the integral Legendre basis function to construct the shape function, so that there is no hanging point problem, and the numerical implementation is simpler. At the same time, it fully utilizes the advantages of h-subdivision and p-subdivision, allowing high convergence in dealing with problems with high gradient changes. These characteristics are suitable for phase field fracture simulation, where multi-level fine meshes are locally covered at the crack, thereby reducing the consumption of computing resources.

[0006] In the present invention, a multi-level hp finite element method is used as an adaptive strategy to dynamically discretize the phase field fracture model under thermal-mechanical coupling, and three prefabricated crack processing techniques are developed to accurately apply 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 present invention proposes an efficient numerical method for analyzing the fracture failure problem of elastic-brittle materials under thermal shock, which can efficiently capture the crack evolution process, does not have the hanging point problem in the traditional adaptive method, and is simple to implement numerically. It is used to simulate the elastic-brittle fracture mechanical behavior of the structure when it is subjected to thermal shock load, and provides an effective way for the accurate analysis of engineering fracture problems. Summary of the invention

[0007] Aiming at the efficient numerical analysis of the fracture and destruction of elastic-brittle materials under thermal shock, the present invention innovatively proposes an adaptive multi-level grid phase field fracture simulation method (English name: Adaptive multi-level phase-fieldmodel, English abbreviation: AMLPFM), the purpose of which is: First, in order to overcome the deficiency of traditional damage models that it is difficult to accurately capture the crack propagation path, and to explore the fracture behavior under thermal shock, the present invention derives a phase field fracture model for the fracture and destruction behavior of elastic-brittle materials under thermal-mechanical coupling based on the energy functional variational principle, and adopts a staggered iterative algorithm to solve the coupling control equation. Secondly, in view of the fine grid requirements of the fracture phase field and the problem of hanging points in traditional adaptive strategies, the present invention proposes a novel adaptive multi-level grid discretization scheme that can be more easily implemented numerically. In the simulation process, based on the advancement of the phase field contour, the crack tip is tracked to identify the refined area, and the grid is dynamically updated in real time to reduce the computational cost and improve the computational efficiency. In addition, the present invention proposes three prefabricated crack processing techniques (splitting geometric model, weakly applying Dirichlet boundary conditions, and weakening material properties at the crack) to accurately consider the influence of various prefabricated cracks under the framework of multi-level grids. Finally, the present invention aims to solve the shortcomings of the existing finite element-phase field method in analyzing the fracture failure problem under thermal shock, which is caused by the high computational cost due to the small grid, and the limitation of the complex numerical implementation caused by the hanging point problem in the previous adaptive method.

[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 crack simulation method (AMLPFM). In order to accurately capture the high gradient changes of thermal stress during thermal shock and reduce the computational cost in fracture simulation analysis, the present invention combines a new adaptive method based on multi-layer grids (The multi-level hp-method) with a phase field fracture model under thermal-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 crack simulation method, the specific steps are as follows:

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

[0011]

[0012] Where, equations (1), (2) and (3) are respectively the displacement field u = Ω × [0, t a], temperature field θ=Ω×[0,t a ] and phase field d=Ω×[0,t a ]→[0,1], 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 at which temperature changes 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 diffuse crack zone, b is the body force, J is the heat flow, and n is the solid domain boundary The unit external normal vector of is; the boundary conditions of displacement, surface force, temperature and heat flow are and The displacement boundary is applied accordingly Surface force temperature With heat flow

[0013] Considering the strain spectrum decomposition format given by Miehe et al., the strain energy density ψ e Decompose into stretch parts With compression part And the stretched part is weakened, that is,

[0014]

[0015] Among them, ε e is the elastic strain, w(d)=(1-d) 2 +ξ is a degenerate function, 0<ξ<<1 is added for better numerical stability.

[0016] Then, the Cauchy stress σ can be derived as,

[0017]

[0018] In the thermomechanical coupling framework, the elastic strain ε e and thermal strain ε θ There is the following relationship between them:

[0019] ε e =ε-ε θ , (9)

[0020] Where, under the small strain assumption, the total strain ε θ =α(θ-θ0)I is the thermal strain expression, α is the thermal expansion coefficient, θ0 is the initial temperature, and I is the second-order unit tensor.

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

[0022]

[0023] Thus preventing Non-physical crack healing due to decline.

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

[0025]

[0026] In order to ensure that no heat flow passes through the crack surface, the thermal conductivity coefficient k d is also degenerated by the phase field value d,

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

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

[0029] The above is a phase field fracture model under thermal-mechanical coupling. When a mesh-based method is used to solve the problem, the crack needs to be divided into extremely fine meshes, which has a high computational cost. It is necessary to combine an adaptive subdivision method. In the present invention, a novel adaptive strategy (multi-level hp finite element method, The multi-level hp-FEM) is adopted. During the simulation process, the fine mesh is subdivided in real time as the crack expands, and there is no hanging point problem. The numerical implementation is simple, thereby carrying out efficient and effective numerical simulation analysis of the fracture failure behavior under thermal shock.

[0030] Different from the previous adaptive methods, the present invention adopts the superposition idea to realize 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, called the covering grid, as shown in Figure 2 As shown. According to the superposition concept, the approximate physical value at any x0 can be expressed as,

[0031]

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

[0033] In the present invention, the above-mentioned adaptive subdivision strategy with superposition idea is introduced into the control equations (1)-(6) of the phase field fracture model under thermal-mechanical coupling, and the overall calculation domain Ω is discretized into the base grid Ωb With 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 value θ in each layer of the grid after discretization b ,θ on The interpolation value is,

[0036]

[0037] Among them, η · is the local coordinate of each layer of grid, N · N 、N · E With N · F is the shape function attached to the point, line and surface, and the corresponding temperature value to be determined on these shape functions is θ · N ,θ · E and θ · F For displacement value u b 、u on and the phase field value d b ,d on A consistent interpolation format can be used.

[0038] Consider any dummy variable δθ b ,δu b ,δd b ,δθ on ,δu on and δd on , the following weak form can be obtained from the control equations (1)-(6):

[0039]

[0040] Among them, the virtual forces of temperature field and displacement field are:

[0041]

[0042] According to the above theory, the basic format of coupled temperature-displacement-phase field under multi-level hp grid is derived. Through quasi-static slow loading and the Newton-Raphson staggered solution strategy, the crack initiation and propagation process is simulated, and the adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock proposed in the present invention can be realized.

[0043] The present invention is implemented by programming with MATLAB software, and the result visualization is realized by using ParaView software. Figure 3 The AMLPFM method calculation cycle diagram is shown, and its specific implementation steps are as follows:

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

[0045] Step 2: Initialize the phase field value according to the location of the prefabricated crack using a split geometry model, a weakly imposed phase field Dirichlet boundary condition, or weakened crack material properties;

[0046] Step 3: In the current time step, the Newton-Raphson iteration method is used to implement the staggered solution strategy, and the temperature field, displacement field and phase field are updated in sequence. When the iteration converges, the temperature value, displacement value and phase field value in the current time step are obtained;

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

[0048] In step 1, quadrilateral elements are used for discretization, and the prefabricated crack tip or path is subdivided into multiple layers of grids, such as Figure 4 For simplicity, the symbol m is used to represent the base grid b and the cover grid o. Then, the temperature field θ of each layer of grid after discretization is m , displacement field u m With phase field d m It can be expressed as,

[0049]

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

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

[0052]

[0053] in, and for

[0054]

[0055] in, and Represent the partial derivatives with respect to the x-axis and y-axis respectively.

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

[0057]

[0058] in, and is the residual of temperature, displacement and phase field in each layer of grid. The subscripts b and o are replaced by m. Then the residual of temperature field in each layer of grid is Displacement field residual and phase field residual The expression is,

[0059]

[0060] At the same time, the temperature field force External force with displacement field for,

[0061]

[0062] In the framework of multi-level hp finite element, linear independence and compatibility between shape functions must be guaranteed. Linear independence requires that the shape functions of coarse-scale meshes cannot be linearly combined with the shape functions of fine-scale meshes. Compatibility requires that the approximate solution after superposition has C 0 Continuity can be simply achieved in numerical implementation by deactivating (i.e. removing) the corresponding mesh topology components.

[0063] In step 2, the present invention selects a suitable method to initialize the phase field according to the position of the prefabricated crack, such as Figure 4 As shown in the figure, there are two kinds of prefabricated cracks. Crack 1 coincides with the unit edge of each level of mesh, and crack 2 passes through the unit obliquely. Below, the influence of prefabricated cracks will be considered through three methods:

[0064] Method I: Split the geometric model. The prefabricated crack 1 is only locally subdivided in the crack tip area, and the mesh on the crack path is not processed otherwise. The 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 prefabricated crack.

[0065] Method II: Applying phase field Dirichlet boundary conditions. By ensuring that the phase field value d=1 at the prefabricated crack in the simulation process, the influence of any crack can be considered. However, for crack 2 that passes obliquely through the unit, the Dirichlet boundary conditions cannot be directly applied to the corresponding topological components to constrain the phase field value inside the unit. Here, the penalty function method is used to weakly apply the phase field Dirichlet boundary conditions to the unit. 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 constraints are as follows:

[0067]

[0068] Where G is the shape function matrix of the selected point in the base grid and the cover grid, and V is a vector filled with element 1;

[0069] Then, taking the above constraint equation into account in 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, K d is the phase field global stiffness matrix;

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

[0073] Method III: Weakening material properties. For prefabricated crack 2, the material properties at the prefabricated crack can also be weakened, which is equivalent to making the material near the crack zone particularly soft, and the effect of method I can be achieved. Here, a minimum scalar β is introduced to weaken the elastic matrix of the crack zone, that is,

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

[0075] Among them, 0<<β<1;

[0076] Both Method II and Method III can handle various types of 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 performs the following steps at the current time step [t l ,t l+1 ] adopts the Newton-Raphson iteration method to implement the staggered solution strategy, and updates the temperature field, displacement field and phase field in turn; among them, the temperature field residual equation (24) Using the backward difference method, it can be expressed 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, where the solution format of the i-th iteration step is as follows:

[0080] First, fix the displacement field of the previous iteration step Phase Field According to equation (24), solve the temperature field a of the current iteration step (i) ,

[0081]

[0082] in, and are the temperature stiffness matrices of the base grid, cover grid and two-level grid coupling, respectively, and the expression is:

[0083]

[0084] Among them, the thermal conductivity Degradation due to material damage;

[0085] Next, fix the temperature field a of the current iteration step (i) and the phase field of the previous iteration According to equation (25), the displacement field of the current iteration step is solved

[0086]

[0087] Among them, the elasticity matrix Stiffness Matrix and The expression of It is worth noting that when method III is used to treat the prefabricated crack, the elastic matrix near the crack zone is weakened, that is, D * =βD;

[0088] Finally, fix the temperature field a of the current iteration step (i) With displacement field According to equation (26), the phase field is solved

[0089]

[0090] Among them, the phase field stiffness matrix and The expression of The expression is similar to;

[0091] It should be noted that when using method II to deal with the influence of prefabricated cracks, there are constraints The constraint equation is taken into account in the phase field solution format (39) by the penalty function method, that is,

[0092]

[0093] Combination 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 in the present invention will be presented in the form of the following pseudo code:

[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 the appropriate method (splitting geometry model, weakly applied phase field Dirichlet boundary conditions or weakened material properties) according to the prefabricated crack location to initialize the phase field value distribution;

[0096] 3) According to the discrete finite element model, determine the activation state of the grid topology component;

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

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

[0099] 4.2) Iterative solution starts;

[0100] 4.2.1) i←i+1;

[0101] 4.2.2) According to equation (33), the fixed displacement field and phase field Update temperature field

[0102] 4.2.3) According to equation (37), the fixed temperature field a (i) and phase field Update 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)), the fixed temperature field a (i) and displacement field Update Phase Field

[0104] 4.2.5) If and Go to step 4.3), otherwise, return to step 4.2) and continue iterating;

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

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

[0107] 5) Return to step 4) until the calculation is completed;

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

[0109] Beneficial effects 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 element in the field of fracture analysis. The phase field model can effectively overcome the deficiency of the traditional continuous medium model in accurately capturing the crack propagation path, and can relatively simply handle complex crack propagation problems, such as crack crossing, bifurcation and free crack propagation in three-dimensional space, etc., and the present invention can further extend the material constitutive law to the fracture failure analysis of other materials, such as large deformation fracture analysis, plastic material fracture analysis, etc.;

[0111] (2) The present invention provides an efficient adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock. In this method, multi-level hp finite element is used as an adaptive strategy, which can reduce computing resources and improve computing efficiency. Compared with previous adaptive methods, multi-level hp finite element achieves local subdivision by superimposing grids, and there is no hanging point problem, which makes the numerical implementation simple. At the same time, this method makes full use of the advantages of h subdivision and p subdivision, 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, topological optimization, and aerodynamics research;

[0112] (3) The present invention provides an efficient adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock. Three multi-level crack processing technologies are developed under the multi-level grid framework. Method II uses prefabricated cracks as boundary conditions and modifies the solution format through the penalty function method, thereby weakly applying boundary conditions; Method III weakens the material properties at the crack to achieve an effect similar to that of Method I. Its essence is to regard the crack as a virtual domain, which can be further extended to the virtual domain embedding technology of complex geometric bodies, and carry out efficient numerical simulation analysis for complex structures such as lattices and heat sinks;

[0113] (4) The efficient adaptive multi-level grid phase field method for brittle fracture analysis under thermal shock provided by the present invention adopts the Newton-Raphson iteration method to perform an alternating solution of the temperature field, displacement field and phase field, which greatly reduces the difficulty of numerical implementation and improves the calculation efficiency, providing a feasible solution for the parallel calculation research of large-scale complex fracture failure problems. BRIEF DESCRIPTION OF THE DRAWINGS

[0114] Figure 1 Schematic diagram of the fracture problem under thermal-mechanical coupling of an arbitrary solid domain Ω containing a sharp crack interface Γ and a diffuse crack according to the present invention, wherein (a) is a sharp crack and (b) is a diffuse crack;

[0115] Figure 2 Schematic diagram of the multi-level hp finite element method of the present invention, wherein (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 a high-efficiency 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 multi-level grids at two prefabricated cracks of the present invention;

[0118] Figure 5 is the Gaussian integral point distribution in the prefabricated crack processing method III of the present invention;

[0119] Figure 6 Schematic diagram of the structure and boundary conditions of Example 1 of the quasi-static tensile fracture of a prefabricated rectangular plate with double oblique cracks according to the present invention;

[0120] Figure 7 The crack phase field evolution cloud diagram of the prefabricated cracks treated by method II and method III in Example 1 of the present invention, wherein (a) is method II and (b) is method III;

[0121] Figure 8 This is a comparison diagram of force-displacement curves of treating prefabricated cracks using method II and method III in Example 1 of the present invention;

[0122] Fig. 9 These are the phase field evolution cloud diagrams and experimental result diagrams of the quenching fracture of the ceramic plate under thermal shock in Example 2 of the present invention, where (a) is t=0.5ms, (b) is t=5ms, (c) is t=30ms, (d) is t=75ms, (e) is t=200ms, and (f) is the experimental result. DETAILED DESCRIPTION

[0123] The specific implementation of the present invention is further described below in conjunction with the accompanying drawings and technical solutions.

[0124] Example

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

[0126] First, a quasi-static tensile brittle fracture simulation of a square plate with an oblique initial crack will be implemented to illustrate the reliability of the developed prefabricated crack processing technology, as well as the effectiveness and efficiency of the adaptive grid-phase field method proposed in the present invention in solving fracture problems; then, the quenching thermal cracking of the ceramic plate is simulated, and the crack path is compared with the experimental results to verify the accuracy of AMLPFM in 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 the grid size h e / 2 is used to satisfy the rationality of the phase field fracture numerical simulation results. All embodiments are considered as plane strain assumptions, and the gravity effect is not considered.

[0127] 1) Example 1: Quasi-static tensile fracture of a prefabricated rectangular plate with double oblique cracks (see Appendix Figure 6-8 )

[0128] Example 1 simulates the following Figure 6 The tensile fracture of a rectangular plate with double oblique cracks is shown. The mesh is subdivided into three levels at the crack, and the prefabricated cracks are treated using method II (weakly applied phase field boundary conditions) and method III (weakened material properties at the 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.

[0129] Table 1. Quasi-static tensile brittle fracture material parameters of prefabricated double-oblique crack rectangular plate

[0130]

[0131] Figure 7 The phase field evolution cloud diagrams of method II and method III for processing prefabricated cracks are given. Method II ensures that the phase field value d=1 at the prefabricated crack before displacement loading, so that the material at the crack has no bearing capacity; method III does not initialize the phase field value, and its phase field value at the prefabricated crack gradually changes from 0 to 1 with displacement loading. Both methods obtain consistent crack propagation paths, and dynamically update the mesh by tracking the crack tip during the simulation process, further reducing the total number of degrees of freedom, making the AMLPFM method have great potential to solve large-scale complex fracture problems. Figure 8 The force-displacement curves of the two methods are given in the figure, 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. In summary, Example 1 proves the effectiveness and efficiency of the adaptive grid-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 plate under thermal shock (see Appendix) Fig. 9 )

[0133] Example 3 simulates the quenching fracture phenomenon of ceramic plates and compares it with the experimental results. When the ceramic plate is heated to 300°C and thrown into a 20°C water pool, the ceramic plate cools down rapidly, causing cracks to initiate and expand. After being fished out and stained with ink, the internal crack expansion pattern can be observed. The results are as follows: Fig. 9 Considering the symmetry of the geometry, only one quarter of the model is used for numerical simulation, and its material parameters are shown in Table 2.

[0134] Table 2. Material parameters of ceramic plate quenching fracture under thermal shock

[0135]

[0136] Fig. 9 The evolution cloud diagram of the quenched fracture of the ceramic plate is given and compared with the experimental results. It can be seen that the grid will be dynamically and adaptively subdivided as the crack expands, and many cracks are generated inside the ceramic plate, which is close to the experimental results. This verifies that the AMLPFM proposed in this invention can carry out efficient and effective numerical simulation research on the fracture and failure behavior of elastic-brittle materials under thermal shock. Fig. 9 As shown in (a), at the beginning of the simulation, a consistent damage zone is generated on the outer boundary of the ceramic plate, and no cracks initiate. Then, in 9(b), tiny cracks of equal length and distance quickly initiate on the boundary. Then, as shown in Fig. 9 As shown in (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. Fig. 9(d) and then the simulation stops after the crack ③ is extended. The final crack path is shown in Fig. 9 (e) shows that the result is close to the experimental result. The fracture process is extremely short, and it is difficult to observe the crack evolution process through experiments. However, numerical simulation can more easily explore the crack propagation process under thermal shock, which further illustrates 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 respectively verify the accuracy and effectiveness of the AMLPFM method proposed in the present invention in detail from different levels, and demonstrate the significant advantages of the method in terms of computational efficiency and computational scale, proving the necessity of its invention. At the same time, a multi-level hp finite element is used as an adaptive strategy, and its grid is dynamically subdivided as the crack propagates. Compared with other adaptive strategies, there is no hanging point problem, and the numerical implementation is simple, and it can also be used in other numerical simulations of high gradient changes. Therefore, the efficient adaptive multi-level grid phase field method (AMLPFM) for elastic-brittle fracture analysis under thermal shock proposed in the present invention is a high-performance numerical method with great development prospects.

[0138] The embodiments of the present invention are given for the purpose of illustration and description, and are not intended to be exhaustive or to limit the invention to the disclosed forms. Many modifications and variations will be apparent to those of ordinary skill in the art. The embodiments are selected 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 thereby design various embodiments with various modifications suitable for specific uses.

Claims

1. An adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock, characterized in that: Here are the steps: Consider an arbitrary solid domain under thermal-mechanical coupling, namely the computational domain Ω, at time t∈[0,t a ] contains a sharp crack Γ; combined with the existing phase field fracture model under thermal-mechanical coupling, the following control equation is obtained: Wherein, equations (1), (2) and (3) are respectively the displacement field u = Ω × [0, t a ], temperature field θ=Ω×[0,t a ] and 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 at which temperature changes 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 diffuse crack zone, b is the body force, J is the heat flow, and n is the solid domain boundary The unit external normal vector of is; the boundary conditions of displacement, surface force, temperature and heat flow are and The displacement boundary is applied accordingly Surface force temperature With heat flow According to the given strain spectrum decomposition format, the strain energy density ψ e Decompose into stretch parts With compression part And the stretched part is weakened, that is, Among them, ε e is the elastic strain, w(d)=(1-d) 2 +ξ is a degenerate function, 0<ξ<<1 is added to obtain better numerical stability; Then, the Cauchy stress σ can be derived as, In the thermomechanical coupling framework, the elastic strain ε e and thermal strain ε θ There is the following relationship between them: e e =e-e θ (9) Where, under the small strain assumption, the total strain ε θ =α(θ-θ0)I is the thermal strain expression, α is the thermal expansion coefficient, θ0 is the initial temperature, and I is the second-order unit tensor; In order to satisfy the irreversibility of crack growth, a historical maximum strain energy H is introduced to replace the Thus preventing Non-physical crack healing due to descent; According to Fourier's law, the heat flow J is, In order to ensure that no heat flow passes through the crack surface, the thermal conductivity coefficient k d is also degenerated 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 capacity of the material is not affected, while when d = 1, the material no longer has the heat transfer capacity. A novel adaptive strategy is used to solve the above-mentioned phase field fracture model under thermal-mechanical coupling, as follows: The superposition idea is used to realize adaptive subdivision. 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, called covering grids. According to the superposition idea, the approximate physical value at any x0 is expressed as, The subscripts b and o represent the base grid and the overlay grid, respectively. b 、u o2 and Represent the approximate physical values ​​on the base grid and the overlay grid respectively; The above-mentioned adaptive subdivision strategy with superposition idea is introduced into the control equations (1)-(6) of the phase field fracture model under thermal-mechanical coupling, and the overall computational domain Ω is discretized into the base grid Ω b With 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, The approximate temperature value θ in each layer of the grid after discretization b ,θ on The interpolation value is, Among them, η is the local coordinate of each layer of grid, N . N 、N . E With N . F is the shape function attached to the point, line and surface, and the corresponding temperature value to be determined on these shape functions is θ . N ,θ . E and θ F ; For displacement value u b 、u on and the phase field value d b ,d on Use a consistent interpolation format; Consider any dummy variable δθ b ,δu b ,δd b ,δθ on ,δu on and δd on , the following weak form is obtained from the governing equations (1)-(6), Among them, the virtual forces of temperature field and displacement field are: Based on the above theoretical deduction, the basic format of coupled temperature-displacement-phase field under multi-level hp grid is obtained. Through quasi-static slow loading and the Newton-Raphson staggered solution strategy, the crack initiation and propagation process is simulated, thus realizing 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 is characterized in that: The adaptive multi-level grid phase field of elastic-brittle fracture analysis under thermal shock is realized by MATLAB software programming, and the results are visualized by ParaView software. The specific steps are as follows: Step 1: Establish a multi-level discrete grid model and finite element interpolation format, define material parameters, consider the linear independence and compatibility of shape functions, and determine the activation status of topological components, i.e. points, lines, and surfaces, in each layer of the grid; Step 2: Initialize the phase field value according to the location of the prefabricated crack using a split geometry model, a weakly imposed phase field Dirichlet boundary condition, or weakened crack material properties; Step 3: In the current time step, the Newton-Raphson iteration method is used to implement the staggered solution strategy, and the temperature field, displacement field and phase field are updated in sequence. When the iteration converges, the temperature value, displacement value and phase field value in the current time step are obtained; Step 4: Store and output relevant variable information, return to step 3, and enter the next time step 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 is characterized in that: In step 1, quadrilateral elements are used for discretization, and the prefabricated crack tip or path is subdivided into multiple layers of grids; the symbol m represents the base grid b and the cover grid o; then, the temperature field θ of each layer of grid after discretization is m , displacement field u m With phase field d m It is expressed as, in, and are the shape functions of temperature, displacement and phase field of the corresponding topological components, and The shape function matrices of the three fields are respectively, and the temperature value θ attached to the topological component mA , displacement value u mA and the phase field value d mA are assembled into vector a m , and Then the spatial derivatives of the three fields are further expressed as, in, and for in, and Respectively represent the partial derivatives of the x-axis and y-axis; Considering the above discretization and combining the weak form (17), the residual equations of the three fields are obtained as follows: in, and is the residual of temperature, displacement and phase field in each layer of grid. The subscripts b and o are replaced by m. Then the residual of temperature field in each layer of grid is Displacement field residual and phase field residual The expression is, At the same time, the temperature field force External force with displacement field for, 4. 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 2, the phase field is initialized according to the position of the prefabricated cracks. There are two types of prefabricated cracks. Prefabricated crack 1 coincides with the unit edge of each level of grid, and prefabricated crack 2 obliquely passes through the unit. The influence of prefabricated cracks will be considered in three ways: Method I: Split the geometric model. The prefabricated crack 1 is only locally subdivided in the crack tip area, and no other processing is done on the mesh on the crack path. The prefabricated crack 1 is processed by splitting the geometric model. Method I is only applicable to the first prefabricated crack. Method II: Applying phase field Dirichlet boundary conditions, by ensuring that the phase field value d=1 at the prefabricated crack in the simulation process, the influence of any crack can be considered; however, for the prefabricated crack 2 that obliquely passes through the unit, the Dirichlet boundary conditions cannot be directly applied to the corresponding topological components to constrain the phase field value inside the unit; the penalty function method is used to weakly apply the phase field Dirichlet boundary conditions inside the unit, 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 constraints are as follows: Where G is the shape function matrix of the selected point in the base grid and the cover grid, and V is a vector filled with element 1; Then, the above constraint equation is taken into account in the pure phase field diffusion model to obtain the initial phase field distribution, and its solution format is: Where α = 1 × 10 5 is the penalty factor, K d is the phase field global stiffness matrix; Finally, when using the Newton-Raphson iteration method in the subsequent time steps, V = 0 to ensure that the phase field value on the prefabricated crack always maintains d = 1; Method III: Weakening material properties. For prefabricated crack 2, the material properties at the prefabricated crack can also be weakened, which is equivalent to making the material near the crack zone particularly soft, to achieve the effect of method I. Here, a minimum scalar β is introduced to weaken the elastic matrix of the crack zone, that is, D * =βD (31) Among them, 0<<β<1; Both Method II and Method III can handle various types of 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, at the current time step [t l ,t l+1 ] adopts the Newton-Raphson iteration method to implement the staggered solution strategy, and updates the temperature field, displacement field and phase field in turn; among them, the temperature field residual equation (24) Using the backward difference method, it can be expressed as: 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, where the solution format of the i-th iteration step is as follows: First, fix the displacement field of the previous iteration step Phase Field According to equation (24), solve the temperature field a of the current iteration step (i) , in, and are the temperature stiffness matrices of the base grid, cover grid and two-level grid coupling, respectively, and the expression is: Among them, the thermal conductivity Degradation due to material damage; Next, fix the temperature field a of the current iteration step (i) and the phase field of the previous iteration According to equation (25), the displacement field of the current iteration step is solved Among them, the elasticity matrix Stiffness Matrix and The expression of The expression of is similar to that of; when method III is used to treat the prefabricated crack, the elastic matrix near the crack zone is weakened, that is, D * =βD; Finally, fix the temperature field a of the current iteration step (i) With displacement field According to equation (26), the phase field is solved Among them, the phase field stiffness matrix and The expression of The expression is similar to; It should be noted that when using method II to deal with the influence of prefabricated cracks, there are constraints The constraint equation is taken into account in the phase field solution format (39) by the penalty function method, that is, 6. The adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock according to claim 2, characterized in that: The specific implementation process of the adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock will be presented in the following pseudo code form: 1) Definition parameters include density ρ, elastic modulus E, Poisson's ratio v, thermal conductivity k0, specific heat capacity c, critical fracture energy density G c , crack characteristic width l0, time increment Δt, penalty factor α, weakening coefficient β; 2) Discretize the finite element model of the multi-level hp grid, and according to the prefabricated crack location, select the splitting geometry model, weakly apply the phase field Dirichlet boundary condition or weaken the material properties, and initialize the phase field value distribution; 3) According to the discrete finite element model, determine the activation state of the grid topology component; 4) Time step cycle l = 0, 1, 2, ..., N a ; 4.1) Initialize the iteration variable i = 0, and a (i) ←a (l) , and 4.2) Iterative solution starts; 4.2.1) i←i+1; 4.2.2) According to equation (33), the fixed displacement field and phase field Update temperature field 4.2.3) According to equation (37), the fixed temperature field a (i) and phase field Update displacement field 4.2.4) According to equation (39) (if equation (41) is used to process the prefabricated crack, then according to equation (41)), the fixed temperature field a (i) and displacement field Update Phase Field 4.2.5) If and Go to step 4.3), otherwise, return to step 4.2) and continue iterating; 4.3) Update temperature, displacement and phase field 4.4) Output the calculation file of the current time step for post-processing; 5) Return to step 4) until the calculation is completed.

7. The adaptive multi-level grid phase field method for elastic-brittle fracture analysis under thermal shock according to claim 6, characterized in that: Where ∈ is the convergence tolerance, and ∈=1×10 -5 .

Images