Inelastic heat fraction estimation from two successive mechanical and thermal analyses and full-field measurements

European Journal of Mechanics A/Solids 38 (2013) 1e11

Contents lists available at SciVerse ScienceDirect

European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol

Inelastic heat fraction estimation from two successive mechanical and thermal analyses and full-field measurements T. Pottier a,1, F. Toussaint a, *, H. Louche b, P. Vacher a a b

Laboratoire SYMME, Polytech’Savoie, BP 80439, 74944 Annecy le Vieux Cedex, France Laboratoire de Mécanique et de Génie Civil (LMGC), Université Montpellier 2, CNRS, Place Eugène Bataillon, 34095 Montpellier Cedex, France

a r t i c l e i n f o

a b s t r a c t

Article history: Received 6 January 2012 Accepted 4 September 2012 Available online 21 September 2012

A new method to identify inelastic heat fraction (b factor) evolution models is proposed in this paper. It is based on: (i) simultaneous kinematic and thermal field measurements on heterogeneous tensile tests and (ii) Finite Element Updating (FEU) inverse method. Two inverse calculations, that involve a LevenbergeMarquardt optimization algorithm, are successively used to identify on the one hand the material parameters of a mechanical constitutive model (anisotropic plasticity coefficients and hardening parameters) and on the other hand the plastic power ratio converted into heat. The power balance of the thermomechanical problem is calculated and presented, and the assessment of thermal heat sources is detailed. Finally, six mechanical parameters and four evolution models of b are successively identified for commercially pure titanium. Results confirm a strain dependency of b. Ó 2012 Elsevier Masson SAS. All rights reserved.

Keywords: Full-field measurement Inverse problem Inelastic heat fraction Infrared thermography Plastic power

1. Introduction Heat generation in a material during a mechanical load offers important knowledge in order to predict temperature variations, mainly in dynamic conditions. In elastoplastic materials, heat generation is associated with intrinsic dissipation due to plasticity. Part of the mechanical power given to an elementary volume of the material is converted into heat power while the remainder is stored (see for example the thorough review of Bever et al. (1973)). The knowledge of the dissipated heat part is both interesting: (i) to estimate the amount of heat power generated by the irreversible behaviours and (ii) to estimate the stored power (or energy) in the material. Heat generation during the deformation can be estimated either by an experimental approach or by a modelling one (Chaboche, 1993a, 1993b; Chrysochoos et al., 1989; Håkansson et al., 2005; Helm, 2006; Longère and Dragon, 2008; Vincent, 2008; Dumoulin et al., 2010; Stainier and Ortiz, 2010). In these two approaches, the intrinsic dissipation is often computed through Inelastic Heat Fraction (IHF), generally noted b. As mentioned by Rittel (1999), the

* Corresponding author. Tel.: þ33 (0)450 096 576; fax: þ33 (0)450 096 543. E-mail address: [email protected] (F. Toussaint). URL: http://www.symme.univ-savoie.fr 1 Present address: Department of Mechanical Engineering, Faculty of Engineering, Chiang Mai University, 239 Huaykaew Rd., Suthep Amphur Muang, 50200 Chiang Mai, Thailand. 0997-7538/$ e see front matter Ó 2012 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.euromechsol.2012.09.002

b factor, introduced and calculated in the pioneer work of Taylor and Quinney (1937), was an integrated factor defined by the ratio of the dissipated energy to the plastic work. This factor is noted bint ¼ b in Rittel (1999). The power ratio, obtained by dividing the heat power (intrinsic dissipation) by the plastic power was denoted bdiff by Rittel (1999). In many Finite Element codes, the b factor was assumed to be constant and typically about 0.9 for metals (Needleman and Tvergaard, 1995; Zhou et al., 1996; Kappor and Nemat-Nasser, 1998). In these particular cases, b ¼ bint ¼ bdiff. However, several experimental studies have shown that this assumption does not hold in many material types (Taylor and Quinney, 1937; Bever et al., 1973; Chrysochoos et al., 1989; Macdougall, 2000; Hodowany et al., 2000; Oliferuk et al., 2004; Chrysochoos et al., 2009; Dumoulin et al., 2010). In these studies, b was different than 0.9 and was strain (and sometimes strain rate) dependent. A summary of such variations in b values can be found in Macdougall (2000). Several experimental techniques have been proposed to measure, at a macroscopic scale, the IHF (b or bdiff) or the stored energy ratio. These techniques aim at measuring the dissipated part of the plastic work (or plastic power) through the following methods: (1) Global energy measurements during quasi-static tests by: microcalorimetry (Bever et al., 1973; Chrysochoos et al., 1989); comparison with a dummy specimen heated electrically (Oliferuk et al., 1995).

2

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11

