An anisotropic linear thermo-viscoelastic constitutive law: Elastic relaxation and thermal expansion creep in the time domain.

A constitutive material law for linear thermo-viscoelasticity in the time domain is presented. The time-dependent relaxation formulation is given for full anisotropy, i.e., both the elastic and the viscous properties are anisotropic. Thereby, each element of the relaxation tensor is described by its own and independent Prony series expansion. Exceeding common viscoelasticity, time-dependent thermal expansion relaxation/creep is treated as inherent material behavior. The pertinent equations are derived and an incremental, implicit time integration scheme is presented. The developments are implemented into an implicit FEM software for orthotropic material symmetry under plane stress assumption. Even if this is a reduced problem, all essential features are present and allow for the entire verification and validation of the approach. Various simulations on isotropic and orthotropic problems are carried out to demonstrate the material behavior under investigation.

Introduction

The apparent effects of viscoelastic material behavior manifest themselves in a time-dependent response to loading and accompanying energy dissipation. Relaxation or creep occurs when a viscoelastic material is exposed to quasi-static loads and load changes. Under cyclic excitation damping is exhibited. These phenomena are treated typically in the time and in the frequency domain, respectively.

Such viscoelastic effects are widespread in natural as well as in engineering materials. Among them are almost all biological tissues and most polymers, in particular thermoplastic materials. Biological materials are composites “by nature”, whereas engineering polymers are often mixed with other constituents to improve their performance. These materials are likely to exhibit thermal expansion relaxation/creep behavior, i.e. time-dependent changes of unconstrained thermal strains as a response to temperature changes. Such composites, natural and man made ones, often show direction dependent properties and the consideration of anisotropy becomes inevitable.

A general introduction to viscoelasticity can be found, e.g. in Lakes (2009), Gurtin and Sternberg (1962), Schapery (2000), which focus predominantly on isotropic behavior and treat the time as well as the frequency domain. Transversely isotropic (and orthotropic) viscoelastic models have been presented, e.g. in Holzapfel (2000), Bergström and Boyce (2001), Pandolfi and Manganiello (2006), Puso and Weiss (1998), Taylor et al. (2009), Nedjar (2007), based on invariants of the strain representation, which are commonly used for biological soft tissues. Typically, they combine non-linear orthotropic elasticity with linear isotropic viscosity. Linear viscoelasticity of orthotropic media is presented in Carcione (1995), Dong and McMechan (1995) in the context of geo-materials. Linear and non-linear orthotropic viscoelasticity adopting not the full set of orthotropic viscous effects is presented in Melo and Radford (2004), Kaliske (2000) and applied to composites and foam materials, respectively. Full elastic anisotropy is treated in Lévesque et al. (2008), but with one single relaxation function operating on the entire tensor. In Schapery (1964) full anisotropy is modeled with various relaxation functions applied to elasticity tensor contributions. A review on damping of composites is given in Chandra et al. (1999). Anisotropic non-linear viscoelastic constitutive laws have been introduced in Schapery (1969), Zocher et al. (1997), Sawant and Muliana (2008), Poon and Ahmad (1998) which found widespread application to study composite materials. Implementation aspects are discussed in Crochon et al. (2010), Sorvari and Hämäläinen (2010).

Relaxation type effects occurring for the thermal expansion behavior are mentioned in Lakes (2009), Zocher et al. (1997), Schapery (1964), and observed experimentally for polymers, e.g. Spencer and Boyer (1946). A hygro-thermal expansion relaxation function is presented in Bažant (1970) in the context of concrete behavior. Otherwise, to the knowledge of the authors, there are no modeling attempts in the literature for thermal expansion creep. It is noted that thermal expansion relaxation/creep, here, does not refer to secondary effects arising, e.g. from transient thermal fields, but solely denotes the time-dependent response as inherent material behavior.

In the present work linear thermo-viscoelasticity of anisotropic materials is treated in the time domain. Thereby the thermo-elastic as well as the viscous responses are considered to behave anisotropically. The temperature dependence for thermorheologically simple materials is accounted for by a time–temperature shift function. Time-dependent, anisotropic thermal expansion relaxation/creep is treated.

