Elastoplastic Fracture Analysis of the P91 Steel Welded Joint under Repair Welding Thermal Shock Based on XFEM

: P91 steel is a typical steel used in the manufacture of boilers in ultra-supercritical power plants and heat exchangers in nuclear power plants. For the long-term serviced P91 steel pressurized structures, the main failure mode is the welded joint failure, especially the heat a ﬀ ected zone (HAZ) failure. Repair welding technique is an e ﬀ ective method for repairing such local defects. However, the thermal shock composed of high temperature and thermal stress in the repair welding process will pose a critical loading condition for the existing defects near the heat source which cannot be detected by conventional means. So, the evaluation of structural integrity for the welded joint in the thermal-mechanical coupling ﬁeld is necessary. In this work, the crack propagation law in the HAZ for the P91 steel welded joint was investigated under repair welding thermal loads. The weld repair model of the P91 steel welded joint was established by ABAQUS. The transient temperature ﬁeld and stress ﬁeld in repair welding process were calculated by relevant user subroutines and sequential coupling simulation method. The residual stress was determined by the impact indentation strain method to verify the feasibility of the ﬁnite element (FE) model and simulation method. In order to obtain the crack propagation path, the elastoplastic fracture analysis of the welded joint with initial crack was performed based on the extended ﬁnite element method (XFEM). The inﬂuence of di ﬀ erent welding linear energy on the crack propagation was analyzed. The results show that the cracks in the HAZ propagate perpendicular to the surface and tend to deﬂect to the welding seam under repair welding thermal loads. The crack propagation occurs in the early stage of cooling. Higher welding linear energy leads to larger HAZ and higher overall temperature. With the increase of welding linear energy, the length and critical distance of the crack propagation increase. Therefore, low welding linear energy can e ﬀ ectively inhibit the crack propagation in the HAZ. The above calculation and analysis provide a reference for the thermal shock damage analysis of repair welding process, which is of great signiﬁcance to improving the safety and reliability of weld repaired components.


Introduction
P91 steel is widely used in boiler components of ultra-supercritical power plants due to its excellent creep strength and steam corrosion resistance at high temperature [1]. It is also one of the common choices of steel pipe material for heat exchangers in nuclear power plants [2]. In practical application, P91 steel is usually connected by welding to form the pressure equipment components. During long-term service, local defects will occur inevitably in pressure equipment due to the effect of harsh service environment or improper operation procedure, especially in vulnerable parts, such as Repair welding is a rapid welding heat transfer process, and the arc temperature during welding process may reach thousands of degrees centigrade [8]. The welding heat source loads a large amount of heat flow into the welding seam and nearby area in a short time. The high temperature and the unsteady thermal stress act on the weldment in the form of thermal shock. On the one hand, it will cause the formation of hot cracks in the welding seam; on the other hand, the interaction between existing cracks and thermal-mechanical coupling field will lead to the crack propagation and even the failure of the whole component. In order to ameliorate the quality of weldments and weld repaired components, it is necessary to study the effect of the welding process on the cracking behavior. Alvarez et al. [9] compared the hot cracking sensitivity of tungsten inert gas welding (TIG) and laser beam welding (LBW) by analyzing the microstructure and chemical composition of the welding seam of 718 alloy. Chelladurai et al. [10] studied the solidification cracking behavior of 316 stainless steel under different energy transfer modes by adjusting the pulse parameters in pulsed laser welding (PLW). Hosseini et al. [11] investigated the effect of heat input and welding speed in electron beam welding (EBW) on the hot cracking sensitivity of AA2024-T351 alloy experimentally. Coniglio et al. [12] explored the initiation mechanism of solidification crack in the arc welding seam for 6060 aluminum alloy theoretically and numerically. Wei et al. [13] developed a software which can automatically plot the driving force and resistance curves of solidification cracks according to the numerical and experimental results. Bordreuil et al. [14] established a solidification cracking model which can predict microstructure and pore nucleation in the welding seam by combining cellular automata model with intergranular fluid flow model. Agarwal et al. [15] analyzed the influence of metallurgical factors and molten pool shape on the solidification cracking in laser beam welding of advanced high strength steels numerically and experimentally. Jiang et al. [16] researched the effect of repair width on residual stress for the composite plate to reduce the probability of interface cracking. Hyde et al. [17] studied the creep crack propagation behavior of P91 steel weldment under constant loads by both experimental method and finite element (FE) method. Pandey et al. [18] investigated experimentally the effect of different diffusible hydrogen concentrations in weld metal Repair welding is a rapid welding heat transfer process, and the arc temperature during welding process may reach thousands of degrees centigrade [8]. The welding heat source loads a large amount of heat flow into the welding seam and nearby area in a short time. The high temperature and the unsteady thermal stress act on the weldment in the form of thermal shock. On the one hand, it will cause the formation of hot cracks in the welding seam; on the other hand, the interaction between existing cracks and thermal-mechanical coupling field will lead to the crack propagation and even the failure of the whole component. In order to ameliorate the quality of weldments and weld repaired components, it is necessary to study the effect of the welding process on the cracking behavior. Alvarez et al. [9] compared the hot cracking sensitivity of tungsten inert gas welding (TIG) and laser beam welding (LBW) by analyzing the microstructure and chemical composition of the welding seam of 718 alloy. Chelladurai et al. [10] studied the solidification cracking behavior of 316 stainless steel under different energy transfer modes by adjusting the pulse parameters in pulsed laser welding (PLW). Hosseini et al. [11] investigated the effect of heat input and welding speed in electron beam welding (EBW) on the hot cracking sensitivity of AA2024-T351 alloy experimentally. Coniglio et al. [12] explored the initiation mechanism of solidification crack in the arc welding seam for 6060 aluminum alloy theoretically and numerically. Wei et al. [13] developed a software which can automatically plot the driving force and resistance curves of solidification cracks according to the numerical and experimental results. Bordreuil et al. [14] established a solidification cracking model which can predict microstructure and pore nucleation in the welding seam by combining cellular automata model with intergranular fluid flow model. Agarwal et al. [15] analyzed the influence of metallurgical factors and molten pool shape on the solidification cracking in laser beam welding of advanced high strength steels numerically and experimentally. Jiang et al. [16] researched the effect of repair width on residual stress for the composite plate to reduce the probability of interface cracking. Hyde et al. [17] studied the creep crack propagation behavior of P91 steel weldment under constant loads by both experimental method and finite element (FE) method. Pandey et al. [18] investigated experimentally the effect of different diffusible hydrogen concentrations in weld metal on hydrogen-induced cracking features of P91 steel. Zhang et al. [19] applied the high energy spark deposition (HESD) method to weld repair and tested the fracture properties of the repair welding seam. He et al. [20] performed the elastic fracture analysis of cracked aluminum alloy plate during metal inert gas welding (MIG) and cooling processes.
According to the existing literature, a lot of theoretical and experimental studies have been carried out on the hot cracks in welding processes. Almost all the current numerical simulations of cracking behavior for welding focus on the structure without crack or with stationary crack. The thermal effect of welding, especially repair welding, will make the existing cracks which are neglected by conventional detection methods further extend and even penetrate through the whole structure. However, there is little research on the crack propagation behavior under repair welding thermal shock. Therefore, it is necessary to analyze the structural integrity of weldments under repair welding thermal loads.
Combined with the specific damage model, the extended finite element method (XFEM) can simulate the ductile crack growth behavior accurately [21]. Belytschko [22] first proposed the embryonic form of the XEFM. It was independent of mesh generation and enhanced the finite element approximation by adding discontinuous enrichment functions to the displacement field near the crack tip. Moës et al. [23] introduced the step function and crack tip asymptotic function to describe the crack surface and crack tip, respectively, making it successfully used in the analysis of fracture mechanics, and called this technique "extended finite element method". Daux et al. [24] added several different asymptotic functions and step functions to the crack tip and crack surface to simulate the crack branching. Chessa et al. [25] improved the convergence of the blending element in the XFEM (shown in Figure 2) by utilizing the extended strain method. Stolarska et al. [26] used the level set method to locate the crack location, which further improved the theory of the XFEM. For completing the modelling of dynamic crack, Song et al. [27] applied the phantom node method to the XFEM, which implemented the definition of the crack propagation behavior within the framework of the conventional finite element method. Until now, FE software such as ABAQUS, ANSYS and LS-DYNA have implemented the XFEM program into the function module of fracture analysis. This paper uses the XFEM to predict how the repair welding thermal shock affects the cracking behavior of the heat affected zone (HAZ) in a P91 steel welded joint. The influence of welding linear energy on the cracking features has been studied, which provides a reference for repair welding of structures containing defects. This study has guiding significance for the life extension of pressure equipment under long-term service or extended service.
Metals 2020, 10, x FOR PEER REVIEW 3 of 26 on hydrogen-induced cracking features of P91 steel. Zhang et al. [19] applied the high energy spark deposition (HESD) method to weld repair and tested the fracture properties of the repair welding seam. He et al. [20] performed the elastic fracture analysis of cracked aluminum alloy plate during metal inert gas welding (MIG) and cooling processes. According to the existing literature, a lot of theoretical and experimental studies have been carried out on the hot cracks in welding processes. Almost all the current numerical simulations of cracking behavior for welding focus on the structure without crack or with stationary crack. The thermal effect of welding, especially repair welding, will make the existing cracks which are neglected by conventional detection methods further extend and even penetrate through the whole structure. However, there is little research on the crack propagation behavior under repair welding thermal shock. Therefore, it is necessary to analyze the structural integrity of weldments under repair welding thermal loads.
Combined with the specific damage model, the extended finite element method (XFEM) can simulate the ductile crack growth behavior accurately [21]. Belytschko [22] first proposed the embryonic form of the XEFM. It was independent of mesh generation and enhanced the finite element approximation by adding discontinuous enrichment functions to the displacement field near the crack tip. Moës et al. [23] introduced the step function and crack tip asymptotic function to describe the crack surface and crack tip, respectively, making it successfully used in the analysis of fracture mechanics, and called this technique "extended finite element method". Daux et al. [24] added several different asymptotic functions and step functions to the crack tip and crack surface to simulate the crack branching. Chessa et al. [25] improved the convergence of the blending element in the XFEM (shown in Figure 2) by utilizing the extended strain method. Stolarska et al. [26] used the level set method to locate the crack location, which further improved the theory of the XFEM. For completing the modelling of dynamic crack, Song et al. [27] applied the phantom node method to the XFEM, which implemented the definition of the crack propagation behavior within the framework of the conventional finite element method. Until now, FE software such as ABAQUS, ANSYS and LS-DYNA have implemented the XFEM program into the function module of fracture analysis. This paper uses the XFEM to predict how the repair welding thermal shock affects the cracking behavior of the heat affected zone (HAZ) in a P91 steel welded joint. The influence of welding linear energy on the cracking features has been studied, which provides a reference for repair welding of structures containing defects. This study has guiding significance for the life extension of pressure equipment under long-term service or extended service.

