## Abstract

The behaviour of concrete at elevated temperatures is important for an assessment of integrity (strength and durability) of structures exposed to a high-temperature environment, in applications such as fire exposure, smelting plants and nuclear installations. In modelling terms, a coupled thermomechanical analysis represents a generalization of the computational mechanics of fracture and damage. Here, we develop a fully coupled anisotropic thermomechanical damage model for concrete under high stress and transient temperature, with emphasis on the adherence of the model to the laws of thermodynamics. Specific analytical results are given, deduced from thermodynamics, of a novel interpretation on specific heat, evolution of entropy and the identification of the complete anisotropic, thermomechanical damage surface. The model is also shown to be stable in a computational sense, and to satisfy the laws of thermodynamics.

## 1. Introduction

There is considerable experimental evidence in the literature for the damage which occurs in concrete at elevated temperatures (e.g. Ehm & Schneider 1985; Schneider 1988). It is clear that the elastic modulus, *E*_{0}, ultimate strength (and failure surface in biaxial stress states) and plastic strain are reduced by exposure to high temperatures. A striking feature, which has been shown to be crucial in numerical studies (Khennane & Baker 1992*c*), is that the imposition of a moderate compression during heating actually inhibits the thermal damage, presumably by constraining the micro-crack formation due to differential expansion of paste and aggregates.

Schneider (1988) shows the loss of the Young modulus for concrete tested at different elevated temperatures. While dependent on aggregate type, for most normal concretes the response can be idealized to a parabolic dependence of temperature. For the case of unsealed specimens, the material appears to be completely damaged at *θ*=1000°C, though melting temperature is higher, say, 1300°C. Since for lightweight concretes the reduction is less dramatic, Khennane & Baker (1993) assumed the linear reduction, and obtained reasonable results using a classical plasticity approach.

It is also reasonable to interpret this degradation as damage since it is irreversible. Dias *et al*. (1990) shows a similar plot for mortar, where additional samples were tested after a pre-heat to a given temperature, cooling and heating again to the test temperature. Within reasonable bounds of experimental accuracy, one could draw horizontal lines on the figure intersecting the curve at the pre-heat temperature, i.e. there was no damage recovery. It has also been reported that tests conducted at room temperature after the samples were pre-heated to a given temperature showed similar (reduced) *E*_{0} values to those samples tested while at the temperature corresponding to the pre-heat value: thermal damage is irreversible.

Since thermal damage is likely to be a result of differential expansion of mortar and aggregates—certainly aggregate type is vitally important to the overall measured response—it seems reasonable to develop a two-phase theory of mixtures in order to quantify that damage, akin to the analysis employed by Ortiz (1985) for the inelastic response of concrete at room temperatures. An alternative, exploited in the work of Bazant & Prat (1988), is to consider an activation energy approach to the development of micro-cracks. Again, if carried out on the phases, one would need a reduction in aggregate stiffness to achieve the appropriate concrete damage, but if applied directly to concrete, as Bazant & Prat (1988) did in defining a size-independent fracture energy, one obtains an exponential decay in *E*_{c} with temperature; this will be quantified later when considering thermal softening.

The micro-mechanical (i.e. two phase) formulation following Ortiz (1985) holds considerable appeal from the phenomenological point of view, in that known characteristics can be defined for each phase in developing the macroscopic model. Such features would include: micro-cracking in the mortar; plasticity and thermal softening in the aggregates; different thermal expansions and heat conduction; the effects of phase changes; and the effects of specific aggregate types which are known to have a major influence.

Yazdani & Schreyer (1988), who follow Ortiz (1985), add further inelastic components, and then fit the expressions directly to a concrete model, thus finding a macro-scopic model immediately. This notion will be viable with larger scale problems, with regard to computation, and can be fitted to thermal softening of concrete through a constitutive definition of entropy. Of course, it loses the micro-mechanical spirit, but the identification of properties like fluxes at the global level is more straightforward with global damage. It is the extension of these concepts to a thermomechanical setting that we pursue here.

This paper presents the details of an anisotropic thermomechanical damage model developed for concrete at transient elevated temperatures. The model has been derived within the consistent framework of thermodynamics, drawing on the temperature-independent damage models of Ortiz (1985) and Yazdani & Schreyer (1988) and the thermomechanical coupling aspects of Simo & Miehe (1992). In addition, account has been taken of the known stress–temperature dependence of concrete through the descriptions of thermal and thermomechanical damage, and thermal softening. Explicit expressions have been derived for the free energy including elastic energy, damage due to micro-crack formation, thermal–mechanical coupling (thermal expansion and creep) and thermal energy.

Carol & Willam (1994) suggest a basic deficiency with the stress-induced anisotropy in these models with regard to crack closure. In such cases, they have shown a spurious dissipation under crack closure and rotation of principal stress. Simo & Ju (1987*a*,*b*) suggest an alternative spectral decomposition, but this too is not entirely freed from the problem. Carol & Willam (1994) recommend a non-local approach which alleviates the problems induced by the directionality of the response functions. At this stage, however, it is not clear how significant the problem is. Nonetheless, we believe that a ‘rotating damage model’ is appropriate, and so we utilize the spectral decomposition and response function concepts. Micro-crack closure would occur in local axes, and though the damage tensor may rotate, the irreversibility condition on damage evolution should ensure that no stiffness recovery occurs.

