Instability analysis of non-homogeneous materials under biaxial loading

International Journal of Plasticity 21 (2005) 1970–1999 www.elsevier.com/locate/ijplas

Instability analysis of non-homogeneous materials under biaxial loading George Chatzigeorgiou, Nicolas Charalambakis

*

Department of Civil Engineering, Aristotle University of Thessaloniki, Thessaloniki, GR 541 24, Greece Received in final revised form 22 June 2004 Available online 13 February 2005

Abstract In this paper, we formulate the instability problem of non-homogeneous materials under biaxial loading, in the sense that the mechanical parameters, and more specifically the strain hardening or softening and the strain-rate sensitivity, are spatially dependent functions. We first study the time behavior of the process and present the evolution of strain non-uniformities before the critical time, to show that it depends on the interplay between the spatial distribution of mechanical parameters and the initial non-uniformities. We next study the instability modes of a material exhibiting inhomogeneities along one direction by an effective instability analysis [Dudzinski, D., Molinari, A., 1991. Perturbation analysis of thermoviscoplastic instabilities in biaxial loading. International Journal of Solids and Structures 27 (5), 601–628] adapted to the non-homogeneous case. This method is based on a suitable nonuniform reference solution and allows to select the localization modes activated at different deformation levels for different ‘‘deformation paths’’ of the material zones.
2005 Elsevier Ltd. All rights reserved. Keywords: Non-homogeneity; Functionally graded materials; Non-uniformity; Necking

*

Corresponding author. Tel.: +30 231 995931; fax: +30 231 995679. E-mail address: [email protected] (N. Charalambakis).

0749-6419/$ - see front matter
2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijplas.2005.01.002

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1971

1. Introduction In the context of theories applicable to experiments on the necking under uniaxial and biaxial loading of thin strips cut from rolled sheet (Hill, 1948, 1952; Rudnicki and Rice, 1975; Hill and Hutchinson, 1975; Storen and Rice, 1975), the onset of localized necking is explainable as a bifurcation from a state of uniform deformation, without the assumption of pre-existing macroscopic inhomogeneities, such as those due to a position-dependent strain hardening or softening of the material. The neck forms normal to the largest principal strain for biaxial stretching, or along the direction of zero extension for negative ratio of principal strains (Storen and Rice, 1975). Since the assumption of initial imperfections is necessary for the prediction of necking under biaxial stretching, the MK defect theory of Marciniak and Kuczynki (1967) has been developed. A second, equally popular theory, based on linearized stability analysis, was proposed by Rice (1976), Bai (1982), Anand et al. (1987). Important perturbation techniques have been developed by Molinari (1985), Dudzinski and Molinari (1991) and Toth et al. (1996) to deal with the problem of predicting strain localization and forming limit diagrams. The sheet is assumed under homogeneous deformation caused by proportional straining rate at the boundaries. A perturbation is superimposed to the reference homogeneous solution. It is well known that instability or stability is characterized by the fact that the perturbation is increasing or decreasing. The new ingredient of these techniques relies on the fact that they include an analysis of the rate of growth of the perturbations, providing the link between the bifurcation results of Hill (1952), Storen and Rice (1975) and the defect theory of Marciniak and Kuczynki (1967). This analysis, herein named effective instability analysis, characterizes the overall strains for which a certain intensity in the instability growth is developed. The main advantage of the perturbation analysis relies on the fact that it is suitable for more general constitutive laws, such as rate-dependent and thermal softened materials. It is worth noticing that Barbier et al. (1998) proved that the bifurcation analysis for a non-viscous material should be considered as a limit case of the linear perturbation analysis. In this paper, the limits of sheet metal ductility caused by macroscopic material heterogeneities, are discussed. Material heterogeneities may occur, for instance, as a result of the specific material structure or incomplete production or phase transition. The role of spatially variable strain hardening or softening and strain rate sensitivity in the onset and evolution of strain non-uniformities is not well elucidated by the perturbation or the bifurcation analysis. Hill and Hutchinson (1992) proposed a constitutive framework to take account of progressive changes in a stress-based criterion of yielding under biaxial loading. A new concept of engineering the materialÕs microstructure is related to the functionally graded materials (FGMs), wherein the microstructural details are spatially varied through non-uniform distribution of the reinforcement phases. In Fig. 1 (Hirai, 1996; Aboudi et al., 1999), we see an example of continuously graded microstructure. The result is a microstructure that produces continuously changing mechanical properties at the microscopic level. The complex functionality of

1972

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

Fig. 1. Continuously graded material.

materials, such as aerospace engines, circuit boards, artificial tissues, ceramic-tometal joints, includes a continuous change from one structure or composition to another (Pompe et al., 2003). By spatially varying the microstructure, the material can be tailored to yield optimal mechanical behavior, characterized by improved fatigue resistance and interlaminar stress reduction. Two commonly used methods for modeling functionally graded materials are the uncoupled approach, in which spatially varying material properties are employed, and the coupled approach, in which the material microstructure is explicitly taken into account (Collier et al., 2002). Recently, Chaboche et al. (2004) considered mean-field approximation methods in order to transfer any arbitrary local plasticity constitutive description to an overall analytical continuum description. The numerical solution is based on unit cell FE analysis in the context of periodic homogenization and reference solutions for the overall stress–strain responses are presented. Functionally graded materials having continuously varying properties have been studied by Ozturk and Erdogan (1993, 1995), Jin and Batra (1996, 1998), Wu et al. (1999, 2000), Cheng and Batra (2000), Vel and Batra (2002, 2003), Berezovski et al. (2003), Qian et al. (2003), Batra and Love (2004), Zhang and Paulino (2005). Different issues concerning the application of effective properties in the analysis of heterogeneous materials have also been examined from the micromechanical point of view (Finot et al., 1994; Dolbow and Nadeau, 2002; Aboudi et al., 2003; Majta and Zurek, 2003). The problem of shear linear stability of viscoplastic stratified composites was analyzed in Georgievskii (2000). Recently, a damage model for metal-matrix composite was proposed in Naboulsi and Palazotto (2001, 2003) while in Charalambakis and Baxevanis (2004) the adiabatic shearing of non-homogeneous thermoviscoplastic materials was studied. Havner (2004) considered the inhomogeneity

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1973

represented by a single modulus varying slightly with position, with the same stress– strain curve throughout. Ostoja-Starzewski (2004) studied the scale effects in plasticity of random media. In this paper, we adopt the concept of a non-homogeneous model in which the microstructure produces continuously changing mechanical properties at the macroscopic level. The main difficulty is to modelize the macroscopic properties, such as the strain rate sensitivity and the strain hardening or softening, in terms of functions of the spatial variables. In many works, the problem of thermoelastic crack has been analyzed, with the assumption of an exponential spatial variation of the elastic modulus and the thermal properties. Local macroscopic properties, within the functionally graded materials, were first obtained through homogenization techniques based on a chosen scheme and then subsequently used in a global analysis (Hirai, 1996; Aboudi et al., 1999). Recently, Yoshida et al. (2003) presented the experimentation, fundamentals and technique of identification of the mechanical properties of individual component layers by minimizing the difference between the experimental results and the corresponding results of numerical simulation. The method concerns two types of constitutive models (Chaboche–Rousselier and Prager). It is worth noticing that the method is based on experimental data (tensile load versus strain curve in uniaxial tension) for a whole bimetallic sheet but not for individual component layers. Moreover Harth et al. (2004) proposed stochastic methods to describe the influence of scattering test data on the identification of material parameters for the viscoplastic constitutive model of Chan, Bodner, Lindholm in its uniaxial form. Based on existing homogenization relations between the local state distribution function and the elastic yield properties, Adams et al. (2004) presented allowable combinations of properties. A crucial assumption in this work is that the above macroscopic properties can be expressed as smooth functions of the position, approximating the given properties of the material components. This assumption, leading to continuous stress between the material components, is justified only under zero stress joining methods, excluding tensile stress concentration and debonding of the joint. This is achieved by gradually changing the volume fraction of the constituent materials usually only in one direction to obtain a smooth variation of material properties (see Fig. 1 and also Cheng and Batra (2000)). In the above context, we consider a rigid-viscoplastic material with strain and strain rate sensitivity approximated by smooth functions of the spatial variables. The perturbation approach shows the conditions under which an unstable behavior occurs. The instability for homogeneous materials is activated for a uniformly strained sheet under plane stress. However, for a heterogeneous material, the concept of a uniform reference strain is no more applicable. In this paper, we present a non-linear analysis formulation in order to follow the spatial distribution and evolution with time of strain non-uniformities, caused by material heterogeneities, such as position dependent strain hardening or softening and strain-rate sensitivity (Section 2). We focus on the interplay between the necking effect caused by initial imperfections and the strain non-uniformities induced by the above inhomogeneities. In Section 3, the effective instability analysis, proposed by Dudzinski and Molinari (1991) for a homogeneous material, is presented for the case