Model Geometry
Two P91 steel plates with a size of 200 mm × 100 mm × 12 mm (length, width and thickness, respectively) are butt welded with four welding layers, including two layers of root welding (welding layers 1 to 2), one layer of filler welding (welding layer 3) and one layer of cover welding (welding layer 4). The selection of the form and size of welding groove is based on the criteria of Yang et al. [28]. A V-shaped groove is adopted, the root gap and blunt edge are both 2 mm, the groove angle is

Model Geometry
Two P91 steel plates with a size of 200 mm × 100 mm × 12 mm (length, width and thickness, respectively) are butt welded with four welding layers, including two layers of root welding (welding layers 1 to 2), one layer of filler welding (welding layer 3) and one layer of cover welding (welding layer 4). The selection of the form and size of welding groove is based on the criteria of Yang et al. [28]. A V-shaped groove is adopted, the root gap and blunt edge are both 2 mm, the groove angle is 60 • . The area containing a type IV crack in the HAZ will be removed and then refilled by repair welding. The size of weld repair area is 80 mm × 8 mm × 3 mm (length, width and thickness, respectively), and the angle of groove surface for repair welding is 30 • . Figure 3 shows the schematic drawing of the model geometric dimensions.
Metals 2020, 10, x FOR PEER REVIEW 4 of 26 60°. The area containing a type IV crack in the HAZ will be removed and then refilled by repair welding. The size of weld repair area is 80 mm × 8 mm × 3 mm (length, width and thickness, respectively), and the angle of groove surface for repair welding is 30°. Figure 3 shows the schematic drawing of the model geometric dimensions. The finite element modelling tool of thermal-mechanical coupling process and fracture behavior is ABAQUS 6.14-5 (Dassault Systemes, Paris, France). In order to save the calculating work and avoid the convergence problem in 3-D computational fracture mechanics, a 2-D plane strain analysis model is set up. The plane studied in this paper is the middle cross section of the plate, as shown in Figure  4. Mesh refinement is carried out in the welding seam and HAZ, where exist high heat flux loading, to overcome the mesh-sensitivity problem. The area far away from the welding seam and HAZ is coarsened appropriately to reduce the calculation time. The transitional mesh is adopted between the refined mesh and the coarsened mesh. The FE model and meshing are shown in Figure 5. A total of 9063 elements and 9322 nodes are meshed. The element edge size for the refined mesh is 0.25 mm, and that for the coarsened mesh is 2 mm.  The finite element modelling tool of thermal-mechanical coupling process and fracture behavior is ABAQUS 6.14-5 (Dassault Systemes, Paris, France). In order to save the calculating work and avoid the convergence problem in 3-D computational fracture mechanics, a 2-D plane strain analysis model is set up. The plane studied in this paper is the middle cross section of the plate, as shown in Figure 4. Mesh refinement is carried out in the welding seam and HAZ, where exist high heat flux loading, to overcome the mesh-sensitivity problem. The area far away from the welding seam and HAZ is coarsened appropriately to reduce the calculation time. The transitional mesh is adopted between the refined mesh and the coarsened mesh. The FE model and meshing are shown in Figure 5. A total of 9063 elements and 9322 nodes are meshed. The element edge size for the refined mesh is 0.25 mm, and that for the coarsened mesh is 2 mm.
Metals 2020, 10, x FOR PEER REVIEW 4 of 26 60°. The area containing a type IV crack in the HAZ will be removed and then refilled by repair welding. The size of weld repair area is 80 mm × 8 mm × 3 mm (length, width and thickness, respectively), and the angle of groove surface for repair welding is 30°. Figure 3 shows the schematic drawing of the model geometric dimensions. The finite element modelling tool of thermal-mechanical coupling process and fracture behavior is ABAQUS 6.14-5 (Dassault Systemes, Paris, France). In order to save the calculating work and avoid the convergence problem in 3-D computational fracture mechanics, a 2-D plane strain analysis model is set up. The plane studied in this paper is the middle cross section of the plate, as shown in Figure  4. Mesh refinement is carried out in the welding seam and HAZ, where exist high heat flux loading, to overcome the mesh-sensitivity problem. The area far away from the welding seam and HAZ is coarsened appropriately to reduce the calculation time. The transitional mesh is adopted between the refined mesh and the coarsened mesh. The FE model and meshing are shown in Figure 5. A total of 9063 elements and 9322 nodes are meshed. The element edge size for the refined mesh is 0.25 mm, and that for the coarsened mesh is 2 mm.  The sequential coupling method is used in the thermal-mechanical coupling simulation. Firstly, the thermal analysis is performed, and then the temperature field results are applied to the nodes of the mechanical model to carry out the mechanical analysis. The removal and deposition of the metal are implemented by the method of Birth and Death Element. A 4-node linear heat transfer quadrilateral (DC2D4) element is adopted for the thermal analysis, and a 4-node bilinear plane strain quadrilateral reduced-integration (CPE4R) element is adopted for the mechanical analysis. The topological relationship for meshing between the thermal model and the mechanical model is consistent. The sequential coupling method is used in the thermal-mechanical coupling simulation. Firstly, the thermal analysis is performed, and then the temperature field results are applied to the nodes of the mechanical model to carry out the mechanical analysis. The removal and deposition of the metal are implemented by the method of Birth and Death Element. A 4-node linear heat transfer quadrilateral (DC2D4) element is adopted for the thermal analysis, and a 4-node bilinear plane strain quadrilateral reduced-integration (CPE4R) element is adopted for the mechanical analysis. The topological relationship for meshing between the thermal model and the mechanical model is consistent.

