Modeling mechanical response of intumescent mat material at room temperature

Available online at www.sciencedirect.com

International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

Modeling mechanical response of intumescent mat material at room temperature J.S. Kim1 , M.E. Walter∗ , J.K. Lee Mechanical Engineering Department, The Ohio State University, 2075 Robinson Laboratory, 206 W. 18th Ave., Columbus, OH 43210, USA Received 25 October 2001; received in revised form 23 October 2002; accepted 3 November 2002

Abstract Intumescent mat material is widely used to support ceramic substrates in catalytic converters and behaves very much like hyper-foam material under compressive loading. Experiments show that compressive loading curves depend on the ram speed and the number of cycles. The unloading curves show di2erent slopes and paths that depend less on the ram speed and number of cycles. The slopes of the unloading curves decrease as the plastic strain increases; this is referred to as “softening” in this study. The e2ects of rate, softening, and plastic deformation must be considered to model the mechanical response of intumescent mat material. Finite deformation theory is applied with a multiplicative decomposition of the deformation gradient tensor. The developed theory is implemented as an implicit 6nite element algorithm in ABAQUSTM =STANDARD. The necessary material parameters are extracted from experiments. Numerical simulations show good agreement with experiments. ? 2002 Elsevier Science Ltd. All rights reserved.

1. Introduction Energy absorbing support structures are often composed of compressible hyper-foam materials that exhibit viscous reaction and high capacity to retain elastic strain energy for wide ranges of deformation. A hyper-foam material can be de6ned as a highly porous material that can undergo large deformations. One application of a hyper-foam support material is intumescent or “swelling” mat in catalytic converters. The intumescent mat is placed between the catalytic converter substrate and the outer can. The intumescent mat has very little tensile strength and is primarily subjected to compression. While the outer can of a catalytic converter expands and shrinks during the assembly ∗

Corresponding author. 2075 Robinson Laboratory, 206 W. 18th Ave., Columbus, OH 43210, USA. E-mail address: [email protected] (M.E. Walter). 1 Present address: ArvinMeritor, 950W 450S Bldg. 2, Columbus, Indiana, 47201, USA.

0020-7403/02/$ - see front matter ? 2002 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0020-7403(02)00176-5

2286

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

Nomenclature h F  W J p Xi Ei F C e E U   y S D L (4) E∞ Eˆ (4) D(4) (4) T∞ (4) Tem N n PK2

softening parameter yield function density Helmholtz free energy strain energy function determinant of the deformation gradient tensor F principal stretch position vector unit direction vector in the reference con6guration deformation gradient tensor right Cauchy–Green metric tensor Eulerian strain tensor Lagrangian strain tensor right stretch tensor Cauchy stress Kirchho2 stress yield stress 2nd Piola–Kirchho2 stress rate of deformation gradient tensor velocity gradient tensor spatial elasticity tensor without rate e2ects total spatial elasticity tensor material elasticity tensor 4th order transformation tensor between De and De 4th order transformation tensor between De and Dm principal direction vector in the reference con6guration principal direction vector in the current con6guration 2nd Piola–Kirchho2 stress

operation and under service conditions, the compressed mat supports the substrate. In this application, predicting mat holding force for both mechanical loading and unloading is very important in catalytic converter design [1–4]. The intumescent mat material contains aluminosilicate glass 6bers, vermiculite mineral particles, and an organic binder. The vermiculite is a micaceous hydrated magnesium aluminum-silicate mineral which, when heated, looses water and increases the mat layer thickness. At temperatures between 400 and 600◦ C, the increase in thickness can be as high as 10 times. Fig. 1 shows a material structure and Table 1 shows its composition based on mass fraction. Although the mechanical behavior at elevated temperature is quite di2erent from that at room temperature, the process of assembling the converter part takes place at room temperature, and as a 6rst step, it is worth developing a room temperature material response model.

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315 Vermiculite Platelet Vermiculite Platelet

2287

Aluminosilicate Fibers

Aluminosilicate Fibers Thickness (y)

Lateral (x) Before Heating

After Heating

(a)

(b)

Fig. 1. A schematic side view (a) of intumescent mat before heating and after heating. A top view (b) (view from the thickness direction) of a 29 mm diameter mat specimen, where black specks are vermiculite particles. Table 1 Mass based composition and density of intumescent mat Vermiculite

Aluminosilicate

Organic binder

Density (Kg=m3 )

40 –50%

3–10%

31– 43%

632

A great deal of progress has been made on constitutive model development for foam materials e.g., [5–10]. Simo [11] and Papoulia et al. [12] have proposed visco-hyperelastic models to account for viscous e2ects, and rate dependent plasticity formulations are employed to represent plastic deformations [13–16]. Plastic strains are obtained incrementally by satisfying assumed consistency conditions. However, most of these models assume constant unloading slopes and will not provide accurate results when applied to intumescent mat materials under cyclic loading conditions. Fig. 2, a schematic of a quasi-static cyclic compression experiment for intumescent mat, shows a continuous change in the slope of the unloading and reloading curves. The unloading curve shows more softening (a smaller tangent slope) at a higher compressive strain. Although this softening is often referred to as “material damage,” in order to distinguish between similar softening of the long term elastic response, this paper will use the term “softening” to describe the observed decrease in

2288

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315 3 Monotonic loading 1st unloading 1st reloading 2nd unloading 2nd reloading

Compressive stress (MPa)

2.5

2

1.5

1 Material softening

0.5

0 0

0.2

0.4 0.6 0.8 Compressive true strain

1

Fig. 2. The static, cyclic compression of intumescent mat material. At higher compressive strains the unloading curves show increased softening.

tangent slope. The loading and unloading curves also show very di2erent viscous properties when the ram speed is varied in cyclic compression experiments. Therefore, for realistic modeling of intumescent mat materials, in addition to including hyperelastic behavior, it is critical to incorporate the softening and visco-plastic behavior. So far, relatively little research has been done regarding unloading phenomena. The hyper-foam theory cited above typically uses linear elastic or hyperelastic functions for the unloading response. The elastic unloading response function is assumed to be independent of plastic strain and therefore remains constant for all loading histories. With the theories cited above, it would not be possible to model the material response shown in Fig. 2. This study is limited to mat type foam that has a small thickness compared to other dimensions. In addition, the deformation in the thickness direction is much greater than the deformation in the lateral direction. Furthermore, no signi6cant in-plane loading is expected for the application at hand (catalytic converter packaging). Therefore, a 1-D plasticity model in the thickness direction can be created for this material. In this study, to model intumescent mat material response, several experimental results are 6rst reviewed. Next, a constitutive model is extended from hyperelasticity to visco-hyperelasticity. The 1-D plasticity is established based on mechanical strain, and a new elastic response function is introduced to account for material softening due to the strain. The complete spatial elasticity tensor is derived, and the Newton–Raphson method is employed for the plastic strain calculation. The theory is then coded in ABAQUSTM =Standard through the user subroutine UMAT. To verify the theory, the results of several numerical simulations are compared to experiments.

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2289