The constitutive formulations are developed and the pertinent equations are given, including the consistent material tangent operator. The state of the material is handled by internal variables. The algorithms are implemented into an implicit Finite Element Method (FEM) software package. This is done under the assumption of plane stress states and orthotropic materials, but the approach is general and the implementation can be extended to tri-axial problems and general anisotropy in a straightforward manner just by adding additional terms of the same form.

For examples of applications, homogenized orthotropic material data is computed from a simple unit cell approach which serves as input to the developed constitutive laws. This allows for a consistent verification and validation of the developments. Not only the principal load cases are simulated to extract the input data, but also general loading scenarios and loading histories.

Linear viscoelastic constitutive model

The constitutive model for considerations in time domain has to give the stress response as a function of the strain (rate) history and the appropriate instantaneous material behavior. For isotropic material symmetry the stress and strain tensors can be split into their volumetric and deviatoric part. The corresponding scalar valued material properties are the bulk and shear relaxation functions. A similar decomposition can be performed for cubic material symmetry (e.g. used in Pettermann and Hüsing 2012) and for the transversely isotropic case; see e.g. Suquet and Bornert (2001). Alternatively, one can deal with the tensorial equations, as is pursed in the present work which allows for the treatment of the most general anisotropic case. Both the elastic and the viscous contributions to the response are allowed to be of general anisotropy.

Adopting Voigt notation, the hereditary integral in tensorial form reads Stress and strain are given by with . Note that in Voigt notation, coordinate transformation operators must be formulated correspondingly. The time-dependent material tensor for anisotropic materials is given by the symmetric matrix, . In the most general case it is composed of 21 individual relaxation functions of the form, which are given by Prony series expansions with being the instantaneous values, i.e., the elasticity tensor elements giving the short term behavior. The sums over Prony terms contain the relative relaxation elements, , and the corresponding characteristic times, . Note that the relaxation tensor, , possesses off-diagonal terms with , implicitly containing some Poisson relaxation which is not necessarily monotonic; see e.g. Lakes and Wineman (2006). The stress components can be expressed by combining Eqs. (1) and (3) as, where the sum over gives the contributions of the individual strain components. The “relative creep contribution” for each Prony term, is introduced and summed to the “total creep contributions” by where the superscript indicates the cause (i.e., the applied strain) and the subscript denotes the stress component in Eq. (4) which is affected. This unusual super/subscript notation is chosen to discern different types of “strain like” quantities unambiguously. Finally, the stress components read,

Note that the expressions and cannot be interpreted as creep strains, and . Consequently, they are rather property type quantities carrying information on how much of the relaxation potential has been consumed until time . Nevertheless, the for every Prony term represent the internal state variables.

The proposed approach may be viewed as an ad hoc extension of common linear viscoelastic models. Its thermodynamic consistency, however, is not ensured and could be violated by a unsuited set of material parameters.

The temperature dependence of viscoelastic thermorheologically simple materials is widely treated by introducing a reduced time; see e.g. Lakes (2009), with the shift function, , depending on the temperature, . In the present work the empirical Williams–Landel–Ferry (WLF) function, see e.g. Lakes (2009), is used with being the reference temperature and and are constants pertaining to the particular material under consideration. In general cases, each relaxation element, , can have its own time–temperature shift function, .

Thermal expansion relaxation/creep constitutive model

In linear thermo-(visco)elasticity assuming small strains, additive decomposition can be applied, which gives the total strain (in Voigt notation), as the sum of the viscoelastic ve and the thermal th contributions.

Now, materials may show thermal expansion which exhibits time dependence as an inherent material behavior (e.g. composites). Following Lakes (2009), Zocher et al. (1997), the time-dependent thermal expansion relaxation is given by with being the time-dependent coefficients of thermal expansion and, here, temperature independent behavior is assumed. is the temperature change from some (stress free) starting temperature. For the time independent case, the usual thermal expansion results from integration of Eq. (11) as . Note that Eq. (11) exactly resembles the hereditary integral in Eq. (1) for the time-dependent stresses under strain loading. Thus, the further treatment goes in perfect analogy to linear viscoelasticity as studied before.