With regard to previous computational modelling in thermomechanics, there have been a number of works investigating thermo-plasticity for concrete at elevated temperatures, e.g. de Borst & Peeters (1989), Khennane & Baker (1992*a–c*, 1993) and more recently, Heinfling *et al*. (1997, 1998). Schrefler *et al*. (1997), Ulm *et al*. (1998) and Consolazio & Chung (2004) have appropriately coupled considerations of vapour diffusion and/or chemical degradation, but without necessarily the strong thermomechanical interaction at high stress and temperature. Other works have concentrated on macroscopic models for thermomechanical damage (or coupled plasticity and damage), albeit generally isotropic, e.g. Baker & Stabler (1998), Stabler & Baker (2000*a*), and general computational procedures, e.g. Simo & Miehe (1992) or Stabler & Baker (2000*b*).

This paper is concerned with anisotropy in a fully coupled thermomechanical setting, and particularly macroscopic models for extreme conditions of high stress and transient temperature. Our interest centres on the formulation of the models and, in accord with the theme of this special issue, their adherence to the laws of thermodynamics.

As such, we begin with an overview of a general framework of thermodynamics, and follow this with a discussion of isothermal (and room temperature) models for stress-induced anisotropy in concrete. We then present details of a coupled anisotropic thermomechanical damage model, and investigate the model's behaviour with regard to dissipation and the second law of thermodynamics. Brief details are given for the stability of computational schemes, again with regard to the dissipation inequality.

## 2. Dissipation and coupled thermomechanics

### (a) General framework of thermodynamics

Complete details of a general framework of thermodynamics can be found in the original internal variable representations of Coleman & Gurtin (1967) and Lubliner (1972), and from more recent developments in damage mechanics (Krajcinovic 1989; Mazars & Pijaudier-Cabot 1989), plasticity (Maugin 1992), coupled plasticity and damage (Hansen & Schreyer 1994), and thermoplasticity with a distinct computational bias (Simo & Miehe 1992; Armero & Simo 1993). Here we wish to outline that framework, motivated by our objective to develop a fully coupled damage model for the thermomechanical analysis of a two-phase composite like concrete.

Armero & Simo (1993) summarize the general form of the local balance laws (momentum and energy balance) for small stress and strain in a domain *Ω*×[0,*t*] as(2.1)where (., *t*) is the motion of particles in a continuum defined on a reference configuration *Ω* in a time interval [0,*t*], with velocity field * V*=∂/∂

*t*, deformation gradient defined by , nominal stress tensor

**, associated temperature field**

*σ**Θ*(.,

*t*), entropy

*η*, density

*ρ*

_{0}, heat flux vector

*, with*

**Q***and being body forces and heat source, respectively.*

**B**The internal dissipation is given by(2.2)where the three terms would be interpreted, respectively, as stress power, thermal power and change in internal energy; this latter should here include contributions from the free energy to create new micro-crack surfaces (Ortiz 1985), as well as for the creation of phase changes. Now we can write both the internal energy and heat flux in the form(2.3)where * D*=

*(*

**D****,**

*σ**Θ*) is a set of micro-structural variables, such as plastic strains, hardening parameters or damage variables, which define the mechanical state of the material. These are typically described in constitutive models in the form of rate equations (i.e. evolution laws)(2.4)

For an analysis of concrete, where thermal damage is known to be partially inhibited by compressive stress, the definition of this latter evolution law is very important. Within the normal context of thermodynamics, we find the nominal stress, absolute temperature and thermodynamic ‘forces’ (Maugin 1992) through derivatives of the internal energy(2.5)

The chain rule applied to internal energy yields(2.6)

Accordingly, equation (2.6) with (2.5) in equation (2.2) yields(2.7)

### (b) Free energy and conjugacy

Free energy—that which can be released to the system—is found from internal energy by subtracting the heat lost due to entropy production. The functional form of the Helmholtz free energy, *ψ*, is thus obtained from the internal energy through the Legendre transformation(2.8)

From equation (2.8), the rate change of free energy is(2.9)which when introduced into equation (2.2) gives(2.10)

For completeness, we deduce internal dissipation from constitutive laws like equation (2.5) based on the free energy function. Firstly, from equation (2.8) it is clear that(2.11)where relation equation (2.11*b*) is conjugate to relation equation (2.5*b*). The chain rule applied to free energy yields(2.12)which with constitutive laws equation (2.11) substituted into equation (2.10) yields(2.13)as in equation (2.7).

### (c) Gibb's energy and complementarity

To date, we have outlined thermodynamic relations based on internal and free energies written as functions of strain. For later developments, we now consider the complementary problem where energy is written as a function of stress. The functional forms of complementary internal energy and Gibb's free energy are defined as(2.14)

Beginning with the definitions of stress and strain, based on the stored energy functions(2.15)we write the complementary energies per unit volume as(2.16)

From this relationship, we find transformations like equation (2.8)(2.17)and the conjugate relations(2.18)

We reflect briefly on a comparison between equations (2.11*b*) and (2.18*b*) for concrete. For a problem described by Helmholtz energy and stiffness, which decreases with temperature in concrete, then the negative sign in equation (2.11*b*) ensures an increase in entropy with temperature. For a definition based on Gibb's free energy and compliance, which correspondingly increases with temperature, one again finds an increase in entropy from equation (2.18*b*).

From equations (2.8) and (2.16), we find the relationship between Gibbs energy and internal energy(2.19)

This reflects the origins of Gibbs energy in enthalpy, and its corresponding reliance on pressure–volume. In solid mechanics, pressure–volume is interpreted as stress–strain, with appropriate definitions on representative volumes. With proper sign conventions for stress, Gibbs energy in equation (2.19) can in fact be derived directly from enthalpy, and the Legendre transformation equation (2.17).