Material Property
Temperature-dependent nonlinear features are considered for the thermal physical properties and mechanical properties of P91 steel, as shown in Figure 6. The thermal physical properties include thermal conductivity, specific heat and thermal expansion coefficient, while the mechanical properties include Young's modulus, Poisson's ratio and yield strength. The density of P91 steel is assumed to be constant at a value of 7780 kg/m 3 . The material properties of the weld metal are also considered in the FE simulation for obtaining higher accuracy. The solidus and liquidus temperatures of both base metal and weld metal are set at 1420 °C and 1500 °C, respectively. The melting heat and solidification heat of the material are both set at 260 kJ/kg.

Material Property
Temperature-dependent nonlinear features are considered for the thermal physical properties and mechanical properties of P91 steel, as shown in Figure 6. The thermal physical properties include thermal conductivity, specific heat and thermal expansion coefficient, while the mechanical properties include Young's modulus, Poisson's ratio and yield strength. The density of P91 steel is assumed to be constant at a value of 7780 kg/m 3 . The material properties of the weld metal are also considered in the FE simulation for obtaining higher accuracy. The solidus and liquidus temperatures of both base metal and weld metal are set at 1420 • C and 1500 • C, respectively. The melting heat and solidification heat of the material are both set at 260 kJ/kg. The sequential coupling method is used in the thermal-mechanical coupling simulation. Firstly, the thermal analysis is performed, and then the temperature field results are applied to the nodes of the mechanical model to carry out the mechanical analysis. The removal and deposition of the metal are implemented by the method of Birth and Death Element. A 4-node linear heat transfer quadrilateral (DC2D4) element is adopted for the thermal analysis, and a 4-node bilinear plane strain quadrilateral reduced-integration (CPE4R) element is adopted for the mechanical analysis. The topological relationship for meshing between the thermal model and the mechanical model is consistent.

Material Property
Temperature-dependent nonlinear features are considered for the thermal physical properties and mechanical properties of P91 steel, as shown in Figure 6. The thermal physical properties include thermal conductivity, specific heat and thermal expansion coefficient, while the mechanical properties include Young's modulus, Poisson's ratio and yield strength. The density of P91 steel is assumed to be constant at a value of 7780 kg/m 3 . The material properties of the weld metal are also considered in the FE simulation for obtaining higher accuracy. The solidus and liquidus temperatures of both base metal and weld metal are set at 1420 °C and 1500 °C, respectively. The melting heat and solidification heat of the material are both set at 260 kJ/kg.

Heat Source Definition
The heat source loads should reflect the actual temperature change during welding as accurately as possible. The distributed heat flux (DFLUX) has been widely used because of its mature theory. Among the distributed heat flux, Gaussian heat flux distribution, uniform body heat flux distribution and double ellipsoidal heat flux distribution are commonly used [30]. The thermal exchange during welding is mainly the thermal conduction inside the weldment, which follows Fourier's law (Equation (1)), and the governing equation of temperature field (Equation (2)) follows the law of conservation of energy. where: R is the heat flux (W/m 2 ); λ is the thermal conductivity (W/m·k); T is the distribution function of temperature field (K); ρ is the density (kg/m 3 ); c is the specific heat (J/kg·K); t is the transient time (s); Q is the intensity of thermal energy (W), including the thermal energy generated by the heat source and the thermal energy generated by the solid-liquid phase change.
The conventional uniform body heat flux distribution is defined by Equation (3) [29].
where q is the heat generation rate by heat source; η is the arc efficiency factor; U is the arc voltage; I is the welding current; V act is the action volume of the heat source. The conventional uniform body heat flux distribution adopts amplitude curves to control the piecewise changes of the heat generation rate with time. The forms of the changes are linear in each time segment, which cannot reflect the real situation of the heat source passing through the plane. In addition, the value of V act needs to be estimated first and then calibrated according to the experimental results. Continuous adjustment must be made until the simulated results are consistent with the experimental results.
An improved uniform body heat flux distribution is used to simulate the welding process. The action volume of the heat source is calculated according to the groove size and welding process parameters. As shown in Figure 7, assuming that all thermal energy is concentrated in the filler metal, V act is determined by the product of the cross-sectional area of the welding layer (A) and the action length of the heat source (C). C is estimated by Equation (4) [31].
where τ is the coefficient determined by the welding method and welding process parameters, and in this paper τ = 3 mm·kW −1 . Therefore, the total action time of the heat source is obtained by Equation (5). Therefore, the total action time of the heat source is obtained by Equation (5).
where v is the welding speed.
Considering the effect of the heat source approaching and leaving the studied cross section during the welding process, the improved uniform body heat flux distribution changes with time instantaneously. The distributed heat flux is calculated by Equation (6) [32].
where t is the transient loading time; t 0 is the time taken for the center of the heat source to move to the studied cross section, equal to half of t act . The modelling of the improved uniform body heat flux distribution is done by FORTRAN language written in DFLUX, one of user subroutines in ABAQUS. The arc efficiency factor is assumed to be 0.8 for the shielded metal arc welding (SMAW), and 0.6 for the gas tungsten arc welding (GTAW).