(2) Local temperature measurements by: embedded thermocouples during dynamic tests (Rittel, 1999); infrared radiometers during dynamic tests (Hodowany et al., 2000; Macdougall, 2000); infrared thermography in quasi-static tests (Chrysochoos et al., 1989, 2009; Dumoulin et al., 2010. (3) Inverse method: analytical 3D heat diffusion modelling (Bouffera et al., 2005). The method proposed in this work is in the third type. Different models of bdiff evolution with the strain are chosen a priori and identified through an inverse method. The chosen method is the Finite Element Updating (FEU), broadly used in several applications such as elasticity (Lecompte, 2007), hardening (Kajberg and Lindkvist, 2004) and anisotropy parameters identification (Meuwissen, 1998), but also many others. An overview of the proposed method is given in the second section and the constitutive models to be identified and the heat sources assessment procedure are discussed in Section 3. Then, the experimental setup is presented in Section 4. The numerical model and the boundary conditions are detailed in Section 5. Finally, results and discussion are presented in the last section. 2. Overview of the method In thermomechanical analysis, the material self-heating is often simulated in the FE calculations as a conversion rate of the volume plastic power w0p into heat:

b¼

w0d w0 ¼ 1  0s w0p wp

(1)

where w0d is the dissipated power generated by irreversibilities (plasticity, damage.), w0s is the stored power, and b ¼ bdiff to simplify the notation. This paper aims at applying the FEU method to the identification of four a priori models of the evolution of b when strain varies. For this purpose, two FEU inverse methods are successively run from a restricted number of tests (Fig. 1): (i) The first one denoted FEU-M is a mechanical inverse identification and aims to determine the parameters of the mechanical model. Global measurement (reaction force of the sample F exp(t)) and local data obtained by Digital Image Correlation (DIC) are used to identify both Ludwick’s hardening parameters and Hill’s anisotropy parameters. DIC provides full-field measurements of longitudinal ~ ; tÞ and uexp ~ ; tÞ at the and transverse displacement fields uexp x ðx y ðx ~ . Note that in the present paper superscript exp surface location x denotes quantities assessed from experimental measurement and num those computed through the choice of constitutive equations. (ii) The second one, denoted FEU-T, identifies the parameters of different inelastic heat fraction b models. It takes into account the plastic power w0p, the thermomechanical coupling power (denoted w0tmc ) computed in the previous FEU-M analysis and the experimental thermal field T exp measured by the infrared camera. The present uncoupled approach relies on the assumption that temperature raises measured during the experiment (less than 12 K) are not sufficient to imply a significant material softening. The works of Khan et al. (2004) and Ogi et al. (2004) have shown that elastic and plastic constitutive parameters after such self-heating remain close to their nominal value (estimated at ambient

Fig. 1. Mechanical parameters and inelastic heat fraction identification flowchart.

temperature). Hence, for the sake of simplicity, the thermal and mechanical problems are assumed to be independent and the influences of temperature on the mechanical behaviour and both elastic and plastic moduli are discarded. As shown in Fig. 1, the mechanical modelling provides an estimation of the local stress state through the Cauchy stress tensor (noted Tðx; tÞ) and displacement fields (unum(x,t)) which are used to assess the plastic work. In the FEU-M procedure, the imposed displacements at the domain boundaries are repeated throughthickness. Subsequently, output data (T and unum) are known ~ . Then, the over the sample volume x and not only over its surface x knowledge of the power balance and a value of b are sufficient to assess internal heat sources w0ch according to Eq. (1) and the split of w0ch in a part associated with the thermomechanical coupling w0tmc and a dissipated part w0d (Chrysochoos et al., 1989):

w0ch ¼ w0d þ w0tmc :

(2)

Finally, in the FEU-T procedure, the value of b is optimized in order to minimize the discrepancy between the calculated ~ ; tÞÞ temperature fields. Even (Tnum(x,t)) and the measured ðT exp ðx

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11

though the temperature field is calculated over the sample volume ~ ; tÞÞ are used for comparison x, only the surface data ðT num ðx purpose. 3. Constitutive models 3.1. Thermal model The volume heat sources w0ch , which relate the mechanical to the thermal behaviour, are set in the right hand side of the classical heat diffusion equation. In the current configuration, it therefore comes:

rC q_  kDq ¼ w0ch ;

(3)

where r is the mass density, C the specific heat, k the thermal conductivity and q is the temperature variation defined as q ¼ T  T0 with T0, the initial temperature. As detailed in Rittel (1999), the inelastic heat fraction b can be defined in terms of energy or power ratio and is then associated with an integrated or an instantaneous quantity. In the present study, the power formalism is chosen. According to Eqs. (1) and (2) the heat diffusion equation becomes:

rC q_  kDq ¼ w0ch ¼ bw0p þ w0tmc :

(4)

Due to the strain magnitude in the sample (the maximum longitudinal GreeneLagrange strain is Emax xx z0:8), the hypothesis of small perturbations no longer holds. Hence, both mechanical and thermal constitutive models must be written in the finite deformation framework. For this purpose, a choice upon the gradient tensor decomposition needs to be made. In the present paper, Lee’s decomposition is assumed (Lee, 1969, 1981). It consists in a multiplicative split of the gradient tensor in an elastic and a plastic part as following:

F ¼ Fe Fp :

(5)

This assumption introduces an unstressed configuration (a bar will be added when variables are expressed in this configuration) obtained after the transformation Fp , which is neither Eulerian nor Lagrangian. A significant consequence of the chosen split of the gradient tensor is the subsequent split of the strain rates tensor as: _ 1 ¼ Ge þ Fe Gp F1 . This relation allows to write the G ¼ FF e external volume power w0ext in the unstressed configuration as follows (Chrysochoos, 1985; Sidoroff and Dogui, 2001):

 r
r w0ext ¼ S : FTe Fe Gp þ S : E_ e ; r r |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflffl{zfflfflfflffl} w0p

(6)

w0e

where “:” denotes the double contraction of two second order tensors, w0e is the elastic volume power, Ee ¼ 1=2ðFTe Fe  IÞ is the elastic GreeneLagrange strain tensor, r=r is the volume variation over the elastic transformation. Gp is the velocity gradient tensor associated to Fp such as Gp ¼ F_ p F1 p . Finally, the second Piolae Kirchoff stress tensor S related to the unstressed configuration is estimated, in the matlab environment through:

r S ¼ F1 TF1T ; e r e

(7)

where T is the Cauchy stress tensor. The expressions of w0ext (Eq. (6)) and S (Eq. (7)) provide a way to assess w0p from mechanical outputs. Indeed, the mechanical outputs are reset in the unstressed configuration within Matlab environment and allow to assess w0ext prior to the thermal finite element calculations (performed by the means of Abaqus FE code).

3

For this purpose, the thermomechanical coupling power w0tmc must also be estimated (Eq. (4)). This power is the sum of different terms, including the variation of hardening parameters (thermodynamical forces) with the temperature (Chrysochoos et al., 1989. In the case of titanium, this dependence is supposed negligible and the thermomechanical coupling power is reduced to the only thermoelastic power denoted w0thel . The thermoelastic coupling is related to the elastic transformation Fe (defined from the unstressed to the Eulerian configurations), thus the infinitesimal deformation framework can be satisfied if small elastic strains are assumed. Therefore, the choice of the stress tensor used to evaluate thermoelastic coupling does not significantly influence the value of w0thel . It therefore comes:


 _ ; w0tmc ¼ w0thel z  lT0 tr S

(8)

where l stands for the linear thermal expansion coefficient. According to Boulanger et al. (2004), this assumption holds as long as T0El2/rC  1, and such a condition is verified if the values given in Table 1 are used (T0El2/rC z 1.6  104). As discussed in paragraph 1, b is often chosen to be constant; in this case, only one parameter is to be identified. In the present study, four models that take into account the evolution of b during the deformation process are investigated. The first model is chosen to be the constant model. In accordance with experimental studies (Oliferuk et al., 1995; Chrysochoos et al., 1989; Macdougall, 2000; Hodowany et al., 2000), the three other models are assumed to depend on the plastic strain Ep ¼ 1=2ðFTp Fp  IÞ, the fourth model is the model proposed by Zehnder (1991) and discussed in Appendix A:

8 > model no: 1 : > > > > > > model no: 2 : > > > < model no: 3 : > > > > > > > > model no: 4 : > > :

b ¼ a; 







b Epxx ¼ bEpxx þ c; 



b Epxx ¼ d Epxx e þf ;

(9)



1n ~ h Exx   b Epxx ¼ p 1n ; xx Ep

where (a,b,.,f) are the constant material parameters to be identified. In the first model, a is supposed to be included in [0,1]. In the ~ is the constant Zehnder’s model, n is the hardening exponent and h to be determined. Subsequently, the number of parameters to identify varies between 1 (models no. 1 and no. 4) and 3 (model no. 3). The choice of a dependency of b to a single component of the strain tensor constitutes a strong simplification that relies on the fact that, in the presented experimental case, the tensile loading rules the sample deformation. It was thus assumed and verified that in this specific case only, (i) Epxx zp where p stands for the cumulated plastic strain and that (ii) the time evolution of this component is in accordance with the evolution of the cumulated plastic strain. This latter must be preferred in the general case. 3.2. Mechanical model The solving of the mechanical finite element problem is performed using Abaqus FE code. The calculations output is the Cauchy Table 1 Constant parameters of the study. Parameters

Units Value

Mechanical parameters

Thermal parameters

r

E

n

S0

k

C

l

kg m3 4500

GPa 104

e 0.34

MPa 396

W m1 K1 17

J kg1 K1 520

K1 8.4$106

4

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11

stress tensor T which is used to estimate the principal extension ratios (l1,2,3). T is defined in the spatial configuration and therefore needs to be reset in the unstressed configuration in order to perform the power estimation presented above (Eq. (6)). For this purpose the displacement fields are also used as outputs in order to estimate the strain gradient tensor F. For the sake of consistency, in the present work, the mechanical model is expressed with respect to the unstressed configuration. An anisotropic elasticeplastic mechanical model is chosen. The elastic behaviour is modelled by the Hooke’s law. Elastic isotropy is assumed and leads to a definition of a linear model ruled by two parameters E the Young’s modulus and n the Poisson’s ratio. The plastic behaviour is thus defined in the unstressed configuration. Associated plasticity is assumed and implies the definition of the plastic strain according to the normality rules (Sidoroff and Dogui, 2001; Simo and Ortiz, 1985) such as:

_ 1 F1T p Ep Fp ¼ l


vf S; R; S0 vS

(10)

where E_ p and S are the plastic GreeneLagrange strain rate (reference configuration) and the second PiolaeKirchoff (unstressed configuration). Finally, l is the plastic factor and f and R are respectively the yield function and the isotropic hardening variable. These two latter quantities are now to be defined. i) The yield surface f ðS; R; S0 Þ is here settled according to the Hill’s criterion, it allows the assessment of the plastic strain rate at every stage of the plastic deformation. According to Simo and Ortiz (1985) the yield surface should include the elastic right CauchyeGreen strain tensor relative to the unstressed configuration. However, under the assumption of small elastic strains, this tensor should be approximated by the identity tensor. Finally, the yield function is:


  2 2 2 1
 f S; R; S0 ¼ pffiffiffi F S22  S33 þG S33  S11 þH S11  S22 2 12
2 2 2 þ 2 LS23 þ MS31 þ NS12 R  S0 ; ð11Þ where Sij are the components of the second PiolaeKirchoff stress tensor in the unstressed configuration S. S0 is the initial yield stress and is chosen to equal the yield stress determined from a tensile test performed at 0 from the rolling direction. Such an assumption classically allows to fully constraint the Hill’s criterion and therefore imposes that G þ H ¼ 1 (Toussaint et al., 2008b; Pottier et al., 2011). ii) A Ludwick’s hardening law is finally chosen to complete the plastic model definition and is written as:

R ¼ KðpÞn ;

(12)

where R is the isotropic hardening variable, K and n are the material parameters to identify. p is the cumulated plastic strain. However, under the assumption of monotonous load cases, which are assumed in the present work, it can be stated that:




pz

eq Ep



Z t rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2_ ¼ Ep ðtÞ : E_ p ðtÞdt; 3

(13)

0 eq

where Ep is the equivalent GreeneLagrange plastic strain. The plastic strain is defined by Ep ¼ 1=2ðFTp Fp  IÞ. The isotropy of the hardening behaviour has been experimentally verified from standard identification carried out on three tensile samples cut at 0 , 45 and 90 from the rolling direction (such an identification procedure is presented in Toussaint et al. (2008a) and Pottier et al. (2011)). The complete material model basically requires eleven parameters to be identified (K, n, F, G, H, L, M, N, E, n and S0). However, Young’s modulus E, Poisson’s ratio n and yield stress S0 are assumed to be known and constant (Table 1). Moreover, it has been verified that the sensitivity of the experiment to L and M is too small to enable their identification. Thus, they have been fixed to their isotropic values. Therefore, only six parameters remain to be identified: two for the hardening law (K and n) and four for the yield criterion (F, G, H, and N). 4. Experimental setup and full-field measurements 4.1. Experimental setup The tested samples are cut from CP titanium rolled sheets (thickness: 1.6 mm). The rolling and the transverse directions are respectively denoted 0 and 90 . A notched geometry (Fig. 2) is chosen, according to the work of Meuwissen (1998), so that the obtained strain fields exhibit a shear band in the central zone of the sample. Tests are performed on an Instron 50 kN tensile device that undergoes a constant cross-head velocity of 5 mm/min. 4.2. Kinematic field measurements Digital Image Correlation (DIC) is used for kinematic measurement purpose. The exposed face of the sample is speckled using matte black and white painting sprays. A Nikon D200 camera captures frames at 0.33 Hz with a resolution of 3872  2592 pixels. DIC is processed using 7D software (Vacher et al., 1999). The dimension of the grid is set to 16  16 pixels and the dimension of the pattern used to compare sub images is 16  16 pixels. The

