A novel thermo-mechanical coupling approach for thermal fracturing of rocks in the three-dimensional FDEM

Abstract

This paper presents a new three-dimensional thermo-mechanical (TM) coupling approach for thermal fracturing of rocks in the finite–discrete element method (FDEM). The linear thermal expansion formula is implemented in the context of FDEM according to the concept of the multiplicative split of the deformation gradient. The presented TM formulation is derived in the geo-mechanical solver, enabling thermal expansion and thermally induced fracturing. This TM approach is validated against analytical solutions of the Cauchy stress, thermal expansion and stress distribution. Additionally, the thermal load on the previously validated configurations is increased and the resulting fracture initiation and propagation are observed. Finally, simulation results of the cracking of a reinforced concrete structure under thermal stress are compared to experimental results. Results are in excellent agreement.

1 Introduction

Thermally induced deformations and fracturing are significant concerns in many engineering fields such as engine design, nuclear fission reactors, geothermal energy, hydrocarbon production and radioactive waste disposal. When a material is exposed to a temperature gradient, it will deform and cracks may appear, the prediction of those cracks is key in making robust, efficient and safe designs. In geothermal energy exploitation, the injection cycles of cool water in the hot reservoir are likely to induce fractures [37], and in hydrocarbon production, thermal shocks are found to enhance the hydraulic fracturing efficiency [6]. In such applications, the fractures increase the permeability of the reservoirs and have a positive influence on production. On the contrary, in the context of radioactive waste disposal fractures are undesirable as initiating new cracks in the rock is potentially creating new flow pathways for radionuclides to be transported into the biosphere. As high temperatures are expected in the vicinity of a geological disposal facility (GDF), predicting thermally induced deformations and failure is of importance for safe disposal of the radioactive waste. In radioactive repository applications, the major heat-driven process is the thermal expansion of the rock. Thermal expansion is a strongly coupled process because it requires the mechanical equations to be modified. A typical granite has a coefficient of thermal expansion of $8 \times 10^{-6}/^\circ \hbox {C}$. For a temperature rise of $50\,^\circ \hbox {C}$, this is equivalent to a $0.015 \%$ volume change. Although a ground surface uplift of up to a few tens of centimetres is expected above the GDF [13, 31], large thermal stresses will be localised near the deposition holes. This is because significant thermal gradients will only be found within a few tens of metres from the heat source [13]. In the vicinity of a deposition hole, thermal stresses will influence the porosity of the media and the fracture apertures. For a horizontal deposition hole, the generated horizontal thermal stress is expected to close up vertical fractures and thus to reduce the overall hydraulic transmissivity of the rock [3, 24, 25].

The thermo-mechanical coupling is a common feature among numerical simulators but thermally induced fracturing can only be performed with methods that represent fractures explicitly. Some research groups developed TM models for thermal fracturing based on: the extended finite element method (X-FEM) [5, 17, 38], the boundary element method (BEM) [8, 23, 26] and the discontinuous deformation analysis (DDA) [14]. However, these methods do not consider the effect of the thermal contacts in a multibody system, e.g. contact between two fracture walls, or heat resistance due to contact of two solid bodies, etc. In the last two decades, thermo-mechanical coupling and thermal fracturing capabilities were increasingly developed within the particle-based method, a sub-class of distinct element methods (explicit DEM), originally developed to model the behaviour of granular materials [4]. In particle-based methods, the discrete elements are typically discs in 2D and spheres in 3D. They are considered rigid; thus, only the element interaction needs to be solved. As a result, this method benefits from a low computational cost and large number of particles can be used. This makes particle-based methods advantageous in micro-fracturing where the particle approach is a good approximation of the microstructure of rocks [16]. To model fracture mechanics in particle-based methods, Potyondy et al. [21] developed a bonded particle model (BPM) where intact materials are represented by an assembly of bonded particles which may break under certain stress intensity. Within this framework, Wanne et al. [30] proposed an extension of the BPM with thermo-elastic bonds for the thermal fracturing of granites. To resolve heat conduction in the particle system, Feng et al. [7] presented a 2D model where contacting circular particles share heat flux bonds with their neighbours; such models also exist in 3D using spherical particles [22]. When combining the thermal and bonded particle type of models, thermal fracturing was achieved by Xia et al. [32, 33] for circular particles and by Andre et al. [2] for spherical particles. Among particle-based methods, we also find the peridynamic method [27] which was the principal advantage that its governing equations stay valid over discontinuities. Within the peridynamics framework, Wang et al. [29] recently developed a thermo-mechanical model for the thermal cracking of rocks. However, particle-based methods rely on the assumption that the temperature is uniform within each particle which is not valid in some configurations. To compute the heat distribution within each particle, particle-based methods rely, when available, on extra analytical solutions.