2. Experimental investigation All the experiments described in this study used 10 mm thick intumescent mat specimens. The specimens were punch cut from a large sheet into 29 mm diameter buttons. All the experiments were performed under compression. They can be divided into the following two main categories: load control (dead weight) compression experiments and displacement control compression experiments. 2.1. Procedures for load control (dead-weight) compression experiments The load control experiment was performed by applying constant pressure on the specimen and recording specimen thickness changes. A custom load frame [17] was built to perform the load control experiments. During the loading process, dead weights were added after the mat was fully relaxed at each weight. For small loads, thickness readings were taken after waiting 30 min. This was determined to be suOcient time for the displacement rate of change to be small enough. In the same way, for the unloading cases, stacked dead weights were removed one by one after no additional mat thickness change was found. Thickness was measured with an LVDT, and the actual load was measured with a load cell. The data was used to construct static stress versus stabilized strain curves. 2.2. Procedures for displacement control compression experiments The displacement control experiments consisted of compressing the specimen to a certain thickness and monitoring how the stress changes as a function of strain and time. A servo-hydraulic load frame was used to compress specimens from 10 mm to a speci6ed position at a 6xed displacement rate. The distance between the loading rods was then held constant for a certain hold time. It was observed that after a hold time of 120 s, the stress was decreasing at a very small rate and that unloading could begin. After the hold time, the upper rod was moved upward to the original 10 mm gap position. Specimens were tested to maximum strains of 0.3, 0.4, 0.5, and 0.6 engineering strain. Experiments were also performed at various ram speeds to investigate rate e2ects. The cyclic experiments employed the same loading pro6le repeatedly. See Kim et al. [17] and Kim [18] for additional details. 2.3. Experimental results Fig. 3 shows stress versus strain curves for di2erent compression speeds up to 60% compressive engineering strain. The static stress versus strain curve comes from the load control experiment with 30 min stress relaxation at each loading. And stress versus strain curves with di2erent ram speeds come from the displacement control experiments. The curves in the 6gure show that the mat becomes sti2er as the compression speed increases and that the mat is highly rate dependent. Fig. 4 shows loading and unloading curves during the 6rst 10 cycles of a 40% cyclic displacement control experiment. There is a signi6cant di2erence between the 1st and 2nd cycle loading curves. Only minor di2erences are seen between the 1st and 2nd cycle unloading curves. Also, hysteresis loops are distinguishable between the 1st cycle and 2nd cycle of loading and unloading. The size of the hysteresis loop decreases from the 1st to 2nd cycle and there is a little further size reduction

2290

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

Compressive stress (MPa)

8

T

v=1.6mm/s v=0.4mm/s v=0.005mm/s Static test

1

6

4

2

T +T 1

2

0 0

0.2

0.4

0.6

0.8

1

Compressive true strain

Fig. 3. Stress versus strain curves to 60% compression for di2erent compression speeds. T1 is the time to reach 60% compression. T2 is the hold time at the maximum strain. T1 + T2 is the time at which displacement unloading begins.

1.2

Compressive stress (MPa)

1

0.8

1st cycle loading

0.6 2nd cycle loading 0.4 3rd cycle loading 0.2 Unloading 0 0

0.1

0.2 0.3 0.4 0.5 Compressive true strain

0.6

Fig. 4. Experimental cyclic compression results for 40% compression (v = 1:6 mm=s).

after the 2nd cycle. The peak stresses from the cyclic 40% displacement control experiments are plotted in Fig. 5. Signi6cant stress reduction takes place from the 1st to 2nd cycle and only gradual stress reductions continue after the 3rd cycle.

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2291

Peak compressive stress (MPa)

1.2

1 v=1.6mm/s v=0.4mm/s

0.8

0.6

0.4

0.2

0 2

4 6 Number of cycles

8

10

Fig. 5. Peak stress variations versus number of cycles for two di2erent compression speeds (40% compression). 3

Monotonic loading Cyclic loading

Compressive stress (MPa)

2.5

2

1.5

1

0.5

0 0.3

0.4

0.5 0.6 0.7 0.8 0.9 Compressive true strain

1

1.1

Fig. 6. Comparison of monotonic static loading and cyclic loading from load control compression experiments.

Fig. 6 shows monotonic and cyclic static stress versus strain curves from the constant load experiments. During the monotonic loading process, load is applied continuously without unloading. The cyclic static loading experiment is performed to investigate the plastic behavior. Load is applied to a certain stress level, and after removing the load, the material is then reloaded to a higher stress level. The direction of loading and unloading is shown in Fig. 6 with arrows. The results show that the unloading curves do not follow the monotonic loading curve in the range of strains we are interested

2292

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315 3.5

Compressive stress (MPa)

3

30% compression 40% compression 50% compression

2.5 2 1.5 1 0.5 0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Compressive true strain

Fig. 7. First cycle compression results for di2erent maximum compressive strains (v = 1:6 mm=s).

in. However, unloading and reloading curves are nearly the same. Note that in Fig. 6, the reloading curve at higher compressive strain is slightly softer than that at lower strain. This softening e2ect is discussed in Section 6.1. All reloading curves appear to rejoin the monotonic loading curve. The displacement control experiment results in Fig. 7 show the 1st cycle of loading and unloading curves for three di2erent maximum compressive strains. All loading curves follow the same loading path. Fig. 8 shows a comparison between the 1st and 2nd cycle loading curves when di2erent compressive strain is applied. The 1st loading path is the same regardless of the maximum applied strain. The path of the 2nd loading is very dependent on the maximum strain reached during the 1st loading. Also, note that the peak stress at the maximum strain level is much less for the second loading. These di2erences between the 1st and 2nd loading cycles suggest that, for a given maximum strain level, the rate dependencies during the 1st loading are di2erent from those during subsequent loadings. The di2erences may be due to damage and plastic deformation incurred during the 1st loading cycle. From these experiments, it is concluded that intumescent mat is a highly nonlinear material with distinct loading and unloading characteristics that depend on the maximum displacement or load, the cycle number, and the rate of loading.

3. Model development Based on experimental results, the four di2erent con6gurations shown in Fig. 9 will be used to account for rate e2ects, material softening, and plastic strain. For a given mechanical deformation Fm , it is assumed that the material reference con6guration 0 6rst changes to the intermediate con6guration 1 due to plastic deformation FP . When the material softens, then the intermediate con6guration 1 changes to another intermediate con6guration 2 by the softening deformation Fh . In this case, the

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

1st loading (40%) 1st loading (50%) 1st loading (60%) 2nd loading (40%) 2nd loading (50%) 2nd loading (60%)

8 Compressive stress (MPa)

2293

6

4

2

0 0

0.2

0.4

0.6

0.8

1

Compressive true strain

Fig. 8. Comparison of the 1st and the 2nd loading curves for di2erent maximum compressive strains (v = 1:6 mm=s).

Current Configuration β Position vector X

F'

Fm

Reference Configuration β0

e

F

Virtual Intermediate Configuration β2 Position vector X2

Position vector X0

Fp

Fh

Intermediate Configuration β1 Position vector X10

Fig. 9. Con6gurational changes due to plastic strain and material softening.

2294

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

intermediate con6guration 1 is a state after plastic deformation and the intermediate con6guration 2 is a state after material softening, that is, dx = F Fh Fp dX0 = Fm dX0 . A visco-hyperelastic model is developed based on the con6gurations and 1 . The material softening model is based on the con6gurations , 1 , and 2 . The elasto-plastic model is based on the con6gurations ; 0 and 1 . Then, through the 4 con6gurations shown in Fig. 9, each model is combined to form a complete visco-hyper-elasto-plastic model with material softening. 3.1. Hyperelastic model development Hyperelastic models have been in existence since the 1940s and are well documented [6–8,10, 19–22]. For the elastic part of the current development, we follow the work reported by Ogden [20] and by Simo and Taylor [21], in which the strain energy is expressed in terms of principal stretches. From the virtual work principal, the internal energy variation is expressed as


