A thermomechanical finite element tool for Leak‐before‐Break analysis

A new finite element tool is presented, which utilises the extended FEM (XFEM) to model leaks through cracks. Heat flux and pressure boundary conditions are imposed on the crack in the form of jump terms. Enrichments are chosen to model either strong or weak discontinuities across the crack, as well as singularities at the crack tips. Excellent convergence rates are achieved for both the thermal and mechanical models, where errors are calculated relative to analytical solutions derived for this specific problem. A more general temperature approximation is also presented, which makes no assumptions about the continuity of temperature or heat flux across the crack. Results indicate that this is a robust way of modelling the temperature of a plate containing a crack with or without a leaking fluid. Thermomechanical simulations were then carried out to demonstrate the applicability of the FEM for analysing leak rates in nuclear reactor primary pipework. A two‐phase flow model based on the Henry–Fauske model is chosen for the fluid aspect, and this is coupled to the structure through a convection law. Crack closure is shown to reduce the leak rate by up to 40%. © 2014 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons, Ltd.


INTRODUCTION
In nuclear power plants, the low probability of failure is justified by multilegged safety cases. There are the standard manufacturing quality assurance requirements, as well as inspection methods and defect tolerance assessments for small flaws. In addition to this, there are defence in depth arguments such as the Leak-before-Break (LbB) philosophy. It says that if there is a crack in a component, then it will be possible to detect a leak within a short enough time frame such that any corrective action may be taken before a catastrophic failure occurs [1]. Essentially, this eliminates the possibility of a guillotine type fracture, which would normally be mitigated against with pipe whip restraints and jet impingement shields. The economic benefits of not imposing these protective measures are made possible by safety arguments such as Leak-before-Break.
The full Leak-before-Break procedure requires an assessment of the critical crack size of a component. This is the size of crack at which a guillotine type fracture will occur. This can be calculated using a failure assessment diagram methodology such as R6 [2]. It is then possible to calculate the crack opening area (COA) for a crack size much smaller than the critical crack size. Analytical equations can be used to obtain COA, and if further accuracy is required, the FEM could be used. Using this COA, it is then possible to obtain the leak rate. If this leak rate is greater than a specified detection limit, a Leak-before-Break case can be made. The biggest uncertainty in a Leak-before-Break assessment is the calculation of leak rate. This paper focusses on the calculation of leak rates 679 using a new finite element tool. This tool is designed to reduce error by incorporating crack opening and fluid calculations within elements. The fluid mechanics of the leaking fluid will be coupled to the solid via mechanical and thermal boundary conditions. In particular, the effect of heat transfer from the fluid to the crack face is investigated as it is hypothesised to reduce the COA. Therefore, the focus of this work is on the thermal exchange between the fluid and the structure, in contrast to [3], which focusses more on the mechanical aspects of a similar coupled problem.
The high inlet-to-outlet pressure ratios that occur when a through wall crack appears generally necessitate the use of critical flow models. These are either single or two-phase flow, depending on whether a gas or water coolant is used. CFD gives leak rates accurately but are impractical for LbB assessments due to the complexity of the models. Therefore, analytical or simplified 1D finite difference models are preferred as they give good accuracy relatively cheaply. The single phase mass flow rate can be derived by considering the mass, energy and momentum equations, as well as an equation of state for an ideal gas. This leads to closed-form solutions that are straightforward to solve. For the two-phase case, a model based on Henry-Fauske is employed [4]. This treats the mixture as a single fluid with properties based on the mixture quality. There is also a consideration of the thermal non-equilibrium between phases, which occurs because of the very short time the fluid is in the crack. The error is exaggerated by the fact that there is uncertainty in the crack shape [5].
It is apparent therefore that Leak-before-Break has both structural and fluid mechanics considerations. It is a present practice to analyse these two aspects independently [1]. However, modern computational techniques are making a coupled approach more feasible. The motivation for a coupled approach stems from the fact that thermal heat transfer from the fluid to the crack walls is known to affect leak rate [6]. Two-phase flow is particularly sensitive to wall heat transfer [7], this is due to the steam evolution as pressure and temperature change through the thickness. Also, the additional heating of the crack walls due to a leaking fluid causes thermal expansion, which changes the crack opening displacement (COD). So it is clear that there should be a link between fluid and structure through a thermal coupling. In addition to this, there should be a consideration for the end user, so a model that can be efficiently used in an industrial context is essential. To facilitate this requirement, a new finite element tool has been developed, which accounts for all the physics that affects leak rate. The approach outlined is founded on the extended FEM (XFEM), as this provides a framework to build an element that can capture the essential behaviour, as well as being simple to use once implemented into a finite element code.