For the most general anisotropic case, the time-dependent thermal expansion tensor (in correspondence to the strain in Eq. (2)) reads The components are given by the six individual relaxation functions, with the instantaneous thermal expansion coefficients and the relative relaxation values as well as characteristic times for the Prony series with terms. For common materials, the relative relaxation values, , are negative which give rise to creep type thermal expansion behavior.

Applying Eqs. (4) to (7) analogously, one obtains the thermal strain: The “relative creep temperatures”, , have no direct physical meaning, however, they carry the information on the amount of relaxation and, consequently, are used as internal state variables. The concept of reduced time and a shift function, Eq. (8), can be applied equivalently.

The viscoelasticity and the thermal expansion relaxation have to be solve simultaneously, fulfilling Eq. (10).

Time integration and FEM implementation

The scheme of the implicit incremental time solution for arbitrary strain histories in the context of FEM can be found e.g. in Crochon et al. (2010), Sorvari and Hämäläinen (2010), Poon and Ahmad (1998), Muliana and Khan (2008). The pertinent equations are given briefly in the following. Key to the stress update and the consistent Jacobian formulation is the time integration of Eq. (1) with the material properties given in Eq. (3). Besides the strain rate, the only time-dependent contributions are the exponential functions in Eq. (3). They appear in all elements of , just with different characteristic times. With focus on these terms one can extract an integral of the form, where is the appropriate strain rate. Here, is the “relative creep contribution” as introduced in Eq. (5). The above equation is solved for a finite increment of time, , added to time at which the state is completely known and under the assumption that the strain increases linearly with time. This yields the increment of the relative creep contribution as, with the superscript n denoting the state at the beginning of the increment. With this expression the stress update can be given as The role of the solution dependent state variables is taken by , accumulated over the time increments. The equation for computing the consistent tangent moduli elements depends on the time increment only, and does not include any matrix summation. In the implementation of this solution scheme, “symmetry” on the creep contributions (and ) must not be enforced; see Eq. (7). However, the material input requires , , and , to maintain symmetry in the relaxation tensor at any time.

The treatment of thermal expansion relaxation/creep is equivalent to the stress update, because of the identical structure of the governing equations, and it yields the increment of thermal strain components, Similar to the consistent tangent moduli, the coefficients of thermal expansion at the end of the time increment are required for effective convergence. They are computed in perfect analogy with Eq. (18).

Orthotropic plane stress implementation

The implementation of the viscoelastic model is carried out with the UMAT option, the thermal expansion with the UEXPAN option, for the implicit FEM software package ABAQUS/Standard v6.14 (Dassault Systèmes Simulia Corp., Providence, RI, USA). The present implementation allows for 13 Prony terms maximum for any of the relaxation type properties. It is carried out for the case of orthotropic material symmetry under plane stress conditions. Even if this poses a reduced problem, all essential features of the proposed constitutive law are included, and the model formulation and the software realization can be tested and verified entirely. Extensions to more general cases just need the consideration of additional terms in the summation loops of exactly the same type.

For constant temperatures which are not equal to the WLF-reference temperature, Eqs. (8) and (9) can be applied directly to compute the reduced time and likewise the reduced time increment. If the temperature changes within a time increment, an approximation for small temperature increments is utilized, assuming one unique time–temperature shift function for all elements. Extension to individual shift functions is straightforward by modifying Eqs. (16) and (18). The time–temperature shift option is implemented for the viscoelastic behavior only, but not (yet) for the time-dependent thermal expansion relaxation/creep.

To control the accuracy of the solution, i.e. controlling in Eq. (16), an automatic time stepping procedure is implemented.

Material parameters

For the measurement of relaxation data which can serve as input to a viscoelastic constitutive law there are standard procedures for simple loading scenarios. For intricate material response as given by , but in particular by , a direct measurement can be very difficult or impossible. So, to calibrate the input from a series of (feasible) experiments some kind of optimization procedure might be useful.

In the present work the input data is generated by micromechanical type unit cell predictions which directly yields the viscoelastic relaxation matrix elements and thermal expansion relaxation coefficients. These functions of time are fitted by Prony series expansions.