1974

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

of a material with gradually changing properties only in one direction. The analysis is based on the perturbation imposed on an appropriate non-uniform reference solution, which is characterized by non-uniform strain rate, stress components and thickness in the direction of inhomogeneity. This reference solution can be obtained numerically. The non-homogeneous effective instability method permits to study the critical conditions and select the localization modes by studying the overall strains for which a certain level of instability growth is developed at a material zone. The optimization of the critical orientation is based not only on the measure of the mode growth, but also on the ‘‘deformation path’’ at the material point considered. Finally, we present numerical examples and compare with the results of Dudzinski and Molinari (1991) for the homogeneous case.

2. Formulation of the problem and non-uniformity evolution equations We consider a thin rectangular plate under biaxial loading in the midsurface at remote boundaries. The plate is symmetric about the midsurface Oxy (Fig. 2), with material properties independent of z, where the axis Oz is perpendicular to Oxy, and continuously varying along any direction in Oxy. A special case of this model is a plate made of a functionally graded material in which the volume fraction of the constituent materials change only in the y-direction. It seems that the effect of the gradation on the stress distribution, combined with the joining techniques, may be important (Saraev and Schmauder, 2003; Berezovski et al., 2003). In the sequel, we assume that the smoothness assumption on the material parameters and an appropriate method of joining dissimilar materials with very low residual stresses result to the continuity of traction/strain/displacements across the joints and eliminate the through thickness stress, since, as large strains are involved at the limiting stage, elasticity is neglected. Therefore, we assume that the plate is subjected to plane stress conditions. Moreover, it is assumed Von Mises isotropic behavior.

Fig. 2. Coordinate systems.

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1975

We denote by rxx(x,y,t), ryy(x,y,t) and rxy(x,y,t) the stress components at the point M(x,y) of the plate with respect to the fixed system Oxy (Fig. 2). Neglecting the inertia effects, the stress components satisfy the equilibrium conditions: oðrxx hÞ oðrxy hÞ þ ¼ 0; ox oy

ð1Þ

oðrxy hÞ oðryy hÞ þ ¼ 0; ox oy

ð2Þ

where h(x,y,t) the thickness of the plate. We denote by r11, r22, r12 the stress components of the material point M(x,y) with respect to the system Ox1x2 (Fig. 2), dependent on the angle W. The material behavior is stated to be rigid viscoplastic, obeying the flow law e_ 11 ¼ k_

or ; or11

e_ 22 ¼ k_

or ; or22

e_ 12 ¼ k_

or ; or12

ð3Þ

where e_ 11 ; e_ 22 and e_ 12 the strain-rates at M in the sheet plane, related to the forth component by the incompressibility condition e_ zz ¼ e_ 33 ¼
_e11
e_ 22 :

ð4Þ

The effective stress r is given by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ r211 þ r222
r11 r22 þ 3r212

ð5Þ

and it is work-conjugated to the equivalent strain-rate e_ , where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e_ ¼ ð2=3Þð_e211 þ e_ 222 þ e_ 233 þ 2_e212 Þ:

ð6Þ

Using (5) and (3) 2r11
r22 ; e_ 11 ¼ k_ 2r

2r22
r11 e_ 22 ¼ k_ ; 2r

3r12 e_ 12 ¼ k_ : 2r

ð7Þ

Substituting (7) and (4) into (6) we obtain e_ ¼ k_

ð8Þ

and consequently (3) can be written in the form e_ 11 ¼ e_

or ; or11

e_ 22 ¼ e_

or ; or22

e_ 12 ¼ e_

or or12

ð9Þ

or, recalling (5) e_ 11 ¼ D1 e_ ;

e_ 22 ¼ D2 e_ ;

e_ 12 ¼ D3 e_ ;

ð10Þ

where D1 ¼

2r11
r22 ; 2r

D2 ¼

2r22
r11 ; 2r

D3 ¼

3r12 : 2r

ð11Þ

1976

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

The initial thickness of the plate is h0 = h(x,y,t0) and the initial effective strain and effective strain-rate are, respectively, e0 ¼ eðx; y; t0 Þ and e_ 0 ¼ e_ ðx; y; t0 Þ, both positive functions. A difficult question is how to quantify the non-homogeneous constitutive behavior of the material, needed for this formulation, i.e., how to define the physical grounds in order to have meaningful material variations (see Aboudi et al., 1999; Collier et al., 2002; Berezovski et al., 2003 and related references in Section 1). In this first step work, we ‘‘solve’’ this problem by smoothly approximating the piecewise continuous parameters of the material components by continuous functions. The smooth functions between the values of the material parameters are supposed to correspond to the transition phases of Fig. 1. Then the material is rigid-viscoplastic, exhibiting strain softening (or hardening) characterized by the function n(x,y) < 0 (or n(x,y) > 0) and strain-rate sensitivity m(x,y) > 0. The softening (or hardening) law is described by the following empirical law, adopted by Molinari and Clifton (1983) r ¼ Gðx; yÞenðx;yÞ e_ mðx;yÞ ;

ð12Þ

where
1 < n < 1, 0 < m < 1, for every x,y. G is a material function and e is the effective strain, given by Z t eðx; y; tÞ ¼ e0 ðx; yÞ þ e_ ðx; y; sÞ ds: ð13Þ t0

Using the plane stress assumption, the compatibility condition reads o2 e_ 11 o2 e_ 22 o2 e_ 12 þ 2 ¼2 : 2 ox2 ox1 ox1 ox2

ð14Þ

The current thickness h is related to its initial value by hðx; y; tÞ ¼ h0 ðx; yÞeezz ðx;y;tÞ ;

ð15Þ

where ezz is the logarithmic deformation in the direction z, orthogonal to the plate Z t ð16Þ ezz ðx; y; tÞ ¼ e0zz ðx; yÞ þ e_ zz ðx; y; sÞ ds t0

and e_ zz is given by (4). Differentiating (15) with respect to x or y, we find expressions for the necking gradient in the directions x or y, respectively, in terms of the initial defects and the corresponding strain non-uniformities, for instance   oh oh0 ezz oe11 oe22 ¼ e
h þ : ð17Þ ox ox ox ox Finally, the equilibrium equations with respect to Ox1x2 read: of1 of3 þ ¼ 0; ox1 ox2

ð18Þ

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

of3 of2 þ ¼ 0; ox1 ox2

1977

ð19Þ

where f1 ðx1 ; x2 ; tÞ ¼ r11 h;

ð20Þ

f2 ðx1 ; x2 ; tÞ ¼ r22 h;

ð21Þ

f3 ðx1 ; x2 ; tÞ ¼ r12 h:

ð22Þ

We note that (6) is equivalent to (10) and (5), so we will omit it in the sequel and e_ will be taken as independent variable. Using (5) r ¼ gðx1 ; x2 ; tÞh
1 ;

