InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Chemical Engineering » "Thermodynamics - Physical Chemistry of Aqueous Systems", book edited by Juan Carlos Moreno-Piraján, ISBN 978-953-307-979-0, Published: September 15, 2011 under CC BY-NC-SA 3.0 license. © The Author(s).

Chapter 12

Three-Dimensional Constitutive Viscoelastic Model for Isotropic Materials

By Donald Picard and Mario Fafard
DOI: 10.5772/22485

Article top


Scheme of creep strain (a) and creep strain rate (b) during a constant loading test.
Figure 1. Scheme of creep strain (a) and creep strain rate (b) during a constant loading test.
Diagram of the Hall-Héroult electrolysis cell-from (Richard, 2004).
Figure 2. Diagram of the Hall-Héroult electrolysis cell-from (Richard, 2004).
Diagram of the electrical preheating stage of the Hall-Héroult electrolysis cell- from (Goulet, 2004).
Figure 3. Diagram of the electrical preheating stage of the Hall-Héroult electrolysis cell- from (Goulet, 2004).
Diagram of the three-dimensional visco-elastic rheological model.
Figure 4. Diagram of the three-dimensional visco-elastic rheological model.
Simple finite element model for the carbon cell lining materials – Taken from (D'Amours, 2004).
Figure 5. Simple finite element model for the carbon cell lining materials – Taken from (D'Amours, 2004).
Strain evolution at the center of each component of the model (Reference case).
Figure 6. Strain evolution at the center of each component of the model (Reference case).
Normal stresses at the center of each component of the model (creep case)
Figure 7. Normal stresses at the center of each component of the model (creep case)

Three-Dimensional Constitutive Viscoelastic Model for Isotropic Materials

Donald Picard1, and Mario Fafard2

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):


where E,σ,ε,F,A0,Lf and L0are 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):


where c,ρandν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):


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σufor the concrete, whereσuis 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ν have to be identified.


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/cm2 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 to0.2σu. Linear creep is thus mainly observed at stress levels under0.2σ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:

a,A:Scalars{a},a,a:Vectors[A]:Matrix[A]:Second-order tensor[A]¯¯:Fourth-order tensor

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


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:


whereeis the specific internal energy, Eis the internal energy, Cis the kinematic energy, Pextis the power of external forces and otheris related to all other kinds of energy (thermal, electrical, magnetic, etc.). Also, in equations (7) and (8), vis 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 energye, as described in Eq. (8), and the specific entropy densitys:

while considering

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


where Q/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:


where f is a body load, tis a prescribed stress, ris a heat source, qis a prescribed heat flux and n is a normal vector to a surface, considering that:

dρdt+ρdiv(v)=0mass conservation

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




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


whereSwas already defined in Eq. (10). After some mathematical manipulations, equation (20) becomes:


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:


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ψ. The parameters will be identified later in laboratory. Different approaches can be used to establish the form of ψ 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ψ. 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 [qα] (α=1,...,n). So, Eq. (22) becomes:


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


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:


The Clausius-Duhem inequality thus becomes:


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:


Where the coefficients of [bα]¯¯must be identified. In the present case, the coefficients of [bα]¯¯ will be considered constant, i.e., independent of timet, temperatureT, 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


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:


where [E]¯¯,[Bα]¯¯ and [Aα]¯¯ 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:


The equations (27) to (30) lead to a set of equations for each value ofα:


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:


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 (σ11=σ=constant) of a cylindrical sample is




Coefficients G and K represent the instantaneous shear and compressibility modulus, respectively. In addition to elastic coefficients, one has to identifyGd,Kd ,λ1α parameters andλ2α. Considering the hypothesis that ξα=μα for 3D isotropic materials leads to λ1α=λ2α=λα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[qα]¯¯, 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 [εe] 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 [ε] is reduced to:


[εe]is the instantaneous elastic strain, [εαa]the anelastic strain associated with the Kelvin-Voigt elementα, 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:


where [He]¯¯ is the Voigt form of the fourth order elasticity tensor and [Hαa]¯¯ is the fourth order tensor related to the spring of the α 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,ϕ , per unit of volume, can be chosen as a quadratic form of its arguments:


where the fourth-order tensor [ηαa]¯¯ is related to the dashpot of the α 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:

where φT represents the thermal terms and φM the mechanical terms. For the visco-elastic problem, and assuming that the thermal part, φT, is satisfied, Eq. (41) can be rewritten using Eq. (39), with the mechanical terms onlyas:


where [σ] 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:


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


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


Assuming that [ηαa]¯¯ 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 ψ and ϕ in Eq. (46) leads to Voigt notation (Mase and Mase, 1999),


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


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:




where [Xα]is a matrix containing the eigenvectors {xα} associated with the eigenvalues λα of the following system:


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[Xα]=[X] and thatλ2=λ3=λ4=λ5=λ6. Also, the eigenvalues of Eq. (49) are


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