Verification and examples

For clarity of the following, the stress and strain components will be given with their full indication. For elements of the plane stress relaxation tensors (elastic and thermal expansion), the compact Voigt notation is maintained. Unless stated otherwise, all simulations are carried out at the WLF-reference temperature. All simulations employ four noded plane stress elements, using reduced integration for single element tests and full integration for the structural simulations. The element thickness is 1 mm.

Isotropic material—single element tests

Single element tests with isotropic linear viscoelastic material are carried out first. The element size is . The ABAQUS built in material is used to obtain reference solutions, the generic material data (of some glassy polymer) is given in Table 1, which show shear relaxation only, plotted in Fig. 1 (left) as “R33”. Simulations are performed to compute the material parameters pertaining to the developed material model, i.e. the response to “uni-axial strain” loading is sought for. Here, uni-axial strain denotes the in-plane state, in the perpendicular direction plane stress is maintained throughout this work. A Heaviside step function is applied at and kept constant afterwards as . The other in-plane strain components are and at any time. From this load case the components of the relaxation matrix can be deduced directly, see Fig. 1 (thin solid lines). Note that the response which contains the Poisson effect, , is non-monotonic. For completeness, the shear behavior is added which exactly resembles the input, .

Table 1

Isotropic linear thermo-viscoelastic plane stress material data as input to the ABAQUS material law; instantaneous elastic moduli, bulk and shear relaxation by one Prony term each, WLF temperature shift data, and coefficient of thermal expansion

E_0=2500MPaν_0=0.25(G_0=1000MPa)

b=0(τ^b=1.0s)g=0.5τ^g=1.0s

T_ref=0^∘CC_1=17C_2=50

α_0=1×10^−4/^∘C

Fig. 1

Single element predictions of the relaxation functions as response to Heaviside step strains by the ABAQUS built in material law (“abamat”) and the developed constitutive material law (“umat”) calibrated to the former one

Based on this data, the material parameters for the UMAT are calibrated. The (sum of the) relative relaxation elements can be computed directly from the long term and short term response. The characteristic times can be deduced from the initial slope. If two Prony terms are required, a manual “trail and error” procedure is used.

The diagonal element is well captured by one Prony term, for the off-diagonal element, , at least two Prony terms are required to fit the non-monotonic behavior. The resulting material parameters for the UMAT are listed in Table 2. Now, these parameters and the developed UMAT are employed to simulate the same load case as before; see Fig. 1 (bold dashed lines). Comparison to the input shows some deviations but the resulting material behavior is considered sufficiently accurate. Note the difference in relaxation for various term.

Table 2

Isotropic linear viscoelastic plane stress material data as input to the UMAT; instantaneous elasticity matrix elements, relaxation matrix elements by two Prony terms each, and WLF temperature shift data

R_110=2666MPaR_220=2666MPaR_120=R_210=666MPaR_330=1000MPa

k=1r_11=0.41r_22=0.41r_12=r_21=0.25r_33=0.25τ^r_11=1.25 sτ^r_22=1.25 sτ^r_12=τ^r_21=2.0 sτ^r_33=1.0 s

k=2r_11=0r_22=0r_12=r_21=−0.11r_33=0.25τ^r_11=1.25 sτ^r_22=1.25 sτ^r_12=τ^r_21=0.6 sτ^r_33=1.0 s

T_ref=0^∘CC_1=17C_2=50

A series of additional tests and verifications are run based on single elements. Among them are classical relaxation test (under uni-axial stress), bi-axial loading, and arbitrary load cases. Besides the response to step functions, arbitrary ramp loading and load reversals are realized. All comparisons give satisfactory results. Also, the time–temperature shift function for temperature dependent viscoelastic material data is tested successfully as well as the automatic time stepping algorithm.

Isotropic material—structural simulations