Recently, Yan and Zheng [35, 36] developed a coupled thermo-mechanical model based on the FDEM with demonstration of thermal fracturing in both two and three dimensions. They demonstrated that simulation results were in good agreement with the analytical solutions and experimental results. However, the expression of the thermal stress used in their paper (equation 50) belongs to the infinitesimal strain theory in which the displacements of the solid body are assumed to be much smaller than any relevant dimension of the body. In order to deal with large strains and/or rotations (as often arise in FDEM application fields) which invalidate assumptions inherent in infinitesimal strain theory, this paper presents a novel thermo-mechanical coupling formulation in three dimensions for the large strain FDEM with thermally induced fracturing. This model is also consistent with a considerable body of FDEM theory [12, 28], i.e. the multiplicative decomposition of the deformation gradient tensor, which specifically adopts finite strain theory based on the left Cauchy–Green tensor for a consistent Cauchy stress calculation that takes into account the dissipative part of the thermal stress. The thermal model introduced in [15] is coupled with the existing finite strain model of the FDEM approach [34] and its fracture model [9]. The thermo-elastic theory for large deformations is presented first, then the thermal stress is validated against analytical solutions and finally a three-dimensional validation of thermally induced fracturing is presented.

2 Thermo-elastic theory for finite strain FDEM

2.1 Deformation gradient F

We define $\mathbf{X}$ as the vector of the initial positions in the element and $\mathbf{x}$ the vector of current positions [34], and we write:

$$\begin{aligned} \mathbf{X }= & {} \begin{bmatrix} X\\ Y\\ Z\\ \end{bmatrix} =\mathbf{X} _{on}{} \mathbf{N} \end{aligned}$$

(1)

$$\begin{aligned} \mathbf{x}= & {} \begin{bmatrix} x\\ y\\ z\\ \end{bmatrix} =\mathbf{x} _{cn}{} \mathbf{N} . \end{aligned}$$

(2)

$\mathbf{N} =[N_1 \ N_2 \ ...\ N_{n_e}]$ is the shape function or interpolation function, $n_e$ is the number of nodes in the element and $\mathbf{X} _{on}$ and $\mathbf{x} _{cn}$ are, respectively, the original and current position arrays

$$\begin{aligned} \mathbf{X} _{on}= & {} \begin{bmatrix} X_1 &{}\quad X_2 &{}\quad \ldots &{}\quad X_{n_e}\\ Y_1 &{}\quad Y_2 &{}\quad \ldots &{}\quad Y_{n_e}\\ Z_1 &{}\quad Z_2 &{}\quad \ldots &{}\quad Z_{n_e}\\ \end{bmatrix} \end{aligned}$$

(3)

$$\begin{aligned} \mathbf{x} _{cn}= & {} \begin{bmatrix} x_1 &{}\quad x_2 &{}\quad \ldots &{}\quad x_{n_e}\\ y_1 &{}\quad y_2 &{}\quad \ldots &{}\quad y_{n_e}\\ z_1 &{}\quad z_2 &{}\quad \ldots &{}\quad z_{n_e}\\ \end{bmatrix}. \end{aligned}$$

(4)

The deformation gradient tensor links the current position to the initial element positions:

$$\begin{aligned} \mathbf{F} =\frac{\partial \mathbf{x} }{\partial \mathbf{X} }= \mathbf{x} _{cn} \frac{\partial \mathbf{N} }{\partial \mathbf{X} }. \end{aligned}$$

(5)

The determinant of the matrix F is called the Jacobian

$$\begin{aligned} J = {\hbox {det}}(\mathbf{F} ). \end{aligned}$$

(6)

From here, we define the left Cauchy–Green tensor:

$$\begin{aligned} \mathbf{B} =\mathbf{F} {} \mathbf{F} ^\top . \end{aligned}$$

(7)

2.2 Velocity gradient L

We introduce the velocity vector as

$$\begin{aligned} \mathbf{v} = \begin{bmatrix} v_x\\ v_y\\ v_z\\ \end{bmatrix} = \mathbf{v} _{cn} \mathbf{N} . \end{aligned}$$