Boundary Conditions
The ambient temperature is assumed as 20 • C. The convection heat transfer exists between air and weldment surface during welding due to the temperature difference. Additionally, the temperature difference between weldment surface and surrounding environment will result in continuous radiation heat transfer. The convection heat transfer follows Newton's law of cooling (Equation (7)).
where q c is the heat flux of convection heat transfer; α c is the convection heat transfer coefficient; ∆T is the value of temperature difference between air and weldment surface. The convection heat transfer coefficient is calculated by Equation (8), which is introduced into the FE model by the user subroutine FILM in ABAQUS 6.14-5.
where T is the transient temperature of the weldment surface. The heat flux radiated outward by weldment follows the Stefan-Boltzmann law (Equation (9)).
where ε is the emissivity, assumed as 0.85; σ is the Stefan-Boltzmann constant with the value of 5.67 × 10 −8 W/m 2 ·K 4 . T is the temperature of the weldment surface. The radiation heat transfer between weldment surface and surrounding environment is calculated by Equations (10) and (11).
where q r is the heat flux of radiation heat transfer; α r is the radiation heat transfer coefficient; T f is the ambient temperature.
During the mechanical analysis, degrees of freedom for node A and node B at two ends of the bottom surface of the FE model (shown in Figure 5) are constrained in the X direction and Y direction to prevent the rigid body displacement.

Thermal Loading Patterns
Welding linear energy is the heat energy input by welding heat source to unit length welding seam, which influences the surface forming of welding seam and the formation of welding defects. The magnitude of welding linear energy is calculated by Equation (12).
Three different repair welding linear energy levels will be considered in the fracture analysis. The values of three selected linear energies are 10 kJ/cm, 16 kJ/cm and 25 kJ/cm, respectively. According to the heat flux distribution defined by Equation (6), the time distribution curves of heat generation rate in the weld repair area under three different linear energy rates is plotted in Figure 8. The heat generation rate presents a Gaussian profile with time, which can reflect the moving features of the heat source. With the increase of linear energy, both the duration of the thermal loading and the generated total heat energy increase. The duration of the thermal loading is 3.96 s for the linear energy of 10 kJ/cm, 6.34 s for the linear energy of 16 kJ/cm and 9.90 s for the linear energy of 25 kJ/cm. All three thermal loading patterns have the same peak heat generation rate in the intermediate instant of the corresponding loading time, which is 14.8 GW/m 3 . The intermediate instant is the time when the center of the heat source passes through the studied cross section. 3 4 273 273 100 100 where qr is the heat flux of radiation heat transfer; αr is the radiation heat transfer coefficient; Tf is the ambient temperature. During the mechanical analysis, degrees of freedom for node A and node B at two ends of the bottom surface of the FE model (shown in Figure 5) are constrained in the X direction and Y direction to prevent the rigid body displacement.

Thermal Loading Patterns
Welding linear energy is the heat energy input by welding heat source to unit length welding seam, which influences the surface forming of welding seam and the formation of welding defects. The magnitude of welding linear energy is calculated by Equation (12).
Three different repair welding linear energy levels will be considered in the fracture analysis. The values of three selected linear energies are 10 kJ/cm, 16 kJ/cm and 25 kJ/cm, respectively. According to the heat flux distribution defined by Equation (6), the time distribution curves of heat generation rate in the weld repair area under three different linear energy rates is plotted in Figure 8.

Damage Model
The XFEM used in this paper is based on the cohesive crack model [33]. The damage response of a crack takes the traction-separation law as the constitutive relation. As shown in Figure 9, the

Damage Model
The XFEM used in this paper is based on the cohesive crack model [33]. The damage response of a crack takes the traction-separation law as the constitutive relation. As shown in Figure 9, the damage development includes three stages: damage initiation, damage evolution and element failure.  The damage initiation of a crack is the starting point at which the cohesive stiffness between the crack surfaces begins to degrade. It occurs when the stress or strain of the element reaches a critical value. Here, the maximum principal stress (MPS) criterion is used to judge whether a crack reaches the damage initiation condition. Meanwhile, the crack propagation direction is also controlled by the The damage initiation of a crack is the starting point at which the cohesive stiffness between the crack surfaces begins to degrade. It occurs when the stress or strain of the element reaches a critical value. Here, the maximum principal stress (MPS) criterion is used to judge whether a crack reaches the damage initiation condition. Meanwhile, the crack propagation direction is also controlled by the MPS criterion. The cohesive crack will propagate along the direction of maximum principal stress. The damage initiation criterion satisfies Equation (13).
where < σ max > is the actual maximum principal stress of the element, the symbol "< >" means the damage initiation does not exist in a pure compression state; σ 0 max is the allowable maximum principal stress defined; f tol is the tolerance with a value of 0.05, which is a default value recommended by ABAQUS.
The damage evolution defines the degradation patterns of cohesive stiffness after damage initiation. The nonlinear exponential softening response is adopted to describe the degradation of cohesive stiffness for the analysis of elastoplastic fracture mechanics (EPFM). As the fracture energy is one of the fracture properties of materials, the evolution law based on energy method is selected. Once the crack propagation driving force exceeds the equivalent critical energy release rate (equivalent to fracture energy), the crack will extend. The equivalent critical energy release rate is calculated by the Benzeggagh-Kenane law (Equation (14)) [34].
where G eqC is the equivalent critical energy release rate; G I , G II and G III are the energy release rates of mode I, II and III cracks, respectively; G IC and G IIC are the critical energy release rates of mode I and II cracks, respectively; η is the exponent, to which the response is insensitive for isotropic failure. Fracture toughness is an index to measure the capacity of materials to prevent crack propagation, which is one of the inherent properties of materials. It can be expressed by a single parameter such as energy release rate G, stress intensity factor K, crack tip opening displacement CTOD and J-integral. The critical energy release rate is used as the fracture toughness here. The fracture failure mode is assumed to be isotropic, so G IC = G IIC . The fracture energy in the process of damage evolution is calculated on the basis of the fracture toughness of P91 steel. The tensile strength of P91 steel is taken as an approximation for the MPS. For both fracture toughness and tensile strength, temperature-dependent data are used to simulate the actual fracture behavior. The fracture parameters under high temperature are extrapolated by a linear interpolation method. The fracture properties of P91 steel are shown in Figure 10. With the development of damage evolution, the cohesive traction between crack surfaces decreases. When the cohesive traction is reduced to zero, the element fails, which means that the crack surfaces have been completely opened.
The fracture problems with damage definition have strong nonlinearity and discontinuity, which makes numerical solutions difficult to converge. ABAQUS provides a viscous regularization  Figure 10. Temperature-dependent fracture properties for P91 steel, used in the FE simulation, data from [35,36].
With the development of damage evolution, the cohesive traction between crack surfaces decreases. When the cohesive traction is reduced to zero, the element fails, which means that the crack surfaces have been completely opened.
The fracture problems with damage definition have strong nonlinearity and discontinuity, which makes numerical solutions difficult to converge. ABAQUS provides a viscous regularization method to solve this problem. By setting an appropriate small viscosity coefficient, the convergence of the fracture model will be significantly improved. Here, the viscosity coefficient is set to be 1.0 × 10 −5 .
The type IV crack with a size of 0.5 mm perpendicular to the surface is prefabricated in the HAZ of repair welding, as shown in Figure 11a. d is the distance between the crack and the fusion line. The initial crack is introduced into the FE model by the XFEM technique, and the contact attributes between the crack surfaces are set as frictionless for tangential behavior and "hard" contact for normal behavior, as shown in Figure 11b. Figure 10. Temperature-dependent fracture properties for P91 steel, used in the FE simulation, data from [35,36].
With the development of damage evolution, the cohesive traction between crack surfaces decreases. When the cohesive traction is reduced to zero, the element fails, which means that the crack surfaces have been completely opened.
The fracture problems with damage definition have strong nonlinearity and discontinuity, which makes numerical solutions difficult to converge. ABAQUS provides a viscous regularization method to solve this problem. By setting an appropriate small viscosity coefficient, the convergence of the fracture model will be significantly improved. Here, the viscosity coefficient is set to be 1.0 × 10 −5 .
The type IV crack with a size of 0.5 mm perpendicular to the surface is prefabricated in the HAZ of repair welding, as shown in Figure 11a. d is the distance between the crack and the fusion line. The initial crack is introduced into the FE model by the XFEM technique, and the contact attributes between the crack surfaces are set as frictionless for tangential behavior and "hard" contact for normal behavior, as shown in Figure 11b.