εaxial=ε11=σ1EHeElastic strainσ1(19α=1NKαa11eλα1tλα1+13α=1NGαa11eλα2tλα2)Creep strain
εradial=ε22=ε33=σ1νEHeElastic strainσ1(19α=1NKαa11eλα1tλα116α=1NGαa11eλα2tλα2)Creep strain

where σ1 is the constant axial stress, ε1is the axial strain and εrthe hoop strain. All other strain components are null. Thecoefficients Kαa1and Gαa1are related to the hydrostatic and deviatoric creep mechanisms, respectively. Moreover, all λαi coefficients are positive and 1/λαiare related to a relaxation time. The elastic parameters E and ν(elastic properties) must be identified through standard compression tests, while the anelastic parametersEHαa,μα ,Eηαa andζα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:


Assuming a forward difference Euler scheme (Reddy, 1993)




Moreover, it can be assumed that


Equation (57) and (59) lead to a system of equation of dimension N×Ndescribed as follows:


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


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 toN=2, the tangent matrix is:




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” ξ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:


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=500W/m2/K) was imposed in this “Anode shadow” region, where the temperature Twas 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=40hours) is almost negligible compared to the other strains (e.g., plastic, thermal, etc.). This result directly ensues from the assumption that the baked (ξ=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α]¯¯:fourth-order tensor

[bα]¯¯:fourth-order tensor

[Bα]¯¯:fourth-order tensor

c:speed of sound

C:kinematic energy

[d]:second-order strain rate tensor

e:specific internal energy

E:Young’s modulus

EAα:constitutive parameter of [Aα]¯¯

Ebα:constitutive parameter of [bα]¯¯

EBα:constitutive parameter of [Bα]¯¯

EE:constitutive parameter of [Eα]¯¯

EHαa:constitutive parameter of [Hαa]¯¯

Eηαa:constitutive parameter of [ηαa]¯¯

[E]¯¯:fourth-order tensor

[Eα]¯¯:fourth-order tensor

f:body load


G:instantaneous shear modulus

Gαa1:deviatoric coefficent

[He]¯¯:fourth-order elastic tensor

[Hαa]¯¯:fourth-order anelastic tensor

K:instantaneous compressive modulus

Kαa1:hydrostatic coefficient


n:normal vector

q:heat flux

[qα]:second-order internal state variable tensor


r:heat source

s:specific entropy density



t:stress vector




Wext:external work





[ε˙]:second-order strain rate tensor

[εe]:second-order elastic strain tensor

[εαa]:second-order anelastic strain tensor

ζα:constitutive parameter of [Bα]¯¯

ςα:constitutive parameter of [ηαa]¯¯

[ηαa]¯¯:fourth-order anelastic tensor


μα:constitutive parameter of [Aα]¯¯

μα:constitutive parameter of [Hαa]¯¯

ν:Poisson’s ratio

ξα:constitutive parameter of [bα]¯¯

ρ:mass density


[σ]:second-order stress tensor

φ:dissipative potential

[χ]:matrix of eigenvectors

ψ:Helmholtz free energy