Finally, we record the constitutive laws(2.20)

The chain rule applied to Gibb's energy yields(2.21)

The alternate forms of the dissipation rate (from equations (2.1), (2.16) and (2.17))(2.22)both of which again reduce to(2.23)

### (d) Entropy, dissipation and the second law

The second law of thermodynamics requires that the rate of change of entropy should not be less than the production of entropy from heat (Maugin 1992), which gives rise to the dissipation inequality(2.24)where the dissipation due to conduction is(2.25)

Further, we split the internal dissipation into parts associated with the thermal production of entropy, and the contribution to dissipation of the purely mechanical theory: . Since the second law must be satisfied under non-conductive and isothermal processes, one requires separately that and , where the mechanical dissipation under isothermal conditions follows from Clausius–Duhem form (equation (2.2) with (2.8))(2.26)

The general problem of maximum dissipation leads to the evolution equations for an associative process in plasticity or damage (Simo & Miehe 1992; Hansen & Schreyer 1994). Thus given a smooth damage function, , extension of the isothermal principle of maximum dissipation to a thermomechanical version gives the evolution of entropy as(2.27)where the consistency parameter, , is the Lagrange multiplier satisfying the Kuhn–Tucker conditions: .

The flow rule equation (2.27) is properly interpreted as that controlling the evolution of thermal softening. The consequence is that the thermal dissipation can be written(2.28)

We note from equation (2.27) that for a purely mechanical (temperature independent) process. Generally, of course, is only zero for an *isentropic* process not an *isothermal* process , because entropy can be increased when the mechanical process is temperature dependent even when temperature does not change, as evident from equation (2.27).

### (e) Interpretation of entropy and internal variables

In order to account for the additional configurational entropy associated with dislocation and lattice defect motion, Armero & Simo (1992) split the total entropy into ’elastic’ and ’plastic’ parts(2.29)with micro-structural variables correspondingly being plastic strain, hardening variables, and plastic entropy. Further, with the internal energy equation (2.3) written with a functional dependence on the elastic strains and *η*^{e}, and a hardening potential, one finds more specific constitutive equations from equation (2.4)(2.30)where the * D* include

*η*

^{p}. As a result, the meaning of thermal dissipation would be refined to(2.31)with evolution of thermal softening as(2.32)

In a damage model, one could use *η*^{p} conceptually in an identical fashion to govern thermal softening of the damage surface, for which we will reserve the same notation . Indeed, for a damage model of concrete, where mechanical damage is promoted through micro-cracking, we will note later a clear relationship between plastic entropy and the (free) surface energy of micro-crack formation, and thence softening as this energy changes with temperature.

### (f) Temperature evolution

The first step in seeking a solution to the governing equation (2.1) is to transform the energy balance equation into a first order equation in temperature evolution, with displacement and temperature as state variables.

Based on a Gibb's free energy formulation, it is easily shown from an application of the chain rule that(2.33)where the last follows from equation (2.22). Then, after substituting equation (2.33*c*) into equation (2.1*c*), one finds(2.34)where(2.35)

For specific damage models the first term on the right hand side of equation (2.34) reduces to explicit derivatives of the damage surface with respect to temperature and internal damage variables. Of course, in the absence of mechanical effects, equation (2.34) reduces to the standard heat conduction equation. Simo & Miehe (1992) use the entropy split equation (2.29) and evolution equation (2.32) to arrive at the same form based on a Helmholtz free energy.

The thermal–mechanical coupling is completed once a functional form for the flux, * Q*, is defined. For a purely isotropic process, one would use the classical form(2.36)

Even though thermal damage, in the absence of load, may be considered isotropic, there is an induced anisotropy in the conductivity coefficients as a result of the coupling between temperature and damage induced anisotropy. We thus write a general tensorial form(2.37)where * Q* is the flux vector and

*k*is the second order heat conduction tensor.

## 3. Isothermal damage models

### (a) Two-phase model of Ortiz

Ortiz (1985) developed a damage model for the inelastic response of concrete where damage was deemed to accrue through micro-crack formation in the mortar. Damage was applied directly to the compliance of the mortar, with the flow of damage related to response functions for tensile and compressive modes of fracture (i.e. cleavage and splitting modes) which are dictated by the orientation of micro-cracks. The aggregate phase is treated with a relatively simple plasticity model, with the composite response being determined through a theory of mixtures.

Beginning with the Gibb's free energy, for isothermal processes, Ortiz (1985) writes(3.1)where ** σ** is the small deformation stress tensor, is the fourth order compliance tensor, and

*A*

^{c}is the free energy to form micro-cracks. Then one finds strains and strain rates as(3.2)and, in the normal fashion(3.3)where the two terms correspond to rate of elastic and inelastic deformations, respectively. Flexibility is then written in additive form(3.4)where the first term represents the elasticity tensor in the absence of micro-cracks, the second the flexibility due to the active micro-cracks, and is regarded as the damage tensor, with

*μ*being the damage variable. Thus, the flow rule for the inelastic strain, , in equation (3.3) requires the identification of an evolution law for the damage tensor in the form(3.5)where is a fourth order tensor governing the direction of damage.

To evaluate , Ortiz (1985) initially finds the flexibility assuming all micro-cracks are open, , and projects this onto the positive cone of the stress tensor. Noting the cross-effects on micro-cracks under compressive splitting stresses, Ortiz (1985) uses an additive form for the cracked flexibility, with the compressive part being projected onto the negative cone, thus making use of the full spectral decomposition of the stress tensor into its positive and negative parts(3.6)