(8)

With $\mathbf{v} _{cn}$,

$$\begin{aligned} \mathbf{v} = \begin{bmatrix} v_{x,1} &{}\quad v_{x,2} &{}\quad \ldots &{}\quad v_{x, n_e}\\ v_{y,1} &{}\quad v_{y,2} &{}\quad \ldots &{}\quad v_{y, n_e}\\ v_{z,1} &{}\quad v_{z,2} &{}\quad \ldots &{}\quad v_{z, n_e}\\ \end{bmatrix} = \mathbf{v} _{cn} \mathbf{N} . \end{aligned}$$

(9)

The velocity gradient is

$$\begin{aligned} \mathbf{L} = \nabla \mathbf{v} = \frac{\partial \mathbf{v} }{\partial \mathbf{x} } =\mathbf{v} _{cn} \frac{\partial \mathbf{N} }{\partial \mathbf{x} }. \end{aligned}$$

(10)

We can also write the velocity gradient in terms of the deformation gradient with

$$\begin{aligned} \dot{\mathbf{F }}= & {} \frac{\partial \mathbf{F} }{\partial t} = \frac{\partial }{\partial t} \left( \frac{\partial \mathbf{x} }{\partial \mathbf{X} } \right) = \frac{\partial }{\partial \mathbf{X} } \left( \frac{\partial \mathbf{x} }{\partial t} \right) = \frac{\partial \mathbf{v} }{\partial \mathbf{X} } \end{aligned}$$

(11)

$$\begin{aligned} \dot{\mathbf{F }}= & {} \frac{\partial \mathbf{v} }{\partial \mathbf{X} } = \frac{\partial \mathbf{v} }{\partial \mathbf{x} } \frac{\partial \mathbf{x} }{\partial \mathbf{X} }. \end{aligned}$$

(12)

Hence,

$$\begin{aligned} \dot{\mathbf{F }} = \mathbf{L} . \mathbf{F} . \end{aligned}$$

(13)

And

$$\begin{aligned} \mathbf{L} = \dot{\mathbf{F }} . \mathbf{F} ^{\text {-}1}. \end{aligned}$$

(14)

Finally, we write the rate of deformation matrix as:

$$\begin{aligned} \mathbf{D} =\frac{1}{2} (\mathbf{L} + \mathbf{L} ^\top ). \end{aligned}$$

(15)

2.3 Thermal expansion

The law of heat conduction, also known as Fourier’s law, is implemented to simulate the temperature distribution within the solid (see detailed implementation in [15]). The linear thermal expansion coefficient corresponds to the one-dimensional change in length for a given temperature change and is noted $\alpha$. According to the concept of multiplicative split of the deformation gradient [28], the deformation gradient may be decomposed into a thermal component $\mathbf{F} _T$ and an elastic $\mathbf{F} _e$ component:

$$\begin{aligned} \mathbf{F} =\mathbf{F} _{T}{} \mathbf{F} _{e} \end{aligned}$$

(16)

with

$$\begin{aligned} \mathbf{F} _{T}=\varUpsilon (T) \mathbf{I} . \end{aligned}$$

(17)

For an isotropic material, $\varUpsilon =\varUpsilon (T)$ is the ratio linking the variation of temperature to a deformation in a principal direction and $\mathbf{I}$ is the identity matrix. We also deduct from the expression of the velocity gradient (14) that

$$\begin{aligned} \mathbf{L} _T= \frac{\dot{\varUpsilon }}{\varUpsilon } \mathbf{I} = \frac{1}{\varUpsilon } \frac{\partial {\varUpsilon }}{\partial t} \mathbf{I} = \frac{1}{\varUpsilon } \frac{\partial {\varUpsilon }}{\partial T} \frac{{\partial T}}{\partial t} \mathbf{I} . \end{aligned}$$

(18)

Consider now the volume $V_T$ resulting from a thermal expansion event of an initial volume $V_0$ from a temperature $T_0$ to T

$$\begin{aligned} dV_T = det (\mathbf{F} _T) dV_0. \end{aligned}$$

(19)

According to the relationship $\dot{J} = J \ . \ tr(\mathbf{L} )$ [20], the time derivative of the above expression is:

$$\begin{aligned} \frac{d}{dt}(dV_T) = tr(\mathbf{L} _T) dV_0. \end{aligned}$$

(20)

Integrating equation 18 yields