Weld Repair Specimens
Weld repair experiments of the P91 steel welded joint were carried out to verify the feasibility of the FE model. The size of the specimens was the same as that of the FE model. The chemical composition of experimental base metal is shown in Table 1, which measures up to ASME BPVC.II.A-2019 [37]. The selected weld metal was particularly suited for matching P91 steel. ER90S-B9 welding rod was used as a filler metal for GTAW, while E9015-B9 stick electrode was used as a filler metal for

Weld Repair Specimens
Weld repair experiments of the P91 steel welded joint were carried out to verify the feasibility of the FE model. The size of the specimens was the same as that of the FE model. The chemical composition of experimental base metal is shown in Table 1, which measures up to ASME BPVC.II.A-2019 [37]. The selected weld metal was particularly suited for matching P91 steel. ER90S-B9 welding rod was used as a filler metal for GTAW, while E9015-B9 stick electrode was used as a filler metal for SMAW. The chemical composition of the weld metal is shown in Table 2, which measures up to ASME BPVC.II.C-2019 [38]. Firstly, the base metal plates were connected by multilayer welding. GTAW was used for root welding, followed by SMAW used for filler welding and cover welding. The shielding gas for GTAW was argon. Then the welded joint was air-cooled to the ambient temperature. The metal near the FGHAZ of the welded joint, where it was prone to IV type cracking, was removed. Finally, repair welding technique was utilized to refill it by SMAW with the welding linear energy of 10 kJ/cm. The specific welding parameters for each layer and weld repair are shown in Table 3. The interlayer temperature was maintained at 200-300 • C during welding. The morphology of welding seams after initial welding and repair welding is shown in Figure 12. Firstly, the base metal plates were connected by multilayer welding. GTAW was used for root welding, followed by SMAW used for filler welding and cover welding. The shielding gas for GTAW was argon. Then the welded joint was air-cooled to the ambient temperature. The metal near the FGHAZ of the welded joint, where it was prone to IV type cracking, was removed. Finally, repair welding technique was utilized to refill it by SMAW with the welding linear energy of 10 kJ/cm. The specific welding parameters for each layer and weld repair are shown in Table 3. The interlayer temperature was maintained at 200-300 °C during welding. The morphology of welding seams after initial welding and repair welding is shown in Figure 12.

Verification of the FE Model
After the fabrication of weld repair specimens, the impact indentation strain method [39] was used to measure the residual stress. A spherical indenter was used to generate indentations in the form of impact loads at the location of measuring points. On the one hand, the impact action makes the material appear to have plastic flow deformation. On the other hand, the elastoplastic deformation caused by indentations itself changes under the action of residual stress. The total deformation after superposition of the two kinds of deformation is called strain increment. The relationship between the strain increment generated by the impact indentation strain method and elastic strain satisfies Equation (15).
where ∆ε is the strain increment; ε e is the elastic strain; a 1 , a 2 , a 3 are the coefficients obtained from the calibration curve; B is the strain increment in a zero-stress state. The biaxial resistance strain gauge containing two sensing grid elements with different axial directions, namely X-axis direction and Z-axis direction, was used to measure residual strain. A certain size of indentation was produced at the midpoint of the grid axis by the impact load. The value of strain increment was recorded by the strain recording instrument. The value of residual strain was determined according to the relationship between the calibrated elastic strain and the strain increment. The value of residual stress was calculated by Hooke's law.
The surface residual stress of specimens after both initial welding and repair welding was measured. Five measuring points were taken from each specimen, and the distance between adjacent measuring points on each side of the welding seam was 2.5 mm. The specific position of each measuring point is shown in Figure 13. The measurement of surface residual stress was finished by the KJS-3\4 indentation stress measurement system (Developed by Institute of Metal Research, Chinese Academy of Sciences, Shenyang, China). Before measurements, the accuracy of the indentation stress measurement system was calibrated with conventional residual stress measurement methods, such as X-ray diffraction, which measures up to GB/T 24179-2009. As shown in Figure 14a, the measurement system consists of two parts: portable intelligent stress tester and indentation generating device. The indentation generating device includes a striking rod (for generating indentation), a permanent magnet fixed base (for restricting the displacement of the measurement device) and a centering microscope (for centering the grid axis of strain gauge). The residual stress on-site measuring is shown in Figure 14b.
relationship between the strain increment generated by the impact indentation strain method and elastic strain satisfies Equation (15).
where Δε is the strain increment; εe is the elastic strain; a1, a2, a3 are the coefficients obtained from the calibration curve; B is the strain increment in a zero-stress state. The biaxial resistance strain gauge containing two sensing grid elements with different axial directions, namely X-axis direction and Z-axis direction, was used to measure residual strain. A certain size of indentation was produced at the midpoint of the grid axis by the impact load. The value of strain increment was recorded by the strain recording instrument. The value of residual strain was determined according to the relationship between the calibrated elastic strain and the strain increment. The value of residual stress was calculated by Hooke's law.
The surface residual stress of specimens after both initial welding and repair welding was measured. Five measuring points were taken from each specimen, and the distance between adjacent measuring points on each side of the welding seam was 2.5 mm. The specific position of each measuring point is shown in Figure 13. The measurement of surface residual stress was finished by the KJS-3\4 indentation stress measurement system (Developed by Institute of Metal Research, Chinese Academy of Sciences, Shenyang, China). Before measurements, the accuracy of the indentation stress measurement system was calibrated with conventional residual stress measurement methods, such as X-ray diffraction, which measures up to GB/T 24179-2009. As shown in Figure 14a, the measurement system consists of two parts: portable intelligent stress tester and indentation generating device. The indentation generating device includes a striking rod (for generating indentation), a permanent magnet fixed base (for restricting the displacement of the measurement device) and a centering microscope (for centering the grid axis of strain gauge). The residual stress on-site measuring is shown in Figure 14b.  Figure 15 presents a comparison of the residual stress within and around the welding seam zone obtained from experiments and FE simulations. It is shown that the simulated results are in good agreement with the experimental results. The maximum difference between the simulated value and the experimental value is less than 10%. Therefore, the FE model and program developed in this paper are suitable for the thermal-mechanical coupled simulation of repair welding process for the welded joint. The reason for the existing difference between the simulated results and the experimental results is that the 2-D numerical model is a simplification of the 3-D experimental model. The 2-D plane in the FE simulations actually corresponds to the middle plane of an infinite plate. Thus, only when the size of specimens in the length direction is infinitely large, the simulated  Figure 15 presents a comparison of the residual stress within and around the welding seam zone obtained from experiments and FE simulations. It is shown that the simulated results are in good agreement with the experimental results. The maximum difference between the simulated value and the experimental value is less than 10%. Therefore, the FE model and program developed in this paper are suitable for the thermal-mechanical coupled simulation of repair welding process for the welded joint. The reason for the existing difference between the simulated results and the experimental results is that the 2-D numerical model is a simplification of the 3-D experimental model. The 2-D plane in the FE simulations actually corresponds to the middle plane of an infinite plate. Thus, only when the size of specimens in the length direction is infinitely large, the simulated results may be very close to the experimental results.  Figure 15 presents a comparison of the residual stress within and around the welding seam zone obtained from experiments and FE simulations. It is shown that the simulated results are in good agreement with the experimental results. The maximum difference between the simulated value and the experimental value is less than 10%. Therefore, the FE model and program developed in this paper are suitable for the thermal-mechanical coupled simulation of repair welding process for the welded joint. The reason for the existing difference between the simulated results and the experimental results is that the 2-D numerical model is a simplification of the 3-D experimental model. The 2-D plane in the FE simulations actually corresponds to the middle plane of an infinite plate. Thus, only when the size of specimens in the length direction is infinitely large, the simulated results may be very close to the experimental results. (c) (d) Figure 15. Comparison of the residual stress obtained from FE simulations and experiments: (a) S11 stress after initial welding; (b) S33 stress after initial welding; (c) S11 stress after repair welding; (d) S33 stress after repair welding.

