Three-Dimensional Constitutive Viscoelastic Model for Isotropic Materials

1. Introduction

1.1. Content

The characterization of materials behaviour is of importance in many research fields. From an analytical point of view, the mechanical behaviour of materials could be described either by using empirical laws based on experimental observations or by using a framework to develop constitutive laws. In the latter case, the thermodynamic of irreversible processes could be used as the framework. These constitutive laws are generally dedicated to complex problems and are thus developed in a three-dimensional context. Materials behave in different ways under loading but, under specific conditions, they all will generally exhibit instantaneous and time-dependant deformations. Instantaneous deformation could be elastic, plastic and so forth, while time-dependant deformation generally refers to the viscosity of the material.

Moreover, as for the Navier-Stokes equation in fluid mechanics which is only valid for Newtonian fluids, constitutive laws are developed considering the nature of the materials. Hence, in this chapter, targeted materials are isotropic granular materials such as concrete, carbon block materials and ramming paste (Sørlie and Øye, 2010). Contrary to metals, this implies that the tensile behaviour is different from the compressive one. Models are also meant to represent various phenomena such as creep and relaxation, which are mainly those of interest in this chapter. The upcoming sections will give a brief overview of the mechanical deformations such as creep/relaxation that are involved in viscoelasticity. Some basic principles governing the development of constitutive law modelling within a thermodynamic framework will be discussed, as well.

1.2. Viscoelasticity

Viscoelasticity of materials could be expressed in many ways. In this chapter, it will be described as the ability of the material to deform elastically, viscously and/or a combination of those. The interaction between the elastic and viscous behaviours could be explained with the help of the rheology of materials, as discussed in section 3.4.

For solids, the elastic behaviour is related to the instantaneous deformation of a material and is expressed using Young’s modulus. This is generally measured by two different methods. The first one consists of calculating the slope of the stress-strain curve resulting from a uniaxial loading test as described by the one-dimensional (1D) Hooke’s law (1):

$$\begin{array}{ccc}E=\frac{\sigma }{\epsilon }=\frac{F/A_0}{\Delta L/L_0}=\frac{F{L}_0}{A_0\Delta L}& where& \Delta L=\left|L_f-L_0\right|\end{array}$$E1

where $E,\sigma ,\epsilon ,F,A_0,L_f$ and $L_0$are respectively the Young’s modulus, the engineering stress, the engineering strain, the maximum load applied, the initial surface area on which the force is applied, the final length and the initial length of the sample. In the case of hyper elasticity, the true stress and true strain should be used but won’t be discussed here. Both static and dynamic loading tests could be performed and thus leading to the static and dynamic Young’s modulus. The second method calculates the Young’s modulus using acoustic wave technique (Kinsler, 2000). A relation is established between the elasticity of the material, i.e., the Young’s modulus, and the speed of sound travelling through it, as described by Eq. (2):

$$E=c^2\rho \frac{(1+\nu )(1-2\nu )}{(1-\nu )}$$E2

where $c,\rho$and$\nu$are respectively the speed of sound in the material, the density and the Poisson’s ratio. For a cylindrical sample, the Poisson’s ratio is defined as the ration of the radial strain (perpendicular to the applied load) over the axial strain (in the direction of the applied load):

$$\nu =-\frac{{\epsilon }_{radial}}{{\epsilon }_{axial}}$$E3

Young’s modulus should therefore be calculated adequately with either Eq. (1) or (2) according to the context. Materials discussed in this chapter generally undergo static loading and thus elasticity of materials will be defined as the slope of the stress-strain curve resulting from a uniaxial static loading test, as defined in Eq. (1).

The viscous behaviour of material could be rather complex. As aforementioned, it may be described with time-dependant functions. Time-dependant deformations can be both reversible and/or irreversible but in either case dissipation of energy is involved. This point will be clarified in section 3. According to material rheology, a reversible mechanism means that the deformations would be recovered once the material is unloaded. For example, the elastic deformation is fully reversible and instantaneous. On the other hand, an irreversible mechanism means that the deformations will not be recovered once the material is unloaded. However, emphasis will be put on a reversible mechanism such as the creep/relaxation phenomenon.

Creep in materials is generally defined as a three stage phenomenon, as shown in Figure 1. Under constant loading, the creep strain rate will decrease (stage I) up to a steady state (stage II) before increasingup to the failure (stage III). The three stages are also called the primary creep, secondary creep and tertiary creep. Tertiary creep may or may not happen depending on the stress level applied. If the applied stress is low (e.g. below 0.4${\sigma }_u$for the concrete, where${\sigma }_u$is the compressive strength), tertiary creep won’t happen and a steady-state could be reached instead. Creep mechanisms are also function of the material. For metals, the mechanisms could be related to the movements of dislocations. For granular materials, mechanisms causing the creep are very different. In carbon cathode block for example, creep may originate from the cleavage and slip of basal planes within the carbon based filler particles (Picard et al., 2008).

Figure 1.

Scheme of creep strain (a) and creep strain rate (b) during a constant loading test.

In most engineering applications, however, models are developed considering the macroscopic behaviour of materials and thus the creep mechanisms (movement of dislocations, cleavage and slip of basal planes, etc.) are not directly taken into consideration. Nevertheless, the creep mechanisms could be used to give a physical meaning to parameters of the constitutive laws, as discussed in section 3.4.

Also, with the increase of numerical power calculation, complex problems are more and more solved within a three-dimensional context. Simplification of 3D complex problems to 1D ones may then be irrelevant in some cases to get approximated solutions. This may however increase the complexity of the experimental characterisation in order to identify model parameters, as discussed in the following section. This requires that the way models are established be carefully performed in order to make parameters easily identifiable in laboratory.

1.3. Three-dimensional context

As discussed in the previous section, strains and stresses are closely linked together through material characteristics and deformation mechanisms. In solid mechanic, it is useful to mathematically express this relation through constitutive laws. Also, as aforementioned, numerical modelling now requires three-dimensional constitutive models. One of the simplest expressions of a constitutive law is the Hooke’s law (4) in Voigt notation (Mase and Mase, 1999), for isotropic materials. In this case, assuming the small strain theory, only two parameters which are the Young’s modulus $E$ and the Poisson’s ratio$\nu$ have to be identified.

$$\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ 2{\epsilon }_{23}\\ 2{\epsilon }_{31}\\ 2{\epsilon }_{12}\end{array}\right]=\left[\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{23}\\ {\gamma }_{31}\\ {\gamma }_{12}\end{array}\right]=\frac{1}{E}\left[\begin{array}{cccccc}1& -\nu & -\nu & 0& 0& 0\\ -\nu & 1& -\nu & 0& 0& 0\\ -\nu & -\nu & 1& 0& 0& 0\\ 0& 0& 0& 2(1+\nu )& 0& 0\\ 0& 0& 0& 0& 2(1+\nu )& 0\\ 0& 0& 0& 0& 0& 2(1+\nu )\end{array}\right]\left[\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\sigma }_{23}\\ {\sigma }_{31}\\ {\sigma }_{12}\end{array}\right]$$E4

Constitutive laws for viscoelastic materials are however more complex than the Hooke’s law and consequently the number of parameters to identify increases. Frameworks could also be used to ensure that all parameters of the proposed laws are physically admissible. To do so, it is possible to use the thermodynamic of irreversible processes as the framework. Based on the concepts of continuum mechanics and irreversible thermodynamics, the Clausius-Duhem inequality is obtained for given problems where dissipation mechanisms are of importance, e.g., viscous deformation. Fundamental equations leading to a generic form of the Clausius-Duhem inequality have been well covered by many authors (Bazarov, 1964; Coussy, 2010; Lemaitre and Chaboche, 1990; Mase and Mase, 1999) and thus will be only summarized later in section 3. Based on the generic form of the Clausius-Duhem inequality, models or constitutive laws are further developed considering various assumptions closely related to materials of interest.