e W dV1 = S · E dV1 =  · ee dV; (1) where W is a strain energy function, S is the 2nd Piola–Kirchho2 stress (PK2 stress),  is the Kirchho2 stress, Ee is the Lagrangian elastic strain tensor, ee is the Eulerian elastic strain tensor, V1 is the reference volume, and V is the current volume. From the relation W = S · Ee in Eq. (1), the PK2 stress can be expressed as 3  1 @W 1 1 S= (2) e @e (Ni ⊗ Ni );  i i i=1 where ie is the principal stretch associated with the unit vector Ni1 in the reference con6guration and ⊗ represents the tensor outer product. The PK2 stress can be related to the Cauchy stress  by the elastic deformation gradient tensor Fe by 1 T  = e Fe SFe ; (3) J where J e is the determinant of Fe . The PK2 stress rate can be related to the Lagrangian elastic strain rate by S˙ = D(4) E˙ e . Hence, the fourth order material elasticity tensor D(4) can be obtained by di2erentiating Eq. (2) with respect to time. 3.2. Visco-elastic model development Classical linear viscoelastic models have been incorporated into non-linear large deformation models [11,12,23,24]. The visco-elastic stress response can be viewed as a superposition of the long-term elastic response ∞ and a Prony series (see Papoulia et al. [12] for details): N  H ; (4)  = ∞ + =1

where the long-term elastic response can be derived from the virtual work principle with an appropriately de6ned strain energy function W∞ : @W∞ eT 1 ∞ = e Fe F : (5) J @Ee

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

The Prony series terms are de6ned as  
t 1 e @W∞ eT −1= (t −t
) d H = dt  : a e F F dt  J e @Ee 0

2295

(6)

The constants with subscript  are material constants. Through the above equations, the stress response is divided into static and transient terms. Depending on the complexity of the energy function, the time integral in Eq. (6) may not be analytically tractable. In general, a numerical integration scheme is required to implement these equations in the 6nite element method. Time t in Eq. (6) can be divided into n + 1 time steps so that t n+1 = t n + St. Then, based on the known stress and history at time step n, stress at time step n + 1 can be computed by assuming that ∞ is linear during the time increment St: n+1 n+1 n+1 n+1 = ∞ + vis = ∞ +

N 

Hn+1

(7)

n+1 n (∞ − ∞ ) a  (1 − e−St= ): St

(8)

=1

with Hn+1 = e−St= Hn +

n+1 can be determined from Eq. (5). Since the deformation gradient tensor is known at all times, ∞ All other terms in Eq. (8) are material constants or are known at time step n and therefore Hn+1 can be calculated. The Euler forward integration scheme requires only information at step n to calculate the stress at step n + 1. This is a signi6cant advantage in terms of saving computer memory and computing time. The spatial elasticity tensor Eˆ (4) that includes visco-elastic terms is obtained as (4) ; Eˆ (4) = vis E∞

(9)

where vis = 1 +

N  1 a  (1 − e−St= ): St =1

(10)

(4) maps the elastic rate of stretch tensor to the Jaumann The fourth order spatial elasticity tensor E∞ rate of Cauchy stress, which can be obtained by taking the Jaumann rate of Eq. (5). This procedure is described in Appendix A, and the result is as follows: (4) E∞ = −∞ ⊗ I +

1 (4) P + H(4) : Je

(11)

The components of H(4) and P(4) are related to the Cauchy stress, elastic deformation gradient Fe and material form of the elasticity tensor D(4) , in the following manner (see Appendix A for details): (4) = 12 (ia ∞bj + ∞ia bj + ib ∞aj + ∞ib aj ); Hijab

(12)

(4) (4) e e e e = FiP FjQ FkR FmS DPQRS : Pijkm

(13)

2296

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

4. Softening model development for the elastic unloading function 4.1. Assumptions for the elastic unloading function The experimental results in Fig. 6 show that the unloading stress response softens as compressive strain increases. Constitutive modeling must account for di2erent unloading functions for di2erent compressive strains. Ordinarily this would require many material constants, however, for intumescent mat, experimental observations indicate that the unloading stress responses are self-similar. Let the stress be a function of strain at an arbitrary base state,  = f(U):

(14)

The similarity can be treated by scaling the strain by a factor h, while keeping stress constant:  = f(hU) = f(U ):

(15)

The scaling factor h is called the softening parameter and is assumed to be a function of the 3rd m m component of the right Cauchy–Green metric tensor, i.e., h = h(C33 ). Physically, C33 is the thickness direction stretch ratio which can be used as a thickness strain measure. The stress-strain curve obtained in this manner is called the “master curve”. 4.2. Material softening in large deformation theory Let us imagine two di2erent bodies that deform independently to exactly the same shape and size. Furthermore, assume that one body deforms to the 6nal state by the deformation gradient tensor Fe and the other body by F . If the two di2erent bodies generate the same stress 6eld, then F and Fe are called equivalent deformation gradient tensors. Material softening phenomena for unloading functions can be understood by introducing the three di2erent con6gurations and the respective deformation gradient tensors shown in Fig. 9. It is assumed that the reference con6guration is a state before material softening shown as 1 in Fig. 9 and the intermediate con6guration is a state after material softening shown as 2 in Fig. 9. The transition from the reference con6guration to the intermediate con6guration takes place without changing mechanical properties. Stress is calculated from F associated with the transition from the intermediate to the current con6guration. Similar methods are used in large deformation plasticity theory. From the con6gurations shown in Fig. 9, the deformation gradient tensor Fe can be de6ned as follows: Fe = F Fh :

(16)

It is necessary to relate Fe and F to determine material properties at di2erent states. This relationship is established for large deformation theory with concepts similar to those used with the in6nitesimal elastic functions in Eq. (15). In large deformation theory, suppose that the right Cauchy–Green metric tensor C is used to represent the stress function. The metric tensors, Ce and C , can be expressed through the spectral theorem: Ce = (1e )2 N11 ⊗ N11 + (2e )2 N21 ⊗ N21 + (3e )2 N31 ⊗ N31 ;

(17)

C = (1 )2 N1 ⊗ N1 + (2 )2 N2 ⊗ N2 + (3 )2 N3 ⊗ N3 :

(18)

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2297

Since, in Eq. (15), the material softening occurs in the principal strain direction and since the direction does not change after material softening, it is assumed that the same physical hypothesis can be applied to material softening in large deformation. In this case, the following relations for principal stretches and directions between C and Ce are assumed, i = (ie )h

(19a)

Ni1 = Ni :

(19b)

and Therefore, Eqs. (19a) and (19b) lead to the metric tensor relation, C = (Ce )h where Ce is raised the power h (the softening parameter). By di2erentiating Eq. (16) with respect to time, the relationship between velocity gradients L and Le can be obtained: L = Le − GL(4) Le ;

(20)

for which GL(4) is de6ned with constant h as follows (see Appendix B for details):   3   je ie j − i je ie (4) (1 − h)niiii + (nijji + nijij ) +  (njiji + njiij ) : GL = (ie )2 − (je )2 i j i=1 ij=12;13;23

(21)