Thermal Analysis
The transient temperature contours at the intermediate instant and the end instant of repair welding for different linear energy are shown in Figure 16. The variation of linear energy has a significant influence on the temperature profile of the welding seam and HAZ. With the increase of linear energy, the peak temperature at the intermediate instant of repair welding increases. The peak temperature is 1819 °C for the linear energy of 10 kJ/cm, 2088 °C for the linear energy of 16 kJ/cm and 2320 °C for the linear energy of 25 kJ/cm. This is because higher linear energy involves longer loading times, thus more thermal energy is poured into the bulk metal during repair welding. After the intermediate instant of repair welding, due to the fact that the heat source starts to depart from the cross section, the heat generation rate decreases continuously. At the end instant of repair welding, the peak temperature for different linear energy has no obvious difference, all are around 1500 °C. In addition, with the increase of repair welding linear energy, the range of the HAZ expands. The higher linear energy brings higher overall temperature of the HAZ and poorer fracture properties.

Intermediate instant of repair welding
End instant of repair welding

Thermal Analysis
The transient temperature contours at the intermediate instant and the end instant of repair welding for different linear energy are shown in Figure 16. The variation of linear energy has a significant influence on the temperature profile of the welding seam and HAZ. With the increase of linear energy, the peak temperature at the intermediate instant of repair welding increases. The peak temperature is 1819 • C for the linear energy of 10 kJ/cm, 2088 • C for the linear energy of 16 kJ/cm and 2320 • C for the linear energy of 25 kJ/cm. This is because higher linear energy involves longer loading times, thus more thermal energy is poured into the bulk metal during repair welding. After the intermediate instant of repair welding, due to the fact that the heat source starts to depart from the cross section, the heat generation rate decreases continuously. At the end instant of repair welding, the peak temperature for different linear energy has no obvious difference, all are around 1500 • C. In addition, with the increase of repair welding linear energy, the range of the HAZ expands. The higher linear energy brings higher overall temperature of the HAZ and poorer fracture properties.
The transient temperature contours at the intermediate instant and the end instant of repair welding for different linear energy are shown in Figure 16. The variation of linear energy has a significant influence on the temperature profile of the welding seam and HAZ. With the increase of linear energy, the peak temperature at the intermediate instant of repair welding increases. The peak temperature is 1819 °C for the linear energy of 10 kJ/cm, 2088 °C for the linear energy of 16 kJ/cm and 2320 °C for the linear energy of 25 kJ/cm. This is because higher linear energy involves longer loading times, thus more thermal energy is poured into the bulk metal during repair welding. After the intermediate instant of repair welding, due to the fact that the heat source starts to depart from the cross section, the heat generation rate decreases continuously. At the end instant of repair welding, the peak temperature for different linear energy has no obvious difference, all are around 1500 °C. In addition, with the increase of repair welding linear energy, the range of the HAZ expands. The higher linear energy brings higher overall temperature of the HAZ and poorer fracture properties.

Mechanical Analysis
The calculated temperature field loads are used as the predefined field of stress calculation to perform the mechanical analysis. In Figure 18 seam is greater than that of the HAZ. This is because of the different material properties considered between the weld metal and the base metal. Figure 17 shows the temperature distribution at a depth of 0.5 mm in the HAZ at the end instant of repair welding for different linear energy. The overall temperature in the HAZ increases with the increment of linear energy. The peak temperature of the HAZ is 961 • C for the linear energy of 10 kJ/cm, 1130 • C for the linear energy of 16 kJ/cm and 1258 • C for the linear energy of 25 kJ/cm. The temperature distribution of the HAZ is gentler for high linear energy. Q = 25 kJ/cm  Figure 17 shows the temperature distribution at a depth of 0.5 mm in the HAZ at the end instant of repair welding for different linear energy. The overall temperature in the HAZ increases with the increment of linear energy. The peak temperature of the HAZ is 961 °C for the linear energy of 10 kJ/cm, 1130 °C for the linear energy of 16 kJ/cm and 1258 °C for the linear energy of 25 kJ/cm. The temperature distribution of the HAZ is gentler for high linear energy.