This latter is carried out by finding the eigenvalues, *σ*_{α}, and eigenvectors, , of the stress tensor, and summing for all positive eigenvalues(3.7)

The operator, *P*^{+}, which projects the stress onto its positive cone is given by the relations(3.8)(3.9)and projection to the negative cone given by *σ*^{−}=*P*^{−} : ** σ** with

*P*^{−}=

*I*^{4}−

*P*^{+},

*I*^{4}being the fourth order identity tensor. Hence, the rate change of damage tensor is given by(3.10)where are the parts of the evolution for the damage tensor related to the tensile and compressive failure modes, respectively, assuming all micro-cracks open. The free energy can then be conveniently rewritten in terms of the two modes of micro-cracking as(3.11)

The rate independent damage rules are written for the two modes of cracking in terms of the response functions mentioned earlier(3.12)which for an associated process are shown to be(3.13)where *ω* is a parameter representing the relative strengths in uniaxial compression and tension. With equations (3.5), (3.12) and (3.10), we also find(3.14)

We restate that in evaluating the ‘damage flexibility’, one actually evaluates the damage tensor of equation (2.4) directly. Given, then, the specific notation of this section we observe the equivalence relations on the notations(3.15)

Finally, after applying the principle of maximum dissipation, one finds from equations (3.11), (3.12) and (3.13) a damage function (Ortiz 1985)(3.16)where *t*(*μ*) is denoted the critical stress for the onset of damage given by(3.17)

Any legitimate tension softening stress–strain relationship can be used to determine an explicit form for critical stress in terms of cumulative damage. For a given value of critical stress, the equation (*μ*)=0, given by equation (3.15), describes a circle in both biaxial tension and compression, and an ellipse in the tension–compression. Since the limiting value of critical stress is the tension strength, *f*_{t}, it is clear from the conditions of uniaxial compression that: .

Ortiz (1985) considers flexibility along the axis as 1/*E*_{0}+*μ*. Then by equating the general stress–strain relation to a specific form such as that from Smith & Young (1955)(3.18)where is the strain at peak stress, one finds the tangent stiffness and critical stress to be(3.19)

For later use, we note also the derivative(3.20)which is clearly negative beyond peak strain.

### (b) Composite model of Yazdani & Schreyer

Schreyer (1989) and Yazdani & Schreyer (1988, 1990) developed a coupled plasticity and damage model for concrete at room temperatures, where the damage part uses a similar formulation to that which Ortiz (1985) used for mortar, but applied to the composite concrete. Their damage surface thus shows a number of modifications from equation (3.16) to account for the observed macroscopic behaviour of concrete. For example, they include a direct pressure term to account for positive dilatancy, modified stress cones to inhibit damage under hydrostatic compression, and deviatoric stress dependency to account for inelastic flow strains.

Those authors show a strain rate decomposition, which follows an integration of the Gibb's free energy, into elastic, damage and inelastic strain rate components(3.21)where(3.22)

The damage strains, , thus correspond to the recoverable strain due to micro-crack formation, whereas the inelastic strains, , represent the irrecoverable part. Depending on the relative magnitude of these two components, one can vary between purely brittle and purely plastic processes: this relative magnitude of strain components as a function of temperature is unfortunately almost completely undocumented, except that concrete does appear to be very ductile at high temperatures with a commensurate increase in flow strains.

Schreyer (1989) used the same rationale as Ortiz (1985), except that the response function for micro-cracking under compressive load is comprised of two parts. The first, which corresponds to the cross-effect, uses the shifted negative cone given by(3.23)where *λ*_{max} is the maximum eigenvalue of *σ*^{−} introduced in order to preclude damage under hydrostatic compression (Schreyer 1989), and *i* is the second order identity tensor. That is,(3.24)where *ω* takes the same meaning as before. The second part, , is a semi-empirical term included to account for the effect of lateral pressure(3.25)where *α* is a factor which accounts for the effect of lateral pressure, and and *H*(.) are the smallest eigenvalue of and the Heaviside function, respectively.

Yazdani & Schreyer (1988, 1990) postulate that imperfect fracture only occurs with the compressive mode of cracking, and thus relate the inelastic strain rate to the positive and negative cones of the deviatoric component of the stress tensor. In the spirit of equation (3.3) with flow rules equation (3.12*b*), they write an evolution equation for permanent strains as(3.26)where the material parameter *β*>1 for permanent deformation, and the tensors, *s*^{+} and *s*^{−} are the projections of the deviatoric stress tensor on the positive and negative cones.

The resulting damage surface takes the form(3.27)where for convenience we have defined(3.28)

The work of Yazdani & Schreyer (1990) is general for three dimensional states of stress where the pressure is the usual hydrostatic term: . If one uses the form of equation (3.27) for two-dimensional problems, with the corresponding trace function for pressure, , then it can be shown that there is no solution to (*μ*)=0, this being complex.

The use of the shifted cone is certainly effective, but produces a sharp point in hydrostatic biaxial stress (Baker & de Borst 1995) which is not always representative; we note from Yazdani & Schreyer (1990) that this is not so in three dimensions.

Baker & de Borst (1995) show a more rounded surface obtained by using the response function in equation (3.13*b*) with the negative cone of *σ* as in Ortiz (1985), and not the shifted cone in equation (3.23). This is more appropriate in shape for normal concretes over a wide temperature range, though for lightweight concrete, the beneficial effect of confinement as with equation (3.23) is applicable. Indeed, since we are not here using phase descriptions following Ortiz (1985), one might to an extent control the influence of aggregate type by using the generalized shifted cone(3.29)

Furthermore, one notes that the biaxial failure surface for concrete obtained isothermally at 300°C is certainly not convex (Khennane & Baker 1992*c*), and that the maximal benefit of confinement is obtained between the uniaxial and hydrostatic compression states. That is, the addition of some lateral pressure to a uniaxial compression test inhibits damage, but too much increases it relatively. This feature can be retrieved from equation (3.27) by exploiting the term: *ω*(** τ**:

*σ*).

The factor *β* was included to generate permanent deformations. By changing *β*, damage surfaces using equations (3.13) and (3.14) can be modified to generate the shape corresponding to the 300°C biaxial data of Ehm & Schneider (1985). We note that since *s*^{+} and *s*^{−} represent decompositions of the deviatoric stress, they are null under hydrostatic compression, so that the value of *β* is irrelevant. An analytical expansion of this term as a function of principal stresses, *σ*_{1} and *σ*_{2}, reveals(3.30)which is quadratic between the points of uniaxial compression (*σ*_{1}=0) and hydrostatic compression (*σ*_{1}=*σ*_{2}<0), and *β* controls the amplitude. It can thus be seen that for uniaxial compression, the terms in *β* are again irrelevant, being multiplied by the zero principal stress.

## 4. A thermomechanical damage model

There are many models for thermal damage which can be embedded within the basic framework of thermodynamics and damage. Here we outline just one model for degradation of elastic constants and thermal softening of the limiting damage.

### (a) Gibb's free energy

The formulation of the constitutive equations could begin by proposing that strain rate be decomposed, and then invoking a Hookian law for the details. However, following Yazdani & Schreyer (1988), we begin with the relationship between compliance and stress tensors, and the Gibb's free energy(4.1)

Integrating once with respect to ** σ** yields the total strain decomposition(4.2)where is the sum of the free thermal strain and creep strains, both isothermal and transient. Written in rate form we recover an expression like equation (3.21)(4.3)where the total inelastic strain rates are given by(4.4)

A second integration yields the free energy function(4.5)where the Gibb's (free) energy is shown as functionally dependent on the stress, damage tensor and temperature. Flexibility is given by equation (3.4), represents the sum of thermal strains, and *T*(*Θ*) is a function of temperature only, being a constant of integration.

### (b) Elastic entropy

From the form of the Gibb's free energy equation (4.5), we find an explicit form for the elastic entropy(4.6)

Intuitively, one might think of the change in micro-crack surface energy with temperature as indicative of thermal softening, in which case the last term would relate to plastic entropy. This concept is now explored a little further. From the damage models of §3*a* we know that the damage surface, =∂*G*/∂*μ*, includes a critical stress term: *t*^{2}(*μ*)−2∂*A*^{c}/∂*μ*. Thus the flow rule for plastic entropy equation (2.32) will include a cross derivative term(4.7)from which the explicit relation to thermal softening is evident: given the relationship between the softening law for critical stress and material properties, *f*_{t} and *E*_{0}, the temperature derivative above represents the thermal softening of those properties. The actual form of the softening will be investigated later.