Next, simulations of structural problems with isotropic material are carried out which allow direct comparison with reference solutions. For this purpose a quadratic patch of unit length is discretized with plane stress elements. The patch contains two parallel rectangular voids with ; see Fig. 2. The boundary conditions are “unit cell” like to be used later for homogenization; see Pettermann and Suresh (2000), Böhm (2004). The deformations of the nodes at the edges are coupled to the adjacent corner nodes of the model to remain on a straight line. The deformation of the model is controlled by the corner nodes, whereas the upper right corner is coupled to the others to maintain a parallelepipedic shape. This way, a unit cell type model is set up with three master nodes to apply loads and to read the response. Note that this does not give correct periodic boundary conditions (see e.g. Pettermann and Suresh 2000; Böhm 2004) under shear loads. However, it gives some homogenized behavior which sufficiently serves for the purpose of comparison in the present study.

Fig. 2

Structure with voids (light gray), also representing a unit cell of an orthotropic material

The master nodes are used to apply a variety of loading scenarios as discussed above. The material model employed is either the ABAQUS built in one or the developed UMAT, whereas the voids are realized by setting the Young’s modulus to . The predicted master node responses and the local, time-dependent stress fields are compared, and for all studied load cases almost one to one agreement is found.

The implemented automatic time stepping as well as the consistent tangent matrix, i.e., the Jacobian, work correctly and the performance of the implemented material law is very satisfactory.

Orthotropic porous material

Now, the structure investigated in the previous section can be interpreted as porous material with orthotropic properties. By homogenization, the effective linear viscoelastic orthotropic plane stress behavior can be predicted. Applying the load cases of Heaviside step uni-axial strain in the two principal directions and for shear, the entire constitutive behavior can be predicted (as introduced in Sect. 5.1). The relaxation response is presented in Fig. 3 (denoted by “unit cell”), with . All elements of are approximated by one Prony term, the material parameters are listed in Table 3. These material parameters together with the developed UMAT are now employed to run single element simulations. Comparison to the unit cell response is shown in Fig. 3 and excellent agreement can be seen for all elements. In addition to the direct comparison of the fitted response, a series of general load cases are investigated and compared successfully. The assumption of taking the same time–temperature shift as for the base material is confirmed to hold.

Fig. 3

Orthotropic relaxation functions as response to Heaviside step strains; homogenization of a voided structure (“unit cell”) and single element predictions by the developed constitutive material law (“umat”) calibrated to the homogenized behavior

Table 3

Orthotropic linear viscoelastic plane stress material data as input to the UMAT; instantaneous elasticity matrix elements, relaxation matrix elements by one Prony term each, and WLF temperature shift data (taken to be equal to that of the input material)

R_110=1693MPaR_220=1124MPaR_120=R_210=227MPaR_330=242MPa

k=1r_11=0.440r_22=0.439r_12=r_21=0.234r_33=0.453τ^r_11=1.09 sτ^r_22=1.09 sτ^r_12=τ^r_21=1.98 sτ^r_33=1.08 s

T_ref=0^∘CC_1=17C_2=50

Note the orthotropic behavior of the relaxation, too, which is fully captured by the present approach.

Thermo-viscoelastic orthotropic composite

In this section an example is presented to study the thermal expansion relaxation. To this end the geometrical model of the structure from the previous examples is used, but the regions of the voids are now filled with rigid material to obtain a model composite. The reinforcement material is taken to be quasi-rigid with a Young’s modulus of , the Poisson ratio is set to 0.25. The isotropic coefficient of thermal expansion is set to . For the matrix material the data from Table 1 is used, disregarding temperature dependence (i.e., no WLF shift is applied).

As in the previous examples, the unit cell is employed to perform material characterization, i.e., to compute the effective linear thermo-viscoelastic properties of the composite. The necessary (mechanical) load cases are simulated and the input data to the UMAT is fitted. Again, the off-diagonal term, , exhibits some curvature change around and is fitted by two Prony terms. The calibrated homogenized material data is listed in Table 4, graphical presentation is omitted.

Table 4

Orthotropic linear thermo-viscoelastic plane stress material data of a rigid inclusion reinforced model composite as input to the UMAT and the UEXPAN; instantaneous elasticity matrix elements, relaxation matrix elements, instantaneous coefficients of thermal expansion, and thermal expansion relaxation data

R_110=5301MPaR_220=3857MPaR_120=R_210=817MPaR_330=714MPa