Mechanical Analysis
The calculated temperature field loads are used as the predefined field of stress calculation to perform the mechanical analysis. In Figure 18 The S11, S33 and Von Mises stresses at a depth of 0.5 mm in the HAZ at the end instant of repair welding for different linear energy are shown in Figure 19. With the increase of linear energy, the range of high stress area of S11, S33 and Von Mises expands and moves towards the direction far away from the fusion line. The peak S33 stress and Von Mises stress has little change with the increase The S11, S33 and Von Mises stresses at a depth of 0.5 mm in the HAZ at the end instant of repair welding for different linear energy are shown in Figure 19. With the increase of linear energy, the range of high stress area of S11, S33 and Von Mises expands and moves towards the direction far away from the fusion line. The peak S33 stress and Von Mises stress has little change with the increase of linear energy, while the peak S11 stress increases slightly with the increment of linear energy. The S11 stress and S33 stress are tensile stress near the fusion line and compressive stress far away from the fusion line.  Figure 20 shows the S11, S33 and Von Mises stress distribution at a depth of 0.5 mm in the HAZ at the end instant of cooling for different linear energy. Both the S33 stress and Von Mises stress have similar peaks for different linear energy. The peak S33 stress is about 525 MPa, and the peak Von Mises stress is about 480 MPa. In the cooling process, the S11 stress near the fusion line in the HAZ changes from tensile stress to compressive stress, and the S11 stress far away from the fusion line changes from compressive stress to tensile stress. Almost all the S33 stress in the HAZ changes from compressive stress to tensile stress, and the peaks change from −500 MPa to +525 MPa, which is caused by the shrinkage of welding seam during cooling. The remarkable tensile stress may lead to the propagation of initial cracks in the HAZ. Moreover, with the increase of linear energy, the range of the high stress area in the HAZ expands. The range of the high stress area in the HAZ is 4 mm for the linear energy of 10 kJ/cm, 6 mm for the linear energy of 16 kJ/cm and 10 mm for the linear energy of 25 kJ/cm.  Figure 20 shows the S11, S33 and Von Mises stress distribution at a depth of 0.5 mm in the HAZ at the end instant of cooling for different linear energy. Both the S33 stress and Von Mises stress have similar peaks for different linear energy. The peak S33 stress is about 525 MPa, and the peak Von Mises stress is about 480 MPa. In the cooling process, the S11 stress near the fusion line in the HAZ changes from tensile stress to compressive stress, and the S11 stress far away from the fusion line changes from compressive stress to tensile stress. Almost all the S33 stress in the HAZ changes from compressive stress to tensile stress, and the peaks change from −500 MPa to +525 MPa, which is caused by the shrinkage of welding seam during cooling. The remarkable tensile stress may lead to the propagation of initial cracks in the HAZ. Moreover, with the increase of linear energy, the range of the high stress area in the HAZ expands. The range of the high stress area in the HAZ is 4 mm for the linear energy of 10 kJ/cm, 6 mm for the linear energy of 16 kJ/cm and 10 mm for the linear energy of 25 kJ/cm.

XFEM Analysis
Under repair welding thermal shock, the welding seam and HAZ deform plastically, so the fracture mode is elastoplastic fracture. The whole fracture process has strong nonlinearity, including not only the nonlinearity of material properties, but also the geometric nonlinearity consisting of large thermal deformation and crack propagation.
The crack at the position of d = 2 mm under the repair welding linear energy of 10 kJ/cm is selected to analyze and the simulated crack propagation paths are shown in Figure 21. It is noteworthy that the crack propagation occurs during the cooling process, which is in agreement with the mechanical analysis. The crack propagation direction is almost perpendicular to the surface due to the control of MPS criterion and there is a certain tendency for the crack to deflect to the welding seam during the propagation process. A total of two elements are extended and the extracted crack propagation length is 0.52 mm. Figure 21a shows the MPS contours for the crack zone at different transient times. t = 3.96 s is the end instant of repair welding and the start instant of cooling. At the start instant of cooling, the MPS around the crack tip is lower than the allowable MPS at the corresponding temperature according to the damage initiation criterion, so the crack does not propagate. With the cooling process going on, the tensile stress increases continuously, which leads to the increment of the MPS around the crack tip. At t = 4.43 s, the crack damage initiation condition is reached, and the new crack surface begins to form. After the formation of the new crack surface,