The temperature derivative of inelastic strains in equation (4.6) will also influence the plastic entropy, so that the entropy split equation (2.29) could be identified if required. Equally, once the details of the parameters in equation (4.6) have been defined, the explicit constitutive law for *η*^{e} would follow. The only bar to an explicit evaluation of equation (4.6) is the expansion of the final term in *A*^{c}. One can begin with the relationship between surface energy and critical stress and integrate equation (3.19*b*) to find a functional form for *A*^{c}, which could be differentiated with respect to *Θ* once the softening law has been defined. That process leads to, *A*^{c}=*A*^{c}(0, *Θ*)+*A*^{c}(*μ*, *Θ*), where the first term is a constant of integration at zero damage. It is, however, a function of temperature and without experimental evidence to guide a functional form, there seems little point in pursuing that line.

For a fully coupled solution, it is more realistic to integrate the evolution law for entropy in an incremental fashion, since that flow rule will be well defined.

### (c) Evolution of thermal damage

Given our understanding of damage under increasing temperature, we propose a temperature dependence of the compliance which is incorporated within the response functions of §3*a*. We will evaluate the critical stress as Ortiz (1985), but assume all increase in flexibility, thermal or mechanical, to be attributed to the damage variable *μ*. Thus equation (3.4) reads(4.8)

The central concept is to increase the damage compliance by a temperature dependent term. One option would be to scale the elastic compliance by a function of temperature, determined empirically to satisfy thermal degradation. However, damage would then occur under any state of stress and ignores the restraining effect of compression. Alternatively, one might scale flexibility by a function of volumetric strain, accounting for stress and pore pressure. However, while this model should be studied in detail, there are a few deficiencies: (i) damage would be isotropic under any stress; (ii) it would be difficult to quantify the effect from current data, though a viable approach might be to re-cast the data for free thermal strain for different initial stresses, and relate these to the decrease in the Young modulus; (iii) the partial derivatives with respect to temperature, in equation (4.6) for example, would be awkward. It seems better then to propose an evolution law like those in equation (3.12). Given the inhibiting effect of compression, and the desired anisotropy under mixed states of stress, we suggest an additive term in the same form as equation (3.12*a*)(4.9)where the function *f*(*Θ*) should capture experimental results. Note that even when *σ*^{+}=**0**, damage still occurs mechanically under large compression and thermal softening so that the observed beneficial effect of reasonable initial stress is captured.

To define *f*(*Θ*), we might suggest a linear decrease in the Young modulus, though this is difficult to expand into the additive form of equation (3.4). Instead we approximate(4.10)