$$\begin{aligned} \frac{d}{dt}(dV_T) = \frac{3}{\varUpsilon } \frac{\partial {\varUpsilon }}{\partial T} \frac{{\partial T}}{\partial t} dV_0. \end{aligned}$$

(21)

The traditional temperature-dependent thermal expansion coefficient is:

$$\begin{aligned} \alpha (T)= \frac{1}{\varUpsilon } \frac{\partial \varUpsilon }{\partial T}. \end{aligned}$$

(22)

Hence,

$$\begin{aligned} \frac{{\hbox {d}}}{{\hbox {d}}t}(dV_T) = 3\alpha (T) \frac{\partial T}{\partial t} dV_0. \end{aligned}$$

(23)

Integrating equation 22 over the temperature change gives

$$\begin{aligned} \varUpsilon (T) = {\hbox {exp}} \left( \int ^T_{T_0}\alpha (T) dT \right) . \end{aligned}$$

(24)

If the thermal coefficient is considered independent of temperature i.e. $\alpha (T)= \alpha$ and $\alpha |T -T_0 |<<1$, we can admit the following approximation [28]

$$\begin{aligned} \varUpsilon (T) = 1 + \alpha \Delta T, \quad \Delta T = T - T_0. \end{aligned}$$

(25)

Fig. 1

Cohesive zone model. (a) Schematic model of the transition between elastic, plastic and broken zones, (b) representation of the cohesive zone model in FDEM [11]

Fig. 2

Insertion of joint elements in 3D: a continuous FEM formulation, b discontinuous FEM formulation. From [10]

2.4 Thermo-elastic tensors

Let us now rewrite the previous tensors in terms of the thermo-elastic decomposition. We start with developing equation 25 into equation 16

$$\begin{aligned} \mathbf{F} =(1 + \alpha \Delta T) \mathbf{F} _e. \end{aligned}$$

(26)

The Jacobian becomes

$$\begin{aligned} J = {\hbox {det}}(\mathbf{F} ) = {\hbox {det}}(\mathbf{F} _T) \ . \ {\hbox {det}}(\mathbf{F} _e) = (1 + \alpha \Delta T )^3 {\hbox {det}}(\mathbf{F} _e). \end{aligned}$$

(27)

We may also write

$$\begin{aligned} J = J_{T} J_e. \end{aligned}$$

(28)

The velocity gradient tensor becomes

Table 1 Mechanical and thermal parameters for validation, material properties of concrete taken from [36]

Fig. 3

Finite element model of the cube with two different mesh sizes

$$\begin{aligned} \mathbf{L} = \dot{\mathbf{F }} . \mathbf{F} ^{-1} = \big [ \mathbf{F} _T \dot{\mathbf{F }}_e + \mathbf{F} _e \dot{\mathbf{F }}_T \big ] \big [\mathbf{F} _T \mathbf{F} _e \big ]^{-1} = \dot{\mathbf{F }}_T \mathbf{F} _T^{\text {-}1} + \dot{\mathbf{F }}_e \mathbf{F} _e^{-1}. \end{aligned}$$

(29)

Note that there is commutativity of the matrix product in the above development because both $\mathbf{F} _T$ and $\dot{\mathbf{F }}_T$ are a product between a scalar and the identity matrix. We obtain

$$\begin{aligned} \mathbf{L} = \mathbf{L} _{T} + \mathbf{L} _e. \end{aligned}$$

(30)

The left Cauchy–Green tensor becomes

$$\begin{aligned} \mathbf{B} =\mathbf{F} {} \mathbf{F} ^\top = \mathbf{F} _T\mathbf{F} _e (\mathbf{F} _T\mathbf{F} _e)^\top . \end{aligned}$$

(31)

And

$$\begin{aligned} \mathbf{B} = \mathbf{B} _{T} \mathbf{B} _e. \end{aligned}$$

(32)

Finally, the rate of deformation matrix becomes

$$\begin{aligned} \mathbf{D} =\frac{1}{2} (\mathbf{L} + \mathbf{L} ^\top ) = \mathbf{D} _T + \mathbf{D} _e. \end{aligned}$$

(33)

Fig. 4

Temperature versus volume change of a 1-m cube for two different mesh sizes

2.5 Cauchy stress

The Cauchy stress is defined as force per unit area and is necessary to solve the momentum equation. As given by [34]:

$$\begin{aligned} \mathbf{C} =\frac{\mu }{J}(\mathbf{B} -\mathbf{I} ) + \frac{\lambda }{J}(\ln J)\mathbf{I} + \mathbf{C} _D \end{aligned}$$

(34)

with $\lambda$ and $\mu$ being the Lamé coefficients and $\mathbf{C} _D$ the dissipative part of the stress defined as:

$$\begin{aligned} \mathbf{C} _D=2\eta \mathbf{D} \end{aligned}$$

(35)

with $\eta$ being the viscosity of the material considered.

Fig. 5

3D model of the hollow cylinder

2.6 Analysis of stress response

Now considering only internal forces caused by the thermal expansion, we have

$$\begin{aligned} \mathbf{C }_T=\frac{\mu }{J_T}(\mathbf{B} _T-\mathbf{I} ) + \frac{\lambda }{J_T}(\ln J_T)\mathbf{I} + \mathbf{C} _{D_T} \end{aligned}$$

(36)

with $J_T$ being the volume change induced by thermal expansion:

$$\begin{aligned} J_T = det(\mathbf{F} _T)= (1 + \alpha \Delta T)^3 \end{aligned}$$

(37)

which yields

$$\begin{aligned} \begin{aligned} \mathbf{C} _T&=\left[ \frac{\mu }{(1 + \alpha \Delta T)^3}((1 + \alpha \Delta T)^2-1)\right. \\&\quad + \frac{\lambda }{(1 + \alpha \Delta T)^3}\ln (1 + \alpha \Delta T)^3 \\&\quad \left. + \frac{2\eta \alpha }{1 + \alpha \Delta T} \frac{\partial {T}}{\partial t} \right] \mathbf{I} . \end{aligned} \end{aligned}$$

(38)

We can use the two above expressions to perform a validation test, controlling the volumetric expansion with $J_T$ and the associated stress with $\mathbf{C} _T$.

2.7 Fracture model

The model used to simulate thermally induced fracturing in this paper was developed by Guo el al. [10], a novelty in the three-dimensional FDEM. In this section, this fracture model is introduced briefly. For the nonlinear fracturing process, a cohesive zone model [18] handles the transitional elasto-plastic behaviour of joint elements (Fig. 1).

Table 2 Mechanical and thermal parameters for validation

When using the fracture model joint, elements are inserted between triangular (2D) or tetrahedral finite elements (3D) pre-simulation as shown in Fig. 2.

The joint elements are four-noded in 2D and six-noded in 3D. They may be cohesive (i.e. unbroken) to represent the intact rock or broken to represent new or pre-existing fractures. In tension, joint elements act as springs until they are broken which is analogous to the approach used in the explicit DEM codes, while in compression, joint elements will not penetrate each other significantly because the contact force is calculated with the penalty function method [9].

3 Validation work

The TM model has been implemented in the framework of a solid modelling code: Solidity (Solidityproject.com) which is an open-source, multi-purpose explicit mechanical solver. Solidity can solve the mechanical equations with a hybrid continuum (FEM)–discontinuum (DEM)—FDEM approach; this reduces to the FEM when dealing with a single solid; FEM simulations will be referred as solving a ‘continuum’. Solidity can also solve with a fully discontinuum explicit DEM (discrete element method) approach, i.e. the fracture model introduced previously which will be referred as solving for a ‘discontinuum’. In this section, the thermo-mechanical coupling implemented in Solidity will be verified for both continuum and discontinuum problems.

3.1 Cauchy stress validation

To make sure that the thermo-mechanical coupling model presented in this paper is implemented correctly, the Cauchy stress response to a temperature change is verified. A single finite element is fixed in all directions and subjected to a temperature change. The Cauchy stress generated by thermal expansion at the element level is compared to the analytical result (Eq. 38). With the parameters given in Table 1 [36], the Cauchy stress in all three directions is 3.332 MPa and the value computed with Solidity for one element is 3.332 MPa. Note that we do not consider here the transient term of the Cauchy stress which is associated with the dissipative part of the stress and depends on how fast the temperature changes.

3.2 Thermal expansion validation

