- Email: [email protected]
Second-order estimate of the macroscopic behavior of periodic hyperelastic composites: theory and experimental validation
Journal of the Mechanics and Physics of Solids 52 (2004) 27 – 49
www.elsevier.com/locate/jmps
Second-order estimate of the macroscopic behavior of periodic hyperelastic composites: theory and experimental validation N. Lahellec∗ , F. Mazerolle, J.C. Michel L.M.A./C.N.R.S., 31, Chemin Joseph Aiguier, Marseille 13402, Cedex 20, France Received 15 August 2002; received in revised form 30 June 2003; accepted 30 June 2003
Abstract This paper deals with some theoretical and experimental aspects of the behavior of periodic hyperelastic composites. We focus here on composites consisting of an elastomeric matrix periodically reinforced by long 6bers. The paper is composed of three parts. The 6rst part deals with the theoretical aspects of compressible behavior. The second-order theory of Ponte Casta˜neda (J. Mech. Phys. Solids 44 (1996) 827) is considered and extended to periodic microstructures. Comparisons with results obtained by the 6nite element method show that the composite behavior predicted by the present model is much more accurate for compressible than for incompressible materials. The second part deals with the extension of the method to incompressible behavior. A mixed formulation (displacement–pressure) is used which improves the accuracy of the estimate given by the model. The third part presents experimental results. The composite tested is made of a rubber matrix reinforced by steel wires. Firstly, the matrix behavior is identi6ed with a tensile test and a shear test carried out on homogeneous samples. Secondly, the composite is tested under shearing. The experimentally measured homogenized stress is then compared with the predictions of the model. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Homogenization; Hyperelasticity; Rubber material; Experiments
1. Introduction The aim of this paper is to predict the behavior of hyperelastic composites using the homogenization procedure. Knowing the hyperelastic behavior of the composite’s constituents, it is proposed to predict the homogenized behavior as analytically as possible. ∗
Corresponding author. Fax: +33-491-1644-81. E-mail address: [email protected] (N. Lahellec).
0022-5096/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0022-5096(03)00104-2
28
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
When the strain energy is a convex function of the deformation gradient, the homogenized behavior is given by the homogenized strain energy. This function is computed by solving the in6mum strain energy of a single unit cell subjected to the macroscopic deformation gradient, see Marcellini (1978). Unfortunately, hyperelastic strain energies are not convex and, therefore, the in6mum problem cannot be studied in a single unit cell and has to be considered in the structure as a whole, see MBuller (1987). These calculations are very long and in the present study the homogenized strain energy is computed in a single unit cell, based on the assumption that no instability is observed in the range of loads studied. The nonlinearity of the behavior is an additional diCculty for hyperelastic materials. The relationship between stress and strain is strongly nonlinear, which complicates the estimation of their macroscopic behavior. Previous authors, Dumontet (1983, 1989), Geymonat et al. (1993) and Pruchnicki (1998) have studied the case of laminates where the problem is one dimensional. In the compressible case, Dumontet (1983) gave an analytical expression for the homogenized strain energy and in the incompressible case, Pruchnicki (1998) shows that the homogenization problem can be reduced to a simple nonlinear problem. Other authors have linearized the initial problem using an incremental method (Abeyaratne and Triantafyllidis, 1984; Devries, 1996) to approach more complex microstructures, and solved the linear problem using the periodic homogenization method. These methods always lead to a 6nite element method (FEM) resolution. The last group of studies dealt with developing analytical methods for estimating the macroscopic behavior of nonlinear composites, see Ponte Casta˜neda and Suquet (1998), Talbot and Willis (1997). Ponte Casta˜neda (1989) studied Voigt and Reuss bounds for hyperelastic materials and adapted the Hashin–Shtrikman procedure to get a nonlinear bound. Ponte Casta˜neda and Tiberio (2000) applied the second-order theory, 6rst proposed by Ponte Casta˜neda (1996), to hyperelastic composites. Here, we propose to extend this second-order method to estimate the macroscopic behavior of periodic and incompressible composites. We model the periodic microstructure using the lower periodic Hashin–Shtrikman bound proposed by Suquet (1990) and Iwakuma and Nemat-Nasser (1983). This bound is given by a simpli6ed computation of the linear problem in the Fourier space. To overcome the incompressibility constraint, we use a mixed formulation, displacement–pressure, before applying the second-order method. This paper concludes with an experimental part. The method developed is compared with experimental results obtained with a composite consisting of a rubber matrix reinforced by steel wires. Firstly, the matrix behavior is identi6ed with an original strain energy function. Secondly, the estimate of the homogenized stress is compared with the experimental data obtained in a shearing experiment. 2. Hyperelastic compressible and periodic composites 2.1. Second-order method P. Ponte Casta˜neda proposed the second-order method for estimating the macroscopic behavior of nonlinear composites, initially for power-law materials Ponte Casta˜neda
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
29
(1996), Ponte Casta˜neda and Suquet (1998) and more recently in the case of hyperelastic isotropic macroscopic behavior Ponte Casta˜neda and Tiberio (2000). When W denotes the strain energy and F the deformation gradient, the homogenized strain energy is given by the variational principle 1 K ˜ W (F) = inf W (X; F) dV ; (1) K |V0 | V F∈Kper (F) 0 K the set of kinematically where V0 is the unit cell in initial con6guration and Kper (F) admissible deformation gradients with the average deformation gradient FK K K u−(F−I)·X K Kper (F)={F | F=I+∇u; F=F; with : = (1=|V0 |) V0 : dV . The proposed method involves two stages:
V0 –periodic; det F¿0}
(2)
1. Linearization of the initial problem by the second-order procedure, which gives a linear but heterogeneous problem. 2. Resolution of this problem using the Hashin–Shtrikman periodic lower bound. The second-order theory consists of replacing W (F) by its second-order Taylor expansion about FK r which is the average of F over phase r. W (F) = W (FK r ) + (F − FK r ) : r + 12 (F − FK r ) : Lr : (F − FK r ); Lr =
r =
9W K r (F ); 9F
92 W K r (F ): 9F2
Using this expansion, Ponte Casta˜neda (1996) established that the second-order expansion of the homogenized strain energy, W˜ , can be written as K = W˜ (F) (3) cr W (FK r ) + 12 (FK − FK r ) : r : r
K The homogenized stress is then given by the derivative of (3) with respect to F: 1 ˜ K 9FK r ˜ F) K = 9W (F) = ( (4) + r : cr [r + (FK − FK r ) : Lr ] : 2 9FK 9FK r The average deformation gradient in each phase is still present in both expressions (3) and (4) and must be calculated. This is done by solving the following variational problem for a linear heterogeneous material: r inf ∇u : (r − Lr : ∇u ) + 12 ∇u : Lr : ∇u : (5) K F∈Kper (F)
The solution of (5) is studied in the next section.
30
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
2.2. Solving the linear heterogeneous problem 2.2.1. Hashin–Shtrikman lower bound In this section, we solve the linear heterogeneous problem (5) using the Hashin– Shtrikman procedure of Willis (1981). This estimate is obtained with a linear homogeneous reference material, having an elastic strain energy de6ned by W 0 = 12 ∇u : L0 : ∇u:
(6)
The Hashin–Shtrikman estimate is given by r
W HS (∇u) = 12 ∇u : L0 : ∇u + ∇u : (Lr − L0 )∇u :
(7)
When L0 is the matrix modulus (the most compliant constituent), this estimate gives a lower bound. The original linear heterogeneous problem is transformed into a linear homogeneous problem de6ned by
1 r ∇u : (r − L0 : ∇u ) + ∇u : L0 : ∇u : inf (8) K 2 F∈Kper (F) 2.2.2. Solving the periodic and homogeneous problem Ponte Casta˜neda and Tiberio (2000) proposed to use the Green operator in an in6nite medium for random composites. The composites considered here are periodic and we use the periodic Green operator, which can be determined by carrying out a Fourier transform of (8) in line with the procedure used in Suquet (1990) or Iwakuma and Nemat-Nasser (1983). In indicial notation, (8) reads pr L0ijkh uk; hj + ij; j =0
with uk; ij =
9 2 uk 9Xi 9Xj
and
r
pr = r − L0 : ∇u :
(9)
The Fourier transform of (9) is carried out with spatial frequencies de6ned by i = 2P=Li in which P is an integer and Li the unit cell dimension in the i direction, with orthogonal symmetry axes. With V0 the unit cell in its reference con6guration, f(X) a V0 -periodic vector valued function, the Fourier transform of f is de6ned as i·X ˆ ˆ f() = f(X)e−i·X ; f(X) = f + f()e and ∈R∗ −{0}
[ 9f(X) 9Xj
ˆ = ij f():
With these relations, (9) becomes pr −L0ijkh h j uk () + ij ij () = 0:
(10)
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
31
Inverting the acoustic tensor · L0 · and diPerentiating with respect to h yields the following expression for u k; h : pr −1 0 u k; h () = −h [Lijkh h j ]ki j ij ():
(11)
pr We now apply relations (10) to u k; h () and ij () which gives the local transformation gradient value pr −i·X i·X uk; h (X) = uk; h − h [L0ijkh h j ]−1 e : (12) ki j ij (X)e ∈R∗ −{0}
∇u given by (12) is the solution of (9). Strain energy expression (3) contains the average displacement gradient in each phase which is obtained by averaging (12) s Psr : pr ∇u = ∇u + (13) r
with Psr the fourth-order tensor associated with the periodic Green operator and de6ned by −i·X Psr r ei·X s : (14) h [L0ijkh h j ]−1 khij = − ki j cr e ∈R∗ −{0}
This operator depends on the microstructure through e−i·X r ei·X s and on the mas terial nonlinearity through L0 . At this stage, a relation between ∇u , Psr and pr has been established. Since these two terms depend on FK r , (13) is an implicit equation for r the unknown ∇u . K Expression (4) for the homogenized stress contains the fourth-order tensor 9FK r =9F, K which is calculated by taking the derivative of Eq. (13) with respect to F: pr 9FK s 9Psr 9FK r 9FK r pr sr 9 P : : (15) =I+ + : : : 9FK 9FK r 9FK 9FK r 9FK r 2.3. Comparison with 4nite element results In this section, the case of a periodic two-phase material composed of a matrix reinforced with long 6bers with a circular section is investigated. In the initial con6guration, the 6bers were assumed to be arranged in a square array. Two 6ber volume fractions are considered: 19 and 14 . In what follows, the unit cells associated with these microstructures are named S9 and S4. Both the 6bers and the matrix are modeled as Mooney–Rivlin hyperelastic material with strain energy (see ABAQUS, 2001): W (F) = C10 (I1 − 3) + C01 (I2 − 3) +
1 (det(F) − 1)2 ; d
(16)
where F =
F ; det(F)1=3
C = F T · F ;
I1 = tr(C );
I2 = 12 (tr(C )2 − tr(C · C )):
32
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Table 1 De6nitions of the studied loadings Name
Macroscopic deformation gradient
Uniaxial tension Incompressible tension Simple shear Compression
FK = e1 ⊗ e1 + e2 ⊗ e2 + e3 ⊗ e3 ; 1 6 6 3 FK = e1 ⊗ e1 + 1 e2 ⊗ e2 + e3 ⊗ e3 ; 1 6 6 3 FK = e1 ⊗ e1 + e2 ⊗ e2 + e3 ⊗ e3 + e1 ⊗ e2 ; 1 6 6 2 FK = !e1 ⊗ e1 + e2 ⊗ e2 + e3 ⊗ e3 ; 0:3 6 ! 6 1
Fig. 1. Typical mesh (4-node bilinear elements), S4 unit cell.
Note that quasi-incompressible behavior can be modeled by increasing the ratio 1=d (penalty method). Note also that in the limiting case of small strains, these materials can be regarded as linear elastic and isotropic with bulk modulus k = 2=d and shear modulus = 2(C10 + C01 ). The contrast between the 6bers and the matrix is de6ned as W 6bers (F) = contrast × W matrix (F): Only in-plane macroscopic strains are considered here, (O; X1 ; X2 ) parallel planes and loading are de6ned in Table 1, and the problem is solved in plane strain. The homogenized strain energy and stress obtained using the second-order procedure are compared with those resulting from 6nite element calculations (ABAQUS, 2001). A typical mesh using four-node bilinear elements is shown in Fig. 1. 2.3.1. Compressible case The compressible case is solved 6rst. The matrix parameters are matrix matrix C10 = C01 = (1=d)matrix = 1 Pa:
Uniaxial tension: The uniaxial tension case corresponds to an average deformation gradient FK given in Table 1. The maximum strain is 200% ( = 3). In both
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
33
Fig. 2. Homogenized stress, ˜ 11 in (a) and ˜ 22 in (b), vs. load factor during simulated uniaxial tension applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively). The contrast between 6bers and matrix was 20.
microstructures, the 6bers are 20 times stiPer than the matrix. The stress estimate given by the second-order procedure (SOE) is compared with that given by the FEM in Fig. 2. With S9 and S4 microstructures, it can be seen that the diPerence between the two estimates is very small. The greatest error in the stress is less than ˜ FEM ˜ FEM 2%, (error = (˜ SOE 11 − 11 )= 11 ). In all the cases studied, the maximum error was reached with the greatest 6ber ratio, 25% volume fraction and with the largest strain. Compression: We simulated the compression test (Table 1) with microstructures S9 and S4, with two contrasts, ! is the load factor: • Contrast equal to 20, see Figs. 3(a) and (b): the greatest error in the estimate of 11 was nearly 1%. With 22 , the error was of the same order, above ! ≈ 0:55 in the S4 case and ! ≈ 0:38 for S9. Below these values, the error increased to 10% with S4. • Contrast equal to 100: Figs. 3(c) and (d) show that the estimates given by the second-order procedure are less satisfactory. The error was still low, above ! ≈ 0:65 with S4 and ! ≈ 0:5 with S9, but increased below these values to 25% with 11 and more than 50% with 22 . Simple shear: We simulated the simple shear test de6ned by FK in Table 1 up to 200% shear strain ( = 2). Microstructures S4 and S9 were taken to have a contrast equal to 20. In Fig. 4, the SOE and the FEM are compared. In the case of S9, the 12 and 21 estimates were accurate, and the greatest error was less than 3%. The estimate was also accurate in the case of S4 and at a shear strain less than 100%. Above this strain value, the error amounted to 15% at the maximum strain ( = 2). Conclusion for the compressible case: The nonlinearity of the behavior of hyperelastic composites makes it diCcult to calculate their homogenized or ePective behavior. This nonlinearity is due to the large strains involved and the resulting changes in the microstructure. The present model takes into account the nonlinearity of the material with the inTuence of the average deformation gradient in the matrix phase, in L0 and
34
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 3. Homogenized stress vs. load factor ! during simulated uniaxial compression applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively). The contrast between 6bers and matrix was 20 in (a) and (b) and 100 in (c) and (d).
Fig. 4. Homogenized stress, ˜ 12 in (a) and ˜ 21 in (b), vs. load factor during simulated simple shear applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively). The contrast between 6bers and matrix was 20.
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
35
Fig. 5. Pre-compressed shear simulation and description of microstructure evolution, with microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively). The contrast between 6bers and matrix was 20.
the various eigenstresses. The simulation, which gave very accurate estimates, uniaxial tension simulation, showed that this model accounts accurately for the nonlinearity when the strain 6eld is not localized. In the simulation, the 6ber interaction zone, namely the area in between two 6bers, is large and does not shrink with the loading factor. In the case of compression, the interactions between 6bers evolve. By increasing the loading factor, the 6bers move closer and the matrix between them is compressed and becomes stiPer, because the volume cannot be compressed to zero (det F 9 0). With small contrasts, this ePect is not very important because the interaction zone is large enough to inTuence the average matrix behavior in FK 2 . On the other hand, increasing the contrast decreases the interaction zone as well as its ePect on the average matrix behavior. In the limiting case of rigid 6bers, in6nite contrast, the composite locks up when the 6bers come into contact. With the microstructures modeled here, max 1 = 2 c1 =, where square unit cells with circular 6ber, this locking occurs when F11 max 1 max 1 (S4) = 0:56; F11 (S9) = 0:38). The model c1 is the 6ber volume fraction (F11 2 max 2 K = c1 , which is smaller proposed here locks up when det F = 0, which occurs at F11 1 max . Simulation of simple shear shows changes in the microstructure shape than F11 which are not very well taken into account by the proposed model. In this case, the microstructure alternates between a con6guration in which the 6bers are aligned and one in which they are arranged crosswise. This tendency is even more pronounced in a pre-compressed shear simulation with the deformation gradient de6ned as FK = e1 ⊗ e1 + 12 e2 ⊗ e2 + e3 ⊗ e3 + e1 ⊗ e2 ;
: 0 → 2:
(17)
During this simulation, every 6ber is lined up with either another 6ber or between two others (see Fig. 5). The macroscopic behavior depends strongly on the microstructure
36
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 6. Homogenized stress ˜ 11 vs. load factor during simulated incompressible tension in (a) and ˜ 21 vs. load factor during simulated simple shear in (b), applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 11% and 25%, respectively), incompressible behavior. The contrast between 6bers and matrix was 20.
changes. In the case of S9, SOE accurately predicts the two 6rst 6bers alignments, whereas in the case of S4, the 6rst alignment is well predicted but then the error becomes very large. 2.3.2. Quasi-incompressible case The aim of this study was to model composites made from incompressible rubber. In order to check the accuracy of the SOE in this case, simulations were performed with the same values for C10 and C01 as in the previous case, increasing the ratio 1=d so as to reach the incompressibility limit. The simulations shown here were carried out with composites de6ned by 6bers 20 times stiPer than the matrix (contrast equal to 20), which are also de6ned by C10 = C01 = 1 Pa and 1=d = 500 Pa. Fig. 6 shows the evolution of the tensile stress and of the shear stress corresponding to incompressible tension and shear as de6ned in Table 1. In both cases, the results of the SOE deviate signi6cantly from the FEM values. Conclusion for quasi-incompressible composites: The above simulations show that as we get closer to the incompressibility limit, the estimate given by the classic second-order formulation is no longer satisfactory. For both microstructures, we varied the 1=d values from 0.1 to 1000 Pa and the accuracy of the estimate was found to decrease sharply at around 500 Pa. This can be explained by the strong nonlinearity introduced by the incompressibility constraint. With materials of this kind, the strain energy can be written as follows: C10 (I1 − 3) + C01 (I2 − 3) if det(F) = 1; (18) W (F) = +∞ if det(F) = 1: This strong nonlinearity of det(F) is not accounted for by a second-order Taylor expansion, which explains these poor results. This was noticed by Ponte Casta˜neda and Tiberio (2000), who established, in the case of rigid inclusions that the straightfor-
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
37
ward use of the second-order approximation led to an approximate incompressibility constraint which is inconsistent with the exact one. 2.4. Conclusion The above simulations (in the compressible and quasi-incompressible cases) show that SOE is an accurate tool for dealing with compressible periodic composites in loading cases where the macroscopic behavior does not depend too strongly on 6ber interactions. In the next section, the second-order method is adapted to model incompressible composites.
3. Modi ed method for incompressible hyperelastic materials 3.1. New formulation It was shown in Section 2.3 that the traditional second-order formulation is not appropriate for dealing with incompressible materials. Ponte Casta˜neda and Tiberio (2000) developed a method of computing the incompressible limit for rigid inclusions, but it cannot be applied here (because of the 6nite contrast). Therefore, an alternative formulation has to be developed. The initial variational principle reads K = W˜ (F)
inf
K F∈Kper; inc (F)
W (F);
(19)
K denotes the set of kinematically admissible deformation gradients where Kper; inc (F) with the average deformation gradient FK and the incompressibility constraint K = {F | F = I + ∇u; F Kper; inc (F) K u − (FK − I) · X V0 –periodic and det(F) = 1}: = F; We replace this constrained in6mum problem by a saddle-point problem, where the incompressibility constraint is relaxed by introducing its associated Lagrange multiplier p, the pressure. K p) ˜ F; L( K = sup
inf
K p∈P(p) K F∈Kper (F)
L(F; p);
L(F; p) = W (F) − p(det(F) − 1); (20)
with P as the set of admissible pressures p. P(p) K = {p | p V0 –periodic; p = p}: K When W is a convex function of the deformation gradient, formulas (19) and (20) are equivalent, but a hyperelastic strain-energy cannot be convex and it is not generally known whether the “sup inf ” inversion is legitimate (sup inf 6 inf sup). Saddle-point problem (20) is nonlinear and we simpli6ed it by applying the second-order method.
38
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
3.2. Second-order expansion of the Lagrangian In Section 2, we applied the in6mum problem to the second-order expansion of W with respect to the deformation gradient F. In this section, we apply the second-order procedure to saddle-point problem (20) by expanding the Lagrangian L to second-order with respect to (F; p) about a piecewise constant reference deformation gradient Fr and a piecewise constant pressure pr 9L r 9L + (p − p ) L(F; p) ≈ L(Fr ; pr ) + (F − Fr ) : 9F Fr ; pr 9p Fr ; pr 1 92 L + (F − Fr ) : : (F − Fr ) 2 9F2 F˜r ; p˜r 92 L r + (F − F ) : (p − pr ): (21) 9F9p F˜r ; p˜r Let us de6ne 9L = r ; 9F Fr ; pr
92 L = L˜ r ; 9F2 F˜r ; p˜r
92 L = −A˜ r : 9F9p F˜r ; p˜r
(22)
Then applying expansion (21) to L in (20) yields K p) ˜ F; L( K ≈ sup infper L(Fr ; pr ) + (F − Fr ) : r − (p − pr )(det Fr − 1) p∈P u∈K
1 r r r r r r ˜ ˜ + (F − F ) : L : (F − F ) − (F − F ) : A (p − p ) : 2
(23)
To determine the piecewise constant reference 6elds (Fr ; pr ; F˜ r ; p˜ r ), we study the stationarity of the homogenized Lagrangian de6ned by (23), in line with Ponte Casta˜neda and Tiberio (2000). Since one cannot cancel the derivatives with respect to F˜ r and p˜ r , the 6rst and the second derivatives about FK r and pK r are cancelled, which gives Fr = F˜ r = FK r
and
pr = p˜ r = pK r :
(24)
We note F = I + ∇u and we rewrite (23) as r K p) ˜ F; L( K = L(FK r ) − ∇u : r + pK r (det FK r − 1) 1 r r r r r r + ∇u : L : ∇u − ∇u : A pK 2
1 r + sup infper ∇u : [r − ∇u : Lr + pK r Ar ] + ∇u : Lr : ∇u 2 p∈P F∈K r (25) − p(det FK r − 1 + Ar : (∇u − ∇u )] :
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
39
(25) is a saddle-point problem for a linear and heterogeneous solid, the Euler–Lagrange equations of which are r
Div(Lr : ∇u + r − ∇u : Lr + pK r Ar − pAr ) = 0 r det FK r − 1 + Ar : (∇u − ∇u ) = b
∇u = ∇u
and
p = pK
F and p V0 —periodic;
in V0 ;
(26a)
in V0 ;
(26b)
in V0 ;
(26c)
· n V0 —antiperiodic:
(26d)
Eq. (26b) stands for the linearized incompressibility constraint in which b denotes the Lagrange multiplier associated with the constraint p = p, K which has to be taken into account when the variations with respect to p are calculated. By averaging (26b) in each phase one obtains b= cr (det FK r − 1); (27) r
Hill’s lemma, see Ponte Casta˜neda and Suquet (1998) and references therein, according to which the average of the internal work equals the macroscopic internal work ( : F = : F), applied to (26) yields the following formula for the homogenized Lagrangian: 1 9W K r K ˜ L(F; p) K = cr W (FK r ) + (FK − FK r ) : (F ) 2 9F r − pK
r
1 det FK r − 1 + Ar : (FK − FK r ) 2
;
(28)
as well as the homogenized stress ˜ K ˜ F) K = 9L(F) : ( 9FK
(29)
3.3. Solving the linear heterogeneous saddle-point problem In Section 3.2, initial nonlinear saddle-point problem (20) was replaced by linear problem (25) or (26). In the present section, this heterogeneous problem is solved by using an estimate of the Hashin–Shtrikman kind. The superscript 0 is used to denote the reference medium taken to be the matrix. Reference operators L0 and A0 are introduced and problem (26) is re-written as Div(L0 : F + (Lr − L0 ) : F + Tr − pA0 − p(Ar − A0 )) = 0 −A0 : F − (Ar − A0 ) : F + k r = b
in V0 :
in V0 ; (30)
We then make the assumption that the terms arising from the heterogeneity of the material are piecewise constant, which simpli6es (30) into a linear homogeneous problem
40
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
with piecewise eigenstress de6ned by Div(L0 : F + pr − pA0 ) = 0 in V0 ; pr = Tr + (Lr − L0 ) : FK r − pK r (Ar − A0 )) A :F=c in V0 ; r c = −b − (Ar − A0 ) : FK r + k r 0
in V0 ;
r
(31)
in V0 :
Problem (31) is solved in the same way as was (8) in Section 2.2.2 with the additional unknown p and an additional equation linearizing the incompressibility constraint. Applying the Fourier transform (de6ned by (10)) to (31) yields pr 0 −L0 h j uk () − ij p()A ˆ ijkh ij = −ij ij (); (32) A0 uˆ () = cr (): ij j i
Eq. (32) is a system of linear equations for the unknown uˆ and p. ˆ Denoting
pr
L0 iA0 () ˆ u() i ; Z= ; D= ; B= p() ˆ T 0T cr () i A 0 system (32) becomes
with B−1 =
Z = B−1 D
H() Q() N() O()
:
(33)
The Fourier transform of the displacement gradient and of the pressure is given as pr [kh () = −Hki ()j h r ∇u ij () + ih Qk ()c (); pr r p() ˆ = iNi ()j ij () + O()c ():
An inverse Fourier transform is now performed to obtain the 6elds of the displacement gradient and of the pressure. Finally, by averaging both expressions over each individual phase, the following expressions are obtained: s ∇u = ∇u + [Hsr : pr + Qsr cKr ]; pK s = pK +
r
[N sr : pr + Osr cKr ]:
r
Introducing the periodic operators Hsr h Hki j cr 'r e−i·Y r ei·Y s ; khij = − ∈R∗ −{0}
(34)
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
sr Qkh =i
41
h Qk cr 'r e−i·Y r ei·Y s ;
∈R∗ −{0}
Nsr ij = i
Ni j cr 'r e−i·Y r ei·Y s ;
∈R∗ −{0}
Osr =
Ocr 'r e−i·Y r ei·Y s ;
∈R∗ −{0} s
(34) is a system of implicit equations, the solution of which gives ∇u and pK s . The expression of the homogenized Lagrangian and of the homogenized stress can then be derived. Remark. The initial constrained in6mum problem has been replaced by a saddle-point problem containing a Lagrange multiplier p associated with the local incompressibility constraint det F=1. The above method gives a homogenized Lagrangian as a function of the macroscopic deformation gradient FK and of p. K One might wonder which constraint is associated with the Lagrange multiplier p. K In the two-dimensional case, a simple K ˜ can be written as p(det calculation shows that the ePects of pK in L K F−1), which shows that pK is the Lagrange multiplier associated with the macroscopic incompressibility constraint. In the three-dimensional case, no such calculation of this kind can be made and we do not know what constraint is associated with p. K To proceed, one can express W˜ as K = L( K 0): ˜ F; W˜ (F)
(35)
This particular choice for pK removes the associated constraint but introduces a new approximation into the model because it reduces the size of the set P. We can then arti6cially introduce the macroscopic incompressibility constraint by K 0) if det FK = 1; ˜ F; L( K = W˜ (F) (36) +∞ if det FK = 1: 3.4. Comparisons with 4nite element results This new model is applied to the composite de6ned above (microstructures S4 and S9), with the incompressible tension and the shear simulation as de6ned in Table 1. On the other hand, here we considered only incompressible Mooney–Rivlin materials de6ned by W (F) = C10 (I1 − 3) + C01 (I2 − 3); where I1 and I2 are the two 6rst invariants of the right Cauchy–Green tensor de6ned by C = FT · F. In Fig. 7, the estimates of the homogenized stress given by the modi6ed formulation (MSOE) are compared with those given by the FEM, in the case of the incompressible tension and of simple shear. In both cases, the maximum error is lower than 7%. These simulations show that the modi6ed formulation de6nitely
42
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 7. Homogenized stress ˜ 11 vs. load factor during simulated incompressible tension in (a) and ˜ 12 vs. load factor during simulated simple shear in (b) applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively), incompressible behavior. The contrast was equal to 20 and the MSOE is also given for the sake of comparison.
improves the accuracy of the estimate given by the second-order procedure in the case of incompressible materials. 4. Comparisons between the model and experimental data The aim of this section is to compare the homogenized behavior predicted by the modi6ed second-order estimate (MSOE) with experimental data for a composite in a quadruple shearing test. As previously seen, the MSOE predicts the macroscopic behavior of hyperelastic composites knowing: the behavior of each phase, through the strain-energy function and the microstructure shape, unit cell in its original con6guration. The specimens tested here were made of a rubbery matrix, the behavior of which is not known a priori, reinforced by steel wires, assumed to be rigid. In the 6rst subsection, we detail the experimental procedure used to identify the matrix behavior. In the second subsection, we describe the experiments performed with the composites and comparisons are made between the experimental data and the MSOE prediction. 4.1. Identi4cation of the elastomer matrix behavior The elastomer tested was a black carbon 6lled and vulcanized rubber manufactured by MICHELIN. Its hyperelastic properties were determined by a tension test and a quadruple shearing test. The aim was to 6nd a strain-energy function capable of predicting the behavior of the sample in these two experiments as accurately as possible. The manufacturer guaranteed that the elastomer was isotropic and incompressible, which would mean that the strain energy function would only depend on the two 6rst invariants I1 and I2 . Consequently, the stress contains only the 6rst derivative of W (strain-energy function) with respect to I1 and I2 (the superscript stands for the
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
43
Fig. 8. Tension test sample (in mm).
diPerence between the modeled stress and the experimentally measured stress) mod (X) = H (I1 ; I2 )
9I1 9I2 + K(I1 ; I2 ) − pF−T 9F 9F
(37)
with H (I1 ; I2 ) =
9W (X; F) 9I1
and
K(I1 ; I2 ) =
9W (X; F) : 9I2
H and K are the material functions to be determined. 4.1.1. Tension test Tension tests were carried out at room temperature using a classical specimen, as shown in Fig. 8, which was loaded with a tension machine, capacity 100 kN. We measured the total load (denoted N ) with a load cell. The strain 6eld on one face of the sample was also analyzed with an image analysis system not shown in the 6gure. Assumptions had to be made about this experiment in order to interpret the results. • First, it was assumed that the stress 6eld was homogeneous in each cross section (perpendicular to the tension direction, St = 20 × 10 mm2 ) so as to be able to easily
44
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
compute the stress as a function of N and St , the initial section: exp =
N e1 ⊗ e 1 St
(superscript exp denoted experimentally measured stress):
(38)
• Secondly, the elastomer was assumed to be isotropic and incompressible. In addition, the strain state was also assumed to be homogeneous in the cross section to determine the 3D strain 6eld with a 2D measurement. Under these assumptions, the deformation gradient is measured 1 F = e1 ⊗ e1 + √ (e2 ⊗ e2 + e3 ⊗ e3 ); I2 ( ) =
I1 ( ) =
2 + 2 ;
1 + 2 : 2
(39)
The strain rate was ˙ = 4 × 10−3 s−1 . 4.1.2. Shear test The second experiment is a shear test on the quadruple shearing device shown in Fig. 9, composed of four elastomer specimens vulcanized between four steel plates. This specimen was tested in the same tension machine used for the tension test, as shown in Fig. 10. The shearing of the elastomeric samples resulted from the traction applied to the two middle steel plates. The total load was measured with the load cell (N ), see Fig. 10, and the total displacement with a displacement transducer (YL with L = 535 mm, see Fig. 10). As in the tension test, we had to make assumptions about this experiment in order to interpret the results: • the stress 6eld was assumed to be homogeneous in each cross section (Ss = 80 × 193 mm2 ). It was then given by exp =
N e1 ⊗ e2 : 2Ss
(40)
• the strain 6eld was also assumed to be homogeneous in each sample and de6ned by (YL is the total displacement and h the sample’s thickness): F = I + e1 ⊗ e2 ;
=
YL ; 2h
I1 ( ) = I2 ( ) = 3 + 2 :
(41)
The strain rate was ˙ = 4 × 10−3 s−1 . 4.1.3. Strain-energy identi4cation With the above described experimental results, tension and shear test, we looked for the strain energy functions H and K. The error between the experimentally measured
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
45
Fig. 9. Quadruple shearing sample (in mm).
LOAD CELL
QUADRUPLE SHEARING SAMPLE DISPLACEMENT TRANSDUCER
Fig. 10. Quadruple shearing experiment.
stress, exp , and the stress from the model, mod , was de6ned as error =
i
exp mod (11 ( i ) − 11 ( i ))2 +
i
exp mod (12 ( i ) − 12 ( i ))2 ;
(42)
where the subscript i refers to each measurement point. We looked for H and K functions minimizing this error. We checked some known hyperelastic functions such
46
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 11. Matrix response and identi6cation with the erf strain energy function (Eq. (44)), tension test in (a) and shear test in (b).
as the generalized Mooney–Rivlin function (ABAQUS, 2001; Ogden, 1997) and the exponential function, see Lambert-Diani and Rey (1998), but none of them 6tted with the experiments well because the material response is highly nonlinear even at small strains. This can be seen from Fig. 11: the response, both in tension and in shear, is nonlinear below a threshold ( = 1:15 for tension and = 0:15 for shear) and becomes nearly linear at higher values. We therefore propose a strain energy de6ned H (I1 ; I2 )=a1 +
a2(1−e−a3 (I1 −3) ) √ ; I1 −3
K(I1 ; I2 )=a4 +
a5(1−e−a6 (I2 −3) ) √ : I2 −3
(43)
These functions H and K are composed of constant terms, a1 and a4 which give the large strain linear response, and nonlinear ones which cancel out when I1 and I2 tend to in6nity (both a3 and a6 are positive numbers). Integrating those derivatives gives the following expression for the strain energy function, called erf strain energy: √ a2 Erf ( a3 (I1 − 3)) W (I1 ; I2 ) = a1 (I1 − 3) + 2a2 I1 − 3 − √ a3 √ a5 Erf ( a6 (I2 − 3)) + a4 (I2 − 3) + 2a5 I2 − 3 − ; (44) √ a6 where Erf is the error function. In Fig. 11, the measured and modeled stresses are compared in tension and shear. The strain-energy function was found to model the response of the elastomer in a very satisfactory manner, the identi6ed parameters being: a1 =2:05 MPa; a2 =−1:46 MPa; a3 =1:13; a4 =0:03 MPa; a5 =−0:2 MPa and a6 =17:5. 4.2. Composite testing and comparisons with the model The composite was tested in the quadruple shearing experiment described below. The shape of the specimen was the same as in Fig. 9, but here, each elastomer was reinforced with 600 steel wires transversally arranged with respect to the shearing
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
47
Fig. 12. Picture of one side of the composite tested.
Fig. 13. Unit cell in its initial con6guration (V0 ).
direction, 60 wires in the horizontal direction and 10 in the vertical direction. The wire radius was 0:72 mm and the microstructure (see the picture of a sample face in Fig. 12) was made of unit cells with a wire volume fraction equal to 0.254, as de6ned in Fig. 13. In Fig. 14, the stress measured during the quadruple shearing experiment is compared with that predicted by the model, where the only inputs were the matrix behavior and the unit cell shape in its reference con6guration. The prediction is seen to be very accurate: the maximum error is less than 7% at about = 0:2. This result looks satisfactory, but two assumptions made have to be kept in mind: • The 6rst concerns the periodicity of the solution. In fact, linear problem (31) was solved by making the assumption that the displacement and the pressure 6eld are both periodic. This assumption should be satis6ed for composites made of an in6nite number of unit cells, which is not the case for the composite tested here, namely 60 unit cells in the horizontal direction but only 10 in the vertical direction. • The second concerns the hyperelastic behavior. Experiments showed that the behavior of the elastomer tested was inTuenced by the strain rate. All the results given here are de6ned by the following loading rate: ˙ = ˙ = 4 × 10−3 s−1 . A viscohyperelastic behavior might give a more realistic description of the response of this elastomer.
48
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 14. Comparison between the homogenized stress predicted by the MSOE and the experimentally measured stress.
5. Conclusion In this paper, based on the second-order method proposed by Ponte Casta˜neda (1996) and Ponte Casta˜neda and Tiberio (2000), two estimates of the macroscopic behavior of hyperelastic composites were obtained. The 6rst was obtained by directly applying the method proposed by Ponte Casta˜neda and Tiberio (2000) to the case of periodic composites, using the periodic Green operators following the lines of Suquet (1990) or Iwakuma and Nemat-Nasser (1983). Two-dimensional 6nite element simulations showed that the estimate was fairly accurate in the case of compressible composites, but not with incompressible ones. A modi6ed method, based on a mixed formulation, namely displacement–pressure of the equilibrium equations was then proposed for incompressible materials. Other FE simulations have shown that this modi6ed method accurately predicts the macroscopic behavior of incompressible composites. Finally, this method was applied to the case of an industrial composite made of a rubber matrix reinforced with steel wires. Comparisons were made between the homogenized stress predicted by the modi6ed second-order method and the experimentally measured data. The response of the composite tested was well predicted using the method proposed. One of the approximations made rests on the modeling of elastomer materials as hyperelastic material: the behavior of these materials seems to be time dependent and the use of a viscohyperelastic model is suggested. This is left for future work. Acknowledgements Authors thank F. Barbarin, Michelin Co., for providing the specimens, and P. Suquet for all his useful advice on some of the theoretical points.
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
49
References ABAQUS, 2001. ABAQUS/Standard Theory Manual Version 6.2, Hibbit, Karlson & Sorensen, Rhodes Island. Abeyaratne, R., Triantafyllidis, N., 1984. An investigation of localization in a porous elastic material using homogenization theory. ASME J. Appl. Mech. 51, 481–486. Devries, F., 1996. Calcul du comportement homog[en[eis[e de composites hyper[elastiques. Rev. Composites Mat[er. Av. 2, 217–248. Dumontet, H., 1983. Homog[en[eisation de mat[eriaux strati6[es de type e[ lastique lin[eaire, non lin[eaire et visco-[elastique. Ph.D. Thesis, Universit[e P. et M. Curie. Dumontet, H., 1989. Homog[en[eisation et ePets de bords dans les mat[eriaux composites. Ph.D. Thesis, Universit[e P. et M. Curie. Geymonat, G., MBuller, S., Triantafyllidis, N., 1993. Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity. Arch. Ration. Mech. Anal. 122, 231–290. Iwakuma, T., Nemat-Nasser, S., 1983. Composites with periodic microstructure. Comput. Struct. 16, 13–19. Lambert-Diani, J., Rey, C., 1998. Elaboration de nouvelles lois de comportement pour les e[ lastom\eres: principe et avantages. C. R. Acad. Sci. Paris S[er. IIb 326 (8), 483–488. Marcellini, P., 1978. Periodic solutions and homogenization of non linear variational problems. Ann. Mat. Pura Appl. 117, 139–152. MBuller, S., 1987. Homogenization of nonconvex integral functionals and cellular elastic material. Arch. Ration. Mech. Anal. 99, 189–212. Ogden, R.W., 1997. Non-Linear Elastic Deformations. Dover, New York. Ponte Casta˜neda, P., 1989. The overall constitutive behaviour of nonlinearly elastic composites. Proc. R. Soc. London A 422, 147–171. Ponte Casta˜neda, P., 1996. Exact second-order estimates for the ePective mechanical properties of nonlinear composite materials. J. Mech. Phys. Solids 44, 827–862. Ponte Casta˜neda, P., Suquet, P., 1998. Nonlinear composites. Adv. Appl. Mech. 34, 171–302. Ponte Casta˜neda, P., Tiberio, E., 2000. A second-order homogenization method in 6nite elasticity and applications to black-6lled elastomers. J. Mech. Phys. Solids 48, 1389–1411. Pruchnicki, E., 1998. Hyperelastic homogenized law for reinforced elastomer at 6nite strain with edge ePects. Acta Mech. 129, 139–162. Suquet, P., 1990. Une m[ethode simpli6[ee pour le calcul des propri[et[es e[ lastiques de mat[eriaux h[et[erog\enes a\ structure p[eriodique. C. R. Acad. Sci. Paris S[er. IIb 311, 769–774. Talbot, D.R.S., Willis, J.R., 1997. Bounds of third order for the overall response of nonlinear composites. J. Mech. Phys. Solids 45, 87–111. Willis, J.R., 1981. Variational and related methods for the overall properties of composites. Adv. Appl. Mech. 21, 1–78.
www.elsevier.com/locate/jmps
Second-order estimate of the macroscopic behavior of periodic hyperelastic composites: theory and experimental validation N. Lahellec∗ , F. Mazerolle, J.C. Michel L.M.A./C.N.R.S., 31, Chemin Joseph Aiguier, Marseille 13402, Cedex 20, France Received 15 August 2002; received in revised form 30 June 2003; accepted 30 June 2003
Abstract This paper deals with some theoretical and experimental aspects of the behavior of periodic hyperelastic composites. We focus here on composites consisting of an elastomeric matrix periodically reinforced by long 6bers. The paper is composed of three parts. The 6rst part deals with the theoretical aspects of compressible behavior. The second-order theory of Ponte Casta˜neda (J. Mech. Phys. Solids 44 (1996) 827) is considered and extended to periodic microstructures. Comparisons with results obtained by the 6nite element method show that the composite behavior predicted by the present model is much more accurate for compressible than for incompressible materials. The second part deals with the extension of the method to incompressible behavior. A mixed formulation (displacement–pressure) is used which improves the accuracy of the estimate given by the model. The third part presents experimental results. The composite tested is made of a rubber matrix reinforced by steel wires. Firstly, the matrix behavior is identi6ed with a tensile test and a shear test carried out on homogeneous samples. Secondly, the composite is tested under shearing. The experimentally measured homogenized stress is then compared with the predictions of the model. ? 2003 Elsevier Ltd. All rights reserved. Keywords: Homogenization; Hyperelasticity; Rubber material; Experiments
1. Introduction The aim of this paper is to predict the behavior of hyperelastic composites using the homogenization procedure. Knowing the hyperelastic behavior of the composite’s constituents, it is proposed to predict the homogenized behavior as analytically as possible. ∗
Corresponding author. Fax: +33-491-1644-81. E-mail address: [email protected] (N. Lahellec).
0022-5096/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0022-5096(03)00104-2
28
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
When the strain energy is a convex function of the deformation gradient, the homogenized behavior is given by the homogenized strain energy. This function is computed by solving the in6mum strain energy of a single unit cell subjected to the macroscopic deformation gradient, see Marcellini (1978). Unfortunately, hyperelastic strain energies are not convex and, therefore, the in6mum problem cannot be studied in a single unit cell and has to be considered in the structure as a whole, see MBuller (1987). These calculations are very long and in the present study the homogenized strain energy is computed in a single unit cell, based on the assumption that no instability is observed in the range of loads studied. The nonlinearity of the behavior is an additional diCculty for hyperelastic materials. The relationship between stress and strain is strongly nonlinear, which complicates the estimation of their macroscopic behavior. Previous authors, Dumontet (1983, 1989), Geymonat et al. (1993) and Pruchnicki (1998) have studied the case of laminates where the problem is one dimensional. In the compressible case, Dumontet (1983) gave an analytical expression for the homogenized strain energy and in the incompressible case, Pruchnicki (1998) shows that the homogenization problem can be reduced to a simple nonlinear problem. Other authors have linearized the initial problem using an incremental method (Abeyaratne and Triantafyllidis, 1984; Devries, 1996) to approach more complex microstructures, and solved the linear problem using the periodic homogenization method. These methods always lead to a 6nite element method (FEM) resolution. The last group of studies dealt with developing analytical methods for estimating the macroscopic behavior of nonlinear composites, see Ponte Casta˜neda and Suquet (1998), Talbot and Willis (1997). Ponte Casta˜neda (1989) studied Voigt and Reuss bounds for hyperelastic materials and adapted the Hashin–Shtrikman procedure to get a nonlinear bound. Ponte Casta˜neda and Tiberio (2000) applied the second-order theory, 6rst proposed by Ponte Casta˜neda (1996), to hyperelastic composites. Here, we propose to extend this second-order method to estimate the macroscopic behavior of periodic and incompressible composites. We model the periodic microstructure using the lower periodic Hashin–Shtrikman bound proposed by Suquet (1990) and Iwakuma and Nemat-Nasser (1983). This bound is given by a simpli6ed computation of the linear problem in the Fourier space. To overcome the incompressibility constraint, we use a mixed formulation, displacement–pressure, before applying the second-order method. This paper concludes with an experimental part. The method developed is compared with experimental results obtained with a composite consisting of a rubber matrix reinforced by steel wires. Firstly, the matrix behavior is identi6ed with an original strain energy function. Secondly, the estimate of the homogenized stress is compared with the experimental data obtained in a shearing experiment. 2. Hyperelastic compressible and periodic composites 2.1. Second-order method P. Ponte Casta˜neda proposed the second-order method for estimating the macroscopic behavior of nonlinear composites, initially for power-law materials Ponte Casta˜neda
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
29
(1996), Ponte Casta˜neda and Suquet (1998) and more recently in the case of hyperelastic isotropic macroscopic behavior Ponte Casta˜neda and Tiberio (2000). When W denotes the strain energy and F the deformation gradient, the homogenized strain energy is given by the variational principle 1 K ˜ W (F) = inf W (X; F) dV ; (1) K |V0 | V F∈Kper (F) 0 K the set of kinematically where V0 is the unit cell in initial con6guration and Kper (F) admissible deformation gradients with the average deformation gradient FK K K u−(F−I)·X K Kper (F)={F | F=I+∇u; F=F; with : = (1=|V0 |) V0 : dV . The proposed method involves two stages:
V0 –periodic; det F¿0}
(2)
1. Linearization of the initial problem by the second-order procedure, which gives a linear but heterogeneous problem. 2. Resolution of this problem using the Hashin–Shtrikman periodic lower bound. The second-order theory consists of replacing W (F) by its second-order Taylor expansion about FK r which is the average of F over phase r. W (F) = W (FK r ) + (F − FK r ) : r + 12 (F − FK r ) : Lr : (F − FK r ); Lr =
r =
9W K r (F ); 9F
92 W K r (F ): 9F2
Using this expansion, Ponte Casta˜neda (1996) established that the second-order expansion of the homogenized strain energy, W˜ , can be written as K = W˜ (F) (3) cr W (FK r ) + 12 (FK − FK r ) : r : r
K The homogenized stress is then given by the derivative of (3) with respect to F: 1 ˜ K 9FK r ˜ F) K = 9W (F) = ( (4) + r : cr [r + (FK − FK r ) : Lr ] : 2 9FK 9FK r The average deformation gradient in each phase is still present in both expressions (3) and (4) and must be calculated. This is done by solving the following variational problem for a linear heterogeneous material: r inf ∇u : (r − Lr : ∇u ) + 12 ∇u : Lr : ∇u : (5) K F∈Kper (F)
The solution of (5) is studied in the next section.
30
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
2.2. Solving the linear heterogeneous problem 2.2.1. Hashin–Shtrikman lower bound In this section, we solve the linear heterogeneous problem (5) using the Hashin– Shtrikman procedure of Willis (1981). This estimate is obtained with a linear homogeneous reference material, having an elastic strain energy de6ned by W 0 = 12 ∇u : L0 : ∇u:
(6)
The Hashin–Shtrikman estimate is given by r
W HS (∇u) = 12 ∇u : L0 : ∇u + ∇u : (Lr − L0 )∇u :
(7)
When L0 is the matrix modulus (the most compliant constituent), this estimate gives a lower bound. The original linear heterogeneous problem is transformed into a linear homogeneous problem de6ned by
1 r ∇u : (r − L0 : ∇u ) + ∇u : L0 : ∇u : inf (8) K 2 F∈Kper (F) 2.2.2. Solving the periodic and homogeneous problem Ponte Casta˜neda and Tiberio (2000) proposed to use the Green operator in an in6nite medium for random composites. The composites considered here are periodic and we use the periodic Green operator, which can be determined by carrying out a Fourier transform of (8) in line with the procedure used in Suquet (1990) or Iwakuma and Nemat-Nasser (1983). In indicial notation, (8) reads pr L0ijkh uk; hj + ij; j =0
with uk; ij =
9 2 uk 9Xi 9Xj
and
r
pr = r − L0 : ∇u :
(9)
The Fourier transform of (9) is carried out with spatial frequencies de6ned by i = 2P=Li in which P is an integer and Li the unit cell dimension in the i direction, with orthogonal symmetry axes. With V0 the unit cell in its reference con6guration, f(X) a V0 -periodic vector valued function, the Fourier transform of f is de6ned as i·X ˆ ˆ f() = f(X)e−i·X ; f(X) = f + f()e and ∈R∗ −{0}
[ 9f(X) 9Xj
ˆ = ij f():
With these relations, (9) becomes pr −L0ijkh h j uk () + ij ij () = 0:
(10)
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
31
Inverting the acoustic tensor · L0 · and diPerentiating with respect to h yields the following expression for u k; h : pr −1 0 u k; h () = −h [Lijkh h j ]ki j ij ():
(11)
pr We now apply relations (10) to u k; h () and ij () which gives the local transformation gradient value pr −i·X i·X uk; h (X) = uk; h − h [L0ijkh h j ]−1 e : (12) ki j ij (X)e ∈R∗ −{0}
∇u given by (12) is the solution of (9). Strain energy expression (3) contains the average displacement gradient in each phase which is obtained by averaging (12) s Psr : pr ∇u = ∇u + (13) r
with Psr the fourth-order tensor associated with the periodic Green operator and de6ned by −i·X Psr r ei·X s : (14) h [L0ijkh h j ]−1 khij = − ki j cr e ∈R∗ −{0}
This operator depends on the microstructure through e−i·X r ei·X s and on the mas terial nonlinearity through L0 . At this stage, a relation between ∇u , Psr and pr has been established. Since these two terms depend on FK r , (13) is an implicit equation for r the unknown ∇u . K Expression (4) for the homogenized stress contains the fourth-order tensor 9FK r =9F, K which is calculated by taking the derivative of Eq. (13) with respect to F: pr 9FK s 9Psr 9FK r 9FK r pr sr 9 P : : (15) =I+ + : : : 9FK 9FK r 9FK 9FK r 9FK r 2.3. Comparison with 4nite element results In this section, the case of a periodic two-phase material composed of a matrix reinforced with long 6bers with a circular section is investigated. In the initial con6guration, the 6bers were assumed to be arranged in a square array. Two 6ber volume fractions are considered: 19 and 14 . In what follows, the unit cells associated with these microstructures are named S9 and S4. Both the 6bers and the matrix are modeled as Mooney–Rivlin hyperelastic material with strain energy (see ABAQUS, 2001): W (F) = C10 (I1 − 3) + C01 (I2 − 3) +
1 (det(F) − 1)2 ; d
(16)
where F =
F ; det(F)1=3
C = F T · F ;
I1 = tr(C );
I2 = 12 (tr(C )2 − tr(C · C )):
32
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Table 1 De6nitions of the studied loadings Name
Macroscopic deformation gradient
Uniaxial tension Incompressible tension Simple shear Compression
FK = e1 ⊗ e1 + e2 ⊗ e2 + e3 ⊗ e3 ; 1 6 6 3 FK = e1 ⊗ e1 + 1 e2 ⊗ e2 + e3 ⊗ e3 ; 1 6 6 3 FK = e1 ⊗ e1 + e2 ⊗ e2 + e3 ⊗ e3 + e1 ⊗ e2 ; 1 6 6 2 FK = !e1 ⊗ e1 + e2 ⊗ e2 + e3 ⊗ e3 ; 0:3 6 ! 6 1
Fig. 1. Typical mesh (4-node bilinear elements), S4 unit cell.
Note that quasi-incompressible behavior can be modeled by increasing the ratio 1=d (penalty method). Note also that in the limiting case of small strains, these materials can be regarded as linear elastic and isotropic with bulk modulus k = 2=d and shear modulus = 2(C10 + C01 ). The contrast between the 6bers and the matrix is de6ned as W 6bers (F) = contrast × W matrix (F): Only in-plane macroscopic strains are considered here, (O; X1 ; X2 ) parallel planes and loading are de6ned in Table 1, and the problem is solved in plane strain. The homogenized strain energy and stress obtained using the second-order procedure are compared with those resulting from 6nite element calculations (ABAQUS, 2001). A typical mesh using four-node bilinear elements is shown in Fig. 1. 2.3.1. Compressible case The compressible case is solved 6rst. The matrix parameters are matrix matrix C10 = C01 = (1=d)matrix = 1 Pa:
Uniaxial tension: The uniaxial tension case corresponds to an average deformation gradient FK given in Table 1. The maximum strain is 200% ( = 3). In both
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
33
Fig. 2. Homogenized stress, ˜ 11 in (a) and ˜ 22 in (b), vs. load factor during simulated uniaxial tension applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively). The contrast between 6bers and matrix was 20.
microstructures, the 6bers are 20 times stiPer than the matrix. The stress estimate given by the second-order procedure (SOE) is compared with that given by the FEM in Fig. 2. With S9 and S4 microstructures, it can be seen that the diPerence between the two estimates is very small. The greatest error in the stress is less than ˜ FEM ˜ FEM 2%, (error = (˜ SOE 11 − 11 )= 11 ). In all the cases studied, the maximum error was reached with the greatest 6ber ratio, 25% volume fraction and with the largest strain. Compression: We simulated the compression test (Table 1) with microstructures S9 and S4, with two contrasts, ! is the load factor: • Contrast equal to 20, see Figs. 3(a) and (b): the greatest error in the estimate of 11 was nearly 1%. With 22 , the error was of the same order, above ! ≈ 0:55 in the S4 case and ! ≈ 0:38 for S9. Below these values, the error increased to 10% with S4. • Contrast equal to 100: Figs. 3(c) and (d) show that the estimates given by the second-order procedure are less satisfactory. The error was still low, above ! ≈ 0:65 with S4 and ! ≈ 0:5 with S9, but increased below these values to 25% with 11 and more than 50% with 22 . Simple shear: We simulated the simple shear test de6ned by FK in Table 1 up to 200% shear strain ( = 2). Microstructures S4 and S9 were taken to have a contrast equal to 20. In Fig. 4, the SOE and the FEM are compared. In the case of S9, the 12 and 21 estimates were accurate, and the greatest error was less than 3%. The estimate was also accurate in the case of S4 and at a shear strain less than 100%. Above this strain value, the error amounted to 15% at the maximum strain ( = 2). Conclusion for the compressible case: The nonlinearity of the behavior of hyperelastic composites makes it diCcult to calculate their homogenized or ePective behavior. This nonlinearity is due to the large strains involved and the resulting changes in the microstructure. The present model takes into account the nonlinearity of the material with the inTuence of the average deformation gradient in the matrix phase, in L0 and
34
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 3. Homogenized stress vs. load factor ! during simulated uniaxial compression applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively). The contrast between 6bers and matrix was 20 in (a) and (b) and 100 in (c) and (d).
Fig. 4. Homogenized stress, ˜ 12 in (a) and ˜ 21 in (b), vs. load factor during simulated simple shear applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively). The contrast between 6bers and matrix was 20.
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
35
Fig. 5. Pre-compressed shear simulation and description of microstructure evolution, with microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively). The contrast between 6bers and matrix was 20.
the various eigenstresses. The simulation, which gave very accurate estimates, uniaxial tension simulation, showed that this model accounts accurately for the nonlinearity when the strain 6eld is not localized. In the simulation, the 6ber interaction zone, namely the area in between two 6bers, is large and does not shrink with the loading factor. In the case of compression, the interactions between 6bers evolve. By increasing the loading factor, the 6bers move closer and the matrix between them is compressed and becomes stiPer, because the volume cannot be compressed to zero (det F 9 0). With small contrasts, this ePect is not very important because the interaction zone is large enough to inTuence the average matrix behavior in FK 2 . On the other hand, increasing the contrast decreases the interaction zone as well as its ePect on the average matrix behavior. In the limiting case of rigid 6bers, in6nite contrast, the composite locks up when the 6bers come into contact. With the microstructures modeled here, max 1 = 2 c1 =, where square unit cells with circular 6ber, this locking occurs when F11 max 1 max 1 (S4) = 0:56; F11 (S9) = 0:38). The model c1 is the 6ber volume fraction (F11 2 max 2 K = c1 , which is smaller proposed here locks up when det F = 0, which occurs at F11 1 max . Simulation of simple shear shows changes in the microstructure shape than F11 which are not very well taken into account by the proposed model. In this case, the microstructure alternates between a con6guration in which the 6bers are aligned and one in which they are arranged crosswise. This tendency is even more pronounced in a pre-compressed shear simulation with the deformation gradient de6ned as FK = e1 ⊗ e1 + 12 e2 ⊗ e2 + e3 ⊗ e3 + e1 ⊗ e2 ;
: 0 → 2:
(17)
During this simulation, every 6ber is lined up with either another 6ber or between two others (see Fig. 5). The macroscopic behavior depends strongly on the microstructure
36
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 6. Homogenized stress ˜ 11 vs. load factor during simulated incompressible tension in (a) and ˜ 21 vs. load factor during simulated simple shear in (b), applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 11% and 25%, respectively), incompressible behavior. The contrast between 6bers and matrix was 20.
changes. In the case of S9, SOE accurately predicts the two 6rst 6bers alignments, whereas in the case of S4, the 6rst alignment is well predicted but then the error becomes very large. 2.3.2. Quasi-incompressible case The aim of this study was to model composites made from incompressible rubber. In order to check the accuracy of the SOE in this case, simulations were performed with the same values for C10 and C01 as in the previous case, increasing the ratio 1=d so as to reach the incompressibility limit. The simulations shown here were carried out with composites de6ned by 6bers 20 times stiPer than the matrix (contrast equal to 20), which are also de6ned by C10 = C01 = 1 Pa and 1=d = 500 Pa. Fig. 6 shows the evolution of the tensile stress and of the shear stress corresponding to incompressible tension and shear as de6ned in Table 1. In both cases, the results of the SOE deviate signi6cantly from the FEM values. Conclusion for quasi-incompressible composites: The above simulations show that as we get closer to the incompressibility limit, the estimate given by the classic second-order formulation is no longer satisfactory. For both microstructures, we varied the 1=d values from 0.1 to 1000 Pa and the accuracy of the estimate was found to decrease sharply at around 500 Pa. This can be explained by the strong nonlinearity introduced by the incompressibility constraint. With materials of this kind, the strain energy can be written as follows: C10 (I1 − 3) + C01 (I2 − 3) if det(F) = 1; (18) W (F) = +∞ if det(F) = 1: This strong nonlinearity of det(F) is not accounted for by a second-order Taylor expansion, which explains these poor results. This was noticed by Ponte Casta˜neda and Tiberio (2000), who established, in the case of rigid inclusions that the straightfor-
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
37
ward use of the second-order approximation led to an approximate incompressibility constraint which is inconsistent with the exact one. 2.4. Conclusion The above simulations (in the compressible and quasi-incompressible cases) show that SOE is an accurate tool for dealing with compressible periodic composites in loading cases where the macroscopic behavior does not depend too strongly on 6ber interactions. In the next section, the second-order method is adapted to model incompressible composites.
3. Modi ed method for incompressible hyperelastic materials 3.1. New formulation It was shown in Section 2.3 that the traditional second-order formulation is not appropriate for dealing with incompressible materials. Ponte Casta˜neda and Tiberio (2000) developed a method of computing the incompressible limit for rigid inclusions, but it cannot be applied here (because of the 6nite contrast). Therefore, an alternative formulation has to be developed. The initial variational principle reads K = W˜ (F)
inf
K F∈Kper; inc (F)
W (F);
(19)
K denotes the set of kinematically admissible deformation gradients where Kper; inc (F) with the average deformation gradient FK and the incompressibility constraint K = {F | F = I + ∇u; F Kper; inc (F) K u − (FK − I) · X V0 –periodic and det(F) = 1}: = F; We replace this constrained in6mum problem by a saddle-point problem, where the incompressibility constraint is relaxed by introducing its associated Lagrange multiplier p, the pressure. K p) ˜ F; L( K = sup
inf
K p∈P(p) K F∈Kper (F)
L(F; p);
L(F; p) = W (F) − p(det(F) − 1); (20)
with P as the set of admissible pressures p. P(p) K = {p | p V0 –periodic; p = p}: K When W is a convex function of the deformation gradient, formulas (19) and (20) are equivalent, but a hyperelastic strain-energy cannot be convex and it is not generally known whether the “sup inf ” inversion is legitimate (sup inf 6 inf sup). Saddle-point problem (20) is nonlinear and we simpli6ed it by applying the second-order method.
38
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
3.2. Second-order expansion of the Lagrangian In Section 2, we applied the in6mum problem to the second-order expansion of W with respect to the deformation gradient F. In this section, we apply the second-order procedure to saddle-point problem (20) by expanding the Lagrangian L to second-order with respect to (F; p) about a piecewise constant reference deformation gradient Fr and a piecewise constant pressure pr 9L r 9L + (p − p ) L(F; p) ≈ L(Fr ; pr ) + (F − Fr ) : 9F Fr ; pr 9p Fr ; pr 1 92 L + (F − Fr ) : : (F − Fr ) 2 9F2 F˜r ; p˜r 92 L r + (F − F ) : (p − pr ): (21) 9F9p F˜r ; p˜r Let us de6ne 9L = r ; 9F Fr ; pr
92 L = L˜ r ; 9F2 F˜r ; p˜r
92 L = −A˜ r : 9F9p F˜r ; p˜r
(22)
Then applying expansion (21) to L in (20) yields K p) ˜ F; L( K ≈ sup infper L(Fr ; pr ) + (F − Fr ) : r − (p − pr )(det Fr − 1) p∈P u∈K
1 r r r r r r ˜ ˜ + (F − F ) : L : (F − F ) − (F − F ) : A (p − p ) : 2
(23)
To determine the piecewise constant reference 6elds (Fr ; pr ; F˜ r ; p˜ r ), we study the stationarity of the homogenized Lagrangian de6ned by (23), in line with Ponte Casta˜neda and Tiberio (2000). Since one cannot cancel the derivatives with respect to F˜ r and p˜ r , the 6rst and the second derivatives about FK r and pK r are cancelled, which gives Fr = F˜ r = FK r
and
pr = p˜ r = pK r :
(24)
We note F = I + ∇u and we rewrite (23) as r K p) ˜ F; L( K = L(FK r ) − ∇u : r + pK r (det FK r − 1) 1 r r r r r r + ∇u : L : ∇u − ∇u : A pK 2
1 r + sup infper ∇u : [r − ∇u : Lr + pK r Ar ] + ∇u : Lr : ∇u 2 p∈P F∈K r (25) − p(det FK r − 1 + Ar : (∇u − ∇u )] :
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
39
(25) is a saddle-point problem for a linear and heterogeneous solid, the Euler–Lagrange equations of which are r
Div(Lr : ∇u + r − ∇u : Lr + pK r Ar − pAr ) = 0 r det FK r − 1 + Ar : (∇u − ∇u ) = b
∇u = ∇u
and
p = pK
F and p V0 —periodic;
in V0 ;
(26a)
in V0 ;
(26b)
in V0 ;
(26c)
· n V0 —antiperiodic:
(26d)
Eq. (26b) stands for the linearized incompressibility constraint in which b denotes the Lagrange multiplier associated with the constraint p = p, K which has to be taken into account when the variations with respect to p are calculated. By averaging (26b) in each phase one obtains b= cr (det FK r − 1); (27) r
Hill’s lemma, see Ponte Casta˜neda and Suquet (1998) and references therein, according to which the average of the internal work equals the macroscopic internal work ( : F = : F), applied to (26) yields the following formula for the homogenized Lagrangian: 1 9W K r K ˜ L(F; p) K = cr W (FK r ) + (FK − FK r ) : (F ) 2 9F r − pK
r
1 det FK r − 1 + Ar : (FK − FK r ) 2
;
(28)
as well as the homogenized stress ˜ K ˜ F) K = 9L(F) : ( 9FK
(29)
3.3. Solving the linear heterogeneous saddle-point problem In Section 3.2, initial nonlinear saddle-point problem (20) was replaced by linear problem (25) or (26). In the present section, this heterogeneous problem is solved by using an estimate of the Hashin–Shtrikman kind. The superscript 0 is used to denote the reference medium taken to be the matrix. Reference operators L0 and A0 are introduced and problem (26) is re-written as Div(L0 : F + (Lr − L0 ) : F + Tr − pA0 − p(Ar − A0 )) = 0 −A0 : F − (Ar − A0 ) : F + k r = b
in V0 :
in V0 ; (30)
We then make the assumption that the terms arising from the heterogeneity of the material are piecewise constant, which simpli6es (30) into a linear homogeneous problem
40
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
with piecewise eigenstress de6ned by Div(L0 : F + pr − pA0 ) = 0 in V0 ; pr = Tr + (Lr − L0 ) : FK r − pK r (Ar − A0 )) A :F=c in V0 ; r c = −b − (Ar − A0 ) : FK r + k r 0
in V0 ;
r
(31)
in V0 :
Problem (31) is solved in the same way as was (8) in Section 2.2.2 with the additional unknown p and an additional equation linearizing the incompressibility constraint. Applying the Fourier transform (de6ned by (10)) to (31) yields pr 0 −L0 h j uk () − ij p()A ˆ ijkh ij = −ij ij (); (32) A0 uˆ () = cr (): ij j i
Eq. (32) is a system of linear equations for the unknown uˆ and p. ˆ Denoting
pr
L0 iA0 () ˆ u() i ; Z= ; D= ; B= p() ˆ T 0T cr () i A 0 system (32) becomes
with B−1 =
Z = B−1 D
H() Q() N() O()
:
(33)
The Fourier transform of the displacement gradient and of the pressure is given as pr [kh () = −Hki ()j h r ∇u ij () + ih Qk ()c (); pr r p() ˆ = iNi ()j ij () + O()c ():
An inverse Fourier transform is now performed to obtain the 6elds of the displacement gradient and of the pressure. Finally, by averaging both expressions over each individual phase, the following expressions are obtained: s ∇u = ∇u + [Hsr : pr + Qsr cKr ]; pK s = pK +
r
[N sr : pr + Osr cKr ]:
r
Introducing the periodic operators Hsr h Hki j cr 'r e−i·Y r ei·Y s ; khij = − ∈R∗ −{0}
(34)
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
sr Qkh =i
41
h Qk cr 'r e−i·Y r ei·Y s ;
∈R∗ −{0}
Nsr ij = i
Ni j cr 'r e−i·Y r ei·Y s ;
∈R∗ −{0}
Osr =
Ocr 'r e−i·Y r ei·Y s ;
∈R∗ −{0} s
(34) is a system of implicit equations, the solution of which gives ∇u and pK s . The expression of the homogenized Lagrangian and of the homogenized stress can then be derived. Remark. The initial constrained in6mum problem has been replaced by a saddle-point problem containing a Lagrange multiplier p associated with the local incompressibility constraint det F=1. The above method gives a homogenized Lagrangian as a function of the macroscopic deformation gradient FK and of p. K One might wonder which constraint is associated with the Lagrange multiplier p. K In the two-dimensional case, a simple K ˜ can be written as p(det calculation shows that the ePects of pK in L K F−1), which shows that pK is the Lagrange multiplier associated with the macroscopic incompressibility constraint. In the three-dimensional case, no such calculation of this kind can be made and we do not know what constraint is associated with p. K To proceed, one can express W˜ as K = L( K 0): ˜ F; W˜ (F)
(35)
This particular choice for pK removes the associated constraint but introduces a new approximation into the model because it reduces the size of the set P. We can then arti6cially introduce the macroscopic incompressibility constraint by K 0) if det FK = 1; ˜ F; L( K = W˜ (F) (36) +∞ if det FK = 1: 3.4. Comparisons with 4nite element results This new model is applied to the composite de6ned above (microstructures S4 and S9), with the incompressible tension and the shear simulation as de6ned in Table 1. On the other hand, here we considered only incompressible Mooney–Rivlin materials de6ned by W (F) = C10 (I1 − 3) + C01 (I2 − 3); where I1 and I2 are the two 6rst invariants of the right Cauchy–Green tensor de6ned by C = FT · F. In Fig. 7, the estimates of the homogenized stress given by the modi6ed formulation (MSOE) are compared with those given by the FEM, in the case of the incompressible tension and of simple shear. In both cases, the maximum error is lower than 7%. These simulations show that the modi6ed formulation de6nitely
42
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 7. Homogenized stress ˜ 11 vs. load factor during simulated incompressible tension in (a) and ˜ 12 vs. load factor during simulated simple shear in (b) applied to microstructures S9 and S4 (6bers in square array with 6ber volume fraction equal to 19 and 14 , respectively), incompressible behavior. The contrast was equal to 20 and the MSOE is also given for the sake of comparison.
improves the accuracy of the estimate given by the second-order procedure in the case of incompressible materials. 4. Comparisons between the model and experimental data The aim of this section is to compare the homogenized behavior predicted by the modi6ed second-order estimate (MSOE) with experimental data for a composite in a quadruple shearing test. As previously seen, the MSOE predicts the macroscopic behavior of hyperelastic composites knowing: the behavior of each phase, through the strain-energy function and the microstructure shape, unit cell in its original con6guration. The specimens tested here were made of a rubbery matrix, the behavior of which is not known a priori, reinforced by steel wires, assumed to be rigid. In the 6rst subsection, we detail the experimental procedure used to identify the matrix behavior. In the second subsection, we describe the experiments performed with the composites and comparisons are made between the experimental data and the MSOE prediction. 4.1. Identi4cation of the elastomer matrix behavior The elastomer tested was a black carbon 6lled and vulcanized rubber manufactured by MICHELIN. Its hyperelastic properties were determined by a tension test and a quadruple shearing test. The aim was to 6nd a strain-energy function capable of predicting the behavior of the sample in these two experiments as accurately as possible. The manufacturer guaranteed that the elastomer was isotropic and incompressible, which would mean that the strain energy function would only depend on the two 6rst invariants I1 and I2 . Consequently, the stress contains only the 6rst derivative of W (strain-energy function) with respect to I1 and I2 (the superscript stands for the
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
43
Fig. 8. Tension test sample (in mm).
diPerence between the modeled stress and the experimentally measured stress) mod (X) = H (I1 ; I2 )
9I1 9I2 + K(I1 ; I2 ) − pF−T 9F 9F
(37)
with H (I1 ; I2 ) =
9W (X; F) 9I1
and
K(I1 ; I2 ) =
9W (X; F) : 9I2
H and K are the material functions to be determined. 4.1.1. Tension test Tension tests were carried out at room temperature using a classical specimen, as shown in Fig. 8, which was loaded with a tension machine, capacity 100 kN. We measured the total load (denoted N ) with a load cell. The strain 6eld on one face of the sample was also analyzed with an image analysis system not shown in the 6gure. Assumptions had to be made about this experiment in order to interpret the results. • First, it was assumed that the stress 6eld was homogeneous in each cross section (perpendicular to the tension direction, St = 20 × 10 mm2 ) so as to be able to easily
44
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
compute the stress as a function of N and St , the initial section: exp =
N e1 ⊗ e 1 St
(superscript exp denoted experimentally measured stress):
(38)
• Secondly, the elastomer was assumed to be isotropic and incompressible. In addition, the strain state was also assumed to be homogeneous in the cross section to determine the 3D strain 6eld with a 2D measurement. Under these assumptions, the deformation gradient is measured 1 F = e1 ⊗ e1 + √ (e2 ⊗ e2 + e3 ⊗ e3 ); I2 ( ) =
I1 ( ) =
2 + 2 ;
1 + 2 : 2
(39)
The strain rate was ˙ = 4 × 10−3 s−1 . 4.1.2. Shear test The second experiment is a shear test on the quadruple shearing device shown in Fig. 9, composed of four elastomer specimens vulcanized between four steel plates. This specimen was tested in the same tension machine used for the tension test, as shown in Fig. 10. The shearing of the elastomeric samples resulted from the traction applied to the two middle steel plates. The total load was measured with the load cell (N ), see Fig. 10, and the total displacement with a displacement transducer (YL with L = 535 mm, see Fig. 10). As in the tension test, we had to make assumptions about this experiment in order to interpret the results: • the stress 6eld was assumed to be homogeneous in each cross section (Ss = 80 × 193 mm2 ). It was then given by exp =
N e1 ⊗ e2 : 2Ss
(40)
• the strain 6eld was also assumed to be homogeneous in each sample and de6ned by (YL is the total displacement and h the sample’s thickness): F = I + e1 ⊗ e2 ;
=
YL ; 2h
I1 ( ) = I2 ( ) = 3 + 2 :
(41)
The strain rate was ˙ = 4 × 10−3 s−1 . 4.1.3. Strain-energy identi4cation With the above described experimental results, tension and shear test, we looked for the strain energy functions H and K. The error between the experimentally measured
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
45
Fig. 9. Quadruple shearing sample (in mm).
LOAD CELL
QUADRUPLE SHEARING SAMPLE DISPLACEMENT TRANSDUCER
Fig. 10. Quadruple shearing experiment.
stress, exp , and the stress from the model, mod , was de6ned as error =
i
exp mod (11 ( i ) − 11 ( i ))2 +
i
exp mod (12 ( i ) − 12 ( i ))2 ;
(42)
where the subscript i refers to each measurement point. We looked for H and K functions minimizing this error. We checked some known hyperelastic functions such
46
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 11. Matrix response and identi6cation with the erf strain energy function (Eq. (44)), tension test in (a) and shear test in (b).
as the generalized Mooney–Rivlin function (ABAQUS, 2001; Ogden, 1997) and the exponential function, see Lambert-Diani and Rey (1998), but none of them 6tted with the experiments well because the material response is highly nonlinear even at small strains. This can be seen from Fig. 11: the response, both in tension and in shear, is nonlinear below a threshold ( = 1:15 for tension and = 0:15 for shear) and becomes nearly linear at higher values. We therefore propose a strain energy de6ned H (I1 ; I2 )=a1 +
a2(1−e−a3 (I1 −3) ) √ ; I1 −3
K(I1 ; I2 )=a4 +
a5(1−e−a6 (I2 −3) ) √ : I2 −3
(43)
These functions H and K are composed of constant terms, a1 and a4 which give the large strain linear response, and nonlinear ones which cancel out when I1 and I2 tend to in6nity (both a3 and a6 are positive numbers). Integrating those derivatives gives the following expression for the strain energy function, called erf strain energy: √ a2 Erf ( a3 (I1 − 3)) W (I1 ; I2 ) = a1 (I1 − 3) + 2a2 I1 − 3 − √ a3 √ a5 Erf ( a6 (I2 − 3)) + a4 (I2 − 3) + 2a5 I2 − 3 − ; (44) √ a6 where Erf is the error function. In Fig. 11, the measured and modeled stresses are compared in tension and shear. The strain-energy function was found to model the response of the elastomer in a very satisfactory manner, the identi6ed parameters being: a1 =2:05 MPa; a2 =−1:46 MPa; a3 =1:13; a4 =0:03 MPa; a5 =−0:2 MPa and a6 =17:5. 4.2. Composite testing and comparisons with the model The composite was tested in the quadruple shearing experiment described below. The shape of the specimen was the same as in Fig. 9, but here, each elastomer was reinforced with 600 steel wires transversally arranged with respect to the shearing
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
47
Fig. 12. Picture of one side of the composite tested.
Fig. 13. Unit cell in its initial con6guration (V0 ).
direction, 60 wires in the horizontal direction and 10 in the vertical direction. The wire radius was 0:72 mm and the microstructure (see the picture of a sample face in Fig. 12) was made of unit cells with a wire volume fraction equal to 0.254, as de6ned in Fig. 13. In Fig. 14, the stress measured during the quadruple shearing experiment is compared with that predicted by the model, where the only inputs were the matrix behavior and the unit cell shape in its reference con6guration. The prediction is seen to be very accurate: the maximum error is less than 7% at about = 0:2. This result looks satisfactory, but two assumptions made have to be kept in mind: • The 6rst concerns the periodicity of the solution. In fact, linear problem (31) was solved by making the assumption that the displacement and the pressure 6eld are both periodic. This assumption should be satis6ed for composites made of an in6nite number of unit cells, which is not the case for the composite tested here, namely 60 unit cells in the horizontal direction but only 10 in the vertical direction. • The second concerns the hyperelastic behavior. Experiments showed that the behavior of the elastomer tested was inTuenced by the strain rate. All the results given here are de6ned by the following loading rate: ˙ = ˙ = 4 × 10−3 s−1 . A viscohyperelastic behavior might give a more realistic description of the response of this elastomer.
48
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
Fig. 14. Comparison between the homogenized stress predicted by the MSOE and the experimentally measured stress.
5. Conclusion In this paper, based on the second-order method proposed by Ponte Casta˜neda (1996) and Ponte Casta˜neda and Tiberio (2000), two estimates of the macroscopic behavior of hyperelastic composites were obtained. The 6rst was obtained by directly applying the method proposed by Ponte Casta˜neda and Tiberio (2000) to the case of periodic composites, using the periodic Green operators following the lines of Suquet (1990) or Iwakuma and Nemat-Nasser (1983). Two-dimensional 6nite element simulations showed that the estimate was fairly accurate in the case of compressible composites, but not with incompressible ones. A modi6ed method, based on a mixed formulation, namely displacement–pressure of the equilibrium equations was then proposed for incompressible materials. Other FE simulations have shown that this modi6ed method accurately predicts the macroscopic behavior of incompressible composites. Finally, this method was applied to the case of an industrial composite made of a rubber matrix reinforced with steel wires. Comparisons were made between the homogenized stress predicted by the modi6ed second-order method and the experimentally measured data. The response of the composite tested was well predicted using the method proposed. One of the approximations made rests on the modeling of elastomer materials as hyperelastic material: the behavior of these materials seems to be time dependent and the use of a viscohyperelastic model is suggested. This is left for future work. Acknowledgements Authors thank F. Barbarin, Michelin Co., for providing the specimens, and P. Suquet for all his useful advice on some of the theoretical points.
N. Lahellec et al. / J. Mech. Phys. Solids 52 (2004) 27 – 49
49
References ABAQUS, 2001. ABAQUS/Standard Theory Manual Version 6.2, Hibbit, Karlson & Sorensen, Rhodes Island. Abeyaratne, R., Triantafyllidis, N., 1984. An investigation of localization in a porous elastic material using homogenization theory. ASME J. Appl. Mech. 51, 481–486. Devries, F., 1996. Calcul du comportement homog[en[eis[e de composites hyper[elastiques. Rev. Composites Mat[er. Av. 2, 217–248. Dumontet, H., 1983. Homog[en[eisation de mat[eriaux strati6[es de type e[ lastique lin[eaire, non lin[eaire et visco-[elastique. Ph.D. Thesis, Universit[e P. et M. Curie. Dumontet, H., 1989. Homog[en[eisation et ePets de bords dans les mat[eriaux composites. Ph.D. Thesis, Universit[e P. et M. Curie. Geymonat, G., MBuller, S., Triantafyllidis, N., 1993. Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity. Arch. Ration. Mech. Anal. 122, 231–290. Iwakuma, T., Nemat-Nasser, S., 1983. Composites with periodic microstructure. Comput. Struct. 16, 13–19. Lambert-Diani, J., Rey, C., 1998. Elaboration de nouvelles lois de comportement pour les e[ lastom\eres: principe et avantages. C. R. Acad. Sci. Paris S[er. IIb 326 (8), 483–488. Marcellini, P., 1978. Periodic solutions and homogenization of non linear variational problems. Ann. Mat. Pura Appl. 117, 139–152. MBuller, S., 1987. Homogenization of nonconvex integral functionals and cellular elastic material. Arch. Ration. Mech. Anal. 99, 189–212. Ogden, R.W., 1997. Non-Linear Elastic Deformations. Dover, New York. Ponte Casta˜neda, P., 1989. The overall constitutive behaviour of nonlinearly elastic composites. Proc. R. Soc. London A 422, 147–171. Ponte Casta˜neda, P., 1996. Exact second-order estimates for the ePective mechanical properties of nonlinear composite materials. J. Mech. Phys. Solids 44, 827–862. Ponte Casta˜neda, P., Suquet, P., 1998. Nonlinear composites. Adv. Appl. Mech. 34, 171–302. Ponte Casta˜neda, P., Tiberio, E., 2000. A second-order homogenization method in 6nite elasticity and applications to black-6lled elastomers. J. Mech. Phys. Solids 48, 1389–1411. Pruchnicki, E., 1998. Hyperelastic homogenized law for reinforced elastomer at 6nite strain with edge ePects. Acta Mech. 129, 139–162. Suquet, P., 1990. Une m[ethode simpli6[ee pour le calcul des propri[et[es e[ lastiques de mat[eriaux h[et[erog\enes a\ structure p[eriodique. C. R. Acad. Sci. Paris S[er. IIb 311, 769–774. Talbot, D.R.S., Willis, J.R., 1997. Bounds of third order for the overall response of nonlinear composites. J. Mech. Phys. Solids 45, 87–111. Willis, J.R., 1981. Variational and related methods for the overall properties of composites. Adv. Appl. Mech. 21, 1–78.