In order to preserve the irreversibility assumption, we in effect adopt a simple damage surface for thermal damage (analogous to the basic isotropic surface in Mazars & Pijaudier-Cabot 1989) where damage only accrues once the maximum temperature at a point has been exceeded. Hence we substitute the maximum temperature, *Θ*_{max}, for *Θ* in the above so that no stiffening occurs on local cooling. Finally, since the term like *O*(*Θ*^{2}) may be excessive for a cumulative damage model, we modify the maximum temperature by a factor, 0<*γ*<1, where *γ* should be determined from numerical experiment. Thus, on combining equations (3.12*a*), (4.9) and (4.10), we have a modified compliance for the cleavage mode of micro-cracking(4.11)

It might be argued that since the second response function governs splitting cracks, it should also have a temperature dependence. However, we prefer to induce that dependence through the thermal softening function, *h*(*Θ*), defined later. It should be kept in mind that the softening function will not only control the contraction of the damage surface, but also the growth of damage (and hence compliance) given the softening of the critical stress. That is, damage here is related to strain, and at zero damage *t*=0, so that damage accumulates from the outset. Thus, even in cases of no stress during heating, the critical stress is softening, and when stress is applied damage will accumulate rapidly, giving the desired reduction in stiffness.

For the evolution of damage, we sum the thermomechanical response function in equation (4.11) to the cross effects given by equation (3.24), and a form of equation (3.25) modified to reflect the biaxial nature of this work(4.12)

Returning again to the inhibiting effect of hydrostatic compression, it should be noted that the third response function only exists when *σ*^{−} is not null, and the effect of reducing its magnitude is to expand the damage surface markedly in biaxial compression. The model we suggest includes a beneficial temperature dependence applied to the third response function, written as(4.13)where the Heaviside functions convey the meaning that the benefit only applies under simultaneous heating and compressive loading. In addition, it will be noted that in such cases, the function *g*(*Θ*)=*h*^{2}(*Θ*) so that we essentially negate the softening applied to the critical stress, except that the second response function in biaxial compression has no temperature dependence and hence the overall result is an inhibition of damage, not cancellation.

### (d) Damage function

Both Ortiz (1985) and Yazdani & Schreyer (1988, 1990) define a suitable damage function as the derivative of the Gibb's energy with respect to the damage variable. For the general thermomechanical case, we wish to show that the more formal statement of internal dissipation leads to the same conjugate form here. Given that equation (3.21) includes permanent (flow) strains in addition to damage strains, these should be included in the set of micro-structural variables which govern the mechanical state of the material, in addition to the damage tensor itself; for the flow of permanent strains, we use the form of equation (3.26) precisely. Thus we have(4.14)

From the general statement in equation (2.23), one finds(4.15)whence we write(4.16)

Evidently, we can again argue that the two possible states which also satisfy the Kuhn–Tucker conditions are(4.17)

The evolution of damage in equation (4.16) is found using the incremental form of equation (3.5), Δ*C*=Δ*μ R*(

*σ*), where the three response functions (4.11, 3.24 and 4.13) are added. It should be noted that since the third response function is hydrostatic, the original projection onto the negative cone (as applied to equation (3.25)) becomes the identity:

*P*^{−}:

*I*^{4}:

*P*^{−}=

*I*^{4}. Thus, we write(4.18)

Thence, on combining equations (4.16) and (4.18) with equations (4.11), (3.24), (4.13), we obtain the following explicit form for the damage function(4.19)

### (e) Thermal softening

Thermal softening is here restricted to softening of the critical stress equation (3.19*b*) with temperature. Given the form of equation (3.19*b*) one may seek to vary both *f*_{t} and *E*_{0}. However, the degradation of elastic constants has been accounted for in the response function equation (4.9), and need not be included here; this we recall is consistent with the notion that *μ* is responsible for all departure from the initial elastic tangent at room temperature, *E*_{0}. If the strain at peak stress, , is also held constant, then it is clear from equation (3.19*a*) that the tangent the Young modulus follows the same softening law as *f*_{t}. The empirical variations in strength and stiffness have different forms, and it would therefore be possible to accommodate these by an appropriate functional form for . In so doing, there would seem to be little need for the temperature dependence in equation (4.9), but one would return to the problem of thermal damage being independent of stress state, unless the softening law became very complicated indeed.

Hence, we retain equations (4.9) and (4.10), and define a softening law for *f*_{t} only. Then, the effective stiffness will have a different form from *f*_{t} given the contribution in both equations (3.19*a*) and (4.9). The approach adopted is to form a peak tensile strength that mirrors activation energy concepts, and a compressive strength that follows empirical observations. To begin, we define thermal softening on peak stress, and hence critical stress, in multiplicative form(4.20)where *h*(*Θ*) is the softening function.

The reduction in is essentially parabolic, with little loss up to 200°C (Schneider 1988). Under uniaxial compression, *σ*^{+}=**0**, the damage function for the case takes the form: . Thus, to achieve the appropriate parabolic loss, equation (4.20) suggests that the softening function should take the same form(4.21)given that at 1000 °C all strength has been lost; this curve is shown in Baker & de Borst (1995).

However, because of the temperature dependence of the first term in equation (4.19), the limiting damage in tension will not have the same form. Specifically, for uniaxial tension, the damage function reveals that(4.22)