LITERATURE REVIEW
The extended FEM, first introduced by Belytschko [8], is a development of the FEM, which enables efficient modelling of problems with discontinuities or singularities. The method utilises the partition of unity concept that was presented by Melenk and Babuska [9]. With the addition of enrichment functions, which give additional information in the regions where the conventional finite element may require highly refined meshes, optimum convergence rates can be achieved. Such physical situations where this may arise are cracks in solids, where there is a discontinuity in displacement across the crack, and the singularity occurs in the strain field at the crack tip. The enrichment functions used to describe this behaviour are a Heaviside step function along the crack and a set of functions derived from the Westergaard solution at the crack tip.
Enrichment functions should capture the behaviour in specific regions using exact analytical solutions. For the mechanical model, the well-known elastic behaviour of cracks gives a solution close to the crack tip as a series expansion. The four terms used in XFEM have proven to give good accuracy and convergence rates [10]. Although these expressions were based on the traction free crack assumption, the principle of linear superposition means that they can be used for a constant pressure along the crack faces. This is shown to give good accuracy when comparing the COD from the FE approximation to analytical expressions. The thermal consideration is not as well established, although it is known that there is a singularity in heat flux at the crack tip when considering a homogeneous conducting material. Kit [11] solved the Helmholtz equation for a plate with a crack subject to various boundary conditions. These were either keeping a constant temperature or heat flux along the crack or imposing a jump condition in both quantities. This assumes a heat transfer in the out-of-plane direction,´, which leads to the Poisson equation. However, this work takes no account of the´-direction and hence, an insulated condition is assumed. Therefore, the equation to solve is simply Laplace's equation, which is carried out by using the fundamental solution and applying single layer potential theory [12] (see Appendix). This model provides a means to quantify the convergence behaviour and accuracy of the FE model.
The integration involved with XFEM is complicated by the fact that non-polynomial functions are acting within elements, thereby rendering traditional Gaussian integration techniques inaccurate. Various methods have been suggested to circumvent these problems, a commonly used technique is to decompose the elements into smaller subregions and perform the integration separately within these regions [13]. Other methods have been proposed by Bordas et al. who uses strain smoothing techniques [14] and Schwarz-Christoffel mapping [15] to account for the sharp discontinuity. In [14], it was found that higher levels of accuracy could be achieved with strain smoothing when the enrichment functions were polynomial; however, for non-polynomial functions, the accuracy was less than element subdivision techniques. The results in [15] are shown to be accurate. The integration scheme employed in this study is based on the method of Fries [13], where cracked elements are subdivided into quadrangles and crack tip elements are divided into six triangles each with a vertex at the crack tip. The first use of a decomposition of elements into subquadrilaterals independent of the position of the discontinuity was carried out in [16]. In [17], specific polar integration and geometrical enrichment for singular enrichment is introduced. Elguedj et al. [18] present the use of a fixed subintegration rule for XFEM crack growth. Laborde gives a new polar integration rule and convergence study of the XFEM in [19]. Standard Gauss quadrature in elements containing the discontinuity is used in [20].
Applications of XFEM to problems that are relevant to the one considered here are now discussed. Thermomechanical aspects of XFEM are investigated by Duflot [21]. In this work, the boundary conditions on the crack are constant temperatures using Lagrange multipliers or an insulated condition. A jump in temperature based on a thermal resistance model was considered by Yvonnet et al. [22], which assume a constant heat flux across the crack. This is a special case of the weak form presented in [23], which makes no assumptions about the flux across the crack. The thermal model that will be presented in the paper assumes a heat flux along the crack imposed by a Robin-type boundary condition, with change in heat flux accounted for by the relevant enrichment function.
Extended FEM is applied to coupled fluid solid problems in [24] and [3]. A microscale viscous flow model is coupled to the macroscale finite element model via pressure and displacement requirements imposed on the discontinuity. These requirements are satisfied by exploiting the partition-of-unity property of finite element shape functions. The work presented here develops this idea by imposing jump conditions across the discontinuity due to the pressure and heat flux from the leaking fluid.

