Scaled boundary finite element modelling of thermally induced crack propagation

Numerical modelling of thermal shock processes is a very challenging task due to the rapid development of complex crack networks. In the context of linear elastic fracture mechanics, the scaled boundary finite element method (SBFEM) has gained popularity due to its semi‐analytical nature. Here, modelling of stress singularities can be achieved efficiently and accurately. It also facilitates the derivation of generalized polygon elements and thus provides great flexibility in meshing of complex geometries. In crack propagation modelling, such polygon meshes require minimal re‐meshing in the vicinity of the crack.


SBFEM for heat conduction
The scaled boundary finite element method (SBFEM) is a semi-analytical technique [1]. The boundary of a polygon element is discretized using 1D finite elements. The analytical solution is formulated in the radial direction. SB finite elements can have an arbitrary number of edges, which are meshed with higher-order line elements, c.f. Fig. 1. Inside the polygon, there is a point called the scaling center from where the whole boundary of the polygon domain is visible. The boundary of an element is scaled using the radial coordinate ξ (0 ≤ ξ ≤ 1). One of the main features of SBFEM is the ability to model open polygons (i.e. containing cracks). Open polygon elements are characterized by a built-in singularity due to the analytical solution in the radial direction.

Scaling Center
The scaled boundary coordinate transformation is given as equation (1), wherex,ŷ are the Cartesian coordinates, x, y are the coordinates on the boundary and x, y are the vectors of the nodal coordinates.
The governing equation of heat conduction reads (equation 2), where the parameters c ε , ρ, κ, p and θ correspond to specific heat, mass density, anisotropic heat conductivity matrix, source and temperature field, respectively, The semi-discretization of the temperature field is given as θ = N(η)θ(ξ). Assuming steady-state heat conduction and vanishing source terms p, applying the method of weighted residuals to the governing equation (2) yields the SBFE equation in temperature (3), where E 0 , E 1 and E 2 are the coefficient matrices. These matrices only depend on the material properties, geometry and the discretization of the boundary. Equation (3) can be solved analytically by transforming it to first-order ordinary differential equations (4) with 2n unknowns by introducing the nodal flux vector q, where n is the number of degrees of freedom.

SBFEM for thermal stress
The procedure explained above can be extended to linear elasticity. The SBFE equation in displacement (6) is where ξF σ (ξ), ξ 2 P σ (ξ) and ξP σ (ξ) are inhomogeneous terms corresponding to traction, body load and thermally-induced stress, respectively. The homogeneous problem yields the element stiffness matrixK = Ψ f Ψ −1 u , where Ψ u are displacement modes. For plane strain, the initial thermal stress is expressed as where D and α are the elasticity matrix and coefficient of thermal expansion, respectively. The power function solution (5) of the temperature field θ = θ d (η)ξ d is used to construct the particular solution for the inhomogeneous term P σ (ξ) [3]: whereB 1 (η),B 2 (η) result from expressing the differential operator of linear elasticity in SB coordinates and J is the Jacobian matrix. Using the SB coordinates, the stress vector can be expressed as Equation (9) with the stress modes Ψ σi The constants c i are obtained by evaluating the solution at the boundary, c = Ψ −1 u u b . The stress singularities can be modelled very accurately and efficiently by using only the singular modes (i.e.S n > −1). The normal transformation procedure from Cartesian coordinates to polar coordinates is applied. The singular stress term can be written as where L is a characteristic length so that the generalised stress intensity factors are independent of the system of units used.
The generalized stress intensity factors can be constructed directly from their definition [4]: . (11)

Numerical Example
A plate with an initial edge crack is loaded by a constant flux along the vertical direction. Zero flux is assumed at the crack faces and vertical boundaries. The left and right sides of the plate are fixed in horizontal direction, also the top left and right corners are fixed in vertical direction. The convergence of the normalized stress intensity factor at the initial step is shown in Fig. 3. The maximum circumferential stress criterion is used to determine the crack propagation path. Only local re-meshing is required in the vicinity of crack tip (Fig. 2).  Fig. 3: Normalized stress intensity factor. Reference solutions from [5].