Fig. 2. Experimental setup with kinematic and thermal full-field measurement devices facing different sides of the sample.

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11

5

Fig. 3. Interpolation of the mechanical boundary conditions from the DIC analysis at the domain’s borders.

displacement field has a bi-linear form and the grey level interpolation is bi-cubic. Furthermore, DIC uncertainties at small strains (less than 5%) are generally admitted to remain below 0.01 pixel (Bornert et al., 2008). At finite strain, the work conducted by Coudert (2005) studied the behaviour of 7D software by the means of numerically deformed images (up to Eeq z 1). For a pattern size of 10  10 pixels, displacement uncertainties were estimated at 0.008 pixel (or 8  104 in terms of strain).

the temperature field to become static. Subsequently, temperatures are known over the undeformed geometry of the sample. Thus, each pixel provides temperature data of one and only material point, regardless of the deformations undergone by the material. Such a motion compensation technique is detailed in Pottier et al. (2009) and is based on the approach of Sakagami et al. (2008).

4.3. Thermal field measurements

5.1. Numerical model

A fixed infrared camera (CEDIP Jade III MW) is used to capture thermal fields (Fig. 2). The exposed face of the sample is painted using black matte painting spray which ensures the emissivity of the sample to remain high and assumed to be close to 1 (Schlosser, 2008; Saai et al., 2010). The capture frequency is set to 10 Hz and frame size is 320  240 pixels. The detector wavelength range is in between 35 mm and the thermal measurement uncertainties are estimated to remain below 0.02 K in the case of relative temperatures q ¼ T  T0 (CEDIP, 2000). Finally, the camera is previously warmed up for about one hour. The sample motion within the camera frame during the deformation process prevents retrieving the temperature of a material point. Taking into account the displacement fields obtained through the DIC analysis, the temperature field is corrected in order to cancel the sample motion due to deformation. Displacements are interpolated at every pixel location and subtracted for

A 3D Finite Element model is built to duplicate the experimental geometry. The only visible part of the sample is modelled. In other words, the free ends, hidden in the grips, are not considered. The part is meshed using 7056 solids elements (3 layers). C3D8R elements are used in the mechanical analysis and DC3D8 for the thermal problem. In order to improve heat sources calculations, the mesh is refined in the shear band of the sample (Fig. 3). Due to the unavailability of through-thickness measurements (whether from a kinematic or a thermal point of view), the surface data are repeated over the 3 layers of elements. Two distinct problems need to be solved successively: the mechanical one and the thermal one. Hence, two cost-functions are defined. The first is needed in the mechanical problem and is then built using only mechanical data: the two components of the displacement and the global force. Abaqus Explicit plasticity model is used to solve the Finite Element problem. The second, defining

a

5. Numerical model and boundary conditions

b

Fig. 4. (a) Samples used for the characterization of thermal boundary conditions. (b) Temperature evolution at the centre of the sample.

6

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11

assessment of the parameter column using the following updating equation:

Table 2 Initial and identified mechanical parameters. Parameters

K

n

F

G

H

N

Initial values Identified values

550 370.0

0.7 0.4698

0.5 0.1547

0.5 0.2765

0.5 0.7235

3 1.722

the thermal problem, is built from sample temperature and is solved using Abaqus Standard Heat diffusion model. The mechanical cost-function fM ðpÞ is expressed as follows:

2

20 12 0 12 exp exp num Nt ;Ns ;Nn  ux; jk C Buy; jk  unum ðpÞC y; jk 6 X 6Bux; jk ðpÞ
 A þ@
 A fM ðpÞ ¼ 4 4@ max uexp max uexp i;j;k ¼ 1 x; jk y; jk jk

jk

1 2 3 31 2 exp num ðpÞ F  F B j 77
j  C þNn @ ; 5 5 A max Fjexp

(16)

where r(p) is the residual column constituted by the relative field ~ ; tÞ, uy ðx ~ ; tÞ and F(t) for the FEU-M and differences at each node (ux ðx qðx~ ; tÞ for the FEU-T), lLM is a damping parameter introduced by Levenberg (1944) to stabilize the GausseNewton method near the optimum, J is the Jacobian matrix. This latter is assessed by computing the partial derivative, along each constitutive parameter, of every single system response. The Jacobian is here estimated by the means of a forward finite differences scheme (Meuwissen et al., 1998). 5.2. Mechanical boundary conditions

0

j


 
J T J þ lLM I pðkþ1Þ  pðkÞ ¼ J T rðpÞ;

(14)

i

where p is the parameter column, Nt is the number of tests, Ns is the number of considered time steps and Nn is the number of nodes of the FE mesh. In the thermal problem, only one test is considered (Nt ¼ 1). The thermal cost-function fT ðpÞ is expressed as follows:

In the present study, the FE model is solved using an imposed displacement procedure. As shown in Fig. 3, displacements at the domain’s boundaries (upper and lower element line) are imposed and their values are assessed through an interpolation of the DIC analysis results. Both longitudinal and transverse displacement are interpolated and imposed. Finally, due to the low frequency of the visible camera, boundary conditions are only available at 15 time steps over the deformation process. Subsequently, residual fields are computed at these exact time steps.

2

0 12 31 2 exp num Ns ;Nn 6 X Bqjk  qjk ðpÞ C 7 fT ðpÞ ¼ 4 @ A 5; exp j;k ¼ 1

qj

5.3. Thermal boundary conditions

(15)

where qj is the mean of the measured temperature variation over the sample at the jth time step. Whether a mechanical or a thermal problem is considered, a LevenbergeMarquardt (Marquardt, 1963) algorithm is used to minimize the cost-function. Such a method allows the iterative exp

In order to solve the thermal problem, a transient thermal FE model is built. The nodes locations are unchanged from the mechanical analysis. Heat losses due to convection with the ambient air through the front and the back faces of the sample are modelled using a surface heat flux fface ¼ hface $q, where hface is the film coefficient. In addition, thermal losses through the grips are also modelled using a surface heat flux defined as fgrip ¼ hgrip $q,

a

b

Fig. 5. (a) Force versus displacement for experiments performed at 0 and 90 . Reproductions of the identified model. (b) Relative error on longitudinal displacement, after identification, for experiment at 0 .

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11

7

where hgrip stands for an equivalent film coefficient between the sample and the grips. Finally, the global heat losses balance is: 0 Wloss ¼ fface Sface þ fgrip Sgrip ;

(17)

0 where Wloss are the global heat losses power, Sface and Sgrip are respectively the areas of the contact zones with the ambient air and the grips. The values of hface and hgrip are then evaluated using the method presented by Louche (1999). A first sample (denoted no. 1 in Fig. 4a) is cut at both ends so that the obtained exchange surface matches the convection surface of a normal test. This sample is left to hang freely in the ambient air. The sample no. 2 respects the defined geometry such as used during the tests. This sample is fastened in the grip as usual. Both are homogeneously heated up to 38  C by the means of a heat gun and their return to thermal equilibrium is captured using the IR camera (Fig. 4b). Fig. 4b shows the thermal responses of both samples. In this case, where heat sources are null and temperature variation smalls, a simplified version of the heat equation, with a single time constant s for modelling the heat losses, can be used (Chrysochoos et al., 1989). The observed cooling in one point can then be reasonably fitted using:

st

qðtÞ ¼ q0 $e ;

(18)

where q0 is the initial value of q, t the time and s a time constant. Thus, two time constants are assessed by this mean. The first, denoted sface , represents the heat losses due to convection exchanges with the ambient air. The second represents the global thermal losses and is denoted sglobal . Assuming that the global losses can be split as a sum of convective loss and through-grips loss, it comes:

1

sglobal

¼

1

sface

þ

1

sgrip

;

(19)

where sgrip is the time constant associated with the heat losses through the grips. Therefore, sgrip can be assessed using Eq. (19). Moreover, the relation between the time constants (sface , sgrip ) and the convective exchange coefficients (hface and hgrip) gives (Chrysochoos and Louche, 2000):

sface ¼

rCe 2hface

and

sgrip ¼

rCL ; 2hgrip

(20)

where e ¼ 1.6 mm and L ¼ 66 mm are respectively the thickness and the length of the sample. This latter relation provides the values of hface ¼ 13 W m2 C1 and hgrip ¼ 286 W m2 C1.

Fig. 7. Evolution of power balance during the test.

6. Results and discussion 6.1. Mechanical results The six parameters of the model are identified from two tests carried out with samples cut at 0 and 90 from the rolling direction. The convergence criterion is set to stop the iterative process if the highest parameter update over a step is inferior to 1% of the parameter current value. In practice, parameters of Eq. (14) are Nn ¼ 7056, Ns ¼ 15 and Nt ¼ 2. The values of identified mechanical constants are summarized in Table 2. Fig. 5a shows the experimental and calculated reaction forces for the tests at 0 and 90 . It can be noticed that the optimization process leads to a good fitting of the forces. Fig. 5b shows the fields of the relative displacement errors after identification at five time steps along the deformation process. Even though the optimization process cannot be made to reach the zero of the residual field, the error remains low and in the range of other published inverse identifications using the same material constitutive equations (Lecompte, 2007; Meuwissen et al., 1998). As presented in Section 3.1, the displacement and the stress fields calculated with the optimized parameters are used to assess the volume plastic power w0p (through Eq. (6)) and the volume thermoelastic coupling power w0thel (through Eq. (8)). Thus, a given value of b allows the estimations of the heat sources w0ch (Eq. (4)). Fig. 6 shows the time evolution of the volume plastic power in the sample. The heat sources are located in the reduced section of the sample with higher intensities near the edges.

Fig. 6. Evolution, during loading, of the calculated plastic power fields (in W m3). Numbers 1e5 correspond to the time steps defined in Fig. 5a.

8

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11

Table 3 Identified thermal parameters of the four models of b (Eq. (9)).

Parameter

Model no. 1 Model no. 2

Model no. 3

a

d

Initial values 1 Identified values 0.65

0.3 0.79

0.3 0.43

0.3 1 0.3 1 0.99 0.21 0.01 96  103

ffina

1.910

2.199

1.704

f

~ h

c

a

e

Model no. 4

b

1.810

Final value of the cost function.

Integrated powers over the sample volume (denoted W) of the different volume powers (w0e , w0p , w0thel and w0ext ) are computed. The time evolution of these global powers are plotted in Fig. 7 and compared with the experimental external power used to deform the sample: 0exp

Wext

¼ F$u_ mg :

(21)