1 - I. P. Bazarov, 1964 Thermodynamik, Deutscher Verlag der Wissenschaften, Berlin.
2 - O. Coussy, 1995 Mechanics of porous continua, Wiley, 0-471-95267-2, Chichester.
3 - O. Coussy, 2010 Mechanics and physics of porous solids, Wiley, 9780470710388 (online), 9780470721353 (cloth), 0470721359 (cloth), Chichester, West Sussex, U.K.
4 - G. D’Amours, 2004 Développement de lois constitutives thermomécaniques pour les matériaux à base de carbone lors du préchauffage d’une cuve d’électrolyse. Mechanical engineering. Quebec, Canada, Laval University. Ph.D. Thesis, 223 pp.
5 - G. D’Amours, M. Fafard, et al. 2003 "Mechanical behavior of carbon cathode: Understanding, modeling and identification." Light Metals 2003 633-639
6 - M. Fafard, M. T. Boudjelal, et al. 2001 "Three-dimensional viscoelastic model with nonconstant coefficients." (ASCE) Journal of Engineering Mechanics 127 8 808-815 .
7 - P. Goulet, 2004 Modélisation du comportement thermo-électro-mécanique des interfaces de contact d’une cuve de Hall-Héroult. Quebec, Canada, Laval University. Ph.D. Thesis, 265
8 - H. W. Haslach, 2009 "Thermodynamically consistent, maximum dissipation, time-dependent models for non-equilibrium behavior." International Journal of Solids and Structures 462223 3964 EOF 3976.
9 - G. A. Holzapfel, 1996 "On large strain viscoelasticity: Continuum formulation and finite element applications to elastomeric structures." International Journal for Numerical Methods in Engineering 39 22 3903-3926 .
10 - L. E. Kinsler, 2000 Fundamentals of acoustics (4th), Wiley, 0-471-84789-5, New York.
11 - E. Kreyszig, 2006 Advanced engineering mathematics (9th), Wiley, 0-471-72897-7.978-0-471-72897-9, Hoboken, N.J.
12 - G. D. C. Kuiken, 1994 Thermodynamics of irreversible processes- Applications to diffusion and rheology, J. Wiley & sons, 0-471-94844-6, Chichester New York Brisbane.
13 - J. Lemaitre, J. Chaboche, L. , 1990 Mechanics of solid materials, Cambridge University Press, 0521328535, Cambridge.
14 - C. Y. Li, J. Lua, 2009 "A hyper-viscoelastic constitutive model for polyurea." Materials Letters 63 11 877-880 .
15 - A. Lion, C. Liebl, et al. 2010 "Representation of the glass-transition in mechanical and thermal properties of glass-forming materials: A three-dimensional theory based on thermodynamics with internal state variables." Journal of the Mechanics and Physics of Solids 58 9 1338-1360 .
16 - D. Marceau, S. Pilote, et al. 2011 "Advanced numerical simulation of the thermo-electro-chemo-mechanical behaviour of Hall-Héroult cells under electrical preheating." Light Metals 2011 1041-1046 .
17 - G. T. Mase, G. E. Mase, 1999 Continuum mechanics for engineers (2nd), CRC Press, 0849318556 (acid-free paper), Boca Raton, Fla.
18 - R. W. Ogden, 1972 Large Deformation Isotropic Elasticity- On the Correlation of Theory and Experiment for Incompressible Rubberlike Solids, Proceedings of the Royal Society of London, Series A, Mathematical and Physical Sciences, 326 1567 London, February 1972.
19 - D. Picard, 2007 Modélisation et caractérisation du fluage/relaxation de matériaux à base de carbone présents dans les revêtements cathodiques des cuves d’électrolyse de l’aluminium. Quebec, Canada, Laval University. Ph.D. Thesis, 212
20 - D. Picard, G. D’Amours, et al. 2008 "Constitutive laws of carboneous materials of aluminium electrolysis cell: Current knowledge and future development." Light Metals 2008 987-992
21 - D. Picard, M. Fafard, et al. 2008 "Room temperature long-term creep/relaxation behaviour of carbon cathode material." Materials Science and Engineering a-Structural Materials Properties Microstructure and Processing 496(1-2): 366 EOF 375 EOF .
22 - D. Picard, M. Fafard, et al. 2008 "Three-dimensional constitutive creep/relaxation model of carbon cathode materials." Journal of Applied Mechanics-Transactions of the ASME 75, DOI: 10.1115/1.2840044, Artn 031017.
23 - J. N. Reddy, 1993 An introduction to the finite element method (Second edition), McGraw-Hill, 0-07-051355-4, Boston, MA.
24 - D. Richard, 2004 Aspects thermomécaniques de la modélisation par éléments finis du préchauffage électrique d’une cuve de Hall-Héroult : lois constitutives, conception orientée-objet et validation. Quebec, Canada, Laval University. Ph.D. Thesis, 183
25 - D. Richard, G. D’Amours, et al. 2005 "Development and validation of a thermo-chemo-mechanical model of the baking of ramming paste." Light Metals 2005 733-738
26 - D. Richard, P. Goulet, et al. 2006 "Thermo-chemo-mechanical modeling of a Hall-Heroult cell thermal bake-out." Light Metals 2006, 3 Carbon Technology: 669 EOF674 EOF .
27 - R. A. Schapery, 1997 "Nonlinear viscoelastic and viscoplastic constitutive equations based on thermodynamics." Mechanics of Time-Dependent Materials 1 209240 .
28 - M. Sørlie, H. A. Øye, 2010 Cathode in Aluminium Electrolysis (3rd edition), Aluminium-Verlag Marketing & Kommunikation GmbH, 978-3-87017-294-7, Düsseldorf, Germany.
29 - Y. Sun, K. G. Forslund, et al. 2004 "3-D modelling of thermal and sodium expansion in soderberg aluminium reduction cells." Light Metals 2004 587592 .
30 - J. Tessier, C. Duchesne, et al. 2011 "Multiblock monitoring of aluminum reduction cells performance." Light Metals 2011 407412 .
31 - J. Tessier, C. Duchesne, et al. 2010 "Investigation of the Impact of Preheating, Start-up and Early Operation on Potlife." Light Metals 2010 10511055 .
32 - K. C. Valanis, 1972 Irreversible thermodynamics of continuous media- internal variable theory, Springer, 0-387-81127-3, New York.