XFEM Analysis
Under repair welding thermal shock, the welding seam and HAZ deform plastically, so the fracture mode is elastoplastic fracture. The whole fracture process has strong nonlinearity, including not only the nonlinearity of material properties, but also the geometric nonlinearity consisting of large thermal deformation and crack propagation.
The crack at the position of d = 2 mm under the repair welding linear energy of 10 kJ/cm is selected to analyze and the simulated crack propagation paths are shown in Figure 21. It is noteworthy that the crack propagation occurs during the cooling process, which is in agreement with the mechanical analysis. The crack propagation direction is almost perpendicular to the surface due to the control of MPS criterion and there is a certain tendency for the crack to deflect to the welding seam during the propagation process. A total of two elements are extended and the extracted crack propagation length is 0.52 mm. Figure 21a shows the MPS contours for the crack zone at different transient times. t = 3.96 s is the end instant of repair welding and the start instant of cooling. At the start instant of cooling, the MPS around the crack tip is lower than the allowable MPS at the corresponding temperature according to the damage initiation criterion, so the crack does not propagate. With the cooling process going on, the tensile stress increases continuously, which leads to the increment of the MPS around the crack tip. At t = 4.43 s, the crack damage initiation condition is reached, and the new crack surface begins to form. After the formation of the new crack surface, the stress concentration around the crack tip is released and the MPS decreases. A new driving force is required to promote the further crack propagation. At t = 5.22 s, the new driving force makes the crack extend further.
ranging from 0 to 1.0. A value of 0 means that the material has no damage. A value of 1.0 means that the crack is fully opened and the cohesive traction between the crack surfaces is zero. When the value is greater than 0 and less than 1.0, it means that the cohesive crack has already extended and needs extra driving force to fully open. Although the new cracks are not fully opened, in the complex service environment, the cohesive cracks containing damage are likely to open completely after being driven by the extra driving force, which threatens the structural integrity of weld repaired components.   According to the crack propagation paths in Figure 21, two elements in the initial crack front are selected to analyze. The transient S11, S33 and Von Mises stresses of the elements under repair welding thermal shock are shown in Figure 22. The first 3.96 s of transient time is the repair welding process. In this process, the HAZ is compressed due to the outward expansion of the welding seam. The S11 and S33 stresses are compressive stress, and the value of compressive stress increases with the promoting of repair welding. After repair welding, the welding seam begins to cool and shrink, which results in the tensile stress in the HAZ. The S11 and S33 stress changes from compressive stress to tensile stress in the early stage of cooling and the crack propagation occurs in this stage.
Metals 2020, 10, x FOR PEER REVIEW 20 of 26 According to the crack propagation paths in Figure 21, two elements in the initial crack front are selected to analyze. The transient S11, S33 and Von Mises stresses of the elements under repair welding thermal shock are shown in Figure 22. The first 3.96 s of transient time is the repair welding process. In this process, the HAZ is compressed due to the outward expansion of the welding seam. The S11 and S33 stresses are compressive stress, and the value of compressive stress increases with the promoting of repair welding. After repair welding, the welding seam begins to cool and shrink, which results in the tensile stress in the HAZ. The S11 and S33 stress changes from compressive stress to tensile stress in the early stage of cooling and the crack propagation occurs in this stage. In order to investigate the influence of repair welding thermal shock on the cracks at different distances from the fusion line for three different linear energy values, the initial cracks with a size of 0.5 mm perpendicular to the surface were prefabricated at the positions where the distances from the fusion line were 2 mm, 4 mm, 6 mm, 8 mm and 10 mm, respectively. Figure 23 shows the crack propagation predicted by XFEM at different positions for three different linear energy values. All the crack propagation occurs in the early stage of cooling. For cracks at the same location, with the increase of linear energy, the crack propagation length and the number of damaged elements increase. This is because the higher linear energy brings the higher temperature of the HAZ, which leads to poorer crack resistance. The longest crack propagation appears at the position of d = 2 mm under the linear energy of 25 kJ/cm. The extracted crack propagation length is 1.1 mm, which is almost twice of the initial crack length. In this case, the fully opened crack propagation has occurred, and the crack is likely to develop into the macrocrack, even penetrate through the whole wall thickness. For the same linear energy, the crack propagation length decreases with the increment of distance from the fusion line. Both the cracks with d > 6 mm for the linear energy of 10 kJ/cm and the In order to investigate the influence of repair welding thermal shock on the cracks at different distances from the fusion line for three different linear energy values, the initial cracks with a size of 0.5 mm perpendicular to the surface were prefabricated at the positions where the distances from the fusion line were 2 mm, 4 mm, 6 mm, 8 mm and 10 mm, respectively. Figure 23 shows the crack propagation predicted by XFEM at different positions for three different linear energy values. All the crack propagation occurs in the early stage of cooling. For cracks at the same location, with the increase of linear energy, the crack propagation length and the number of damaged elements increase. This is because the higher linear energy brings the higher temperature of the HAZ, which leads to poorer crack resistance. The longest crack propagation appears at the position of d = 2 mm under the linear energy of 25 kJ/cm. The extracted crack propagation length is 1.1 mm, which is almost twice of the initial crack length. In this case, the fully opened crack propagation has occurred, and the crack is likely to develop into the macrocrack, even penetrate through the whole wall thickness. For the same linear energy, the crack propagation length decreases with the increment of distance from the fusion line. Both the cracks with d > 6 mm for the linear energy of 10 kJ/cm and the cracks with d > 8 mm for the linear energy of 16 kJ/cm do not propagate. The reason is that the lower the linear energy is, the smaller the scope of the HAZ is. Once the crack is far away from the fusion line, the driving force is not enough to promote the crack propagation. It should be noted that the area where cracks propagate is the high stress area where the peak stress is located. Almost all the directions of crack propagation are perpendicular to the surface and tend to deflect to the welding seam, which may be related to the high stress in the welding seam.   Since the crack resistance of materials is related to temperature, the temperature at the crack tip is one of the key factors affecting crack propagation. On the one hand, the higher the temperature is, the poorer the crack resistance is, and the cracks are easier to extend under the same driving force. On the other hand, the larger the temperature gradient is, the greater the thermal stress is, so the driving force is sufficient for crack propagation. Figure 24 shows the thermal cycle curves at different crack tips in the HAZ for three different linear energy values. For the same linear energy, with the decrease of the distance from the fusion line, the difference of the peak temperature at the crack tip becomes larger. As a result, the closer the crack is to the fusion line, the greater the temperature gradient and temperature change around the crack tip, the greater the driving force and the longer the crack propagation length. For the same crack tip, the peak temperature increases with the increment of linear energy. Therefore, the greater the linear energy is, the poorer the crack resistance of materials in the HAZ is. In the early stage of cooling, the temperature value in the HAZ is higher and the temperature difference in the HAZ is larger, so crack propagation occurs in the early stage of cooling. With the cooling process going on, the temperature at the crack tip at different positions gradually tends to become consistent. Since the crack resistance of materials is related to temperature, the temperature at the crack tip is one of the key factors affecting crack propagation. On the one hand, the higher the temperature is, the poorer the crack resistance is, and the cracks are easier to extend under the same driving force. On the other hand, the larger the temperature gradient is, the greater the thermal stress is, so the driving force is sufficient for crack propagation. Figure 24 shows the thermal cycle curves at different crack tips in the HAZ for three different linear energy values. For the same linear energy, with the decrease of the distance from the fusion line, the difference of the peak temperature at the crack tip becomes larger. As a result, the closer the crack is to the fusion line, the greater the temperature gradient and temperature change around the crack tip, the greater the driving force and the longer the crack propagation length. For the same crack tip, the peak temperature increases with the increment of linear energy. Therefore, the greater the linear energy is, the poorer the crack resistance of materials in the HAZ is. In the early stage of cooling, the temperature value in the HAZ is higher and the temperature difference in the HAZ is larger, so crack propagation occurs in the early stage of cooling. With the cooling process going on, the temperature at the crack tip at different positions gradually tends to become consistent.  Figure  25. The distribution of heat flux is similar to the thermal cycle curves at crack tips, which indicates that the temperature at crack tips is proportional to the thermal energy density through the crack  Figure 25. The distribution of heat flux is similar to the thermal cycle curves at crack tips, which indicates that the temperature at crack tips is proportional to the thermal energy density through the crack surfaces. For the same crack surface, the peak heat flux increases with the increment of linear energy. For the same linear energy, in repair welding process and early stage of cooling, the closer the crack surface is to the fusion line, the greater the heat flux gradient is and the greater the heat flux change is. surfaces. For the same crack surface, the peak heat flux increases with the increment of linear energy. For the same linear energy, in repair welding process and early stage of cooling, the closer the crack surface is to the fusion line, the greater the heat flux gradient is and the greater the heat flux change is.

Conclusions
In this study, the XFEM based on a cohesive crack model is used to establish the repair welding model of P91 steel welded joint containing initial cracks. The MPS criterion and energy method are used to define the damage initiation and damage evolution, respectively. Considering the temperature-dependent material physical properties and fracture properties, the effect of repair welding thermal shock on the crack propagation behavior in the HAZ is simulated. The influence of different repair welding linear energy and different crack positions on the cracking features is analyzed. The main conclusions are as follows: (1) With the increase of linear energy, the range of the HAZ expands and the overall temperature of the HAZ rises. With the increase of linear energy, the range of S11, S33 and Von Mises high stress area expands. (2) In the repair welding and cooling process, the material in the HAZ undergoes a change from compressed state to tensioned state due to the expansion and shrinkage of the welding seam. It is the direct cause of crack propagation in the HAZ.

Conclusions
In this study, the XFEM based on a cohesive crack model is used to establish the repair welding model of P91 steel welded joint containing initial cracks. The MPS criterion and energy method are used to define the damage initiation and damage evolution, respectively. Considering the temperature-dependent material physical properties and fracture properties, the effect of repair welding thermal shock on the crack propagation behavior in the HAZ is simulated. The influence of different repair welding linear energy and different crack positions on the cracking features is analyzed. The main conclusions are as follows: (1) With the increase of linear energy, the range of the HAZ expands and the overall temperature of the HAZ rises. With the increase of linear energy, the range of S11, S33 and Von Mises high stress area expands. (2) In the repair welding and cooling process, the material in the HAZ undergoes a change from compressed state to tensioned state due to the expansion and shrinkage of the welding seam. It is the direct cause of crack propagation in the HAZ. (3) Crack propagation occurs in the early stage of cooling. The crack propagation direction is roughly perpendicular to the upper surface and has a tendency to deflect to the welding seam. (4) For cracks at the same position, with the increase of linear energy, the crack propagation length increases. For the same linear energy, with the distance from the fusion line increases, the crack propagation length decreases.
(5) After a certain distance from the fusion line, the cracks in the HAZ do not extend. This distance is the critical distance for the crack propagation. With the increase of linear energy, the critical distance of crack propagation increases. The critical distance of crack propagation is consistent with the high stress area after repair welding. Therefore, for repair welding, low linear energy is preferred to constrain the cracking behavior in the HAZ.