The rate of deformation tensor D is the symmetric part of L and is obtained as follows: (4) e D  = D e − GD D

or

(4) e D = Tee D;

(22)

(4) where Tee is a transformation tensor.

4.3. A thermodynamics approach to hyperelastic modeling with softening Suppose that a material point in the intermediate con6guration is characterized by the Helmholtz free energy , Cauchy stress tensor , softening parameter h, and a right Cauchy–Green metric tensor C : = ˆ (C ; h);

(23a)

ˆ  ; h):  = (C

(23b)

Applying the energy conservation equation and Legendre transformation on the 2nd law of thermodynamics leads to the following equation in the absence of temperature [28]: @ˆ ˙ @ ˆ T  F ·L − h ¿ 0; (24)  @C @h where  is density. With velocity gradients obtained by di2erentiating Eq. (16), the 6rst term in Eq. (24) is rewritten as ˆ · Le − 2F

ˆ · Le = ˆ · L + ˆ · F Lh F−1 :

(25)

2298

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

With Eq. (25) and Lh , the inequality in Eq. (24) can be rewritten as follows (see Appendix C for details):     3  ˆ ˆ 1 ˆ @ @ − 2F + 2 ˆ · (26) FT · L −  ln(i )ni ⊗ ni h˙ ¿ 0: h @C @h h i=1 To ensure that the inequality in Eq. (26) is satis6ed for all mechanical processes, the following relationship must hold: ˆ = 2hF

@ ˆ T F : @C

(27)

˙ following relationship Since the terms in the second parenthesis in Eq. (27) are independent of h, holds: 3

 @ˆ 1  ln(i )ni ⊗ ni : = − 2 ˆ · @h h i=1

(28)

For convenience, the Helmholtz free energy ˆ , de6ned as energy per unit mass, can be converted to energy per unit reference volume or strain energy W  : W  = R2 ˆ ;

(29)

where R2 is the density in the equivalent intermediate con6guration 2 . If the strain energy function W  in Eq. (29) is de6ned as Wˆ h then the Cauchy stress can be written as follows: W =

(30)

1  @Wˆ T F F : (31) J  @C Therefore, when there is material softening, Eq. (31) can be used for calculating the stress for a given strain energy function Wˆ . ˆ = 2

5. Visco-hyper-elasto-plastic model development In this section, based on the above theory for modeling visco-hyperelastic materials with softening, the procedure for developing the visco-hyper-elasto-plastic model is presented. The theoretical development starts with a rate independent plastic model and then is extended to address rate e2ects during the loading and unloading process. 5.1. Assumptions for the rate independent plastic model and its material behavior As was previously mentioned, with the assumption that the lateral deformation is insigni6cant compared to the deformation in the thickness direction, a 1-D plasticity model is plausible. Additional

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

Compressive stress

T Trial stress

4

Y Yield stress

2299

σ >y

σ=y

1

3

Shi so that σ = y Shift

ε

p

ε Compressive strain

m

2

Fig. 10. Softening of the elastic unloading response.

assumptions include: • yielding only takes under compressive stress in the thickness direction, • the yield (Tow) stress y∞ is a function of mechanical strain, and • the long-term elastic response (or unloading) ∞ is a function of elastic strain. In 1-D plasticity, the yield function F is simply de6ned as follows: F = e3 · ∞ e3 − e3 · y∞ e3 = 0

or

e3 · ∞ e3 = e3 · y∞ e3 ;

(32)

where e3 is a unit vector in the current thickness direction and is based on the deformation gradient tensor. Fig. 10 is a schematic of the mat material’s 1-D plasticity response during loading and unloading. The thin solid line represents the plastic (yield stress) response during loading. Experimentally it has been observed that the softening parameter h decreases as strain increases. Therefore, in Fig. 10, as the compressive strain increases, the elastic unloading response (thin dotted line) begins to soften and becomes the thick dotted line. For a given strain $m , the trial stress (Point 4) is obtained from the softened elastic curve (thick dotted line), and the yield stress (Point 1) is determined from the yield stress curve. The trial stress is a stress predictor used in the numerical method for the plastic strain calculation. Its value is evaluated from the elastic unloading curve for the given current deformation based on the between plastic strain at the previous step and is subsequently compared with the yield condition. Until the trial stress satis6es the yield condition, the plastic strain is modi6ed and updated from the previous step. Since the magnitude of the trial stress is higher than the magnitude of the yield stress, the material has yielded and the trial stress must be modi6ed to satisfy the yield condition. Satisfying the yield condition requires shifting the thick dotted line to the thick solid line by the plastic strain magnitude. If the loading is removed from Point 1, the unloading will follow the thicksolid line. Point 2 is the permanent plastic strain when the load is completely removed.

2300

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

5.2. The plastic deformation gradient tensor and yield conditions In large deformation theory, the deformation gradient tensor can be multiplicatively decomposed into elastic and plastic parts as follows [25,26]: Fm = Fe Fp ; Fe

(33) Fp

is due to elastic, and is due to plastic deformation. For analytical convenience, the where material rigid body rotation is assumed to come only from the elastic deformation gradient tensor, and plastic deformation does not contribute any material rotation during the deformation. Then, the plastic deformation gradient tensor can be simpli6ed to Up , a right stretch tensor: Fp = R p U p = U p :

(34)

Since plastic deformation is assumed to take place only in the material thickness direction, the right stretch tensor contains only one component and can be written as: U p = E 1 ⊗ E 1 + E 2 ⊗ E 2 +  p E3 ⊗ E 3 ;

(35) p

is a plastic stretch in the thickness where E1 , E2 and E3 are the material coordinate directions and direction E3 . During the plastic straining, p begins to decrease and the plastic deformation gradient tensor will satisfy the yield condition (Eq. (32)). Therefore, the plastic loading and unloading conditions are de6ned through p as follows: Loading condition ˙p ¡ 0 and p ¡ 1: Unloading condition ˙p = 0: 5.3. De:nitions of loading and unloading functions Fig. 9 shows the four di2erent con6gurations that are needed to de6ne the loading and unloading with material softening and plastic strain. For a given deformation Fm , due to plastic deformation, the material reference con6guration 0 6rst changes to the intermediate con6guration 1 . When the material softens, then the intermediate con6guration 1 changes to another con6guration 2 . Therefore, the unloading stress can be calculated from F . The elastic unloading stress ∞ is de6ned between the current con6guration and the intermediate con6guration 2 . Cauchy stress is obtained from the unloading potential Wˆ U as follows: 3 1   @Wˆ U (1 ; 2 ; 3 ) ∞ =   (ni ⊗ ni ): (36) J i=1 i @i Recall that any point on the loading stress curve is regarded as the yield stress y∞ , which is a function of mechanical strain. Therefore the yield stress is derived from the loading potential WL for a given deformation gradient tensor Fm as follows: 3 1  m @WL (1m ; 2m ; 3m ) y∞ = m  (ni ⊗ ni ): (37) J i=1 i @im

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2301

5.4. Elastic rate of deformation tensors The relation between the total mechanical and elastic rate of deformation tensors (Lm and Le ) can be obtained by di2erentiating Eq. (33) with respect to time: Le = Lm −

˙p ML ; p

(38)

where ML =

3  3  e i

i=1 j=1

je

(Ni1 · E3 )(Nj1 · E3 )ni ⊗ nj :

(39)

The rate of deformation gradient tensor De , which is the symmetric part of Le , can be obtained as follows: De = Dm −

˙p MD ; p

(40)

where MD is the symmetric part of ML . The unknown rate of plastic stretch ˙p can be obtained by taking a time derivative of the yield function. This time derivative results in the following equation: ·

·

e3 · ∞ e3 − e3 · y∞ e3 = B1 ˙p + B2 = 0;

(41)

where the constants B1 and B2 are given by (see Appendix D for details):  1 1 (4) (4) B1 = p (Tee MD · I)(e3 · ∞ e3 ) −  e3 · P(4) Tee M D e3  J  (4) − 2e3 · (ML − GL ML )∞ e3 ; (4) m B2 = 4e3 · (∞ − y∞ )Dm e3 − (Tee D · I)(e3 · ∞ e3 ) +

+ (Dm · I)(e3 · y∞ e3 ) −

(42)

1 (4) m e3 · P(4) Tee D e3 J

1 e3 · PY(4) Dm e3 − 2e3 · GL(4) Dm ∞ e3 : Jm

(43)

(4) e F e F e F e D(4) (4) is de6ned in the = FiP In Eqs. (42) and (43), P(4) is de6ned as Pijkm jQ kR mS PQRS and D e m F m F m F m D(4) and D(4) is de6ned in the = FiP equation S˙ = D(4) E˙ . Also, PY(4) is de6ned as PY(4) jQ kR mS YPQRS Y ijkm (4) m ˙ equation Y = DY E . Y is the second Piola–Kirchho2 stress for the yield stress. Therefore, ˙p can be substituted into Eq. (40) and De is obtained with respect to Dm :   1 e (4) (4) m D = I − B2 Dm = Tem D ; (44) B1  p (4) where the fourth order transformation tensor Tem is introduced for convenience.

2302

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

5.5. Including rate e;ects for loading processes From the experimental results, it is seen that there are two kinds of rate e2ects in intumescent mat material. One is a rate e2ect during the loading process and the other is a rate e2ect during the unloading process. To simplify the formulation of the rate dependency for loading processes, the stress response during the loading process is assumed to be composed of a static yield stress due to plastic strain and a rate dependent stress due to elastic strain. Thus, two rate e2ects are incorporated: one associated with elastic unloading and a di2erent one associated with the elastic portion of the loading. The viscous property constants associated with the unloading function are determined by whether or not plastic deformation has occurred. The existence of plastic deformation is determined by evaluating the static yield condition in Eq. (32). The viscoelastic material properties for the unloading function obtained from the loading process are used when there is plastic deformation. When there is no plastic deformation, then the viscoelastic material properties are switched to properties obtained from the unloading function. The 6nal stress is obtained by adding the rate dependent stress to the static stress. Considering that the mat material has virtually no tensile strength and is used almost exclusively in compression, it is assumed that no rate e2ects are present if the static stress in the thickness direction becomes slightly positive. ◦ Referring to Fig. 9, the Jaumman rate of Cauchy stress  from the 2 to con6gurations can be written as follows:  = E(4) De ; ◦

(45)

where De is the rate of deformation gradient tensor and E(4) is the spatial elasticity tensor. It is necessary to convert the spatial elasticity tensor de6ned for the to 2 con6gurations to one de6ned for the to 1 con6gurations. By using the transformation tensor de6ned in Eq. (22), Eq. (45) can be written in terms of De : (4) e  = E(4) Tee D: ◦

(46)

If the spatial elasticity tensor is de6ned for the and 0 con6gurations, the following relation is obtained through Eq. (44): (4) (4) m  = E(4) Tee Tem D : ◦

(47)

For the proposed visco-elasticity model, (4) m  = Etotal D ◦

and

(4) (4) (4) Eˆ (4) Tee Tem : total = vis E

(48)

The viscous term vis is de6ned in Eq. (10). If there is no plastic strain, then ˙p is zero, and in Eq. (4) (48), Tem becomes the Identity tensor. 5.6. Numerical procedures for the plastic stretch ratio and the stress update When the deformation gradient tensor is determined at each time increment, it is necessary to decompose the strain into elastic strain and plastic strain by satisfying the yield function. Suppose

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315 n

that the plastic stretch ratio p has been obtained after a time step n and that Sp at time step n + 1. Then, to satisfy the yield criterion at time step n + 1: n+1

n

F(p + Sp ) = 0

n+1

2303

is the increment (49)

and p

n+1

n

n+1

= p + Sp :

(50)

Performing a Taylor series expansion of Eq. (49) and neglecting higher order terms results in the following equation: n −F(p ) n+1 ; (51) Sp = (@F(pn ))=@p which is not a linear equation and will require linearization and iteration. The elastic trial stress ˆ∞ at time step n+1 can be calculated from the elastic deformation gradient tensor by using the previous plastic deformation gradient tensor at time step n. The procedure is summarized as follows: n+1 n+1 n (i) assume the elastic deformation gradient tensor to be Fˆ e =Fm Fp−1 with p at the previous time step n, (ii) calculate the elastic trial stress ˆ∞ by Eq. (36), (iii) calculate the yield stress y∞ by Eq. (37), and (iv) check the yield condition by Eq. (32).

If yield occurs, (i) (ii) (iii) (iv)

update Fp by the Newton–Raphson method until the yield condition Eq. (32) is satis6ed, calculate ∞ by Eq. (36), p p p calculate vis by using Eq. (8) with a = a and  =  , and p write the stress as  = ∞ + vis .

If no yield occurs, e by using Eq. (8) with a = ae and  = e and (i) calculate vis   e . (ii) write the stress as  = ˆ∞ + vis

6. Model parameters and comparisons with experiments The Ogden strain energy function employed by [27] is used as a base function for both loading and unloading. The loading function, WL in Eq. (37) and the unloading function, Wˆ U in Eq. (36) are therefore de6ned as follows:

2)L L 1 − L L L L WL (1 ; 2 ; 3 ) = 2 1 + 2 + 3 − 3 + (J − 1) ; (52) L L