The thermal expansion of a material for a given temperature change is verified as follows. Consider a cubical solid with 1-metre edges meshed with two different element sizes, (a) 1 m and (b) 0.1 m (Fig. 3). A simulation is performed with the thermo-mechanical properties and temperature change listed in Table 1. The temperature is increased uniformly from $0 \ ^\circ \hbox {C}$ to $100 \ ^\circ \hbox {C}$ over $0.1\,\hbox {s}$ with an integration time step of $1.0\times 10^{-5} \, \hbox {s}$ for (a) and of $1.0\times 10^{-6} \, \hbox {s}$ for (b). For the given temperature increase, the $1\,{\hbox {m}}^3$ cube expands to reach a volume of $1*\mathbf{J} _T = (1 + \alpha \Delta T)^3 = 1.003003\,\hbox {m}^3$. This value is in accord with the simulation results for models (a) and (b), as shown in Fig. 4; results are in agreement with the analytical solution. Note that no significant thermal stress is generated in this configuration because the temperature is increased uniformly and the cube is free to expand.

Fig. 6

Steady-state temperature profile in the hollow cylinder

Fig. 7

Radial thermal stress $\sigma _{rr}$ in the hollow cylinder

Fig. 8

Transverse thermal stress $\sigma _{\theta \theta }$ in the hollow cylinder

Fig. 9

Transverse thermal stress $\sigma _{\theta \theta }$ in the hollow cylinder

3.3 Thermal stress validation: hollow cylinder

Now that the implementation of the TM model is verified; a validation of the thermal stress field is performed. A thermal gradient must be present in order to generate differential thermal stress in the material. Consider a hollow and thin cylinder as presented in Fig. 5, with an inner and outer radii, respectively, a and b and a height h. A Dirichlet boundary condition is applied on the inner boundary of the thin cylinder with a temperature $T_a$ and on the outer boundary with a temperature $T_b$. The faces normal to the Z direction are considered adiabatic. For this problem, the temperature and stress solution are given by [19]:

Fig. 10

Thermal fracturing of the hollow cylinder at different times (a–c). First row is the maximal principal stress ($\sigma _1$), second row is the crack pattern, and third row is the temperature field displayed on the continuum mesh. Note that the engineering sign conventional is adopted where tensile stress is positive and hot red colour indicates tensile stress is generated at the outer surface in (a) top row. (Color figure online)

Fig. 11

Thermal cracking due to differential thermal expansion of a reinforcement embedded in concrete [1]

$$\begin{aligned}&T(r) = T_a + (T_b - T_a)\frac{\ln {a}}{\ln {b}} \nonumber \\&\sigma _{rr} = \frac{\alpha E}{2}(T_b-T_a)\left[ -\frac{\ln {(r/a)}}{\ln {(b/a)}}+\left( 1 - \frac{a^2}{b^2}\right) \frac{b^2}{b^2-a^2}\right] \nonumber \\&\sigma _{\theta \theta } = \frac{\alpha E}{2}(T_b-T_a)\left[ -\frac{1+\ln {(r/a)}}{\ln {(b/a)}}+\left( 1 + \frac{a^2}{b^2}\right) \frac{b^2}{b^2-a^2}\right] .\nonumber \\ \end{aligned}$$

(39)

To validate the above analytical solution, the hollow cylinder is meshed with 0.005m finite elements, as shown in Fig. 5. Note that when employing the hybrid FDEM on a single body, the discontinuum part is not recruited and the solver only uses the FEM. Thereby, two simulations are performed, one with the FEM solver (referred as ‘continuum’) and other with the explicit DEM fracture model (referred as ‘discontinuum’). Simulations are performed with an integration time step of $\ 5.0\times 10^{-9} s$, and the total simulation time is of $\ 5.0\times 10^{-3} s$. Material properties and boundary conditions are listed in Table 2. Temperature, radial and transverse thermal stress results are, respectively, presented in Figs. 6, 7 and 8, and good agreement is observed between simulation results and analytical solution.

Table 3 Mechanical and thermal parameters for validation

Fig. 12

3D model of the concrete reinforcement

3.4 Thermal fracturing validation: hollow cylinder

Now that the stress field has been validated against Equation 39, and this analytical solution is used again to dimension a simulation that will generate a fracture. For a tensile fracture to appear, the thermal stress in the transverse direction ($\sigma _{\theta \theta }$) must exceed the tensile strength of the material. The thermal expansion coefficient of the material is increased sixfold; the maximal $\sigma _{\theta \theta }$ that will be generated in the cylinder is now higher than the tensile strength of the material, as highlighted in Fig. 9; fractures are expected to initiate on the outer boundary of the cylinder. Note that in Fig. 9, only the continuum stress profile is compared to the analytical solution because when fractures are generated in the discontinuum simulation, the stress field is perturbed. Graphical results of the discontinuum simulation are presented in Fig. 10. The first two cracks appear on the outer boundary of the cylinder (Fig. 10b) and propagate the centre of the cylinder (Fig. 10c). Then, several cracks initiate in the inner boundary of the cylinder, three of them propagate outward (Fig. 10d) and two of them cut through to the outer boundary of the cylinder (Fig. 10d).