k=1r_11=0.419r_22=0.413r_12=r_21=0.203r_33=0.434τ^r_11=1.20 sτ^r_22=1.20 sτ^r_12=τ^r_21=2.00 sτ^r_33=1.00 s

k=2r_11=0.000r_22=0.000r_12=r_21=−0.034r_33=0.000τ^r_11=1.20 sτ^r_22=1.20 sτ^r_12=τ^r_21=0.35 sτ^r_33=1.00 s

α_10=5.16×10^−5/^∘Cα_20=7.58×10^−5/^∘C

k=1a_1=−0.019a_2=−0.045τ^a_1=1.00 sτ^a_2=1.30 s

In addition, the effective thermal expansion behavior of the composite can be predicted by the unit cell simulations. A Heaviside temperature step, , from an undeformed, stress free configuration is applied and the thermal expansion response is evaluated. The constituents of the model composite posses different coefficients of thermal expansion, thus, a uniform temperature change of the composite induces stresses and strains in the constituents (the stresses, of course, are self-equilibrated). Since these stresses and strains show relaxation and creep, the effective thermal expansion shows relaxation/creep, too. Note that no thermal diffusion is involved here. Due to the orthotropic material symmetry there are two independent thermal expansion relaxation/creep functions. The corresponding unit cell predictions are presented in Fig. 4 (thin lines).

Fig. 4

Orthotropic thermal expansion relaxation/creep functions as response to unit temperature Heaviside step; homogenization of a rigid inclusion reinforced model composite (“unit cell”), and single element predictions by the developed thermal expansion law (“uexpan”) calibrated to the homogenized behavior

To obtain input data to the UEXPAN subroutine, the thermal expansion response is fitted, and the material data is listed in Table 4. It comprises the two instantaneous coefficients of thermal expansion and the corresponding relaxation/creep functions with one Prony term each. The fitted behavior is presented in Fig. 4 (thick dashed lines) for comparison with the unit cell predictions.

So far, the material data input is calibrated and other thermo–viscoelastic load cases are studied. Temperature changes under fully and uni-axial in-plane constraint are simulated which give rise to thermally induced stresses. For both cases excellent agreement is found for the unit cell response and the UMAT/UEXPAN prediction. For comparison, the same load cases are simulated with coefficients of thermal expansion taken to be time independent, i.e. suppressing thermal expansion relaxation. As expected, the stress responses deviate for increasing time.

Summary

A constitutive material law for linear thermo-viscoelasticity in the time domain with full anisotropy is presented. For the linear viscoelastic material behavior the anisotropic material symmetry is not only considered for the elastic part. Also for the relaxation response, full anisotropy is realized with the appropriate number of independent material values. The formulation is based on a time-dependent elasticity tensor for which each element possesses its own relaxation function. This way, the mutual coupling of the normal components, i.e., Poisson type effects, is accounted for. The relaxation function for each tensor element is prescribed by its individual Prony series expansions.

An incremental algorithm is set up which gives the time-dependent stress response to a prescribed incremental strain history. State dependent internal variables are utilized to account for the loading history, and the consistent material Jacobian is formulated. Temperature dependence is realized by a time–temperature shift function of the Williams–Landel–Ferry type. It is implemented for orthotropic material symmetry under plane stress assumption into a FEM package. The extension to tri-axial stress and strain states and to general anisotropy is straightforward.

The thermal expansion relaxation/creep is considered as inherent material behavior. It is treated in analogy to the mechanical, i.e., viscoelastic, model. The material behavior is described by instantaneous coefficients of thermal expansion and their individual relaxation functions. The time-dependent thermal strain is modeled as a function of the applied temperature change. Again, internal state variables are used to account for the temperature history.

Various tests on isotropic and orthotropic plane stress problems are carried out successfully for verification and validation. The performance is found to be very satisfactory. Orthotropic test cases are set up by considering heterogeneous model composites and employing a periodic unit cell type approach. Thereof, homogenized material data is computed, serving as input to the developed material law. Also the thermal expansion relaxation behavior is investigated.

Finally, the developed constitutive laws provide tools for structural analyses of components made from orthotropic linear thermo–viscoelastic materials.