2)U U 1 U U − U U ˆ (J − 1) ; (53) WU (1 ; 2 ; 3 ) = 2 1 + 2 + 3 − 3 + U U

2304

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315 Displacement loading with specified speed Rigid plate Contact between mat and upper die

Z X

Fig. 11. A schematic diagram of the cyclic compression simulation.

where J is the Jacobean of the deformation gradient tensor,  are the principal stretches, and , ), and are material parameters to be determined from experiments. The constitutive model is coded in the ABAQUSTM user subroutine called UMAT. A single 8 node solid element is used for comparison with the monotonic and cyclic compression experiments. The size of the element is 10 mm × 10 mm × 10 mm. The bottom surface is 6xed in the normal direction, lateral surfaces are free, and compressive (normal) displacements are applied to the top surface nodes of the element. Di2erent compression speeds are used to compare with experimental results. During cyclic displacement control compression, if the mat response is very slow compared to the machine response and/or if there is plastic deformation during the loading process, the specimen may loose contact with the upper loading rod during unloading. This e2ect is simulated by employing a contact surface between the loading rod and top surface of the mat element. Fig. 11 shows a schematic diagram for the cyclic compression simulation. Cyclic displacement is applied to the rigid plate in the Z direction. 6.1. Determination of material parameters In the proposed visco-hyper-elasto-plastic model, two di2erent loading and unloading functions are de6ned. As was explained in Section 5.7, the loading function has only a static term and the unloading function has a static and a transient term. There are two material parameters in the static term and two material parameters to de6ne the transient term. The loading function material parameters are determined from the static stress versus strain curve shown in Fig. 3. The material parameters for the loading function (L and )L ) are obtained from non-linear curve 6tting of the loading curve. For the unloading function, as shown in Fig. 6, each static curve has a di2erent strain at zero stress and a di2erent slope. From the idea of material similarity for each unloading function, all unloading curves in Fig. 6 can be represented with one master curve by changing the softening parameter h and by shifting the strain at zero stress to zero strain. Fig. 12 shows the master curve and the unloading curves from di2erent strains. As shown in Table 2, the softening parameter h decreases with increasing strain. From the master curve shown in Fig. 12, the static material parameter U

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2305