2. Materials of interest

The needs of the constitutive laws presented in this chapter originated from projects with the aluminium industry. Those laws have been used in numerical simulations (D'Amours et al., 2003; Picard et al., 2008; Richard et al., 2005; Richard et al., 2006) of the Hall-Héroult electrolysis cell (Figure 2) used for aluminum production (Sørlie and Øye, 2010). The materials of interest in these projects were mainly carbon cathode blocks, carbon anode blocks and the ramming paste. All those materials are enclosed in the potshell, thereby restraining their thermal expansion and/or chemical expansion (swelling) in some cases. In fact, the Hall-Héroult electrolysis cell is built at room temperature while its operational temperature is near 1000°C. A preheating phase prior to the aluminum production phase is thus performed to avoid or minimize thermal shocks in materials. The preheating of Hall-Héroult electrolysis cell is achieved using various methods (Sørlie and Øye, 2010). The one of interest in this work is commonly named as the electrical preheating that uses the Joule effect to heat the cell lining materials. This method is performed by passing a high density current (up to 1 A/cm² at the cathode) from the anodes to the cathode through a coke bed as shown in Figure 3.

This preheating phase is critical for the industry since it can have a major impact on the cell lifespan (Tessier et al., 2011; Tessier et al., 2010) and was thus closely investigated over the past years through numerical simulations (D'Amours et al., 2003; Marceau et al., 2011; Richard et al., 2006; Sun et al., 2004). It is therefore required to develop constitutive laws for the three aforementioned materials to feed numerical models. The targeted materials are briefly presented below.

Figure 2.

Diagram of the Hall-Héroult electrolysis cell-from (Richard, 2004).

Figure 3.

Diagram of the electrical preheating stage of the Hall-Héroult electrolysis cell- from (Goulet, 2004).

The carbon cathode consists of several carbon blocks joined together by ramming paste. These blocks are located in the liningof the electrolysis cell (Sørlie and Øye, 2010). The cathode block, as well as the ramming paste, consists of filler particles (derived from coke and graphite) mixed up with a binder (e.g. pitch). However, before the cells are started, the cathode blocks are baked while the interblock ramming paste remains unbaked (green) and will only be partially baked during the pot start-up heating (Sørlie and Øye, 2010). Depending on the block type used, heat treatment can be applied to the filler particles and/or blocks at different stages in the manufacturing process.

As aforementioned, during the preheating phase, temperature in the electrolysis cell starts from ambient and slowly increases to the desired operational value, usually around 960°C. Thermal expansion of all materials of the lining (steel shell, carbon cathode, ramming paste, refractory concrete, etc.) will thus induce mechanical stresses that could lead to material failures. Elasticity, plasticity, damaging, hardening, softening, creep/relaxation, etc. are different mechanisms inducing stresses/strains that may eventually lead to mechanical failure of the carbon materials (ramming paste and cathode). This implies that constitutive laws must take into account temperature evolution as well as many other phenomena (e.g. chemical contamination). Also, the carbon cathode block and ramming paste may exhibit damaging at very low stresses down to$0.2{\sigma }_u$. Linear creep is thus mainly observed at stress levels under$0.2{\sigma }_u$ (Picard et al., 2008) where no other phenomenon takes place.

Even though carbon anodes are not confined within a steel shell like cathode, thermal expansion of the anodic beams as well as thermal shocks may lead to mechanical failure of anodes. Creep/relaxation could also play an important role, mainly in the stub-hole region where high localized stresses could be induced by irregular geometries. Anode properties are very similar to carbon cathode and ramming paste, therefore, constitutive laws developed for one of the materials can be easily adapted to the others.

3. Constitutive laws

3.1. Mathematical notation

In this work, the following mathematical notation is used for scalars, vectors and tensors:

$$\begin{array}{cc}a,A:& \text{Scalars}\\ \left\{a\right\},⟨a⟩,\overrightarrow{a}:& \text{Vectors}\\ \left[A\right]:& \text{Matrix}\\ \left[A\right]:& \text{Second-order tensor}\\ \underset{¯}{\underset{¯}{\left[A\right]}}:& \text{Fourth-order tensor}\end{array}$$

The double contraction product will also be used in this work and for second-order tensors is defined as:

$$\left[A\right]:\left[B\right]=\text{tr}(A^{T}B)=\text{tr}(B{A}^T)=\text{tr}(A{B}^T)=\text{tr}(B^{T}A)$$E5

3.2. Thermodynamic framework

The methodology presented in this section is well covered in the literature (Bazarov, 1964; Coussy, 1995; Lemaitre and Chaboche, 1990; Mase and Mase, 1999). So, it will only be summarized here.

The thermodynamic framework implies that all systems must be based on the first and second thermodynamic principles. In continuum mechanic, the first principle which is the energy conservation can be written as:

$$\frac{dE}{dt}+\frac{dC}{dt}=P_{ext}+{\displaystyle \sum other}$$E6

$$C={\displaystyle \underset{v}{\int }\rho \overrightarrow{v}\cdot \overrightarrow{v}}dv$$E7

$$E={\displaystyle \underset{v}{\int }\rho edv}$$E8

where$e$is the specific internal energy, $E$is the internal energy, $C$is the kinematic energy, $P_{ext}$is the power of external forces and ${\displaystyle \sum other}$is related to all other kinds of energy (thermal, electrical, magnetic, etc.). Also, in equations (7) and (8), $\overrightarrow{v}$is the velocity and $v$ refers to the volume element. All materials also exhibit energy dissipation under external forces. If the dissipation level is of importance for a given problem, the second principle must then be used to obtain the Clausius-Duhem inequality. The inequality is obtained by first establishing a relation between the specific internal energy$e$, as described in Eq. (8), and the specific entropy density$s$:

$$e=\psi +Ts$$E9

while considering

$$S={\displaystyle \underset{v}{\int }\rho sdv}$$E10

where$T$is the temperature in Kelvin and$\psi$ is the Helmholtz free energy. In solid mechanic while considering a thermomechanical problem, the first principle, i.e. Eq. (6), can be rewritten as:

$$\begin{array}{cc}\frac{dE}{dt}+\frac{dC}{dt}=\frac{\partial W_{ext}}{\partial t}+\frac{\partial Q}{\partial t},& \frac{\partial W_{ext}}{\partial t}=P_{ext}\end{array}$$E11

where $\partial Q/\partial t$ is the heat rate received by the domain considered. Starting from equation (11) the Clausius-Duhem inequality has to be defined either in initial configuration or deformed configuration. In the latter case, i.e., deformed configuration, by using equations (7) and (8), equation (11) becomes:

$$\frac{d}{dt}\left(\frac{1}{2}{\displaystyle \underset{v}{\int }\rho \overrightarrow{v}\cdot }\overrightarrow{v}+\rho edv\right)={\displaystyle \underset{v}{\int }(\overrightarrow{f}\cdot \overrightarrow{v}+r)}dv+{\displaystyle \underset{s}{\int }(\overrightarrow{t}\cdot \overrightarrow{v}-\overrightarrow{q}\cdot \overrightarrow{n})}ds$$E12

with

$$P_{ext}={\displaystyle \underset{v}{\int }\overrightarrow{f}\cdot \overrightarrow{v}dv+{\displaystyle \underset{s}{\int }\overrightarrow{t}\cdot \overrightarrow{v}ds}}$$E13

and

$$\frac{\partial Q}{\partial t}={\displaystyle \underset{v}{\int }rdv-{\displaystyle \underset{s}{\int }\overrightarrow{q}\cdot \overrightarrow{n}ds}}$$E14

where $\overrightarrow{f}$ is a body load, $\overrightarrow{t}$is a prescribed stress, $r$is a heat source, $\overrightarrow{q}$is a prescribed heat flux and $\overrightarrow{n}$ is a normal vector to a surface, considering that:

$$\overrightarrow{t}=\left[\sigma \right]\overrightarrow{n}$$E15

and

$$\underset{\text{mass conservation}}{\underbrace{\frac{d\rho }{dt}+\rho div(\overrightarrow{v})=0}}$$E16

where $\left[\sigma \right]$is a second order stress tensor, Eq. (12) can be rewritten, after few other considerations, as:

$$r=\rho \frac{de}{dt}-\left[\sigma \right]:\left[d\right]+div(\overrightarrow{q})$$E17

here

$$\left[d\right]=sym(grad(\overrightarrow{v}))=\left[\dot{\epsilon }\right]$$E18

In order to get the Clausius-Duhem equation, the second principle must then be introduced and take the following form:

$$\frac{dS}{dt}-{\dot{S}}_{ext}\ge 0$$E19

with

$${\dot{S}}_{ext}={\displaystyle \underset{r}{\int }\frac{r}{T}dv-{\displaystyle \underset{s}{\int }\frac{\overrightarrow{q}\cdot \overrightarrow{n}}{T}}}ds$$E20

where$S$was already defined in Eq. (10). After some mathematical manipulations, equation (20) becomes:

$$\rho \frac{ds}{dt}-\frac{r}{T}+\frac{1}{T}div(\overrightarrow{q})-\frac{1}{T^2}\overrightarrow{q}\cdot grad(T)\ge 0$$E21

Finally, combining the two equations obtained with the first and second principles, i.e., equations (17) and (21), and the Helmholtz free energy (9), a general form of the Clausius-Duhem inequality is obtained:

$$\left[\sigma \right]:\left[\dot{\epsilon }\right]-\rho \left(\frac{d\psi }{dt}+s\frac{dT}{dt}\right)-\frac{\overrightarrow{q}}{T}\cdot grad(T)\ge 0$$E22

One of the most important parts in modelling is then to establish a link using the Clausius-Duhem inequality to materials parameters, e.g., by defining the form of Helmholtz free energy$\psi$. The parameters will be identified later in laboratory. Different approaches can be used to establish the form of $\psi$ and the choice is of course function of the problem. In some cases, it is also possible to work with the Gibbs free energy instead of the Helmholtz, e.g. (Haslach, 2009; Lion et al., 2010; Schapery, 1997). For hyperelastic materials, the Ogden model, e.g. (Holzapfel, 1996; Li and Lua, 2009; Ogden, 1972), can also be used. Two different approaches will be presented in this work. The first one will link materials properties to the Clausius-Duhem inequality through generic internal state variables (Fafard et al., 2001) and the second one will use phenomenological approach (Picard et al., 2008) based on rheology of materials. The two approaches have already been published and will thus be only summarized here in sections 3.3 and 3.4.

3.3. Generic internal state variables

The Clausis-Duhem inequality (22) has established a relation between the stress and the strain through the Helmholtz free energy$\psi$. Considering a viscoelastic material and constant and uniform temperature, it can be assumed that the Helmholtz free energy will only be function of the total strain and an undetermined number of second order tensors internal state variables $\left[q_{\alpha }\right]$ ($\alpha =\mathrm{1,...,}n$). So, Eq. (22) becomes:

$$\left[\sigma \right]:\left[\dot{\epsilon }\right]-\dot{\psi }(\left[\epsilon \right],\left[q_{\alpha }\right])\ge 0$$E23

The time derivation of Helmholtz free energy and the substitution into (23) gives

$$\left(\left[\sigma \right]-\frac{\partial \psi }{\partial \epsilon }\right):\left[\dot{\epsilon }\right]-\frac{\partial \psi }{\partial \left[q_{\alpha }\right]}:\frac{\partial \left[q_{\alpha }\right]}{\partial t}\ge 0$$E24

The inequality (24) must be satisfied for any thermodynamic process. According to the local state law (Lemaitre and Chaboche, 1990), the Clausius-Duhem inequality leads to the following state equation:

$$\left[\sigma \right]=\frac{\partial \psi }{\partial \left[\epsilon \right]}$$E25

The Clausius-Duhem inequality thus becomes:

$$-\frac{\partial \psi }{\partial \left[q_{\alpha }\right]}:\frac{\partial \left[q_{\alpha }\right]}{\partial t}\ge 0$$E26

Assuming the linearity of the dissipative mechanism, the so-called dissipative potential ϕ can be chosen as a quadratic form of its arguments (Coussy, 1995; Lemaitre and Chaboche, 1990; Valanis, 1972). Hence, the dissipative potential is:

$$\varphi =\frac{1}{2}{\displaystyle \sum _{}\left[{\dot{q}}_{\alpha }\right]}:\underset{¯}{\underset{¯}{\left[b_{\alpha }\right]}}:\left[{\dot{q}}_{\alpha }\right]$$E27

Where the coefficients of $\underset{¯}{\underset{¯}{\left[b_{\alpha }\right]}}$must be identified. In the present case, the coefficients of $\underset{¯}{\underset{¯}{\left[b_{\alpha }\right]}}$ will be considered constant, i.e., independent of time$t$, temperature$T$, etc. The non-constant coefficients case is well covered by (Fafard et al., 2001). To satisfy the inequality (26), the complementary evolution law can then be written as

$$\frac{\partial \varphi }{\partial \left[{\dot{q}}_{\alpha }\right]}+\frac{\partial \psi }{\partial \left[q_{\alpha }\right]}=0$$E28

Furthermore, the following quadratic form of the Helmholtz free energy is assumed for isotropic, orthotropic or anisotropic materials (Valanis, 1972) for reversible dissipative mechanical problem:

$$\psi =\frac{1}{2}\left[\epsilon \right]:\underset{¯}{\underset{¯}{\left[E\right]}}:\left[\epsilon \right]+{\displaystyle \sum _{\alpha }\left[q_{\alpha }\right]:\underset{¯}{\underset{¯}{\left[B_{\alpha }\right]}}:\left[\epsilon \right]}+\frac{1}{2}{\displaystyle \sum _{\alpha }\left[q_{\alpha }\right]:\underset{¯}{\underset{¯}{\left[A_{\alpha }\right]}}:\left[q_{\alpha }\right]}$$E29

where $$\underset{¯}{\underset{¯}{\left[E\right]}},$$$\underset{¯}{\underset{¯}{\left[B_{\alpha }\right]}}$ and $$\underset{¯}{\underset{¯}{\left[A_{\alpha }\right]}}$$ are symmetric tensors. The first part of (29) is the energy due to the instantaneous deformation. The two other terms represent the change in free energy due to the viscoelastic behaviour of the material. The stress tensor can be obtained from (25) and (29) that takes the following form:

$$\left[\sigma \right]=\frac{\partial \psi }{\partial \left[\epsilon \right]}=\underset{¯}{\underset{¯}{\left[E\right]}}:\left[\epsilon \right]+{{\displaystyle \sum _{\alpha }\underset{¯}{\underset{¯}{\left[B_{\alpha }\right]}}}}^T:\left[q_{\alpha }\right]$$E30

The equations (27) to (30) lead to a set of equations for each value of$\alpha$:

$$\underset{¯}{\underset{¯}{\left[b_{\alpha }\right]}}:\left[{\dot{q}}_{\alpha }\right]+\underset{¯}{\underset{¯}{\left[A_{\alpha }\right]}}:\left[q_{\alpha }\right]=-\underset{¯}{\underset{¯}{\left[B_{\alpha }\right]}}:\left[\epsilon \right]$$E31

Assuming isotropic materials, the following assumption is proposed: the topology of all tensors is similar to that of the elastic one. This means thatfor each tensor there are two unknown coefficients analogous to the Young’s modulus and Poisson’s coefficient (Fafard et al., 2001; Picard et al., 2008). Tensors of Eq. (30) and (31) can then be expressed as:

$$\begin{array}{cc}\underset{¯}{\underset{¯}{\left[A_{\alpha }\right]}}=\left[A_{\alpha }(E_{A_{\alpha }},{\mu }_{\alpha })\right]& \underset{¯}{\underset{¯}{\left[B_{\alpha }\right]}}=\left[B_{\alpha }(E_{B_{\alpha }},{\zeta }_{\alpha })\right]\\ \underset{¯}{\underset{¯}{\left[b_{\alpha }\right]}}=\left[b_{\alpha }(E_{b_{\alpha }},{\xi }_{\alpha })\right]& \underset{¯}{\underset{¯}{\left[E_{\alpha }\right]}}=\left[E_{\alpha }(E_E,\nu )\right]\end{array}$$E32

The analytical solution of Eq. (31) can be obtained using the diagonalization technique or modal projection (Kreyszig, 2006), Thisis fully detailed in (Fafard et al., 2001). The specific solution for the uniaxial creep case (${\sigma }_{11}=\sigma =\text{constant}$) of a cylindrical sample is

$${\epsilon }_{axial}={\epsilon }_{11}=\frac{1}{\gamma \beta }\left[\frac{1}{E}-\frac{1}{3}\frac{1}{G}\left(\frac{K^d}{K}+\gamma \frac{G^d}{G}{e}^{-{\lambda }_{b}t}\right)-\frac{1}{9}\frac{1}{K}\left(\frac{G^d}{G}+\beta \frac{K^d}{K}{e}^{-{\lambda }_{a}t}\right)\right]\sigma$$E33

$${\epsilon }_{radial}={\epsilon }_{22}={\epsilon }_{33}=\frac{1}{\gamma \beta }\left[\frac{-\nu }{E}-\frac{1}{6}\frac{1}{G}\left(\frac{K^d}{K}+\gamma \frac{G^d}{G}{e}^{-{\lambda }_{b}t}\right)-\frac{1}{9}\frac{1}{K}\left(\frac{G^d}{G}+\beta \frac{K^d}{K}{e}^{-{\lambda }_{a}t}\right)\right]\sigma$$E34

where

$$\begin{array}{cc}{\lambda }_a={\lambda }_1\gamma ;& {\lambda }_b={\lambda }_2\beta \\ \gamma =1-\frac{K^d}{K};& \beta =1-\frac{G^d}{G}\end{array}$$E35

$$\begin{array}{cc}G=\frac{E_E}{2\left(1+\nu \right)};& K=\frac{E_E}{3\left(1-2\nu \right)}\end{array}$$E36

and

$$\begin{array}{cc}{\lambda }_{1_{\alpha }}=\frac{E_{A_{\alpha }}}{E_{b_{\alpha }}}\frac{1-2{\xi }_{\alpha }}{1-2{\mu }_{\alpha }};& {\lambda }_{2_{\alpha }}=\frac{E_{A_{\alpha }}}{E_{b_{\alpha }}}\frac{1+{\xi }_{\alpha }}{1+{\mu }_{\alpha }}\end{array}$$E37

Coefficients G and K represent the instantaneous shear and compressibility modulus, respectively. In addition to elastic coefficients, one has to identify$G^d$,$K^d$ ,${\lambda }_{1_{\alpha }}$ parameters and${\lambda }_{2_{\alpha }}$. Considering the hypothesis that ${\xi }_{\alpha }={\mu }_{\alpha }$ for 3D isotropic materials leads to $${\lambda }_{1_{\alpha }}={\lambda }_{2_{\alpha }}={\lambda }_{\alpha }$$and thus reduces the number of parameters to be identified in laboratory. Finally, a link with a viscoelastic rheological model could thereafter be done by defining the rheological parameters as a function of the internal states variables. The alternative solution would be to define rheological parameters as the internal state variables at the very beginning. This corresponds to the second approach (section 3.4) presented in this chapter.

3.4. Strain Based Internal State Variables

Instead of defining the viscoelastic behaviour with generic internal state variable$\underset{¯}{\underset{¯}{\left[q_{\alpha }\right]}}$, the second approach is based on the rheological behaviour of materials. Combined with the thermodynamic of irreversible processes (TIP), rheology gives the possibility to express the mechanical behaviour with simple parameters having a strong physical meaning (Kuiken, 1994). Models can be rather simple or very complex but they are all built using a combination of spring and dashpot elements. The more common elements, also called bodies, are the Kelvin-Voigt and the Maxwell elements. Creep behaviour of the materials studied in this chapter can be expressed using a simple rheological model called the Burger body. It is possible to extend this to another model using a similar approach.

As for the generic internal state variables approach discussed in section 3.3, the rheological approach also needs internal state variables. Instead of choosing generic ISV (Fafard et al., 2001), the choice of these variables is based on the use of a phenomenological approach by assuming a rheological model representing the material behaviour. To represent the creep behaviour of the studied materials, a Kelvin-Voigt rheological model has been chosen (Picard et al., 2008) as shown in Figure 4. As in section 3.3, the small strain theory is assumed.

The parameters of the Kelvin-Voigt model and the internal are fourth order tensors, while the strain and the stress are second order tensors. An undetermined number of Kelvin-Voigt elements give flexibility to the model without increasing the complexity of the constitutive law as it will be discussed later in this section. Assuming a virgin material, having no permanent strain due to earlier where $\left[{\epsilon }^e\right]$ is the instantaneous elastic strain, where

Figure 4.

Diagram of the three-dimensional visco-elastic rheological model.

damage, constant temperature and without other environmental phenomena that can cause additive strains, the total strain $\left[\epsilon \right]$ is reduced to:

$$\left[\epsilon \right]=\left[{\epsilon }^e\right]+{\displaystyle \sum _{\alpha =1}^N\left[{\epsilon }_{\alpha }^a\right]}$$E38

$\left[{\epsilon }^e\right]$is the instantaneous elastic strain, $\left[{\epsilon }_{\alpha }^a\right]$the anelastic strain associated with the Kelvin-Voigt element$\alpha$, and $N$ is the number of Kelvin-Voigt elements.

Now that the ISV is known, it is thus possible to define the Helmholtz free energy form as:

$$\rho \psi =\frac{1}{2}\left[{\epsilon }^e\right]:\underset{¯}{\underset{¯}{\left[H^e\right]}}:\left[{\epsilon }^e\right]+\frac{1}{2}\left[{\epsilon }_{\alpha }^a\right]:\underset{¯}{\underset{¯}{\left[H_{\alpha }^a\right]}}:\left[{\epsilon }_{\alpha }^a\right]$$E39

where $\underset{¯}{\underset{¯}{\left[H^e\right]}}$ is the Voigt form of the fourth order elasticity tensor and $\underset{¯}{\underset{¯}{\left[H_{\alpha }^a\right]}}$ is the fourth order tensor related to the spring of the $\alpha$ Kelvin-Voigt element. The parameters defining those tensors must be identified using appropriate tests. In order to take into account the dissipative mechanism of the model, a dissipative potential,$\varphi$ , per unit of volume, can be chosen as a quadratic form of its arguments:

$$\rho \varphi =\frac{1}{2}{\displaystyle \sum _{\alpha =1}^N\left[{\dot{\epsilon }}_{\alpha }^a\right]:\underset{¯}{\underset{¯}{\left[{\eta }_{\alpha }^a\right]}}:\left[{\dot{\epsilon }}_{\alpha }^a\right]}$$E40

where the fourth-order tensor $\underset{¯}{\underset{¯}{\left[{\eta }_{\alpha }^a\right]}}$ is related to the dashpot of the $\alpha$ Kelvin-Voigt element. The parameters of this tensor must be identified, as well. Based on the thermodynamics of irreversible processes (TIP) and on the choice of the internal state variables (ISV), the Clausius-Duhem inequality can be written as:

$${\phi }_T+{\phi }_M\ge 0$$E41

where $${\phi }_T$$ represents the thermal terms and $${\phi }_M$$ the mechanical terms. For the visco-elastic problem, and assuming that the thermal part, $${\phi }_T$$, is satisfied, Eq. (41) can be rewritten using Eq. (39), with the mechanical terms onlyas:

$${\phi }_M=\left[\sigma \right]:\left[\dot{\epsilon }\right]-\dot{\psi }\ge 0$$E42

$$\left(\left[\sigma \right]-\rho \frac{\partial \psi }{\partial \left[{\epsilon }^e\right]}\right):\left[{\dot{\epsilon }}^e\right]+{\displaystyle \sum _{\alpha =1}^N\left(\left[\sigma \right]-\rho \frac{\partial \psi }{\partial \left[{\epsilon }_{\alpha }^a\right]}\right)}:\left[{\dot{\epsilon }}_{\alpha }^a\right]\ge 0$$E43

where $\left[\sigma \right]$ is the stress tensor. The first term in the Eq. (43) is associated with the elastic strain (Lemaitre and Chaboche, 1990). Since this latter strain is related to a reversible mechanism, it can be assumed that:

$$\left[\sigma \right]=\rho \frac{\partial \psi }{\partial \left[{\epsilon }^e\right]}$$E44

The Clausius-Duhem inequality (43) can then be rewritten:

$${\displaystyle \sum _{\alpha =1}^N\left(\left[\sigma \right]-\rho \frac{\partial \psi }{\partial \left[{\epsilon }_{\alpha }^a\right]}\right)}:\left[{\dot{\epsilon }}_{\alpha }^a\right]\ge 0$$E45

To satisfy the former equation, the complementary evolution law proposed in the following equation is postulated:

$$\left[\sigma \right]-\rho \frac{\partial \psi }{\partial \left[{\epsilon }_{\alpha }^a\right]}=\frac{\partial \varphi }{\partial \left[{\dot{\epsilon }}_{\alpha }^a\right]}$$E46

Assuming that $\underset{¯}{\underset{¯}{\left[{\eta }_{\alpha }^a\right]}}$ are positive, the use of Eq. (40) and Eq. (46) ensures that Eq. (45) is satisfied and thus the Clausius-Duhem equation is also satisfied. Using the definitions of $\psi$ and $\varphi$ in Eq. (46) leads to Voigt notation (Mase and Mase, 1999),

$$\left[{\eta }_{\alpha }^a\right]\left\{{\dot{\epsilon }}_{\alpha }^a\right\}+\left[H_{\alpha }^a\right]\left\{{\epsilon }_{\alpha }^a\right\}=\left\{\sigma \right\}$$E47

for a given$\alpha$. Assuming isotropic materials, tensors of Eq. (39) and (40) could now be expressed as:

$$\begin{array}{l}\underset{¯}{\underset{¯}{\left[H^e\right]}}=\left[H^e(E,\nu )\right]\\ \underset{¯}{\underset{¯}{\left[H_{\alpha }^a\right]}}=\left[H_{\alpha }^a(E_{H_{\alpha }^a},{\mu }_{\alpha })\right]\\ \underset{¯}{\underset{¯}{\left[{\eta }_{\alpha }^a\right]}}=\left[{\eta }_{\alpha }^a(E_{{\eta }_{\alpha }^a},{\varsigma }_{\alpha })\right]\end{array}$$E48

In order to identify the parameters of tensors of Eq. (48), Eq. (47) must be solved. The analytical solution can be obtained using the diagonalization technique (Kreyszig, 2006) and the same assumptions aforementioned in section 3.3 will still be valid. The creep case will also be considered here. The solution of Eq. (47) is thus in Voigt notation:

$$\left\{\epsilon \right\}={\left[H^e\right]}^{-1}\left\{\sigma \right\}+{\displaystyle \sum _{\alpha =1}^N\left[X_{\alpha }\right]\left\{\begin{array}{c}{F}_{{\alpha }_1}{\displaystyle \underset{0}{\overset{\tau }{\int }}{e}^{-{\lambda }_{{\alpha }_1}(t-\tau )}d\tau }\\ F_{{\alpha }_2}{\displaystyle \underset{0}{\overset{\tau }{\int }}{e}^{-{\lambda }_{{\alpha }_2}(t-\tau )}d\tau }\\ ⋮\\ F_{{\alpha }_6}{\displaystyle \underset{0}{\overset{\tau }{\int }}{e}^{-{\lambda }_{{\alpha }_6}(t-\tau )}d\tau }\end{array}\right\}}$$E49

with

$$\left\{F_{\alpha }\right\}={\left[X_{\alpha }\right]}^{-1}{\left[{\eta }_{\alpha }^a\right]}^{-1}\left\{\sigma \right\}$$E50

where $\left[X_{\alpha }\right]$is a matrix containing the eigenvectors $\left\{x_{\alpha }\right\}$ associated with the eigenvalues ${\lambda }_{\alpha }$ of the following system:

$$\left(\left[{\eta }_{\alpha }^a\right]\left[H_{\alpha }^a\right]-{\lambda }_{\alpha }\left[I\right]\right)\left\{x_{\alpha }\right\}=0$$E51

The complete solution is available in (Picard et al., 2008) for further review. Since all tensors have the same topology, it can be shown that$\left[X_{\alpha }\right]=\left[X\right]$ and that${\lambda }_2={\lambda }_3={\lambda }_4={\lambda }_5={\lambda }_6$. Also, the eigenvalues of Eq. (49) are

$$\begin{array}{cc}{\lambda }_{{\alpha }_1}=\frac{E_{H_{\alpha }^a}(1-2{\varsigma }_{\alpha })}{E_{{\eta }_{\alpha }^a}(1-2{\mu }_{\alpha })};& {\lambda }_{{\alpha }_{\mathrm{2,3,4,5,6}}}=\frac{E_{H_{\alpha }^a}(1+{\varsigma }_{\alpha })}{E_{{\eta }_{\alpha }^a}(1+{\mu }_{\alpha })}\end{array}$$E52

Finally, for the uniaxial creep case (constant loading) of a cylindrical sample, Eq. (49) can be rewritten as

$${\epsilon }_{axial}={\epsilon }_{11}=\underset{\text{Elastic strain}}{\underbrace{-\frac{{\sigma }_1}{E_{H^e}}}}-\underset{\text{Creep strain}}{\underbrace{{\sigma }_1\left(\frac{1}{9}{\displaystyle \sum _{\alpha =1}^N{K}_{\alpha }^{a^{-1}}\frac{1-e^{-{\lambda }_{{\alpha }_1}t}}{{\lambda }_{{\alpha }_1}}}+\frac{1}{3}{\displaystyle \sum _{\alpha =1}^N{G}_{\alpha }^{a^{-1}}\frac{1-e^{-{\lambda }_{{\alpha }_2}t}}{{\lambda }_{{\alpha }_2}}}\right)}}$$E53

$${\epsilon }_{radial}={\epsilon }_{22}={\epsilon }_{33}=\underset{\text{Elastic strain}}{\underbrace{{\sigma }_1\frac{\nu }{E_{H^e}}}}-\underset{\text{Creep strain}}{\underbrace{{\sigma }_1\left(\frac{1}{9}{\displaystyle \sum _{\alpha =1}^N{K}_{\alpha }^{a^{-1}}\frac{1-e^{-{\lambda }_{{\alpha }_1}t}}{{\lambda }_{{\alpha }_1}}}-\frac{1}{6}{\displaystyle \sum _{\alpha =1}^N{G}_{\alpha }^{a^{-1}}\frac{1-e^{-{\lambda }_{{\alpha }_2}t}}{{\lambda }_{{\alpha }_2}}}\right)}}$$E54

with

$$\begin{array}{cc}{K}_{\alpha }^{a^{-1}}=\frac{3(1-2{\varsigma }_{\alpha })}{E_{{\eta }_{\alpha }^a}};& G_{\alpha }^{a^{-1}}=\frac{2(1+{\varsigma }_{\alpha })}{E_{{\eta }_{\alpha }^a}}\end{array}$$E55

where ${\sigma }_1$ is the constant axial stress, ${\epsilon }_1$is the axial strain and ${\epsilon }_r$the hoop strain. All other strain components are null. Thecoefficients $$K_{\alpha }^{a^{-1}}$$and $$G_{\alpha }^{a^{-1}}$$are related to the hydrostatic and deviatoric creep mechanisms, respectively. Moreover, all $${\lambda }_{{\alpha }_i}$$ coefficients are positive and $$1/{\lambda }_{{\alpha }_i}$$are related to a relaxation time. The elastic parameters $E$ and $\nu$(elastic properties) must be identified through standard compression tests, while the anelastic parameters$E_{H_{\alpha }^a}$,${\mu }_{\alpha }$ ,$E_{{\eta }_{\alpha }^a}$ and${\zeta }_{\alpha }$must be identified through experimental creep tests.

4. Academic case study

4.1. Finite elements model

Independently of the approach used to get a constitutive law, parameters of the latter must be identified in laboratory. To simplify the identification process, creep tests are preferred (Fafard et al., 2001; Picard et al., 2008) over any other tests to identify the constitutive laws parameters discussed in this work. In such case where the number of parameters to identify is small, methods such as Newton-Raphson algorithm could be used to minimize the error of an objective function through use of the least squares method. If the number of parameters is large, an alternative method like genetic algorithm may be more suitable.

In few cases, parameters identification at a reference state (generally at room temperature) will be sufficient. This implies that in most situations constitutive law parameters must evolve to be representative of the real situation (e.g., the carbon materials have to operate at a temperature near 1000°C while the concrete may be influenced by the relative humidity of the atmosphere).

The analytical solutions of the proposed models in section 3 were obtained by assuming a reference state, i.e., a virgin material at 25°C. As aforementioned, parameters evolution was not taken into account until now. Since these models are dedicated to being used in finite element code (Picard et al., 2008; Richard et al., 2005), each parameter could have its own scalar evolution function. To get this function, the methodology used consists of building datasets of creep measures under different states, e.g., at different temperatures. Then, parameters of the three-dimensional model are identified by using the appropriated parameters identification process, at each different state or temperature, for example. The function obtained will then be used in finite element code to take into account changing environment (temperature, relative humidity, etc.). To do this, the analytical solution of the two proposed approaches must be discretized. Regarding this, in the case of the strain based one (section 3.4), the Eq. (47) can be rewritten as:

$$\left[{\eta }_{\alpha }^a\right]\left\{{\dot{\epsilon }}_{\alpha }^a\right\}+(\left[H_{\alpha }^a\right]+\left[H^e\right])\left\{{\epsilon }_{\alpha }^a\right\}=\left[H^e\right]\left(\left\{\epsilon \right\}-{\displaystyle \sum _{j=\mathrm{1,}j\ne \alpha }^N{\epsilon }_j^a}\right)$$E56

Assuming a forward difference Euler scheme (Reddy, 1993)

$${\left[{\chi }_{\alpha }^a\right]}_{t+\Delta t}{\left\{{\epsilon }_{\alpha }^a\right\}}_{t+\Delta t}=\Delta t{\left[H^e\right]}_{t+\Delta t}{\left\{\epsilon \right\}}_{t+\Delta t}+{\left[{\eta }_{\alpha }^a\right]}_{t+\Delta t}{\left\{{\epsilon }_{\alpha }^a\right\}}_t-{\displaystyle \sum _{j=\mathrm{1,}j\ne \alpha }^N\Delta t{\left[H^e\right]}_{t+\Delta t}{\left\{{\epsilon }_j^a\right\}}_{t+\Delta t}}$$E57

where

$${\left[{\chi }_{\alpha }^a\right]}_{t+\Delta t}={\left[{\eta }_{\alpha }^a\right]}_{t+\Delta t}+\Delta t\left({\left[H_{\alpha }^a\right]}_{t+\Delta t}+{\left[H^e\right]}_{t+\Delta t}\right)$$E58

Moreover, it can be assumed that

$${\left\{{\epsilon }_{\alpha }^a\right\}}_{t+\Delta t}={\left\{{\epsilon }_{\alpha }^a\right\}}_t+\Delta \left\{{\epsilon }_{\alpha }^a\right\}$$E59

Equation (57) and (59) lead to a system of equation of dimension $N\times N$described as follows:

$${\left[\begin{array}{cccc}\left[{\chi }_1^a\right]& \Delta t\left[H^e\right]& \cdots & \Delta t\left[H^e\right]\\ \Delta t\left[H^e\right]& \left[{\chi }_2^a\right]& \ddots & ⋮\\ ⋮& \ddots & \ddots & \Delta t\left[H^e\right]\\ \Delta t\left[H^e\right]& \cdots & \Delta t\left[H^e\right]& \left[{\chi }_N^a\right]\end{array}\right]}_{t+\Delta t}\left\{\begin{array}{c}{\left\{{\epsilon }_1^a\right\}}_{t+\Delta t}\\ {\left\{{\epsilon }_2^a\right\}}_{t+\Delta t}\\ ⋮\\ {\left\{{\epsilon }_N^a\right\}}_{t+\Delta t}\end{array}\right\}=\left\{\begin{array}{c}\Delta t{\left[H^e\right]}_{t+\Delta t}{\left\{\epsilon \right\}}_{t+\Delta t}+{\displaystyle \sum _{\alpha =\mathrm{1,}\alpha \ne 1}^N{\left[{\eta }_{\alpha }^a\right]}_{t+\Delta t}{\left\{{\epsilon }_{\alpha }^a\right\}}_t}\\ \Delta t{\left[H^e\right]}_{t+\Delta t}{\left\{\epsilon \right\}}_{t+\Delta t}+{\displaystyle \sum _{\alpha =\mathrm{1,}\alpha \ne 2}^N{\left[{\eta }_{\alpha }^a\right]}_{t+\Delta t}{\left\{{\epsilon }_{\alpha }^a\right\}}_t}\\ \\ \Delta t{\left[H^e\right]}_{t+\Delta t}{\left\{\epsilon \right\}}_{t+\Delta t}+{\displaystyle \sum _{\alpha =\mathrm{1,}\alpha \ne N}^N{\left[{\eta }_{\alpha }^a\right]}_{t+\Delta t}{\left\{{\epsilon }_{\alpha }^a\right\}}_t}\end{array}\right\}$$E60

From Eq. (59) and (60), the strain increment that will be used in finite element software is:

$$\left\{\begin{array}{c}\Delta \left\{{\epsilon }_1^a\right\}\\ \Delta \left\{{\epsilon }_1^a\right\}\\ ⋮\\ \Delta \left\{{\epsilon }_N^a\right\}\end{array}\right\}={\left[\begin{array}{cccc}\left[{\chi }_1^a\right]& \Delta t\left[H^e\right]& \cdots & \Delta t\left[H^e\right]\\ \Delta t\left[H^e\right]& \left[{\chi }_2^a\right]& \ddots & ⋮\\ ⋮& \ddots & \ddots & \Delta t\left[H^e\right]\\ \Delta t\left[H^e\right]& \cdots & \Delta t\left[H^e\right]& \left[{\chi }_N^a\right]\end{array}\right]}_{t+\Delta t}^{-1}\left\{\begin{array}{c}\Delta t{\left[H^e\right]}_{t+\Delta t}\left({\left\{\epsilon \right\}}_{t+\Delta t}-{\displaystyle \sum _{j=1}^N{\left\{{\epsilon }_j^a\right\}}_t}\right)-\Delta t{\left[H_1^a\right]}_{t+\Delta t}{\left\{{\epsilon }_1^a\right\}}_t\\ \Delta t{\left[H^e\right]}_{t+\Delta t}\left({\left\{\epsilon \right\}}_{t+\Delta t}-{\displaystyle \sum _{j=1}^N{\left\{{\epsilon }_j^a\right\}}_t}\right)-\Delta t{\left[H_2^a\right]}_{t+\Delta t}{\left\{{\epsilon }_2^a\right\}}_t\\ ⋮\\ \Delta t{\left[H^e\right]}_{t+\Delta t}\left({\left\{\epsilon \right\}}_{t+\Delta t}-{\displaystyle \sum _{j=1}^N{\left\{{\epsilon }_j^a\right\}}_t}\right)-\Delta t{\left[H_N^a\right]}_{t+\Delta t}{\left\{{\epsilon }_N^a\right\}}_t\end{array}\right\}$$E61

Other quantities are also needed by finite element software, such as the tangent matrix for the use of Newton-Raphson methods, and are also function of the number $N$ of Kelvin-Voigt elements. In the case of limiting the rheological model to$N=2$, the tangent matrix is:

$$\left[C\right]={\left[H^e\right]}_{t+\Delta t}-\Delta t{\left[H^e\right]}_{t+\Delta t}\left({\left[{\chi }_{12}^a\right]}_{t+\Delta t}^{-1}{\left[{\chi }_{22}^a\right]}_{t+\Delta t}+{\left[{\chi }_{21}^a\right]}_{t+\Delta t}^{-1}{\left[{\chi }_{11}^a\right]}_{t+\Delta t}\right){\left[H^e\right]}_{t+\Delta t}$$E62

where

$$\begin{array}{c}{\left[{\chi }_{12}^a\right]}_{t+\Delta t}={\left[{\chi }_1^a\right]}_{t+\Delta t}-\Delta t^2{\left[H^e\right]}_{t+\Delta t}{\left[{\chi }_2^a\right]}_{t+\Delta t}^{-1}{\left[H^e\right]}_{t+\Delta t}\\ {\left[{\chi }_{21}^a\right]}_{t+\Delta t}={\left[{\chi }_2^a\right]}_{t+\Delta t}-\Delta t^2{\left[H^e\right]}_{t+\Delta t}{\left[{\chi }_1^a\right]}_{t+\Delta t}^{-1}{\left[H^e\right]}_{t+\Delta t}\\ {\left[{\chi }_{11}^a\right]}_{t+\Delta t}=\left[I\right]-\Delta t^2{\left[H^e\right]}_{t+\Delta t}{\left[{\chi }_1^a\right]}_{t+\Delta t}^{-1}\\ {\left[{\chi }_{22}^a\right]}_{t+\Delta t}=\left[I\right]-\Delta t^2{\left[H^e\right]}_{t+\Delta t}{\left[{\chi }_2^a\right]}_{t+\Delta t}^{-1}\end{array}$$E63

The same methodology can be applied to the generic internal state variables approach to get the corresponding strain increment and tangent matrix.

4.2. Case study

The creep/relaxation behaviour of the ramming paste material hasn’t been rigorously investigated yet. However, preliminary numerical results from (D'Amours, 2004) and (Richard, 2004) have clearly shown the importance of this phenomenon. Thus, it is relevant to establish the creep behaviour of this material, even in the form of a qualitative model, in order to take it into account in the pre-heating simulation of the thermo-electro-mechanical behaviour of a Hall-Héroult electrolysis cell. Also, as mentioned in section 2, the ramming paste is initially green and it begins to bake during the pre-heating phase. Thus, the mechanical properties of the paste evolve with the so-called “baking index” $\xi$proposed by (D'Amours, 2004) which is temperature dependant and irreversible. The viscoelastic rheological model presented in section 3.4 has thus been used to numerically simulate the thermo-mechanical behaviour of ramming paste materials (Picard, 2007) as an academic case study. Based on experimental and extrapolated data (Picard, 2007), the following parameters evolution of the viscoelastic model with two Kelvin-Voigt elements ($N=2$) is obtained:

$$\begin{array}{cc}{E}_{{\eta }_1^a}(\text{MPa}\cdot \text{h})& 271263e^{3.1104\xi }\\ {\zeta }_1& 0.1444{\xi }^2-0.3367\xi +0.3922\\ E_{H_1^a}(\text{MPa})& 467276e^{0.2629\xi }\\ {\mu }_1& -0.0444{\xi }^2+.02267\xi +0.2078\\ E_{{\eta }_2^a}(\text{MPa}\cdot \text{h})& 440e^{6.9123\xi }\\ {\zeta }_2& -0.5056{\xi }^2+0.4783\xi +0.0972\\ E_{H_2^a}(\text{MPa})& 45e^{8.7068\xi }\\ {\mu }_2& -0.2889{\xi }^2+0.1733\xi +0.2456\end{array}$$E64

Moreover, it is relevant to note that the creep behaviour of the ramming paste in this case study is function of the baking index only, even at various temperatures. In fact, up to now, no creep/relaxation results are available for the ramming paste at various elevated temperatures. Also, it should be mentioned that the baking and the temperature effects are considered independent. The baking effect is related to the microstructure of the paste which is a function of the highest temperature reached by the paste (irreversible process), while the temperature effect is related to the actual temperature of the paste (reversible process).

The proposed model used by (D'Amours, 2004) was adopted to test evolution of the parameters of the ramming paste. This simple model, illustrated in Figure 5, consists of a quarter carbon block (C) and half a ramming paste seam (P) simulating the electrical preheating phase of an Hall-Héroult electrolysis cell (Sørlie and Øye, 2010). A similar model was also used by (Richard et al., 2005). All the details of that method are available in (Picard, 2007). The four elements on each extremity of the cathode (C) do not contribute to the heat transfer related to the anode since they are not covered by the bed of coke. The elements covered by the anode are represented by the “Anode shadow” region on the model. A forced convection ($h=500$W/m²/K) was imposed in this “Anode shadow” region, where the temperature $T_{\infty }$was ramped linearly from 20°C to 1000°C in a 40 hour interval. All surfaces of the anode except the upper one are insulated. For more information refer to (Picard, 2007).

Figure 5.

Simple finite element model for the carbon cell lining materials – Taken from (D'Amours, 2004).

A first simulation was performed without any creep of the ramming paste to get a reference for future comparison. The strain results after a preheating of 40 hours are plotted in Figure 6.

The important point to note here is that a relatively high stress level in the longitudinal direction (X) is obtained in the materials that consequently initiates a high plastic strain amplitude in the ramming paste. This is in contradiction with the experimental observations (D'Amours, 2004; Richard, 2004) where an important creep/relaxation phenomenon should take place. In fact, the role of the ramming paste in the Hall-Héroult electrolysis cell is to let the cathode bottom block expand during the pre-heating phase without leading to any mechanical failure (e.g., cracking where liquid could leak in). It is then relevant to add a viscoelastic behaviour such as the creep/relaxation to the model. The simulation was then run again with the addition of creep/relaxation in the ramming paste and the results are presented in Figure 7.

These last numerical results show that the plastic strain is obviously greatly influenced by the presence of the creep/relaxation phenomenon. In fact, the level of the plastic strain was considerably reduced from 0.004 to 0.0018, i.e a reduction ratio of 2.2 based on the reference case (Figure 6), which does not take into account creep behaviour of the ramming paste. Also, the anelastic strain level at the end of the simulation ($t=40$hours) is almost negligible compared to the other strains (e.g., plastic, thermal, etc.). This result directly ensues from the assumption that the baked ($\xi =1$) ramming paste creep/relaxation behaviour is similar to that of the carbon cathode block (Richard et al., 2006). This case study thus shows the importance of taking all the relevant phenomena including creep behaviour into account in similar problems. A similar analysis could be done for all other deformations (chemical, thermal, plastic, etc.).

Figure 6.

Strain evolution at the center of each component of the model (Reference case).

Figure 7.

Normal stresses at the center of each component of the model (creep case)

5. Conclusion

The main contribution of this work is to demonstrate how to use thermodynamics framework to develop three-dimensional viscoelastic constitutive laws. Two different approaches using internal state variables have been presented. The first one used generic internal state variables related to dissipative mechanisms. A link with rheological model could be created by defining the rheological parameters as a function of the internal states variables. The second approach was similar to the first one except that the internal state variables were based on phenomenology with the use of rheological models, i.e., a viscoelastic one in the present case. To achieve this, a classical 1D rheological model was extended to the 3D case.

The two constitutive laws presented in this work were obtained considering an isotropic material at a reference state, i.e., a state without any influence of external parameters such as temperature. The constitutive parameters evolutions are instead taking into account within the numerical simulation and thus the constitutive laws have to be discretized using an Euler scheme. A case study using the strain based internal state variables approach has been presented. This study shows how creep/relaxation could influence the results of an industrial problem such as the “baking” of a carbonaceous ramming paste.

Finally, as aforementioned, both approaches have considered isotropic materials. However, the methodology to obtain constitutive laws of anisotropic or orthotropic materials such as wood would be similar.

6. List of symbols

$A:$area

$\underset{¯}{\underset{¯}{\left[A_{\alpha }\right]}}:$fourth-order tensor

$$\underset{¯}{\underset{¯}{\left[b_{\alpha }\right]}}:$$fourth-order tensor

$\underset{¯}{\underset{¯}{\left[B_{\alpha }\right]}}:$fourth-order tensor

$c:$speed of sound

$C:$kinematic energy

$\left[d\right]:$second-order strain rate tensor

$e:$specific internal energy

$E:$Young’s modulus

$E_{A_{\alpha }}:$constitutive parameter of $\underset{¯}{\underset{¯}{\left[A_{\alpha }\right]}}$

$E_{b_{\alpha }}:$constitutive parameter of $$\underset{¯}{\underset{¯}{\left[b_{\alpha }\right]}}$$

$E_{B_{\alpha }}:$constitutive parameter of $\underset{¯}{\underset{¯}{\left[B_{\alpha }\right]}}$

$E_E:$constitutive parameter of $\underset{¯}{\underset{¯}{\left[E_{\alpha }\right]}}$

$E_{H_{\alpha }^a}:$constitutive parameter of $\underset{¯}{\underset{¯}{\left[H_{\alpha }^a\right]}}$

$E_{{\eta }_{\alpha }^a}:$constitutive parameter of $\underset{¯}{\underset{¯}{\left[{\eta }_{\alpha }^a\right]}}$

$\underset{¯}{\underset{¯}{\left[E\right]}}:$fourth-order tensor

$\underset{¯}{\underset{¯}{\left[E_{\alpha }\right]}}:$fourth-order tensor

$\overrightarrow{f}:$body load

$F:$force

$G:$instantaneous shear modulus

$G_{\alpha }^{a^{-1}}:$deviatoric coefficent

$\underset{¯}{\underset{¯}{\left[H^e\right]}}:$fourth-order elastic tensor

$\underset{¯}{\underset{¯}{\left[H_{\alpha }^a\right]}}:$fourth-order anelastic tensor

$K:$instantaneous compressive modulus

$K_{\alpha }^{a^{-1}}:$hydrostatic coefficient

$L:$length

$\overrightarrow{n}:$normal vector

$\overrightarrow{q}:$heat flux

$\left[q_{\alpha }\right]:$second-order internal state variable tensor

$Q:$heat

$r:$heat source

$s:$specific entropy density

$S:$entropy

$t:$time

$\overrightarrow{t}:$stress vector

$T:$temperature

$v:$volume

$\overrightarrow{v}:$speed

$W_{ext}:$external work

$\left\{x\right\}:$eigenvector

$\beta :$constant

$\gamma :$constant

$\epsilon :$strain

$\left[\dot{\epsilon }\right]:$second-order strain rate tensor

$\left[{\epsilon }^e\right]:$second-order elastic strain tensor

$\left[{\epsilon }_{\alpha }^a\right]:$second-order anelastic strain tensor

${\zeta }_{\alpha }:$constitutive parameter of $\underset{¯}{\underset{¯}{\left[B_{\alpha }\right]}}$

${\varsigma }_{\alpha }:$constitutive parameter of $\underset{¯}{\underset{¯}{\left[{\eta }_{\alpha }^a\right]}}$

$\underset{¯}{\underset{¯}{\left[{\eta }_{\alpha }^a\right]}}:$fourth-order anelastic tensor

$\lambda :$eigenvalue

${\mu }_{\alpha }:$constitutive parameter of $\underset{¯}{\underset{¯}{\left[A_{\alpha }\right]}}$

${\mu }_{\alpha }:$constitutive parameter of $\underset{¯}{\underset{¯}{\left[H_{\alpha }^a\right]}}$

$\nu :$Poisson’s ratio

${\xi }_{\alpha }:$constitutive parameter of $$\underset{¯}{\underset{¯}{\left[b_{\alpha }\right]}}$$

$\rho :$mass density

$\sigma :$stress

$\left[\sigma \right]:$second-order stress tensor

$\phi :$dissipative potential

$\left[\chi \right]:$matrix of eigenvectors

$\psi :$Helmholtz free energy