The resulting limiting damage surfaces obtained with no stress during heating can be compared with experimental results in Schneider (1988). The effect of the parabolic softening is to crowd the envelopes together at moderate temperatures, which is generally the case for normal strength concrete; for lightweight concrete a linear softening may be more appropriate. Of course, one can vary the shape of the envelopes as pointed out earlier by varying the factors: *ω*, *β*, *λ*.

Finally, with these definitions, the evolution of (plastic) entropy is completely defined. Given the flow rule for thermal softening equation (2.32), we now examine the required non-negativity of the thermal dissipation equation (2.31), through explicit differentiation of equation (4.19). This yields(4.23)for the case of no initial stress during heating. This is guaranteed non-negative since we ensure that thermal softening of the critical stress decreases monotonically, . In deriving equation (4.23), our model assumes a constant *ω*. If, however, the changing shape of the limiting damage surface with temperature is realized through a definition *ω*=*ω*(*Θ*), then the second, third and fourth terms in equation (4.19) would be represented in equation (4.23), and would require careful definition of in order to satisfy the second law. The same is true for *β* which, as noted earlier, changes the curvature of in stress space.

For the complete expression in equation (4.19), one finds after some manipulation(4.24)which is again non-negative since ≤0 and ∂*h*/∂*Θ*<0. Hence, the dissipation inequality is preserved and the model is thermodynamically consistent.

### (f) Consistent tangent operator

From equation (3.21), we naturally wish to develop the consistent tangent stiffness defined by(4.25)where **d**_{T} is the fourth order tangent stiffness tensor, and the mechanical strain component. In fact, if we cast our expressions for mechanical strain rate in a form like equation (3.21), and particularly if we can write an evolution form for the total inelastic strain rate(4.26)where is an effective stress, then we can follow the derivation exactly as given in Ortiz (1985), without involving the complexities of the damage function through intermediate stages.

From the basic definition of inelastic strain rates in equation (4.4), with equations (3.5) and (3.26), the effective stress is found immediately(4.27)

For an associated model, one would expect the evolution of inelastic strain to be(4.28)which from equation (4.19) does in fact lead to as defined in equation (4.27).

We now invoke the consistency equation, assuming that mechanical damage occurs ‘instantaneously’ with no temperature change(4.29)which also follows from equation (4.19). For the local iterations for flow around the damage surface, we use equations (4.4) and (4.26) to rewrite equation (4.3) as(4.30)or(4.31)with ** d**=

*C*^{−1}. After taking the inner product of with equation (4.30) and substituting into equation (4.29), the damage variable is found from(4.32)whence with equation (4.31) one finds(4.33)

## 5. Thermal evolution

### (a) Specific heat and thermal energy

The thermal energy term in equation (4.5) arose essentially as a constant of integration. For the case of a constant specific heat, *c*, Simo & Miehe (1992) show that the additional term, *T*(*Θ*), must be included in the free energy function in order to satisfy the second derivative(5.1)

In fact, equation (5.1) arises as a consequence of the second law, in that the second variation must be non-negative, so that(5.2)

But we must be careful to interpret the coefficient, , correctly. When the free energy terms apart from *T*(*Θ*) are at most linear in *Θ* (as in Simo & Miehe 1992), then equation (5.2) applies and is a constant (=*c*). A straightforward integration of equation (5.1) then yields a free thermal energy term like(5.3)

For the damage models postulated here, *T*(*Θ*) is still found from equation (5.3) even though the thermomechanical terms are higher order in *Θ*, since it is a function of temperature only. The temperature derivative required by equation (4.6) is thus(5.4)

For the free energy forms used here equation (4.5), and in general, as found from equation (5.2) is not constant. This is essential. Stabler & Baker (2000*a*) showed explicitly for the case of thermo-elasticity with isotropic damage that specific heat could not be constant else the non-negativity of dissipation would be lost. They wrote specific heat as a semi-empirical form of thermal damage parameter, and found that a small variation from a constant *c* was sufficient to restore the dissipation condition.

Hence, it is in equation (5.2) that should be properly interpreted as specific heat. Since the evaluation of requires derivatives of thermal energy *T*(*Θ*) as well as fracture energy *A*^{c} in equation (5.2) with equation (4.6), we can retain the notion of a constant coefficient *c* as the base constant part of specific heat, with a temperature/damage variation coming from the remaining derivatives of *A*^{c}(5.5)

Thence, equation (5.3) is still seen an appropriate form of thermal energy, and equation (5.5) matches the result in Stabler & Baker (2000*a*).

Even in the absence of stress, the evaluation of equation (5.2) yields a non-constant value for , being dependent on thermal damage through the creation of micro-cracks, i.e. *A*^{c}. This seems reasonable for concrete based on empirical evidence. In concrete, heat capacity is not constant (Schneider 1988), except for lightweight concrete, and appears generally to vary parabolically from to 1.2 kJ kg^{−1} K^{−1} in quartzite concrete, albeit the influence of the term does not appear dramatic in simulations Schellekens & Parisch (1994).

### (b) Thermal conductivity

In detail, one would expect the level of damage to influence the coefficients of thermal conductivity. In fact, the influence could be assumed to be small until a macro-scopic crack opens, at which point conductivity across the crack reduces markedly; that parallel to the crack remains unchanged within the concrete but clearly conduction due to convection, either through water diffusion or the air space of the crack, would increase.

For an isotropic damage process, one would write conductivity as a function of effective strain and temperature: . Certain Codes of Practice, e.g. EC2, give conductivity as a function of temperature in the form(5.6)