ð23Þ

where g¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi f12 þ f22
f1 f2 þ 3f32 :

ð24Þ

Recalling (11), we verify that D1 ¼

2f 1
f2 ; 2g

ð25Þ

D2 ¼

2f 2
f1 ; 2g

ð26Þ

D3 ¼

3f 3 : 2g

ð27Þ

Then, by (4) e_ zz ¼
ðD1 þ D2 Þ_e:

ð28Þ

Finally, we note that the stress components with respect to Oxy and Ox1x2 are related by r11 ¼ rxx cos2 W þ ryy sin2 W þ 2rxy sin W cos W;

ð29Þ

r22 ¼ rxx sin2 W þ ryy cos2 W
2rxy sin W cos W;

ð30Þ

r12 ¼ ðryy
rxx Þ sin W cos W þ rxy ðcos2 W
sin2 WÞ:

ð31Þ

We now use (23) to write (12) in the form gðx; y; tÞh
1 ¼ Gðx; yÞenðx;yÞ e_ mðx;yÞ ; from which we obtain by integration over time 1=mðx;yÞ Z t pðx; yÞ gðx; y; sÞ eðx; y; tÞpðx;yÞ ¼ e0 ðx; yÞpðx;yÞ þ ds; Gðx; yÞ1=mðx;yÞ t0 hðx; y; sÞ

ð32Þ

ð33Þ

1978

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

where p¼

mþn : m

ð34Þ

As it is obviously expected, from (33) we obtain immediately that the effective strain and, consequently, the effective strain-rate are non-uniform, if the mechanical parameters are position dependent, even for initially uniform strain. Therefore, using (10), we conclude that the strain-rate components are also non-uniform. The same conclusion of position dependence is valid for the thickness (see (15)) and the stress components (see (1) and (2)). Moreover, since gh
1 > 0, the large time behavior of e at every material zone depends on the sign of m(x,y) + n(x,y). If p(x,y) = (m(x,y) + n(x,y))/m(x,y) < 0 for some (x,y) of the plate, then by (33) there exists a critical time at which the strain at these points becomes infinite. For such a material point (x,y), the critical time tcr is given by 1=mðx;yÞ Z t gðx; y; sÞ e0 ðx; yÞpðx;yÞ Gðx; yÞ1=mðx;yÞ : ds ¼ hðx; y; sÞ jpðx; yÞj t0

ð35Þ

The critical time for the whole plate is obtained as the infimum of all tcr. In the sequel, we assume that for a material exhibiting strain softening, the critical time is not yet reached. For material functions m, n, G and initial data h0, e0, e_0 having physical meaning, this assumption is always valid. This can be explained by the fact that only a finite (and <1) value of strain is needed to cause shear banding. We now use (12), (23) and (15) to write Gh0 en e_ m ¼ ge
ezz ; from which

ð36Þ 

g n ln e þ m ln e_ ¼ ln Gh0


ezz :

Differentiating (37) with respect to time and using (28)  
1 2 n oe oe o e 1 og oe þm þ ðD1 þ D2 Þ ¼ 2 e ot ot ot g ot ot

ð37Þ

ð38Þ

or   2 o2 eðx; y; tÞ nðx; yÞ oeðx; y; tÞ mðx; yÞ
D1 ðx; y; tÞ þ D2 ðx; y; tÞ
ot2 eðx; y; tÞ ot 1 ogðx; y; tÞ oeðx; y; tÞ ¼ 0:
gðx; y; tÞ ot ot

ð39Þ

Taking into account that e0(x,y) > 0, e_ 0 ðx; yÞ > 0, and for strain path satisfying
0:5 6 e_ yy =_exx 6 0, at every point, we apply elementary maximum principle to prove that, if n < 0 or if n > 0 and e0 > n/(D1 + D2), then o2e/ot2 > 0 and e is monotone increasing.

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

Viewing (39) as a BernoulliÕs equation,  2 o2 e oe 1 og oe ¼ 0;
m 2
B ot ot g ot ot

1979

ð40Þ

valid for every x,y, with n B ¼ D1 þ D2
; e

ð41Þ

we obtain the following representation for the strain rate  Z t 1=mðx;yÞ 1=mðx;yÞ
1
1 e_ ðx; y; tÞ ¼ gðx; y; tÞ gðx; y; tÞ ðD1 ðx; y; sÞ e_ 0 ðx; yÞ
m t0

þD2 ðx; y; sÞ
nðx; yÞ=eðx; y; sÞÞgðx; y; tÞ1=mðx;yÞ ds


1 :

ð42Þ

There are two cases in which strain-rate becomes infinite at a finite time tcr: the first corresponds to softened material zone with n(x,y) < 0 and tcr is defined by  Z t jnðx; yÞj gðx; y; tÞmðx; yÞ 1=mðx;yÞ : D1 ðx; y; sÞ þ D2 ðx; y; sÞ
ds ¼ gðx; y; tÞ eðx; y; sÞ e_ 0 ðx; yÞ t0 ð43Þ The second corresponds to material zones exhibiting strain hardening and sufficiently large accumulated effective strain and the critical strain can be obtained in a similar manner. We note that, since only a finite value of the strain rate is needed to create a real shear band, the critical value of time tcr is always unreachable. In any case, we assume that t < tcr. 3. Non-uniform reference solution of continuously graded, along one direction, materials and instability analysis In this section, we first present a non-uniform reference solution appropriate for the linearized instability analysis. For simplicity, we assume that the material exhibits macroscopic inhomogeneity with respect to y-axis. Then a possible reference solution is a solution of the system (10), (32), (14), (18), (19) and (23) (under suitable boundary conditions), which depends only on y and t at most. The deformation is caused by a normal distributed load q(t), parallel to the y-axis, a tangential distributed load r(t) and a uniform longitudinal strain-rate v(t), parallel to the x-axis (Fig. 3). We will verify that there exists a non-uniform time-dependent reference solution S 0inh ¼ f_e011 ; e_ 022 ; e_ 012 ; e_ 0 ; r011 ; r022 ; r012 ; r0 ; h0 g;

ð44Þ

defined by e_ 0xx ¼ e_ 0xx ðx; tÞ  vðtÞ;

ð45Þ

e_ 0yy ¼ e_ 0yy ðy; tÞ;

ð46Þ

1980

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

Fig. 3. Biaxial loading of non-homogeneous plate.

r0xx ðy; tÞh0 ðy; tÞ ¼ pðy; tÞ;

ð47Þ

r0yy ðy; tÞh0 ðy; tÞ ¼ qðtÞ;

ð48Þ

r0xy ðy; tÞh0 ðy; tÞ ¼ rðtÞ

ð49Þ

and (29), (30), (31), (4), (23), (24), (10), 12,13,14,15, (18), (19). The functions v(t), q(t) and r(t) are defined by the boundary conditions, while p(y,t) is an unknown function which will be defined in the sequel. Using (29)–(31) and (20)–(22) we obtain: f10 ¼ r011 h0 ¼ pðy; tÞcos2 W þ qðtÞsin2 W þ 2rðtÞ sin W cos W;

ð50Þ

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1981

f20 ¼ r022 h0 ¼ pðy; tÞsin2 W þ qðtÞcos2 W
2rðtÞ sin W cos W;

ð51Þ

f30 ¼ r012 h0 ¼ ðqðtÞ
pðy; tÞÞ sin W cos W þ rðtÞðcos2 W
sin2 WÞ:

ð52Þ

Since x ¼ x1 cos W
x2 sin W; y ¼ x1 sin W þ x2 cos W; we can use the rule o=ox1 ¼ cos Wo=ox þ sin Wo=oy; o=ox2 ¼
sin Wo=ox þ cos Wo=oy; to see that f10 ; f20 ; f30 satisfy (18) and (19). We now find the strain-rate components with respect to the system Ox1x2: e_ 011 ¼ e_ 0xx cos2 W þ e_ 0yy sin2 W þ 2_e0xy sin W cos W;