MATHEMATICAL FORMULATION
The mathematical derivation of the new finite element tool is presented in this section. Steady-state conditions are assumed in this analysis, as well as small deformation theory for the strains. The crack is considered to be a 1D line discontinuity in a 2D domain, see Figure 1. The crack is represented by a single contour centrally positioned between the two crack faces, see Figure 2. The use of a single contour gives rise to jump conditions, which are necessary to represent sharp changes in surface traction and heat transfer at the crack. The governing equations are the equilibrium equation and heat equation on a domain bounded by @ , respectively, and are given by  where f i is the body force, and Q is the heat source (which in this case will be zero, so the heat equation reduces to the Laplacian T ;jj D 0). The ;j denotes differentiation with respect to the j th coordinate. The constitutive equations in thermoelasticity are given by Á where C is the isotropic fourth-order tensor, is the strain tensor,˛T is the thermal expansion coefficient, k is the thermal conductivity, T is the change in temperature with respect to some reference temperature, and ı ij is the Kronecker delta. The superscripts e and t represent the elastic and thermal components of the total strain and u is the displacement vector. Neumann and Dirichlet boundary conditions in the form of the traction/heat flux and displacement/temperature are imposed on the body, with additional conditions on the crack faces where p and q c are obtained from the pressure and heat flux of the leaking fluid, respectively. I is the second-order identity matrix. The boundary conditions in 3 satisfy T S q D and T T q D ;, u S t D and u T t D ;, c t , N t D pn c for the crack surface. In the case of a thermoelastic crack with convection on the faces c q , q D q c . The jump condition is defined to be OEOE D j 1 C j 2 ‡ . The normals n 1 and n 2 are defined on each face of the crack 1 and 2 , and n is the global normal, as shown in Figure 1. The relationship between the crack face normals is defined to be n 1 D n 2 . The crack is taken to be a single contour c D 1 S 2 as shown in Figure 2. This is physically justified for small crack openings, where the crack is taken to be the limit as the two faces converging, with a jump term across the contour to account for the pressure or heat flux acting on the walls.
Integration of the variational terms f ıu and ıTQ over the domain and on application of the divergence theorem provides the weak form of equation (1) with the prescribed boundary conditions, that is, where D is the conductivity matrix and the following solution spaces u 0 and T 0 of the test functions ıT and ıu are given by The weak form for the thermal problem considered here is derived from the most general form given in [23]. In this analysis, the body contains an interface that describes a curve intersecting with the boundary. The curve is not a feature in the weak form for this application because the physics involved is limited to heat fluxes at surfaces and interfaces.
The formulation is similar to a traditional FE formulation, but in this case, the stiffness matrix has an additional contribution arising from an integral along the crack. The heat vector also has a contribution due to the bulk fluid temperature along the crack. The discretized form is as follows: For the temperature field, it is possible to model a jump in heat flux on either side of the crack. This requires an enrichment function that can account for a change in the derivative of temperature across the crack. Moës' bimaterial formulation [25] is used to account for this where I is the normal distance from the interface to the node. If the crack is assumed insulated, that is, when there is no fluid in the gap, it is more appropriate to use the Heaviside function to give a jump in temperature. This does not account for different fluxes on each side of the crack however. The enrichment function in 7 is used in 10 to account for the jump in heat flux at the crack. So the temperature discretization is i and T c i are additional unknowns that need to be solved for. S c are nodes whose parent element is completely cut by the crack, and S t are nodes whose parent element is partially cut by the crack. Equation 9 is the approximation for a discontinuous temperature across the crack, and Equation 10 is for a discontinuous heat flux across the crack. H is the step function which is C1 on one side of the crack and 1 on the other. i is the crack tip enrichment, and is the level set function employed to locate the crack.
For a more general formulation, the temperature can be approximated using both the weak and strong discontinuity enrichments, that is, An important property of a traditional finite element approximation is the Kronecker delta property of the shape functions, that is, N i .x j / D ı ij . This ensures that the computed unknowns u I are the actual values required of u.x/ at nodes I , that is, u.x I / D u I . It also makes imposing Dirichlet boundary conditions u.x/ simple: u I D N u.x I / for x 2 u . The shifting of the enrichment functions ( j .x// j .x I /) is carried out to ensure that the finite element scheme has this property. The weak discontinuity enrichment is defined in such a way that is has this feature (see Equation 7).
The crack tip enrichments functions describe the singular behaviour at the crack tip and are based on the Westergaard solution.
where r and Â are polar coordinates, with origin at the crack tips. For the thermal field, it is appropriate to use the following enrichment function when a constant temperature is imposed on the crack faces which is suitable because it is constant when switching from to , which are the respective crack faces. Hence, the temperature is continuous across the crack when this enrichment function is used. Taking the derivative with respect to theta gives sinusoidal behaviour, which is discontinuous across the crack, giving a discontinuity in heat flux. When the crack faces are subject to a constant heat flux, the following enrichment function is appropriate: This is discontinuous from to , hence, the temperature can be discontinuous across the crack. The pressure term given in the weak form (5) is now analysed. There is a non-zero traction on each face of the crack, which is equal in magnitude but in opposite directions. The pressure term becomes Z where ıu n is a variation along the crack, and symmetry is assumed on either side of the crack. Similarly, for the thermal case, there is a heat flux on each face of the crack, equal in magnitude but in opposite directions. This flux is defined by a convection law , T bulk and h are known a priori but T must be obtained from the solution; therefore, it is incorporated into the left-hand side.
With the current formulation, it is possible to model either a discontinuity in heat flux or temperature. This is because the temperature approximation can contain the Heaviside or the level set enrichment, which accounts for the strong or weak discontinuity, respectively. Using the aforementioned convection law means that prescribing a temperature along the crack is simply a case of fixing the heat transfer coefficient to be very large. This would give a constant temperature along the crack without the need to use Lagrange multipliers. Currently, this is what is done for imposing Dirichlet conditions in XFEM [21].

FINITE ELEMENT EQUATIONS
Inserting the test functions and trial functions into the weak form leads to the fully coupled mechanical system of equations that can be solved simultaneously with respect to the mechanical and thermal DOFs  . .
where the subscripts M , T and TM correspond to the mechanical, thermal and thermomechanical components, respectively. Mechanical and . The second-order elasticity tensor .C /, for plane stress, is given by where E is the elastic modulus and is the Poisson's ratio. Thermomechanical where 2 ¹0; 1; 2; 3; 4º andˇ2 ¹0; 1º, where 0 corresponds to the weak or strong discontinuity enrichment and 1 is either a cosine or sine enrichment. If the enrichment in Equation (11) is used, 2 ¹0; 1; 2; 3º.

Thermal
The conductivity matrix is made up of the standard matrix plus the contribution from the crack face condition. and where ;ˇ2 ¹0; 1º or ¹0; 1; 2; 3º depending on the approximation used. .D/ contains the conductivity, k, that is, The structure matrix for two dimensions reduces to where K T is a 12 12 matrix, K M ij are 24 24 matrices and K TM ij are 24 12 matrices. In the case of a four-noded quadrilateral linear element, the standard shape functions in a reference domain are as follows: where ; Á 2 OE 1; 1 . The reference domain is mapped to the real element under a smooth map.

FLUID MODELS
The crack is modelled as a 1D channel or pipe with rough walls, and the dimensions can change through the thickness of the vessel, see Figure 3. The crack opening displacement is denoted by w; l is the crack length and t is the length of the channel. The crack opening can either diverge or converge through the channel, current models account for a linear change in this. Choked flow must be considered here as the difference in pressure from the inside of the vessel to the outside is very high. The crack channel shown in 3 is idealised and it is known that in reality, the channel may be more complex. As the fluid flows around corners through the crack length, there are associated pressure losses. The pressure losses can be split up into three distinct categories: surface roughness, changes in direction and recirculation. These are captured in a discharge coefficient C D , which is derived in [5]. The mass flow rate for single phase flow is given by where p 0 and 0 are the inlet pressure and density and A is the average crack opening area. For the two-phase flow, the situation is more complicated and the SQUIRT code is used [26]. SQUIRT uses the Henry-Fauske [4] homogeneous non-equilibrium model, which is essentially the solution of two simultaneous, nonlinear equations. One of these equations accounts for the mass flux, and the other is the pressure balance due to losses along the flow path, including phase change acceleration, area change, friction, entrance and corner losses. The model accounts for the fact that there is thermal non-equilibrium between the two phases during the time the fluid is in the crack. This non-equilibrium mixture quality X c is in the form of an empirical law, which is fitted to experimental data. The mass flux and pressure constraints form a simultaneous nonlinear system, which can be solved by Newton-Raphson iteration, that is, where G c is the critical mass flux, gc and Lc are the specific enthalpies of the vapour and liquid, respectively, is the isentropic exponent and X c and X E are the non-equilibrium and equilibrium mixture quality given by where N D 20 for X E < 0:05 and N D 1 for X E > 0:05. The constant B D 0:0523. The subscripts 0, c, e, a, f , k and aa denote the vessel pressure, critical, entrance, phase change acceleration, frictional, corners and area change acceleration pressure losses, respectively. Properties of steam are obtained from IAPWS IF-97 steam tables [27]. The thermal nonequilibrium is accounted for by a relaxation condition, which is an expression derived from experimental data [28]. This is the most practical way of accounting for heat transfer between the liquid and the vapour. To model this process in detail would be prohibitively expensive. The heat transfer coefficient is calculated from on the flow rate, channel dimensions and fluid properties. The Dittus-Boelter correlation can be used as a good approximation N u D 0:023 Re 0:8 P r 0:4 (33) where where D h is the hydraulic diameter of the channel, and C p , k f and are the specific heat, conductivity and viscosity of the fluid respectively.

INTEGRATION
Numerical integration is complicated somewhat by the fact that there are discontinuous and singular functions in the approximation. These non-polynomial functions cannot be accurately integrated using standard Gaussian quadrature. Therefore, there have been various methods proposed that are designed to overcome these problems, one such method is to subdivide elements into smaller integration regions [13] [29] [8]. Crack tip elements are broken down into triangles, which have a vertex at the crack tip, crack front elements are dissected into subquadrilaterals, which are aligned with the level set. For integration along the level set itself, there is a 1D integral performed and the integration points are distributed along the level set.
The method of subdivision is chosen here because of its robustness and accuracy, as well as it is being relatively simple to implement. Figure 4 shows the integration points in a cracked element, the quadrangles have a side that aligns with the discontinuity. The level set integration points are also shown and the contribution from both of these sets of integration points gives the total integral for the element. Figure 5 shows the integration points in the crack tip. Six triangles decompose the element, all of which have a vertex at the crack tip, within each of these triangles is five integration points. Along the crack, an algorithm populates the segment of the level set that is cracked.

IMPLEMENTATION
The 2D code developed here has the ability to model plates with cracks placed anywhere and with any configuration without any need to change the mesh. The user supplies the dimensions of the plate with a mesh size along with the crack tips and enrichment radius. The code then generates the mesh and level set then finds the enriched elements. After all the initial conditions have been input, the code runs and outputs a leak rate as well as any other quantities of interest. The coupling between the fluid and the structure can be summarised with the help of Figure 6. The initial temperature field is calculated and input to the thermomechanical model and a crack opening displacement is  obtained from the solution. This COD is used in the leak rate calculation to give the leak rate and heat transfer coefficient. A new temperature is then calculated and fed into the thermomechanical model to obtain a new COD. This process is repeated until a steady COD or leak rate is obtained.

One-dimensional case
For the 1D case, a single element is considered, which contains a cohesive crack. This is necessary to maintain the integrity of the element, so a traction separation law that relates the stress to the displacements on each side of the crack [30] is used, that is, Hooke's law is assumed applicable everywhere else. Expressions in the following text are obtained by considering the stresses and strains in the element and following some manipulation give 690 P. GILL AND K. DAVEY Figure 6. Schematic of the code used to find leak rate based on thermomechanical simulation.
where L is the original length of the element, k c is the crack stiffness parameter, p is the internal pressure, x c is the location of the crack and u 1 and u 2 are the displacements at nodes 1 and 2, respectively. The magnitude of the difference between u c and u C c provides an expression for the crack opening displacement, that is, For the thermal case, there can be different temperatures on either side of the crack, the heat transfer is defined by the convection law where n 1;2 are the normals on faces 1 and 2, respectively. Balancing heat flux gives Given the linear nature of this problem, the single elements agree exactly with the analytical models. This is shown in Figures 7 and 8. The 1D study indicates that in principle, an XFEM element can model the thermal and mechanical effects of a leaking fluid within an element. The traction separation law, which is introduced to keep the element held together, is not necessary for a full 2D model. This is because the surrounding material holds the structure together. The temperature model used the Heaviside function, which meant that a jump in crack temperature could occur at the crack. However, if the heat transfer coefficient is sufficiently large, then the temperature on each side of the crack would be the same. For the case of a leaking fluid through a crack, the heat transfer coefficient in most practical situations  will be large enough for the crack temperatures to attain the same temperature. Therefore, it would be more appropriate to use an enrichment that can account for a discontinuity in the derivative of temperature. Both the Heaviside and Moës' enrichments are investigated in the following section for the temperature model.

Two-dimensional case
The mesh used for the 2D model is shown in Figure 9. Enriched elements are those that are cut by the crack and are within a certain radius of the crack tip. These elements are enriched with the Heaviside and crack tip functions, respectively. A visualisation of an angled crack is given in Figure 10. Analytical models (see Appendix) show that in reality, the heat flux singularity is logarithmic, whereas the selected enrichment function gives singular behaviour such as 1= p r. As will be seen in this paper, the p r is a pragmatic choice for the approximation of thermal fields around crack tips. This is because differentiation is more straightforward and high accuracy is achieved.

Convergence
In order to assess the accuracy and convergence performance of the element in a full finite element model, analytical models were derived and energy and temperature error norms calculated. The analytical expression for a prescribed temperature on the crack embedded in an infinite 2D domain is 692 P. GILL AND K. DAVEY Figure 9. Background mesh, squares denote the crack tip nodes and circles are the crack front nodes. A geometric enrichment is used here, where a radius is specified around the crack tip. All elements that fall within this radius are enriched. where is the contour of the crack, s is a parameter, which traverses the length of the crack, a is the crack length and .x; y/ is a point on the plate. The H 1 and L 2 norms are used to show the relative error in the energy and temperature, respectively, where T h is the numerical approximation and T is the exact solution, that is,

693
where the denominator terms are included for normalisation purposes. The singularity in heat flux is obtained in the approximation by choosing a r 1=2 enrichment function for the temperature field. On differentiation, this gives a r 1=2 , which is singular. In order to get a ln r singularity, an enrichment function such as r ln r r would be required. The r 1=2 is chosen as this is straightforward to work with when differentiating and integrating.
For the displacement, an analytical expression is derived on the basis of a method of complex potentials. The pressure term is included via the principle of superposition. The displacements are as follows: where is the complex potential according to Muskhelishvili [31], is a Lamé coefficient and Ä is the Kosolov constant. The variables x 1 , x 2 and u 1 , u 2 represent the coordinates and displacements in 1 and 2 directions. So the COD at the midpoint of the crack is where p and 0 are the pressure and far field stress, respectively. This is then compared with COD obtained from the finite element solution, which is given by the expression where x 0 is the midpoint of the crack. This is because standard shape functions are equal across the crack, and enriched-shape functions are˙1, so taking the difference gives a factor of two. There is an option to modify this crack opening displacement due to plasticity effects, which are known to increase the crack opening area considerably [32]. Currently, the code uses a relationship between elastic and elastic-plastic stress intensity factor .K/ from the R6 failure assessment diagram and the fact that COD is proportional to K 2 . The results of the convergence study are shown in Figures 11-13. Figure 11 shows that the error in temperature is less than 0.001% for a mesh size smaller than 0.01, and the slope of the line is 1.7 when the weak enrichment is used. The mesh size is defined to be the width of an individual element divided by the width of the plate. When a strong enrichment is used, the slope of the line is 1. The optimal convergence rate for a smooth function in standard finite elements is O.h 2 / if the L 2 norm is being used. Therefore, when the correct enrichment is     used, the convergence is getting close to optimal. This comes about because a fixed area or geometric enrichment is used. By utilising such an enrichment, optimal convergence can be achieved; however, this is at the expense of a huge increase in condition number. Ways to avoid such increases in condition number are discussed in [19]. Implementing this method on a larger scale would require more consideration of the condition number to improve solver performance; however, only the accuracy was of interest for this study. The energy error is less than 5% for a mesh size smaller than 0.01 and the slope of the line is 0.9 when weak enrichment is used. This is close to 1, which is the optimal convergence rate of the energy error norm. For a cracked domain, the optimal convergence of standard finite elements is O. p h/ [33]. When the strong enrichment is used for this problem, the slope of the convergence curve is 0.5, so the convergence is no better than standard finite elements. The number of integration points for standard and crack face elements was 3 and for the crack tip elements was 5.

Stress and heat flux
The stress, yy , after application of 155 bar of pressure on the crack faces is given in Figure 14, where the upper and lower surfaces are fixed in the y direction. The stress pattern is very similar to what is observed if a far field stress was applied and the crack faces are traction free. Heat flux and temperature are shown for an angled crack in Figures 15-17. A fluid temperature of 295 ı was taken, and the upper and lower surfaces were fixed at 290 ı . The singular behaviour in Figure 16 can be clearly seen at the crack tips. Figure 18 is a plot of heat flux vectors, this demonstrates how heat flows in a plate containing an angled crack with heat transfer on the faces.

Temperature using all enrichment functions
Using the approximation given in Equation 11, it is possible to exploit the properties of all the enrichment functions so that no assumption is made about the continuity of temperature or heat flux at the crack. A temperature gradient was imposed across a 1 m 1 m plate by fixing temperatures of 292 ı and 290 ı on the upper and lower surfaces. A horizontal crack of length 0.2 m was centred about the origin and the bulk fluid temperature was 295 ı . The heat transfer coefficient was varied between 0 and 1000 W/m 2 K. The results for temperature along the y-axis are shown in Figure 19.

Leak rates
Leakage rates are calculated for water which is slightly subcooled at high pressure. Thermomechanical simulations are performed using approximate values for standard operating conditions of a pressurised water reactor, that is, p 0 D 155 bar, T D 300 ı . Various crack lengths are simulated at three different fluid temperatures of 290 ı , 295 ı and 300 ı . Approximate material properties of stainless steel are taken with E D 200 GPa, D 0:3,˛D 1:282 10 5 and k D 20 W/mK. The dimensions of the plate are 1 m 1 m with a thickness of 25 mm, see Figure 9. Displacements of 0.1 mm and 0.1 mm are imposed on the upper and lower sides of the plate, and the temperatures were fixed at 290 ı , which is also taken as the reference temperature in the thermal strain calculation. The effect of bulk fluid temperature on leak rate for a horizontal crack is shown in Figure 20. Calculations were also carried out for advanced gas cooled reactor (AGR) conditions with a crack of length 0.2 m. Carbon dioxide is the fluid in this case and the temperatures ranged from 800 to 820 K, with an internal pressure was 40 bar. The results for the AGR conditions can be seen in Figure 21.

SUMMARY AND CONCLUSIONS
Presented here is an implementation of XFEM using boundary conditions that are dependent on a fluid model. The heat flux is obtained from Newton's law of cooling, where the heat transfer coefficient is obtained from empirical formulae. This provides a coupling of the fluid to the structure, which is thermomechanically modelled, thereby providing a means of studying changes in crack opening area due to crack wall heating. In order to test the accuracy and hence suitability of the model, convergence studies were undertaken for a 2D test case. Provided that the correct enrichment function is chosen, the convergence rates proved to be optimum for the thermal case. For the mechanical model, the crack opening area was chosen as the parameter on which to base convergence on. When the relative mesh size was about 0.01, the relative error in crack opening displacement was less than 1%.
A more general temperature model was then presented, which made no restrictions on the continuity of temperature or heat flux at the crack. This indicated that when using an approximation containing both the weak and strong discontinuity enrichment functions, situations ranging from adiabatic to high transfer could be modelled with the same model.
Finally, two test cases to demonstrate the features of the new finite element tool were undertaken. The first test case was based on conditions experienced in nuclear reactor primary circuits, specifically those seen in pressurised water reactors. As shown in Figure 20, there is a 30% reduction in leak rate when the fluid is 10 ı hotter than the structure. The next test case is for conditions in an AGR. Figure 21 shows that there is a 40% reduction in leak rate when the fluid is 20 ı hotter than the structure. These results proved the hypothesis that heat flux due to a leaking thermofluid has a noticeable closure effect on the crack opening and ultimately reduces leak rate.

FURTHER WORK
An extension of this work would be to implement this new method in a finite element code to enable 3D simulations. This is work in progress with the open source software Code Aster. Once the thermomechanical XFEM framework is in place, it will be possible to carry out tests to assess the viability for using this in Leak-before-Break assessments. Further to this, there is still much scope to improve thermal hydraulic codes used in this method. Two-phase flow is an area of active research [34,35], particularly for nuclear applications, and the current state of the art could be phased into the scheme presented in this paper.

A.1. Temperature
The problem involves an infinite 2D plate which contains a jump in heat flux or a specified temperature along a line. The temperature field is governed by the Laplacian (r 2 T D 0), or the screened Poisson equation(r 2 T 2 T D 0) if there is heat transfer from the faces of the plate. The jump in heat flux condition is defined to be For the specified temperature, and for a jump in temperature, where s parametrises the crack contour. The solution is obtained by applying Green's third identity to obtain the fundamental solution for the Laplacian, then using the single layer potential theory. The fundamental solution for the Laplacian in 2D is where x 0 is a point along the crack contour and x is any point. Green's identity is x 0 @w .x; x 0 / @n d I w x; x 0 @ .x 0 / @n d where n is the outward pointing normal to the boundary . Application of r 2 D 0 and r 2 w D ı .jx x 0 j/, the temperature can be expressed as .x/ D I .x/ @w .x; x 0 / @n d I w x; x 0 @ .x/ @n d (A. 6) considering the condition D N T D constant on . The density of dipoles is 0, @w @n D 0, so the problem is reduced to finding a function for the density of sources p .x 0 / D @ @n .
x 0 / such that the following integral equation is satisfied The crack contour is a line from a to a on the x-axis. Finally, there is the situation where there is a jump in crack face temperature: Take h.s/ D OEOET , so the density of sources @ @n D 0, and the density of dipoles is @w @n D .x 0 / D 2OEjT j. Therefore, the temperature is the solution of equation

A.2. Displacement
The problem involves an infinite plate under membrane loading with a central crack, which has a pressure acting perpendicular to the crack faces.
One can think of a crack in a plate as the superposition of a plate without a crack under membrane loading, plus the contribution of a compressive stress ( 0 ) along the crack, enforcing the traction free condition. In addition to this, there is the pressure (p) within the crack that acts in the same direction of the compression stress. The stresses and displacements in the plate can be expressed in terms of the complex potentials and