2

Compressive stress (MPa)

Master curve Unloading curve at 1.1MPa (h=1.0) Unloading curve at 1.6MPa (h=0.98) Unloading curve at 2.22MPa (h=0.92)

1.5

1

0.5

0 -0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Compressive true strain

Fig. 12. The master curve and unloading test data at di2erent compressive stresses with material softening removed. Table 2 The softening parameter h and the zero strain point for unloading curves associated with three di2erent maximum applied mechanical strains Mechanical strain Softening parameter h Zero strain

0.8 1 0.532

0.9 0.98 0.6

0.99 0.92 0.65

and )U are determined for the unloading function. In the elastic region, transient material parameters ae and e are obtained from both the relaxation curve at the end of the 1st loading and the rate dependent 2nd loading curves shown in Fig. 8. For the plastic deformation regime, transient material p p parameters a and  are obtained from the transient stress–strain curves in Fig. 6 with the static loading function. 6.2. Simulation results After obtaining the necessary material parameters, simulations are performed to compare with several cycles of the compression experiment results. Fig. 13 shows a comparison of the rate e2ects during the 1st loading process up to 60% compression. It is seen that in all strain ranges except the small strain region, there is good agreement between the experiments and simulations at each compression speed. In the small strain region, the simulations show higher stress than the experimental results, especially as ram speed increases. It appears that the mat material may have di2erent rate e2ects at di2erent compressive strain levels.

2306

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

Compressive stress (MPa)

8

Simulation (static) Analysis(static) Simulation (0.005mm/s) Analysis(0.005mm/s) Analysis(0.4mm/s) Simulation (0.4mm/s) Simulation (1.6mm/s) Analysis(1.6mm/s) Experiment Test(static)(static) Experiment Test(0.005mm/s) (0.005mm/s) Experiment Test(0.4mm/s) (0.4mm/s) Experiment (1.6mm/s) Test(1.6mm/s)

6

4

2

0 0

0.2

0.4 0.6 0.8 Compressive true strain

1

Fig. 13. Comparison of simulation results with experimental data at di2erent compression speeds (60% compression).

Compressive stress (MPa)

2.5

2 Simulation Analysis Experiment Test

1.5

1

0.5

0 0

0.2

0.4 0.6 0.8 Compressive true strain

1

Fig. 14. Comparison of the simulation and experimental static cyclic loading and unloading curves.

Fig. 14 shows a comparison between the static cyclic compression experiment presented in Fig. 6 and the simulation. The simulations show good agreement for the di2erent unloading curves. Fig. 15 compares two cycles from the experimental and simulation results when 60% compression is applied with 1:6 mm=s ram speed. The experimental loading and unloading curves for

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2307

Compressive stress (MPa)

8 1st cycle simulation 2nd cycle simulation 1st 1st cycle cycleexperiment test 2nd 2nd cycle cycleexperiment test

6

4

2

0 0

0.2

0.4

0.6

0.8

1

Compressive true strain

Fig. 15. Comparison of simulations with the 60% cyclic, displacement control compression experiment for 2 cycles (v = 1:6 mm=s). 3.5 1st cycle simulation 2nd cycle simulation 1st cycle experiment 2nd cycletest experiment 2nd cycle

Compressive stress (MPa)

3 2.5 2 1.5 1 0.5 0 0

0.1

0.2 0.3 0.4 0.5 0.6 Compressive true strain

0.7

0.8

Fig. 16. Comparison of two cycles of the 50% cyclic compression simulations and experiments for 2 cycles (v=1:6 mm=s).

the 1st and 2nd cycle are in good agreement with the simulation results. Fig. 16 shows a simulation result for 50% compression with 1:6 mm=s ram speed. During the 1st loading, stress from the simulation is higher than the experimentally determined stress. However, since the operating strain range for this material is generally higher than 50%, the results are considered to be satisfactory for the 1st and 2nd cycles.

2308

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

Compressive stress (MPa)

8 Experiment Simulation 6

4

2

0 0

50

100 150 200 Time (seconds)

250

300

Fig. 17. Comparison of two cycles of stress versus time from the 60% compression simulation and experiments (v = 1:6 mm=s).

Fig. 17 shows two cycles of stress variation with respect to time for 60% compression with 1:6 mm=s ram speed. The 6gure shows the stress jump during loading and the stress relaxation during unloading. The peak value at each cycle and the stress relaxation with respect to time are very well represented by the simulation. However, stress relaxation results from the 2nd cycle show some discrepancies between experiments and simulations. Fig. 18 shows comparisons of peak stresses versus number of cycles from compression experiments with 50% and 60% strain applied under 1:6 mm=s ram speeds. For both cases, good agreement is seen through the 2nd cycle. Simulation results show a constant stress level after the 3rd cycle, but experimental results show continued stress decay with 60% compression data decaying substantially more than for 50% compression. This may be due to cyclic material damage after plastic deformation. Therefore, incorporation of cyclic material damage is believed to be necessary for better simulation of longer term cycling. 7. Concluding remarks This paper presented detailed development of a constitutive model for intumescent mat material. After starting with hyperelasticity theory, the model incorporated rate e2ects, plasticity, and material softening during cyclic loading processes. The model is limited to the room temperature response of materials that have large lateral to thickness dimension ratios and that are deforming primarily in the thickness direction. The results of several experiments on intumescent mat were reviewed. A new material softening model was proposed from the assumption of material similarity. The softening parameter was determined as a function of the strain in the thickness direction. For the 1-D plastic model, a new plastic deformation gradient tensor and yield criteria in the thickness direction were

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2309

Peak compressive stress (MPa)

8 Experiment (60%) Experiment (50%) Simulation (60%) Simulation (50%)

6

4

2

0 2

4 6 Number of cycles

8

10

Fig. 18. Comparison of simulated and experimental peak stress versus number of cycles (v = 1:6 mm=s).

proposed. Also, the plastic yield stress function was established on the basis of strain. Two di2erent rate e2ects in the loading and unloading regions were implemented by employing di2erent viscosity parameters for the plastic and elastic unloading regions. The complete spatial elasticity tensor in the Eulerian frame was derived for implicit 6nite element coding and the Newton–Raphson method was applied to obtain plastic strain. The theory was coded in ABAQUSTM implicit using the user subroutine UMAT. Experimental results and simulations are compared to each other for several experimental cases. For up to two cycles, satisfactory agreement was achieved with the new model. Further study is still required to accurately predict loading and unloading processes for many cycles. It is believed that cyclic damage e2ects after the 2nd loading and di2erent viscous e2ects at di2erent compressive strain levels should be considered in the future. In addition, future versions of the model should attempt to incorporate elevated temperature material responses. Although the new constitutive model developed in this paper has been applied to modeling intumescent mat materials, the same modeling approach could be applied to rubber materials to account for Mullin’s e2ect (i.e., material damage due to cyclic loading). In addition, the softening model can be applied to any anisotropic material to account for material property changes due to temperature. Therefore, numerous applications of the new constitutive model are envisioned. Acknowledgements This work was 6nancially supported by ArvinMeritor Co. We gratefully acknowledge Bill Taylor, Joe Fuehne, and Han-Jun Kim from ArvinMeritor for supporting this project and for their valuable