ð53Þ

e_ 022 ¼ e_ 0xx sin2 W þ e_ 0yy cos2 W
2_e0xy sin W cos W;

ð54Þ

e_ 012 ¼ ð_e0yy
e_ 0xx Þ sin W cos W þ e_ 0xy ðcos2 W
sin2 WÞ;

ð55Þ

from which e_ 011 ¼ vðtÞcos2 W þ e_ 0yy ðy; tÞsin2 W þ 2_e0xy ðy; tÞ sin W cos W;

ð56Þ

e_ 022 ¼ vðtÞsin2 W þ e_ 0yy ðy; tÞcos2 W
2_e0xy ðy; tÞ sin W cos W;

ð57Þ

e_ 012 ¼ ð_e0yy ðy; tÞ
vðtÞÞ sin W cos W þ e_ 0xy ðy; tÞðcos2 W
sin2 WÞ:

ð58Þ

e_ 011 ;

e_ 022 ;

e_ 012

satisfy the compatibility equaUsing (45) and (46), we easily verify that tion (14). We pass to the derivation of all functions in terms of the known functions q(t), r(t), v(t) and the unknown function p(y,t). We apply (10) with respect to the Oxy system to find: e_ 0 ðy; tÞ ¼ vðtÞ=Dx ðy; tÞ;

ð59Þ

e_ 0yy ðy; tÞ ¼ vðtÞDy ðy; tÞ=Dx ðy; tÞ;

ð60Þ

e_ 0xy ðy; tÞ ¼ vðtÞDxy ðy; tÞ=Dx ðy; tÞ;

ð61Þ

where Dx ðy; tÞ ¼ ð2pðy; tÞ
qðtÞÞ=ð2gðy; tÞÞ;

ð62Þ

Dy ðy; tÞ ¼ ð2qðtÞ
pðy; tÞÞ=ð2gðy; tÞÞ;

ð63Þ

Dxy ðy; tÞ ¼ 3rðtÞ=ð2gðy; tÞÞ;

ð64Þ

1982

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

gðy; tÞ ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðpðy; tÞÞ2 þ ðqðtÞÞ2
pðy; tÞqðtÞ þ 3ðrðtÞÞ2 :

ð65Þ

Then, recalling (56)–(58), the strain-rate components at any point with respect to the Ox1x2 system can be expressed only in terms of the unknown function p(y,t), the given control functions q(t), r(t), v(t) and the orientation W. Moreover, by (4) e_ 033 ¼
ð1 þ Dy =Dx ÞvðtÞ;

ð66Þ

from which, integrating over time Z t 0 0 e33 ðy; tÞ ¼ e33 ðy; t0 Þ
½1 þ Dy ðy; sÞ=Dx ðy; sÞvðsÞ ds:

ð67Þ

t0

Finally, integrating (59) over time Z t vðsÞ=Dx ðy; sÞ ds: e0 ðy; tÞ ¼ e0 ðy; t0 Þ þ

ð68Þ

t0

The actual thickness h0 is related to its initial value by 0

h0 ðy; tÞ ¼ h0 ðy; t0 Þee33 ðy;tÞ ;

ð69Þ

where e33 is given by (67). It remains to express the effective stress r, given by (23) and (12), as a function of p(y,t) and the known boundary conditions. Combining (12), (23) and (67) 0

gðy; tÞ ¼ h0 ðy; t0 Þee33 ðy;tÞ GðyÞe0 ðy; tÞ

nðyÞ 0

mðyÞ e_ ðy; tÞ ;

ð70Þ

from which ln gðy; tÞ ¼ ln½h0 ðy; t0 ÞGðyÞ þ e033 ðy; tÞ þ nðyÞ ln e0 ðy; tÞ þ mðyÞ ln e_ 0 ðy; tÞ;

ð71Þ

where g(y,t), e033 ðy; tÞ, e0(y,t), e_ 0 ðy; tÞ are all functions of the unknown function p(y,t). Eq. (71) is valid at every point y and gives the value of p(y,t) at every time t. We note that (71) gives, for t = t0, the following compatibility condition between the initial data, the material functions and the control parameters at t0 ln gðy; t0 Þ ¼ ln½h0 ðy; t0 ÞGðyÞ þ e033 ðy; t0 Þ þ nðyÞ ln e0 ðy; t0 Þ þ mðyÞ ln e_ 0 ðy; t0 Þ: ð72Þ By the numerical solution of (71) we obtain p(y,t), thus all related deformation and stress functions. The instability analysis, based on the above reference solution, is now divided in two parts. In the first part, we obtain from (1), (2), (5), (10), 12), (14) the corresponding linearized system by putting S ¼ S 0inh þ dS;

ð73Þ

where S 0inh is given by (44) and dS ¼ dS 0 egðt
t0 Þ einx1 ;

ð74Þ

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1983

denotes a perturbation in Ox1-direction causing a possible instability in Ox2-direction if Re(g) > 0. In the sequel, by g we will denote the real part of g. The wave number n is chosen so that dS = 0 at the boundaries. The linearized system is presented in Appendix A, (A.1)–(A.3) and the characteristic equation, giving the values of g, for which a non-trivial solution exists, is of the form AX 2 þ BX þ C ¼ 0;

ð75Þ

where A, B, C are given by (A.4)–(A.6) of the Appendix A, respectively, and X ¼

g : e_ 0

ð76Þ

It is important to note that (75) is associated to each material point and that A, B, C and consequently, X are influenced not only by the position dependent reference solution (44), but also by the derivatives of h0 ; r011 ; r022 ; r012 with respect to Ox1 and Ox2 directions (see (A.4)–(A.6)). More specifically, we are concerned with: (a) the dependence of g on the strain path at the critical points; (b) the variation of the critical angle with position; (c) the activation of multiple necking following material inhomogeneities; (d) the relative importance of initial defects and material inhomogeneities to the necking formation. All that are illustrated in Section 4. The second part of the instability analysis is devoted to the application of effective instability technique of Dudzinski and Molinari (1991), adapted to the non-homogeneous case. Our aim is (a) to exhibit the dependence of g on the boundary deformation exx for different loads f2 at the discretized points, (b) to select the dominant instability modes by minimizing exx with respect to W at these points, and finally (c) to give the (exx
eyy) curve at the critical points for different prescribed values of g. We note that in the whole analysis, always the perturbation of the reference solution is carried out. However, once the point of g = 0 is passed, the perturbations grow and make an effect. Their effect on the later deformation is completely neglected in this work, the time is just increased and new perturbations are applied.

4. Numerical example We assume that the plate is made of three rigid-plastic metallic phases with a periodic continuously graded microstructure along the axis y (Fig. 3). Phase I and phase II are made of materials exhibiting strain hardening, while phase III is made of a material exhibiting strain softening. In Fig. 4, we see the type of functionally graded microstructure considered: phase I, phase I with inclusions of phase II, transition region, phase II with inclusions of phase I, phase II, phase II with inclusions of phase III, transition region, phase III with inclusions of phase II, phase III, etc. The mechanical parameters m, n and G can be approximated by periodic functions of the spatial variable y (Fig. 3) with a period of 2 cm. In phase I, the values of m, n and G are 0.015, 0.1 and 270 MPa, respectively, while in phases II and III, the maximum absolute deviations from the above values are 0.01, 0.15 and 130 MPa,

1984

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

Fig. 4. Type of the non-homogeneous material considered.