where umg is the measured mobile grip displacement from the DIC analysis and F is the global force coming from the load cell. As expected, the sum of local behaviour is able to retrieve the global power undergone by the sample during the test. 6.2. Thermal results The current value of the b coefficient, calculated from one of the models presented in Eq. (9), is used to assess the heat sources in each element of the FE mesh using Eq. (4). Boundary conditions are applied following the considerations presented in Section 5.3. Hence, the thermal FE model provides the resulting temperature field that is to be compared with the experimental measurements. The thermal problem is solved for a single test performed at 0 (Nt ¼ 1). As defined for the mechanical analysis, the convergence criterion is set to stop the iterative process if the highest parameter update is inferior to 1%. Identified parameters of the four models are gathered in Table 3. The models no. 3 and no. 4 exhibit the lowest values of the cost function and are thus the most able to predict the measured

temperature. Conversely, model no. 1 is certainly not sufficient to model the observed phenomenon. Moreover, the evolution of the b ratio versus strain is plotted in Fig. 8a. Since the ratio of plastic power converted into heat is not defined for Epxx ¼ 0, the plot has been chosen to start at Epxx ¼ 0:05. The evolution of the global stored power Ws0num, obtained by a spatial integration of w0s over the sample volume, of the 4 identified models plotted in Fig. 8b exhibits a significant behaviour modification along with the chosen models. More specifically, the model no. 1 predicts a constant increasing of the stored power while the three other models predict a decrease of the stored power and hence a reduction of the stored energy evolution at high stains in accordance with physical considerations (Bever et al., 1973). Fig. 8c shows the temperature evolution predicted by the four models and measured at the sample surface for 3 points (P1, P2 and P3) defined in Fig. 6. Fig. 8c confirms the above considerations among which model no. 3 is the most able to predict the observed temperature field. The temperature matching error ematch(x,t) is calculated from the relative difference of calculated and measured temperature field.

ematch ðx; tÞ ¼

T num ðx; tÞ  T exp ðx; tÞ ; T exp ðx; tÞ

(22)

where T num(x,t) is the calculated temperature and T exp(x,t) the measured one. This error is plotted in Fig. 9 for the model no. 3 at five time steps previously defined in Fig. 5a. The value of ematch(x,t) remains in between 5% and 5% during the test, and is mainly located in the shear band vicinity. Fig. 10 presents the matching error evolution along a central vertical line of the sample (namely line L1 defined in Fig. 6). Such a comparison of the four models highlights that models no. 1 and no. 2 provide a poor reproduction of the observed phenomenon while models no. 3 and no. 4 exhibit a better matching error. Thus, considering the presented results, the evolutionary value of b along strain seems obvious. However, the chosen model to represent this evolution is of the utmost importance. For instance, the model no. 2 does not bring major improvements compared to the constant model assumed in many former studies.

Fig. 8. (a) Predicted evolutions of b for the 4 identified models. (b) Evolution of the stored power for the 4 identified models. (c) Comparison with the experimental response of the predicted temperature for the 4 identified models along the central line of the sample (points denoted P1, P2 and P3 are given in Fig. 6).

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11

9

Fig. 9. Temperature matching error in % for 5 time steps of the model no. 3 (test at 0 ).

Up to the maximum load (time around 33 s), where the Considère criterion indicates a strain localization, the best temperature reproductions are obtained for model no. 3. This model reduces the strain increase of b at the beginning of the loading, in comparison to the evolution predicted by model 4 (Fig. 8a). Non-monotonous evolutions of b, as observed for instance in titanium (Hodowany et al., 2000) or steel material (Dumoulin et al., 2010), can’t be predicted by model no. 3. However, one can easily use the same approach, with another strain dependency of b able to generate a non-monotonous strain evolution of b. Moreover, one can also add a stain rate dependency on b. Close to the maximum load, and after, strain localization is observed. The Finite Element solution is mesh dependent and reproducing the local plastic power is difficult. Then, identification of the different models of b has to be limited to the period prior to the maximum load. Since the damage phenomenon becomes significant after the Considere’s point and are not taken into account into the constitutive model, the behaviours reproduced after 33 s are assumed to be non-relevant (this zone has been filled in grey in Figs. 8c and 10).

6.3. Influence of the model definition In the previous section four evolutionary models of b have been identified from a single test carried out at 0 . However, as discussed above, these models depend on the only longitudinal plastic strain component which may not be relevant in the case of non-trivial loadings. The point is here to emphasis the weight of such an assumption. Therefore, the identified model no. 3 is used to estimate the thermal fields of the experiment performed at 90 . The mechanical parameters remain unchanged and the FE thus provides the stress and displacement fields for the test at 90 . These fields are used to evaluate the heat sources that are the input data of a new thermal calculation. The resulting temperature field is then compared to the measured fields (Fig. 11). It can be seen that the matching error remains in the same order of magnitude that in the case of the test carried out at 0 . However, the accuracy of the reproduction decreases in the high localization zone. In this zone, the measured strain fields differ and lead to a mismatch of the model (identified from the test at 0 ) with the thermal behaviour of the test performed at 90 . This point

Fig. 10. Timeespace plot of the absolute value of the matching error for the 4 investigated models. Temperature evolution is plotted along the line L1 versus time. The loading curve is overlaid.

10

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11

Fig. 11. Temperature matching error in % for 5 time steps of the model no. 3 (test at 90 ).

emphasizes the need of defining the evolution models of b from more significant component of strain (equivalent strain or cumulated strain) and possibly from strain rate.

w0s h ¼ ; KM w0p

7. Conclusions

where h is a constant and, in the large strain framework, the modulus ratio KM is settled as:

In the present paper, a new method to identify the inelastic heat fraction (b) during a quasi-static tensile test is proposed. The method is based on an inverse Finite Element Update approach, coupling a Finite Element code (Abaqus) and a Levenberge Marquardt algorithm, and using simultaneous DIC and IR thermal full-field measurements. This method is first used to identify the mechanical parameters of material constitutive equations by the means of DIC measurements. Then, the obtained results are used to assess the inelastic heat fraction for the same test. Four strain evolution models for b are tested. Up to the maximum load, results show a good agreement between measured and predicted temperature fields. The identified evolution models of b, on a CP titanium material, confirm a strain dependency but remain probably too simple to retrieve the whole energetic process. The extensibility of the identified models is weak since only related to specific geometry and constitutive model. However, the presented approach allows identification of b for various hardening models. Such a validation would strengthen the obtained models of b as making them true material functions. Therefore, an interesting perspective of the presented work would be to investigate the correlation between the hardening model and the identified evolution of b. Furthermore, the obtained values of b are very different from the broadly used value (0.9). The best accordance is obtained with a nonlinear evolution of b with the strain, governed by three parameters (model no. 3). The strain evolution is close to the one given by Zehnder model but leads to better temperature estimations at the beginning of the loading. Improvements of b models, in order to obtain non-monotonous evolution with the strain, may also be proposed with the same approach. In addition, a fully coupled thermomechanical analysis may be considered in further researches since it constitutes a more relevant and close-to-physics approach. Appendix A. Zehnder model The evolution of inelastic heat fraction b with strain such as defined by Zehnder (1991) is related to the ratio of stored power w0s to plastic work rate w0p . It is shown that this ratio exhibits a proportional relation with KM: the ratio of the elastic modulus E to the plastic tangent modulus Ep. It therefore comes:

KM ¼

E E ¼ : Ep vR=vEp

(A.1)

(A.2)

The split of the plastic power into a dissipated and a stored part such as w0p ¼ w0d þ w0s and the definition of the inelastic heat fraction (Eq. (1)) gives:

b¼

w0p  w0s w0d K h ¼ ¼ M : 0 KM wp w0p

(A.3)

The Ludwick relation, chosen in the presented formalism states that:

R ¼ KðpÞn ;

(A.4)

Under the assumption of a monotonous load case it can be eq stated that vpzvEp . Therefore, the assessment of the partial derivative of stress to strain gives another expression of KM:

KM ¼

E E 1n ¼ ; E n1 nK p;eq nKEp;eq

(A.5)

which can be used to provide the evolutionary model of the inelastic heat fraction b:

E 1n 1n  h ~ Ep;eq  h Ep;eq   nK b Ep ¼ ¼ ; 1n E 1n Ep;eq Ep;eq nK

(A.6)

~ ¼ nKh=E is another constant. where h References Bever, M., Holt, D., Titchener, A., 1973. The Stored Energy of Cold Work, third ed. Pergamon Press, Oxford. Bornert, M., Brémand, F., Doumalin, P., Dupré, J.-C., Fazzini, M., Grédiac, M., Hild, F., Mistou, S., Molimard, J., Orteu, J.-J., Robert, L., Surrel, Y., Vacher, P., Wattrisse, B., 2008. Assessment of digital image correlation measurement errors: methodology and results. Experimental Mechanics 49, 353e370. Bouffera, R., Pron, H., Henry, J.-F., Bissieux, C., Beaudoin, J.-L., 2005. Study of the intrinsic dissipation associated to the plastic work induced by a ball impact. International Journal of Thermal Sciences 44, 115e119. Boulanger, T., Chrysochoos, A., Mabru, C., Galtier, A., 2004. Calorimetric analysis of dissipative and thermoelastic effects associated with the fatigue behavior of steels. International Journal of Fatigue 26, 221e229.

T. Pottier et al. / European Journal of Mechanics A/Solids 38 (2013) 1e11 CEDIP, 2000. Jade 3 MW IR Technical Specifications. www.flir.com. Chaboche, J.-L., 1993a. Cyclic viscoplastic constitutive equations, part I: a thermodynamically consistent formulation. Journal of Applied Mechanics 60, 813e821. Chaboche, J.-L., 1993b. Cyclic viscoplastic constitutive equations, part II: stored energy, comparison between models and experiments. Journal of Applied Mechanics 60, 822e828. Chrysochoos, A., Louche, H., 2000. An infrared image processing to analyse the calorific effects accompanying strain localisation. International Journal of Engineering Science 38, 1759e1788. Chrysochoos, A., Maisonneuve, O., Martin, G., Caumon, H., Chezeaux, J.-C., 1989. Plastic and dissipated work and stored energy. Nuclear Engineering and Design 114, 323e333. Chrysochoos, A., Wattrisse, B., Muracciole, J.-M., El Kaïm, Y., 2009. Fields of stored energy associated with localized necking of steel. Journal of Mechanics of Materials and Structures 4, 245e262. Chrysochoos, A., 1985. Bilan énergétique en élastoplasticité grandes dédormations. Journal de Mécanique Théorique et Appliquée 4, 589e614. Coudert, T., 2005. Reconstruction tridimensionnelle du volume intérieur d’une chaussure: évaluation du chaussant. Ph.D. thesis, Université de Savoie, France (in french). Dumoulin, S., Louche, H., Hopperstad, O., Børvik, T., 2010. Heat sources, energy storage and dissipation in high-strength steels: experiments and modelling. European Journal of Mechanics e A/Solids 29 (3), 461e474. Håkansson, P., Wallin, M., Ristinmaa, M., 2005. Comparison of isotropic hardening and kinematic hardening in thermoplasticity. International Journal of Plasticity 21, 1435e1460. Helm, D., 2006. Stress computation in thermoviscoplasticity. International Journal of Plasticity 22, 1699e1727. Hodowany, J., Ravichandran, G., Rosakis, A.J., Rosakis, P., 2000. Partition of plastic work into heat and stored energy in metals. Experimental Mechanics 40 (2), 113e123. Kajberg, J., Lindkvist, G., 2004. Characterization of materials subjected to large strains by inverse modeling based on in-plane displacement fields. International Journal of Solids and Structures 41, 3439e3459. Kappor, R., Nemat-Nasser, S., 1998. Determination of temperature rise during high strain rate deformation. Mechanics of Materials 27, 1e12. Khan, A., Suh, Y., Kazmi, R., 2004. Quasi-static and dynamic loading responses and constitutive modeling of titanium alloys. International Journal of Plasticity 20, 233e248. Lecompte, D., 2007. Elastic and Elasto-plastic Material Parameter Identification by Inverse Modeling of Static Tests Using Digital Image Correlation. Ph.D. thesis, Koninklijke Militaire School, Belgique. Lee, E., 1969. Elasticeplastic deformation at finite strain. Journal of Applied Mechanics 36, 1e6. Lee, E., 1981. Some comment on elasticeplastic analysis. International Journal of Solids and Structures 17, 859e872. Levenberg, K., 1944. A method for the solution of certain non-linear problems in least-squares. Quarterly of Applied Mathematics 2, 164e168. Longère, P., Dragon, A., 2008. Evaluation of the inelastic heat fraction in the context of microstructure-supported dynamic plasticity modelling. International Journal of Impact Engineering 35 (9), 992e999. Louche, H., 1999. Analyse par thermographie infrarouge des effets dissipatifs de la localisation dans les aciers. Ph.D. thesis, Université de Montpellier II, France. Macdougall, D., 2000. Determination of the plastic work converted to heat using radiometry. Experimental Mechanics 40, 298e306. Marquardt, D., 1963. An algorithm for least-squares estimation of nonlinear parameters. Journal of the Society for Industrial and Applied Mathematics 11, 431e441. Meuwissen, M., Oomens, C., Baaijens, F., Petterson, R., Janssen, J., 1998. Determination of the elasto-plastic properties of aluminium using a mixed numericale experimental method. Journal of Material Processing Technology 75, 204e211.

