Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 50, 1, pp. 3-22, Warsaw 2012 50th anniversary of JTAM THEORETICAL AND NUMERICAL ASPECTS IN WEAK-COMPRESSIBLE FINITE STRAIN THERMO-ELASTICITY Ahmad-Wahadj Hamkar Stefan Hartmann Clausthal University of Technology, Institute of Applied Mechanics, Clausthal-Zellerfeld, Germany e-mail: stefan.hartmann@tu-clausthal.de In this essay, a constitutive model for nearly incompressible elastic be- havior is extended to the case to thermal effects. First, the use is made of the multiplicative decomposition of the deformation gradient into a thermal and a mechanical part. The thermal part is purely volumetric. Additionally, the mechanical part is multiplicatively decomposed into a volume-preserving and a volume-changing part so that the final stress state shows the influences of the temperature-dependence. The propo- sedmodel is carefully studied in view of the thermo-mechanical coupling effects. Second, the model is implemented into a time-adaptive finite element formulation based on higher-order Rosenbrock-type methods, which is a completely iteration-free procedure so that really fast com- putations are available. The article concludes with a three-dimensional numerical simulation of a representative elastomeric tensile specimen. Key words: thermoelasticity, finite strains, finite elements, Rosenbrock- -type methods 1. Introduction The incorporation of thermo-mechanical coupling effects starts with the de- composition of the deformation gradient F(~X,t) = Grad~χR(~X,t), into a thermal and a mechanical part F=FMFΘ, where ~x = ~χR(~X,t) defines mo- tion of a particle ~X occupying the place ~x at time t. This proposal of Lu and Pister (1975) is taken up in a number of papers (see Miehe, 1988; Holzapfel and Simo, 1996b; Lion, 1997; Heimes, 2005). An additional decomposition of the mechanical part into a volume-changing and volume-preserving part goes 4 A.-W. Hamkar, S. Hartmann back to Flory (1961). There, a particular assumption of the free-energy yields a clear assignment of kinematical quantities to the hydrostatic and deviatoric Cauchy stress state. Since the thermal deformation is purely volumetrical, this approach is of particular interest for thermo-mechanical deformation proces- ses as well. References with similar approaches are given by Simo and Miehe (1992), Miehe (1995), Holzapfel and Simo (1996a). In the case of isothermal problems, it was found out that for Rivlin- Saunders typemodels, see for a detailed investigation inHartmann (2001a,b), it is not possible to show the existence of a solution if the models are exten- ded to the case of nearly incompressibility (see Hartmann and Neff, 2003). Accordingly, Hartmann andNeff (2003) proposed a new class of strain-energy functions,whicharebasedon invariants of theunimodularrightor leftCauchy- Green tensors.Oneadditional aspect lies in the fact that nearly all hyperelasti- city relations, which are based on the decomposition into volume-preserving and volume-changing parts, show non-physical behavior in the tensile and compression test, see Ehlers andEipper (1998) aswell. However, inHartmann and Neff (2003) it is shown that this essentially depends on the choice of the parts in the free-energy. Since thermo-mechanical coupled problems essentially depend on the volume-change, at least resulting from the volumetric tempe- rature expansion, the investigations have to be extended to thermoelasticity, which is done in this paper. In contrast to Miehe (1995) and Holzapfel and Simo (1996a) other func- tions of the free-energy function are applied, and the derivation of the heat equation is based here on the proposed decomposition of the deformation gradient. The latter leads to a consistent representation leading to an exact expression of the thermo-elastic coupling term and the specific heat capacity. This strict derivative has only be carried out by Lion (1997) for the case of incompressibility and without any scopes for the numerical treatment, which results in a different formulation. A second important aspect treats the implementation of the constitutive model into a finite element program.Commercial programs have the difficulty that the heat equation does not offer the possibility to incorporate an arbi- trary process-dependent heat source so that a real thermo-mechanical coupled problem is not able to be computed.Moreover, staggered schemes, which have their advantages in the application of two independent codes (a pure thermal and a pure mechanical one), cannot guarantee to follow the exact solution. Moreover, these are time-integration methods of the order one. Accordingly, an efficient fully monolithic finite element approach has to be developed. In this context, we pick up a very efficient finite element concept developed in Theoretical and numerical aspects in weak-compressible... 5 Hartmann and Wensch (2007), Hartmann and Hamkar (2010) for the case of isothermal problems. Due to the fact that both the discretized heat equation in the finite element framework represents a system of ordinary differential equations and the discretized weak formulation of the equilibrium conditions (or principle of virtual displacements),whichdefinea systemofnon-linear equ- ations, are coupled, Rosenbrock-typemethods are applied to the resulting sys- temof differential-algebraic equations (DAE-system). This is done in thiswork for the first time. Because of the superior behavior of Rosenbrock-type me- thods on the case of smooth problems,which lead after the time-discretization of the DAE-system only to the solution of a linear system of equations within each point in time (see Hartmann andWensch, 2007; Hartmann andHamkar, 2010; Lang, 2000) in the context of finite elements, very efficient and consi- stent computations are obtained. Moreover, the use is made of higher-order time-integration methods and a time-adaptive procedure (step-size control). The paper is organized as follows: first, the thermo-mechanically coupled constitutive model is developed. Afterwards, the resulting DAE-system is de- rived showing the coupling aspects. Finally, a numerical example using the fully iteration-free method for a problemwith different time-scales is conside- red.Thenotation inuse is defined in the followingmanner: geometrical vectors are symbolizedby ~a, second order tensors Abybold-facedRoman letters, and calligraphic letters A define fourth order tensors. Furthermore, we introduce matrices at the global level symbolized by bold-faced italic letters A. 2. Constitutive assumptions In order to develop the constitutive model, the multiplicative decomposition of the deformation gradient F into a thermal FΘ and amechanical part FM according to the proposal of Lu and Pister (1975) is assumed F=FΘFM (2.1) Additionally, the mechanical part FM = F̂MFM (2.2) is decomposed into a volume-preserving part FM and a volume-changing part F̂M according to Flory (1961) F̂M =(detFM) 1/3I detF̂M =detFM FM =(detFM) −1/3FM detFM =1 (2.3) 6 A.-W. Hamkar, S. Hartmann The thermal part is assumed to be purely volumetric FΘ = ϕ 1/3I ϕ = ϕ̂(Θ−Θ0)= ϕ̂(ϑ) (2.4) where ϕ̂(0)= 1 shouldhold and the temperaturedifference ϑ = Θ−Θ0 is defi- ned. Θ is the absolute temperature and Θ0 defines the reference temperature. The determinant of the thermal deformation is obviously given by detFΘ = ϕ̂(Θ−Θ0)= ϕ̂(ϑ) (2.5) describing the volumetric deformation caused by the temperature change ϑ, which is chosen to be linear in view of observations in elastomers ϕ̂(Θ−Θ0) := 1+α(Θ−Θ0)=1+αϑ (2.6) see Treloar (1975). In view of the total deformation detF=det(FΘFM)= (detFΘ)(detFM)= ϕ̂(Θ−Θ0)(detFM) (2.7) holds. Accordingly, we have F= F̂F with { F̂ =(ϕdetFM) 1/3I F =FM (2.8) with F = (detF)−1/3F. Sometimes it is useful to introduce the abbrevia- tion for the determinants J := detF = JΘJM, JΘ := detFΘ = ϕ and JM := detFM = J/ϕ which are used later for short notational purposes. Using the imagination of a fictitious thermal unloading, similar to the case of themultiplicative decomposition of the deformation gradient into an elastic and a plastic state (see Haupt, 1985), the mechanical Green strain tensor EM := lim ϑ→0 E= 1 2 (F⊤MFM − I) (2.9) is defined, where E= 1 2 (F⊤F− I) (2.10) is theGreen strain tensor itself. Thismotivates the thermal part of theGreen- Lagrange-type EΘ =E−EM = 1 2 (F⊤F−F⊤MFM)= 1 2 (C−CM) (2.11) Theoretical and numerical aspects in weak-compressible... 7 or vice versa the additive decomposition E=EM +EΘ (2.12) The push-forward operation F−⊤M EF −1 M yields the decomposition Γ̂= Γ̂M + Γ̂Θ (2.13) with Γ̂=F−⊤M EF −1 M = 1 2 (F⊤ΘFΘ−F−⊤M F −1 M )= 1 2 (ϕ2/3I−B−1M )= ϕ2/3 3 (I−B−1) Γ̂M =F −⊤ M EMF −1 M = 1 2 (I−F−⊤M F −1 M )= 1 2 (I−B−1M ) (2.14) Γ̂Θ =F −⊤ M EΘF −1 M = 1 2 (F⊤ΘFΘ− I)= 1 2 (CΘ− I)= 1 2 (ϕ2/3−1)I where Γ̂, Γ̂M and Γ̂Θ measure the strains relative to the mechanical interme- diate configuration BM. Here, the right and left Cauchy-Green tensors C=F⊤F CΘ =F ⊤ ΘFΘ = ϕ 2/3I CM =F ⊤ MFM = ϕ −2/3C B=FF⊤ BΘ =FΘF ⊤ Θ = ϕ 2/3I BM =FMF ⊤ M = ϕ −2/3B (2.15) are introduced. The definition of strainmeasures (2.14)2,3 have the advantage that they are purely mechanical and purely thermal. For later purposes, the mechanical deformation is decomposed into volume-changing and preserving parts, defined by JM = detFM = (detCM) 1/2 and the unimodular, mecha- nical right Cauchy-Green tensor CM = (detCM) −1/3CM, detCM = 1, is introduced. Additionally, strain-rate tensors relative to the mechanical intermediate configuration BM can bedefinedon the basis of thematerial time derivative of (2.9)-(2.11) by the corresponding push-forward operation F−⊤M (. . .)F −1 M . This yields the strain-rate measures relative to the intermediate configuration BM △ Γ̂=F−⊤M ĖF −1 M = ˙̂ Γ+L⊤MΓ̂+ Γ̂LM △ Γ̂M =F −⊤ M ĖMF −1 M = ˙̂ ΓM +L ⊤ MΓ̂M + Γ̂MLM = 1 2 (LM +L ⊤ M)=:DM (2.16) △ Γ̂Θ =F −⊤ M ĖΘF −1 M = ˙̂ ΓΘ+L ⊤ MΓ̂Θ+ Γ̂ΘLM 8 A.-W. Hamkar, S. Hartmann and implies the additive decomposition △ Γ̂= △ Γ̂M + △ Γ̂Θ (2.17) The strain-rate △ Γ̂M is purely mechanical, whereas the thermal strain-rate re- lative to the intermediate state can be calculated by △ Γ̂Θ = 1 3 ϕ̂′(Θ−Θ0)Θ̇(t)ϕ−1/3I ︸ ︷︷ ︸ ˙̂ ΓΘ +(ϕ2/3−1) △ Γ̂M (2.18) see Eqs. (2.16)3, (2.14)3, and (2.16)2. Next, these strain measures are used within the specific stress power to develop appropriate stressmeasures according to the concept of dual variables (Haupt and Tsakmakis, 1989, 1996) p = 1 ρ T ·D= 1 ρR T̃ · Ė= 1 ρR T̃ · (F⊤M △ Γ̂FM)= 1 ρR (FMT̃F ⊤ M) · △ Γ̂= 1 ρR ŜM · △ Γ̂ (2.19) Here, Eq. (2.16)1 is exploited and aKirchhoff-type stress tensor relative to the mechanical intermediate configuration ŜM =FMT̃F ⊤ M (2.20) is introduced. T̃=(detF)F−1TF−⊤ defines the secondPiola-Kirchhoff stress tensor, T the Cauchy stress tensor and D = (L+L⊤)/2 the spatial strain- rate tensor depending on the spatial velocity gradient L = ḞF−1. Inserting strain-rates (2.17) and (2.18) into specific stress power (2.19), leads to p = 1 ρR ŜM · △ Γ̂= 1 ρR ŜM · ( △ Γ̂M + △ Γ̂Θ) = 1 ρR ϕ2/3ŜM · △ Γ̂M + ϕ̂′(Θ−Θ0)ϕ−1/3Θ̇(t) 3ρR (trŜM) (2.21) In view of thermo-mechanical processes, the Clausius-Duhem inequality has to be fulfilled −ψ̇−Θ̇s+ 1 ρR T̃·Ė− ~q ρΘ ·gradΘ =−ψ̇−Θ̇s+ 1 ρR SM· △ Γ̂− ~q ρΘ ·gradΘ ­ 0 (2.22) Theoretical and numerical aspects in weak-compressible... 9 where ψ defines the specific free-energy, s is the entropy, and ~q the heat flux vector. In the following, the proposals of Lion (2000) and Heimes (2005) are taken up implying the assumption that the free-energy depends on the mechanical deformation EM and the temperature Θ ψ(EM,Θ)= ψM(CM,Θ)+ψΘ(Θ) (2.23) Themechanical part is assumed to be linear in the temperature ψM(CM,Θ)= Θ Θ0 ψM(CM) (2.24) with ψM(CM)= U(JM)+υ(CM)= U(JM)+w(ICM ,II CM ) (2.25) which is decomposed into a free-energy U(JM) depending on the purely vo- lumetric, mechanical deformation and a part depending on the isochoric me- chanical deformation. In the later examples, the use is made of U(JM)= K 50 (J5M +J −5 M −2) w(I B ,II B )= α̂(I3 B −27)+ c10(IB−3)+ c01(II 3/2 B −3 √ 3) (2.26) proposed by Hartmann and Neff (2003) as one possibility in the class of po- lyconvex strain-energy functions. The thermal part of strain-energy (2.23) is defined by ψΘ(Θ)= cp ρR [( (Θ−Θ0)−Θ ln Θ Θ0 ) (1−kpΘ0)− 1 2 kp(Θ 2−Θ20) ] (2.27) resulting from the experimental observation of the linear temperature- dependence of the heat capacity (see Heimes, 2005). The evaluation of thematerial time-derivative of the free-energy ψ defined in Eq. (2.23) ψ̇(EM,Θ)= ( 1 Θ0 ψM(CM) ) Θ̇+ Θ Θ0 dψM(CM) dCM · ĊM +ψ′Θ(Θ)Θ̇ (2.28) is required inClausius-Duheminequality (2.22) yielding,bymeansofdefinition (2.20) and the time derivative of (2.9) expressed bymechanical right Cauchy- 10 A.-W. Hamkar, S. Hartmann Green tensor (2.15)3, ĊM =2F ⊤ M △ Γ̂MFM ( 1 ρR ϕ2/3ŜM −2 Θ Θ0 FM dψM(CM) dCM F ⊤ M ) · △ Γ̂M + [ −s− ( 1 Θ0 ψM(CM)+ψ ′ Θ(Θ) ) + 1 3ρR ϕ−1/3ϕ̂′(Θ−Θ0)trŜM ] Θ̇ (2.29) − ~q ρΘ · gradΘ ­ 0 Exploiting arbitrary mechanical strain and temperature paths, implies the three following expressions ŜM = 2ρR ϕ2/3 Θ Θ0 FM dψM(CM) dCM F⊤M s =− 1 Θ0 ψM(CM)−ψ′Θ(Θ)+ 1 3ρR ϕ−1/3ϕ̂′(Θ−Θ0)trŜM ~q =−κgradΘ κ ­ 0 (2.30) Having particular strain-energy function (2.25), yields for the derivatives in elasticity relation (2.30)1 dU((detCM) 1/2) dCM = 1 2 JMU ′(JM)C −1 M (2.31) and dυ(CM(CM)) dCM = [ dCM dCM ] ⊤ dυ dCM (2.32) with [ dCM dCM ] ⊤ =(detCM) −1/3 [ I− 1 3 (C−1M ⊗CM) ] = J −2/3 M [ I− 1 3 (C −1 M ⊗CM) ] (2.33) In this expression, the identity tensor of fourth order I = [I⊗ I]T23 = δikδjl~ei⊗~ej ⊗~ek⊗~el (2.34) is introduced expressed relative toCartesian coordinates, IA=A.Obviously, C −1 M ⊗CM = C −1 M ⊗CM holds. Caused by the particular dependence on Theoretical and numerical aspects in weak-compressible... 11 the invariants of the mechanical unimodular right Cauchy-Green tensors, the application of the chain-rule leads to dυ dCM = ∂w ∂I CM dI CM dCM + ∂w ∂II CM dII CM dCM =(w1+w2ICM )I−w2CM (2.35) with w1(ICM ,II CM )= ∂w ∂I CM and w2(ICM ,II CM )= ∂w ∂II CM (2.36) In other words, we have dψM(CM) dCM = JM 2 U ′(JM)C −1 M +2J −2/3 M [ I− 1 3 C −1 M ⊗CM ] dυ dCM = J 1/3 M 2 U ′(JM)C −1 M + 2 J 2/3 M ( (w1+w2ICM )I −w2CM − 1 3 (w1ICM +2w2IICM )C −1 M ) (2.37) These expressions are used to express the elasticity relation relative to the current and the reference configuration. First, the push-forward operation of the second Piola-Kirchhoff tensor T̃ leads to S=FT̃F⊤ =F(F−1M SMF −⊤ M )F ⊤ =FMFΘ(F −1 M SMF −⊤ M )F ⊤ ΘF ⊤ M = ϕ 2/3 ŜM (2.38) where the use is made of decomposition (2.1) and Eq. (2.20). S = (detF)T defines the weighted Cauchy tensor (Kirchhoff stresses). Comparing (2.38) with (2.30)1, yields S=2ρR Θ Θ0 FM dψM(CM) dCM F⊤M = ρR Θ Θ0 JMU ′(JM)I +2ρR Θ Θ0 J −2/3 M [FM ⊗FM] T23 [ I− 1 3 C −1 M ⊗CM ] dυ dCM = ρR Θ Θ0 JMU ′(JM)I+2ρR Θ Θ0 ( FM dυ dCM F ⊤ M )D (2.39) AD = A − 1 3 (trA)I defines the deviator operator and trA = ak k de- notes the trace operator of the second-order tensor. Here, the properties [A⊗B]T23[C⊗D] = [ACB⊤⊗D] and [C⊗D][A⊗B]T23 = [C⊗A⊤DB] for the expression [FM ⊗FM]T23 [ I− 1 3 C −1 M ⊗CM ] = [ I− 1 3 I⊗ I ] [FM ⊗FM]T23 (2.40) 12 A.-W. Hamkar, S. Hartmann are exploited. [·]T23 symbolizes the transposition of the second and third in- dex of the fourth-order tensors, see for further properties (Hartman, 2002, p. 1464). In conclusion, under the assumption of isotropy and themechanical unimodular left Cauchy-Green tensor BM = J −2/3 M BM =(J/ϕ) −2/3ϕ−2/3B= J−2/3B=B, (2.41) the Cauchy stress tensor T= ρR ϕ Θ Θ0 U ′(J/ϕ)I+ 2ρR J Θ Θ0 (dυ dB B )D (2.42) is obtained. Finally, the entropy (2.30)2 results in s =− 1 Θ0 ψM(CM)−ψ′Θ(Θ)+ 1 3ρR ϕ−1/3ϕ̂′(Θ−Θ0)trS =− 1 Θ0 ( U (J ϕ ) +υ(C) ) −ψ′Θ(Θ)+ Θ Θ0 ϕ̂′(Θ−Θ0) ϕ2 JU ′ (J ϕ ) (2.43) using relation (2.38). 3. Boundary-value problem and discretization In the following the numerical treatment of the underlying equations are di- scussed. The heat equation in the spatial representation is given by ė(~x,t)=− 1 ρ(~x,t) div~q(~x,t)+r(~x,t)+p(~x,t) (3.1) with the specific internal energy e, the density in the current configuration ρ, the heat flux vector ~q, the volume distributed heat source r, which is assumed to be zero in the following, and the stress power p defined in Eq. (2.19) (see Haupt, 2002, p. 124). Frequently, the specific internal energy e is expressed by e = ψ +Θs requiring the time-derivative ė = ψ̇ + Θ̇s+Θṡ in heat equation (3.1). Accordingly, the time-derivatives of free-energy (2.23), see Eq. (2.28) as well, and entropy (2.43) are necessary. This leads to the concrete form of the heat equation cp(C,Θ)Θ̇ = κ ρ div gradΘ−Θ ( β(C,Θ)C−1−Ξ(C,Θ) ) · Ċ (3.2) Theoretical and numerical aspects in weak-compressible... 13 where Fourier’s model for the heat flux ~q = −κgradΘ is assumed. In this case, the consistent derivation of the heat capacity reads cp = JM Θ2 Θ0 ϕ̂′ ϕ [ 2 ( 1 Θ + 1 2 ϕ̂′′ ϕ̂′ − ϕ̂ ′ ϕ ) U ′(JM)− ϕ̂′JMU ′′(JM) ] −Θψ′′Θ(Θ) (3.3) Additionally, we have the abbreviations β(C,Θ)= 1 2 JM Θ Θ0 ϕ̂′ ϕ [( 1− 1 Θ ϕ̂′ ϕ ) U ′(JM)+JMU ′′(JM) ] Ξ(C,Θ)= 1 J2/3Θ0 [ I− 1 3 (C−1⊗C) ]dυ dC (3.4) Heat equation (3.2) is coupled with the local balance of linear momentum (quasi-static formulation) divT+ρ~k =~0 (3.5) and vice versa, where ~k represents the acceleration of gravity. The spatial discretization of coupled partial differential equations (3.2)-(3.5) using finite elements leads to a system of differential-algebraic equations Cp(t, u,Θ)Θ̇= rΘ(t, u, u̇,Θ) g(t, u,Θ)= 0 (3.6) where Cp represents a time-, displacement- and temperature-dependent heat capacitymatrix.Here, u ∈ Rnuu and Θ∈ RnΘu symbolize the unknownnodal displacements and temperatures, respectively. Three modifications are used to solve the problem: first, we are not only interested in the unknown nodal displacements and temperatures u and Θ, butalso in thenodal reaction forces and the heat fluxes at those nodes where the displacements and temperatures are prescribed. According to Hartmann et al. (2008), and the literature cited therein, the method of Lagrange-multipliers is used in order to obtain the reaction forces. The method of Gear (1986), with a similar idea, however, proposed forODEs, offers the possibility to get the heat fluxes at those nodes, where the temperatures are prescribed. In both cases all nodal displacements and temperatures have to be assumed to be unknown, ua ∈ Rnua and Θa ∈ R nΘa. Here, nua = nuu+nup, where nup are the number of prescribed nodal displacements and nΘp = nΘa−nΘu defines the number of prescribed nodal temperatures. Accordingly, two constraints have to be considered Cu(t, ua(t))= û − u(t)= M⊤u ua(t)− u(t)= 0 CΘ(t,Θa(t))= Θ̂−Θ(t)= M⊤ΘΘa(t)−Θ(t)= 0 (3.7) 14 A.-W. Hamkar, S. Hartmann The sizes of the vectors are given by ua = {u, û}⊤ ∈ Rnua, u ∈ Rnuu, û ∈ Rnup, Θa = {Θ,Θ̂}⊤ ∈ RnΘa, Θ ∈ RnΘu, and Θ̂ ∈ RnΘp. û and Θ̂ are those nodal quantities, where the displacements and the temperatures are prescribed by the functions u(t) and Θ(t). Thematrices M⊤u = [0nup×nuuInup] and M⊤ Θ = [0nΘp×nΘuInΘp] filter the required information fromthe total storage information ua and Θa. The second aspect treats the problem that in Eq. (3.4) the heat capacity matrix Cp is a solution-dependentmatrixwhichmight become singular in cer- tain applications. For this problem, the proposal of Lubich and Roche (1990) is followed, where the substitution u̇a = Du and Θ̇a = DΘ is introduced. Finally, it is well-known that the principle of virtual displacements to solve Eq. (3.6)2 is not appropriate because of volumetric locking effects. Thus, a three-field mixed formulation proposed in Simo et al. (1985) or Simo and Taylor (1991) (see Hartmann, 2002; Hartmann and Hamkar, 2010) as well, is applied. To solve the resulting DAE-system F(t, y(t), ẏ(t)) :=    u̇a− Du Θ̇a− DΘ ga(t, ua,Θa)+ Muλu Cp(Θa, ua)DΘ− rΘ(t, ua, Du,Θa)− MΘλΘ Cu(t, ua) CΘ(t,Θa)    = 0 (3.8) theuse ismadeof a totally iteration-free time-integrationmethod, the so-called Rosenbrock-typemethod. It has been successfully studied in the context of iso- thermalproblems inHartmannandWensch (2007) andHartmannandHamkar (2010). The scheme offers the possibility to apply high-order time-integration methods and a step-size control technique (see Hairer and Wanner, 1996). In each point in time (stage) only a system of linear equations has to be solved. The classical method, either fully coupled or staggered, require iterations to solve the entire system. The proposed procedure is embedded in a 3D-finite element program TASA-FEM, where additional aspects for exploiting line- ar relations of DAE-system (3.8) a priorly considered so that the resulting systems of linear equations have merely the size nua+nΘa. 4. Examples Two examples are studied. First, the aspect of physical relevance in the one- dimensional tensile and compression case if applying the standard volumetric Theoretical and numerical aspects in weak-compressible... 15 part U(JM) is discussed. Second, the three-dimensional finite element simu- lation of an elastomeric specimen using a step-size controlled higher-order Rosenbrock-type method to solve the fully thermo-mechanical coupled pro- blem is shown. 4.1. Uniaxial tensile-compression test for constant temperatures In the first investigation, the uniaxial tensile and compression test is assu- medwith constant temperature. In this case, the deformation gradient has the representation F= λ~e1⊗~e1+λQ(~e2⊗~e2+~e3⊗~e3),where λdefines the prescri- bedaxial stretch and λQ theunknownlateral stretch inabar.Thedeterminant of the deformation gradient reads J =detF= λλ2Q and the unimodular left Cauchy-Green tensor of Eq. (2.41) B= B~e1⊗~e1+BQ(~e2⊗~e2+~e3⊗~e3) with the abbreviations B := (λλ2Q) −2/3λ2 and BQ := (λλ 2 Q) −2/3λ2Q The stress state is given by T = σ~e1 ⊗ ~e1, which has to be inserted into thermo-elasticity relation (2.42). This leads with dυ dB B=(w1+w2IB)B−w2B 2 = a~e1⊗~e1+aQ(~e2⊗~e2+~e3⊗~e3) and the abbreviations a =(w1+w2IB)B −w2B 2 aQ =(w1+w2IB)BQ−w2B 2 Q to the deviator (dυ dB B )D = b~e1⊗~e1− b 2 (~e2⊗~e2+~e3⊗~e3) (4.1) with b = 2/3 ( (w1 +w2IB)(B −BQ)−w2(B 2 −B2Q) ) . Here, strain-energies (2.26) with K = 1000MPa, α̂ = 0.00367MPa, c01 = 0.1958MPa and c10 = 0.1788MPa are applied. Additionally, the first and second invariant of the unimodular left Cauchy-Green tensors have to be calculated I B =I C = trB= B+2BQ =(λλ 2 Q) −2/3(λ2+2λ2Q) II B =II C = 1 2 (I B − trB2)= tr(B−1)= (λλ2Q)2/3(λ−2+2λ−2Q ) (4.2) 16 A.-W. Hamkar, S. Hartmann From Eq. (2.42), expressed by the component representation of the Cauchy stress tensor and (4.1), the two equations σ = ρRΘ ϕΘ0 U ′(J/ϕ)+ 4ρRΘ 3λλ2QΘ0 ( (w1+w2IB)(B −BQ)−w2(B 2−B2Q) ) 0= ρRΘ ϕΘ0 U ′(J/ϕ)− 2ρRΘ 3λλ2QΘ0 ( (w1+w2IB)(B −BQ)−w2(B 2−B2Q) ) (4.3) result. In other words, for a given axial stretch λ and temperature Θ, Eq. (4.3)2 has to be iteratively solved to obtain the lateral stretch λQ. In the second step, the entire true stress, (4.4)1, can be evaluated. Obviously, the experimentally observed “linear” dependence of the tem- perature difference becomes obvious in a stress-stretch diagram, see Fig.1. It has to be remarked, that the function ϕ̂(ϑ) defined in Eq. (2.6) with α =2.06 ·10−4K−1 (see Heimes, 2005), has no essential influence in the range of applications, so that the behavior is close to linear behavior. Fig. 1. Simple tension for constant temperatures (Θ0 =273K). Representation of the linear temperature dependence. σR defines the 1st PK-stress component As mentioned in Hartmann and Neff (2003), the strain-energy part U(JM)= K/2(JM−1)2, commonly applied in commercial finite element com- putations, yields in the case of uniaxial tensile tests an increase of the lateral stretch above a certain amount of the lateral stretch (i.e. the specimen be- comes thicker in the tensile range). Vice versa, in compression, the specimen becomes thinner below a certain compressive stress state. Accordingly, the investigation of the sensitivity of JM = J/ϕ̂(ϑ) is of interest. In Fig.2a the expected non-physical behavior is shown.Again, there is no essential influence Theoretical and numerical aspects in weak-compressible... 17 of the temperature-dependence introduced by ϕ̂(ϑ). All curves are very close to each other. Proposal (2.26)1 does not show such non-physical behavior, see Fig.2b. Unfortunately, there does not exist an analytical proof guaranteeing this property. Only computations in a wide range of stretches and parameters indicate this behavior. Fig. 2. Lateral stretches for various strain-energy functions U(JM)= U(J/ϕ) with K =1000MPa; (a) U(JM)= K/2(JM −1)2, (b) U(JM)= K/50(J5M +J −5 M −2) 4.2. Finite element computations The proposed solution scheme using Rosenbrock-type methods is demon- strated at one-eights of an elastomeric specimen shown in Fig.3. Thematerial parameters are those of the previous example, and for the thermal part we define κ = 0.26W/mK and Θ0 = 293.15K. Additionally, convection is as- sumed at the outer surfaces with q(t) = h0(Θ − Θ∞), h0 = 12.5W/m2K and Θ∞ = 293.15K. Since the heat flux depends essentially on the surface deformation (spatial quantity), the numerical treatment of the deformation- dependent heat flux has to be taken into considerations (this is not discus- sed here). The mechanical and thermal boundary conditions are described in Fig.3. One simplification is the negligence of the first term in Eq. (3.3), which is very small for J ≈ 1.With the thermal part of free-energy (2.27), one obta- ins cp = cp0(1+kpΘ), where cp0 =1.54 ·103J/kgK and kp =3.75 ·10−3K−1 are assumed (see Heimes, 2005). In order to study the step-size behavior, the specimen is loaded sinu- soidally with 200 cycles in the first 200s. Afterwards, the displacements uz(t) = −10mm, t ­ 200s, are held constant, Fig.4a. A step-size control- 18 A.-W. Hamkar, S. Hartmann Fig. 3. Evaluations points, mechanical and thermal boundary conditions for a mesh using eight-nodedQ1P0-elements; (a) evaluation points, (b) mechanical boundary conditions, (c) thermal boundary conditions, q(t)= h0(Θ−Θ∞) Fig. 4. Loading path and the resulting step-size behavior; (a) sine loading path with 200 cycles, (b) step-size behavior of the integrationmethod led third-order method of Rang and Angermann (2008) with an embedded method of order two to estimate the local error is applied to resolve the com- plicated temperature-deformation behavior. The temperature increases in a non-sinusoidal manner, see zoomed picture in Fig.5b, which is a known fact (see Reese, 2001). Moreover, different time-scales (e.g. the final hold-time of 200s) cannot be efficiently solved with constant step-sizes so that a step-size control technique is recommendable. Figure 5a shows the increase of the tem- perature in the specimen during cyclic loading which has the tendency to reach saturation caused by the influence of convection. During the hold-time 200¬ t ¬ 400s, the temperature decreases caused by convection as well. Theoretical and numerical aspects in weak-compressible... 19 Fig. 5. Temperature evolution; (a) temperature evolution at marked nodes, see Fig.3a, (b) temperature evolution at marked node 5 In Fig.6 it is shown, at different times, that the heat evolution requires a good spatial resolution close to the surface, where a finer mesh is used. Fig. 6. Temperature evolution for time t =0,100,200,400s 5. Conclusions A consistent decomposition into thermal and mechanical as well as isochoric and volumetrical effects is applied to the case of thermo-hyperelasticity. A particular poly-convex strain-energy function based on invariants of the uni- 20 A.-W. Hamkar, S. Hartmann modular Cauchy-Green tensors is extended to the thermo-mechanical case. Even in this case, it is demonstrated that a physical meaningful uniaxial ten- sile/compression case requires specific strain-energies. Thesemodels are incor- porated in a fully thermo-mechanical coupled finite element program,where a totally iteration-free third-orderRosenbrock-typemethod to solve the coupled DAE-system given by the heat equation and equilibrium conditions offers a very efficient possibility to solve cyclic loading with different time-scales. Acknowledgment This paper is based on investigations of a grant given by the Simulationswis- senschaftliche Zentrum of the Clausthal University of Technology. We would like to express our sincere thanks. References 1. EhlersW., EipperG., 1998,The simple tension problemat large volumetric strains computed from finite hyperelastic material laws,Acta Mechanica, 130, 17-27 2. Flory P.J., 1961, Thermodynamic relations for high elastic materials,Trans- action of the Faraday Society, 57, 829-838 3. Gear C.W., 1986, Maintaining solution invariants in the numerical solution of ODEs, SIAM Journal on Scientific and Statistical Computing, 7, 3, 734-743 4. Hairer E., Wanner G., 1996, Solving Ordinary Differential Equations II, Springer, Berlin, 2nd revised edition 5. Hartmann S., 2001a, Numerical studies on the identification of the material parameters of Rivlin’s hyperelasticity using tension-torsion tests,Acta Mecha- nica, 148, 129-155 6. Hartmann S., 2001b, Parameter estimation of hyperelasticity relations of ge- neralized polynomial-type with constraint conditions, International Journal of Solids and Structures, 38, 44/45, 7999-8018 7. Hartmann S., 2002, Computation in finite strain viscoelasticity: finite ele- ments based on the interpretation as differential-algebraic equations,Computer Methods in Applied Mechanics and Engineering, 191, 13/14, 1439-1470 8. Hartmann S., Hamkar A.-W., 2010, Rosenbrock-type methods applied to finite element computations within finite strain viscoelasticity, Computer Me- thods in Applied Mechanics and Engineering, 199, 23/24, 1455-1470 9. Hartmann S., Neff P., 2003, Polyconvexity of generalized polynomial-type hyperelastic strainenergy functions fornear-incompressibility, International Jo- urnal of Solids and Structures, 40, 11, 2767-2791 Theoretical and numerical aspects in weak-compressible... 21 10. Hartmann S., Quint K.J., Hamkar, A.-W., 2008,Displacement control in time-adaptive non-linear finite-element analysis, Journal of Applied Mathema- tics and Mechanics, 88, 5, 342-364 11. Hartmann S.,Wensch J., 2007,Finite element analysis of viscoelastic struc- tures using Rosenbrock-typemethods,Computational Mechanics, 40, 383-398 12. HauptP., 1985,On the concept of an intermediate configurationand its appli- cation to representation of viscoelastic-plasticmaterial behavior, International Journal of Plasticity, 1, 303-316 13. Haupt P., 2002, Continuum Mechanics and Theory of Materials, Springer, Berlin, 2 edition 14. Haupt P., Tsakmakis C., 1989, On the application of dual variables in con- tinuummechanics, Journal of ContinuumMechanics and Thermodynamics, 1, 165-196 15. Haupt P., Tsakmakis C., 1996, Stress tensors associated with deformation tensors via duality,Archive of Mechanics, 48, 347-384 16. Heimes T., 2005, Finite Thermoinelastizität, Number 709 in Fortschrittsbe- richte, Reihe 5, Grund- und Werkstoffe/Kunststoffe, VDI-Verlag, Düsseldorf 17. Holzapfel G., Simo J., 1996a, A new viscoelastic constitutive model for continuousmedia at finite thermomechanical changes, International Journal of Solids and Structures, 33, 3019-3034 18. Holzapfel G., Simo J., 1996b, Entropy elasticity of isotropic rubber-like so- lids at finite strains,ComputerMethods in AppliedMechanics andEngineering, 132, 17-44 19. Lang J., 2000, Adaptive multilevel solution of nonlinear parabolic PDE sys- tems. Theory, algorithm, and applications, Springer, Berlin 20. Lion A., 1997, A physically basedmethod to represent the thermomechanical behaviour of elastomers,Acta Mechanica, 123, 1-26 21. LionA., 2000,Thermomechanik von Elastomeren. Experimente undMaterial- theorie, Habilitation, Institute of Mechanics, University of Kassel, Report No. 1/2000 22. Lu S., Pister K., 1975, Decomposition of deformation and representation of the free energy function for isotropic thermoelastic solids, International Journal of Solids and Structures, 11, 927-934 23. Lubich C., Roche M., 1990, Rosenbrock methods for differential-algebraic systems with solution-dependent singular matrix multiplying the derivative, Computing, 43, 325-342 24. Miehe C., 1988, Zur numerischen Behandlung thermomechanischer Prozes- se, Report No. F88/6, University of Hannover, Institut für Baumechanik und NumerischeMechanik 25. Miehe C., 1995, Entropic thermoelasticity at finite strains. Aspects of the formulation and numerical implementation,ComputerMethods in Applied Me- chanics and Engineering, 120, 243-269 22 A.-W. Hamkar, S. Hartmann 26. Rang J., Angermann L., 2008, New Rosenbrock methods of order 3 for PDAEs of index 2, Advances in Differential Equations and Control Processes, 1, 2, 193-217 27. Reese S., 2001, Thermomechanische Modellierung gummiartiger Polymer- Strukturen, Habilitation, Institut für Baumechanik undNumerischeMechanik, Universität Hannover, Report No. F01/4 28. SimoJ.C.,MieheC., 1992,Associativecoupled thermoplasticityatfinite stra- ins: Formulation, numerical analysis and implementation, Computer Methods in Applied Mechanics and Engineering, 98, 41-104 29. Simo J.C., Taylor R.L., 1991, Quasi-incompressible finite elasticity in prin- cipal stretches. Continuumbasis and numerical algorithms,Computer Methods in Applied Mechanics and Engineering, 85, 273-310 30. Simo J.C., TaylorR.L., PisterK.S., 1985,Variational and projectionme- thods for thevolumeconstraint infinitedeformationelasto-plasticity,Computer Methods in Applied Mechanics and Engineering, 51, 177-208 31. Treloar L.R.G., 1975, The Physics of Rubber Elasticity, Clarendon Press, Oxford, 3rd edition Teoretyczne i numeryczne aspekty słabo ściśliwego sformułowania termosprężystości dużych deformacji Streszczenie W pracy przedstawiono model niemal nieściśliwego, sprężystego zachowania się materiału, rozszerzającgona efekty termiczne.Napoczątku rozważańdokonanomul- tiplikatywnej dekompozycji gradientu deformacji na część termiczną i mechaniczną. Część termiczna wykazuje charakter czysto objętościowy. Dodatkowo, część mecha- niczną zdekomponowanona element zachowujący objętość i element o zmiennej obję- tościw ten sposób, żewypadkowystannaprężeńwykazujewrażliwośćna temperaturę. Zaproponowanymodel szczegółowo zbadano w kontekście efektów sprzężenia termo- mechanicznego. W dalszej części pracy, analizowany model zastosowano do czaso- wo adaptacyjnej metody elementów skończonych sformułowanej na podstawiemetod Rosenbrocka wyższych rzędów. Takie sformułowanie umożliwia uzyskanie procedury beziteracyjnej, co z kolei pozwala na wykonanie wyjątkowo szybkich obliczeń nume- rycznych. Artykuł zamyka przykład symulacji numerycznej trójwymiarowej próbki elastomeru poddanej próbie rozciągania. Manuscript received December 1, 2010; accepted for print April 1, 2011