respectively (Fig. 3). The material parameters of the transition phases are supposed to be expressed by the smooth connections of the above values. The thickness of the plate is also a function of y and its initial value is uniform and equal to 1 cm, except at the points 0.9 cm 6 y 6 1.1 cm, where a small defect of 0.2 mm (Fig. 3) exists. The deformation is caused by a time increasing and uniformly distributed load q(t) = 0.1 + 100t M N/m, parallel to y-axis and a constant uniform strain rate t(t) ” 10 s
1. Initially, the effective strain is uniform and equal to e0 = 0.05. The numerical solution of the non-linear equation (71) with respect to p(y,t) is obtained with the help of a FORTRAN code for personal computers, constructed by the authors. Eqs. (71) can be written in the form Z t  Z t UðpÞ ¼
ln ½KðpÞ þ ln ½h0 G
MðpÞ ds þ e0zz þ n ln N ðpÞ ds þ e0 0

0

þ m ln ½N ðpÞ ¼ 0;

ð77Þ

where taking into account that r = 0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pþq KðpÞ ¼ p2 þ q2
pq; MðpÞ ¼ tðtÞ; 2p
q

N ðpÞ ¼

2KðpÞ tðtÞ: 2p
q

The derivative of U(p) with respect to p is equal to Rt Z t RðpÞ ds P ðpÞ RðpÞ U0 ðpÞ ¼

; QðpÞ ds þ n R t 0 þm 0 KðpÞ N ðpÞ N ðpÞ ds þ e 0 0

ð78Þ

ð79Þ

where P ðpÞ ¼

2p
q ; 2KðpÞ

QðpÞ ¼


3q 2

ð2p
qÞ 2½P ðpÞð2p
qÞ
KðpÞ RðpÞ ¼ tðtÞ: 2 ð2p
qÞ

tðtÞ; ð80Þ

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1985

Then, at each time step, we use the Newton–Raphson method to compute p. The integrals are calculated by the trapezoidal rule. By the numerical evolution of the non-homogeneous reference solution (44), we verify that the total time of the analysis is 10 ms. At this moment the effective strain rate at the point d1 (Fig. 3) increases rapidly, causing collapse (Fig. 5). It is important to note that the effective instability analysis is conducted in the framework of a linearized analysis. The instability results are valid only when the linearization process is not too strongly developed (Dudzinski and Molinari, 1991). Therefore, the analysis is consistent only for almost linear effective strain-rate and

effective strain rate (sec-1)

14 12 10 8 6 4 2 0 0

2

4

6

8

10

t (msec) Fig. 5. Evolution of the effective strain rate with time at point d1.

0.4 0.3 0.2 0.1 ρ

0 -0.1 -0.2 -0.3 -0.4 -0.5 0

2

4

6

8

t (msec) Fig. 6. Straining-rate path q versus time at the critical point d1.

10

1986

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

straining rate path and for large values of the instability parameter. Regarding Figs. 5–7, we conclude that the analysis at the points d1 is consistent for straining-rate paths with q < 0.17, instability parameters g P 641 and critical time tcr 6 9 ms. These limit values, in the exception of the critical time, may vary from one material phase to another and must be respected in the sequel. All points d1, situated at the center of the regions occupied by phase III, are critical points. We must underline the differences with the homogeneous case. First of all, the straining-rate path q is no more homogeneous. It depends on the spatial

3

3

critical parameter η (10 )

3.5

2.5 2 1.5 1 0.5 0 -0.5 -0.4 -0.3 -0.2 -0.1

0

0.1

0.2

0.3

0.4

ρ Fig. 7. Critical instability parameter versus straining-rate path q at the critical point d1.

critical angle Ψ (degrees)

35 30 25 ρ<0

20 15 10 5 ρ≥0

0 0.5

1

1.5

2

2.5

3

3.5

3

critical parameter η (10 ) Fig. 8. Critical angle W-critical instability parameter g-curve of critical points following the strain path.

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1987

variable, exhibiting different curves with respect to time and varying between different, for each material phase, limit values. In Fig. 6 we see one of these curves, corresponding to the critical points d1. However, in all regions, the curves q versus t are increasing, due to the increasing loading q(t) at the boundaries and the constant boundary strain-rate t. At the considered time, the reference solution gives ‘‘slow’’ changes of thickness, herein named ‘‘reference’’ changes with respect to the y-direction, due to inhomogeneities. These changes coexist with the perturbations dS, which are suddenly applied changes in any direction, superimposed on the ‘‘slow’’ reference changes. In a forthcoming paper, we introduce a relative perturbation analysis, in which, we examine the behavior of the perturbation divided by the corresponding reference solution.

0 10 20 30 40 50 60 70 80 90 angle Ψ (degrees)

(b)

3

parameter η (10 )

ρ=-0.4628

0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1.2 -1.4 -1.6 -1.8

ρ≈-0.4620

c1 c2 c3 c4 c5 c6 c7 c8

parameter η (103)

0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1

angle Ψ (degrees)

ρ≈-0.4622

b1 b2 b3 b4 b5 b6 b7 b8

0 10 20 30 40 50 60 70 80 90 angle Ψ (degrees) 3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1

ρ=-0.4531

d1 d2 d3 d4 d5 d6 d7 d8

ρ=-0.4621 0 10 20 30 40 50 60 70 80 90

0 10 20 30 40 50 60 70 80 90

(d)

0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 -0.7 -0.8 -0.9 -1

(c)

parameter η (103)

parameter η (103)

(a)

(e)

angle Ψ (degrees)

Fig. 9. (a) Initial thickness of the plate with defect. (b) Initial region a. Perturbation measure versus orientation. (c) Initial region b (8 points). (d) Initial region c with increasing mechanical properties (8 points). (e) Initial region d with decreasing mechanical properties (8 points). The ratio of principal strainrates takes very close negative values everywhere. Instability may be activated only in phase III and transition phase.

1988

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

Using (75), we observe that, at the point d1, the curve of the critical instability parameter g (the parameter of the intensity of the perturbation) versus the ratio q ¼ e_ yy =_exx , consists of three regimes, as it can be seen in Fig. 7. At the beginning, the parameter g decreases from 3100 to 2650. After that it has a slight increase and finally, when q becomes positive, it decreases very fast. An interesting fact is that, in the first two regimes,where the ratio of the principal strain-rates at this point is negative, the critical angle decreases from 34 to 0, while it remains constant, equal to zero, in the final stage of biaxial stretching (q > 0) (Fig. 8). To study the instability parameter in a spatial period of the plate, we need to consider four material regions. The first region is occupied by phase I (points a in Fig. 3),

0.05

0.05

0

0

3

-0.05 -0.1 -0.15 -0.2

parameter η (10 )

parameter η (103)

(a)

ρ=-0.3165

-0.25

-0.05 -0.1 -0.15 -0.2 -0.25 ρ≈-0.314 -0.3

-0.3

0 10 20 30 40 50 60 70 80 90

0 10 20 30 40 50 60 70 80 90

(d)

0.1 0 -0.1 -0.2 c1 -0.3 c2 -0.4 c3 c4 -0.5 c5 -0.6 c6 -0.7 c7 -0.8 ρ≈-0.3255 c8 -0.9 0 10 20 30 40 50 60 70 80 90 angle Ψ (degrees)

(c)

angle Ψ (degrees) 3

3

angle Ψ (degrees)

parameter η (10 )

parameter η (103)

(b)

b1 b2 b3 b4 b5 b6 b7 b8

2.5 ρ≈-0.27 2 1.5 1 0.5

d1 d2 d3 d4 d5 d6 d7 d8

0 -0.5

(e)

0 10 20 30 40 50 60 70 80 90 angle Ψ (degrees)

Fig. 10. (a) Thickness of the plate at t = 4 ms. Activation of multiple necking following inhomogeneities. (b) Region a, in stable situation. (c) Region b, in stable situation. (d) Region c, in stable situation. (e) Region d. Coexistence of unstable and stable situations. The critical angles vary following position. The straining-rate path q obtains different negative values in phase III and transition region.

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1989