Of course, this expression embodies thermal damage, but takes no account of mechanical damage. This might be modified to account for a mechanical damage variable, 0≤*D*_{mech}≤1, so as to recover equation (5.6), or alternatively written in a power law format. However, the difficulty is that conductivity along a crack may indeed increase, whereas equation (5.6) only permits a decrease. It may then be appropriate to restrict the flux to an orthotropic form of equation (2.37), requiring one coupled conductivity coefficient for each coordinate direction, as this may be amenable to experimental verification.

However, we wish to develop a flux linked to anisotropic damage via the decomposition of stress and strain, so that it reflects the directionality of stress-induced damage, while still retaining the orthotropic nature of coefficients involved.

Whether stress or strain is used to control evolution of thermomechanical damage, the thermo-hygral coupling model suggested here could equally be linked to a strain dependence. We write the spectral decomposition in the form(5.7)

Simo & Ju (1987*a*,*b*) write response functions for strain corresponding to the cracking and splitting modes of anisotropy, similar to the stress-based functions in equation (3.13). The projection of these fourth order tensors, again in similar fashion to the analysis in §3*a*, is used to yield second order flux tensors for the positive and negative cones respectively. Writing *ϵ*_{α} and as the eigenvalue and eigenvector of the strain tensor, we write(5.8)

Finally, we write the tensorial form of thermal conductivity as(5.9)where the coefficients *k*_{I} and *k*_{II} refer to conductivity along crack paths and orthogonal to cracks respectively. With regard to the definition of the coefficients, the first, *k*_{I}, should be a decreasing function of positive strain, which reverts to equation (5.6) in the absence of mechanical strain. That is, in cases of zero stress we have only (hydrostatic) thermal expansion, whence **ϵ**^{−}=**0** and **ϵ**^{+}=**ϵ**=*f*(*Θ*)** i**. Then equation (5.9) in equation (2.37) yields the isotropic form like equation (2.36) as required for pure heating. The second coefficient,

*k*

_{II}, should be an increasing function relating to crack width, naturally determined from experiment.

## 6. Stability of the IBVP

Using the governing balance equation (2.1), with the temperature evolution equation (2.34), we write the initial boundary value problem (IBVP) using an operator split rationale (Simo & Miehe 1992, Stabler & Baker 2000*b*)(6.1)where ** B**=

*ρ*

_{0}

**from equation (2.1) and**

*b**u*is the (small) displacement vector.

We defer details of the linearization of weak forms and specific computational details like Stabler & Baker (2000*b*) to another article. Here we wish to derive an a-priori stability estimate for the nonlinear IBVP, analogous to that given in Armero & Simo (1993). We begin with thermal boundary conditions, which are chosen without loss of generality to be *Θ*=*Θ*_{ref} on *Γ* which implies that flux is zero on the boundary: * Q*=−

**k**.∇

*Θ*=0 on

*Γ*. We begin with a total energy functional, which can be written in terms of internal energy (accounting for background heat at the chosen thermal boundary conditions), kinetic energy and potential energy of the loads. From this follows the canonical form of free energy, written in a Gibb's form for consistency with earlier sections(6.2)where

*Π*

_{ext}is the potential energy of the external loads. Given boundary tractions

*=*

**T****.**

*σ***on**

*n**Γ*, where

**is the vector normal to the boundary, and boundary values**

*n*

*u*^{*}=

**;**

*u*

*v*^{*}=

**on**

*v**Γ*, we can write the potential energy and its time derivative as(6.3)

Now, the time derivative of the functional, , can be written(6.4)

To reduce equation (6.4), we make use of the constant thermal boundary conditions, and note that the rate of change of Gibb's energy (from equation (2.21)) is(6.5)

Then, after introducing the momentum balance equation, scaled by ** v**, equation (6.4) becomes(6.6)

We now apply Green's Theorem to the last term in the domain integral, and introduce the energy balance equation, assuming that there is no heat source, to obtain(6.7)using the definition of in equation (2.13).

The integral of div(* Q*) in equation (6.7) vanishes after application of the divergence theorem by virtue of the chosen boundary conditions. For integration of the third term, we note first the identity(6.8)

Now, in integrating equation (6.8) over the domain, the left hand side vanishes, again by virtue of the divergence theorem and the zero flux boundary conditions, leaving(6.9)

Finally, the stability condition in equation (6.7) becomes(6.10)by the non-negativity of the dissipation inequality. Armero & Simo (1993) point out that a stability estimate like equation (6.10) does not preclude the formation of shear bands, for example, since in the nonlinear thermomechanical problem, it is not formally a norm on the displacements; it similarly does not preclude the possibility of local snap-back as found by Yazdani (1993) for an isothermal case. The estimate is however a necessary condition for a stable solution, and should be preserved.

## 7. Concluding remarks

The objective of this article has been to develop and present the underlying formulation of a fully coupled anisotropic thermomechanical model for the analysis of concrete under all regimes of stress and transient temperature. Moreover, the explicit intention has been to study the role and value of thermodynamics both in the detailed formulation and underlying justification for the structure and evolution of the model. The aim has not been to identify the values of all parameters or to undertake large parametric studies. Further study would refine the model characteristics.

To this end, we have presented the background details of a thermodynamic framework, studied interactions of model parameters, developed an anisotropic damage surface and evolution laws, and shown that they all satisfy the fundamental laws of thermodynamics. Formal results are given for the evolution of entropy, the satisfaction of the dissipation inequality, and the stability of the boundary value problem in a computational sense. Thermodynamic requirements for the interpretation of definition of specific heat, which may in general be variable, are also given.

## Footnotes

One contribution of 7 to a Theme Issue ‘Thermodynamics in solids mechanics’.

- © 2005 The Royal Society