2310

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

comments. We would also like to thank John TenEyck from Unifrax Corp. for supplying specimens and for insightful discussions. Appendix A. The Jaumann rate of Kirchho6 stress The 2nd Piola Kirchho2 stress(PK2) can be related to the true stress by  = J e  = Fe SFe : T

(A.1) ◦

The Jaumann rate of Kirchho2 stress  can be obtained by di2erentiating the Kirchho2 stress with respect to time and is described with the rate of stretch tensor De : ◦

˙ e + De + De :  = Fe SF T

(A.2)

By using the following relations: T E˙ e = Fe De Fe ;

(A.3)

J˙e = J e div(C) = J e (De · I);

(A.4)

when S˙ is represented as S˙ = D(4) E˙ e , the Jaumann rate of Cauchy stress is obtained as ◦

 = −(De · I) +

1 e (4) eT e e eT (F D (F D F )F ) + De + De : Je

(A.5)

In tensor component form, Eq. (A.5) can be rewritten as follows. ◦

e ij = −Dmm ij +

1 e e e e (4) e (F F F F D De ) + ik Dkj + Dike kj : J e iP jQ kR mS PQRS km

(A.6)

New 4th order tensors H(4) and P(4) are de6ned as follows: (4) = 12 (ia bj + ia bj + ib aj + ib aj ) Hijab

(A.7)

(4) (4) e e e e = FiP FjQ FkR FmS DPQRS : Pijkm

(A.8) ◦

The Eulerian tangential sti2ness tensor E(4) de6ned in the relation  = E(4) De can be written as E(4) = − ⊗ I +

1 (4) P + H(4) : Je

(A.9)

Appendix B. Intermediate derivations for Eqs. (20)–(21) Time di2erentiation of Eq. (16) is written as follows: L = Le − F Lh F−1 ;

(B.1)

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2311

where Fh , F and Fe are de6ned as: e e e Fh = 1 N11 ⊗ N11 + 2 N21 ⊗ N21 + 3 N31 ⊗ N31 ; 1 2 3

(B.2)

F = 1 n1 ⊗ N11 + 2 n2 ⊗ N21 + 3 n3 ⊗ N31 ;

(B.3)

Fe = 1e n1 ⊗ N11 + 2e n2 ⊗ N21 + 3e n3 ⊗ N31 :

(B.4)

h

Time di2erentiation of F in Eq. (B.2) is rewritten as follows:  ·   e 3 e e      i Ni1 ⊗ Ni1 + i N˙ 1i ⊗ Ni1 + i Ni1 ⊗ N˙ 1i  : F˙ h = i i i i=1

(B.5)