where the mechanical characteristics remain constant. Next is the region with the small defect b. The last two regions c and d are occupied by phase I with inclusions of phase II, phase II with inclusions of phase I and transition phase, and they are characterized by increasing and decreasing material functions. At the beginning, where the thickness is constant everywhere in the exception of the region with the small defect (Fig. 9(a)), it is obvious from Fig. 9(b)–(e) that only the material regions with the decreasing m, n and G may become unstable at an angle between 34 and 35. The strain path q is slightly increasing in the third region and slightly decreasing in the forth region, taking everywhere negative values.

0 10 20 30 40 50 60 70 80 90 angle Ψ (degrees)

parameter η (103)

(b)

(d)

ρ=-0.11

0.05 0 -0.05 -0.1 -0.15 -0.2 -0.25 -0.3 -0.35 -0.4 -0.45

parameter η (103)

0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0

0.1 0.09 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0

ρ≈-0.105

angle Ψ (degrees) 1.4

c1 c2 c3 c4 c5 c6 c7 ρ=-0.145 c8 0 10 20 30 40 50 60 70 80 90 angle Ψ (degrees)

b1 b2 b3 b4 b5 b6 b7 b8

0 10 20 30 40 50 60 70 80 90

(c)

parameter η (103)

3

parameter η (10 )

(a)

1.2 1 0.8 0.6 0.4

ρ=0.008

d1 d2 d3 d4 d5 d6 d7 d8

0.2 ρ=-0.085 0 0 10 20 30 40 50 60 70 80 90 (e) angle Ψ (degrees)

Fig. 11. (a) Thickness of the plate at t = 9 ms. Reference necking due to inhomogeneities ‘‘competing’’ initial defect. (b) Region a. (c) Region b. (d) Region c still in stable situation. (e) Region d. Large differences between the critical angles at different points for straining-rate path q negative. Very large differences between the values of q in phase III and transition region, where the critical angle W = 0 coincides with the direction of minimum extension.

1990

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

At t = 4 ms, where the thickness presents already fluctuations in the non-homogeneous material regions (Fig. 10(a)), there are points with positive g (Fig. 10(e)). The region d is the most unstable, but now there is a small deviation in the critical angle, which starts to decrease as we approach to the critical point and this is due to larger differences in the ratio q (Fig. 10(e)). It is interesting to note that the region c with increasing mechanical parameters is the most ‘‘safe’’. At the critical time t = 9 ms, there is a large reference necking in region d, competing the effect of the initial defect (Fig. 11(a)). In this region, the instability parameter g

critical parameter η (103)

3.5 3

critical points

2.5 2 1.5 1

transition phase

0.5 0 0

0.02

0.04

0.06

phase I 0.08 0.1

0.12

0.14

εxx Fig. 12. Critical parameter g versus principal strain exx at different material zones.

0.14 0.13 0.12

εxx

0.11 0.1 η*=η/1000

0.09 0.08

η*=0.7 η*=1 η*=1.1

0.07 0.06 0

10

20

30 40 50 60 angle Ψ (degrees)

70

80

90

Fig. 13. exx versus angle W in the transition phase for specific levels of parameter g.

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1991

0.138 0.136 0.134 0.132 εxx

0.13 0.128 η*=η/1000

0.126 0.124 0.122

η*=0.01 η*=0.02 η*=0.03

0.12 0.118

0

10

20

30 40 50 60 angle Ψ (degrees)

70

80

90

Fig. 14. exx versus angle W in phase I for specific levels of parameter g.

0.14

A

0.12

εxx

0.1 0.08 0.06 0.04 0.02

B η*=η/1000 η*=1 η*=1.5 η*=2 η*=2.5 η*=3

D C

0 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02 -0.01 εyy

0

Fig. 15. Limit strain curves. exx versus eyy at point d1 (phase III).

starts to decrease. All material points with positive ratio q have a critical angle equal to 0 (direction of minimum extension) (Fig. 11(e)). Moreover, all material regions, in the exception of region c, are now in unstable situation (Fig. 11(b) and (c)), with larger instability parameters and smaller critical angles. The relation of the critical parameter g with the strain exx for the characteristic points is shown in Fig. 12. In this figure, we see a characteristic plateau at the critical points d1, for values of g between 3100 and 2600, where each level of the instability parameter appears for more than one values of exx. After a certain strain level, the

1992

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

curve drops rapidly. At this characteristic strain level, the ratio q becomes positive and the critical points enter to biaxial stretching. It is also clear that, at the points of the ‘‘safe’’ region a (phase I), the parameter g increases with exx. At the points d5 of transition phase, we see that the behavior is similar with other points of phase I, but at a higher level of instability parameters. In Fig. 13 we see that, in the transition phase, and for all values of g tested, the minimum values of strain are observed for the first time for an angle W = 31. The critical angle decreases from 31 to 13 as g increases. The final necking occurs for an angle of 13. At every reference strain, a necking, in a direction depended on the strain level, is superimposed on the reference thinning, which is always parallel to x-axis, forming a kind of variable double necking through this point. The dynamics of multiple neck formation is related to fragmentation in high rate extension of ductile materials (Guduru and Freund, 2002; Mercier and Molinari, 2003). Our result indicates that material inhomogeneities may also create multiple necking along different orientations in the weak regions. Entering to another material phase, the

Fig. 16. Non-homogeneous plate before ‘‘homogenization’’.

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1993

instability band changes orientation. In phase I, higher levels of g are activated for larger strain levels (Fig. 14). Finally, for every level of the instability parameter g, we consider the process described by a boundary load q(t) = 0.1 + at, and strain rate e_ xx ¼ 10 s
1 , with a taking increasing values from 80 to 200. Our aim is to look for the limit curve (exx
eyy), when g becomes large. For g 6 2500 we see in Fig. 15 that there exists a limit curve, decreasing from A to B as a increases from 80 to 200. For g = 3000 a second branch

0.3 0.25

εxx

0.2 ∧ η=30 ∧ η=40

0.15

∧ η=42.5

0.1 0.05 20

25

30 35 angle Ψ (degrees)

40

45

Fig. 17. Critical angle at the transition points of the non-homogeneous plate.

0.28

εxx

0.275

0.27 ∧ η=0

0.265

0.26

∧ η=0.15 ∧ η=0.2

0

2

4 6 angle Ψ (degrees)

Fig. 18. Necking in the homogenized plate.

8

10

1994

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

appears, very close to the line eyy/exx =
1/2, which now increases from C to D as a increases. Next we compare the behavior of a non-homogeneous material with the behavior of a material ‘‘homogenized’’ in the following sense: its material constants are the mean values of the non-homogeneous parameters. The non-homogeneous material is characterized by the functions of Fig. 16. In phase I, the maximum values of m, n and G are 0.03, 0.6 and 400 MPa, respectively, while in phase II the same values are 0.01,
0.1 and 140 MPa, respectively. The average values of m, n and G are 0.02, 0.25 and 270 MPa, respectively. Initially, the thickness and the effective strain are uniform and equal to 1 cm and 0.05, respectively. In the y-direction, we apply a distributed load q(t) = 0.1 + 10t, while in the x-direction we impose a constant uniform strain rate equal to 10 s
1. The homogeneous plate is made of one material with mechanical constants m = 0.02, n = 0.25 and G = 270 MPa. We impose the same initial and boundary conditions with the non-homogeneous plate. By the numerical solution of the non-homogeneous problem we see that at t = 2.3 ms, the effective strain rate of points belonging to phase II increases rapidly, indicating the failure of the plate. For this period of time, at every point of the homogeneous plate, the straining-rate path varies very slightly between
0.44 and
0.36. In the following, we denote by ^ g ¼ g=_e the instability parameter defined in Dudzinski and Molinari (1991). In Fig. 17, we show the relation between exx and angle W at the transition point K of the non-homogeneous plate (Fig. 16) for various critical parameters ^ g. We observe that, the critical angle for ^g P 40 is between 33 and 33.5. In the homogenized plate,under the same boundary conditions and in the same period of time, we can see only the onset of instability for W = 0 at very high strain levels (Fig. 18). In Dudzinski and Molinari (1991), for the same material parameters n = 0.25, m = 0.02 and q =
0.5, the strain takes its minimum value for ^ g ¼ 0 at W = 0. We note that for ^ g ¼ 0, the minimum value of exx is equal to