11

Meuwissen, M., 1998. An Inverse Method for Mechanical Characterization of Metal. Ph.D. thesis, Eindhoven University of Technology, Netherlands. Needleman, A., Tvergaard, V., 1995. Analysis of a brittleeductile transition under dynamic shear loading. International Journal of Solids and Structures 32, 2571e 2590. Ogi, H., Kai, S., Ledbetter, H., Tarumi, R., Hi, M., Takashima, N., 2004. Titanium’s hightemperature elastic constants through the hcpebcc phase transformation. Acta Materialia 52, 2075e2080. Oliferuk, W., Swiatnicki, A., Grabski, W., 1995. Effect of the grain size on the rate of energy storage and microstructure evolution during the tensile deformation of austenitic steel. Materials Science and Engineering: A 197, 49e58. Oliferuk, W., Maj, M., Raniecki, B., 2004. Experimental analysis of energy storage rate components during tensile deformation of polycrystals. Materials Science and Engineering: A 374 (1e2), 77e81. Pottier, T., Moutrille, M.-P., Le-Cam, J.-B., Balandraud, X., Grédiac, M., 2009. Study on the use of motion compensation techniques to determine heat sources. Application to large deformations on cracked rubber specimens. Experimental Mechanics 49, 561e574. Pottier, T., Toussaint, F., Vacher, P., 2011. Contribution of heterogeneous strain field measurements and boundary conditions modelling in inverse identification of material parameters. European Journal of Mechanics e A/Solids 30, 373e382. Rittel, D., 1999. On the conversion of plastic work to heat during high strain rate deformation of glassy polymers. Mechanics of Materials 31, 131e139. Saai, A., Louche, H., Tabourot, L., Chang, H., 2010. Experimental and numerical study of the thermo-mechanical behavior of al bi-crystal in tension using full field measurements and micromechanical modeling. Mechanics of Materials 42, 275e292. Sakagami, T., Nishimura, T., Yamaguchi, T., Kubo, N., 2008. A new full-field motion compensation technique for infrared stress measurement using digital image correlation. Journal of Strain Analysis 43, 539e549. Schlosser, P., 2008. Influence des aspects mécaniques et thermiques sur les mécanismes de déformation d’alliages NiTi. Ph.D. thesis, Université Joseph Fourrier, Grenoble, France. Sidoroff, F., Dogui, A., 2001. Some issues about anisotropic elasticeplastic models at finite strain. International Journal of Solids and Structures 38, 9569e9578. Simo, J., Ortiz, M., 1985. A unified approach to finite deformation elastoplastic analysis based on the use of hyperelastic constitutive equations. Computer Methods in Applied Mechanics and Engineering 49, 221e245. Stainier, L., Ortiz, M., 2010. Study and validation of a variational theory of thermomechanical coupling in finite visco-plasticity. International Journal of Solids and Structures 47, 705e715. Taylor, G., Quinney, H., 1937. The latent heat remaining in a metal after cold working. Proceedings of the Royal Society of London A 163, 157e181. Toussaint, F., Tabourot, L., Ducher, F., 2008a. Experimental and numerical analysis of the forming process of CP titanium scoliotic instrumentation. Journal of Material Processing Technology 197, 10e16. Toussaint, F., Tabourot, L., Vacher, P., 2008b. Experimental study with a digital image correlation (DIC) method and numerical simulation of an anisotropic elastice plastic commercially pure titanium. Archives of Civil and Mechanical Engineering 3, 131e143. Vacher, P., Dumoulin, S., Morestin, F., Mguil-Touchal, S., 1999. Bidimensional strain measurement using digital images. Proceedings of the Institution of Mechanical Engineers 213, 811e817. Vincent, L., 2008. On the ability of some cyclic plasticity models to predict the evolution of stored energy in a type 304l stainless steel submitted to high cycle fatigue. European Journal of Mechanics e A/Solids 27, 161e180. Zehnder, A., 1991. A model for the heating due to plastic work. Mechanics Research Communications 18, 23e28. Zhou, M., Ravichandran, G., Rosakis, J., 1996. Dynamically propagating shear bands in impact-loaded prenotched plates e II. Numerical simulations. Journal of the Mechanics and Physics of Solids 44, 1007e1032.