˙  ) and N˙ 1 are obtained as In Eq. (B.5), (ie = i i ·  e i (1 − h) 1 e = (Ni ⊗ Ni1 ) · E˙ e − i (ln ie )h˙ e   i i i i 3  2 N˙ 1i = (N1 ⊗ Nj1 )E˙ e Ni1 : e 2 (i ) − (je )2 j

(B.6) (B.7)

j=1=i

Lh is obtained as follows from Lh = F˙ h Fh : −1

h

L =

3  (1 − h)

(ie )2

i=1



ie j − i je Niiii E + (ie )2 − (je )2 ij=12;13;23 ˙e



(Njiji + Njiij ) (Nijji + Nijij ) + ie j i je

3  ˙ i1 ⊗ Ni1 ; − (ln ie )hN

 E˙ e

(B.8)

i=1

where E˙ e = Fe De Fe . Therefore, F Lh F−1 in Eq. (B.1) is can be written as T

 h −1

FL F

=

3 

e

(1 − h)niiii D −

i=1

3  i=1



ie j − i je + (ie )2 − (je )2 ij=12;13;23

˙ i ⊗ ni ) (ln ie )h(n 

 je ie (nijji + nijij ) +  (njiji + njiij ) De : i j

(B.9)

Appendix C. Intermediate derivations for Eqs. (24)--(26) From Eq. (B.1), Le is written as follows: Le = L + F Lh F−1 :

(C.1)

2312

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

 Eq. (B.6) can be rewritten with respect to E˙ :



·

ie i



 =



1 −1 h

1 ie ie  1 1 ˙ ˙ (N ⊗ N ) · E − (ln ie )h; i (i )3 i h i

(C.2)

 where E˙ = FT D F . Therefore, from Eq. (C.1), Lh in Eq. (B.8) can be presented as    3    ie j − i je + N ) + N ) (N (N 1 1  jiji jiij ijji ijij h ˙ E˙ e L = + −1 e )2 − (e )2 e   )2 Niiii E +  e h ( (   i i j i j i j i=1 ij=12;13;23

−

3  1 ˙ i1 ⊗ Ni1 ; (ln i )hN 2 h i=1

(C.3)

F Lh F−1 in Eq. (C.1) is obtained as follows:  h −1

FL F

=

3   1



3 1  ˙ i ⊗ ni ) (ln ie )h(n 2 h h i=1 i=1    je ie j − i je ie e + e )2 − (e )2  (nijji + nijij ) +  (njiji + njiij ) D : (  i j i j ij=12;13;23

− 1 niiii L −

(C.4)

Also,  h −1

ˆ · F L F

 =

 3  1 1 ˙ − 1 ˆ · L − 2 ˆ · ln(i )ni ⊗ ni h: h h i=1

(C.5)

In this way Eq. (26) is obtained from Eqs. (C.5) and (24).

Appendix D. Intermediate derivations for Eqs. (41)--(43) The 6rst half of 6rst term in Eq. (41) can be written as ·

e3 · ∞ e3 = e˙3 · ∞ e3 + e3 · ˙∞ e3 + e3 · ∞ e˙3 :

(D.1)

When E3 and e3 are related by the equation Fm E3 = re3 , where r is stretch ratio, e˙3 is obtained as follows: e˙3 = Lm e3 − [(e3 ⊗ e3 ) · Dm ]e3 :

(D.2)

When ∞ is de6ned as ∞ = 1=J  F S∞ FT , ˙∞ can be written as ˙∞ = −(D · I)∞ + L ∞ + ∞ LT +

1 (4)  P D: J

(D.3)

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2313

Therefore, the 1st, 2nd and 3rd terms in Eq. (D.1) are given by e˙3 · ∞ e3 = Lm e3 · ∞ e3 − [(e3 ⊗ e3 ) · Dm ](e3 · ∞ e3 );

(D.4)

e3 · ∞ e˙3 = e3 · ∞ Lm e3 − [(e3 ⊗ e3 ) · Dm ](e3 · ∞ e3 );

(D.5)

and (4) e e3 · ˙∞ e3 = −(Tee D · I)(e3 · ∞ e3 ) + 2e3 · L ∞ e3 +

1 (4) e e3 · P(4) Tee D e3 ; J

(D.6)

and therefore Eq. (D.1) is written as follows: ·

(4) e e3 · ∞ e3 = 2e3 · ∞ Lm e3 − 2[(e3 ⊗ e3 ) · Dm ](e3 · ∞ e3 ) − (Tee D · I)(e3 · ∞ e3 )

+ 2e3 · L ∞ e3 +

1 (4) e e3 · P(4) Tee D e3 : J

(D.7)

In the same way, the 2nd term in Eq. (41) is written as follows: ·

e3 · y∞ e3 = 2e3 · y∞ Lm e3 − 2[(e3 ⊗ e3 ) · Dm ](e3 · y∞ e3 ) − (Dm · I)(e3 · y∞ e3 ) + 2e3 · Lm y∞ e3 +

1 e3 · PY(4) Dm e3 : Jm

(D.8)

(4) e F e F e F e D(4) (4) is de6ned in the equation In Eq. (D.7), P(4) is de6ned as Pijkm = FiP jQ kR mS PQRS and D e m F m F m F m D(4) and D(4) is de6ned in the equation = FiP S˙ = D(4) E˙ . Also, PY(4) is de6ned as PY(4) jQ kR mS YPQRS Y ijkm (4) m ˙ Y=D E . Y is the 2nd Piola Kirchho2 stress for the yield stress. Next, with the following relations: Y

Le = Lm −

˙ ML p

D e = Dm −

˙ MD p

and and

L = Le − GL(4) Le ;

(D.9a)

(4) e D = Tee D;

(D.9b)

Eq. (41) can be expressed as follows: ·

e3 · (∞ − y∞ )e3   1 ˙p (4) (4) = p (Tee MD · I)(e3 · ∞ e3 ) −  e3 · P(4) Tee MD e3 − 2e3 · (ML − GL(4) ML )∞ e3  J (4) m + 4e3 · (∞ − y∞ )Dm e3 − (Tee D · I)(e3 · ∞ e3 ) +

+ (Dm · I)(e3 · y∞ e3 ) −

1 (4) m e3 · P(4) Tee D e3 J

1 e3 · PY(4) Dm e3 − 2e3 · GL(4) Dm ∞ e3 : m J

Therefore the constants B1 and B2 are obtained as in Eqs. (43) and (44).

(D.10)

2314

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

References [1] Fuehne J, Taylor W, Kim J, Lee JK. Characterization of catalytic converter mat material for predictive analysis. SAE International Congress and Exposition, Detroit, MI, 2000. SAE paper #2000-01-0219. [2] Taylor W, Fuehne J, Lyon R, Kim J, Lee JK. FEA simulation and experimental validation of catalytic converter structural integrity. SAE International Congress and Exposition, Detroit, MI, 2000. SAE paper #2000-01-0219. [3] Chen DKS. Integrity of an automotive catalytic converter during the assembly process. ASME Winter Annual Meeting, Dallas, TX, 1990. ASME paper #90-WA/DE-6AC. [4] Muju S, Sager R, Snider B. Catalytic converter durability analysis. IMECE’97, ASME International Congress and Exposition, Atlanta, GA, 1996. [5] Bardenhagen SG, Stout MG, Gray GT. Three-dimensional, 6nite deformation, viscoplastic constitutive models for polymeric materials. Mechanics of Materials 1997;25:235–53. [6] Bonet J, Burton AJ. A simple orthotropic, transversely isotropic hyperelastic constitutive equation for large strain computations. Computer Methods in Applied Mechanics and Engineering 1998;162:151–64. [7] Govindjee S, Simo JC. Mullins’ e2ect and the strain amplitude dependence of the storage modulus. International Journal of Solids and Structures 1992;29:1737–51. [8] Hoger A. A second order constitutive theory for hyperelastic materials. International Journal of Solids and Structures 1999;36:847–68. [9] Moran B, Ortiz M, Shih CF. Formulation of implicit 6nite element methods for multiplicative 6nite deformation plasticity. International Journal for Numerical Methods in Engineering 1990;29:483–514. [10] Weiss JA, Maker BN, Govindjee S. Finite element implementation of incompressible, transversely isotropic hyperelasticity. Computer Methods in Applied Mechanics and Engineering 1996;135:107–28. [11] Simo JC. On a fully three-dimensional 6nite-strain viscoelastic damage model: formulation and computational aspects. Computer Methods in Applied Mechanics and Engineering 1987;60:153–73. [12] Papoulia KD, Kelly JM. Visco-hyperelastic model for 6lled rubbers used in vibration isolation. Journal of Engineering Materials-T ASME 1997;119:292–7. [13] Zhang J, Lin Z, Wong A, Kikuchi N, Li VC, Yee AF, Nusholtz GS. Constitutive modeling and material characterization of polymeric foams. Journal of Engineering Materials-T ASME 1997;119:284–91. [14] Gibson LJ, Ashby MF. Cellular solids: structure & properties. New York: Cambridge University Press, 1997. [15] Puso MA, Govindjee S. A phenomenological constitutive model for polymeric foam. In: Mechanics of plastics and plastic composites, MD-v68/AMD-v215, 1995. ASME International Mechanical Engineering Congress and Exposition, San Francisco, CA, 1995. p. 159 –76. [16] Sherwood JA, Frost CC. Constitutive modeling and simulation of energy absorbing polyurethane foam under impact loading. Polymer Engineering Sciences 1992;32:1138–46. [17] Kim H, Kim JS, Walter ME, Lee JK. The thermo-mechanical stress-strain response of intumescent mat materials through experiments and strain decomposition. Journal of Strain Analysis 2002;37:187–99. [18] Kim HJ. Experimental investigation of the thermo-mechanical response of intumescent mat material. M.S. thesis, The Ohio State University, 2000. [19] Treloar LRG. Stress-strain data for vulcanized rubber under various types of deformation. Transactions of Faraday Society 1944;40:59–70. [20] Ogden RW. Large deformation isotropic elasticity on the correlation of theory and experiment for incompressible rubberlike solids. Proceedings of Royal Society of London Series A 1972;328:565–84. [21] Simo JC, Taylor RL. Quasi-incompressible 6nite elasticity in principal stretches. Continuum basis and numerical algorithms. Computer Methods in Applied Mechanics and Engineering 1991;86:273–310. [22] Swanson SR. A constitutive model for high elongation elastic materials. Journal of Engineering Material Techniques 1985;107:110–4. [23] Johnson AR, Quigley CJ, Freese CE. A viscohyperelastic 6nite element model for rubber. Computer Methods in Applied Mechanics and Engineering 1995;127:163–80. [24] Koeller RC. Applications of fractional calculus to the theory of viscoelasticity. Journal of Applied Mechanics 1984;51:299–307. [25] Lee EH. Elastic-plastic deformation at 6nite strains. Journal of Applied Mechanics 1969;36:1–6.

J.S. Kim et al. / International Journal of Mechanical Sciences 44 (2002) 2285 – 2315

2315

[26] Lubarda VA, Lee EH. A correct de6nition of elastic and plastic deformation and its computational signi6cance. Journal of Applied Mechanics 1981;48:35–40. [27] ABAQUS theory manual, HKS, 1998. [28] Coleman BD, Gurtin ME. Thermodynamics with internal state variables. Journal of Chemical Physics 1967;47: 597–613.