100

Ψ=5° Ψ=20° Ψ=25° Ψ=30° Ψ=35°

90 80 70

∧ η

60 50 40 30 20 10 0 0

0.05

0.1

0.15 εxx

0.2

0.25

0.3

Fig. 19. Non-homogeneous material. Envelope of critical angles at the transition points.

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1995

0.25 0.2

∧ η

0.15 0.1 0.05

Ψ=0° Ψ=1° Ψ=2° Ψ=3° Ψ=4°

0 0.25 0.255 0.26 0.265 0.27 0.275 0.28 0.285 0.29 εxx Fig. 20. Homogenized material. Envelope of critical angles.

0.27 in both cases. The curves exx versus b g for various angles at the transition point are illustrated in Fig. 19. By drawing the envelope of these curves, we see that, as exx increases, the critical angle decreases from 35 to 30. Moreover, we see that for W < 25, ^ g is an increasing function of exx. In Fig. 20, concerning the homogeneous case, we see that, as exx increases, the critical angle is constantly equal to zero. The above results confirm that, as it is expected, the ‘‘homogenization’’ based only on the mean values of the mechanical parameters cannot capture the instability modes of the non-homogeneous materials.

5. Conclusions The main features of the instability modes during the biaxial loading of nonhomogeneous viscoplastic materials were presented. The biaxial loading is caused by a uniform time-dependent strain-rate in x-direction combined by a uniform time-increasing tension in y-direction. The effective instability analysis, based on a non-uniform reference solution, may be an efficient tool to study the biaxial loading of continuously graded materials. The different material phases of the sheet ‘‘absorb’’ different straining rate paths during the deformation process and this fact characterizes the overall behavior of the sheet. For every material phase, absolute and effective instability analysis may provide curves relating the critical angle with the critical instability parameter and the strain path, as well as the critical instability parameter versus principal strain. Moreover, by selecting the curves of the principal strain versus critical angle, for given instability parameters, the method provides the optimal orientation of the band and the corresponding limiting strain. An important feature of the instability mode in the transition phases, is the emergence of double necking in two directions under the same strain level. The critical

1996

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

angle at the critical points with the ‘‘worst’’ material parameters (maximum softening and minimum strain-rate sensitivity) tends to zero, as the straining-rate path at these points tends from negative to positive values. Regarding the absolute instability analysis, the instability parameter at the critical points seems to exhibit a non-monotone behavior for negative straining-rate paths q 6 0, while for q > 0 is monotone decreasing, tending to a limit which is estimated by the overall behavior of the reference solution. Concerning the effective instability analysis, several situations seem to coexist at different deformation levels: loss of thinning capacity, interplay with initial effects, multiple necking, due to the increasing effect of material heterogeneities. The multiple necking seems to result exclusively from the non-homogeneous nature of the material, emerging at the weak material zones. Finally, limit principal strain curves are obtained for increasing boundary loading and large instability parameters.

Appendix A The perturbation method, described in Dudzinski and Molinari (1991), can lead to the following homogeneous system: ~  dS~0 ¼ ~ K 0

ðA:1Þ

with 0

dS~ ¼ fdr011 ; dr022 ; dr012 ; dr0 ; d_e011 ; d_e022 ; d_e012 ; d_e0 g; 0

D01

D02

2D03


1

0

ðA:2Þ 0

0

B 0 0 0
X =r0 0 0 0 B B 0 0 0 0 0 0 0 B D11 =r D12 =r 2D13 =r 0
1=_e 0 0 B 0 B D0 =r0 D0 =r0 2D0 =r0 0 0
1=_e 0 22 23 ~ ¼B K B 12 B D013 =r0 D023 =r0 D033 =r0 0 0 0
1=_ e0 B B 0 0 0 0 0 1 0 B B 0 0 0 @ PX 0 hx2 X 0 he Q 1 he Q 1 0 0 0 0 0 hx2 X PX 0 he Q 2 he Q 2 0

0

1

R=_e0 C C C D1 =_e0 C C D2 =_e0 C C C; D3 =_e0 C C 0 C C C 0 A 0 ðA:3Þ

where 2

D011 ¼ 1
ðD01 Þ ;

2

D022 ¼ 1
ðD02 Þ ;

D012 ¼
ð1=2Þ
D01 D02 ; P ¼ inh0 þ oh0 =ox1 ;

D013 ¼
D01 D03 ;

Q1 ¼ inr011 þ W 1 ;

W 1 ¼ or011 =ox1 þ or012 =ox2 ; R ¼ mX þ n=e0 ;

2

D033 ¼ ð3=2Þ
2ðD03 Þ ; D023 ¼
D02 D03 ; Q2 ¼ inr012 þ W 2 ;

W 2 ¼ or012 =ox1 þ or022 =ox2 ;

h0e ¼
h0 =_e0 ;

h0x1 ¼ oh0 =ox1 ;

h0x2 ¼ oh0 =ox2

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1997

and X is given by (76). The system (A.1) has a nontrivial solution if the determinant ~ is equal to zero. This leads to the equation (75), where: of K 2

2

A ¼
2h0x2 ð1
mÞD023 P
ððD02 Þ þ mD022 ÞP 2 þ ðh0x2 Þ ðD012 ð1
mÞ þ 1=2Þ;

ðA:4Þ

 2 B ¼
ðn=e0 Þ P 2 D022 þ ðh0x2 Þ D012
2D023 h0x2 P þ ðh0 =r0 ÞW 2 ð1
mÞ 2

 ð3r011 =ð2r0 ÞÞD03 P þ ðh0 =r0 ÞPW 1 ð3=4
ð1
mÞðD03 Þ Þ þ 3h0x1 r011 ih0 n=ð4r0 Þ
ðh0 =r0 Þh0x2 ð3=4
ð1
mÞðD03 Þ2 ÞW 2
3ðh0 =r0 Þh0x2 D03 ð1
mÞr022 W 1 =ð2r0 Þ

2
ðh0 =2Þh0x2 D03 in ðr011 þ r022 Þ=r0
3n2 ðh0 Þ2 r011 =ð4r0 Þ h i 2 2
ð3h0 =2Þh0x2 D03 inm ðr012 Þ
r011 r022 =ðr0 Þ ;

ðA:5Þ

C ¼ ð3h0 =2Þðn=e0 ÞD03 ðP =r0 Þð
W 2 r011 =r0 þ W 1 r012 =r0 Þ þ ð3h0 =2Þðn=e0 ÞD03 ðh0x2 =r0 Þð
W 2 r012 =r0 þ W 1 r022 =r0 Þ h i 2 2
ð3h0 =2Þðn=e0 ÞD03 h0x2 in ðr012 Þ
r011 r022 =ðr0 Þ :

ðA:6Þ