4 Application of thermal cracking for concentric cylinders under thermal stress

The cracking of reinforced concrete structures under thermal stress was investigated numerically and experimentally by [1]. This work has been the basis for validation of some of the numerical methods presented in “Introduction” of this paper [29, 36]. Reinforcements in concrete are often made of steel which has a higher thermal expansion coefficient. Upon heating, the steel exerts a pressure on the concrete cover and induces fracturing, as shown in Fig. 11.

Fig. 13

Radial stress $\sigma _{rr}$ in the outer cylinder

A solution of the radial stress is given in [1]:

$$\begin{aligned} \begin{aligned}&\sigma _{\theta \theta }= \frac{a^2(b^2 + r^2)}{r^2(b^2 - a^2)}p\\&p = \frac{(\alpha _a - \alpha _b)\Delta T E_a}{E_a/E_b(\beta + \nu _b)+(1 - \nu _a)}\\&\beta = \frac{b^2 + a^2}{b^2 - a^2}. \end{aligned} \end{aligned}$$

(40)

First, this expression is verified with the continuum and discontinuum models using coefficients of thermal expansion ten times smaller than the experiment: $\alpha _a = 7.2\times 10^{-7}/^\circ \hbox {C}$ and $\alpha _b = 2.2\times 10^{-6} /^\circ \hbox {C}$. Then, the discontinuum model is employed with the thermal expansion coefficients of the experiment (Table 3) and the fracture pattern is compared to the experiment. The model presented in Fig. 12 is meshed with 0.005 m finite elements; for all simulation, an integration time step of $\ 5.0\times 10^{-9}$ s is used. The total simulation time is of $\ 5.5\times 10^{-3}$ s, and temperature is increased uniformly from $0 ^\circ \hbox {C}$ to $100 ^\circ \hbox {C}$ over a time of $5.0\times 10^{-3}$ s. Results are presented in Fig. 13 for $\sigma _{rr}$ and in Fig. 14 for $\sigma _{\theta \theta }$. The stress profiles show that the tangential stress is maximal for $r = a$. This is where fracture is expected to initiate when using the experimental thermal expansion coefficients $\alpha _a = 7.2\times 10^{-6} /^\circ \hbox {C}$ and $\alpha _b = 2.2\times 10^{-5} /^\circ \hbox {C}$. Simulation results of the thermal fracturing simulation performed with the experimental thermal expansion coefficients are shown in Fig. 15. In summary, Fig. 14 shows good agreement between the stress profile and the analytical solution. Moreover, the final fracture pattern (Fig. 15d) is consistent with experimental results (Fig. 11a) in the sense that radial cracks propagate in the concrete, from the interface with steel and up to the outer boundary. In the light of the several successful validations performed in this paper, the thermo-mechanical strong coupling is considered complete.

Fig. 14

Tangential stress $\sigma _{\theta \theta }$ in the outer cylinder

Fig. 15

Thermal fracturing of the reinforced concrete cylinder at different times (a–c). The first row is the maximum principal stress ($\sigma _1$), and the second row is the crack pattern

5 Conclusions

The new three-dimensional TM model in the context of the FDEM is presented in this paper. The thermal expansion formula is implemented and successfully integrated in FDEM based on the multiplicative decomposition of deformation gradient. The first two validation tests: static Cauchy stress and deformation under thermal load, show that the model gives very accurate results and is not sensitive to mesh size. The third and fourth validation tests: distribution of thermal stress and thermal fracturing in a hollow cylinder, also show very good agreement with the analytical solutions. Finally, the implemented finite strain TM model is used for a thermal cracking application, i.e. the cracking of reinforced concrete structures under thermal stress. The results show good agreement between the stress profile and the analytical solution, and the final fracture pattern is consistent with experimental results. The TM method has the potential to tackle complex TM-coupled configurations in fractured rock media. However, it should be noted that enhanced accuracy is achieved at the expense of higher computational costs compared to other continuum- or discontinuum- only numerical methods.