References Aboudi, J., Pindera, M., Arnold, S.M., 1999. Higher-order theory for functionally graded materials. Composites Part B – Engineering 30 (8), 777–832. Aboudi, J., Pindera, M., Arnold, S.M., 2003. Higher-order theory for periodic multiphase materials with inelastic phases. International Journal of Plasticity 19 (6), 805–847. Adams, B.L., Lyon, M., Henrie, B., 2004. Microstructures by design: linear problems in elastic plastic design. International Journal of Plasticity 20 (8–9), 1577–1602. Anand, L., Kim, K.H., Shawki, T.G., 1987. Onset of shear localization in viscoplastic solids. Journal of the Mechanics and Physics of Solids 35 (4), 407–429. Bai, Y.L., 1982. Thermo-plastic instability in simple shear. Journal of the Mechanics and Physics of Solids 30 (4), 195–207. Barbier, G., Benallal, A., Cano, V., 1998. Relation theorique entre la methode de perturbation linaire et lÕ analyse de bifurcation pour la prediction de la localisation des deformations. Comptes Rendus de L Academie Des Sciences Paris, Series II B 326, 153–158. Batra, R.C., Love, B.M., 2004. Adiabatic shear bands in functionally graded materials. Journal of Thermal Stresses 27 (12), 1101–1123. Berezovski, A., Engelbrecht, J., Maugin, G.A., 2003. Stress wave propagation in functionally graded materials. In: WCU 2003, Paris, pp. 507–509. Chaboche, J.L., Kanout, P., Roos, A., 2004. On the capabilities of meanfield approaches for the description of plasticity in metal-matrix composites. International Journal of Plasticity (in press). Charalambakis, N., Baxevanis, T., 2004. Adiabatic shearing of non-homogeneous thermoviscoplastic materials. International Journal of Plasticity 20, 899–914. Cheng, Z.Q., Batra, R.C., 2000. Three-dimensional thermoelastic deformations of a functionally graded elliptic plate. Composites Part B 31, 97–106. Collier, C.S., et al., 2002. Higher order theory-structural/micro analysis code (hot-smac) software for thermo-mechanical analysis of fgms. Metal Powder Industries Federation.

1998

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

Dolbow, J.E., Nadeau, J.C., 2002. On the use of effective properties for the fracture analysis of microstructured material. Engineering Fracture Mechanics 69, 14–16. Dudzinski, D., Molinari, A., 1991. Perturbation analysis of thermoviscoplastic instabilities in biaxial loading. International Journal of Solids and Structures 27 (5), 601–628. Finot, M., Shen, Y., Needleman, A. et al., 1994. Micromechanical modeling of reinforcement fracture in particle-reinforced metal-matrix composites. Metallurgical and Materials Transactions A – Physical Metallurgy and Materials Science 25 (11), 2403–2420. Georgievskii, D.V., 2000. Viscoplastic stratified composites: shear flows and stability. Computers and Structures 76, 205–210. Guduru, P.R., Freund, L.B., 2002. The dynamics of multiple neck formation and fragmentation in high rate extension of ductile materials. International Journal of Solids and Structures 39, 5615–5632. Harth, T., Schwan, S., Lehn, J., Kollmann, F.G., 2004. Identification of material parameters for inelastic constitutive models: statistical analysis and design of experiments. International Journal of Plasticity 20(8–9), to appear. Havner, K., 2004. On the onset of necking in the tensile test. International Journal of Plasticity 20 (4–5), 965–978. Hill, R., 1948. A theory of the yielding and plastic flow of anisotropic metals. Royal Society of London A. Mathematical and Physical Sciences, 281–297. Hill, R., 1952. Mathematical Theory of Plasticity. Clarendon Press, Oxford. Hill, R., Hutchinson, J.W., 1975. Bifurcation phenomena in the plane tension test. Journal of Mechanics and Physics of Solids 23, 239–264. Hill, R., Hutchinson, J.W., 1992. Differential hardening in sheet metal under biaxial loading – A theoretical framework. Journal of Applied Mechanics – Transactions of the ASME 59 (2). Hirai, T., 1996. Functional gradient materials. In: Brook, R.J. (Ed.), Processing of Ceramics, Pt. 2, vol. 17B. Weinheim, New York, NY, pp. 293–341. Jin, Z.H., Batra, R.C., 1996. Some basic fracture mechanics concepts in functionally graded materials. Journal of Mechanics and Physics of Solids 44, 1221–1235. Jin, Z.H., Batra, R.C., 1998. R-curve and strength behavior of a functionally graded material. Materials Science and Engineering A 242, 70–76. Majta, J., Zurek, A.K., 2003. Modeling of ferrite structure after deformation in the two phase region. International Journal of Plasticity 19 (8), 1097–1118. Marciniak, Z., Kuczynki, K., 1967. Limit strains in the processes of stretch-forming sheet metal. International Journal of Mechanical Sciences 9, 609–620. Mercier, R., Molinari, A., 2003. Predictions of bifurcation and instabilities during dynamic extension. International Journal of Solids and Structures 40, 1995–2016. Molinari, A., 1985. Instabilite´ thermoplastique en cisaillement simple. Journal de Me´canique The´orique et Aplique´e 4, 659–684. Molinari, A., Clifton, R., 1983. Localisation de la de´formation viscoplastique en cisaillement simple, re´sultats exactes en the´orie non-line´aire. Comptes Rendus, Acade´mie des Sciences II 296, 1–4. Naboulsi, S.K., Palazotto, A.N., 2001. Thermomicromechanical damage in composites. AIAA Journal 39 (1), 141–152. Naboulsi, S.K., Palazotto, A.N., 2003. Damage model for metal-matrix composite under high intensity loading. International Journal of Plasticity 19 (4), 435–468. Ostoja-Starzewski, M., 2005. Scale effects in plasticity of random media: status and challenges. International Journal of Plasticity 21 (6), 1119–1160. Ozturk, M., Erdogan, F., 1993. The antisymmetric crack problem in a nonhomogeneous medium. Journal of Applied Mechanics 60, 406–413. Ozturk, M., Erdogan, F., 1995. An axisymmetric crack in bonded materials with a nonhomogeneous interfacial zone under torsion. Journal of Applied Mechanics 62, 116–125. Pompe, W., et al., 2003. Functionally graded materials for biomedical applications. Materials Science and Engineering A 362, 40–60.

G. Chatzigeorgiou, N. Charalambakis / International Journal of Plasticity 21 (2005) 1970–1999

1999

Qian, L.F., Batra, R.C., Chen, L.M., 2004. Analysis of cylindrical bending thermoelastic deformations of functionally graded plates by a meshless local Petrov–Galerkin method. Computational Mechanics 33 (4), 263–273. Rice, J.R., 1976. The localization of plastic deformation. In: Koiter, W. (Ed.), Theoretical and Applied Mechanics, pp. 207–220. Rudnicki, J.W., Rice, J.R., 1975. Conditions for the localization of deformation in pressure-sensitive dilatant materials. Journal of the Mechanics and Physics of Solids 23, 371–394. Saraev, D., Schmauder, S., 2003. Finite element modeling of al/sicp metal matrix composites with particles aligned in stripes – a 2d/3d comparison. International Journal of Plasticity 19 (6), 733–747. Storen, S., Rice, J.R., 1975. Localized necking in thin sheets. Journal of Mechanics and Physics of Solids 23, 421–441. Toth, L.S., Dudzinski, D., Molinari, A., 1996. Forming limit predictions with the perturbation method using stress potential functions of polycrystal viscoplasticity. International Journal of Mechanical Sciences 38 (8–9), 805–824. Vel, S.S., Batra, R.C., 2002. Exact solutions for thermoelastic deformations of functionally graded thick rectangular plates. AIAA Journal 40 (7), 1421–1433. Vel, S.S., Batra, R.C., 2003. Three-dimensional analysis of transient thermal stresses in functionally graded plates. International Journal of Solids and Structures 40, 7181–7196. Wang, B.L., Han, J.C., Du, S., 2000. Cracks problem for non-homogeneous composite material subjected to dynamic loading. International Journal of Solids and Structures 37, 1251–1274. Wu, Y., Ling, Z., Dong, Z., 1999. Stress–strain fields and the effectiveness shear properties for three phase composites with imperfect interface. Solids and Structures 37, 1275–1992. Yoshida, F., Urabe, M., Hino, R., Toropov, V., 2003. Inverse approach to identification of material parameters of cyclic elastoplasticity for component layers of a bimetallic sheet. International Journal of Plasticity 19 (12), 2149–2170. Zhang, Z.J., Paulino, G.H., 2005. Cohesive zone modelling of dynamic failure in homogeneous and functionally graded materials. International Journal of Plasticity 21 (6), 1195–1254.