Limit analysis of composite materials with anisotropic microstructures: A homogenization approach

Mechanics of Materials 43 (2011) 574–585

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

Limit analysis of composite materials with anisotropic microstructures: A homogenization approach H.X. Li ⇑ Department of Civil Engineering, University of Nottingham, University Park, Nottingham NG7 2RD, UK

a r t i c l e

i n f o

Article history: Received 14 December 2010 Received in revised form 17 June 2011 Available online 1 July 2011 Keywords: Limit analysis Composites Homogenization theory Finite element method Nonlinear programming

a b s t r a c t A nonlinear mathematical programming approach together with the finite element method and homogenization technique is developed to implement kinematic limit analysis for a microstructure and the macroscopic strength of a composite with anisotropic constituents can be directly calculated. By means of the homogenization theory, the classical kinematic theorem of limit analysis is generalized to incorporate the microstructure – Representative Volume Element (RVE) chosen from a periodic composite/heterogeneous material. Then, using an associated plastic flow rule, a general yield function is directly introduced into limit analysis and a purely-kinematic formulation is obtained. Based on the mathematical programming technique, the finite element model of microstructure is finally formulated as a nonlinear programming problem subject to only one equality constraint, which is solved by a direct iterative algorithm. The calculation is entirely based on a purely-kinematical velocity field without calculation of stress fields. Meanwhile, only one equality constraint is introduced into the nonlinear programming problem. So the computational cost is very modest. Both anisotropy and pressure-dependence of material yielding behavior are considered in the general form of kinematic limit analysis. The developed method provides a direct approach for determining the macroscopic strength domain of anisotropic composites and can serve as a powerful tool for microstructure design of composites. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Limit analysis is a direct approach to the stability analysis of structures subject to quasi-static loads. Comparing with the traditional incremental elastic–plastic method which involves the calculation of the full history of loaddeformation response to find the limit load of a structure, limit analysis can obtain the collapse load in a direct way through searching the critical failure point in a feasible region. Therefore, limit analysis requires much less computational effort and has more stable numerical calculation. It provides a very rigorous method for the stability analysis of structures. ⇑ Tel.: +44 115 9513938. E-mail addresses: [email protected], [email protected] 0167-6636/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.mechmat.2011.06.007

Limit analysis is based on two dual theorems, the static (lower bound) limit theorem (Hill, 1948) and the kinematic (upper bound) limit theorem (Hill, 1951; Drucker et al., 1952; Drucker, 1953). Due to the complexity of engineering problems, numerical methods are often necessary to implement limit analysis. A linear programming (LP) approach was first proposed (see Koopman and Lance, 1965; Maier et al., 1979; Cristiansen, 1981; Turgeman and Pastor, 1982; Pellegrino, 1988; Ponter and Carter, 1997; Ponter et al., 2000). The fundamental technique of this approach is to discretize the nonlinear yield surface and then use the linearized yield condition in the analysis. Accordingly, some additional constraints and computational error are introduced. In the later, a non-linear programming (NLP) approach was developed and widely used in limit analysis (see Gao, 1988; Zouain et al., 1993; Zhang and Lu, 1994, 1995; Liu et al., 1995; Capsoni and

H.X. Li / Mechanics of Materials 43 (2011) 574–585

Corradi, 1997; Capsoni et al., 2001a,b; Li et al., 2001, 2003; Corradi and Vena, 2003; Li and Yu, 2005, 2006a; Corradi et al., 2006; Makrodimopoulos and Martin, 2006, 2007; Pastor et al., 2003, 2009a). However, the above works are mainly focused on macro-structure analysis. As composite and heterogeneous materials are increasingly used in engineering, limit analysis of micro-structures has drawn much attention in recent years. Based on the homogenization technique, Buhan and Taliercio (1991), Taliercio (1992) and Taliercio and Sagramoso (1995) were among the first to introduce limit analysis into the plastic problem of composite materials at the microscopic scale. Some simplified analytical models were proposed to calculate the limit loads of typical microstructures. Francescato and Pastor (1997) presented a linear programming approach to calculate the macroscopic strength of composites with Tresca’s constituents by means of the homogenization technique and finite element method. A convex optimization (Thai et al., 1998; Trillat and Pastor, 2005; Pastor et al., 2008, 2009b) was recently proposed to apply limit analysis to porous materials with the Gurson model, where the homogenization technique was used for heterogeneous materials. Dallot and Sab (2008a,b) also used the homogenization technique to study the limit state of multi-layered plates. Zhang et al. (2004, 2009) applied the homogenization theory into limit analysis and presented a linear programming approach to implement static limit analysis for composite materials. It should be noted that the above microstructure approaches are presented only for isotropic composite materials. Based on a static direct method with the homogenization technique (Weichert et al., 1999a,b), Chen et al. (2010) presented a 3-dimensional finite element approach for limit analysis of periodic metal-matrix composites. The Wilson non-conforming elements were used to increase the accuracy of finite-element solving of an optimization problem. An isotropic microstructure was analyzed and a lower-bound to the limit load of von Mises’ materials can be directly obtained. Based on the homogenization theory and an advanced nonlinear programming approach (Zhang et al., 1991; Liu et al., 1995; Zhang and Lu, 1995; Chen et al., 1998), Li et al. (2001, 2003) and Li and Yu (2006b) developed a nonlinear programming method to apply the kinematic limit analysis for composite and heterogeneous materials obeying the von Mises and ellipsoid yield criteria, where nonlinear yield surface was directly introduced into limit analysis without discretization. The calculation is based on purely-kinematic microscopic velocity fields. The stress field does not need to be calculated and the macroscopic limit state of a composite can be obtained based on its microstructures. Both isotropic and orthotropic materials can be analyzed at the microscopic scale. Even up to now, anisotropic limit analysis still remains a challenge for macrostructure analysis, and therefore it becomes particularly difficult for microstructure analysis. A general method for limit analysis of anisotropic microstructures is still unavailable. The purpose of this paper is to develop a general nonlinear programming approach for kinematic limit analysis of anisotropic composite or heterogeneous materials. First,

575

by means of homogenization theory, the classical kinematic theorem of limit analysis is extended so that it can be applied to a microstructure such as a Representative Volume Element (RVE) chosen from a periodic composite. Then, using an associated plastic flow rule, a general yield function is directly introduced into limit analysis and a purely-kinematic formulation is obtained. The nonlinear yield surface does not need to be linearized, which significantly reduces the number of constraints and the computational effort. Based on nonlinear mathematical programming theory, the finite element model of the kinematic limit analysis for an anisotropic microstructure is finally formulated as a nonlinear programming problem subject to only one equality constraint. A direct iterative algorithm based on the mixed Lagrangean-penalty method is then proposed to solve the resulting nonlinear programming problem. Numerical results are presented to investigate the effect of anisotropic microstructures on the macroscopic strength domain of a heterogeneous/composite material. 2. Homogenization theory As an innovative micromechanics technique, homogenization theory (Suquet, 1987) was developed in the 1980s and has been widely exploited in recent years. The homogenization technique is mostly applied to a periodic composite or heterogeneous material and the term ‘homogenization’ usually qualifies the passage from the micro- to the macro-scale. The periodic feature of a composite can be represented by the microstructure – Representative Volume Element (RVE) as shown in Fig. 1. Then the correlation between the microstructure and macroscopic properties of a composite can be obtained through study of the RVE. Two well-separated scales have to be considered: a microscopic or local scale (x) which is small enough to characterize the microstructures and a macroscopic or overall scale (y) which is large enough to consider the RVE as a macroscopic point of the composite. At the microscopic scale in the RVE, the local strain field e(x) can be split into two parts: the overall strain E which stands for the actual strain field in the RVE if it is homogeneous; and a correction or fluctuation ~eðxÞ which accounts for the presence of heterogeneities. The overall strain E has the average value of e(x) over the RVE. If e(x) derives from a displacement field u(x), the correction part ~eðxÞ derives ~ ðxÞ. Therefore, the from a fluctuation displacement field u average of the correction part over the whole RVE should vanish. Then, the displacement and strain fields in the RVE can be decomposed as:

~ ðxÞ; uðxÞ ¼ E  x þ u ~ eðxÞ ¼ E þ eðxÞ:

ð1Þ ð2Þ

According to the periodic feature of the fluctuation fields, the average of ~eðxÞ over the RVE should vanish and therefore it can get

h~eðxÞi ¼ 0;

E ¼ heðxÞi

ð3a; bÞ

where the operator hi stands for the volume average of a field over the RVE:

576

H.X. Li / Mechanics of Materials 43 (2011) 574–585

x2 P

x1

y2

RVE

y2

y1

y1

Fig. 1. Periodic microstructures and a RVE of a composite.

hf i ¼

1 V0

Z

f ðxÞdx;

ð4Þ

3.1. Kinematic limit theorem

V

where V0 is the volume of the RVE. Similarly, at the microscopic scale in the RVE, the local stress field r(x) can also be split into two parts: the overall/average stress R which stands for the actual stress field in the RVE if it is homogeneous; and a correction or ~ ðxÞ which accounts for the presence of heterfluctuation r ogeneities. The local stress field r(x) is satisfied with the equilibrium condition throughout the entire periodic medium, i.e. both inside RVE and on the boundary of RVE. Then, the local stress field r(x) can be expressed as:

rðxÞ ¼ R þ r~ ðxÞ; div ðrðxÞÞ ¼ 0; ~ ðxÞi ¼ 0; hr

R ¼ hrðxÞi:

ð5Þ ð6Þ ð7a; bÞ

Based on the above decomposition of the local strain and stress fields in the microstructure, the principle of virtual work for the RVE can then be expressed as:

hrðxÞ : eðxÞi ¼ hrðxÞi : heðxÞi;

ð8Þ

which is equivalent to the following expression:

hrðxÞ : eðxÞi ¼ R : E:

ð9Þ

~ It should be emphasized that the fluctuation terms u ~  n are periodic and anti-periodic respectively on and r the boundary of the RVE, whereas n is the outer normal to the boundary. Instead of giving boundary values, the periodic feature of the fluctuation terms on the boundary of the RVE need to be introduced in the microstructure analysis. This represents an important difference between the homogenization technique and other micromechanics approaches (e.g. the self-consistent method or the Mori– Tanaka method). More details about the homogenization technique can be found in the literature (Suquet, 1987).

The collapse load of a structure under quasi-static loading can be found by means of the kinematic theorem of limit analysis (Hill, 1951; Drucker et al., 1952; Drucker, 1953) which states: among all kinematically admissible velocities, the real one yields the lowest rate of the plastic dissipation power. In terms of the kinematically admissible velocities, the kinematic limit theorem can be expressed as follows:

Z  Z Z k t T u_  dC þ f T u_  dv 6 Dðe_ p Þdv ; Ct

V

ð10Þ

V

where k is the limit load multiplier, t is the reference load of external tractions, f is the reference load of body forces, u_ is the displacement velocity, e_ p is the plastic strain rate, Dðe_ p Þ represents a function for the rate of the plastic dissipation power in terms of the admissible strain rate e_ p , the superscript ‘‘*’’ stands for a parameter corresponding to the kinematically admissible velocity, Ct denotes the traction boundary, and V represents the space domain of the structure. If the body force is omitted and the geometry compatible conditions are introduced, the kinematic limit theorem (10) can be re-expressed as the following formulation

8 R T  R p > < k Ct t u_ dC 6 V Dðe_ Þdv ; p _ in V; s:t: e_ ¼ HðuÞ > : u_ ¼ 0 on Cu ;

ð11Þ

where ‘‘s.t.’’ is the abbreviation of ‘‘subject to’’, Cu denotes _ is a linear differenthe displacement boundary, and HðuÞ tial operator which defines the geometry compatibility condition of deformation. The kinematic limit analysis is to find the minimum, optimized limit multiplier k with kt being the limit load of the structure. 3.2. Anisotropic material model

3. Kinematic limit analysis of anisotropic materials The classical kinematic theorem of limit analysis (Hill, 1951; Drucker et al., 1952; Drucker, 1953) was proposed based on the standard material assumption of perfect plasticity with associated plastic flow, and accordingly an upper bound to the limit load of a structure is defined.

In limit analysis, it is assumed that the deformation is small at incipient collapse and the material can be modelled with sufficient accuracy using rigid-perfect plasticity and an associated plastic flow rule. In the paper, the yield criterion of materials (constituents of a composite) is mathematically expressed by a second-order polynomial

577

H.X. Li / Mechanics of Materials 43 (2011) 574–585

which can be used for most of current yield criteria, e.g. isotropic materials (such as von Mises’s criterion), orthotropic materials (such as Hill’s criterion), pressuredependent materials (such as the Mohr–Coulomb or Drucker–Prager criterion), and anisotropic materials (such as the Tsai-Wu criterion). For the purpose of the simplicity of expression when the finite element method is used, the responses of the body to the load, true stress r, true straine, and displacement u, are represented as column vectors, i.e. in the 3-Dimensional model, r ¼ ½r11 ; r22 ; h pffiffiffi pffiffiffi pffiffiffi pffiffiffi pffiffiffi r33 ; 2r12 ; 2r23 ; 2r13 T ; e ¼ e11 ; e22 ; e33 ; 2e12 ; 2e23 ; pffiffiffi 2e13 T , and u = [u1, u2, u3]T. Then, the anisotropic yield function is defined as follows: T

2

6 6  12 6 2rs 6 6 1 6 2r2s P¼6 6 0 6 6 6 6 0 4 0

FðrÞ ¼ r P r þ r Q  1 ¼ 0;

ð12Þ

where F(r) defines a yield function in terms of material strength parameters, P and Q are coefficient matrices related to the strength properties of the material. The material coefficient matrices P and Q can be determined based on the choice of yield criteria. For example, Hill’s yield criterion is frequently used for orthotropic materials and often expressed as: 2

2

þ 2Lr

þ 2M r

þ 2Nr

2 12

 1 ¼ 0;

8 2H ¼ X12 þ Y12  Z12 ; > > < 2G ¼ Z12 þ X12  Y12 ; > > : 2F ¼ 1 þ 1  1 ; 2L ¼

1

Z2

S1

FðrÞ ¼

GþH 6 H 6 6 6 G P¼6 6 0 6 6 4 0 0

1

 2r1 2

0

0

1

0

0

r2s

ð13Þ

ð14a; b; cÞ

1

; 2

S2

2N ¼

1

ð15a; b; cÞ

; 2

S3

3

H

G

0

0

0

HþF F

F FþG

0 0

0 0

0

0

N

0

0

0

0

L

7 7 7 7 7; 07 7 7 05

0

0

0

0 M

T

Q ¼ ½ 0; 0; 0; 0; 0; 0  :

0 0

ð16aÞ

0

7 0 7 7 7 0 7 7 7; 0 7 7 7 7 0 7 5

s

 2r1 2

r2s

0

0

3 2r2s

0

0

0

0

3 2r2s

0

0

0

0

s

3

ð17Þ

3 2r2s

pffiffiffiffi J 2  c0 ¼ 0;

ð18Þ

pffiffiffiffiffiffiffiffiffiffiffiffi rT P r þ rT Q  1 ¼ 0;

ð19Þ

where

2

where X, Y, Z are the uniaxial yield strength of the orthotropic material in the principal anisotropic directions, and S1, S2, S3 are the shear yield strength with respect to the principal axes of anisotropy. The yield criterion (12) can be used to define Hill’s materials when the coefficient matrices P and Q are chosen as

2

0

1 3c2

6 0 6 1 6 6c20 6 6 1 6  6c2 6 0 P¼6 6 0 6 6 6 0 6 4 0

 6c12

1 6c20

0

0

1 3c20

 6c12

0

0

 6c12

1 3c20

0

0

0

0

1 2c20

0

0

0

0

1 2c20

0

0

0

0

0

0

0

0

Q¼

h

u0 c0

;

u0 c0

;

u0 c0

; 0; 0; 0

iT

3

7 07 7 7 7 07 7 7; 07 7 7 07 7 5

ð20aÞ

1 2c20

X2

2M ¼

; 2

0

s

where I1 is the first invariant of stress tensor, J2 is the second invariant of the deviatoric stress tensor, u0 and c0 are strength parameters of the material. If the yield criterion (18) is expressed in terms of stress components, it will become

2

where F, G, H, L, M and N are the material constants defined by the following equations

Y2

 2r1 2

s

where rs is the uniaxial yield stress of the material. For pressure-dependent materials, the Mohr–Coulomb and Drucker–Prager criteria are often used and expressed as

FðrÞ ¼ Fðr22  r33 Þ þ Gðr33  r11 Þ þ Hðr11  r22 Þ 2 31

 2r1 2

FðI1 ; J 2 Þ ¼ u0 I1 þ

T

2 23

1

r2s

:

ð20bÞ

The yield criterion (19) is actually equivalent to the expression (12). In the similar way, it is not difficult to determine the coefficient matrices P and Q for an anisotropic yield material, e.g. obeying the quadratic Tsai-Wu criterion:

FðrÞ ¼ Aij rij þ Bijkl rij rkl  1 ¼ 0;

ð21Þ

where Aij and Bijkl (i, j, k, l = 1, 2, 3) are the material strength parameters. When the stresses are expressed in column vectors, the Tsai-Wu yield criterion (21) has the exact form of Eq. (12). It is noted that expression of Eq. (12) is a second-order, smooth and continuous polynomial, and is convenient to be used in limit analysis. 3.3. Plastic dissipation power

ð16bÞ

The special case of the coefficient matrix P, as expressed in the following Eq. (17), is for the isotropic von Mises yield criterion.

In order to implement the kinematic limit analysis, the rate of the plastic dissipation power needs to be solved. Limit analysis is based on the material model of associated plastic flow. Therefore, when the anisotropic yield function

578

H.X. Li / Mechanics of Materials 43 (2011) 574–585

(12) with the associated plastic flow rule is used, the plastic strain rate can be determined as follows:

e_ p ¼ l_

@wðrÞ @FðrÞ ¼ l_ ¼ 2l_ P r þ l_ Q ; @r @r

ð22Þ

where w(r) denotes a plastic potential function that resembles the yield function F(r), and l_ is a non-negative plastic proportionality factor. If the associated plastic flow is assumed for a material, w(r) = F(r). Then, the stress at the yield surface can be expressed in terms of the plastic strain-rate as

1

1

r ¼ _ P 1 e_ p  P 1 Q : 2l 2

ð23Þ

The material matrix P is required to positive definite or semi-definite, and then P1 can be uniquely determined when P is non-singular. However, the matrix P may be semi-singular when P is positive semi-definite. For this kind of materials, P1 can be approximately determined as (P + cI)1, where I is the identity matrix of the same order with the matrix P and c is a very small positive real number. Since the stress at the yield surface still needs to be satisfied with the yielding condition of Eq. (12). By introducing Eq. (23) into the yield criterion (12), the non-negative plastic flow factor l_ can then be solved as

l_ ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðe_ p ÞT P 1 e_ p 4 þ Q T P 1 Q

ð24Þ

:

Then, the plastic dissipation power for a material obeying the yield criterion (12) with associated plastic flow can be expressed as follows



T

1 1 p 1 1 P e_  P Q e_ p 2l_ 2 1 p T 1 p 1 p T 1 ¼ ðe_ Þ P e_  ðe_ Þ P Q 2l_ 2 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 1 ð4 þ Q T P 1 Q Þ  ðe_ p ÞT P 1 e_ p ¼ 2 1 p T 1  ðe_ Þ P Q : 2

Dðe_ p Þ ¼ rT e_ p ¼

ð25Þ

The rate of the plastic dissipation power for an anisotropic material is finally expressed in terms of a pure strain-rate field without any stress fields. Once the kinematically admissible velocity is obtained, the plastic dissipation power can then be calculated. Moreover, the nonlinear anisotropic yield surface is not linearized. This can avoid the computational error from the linearization of yield surfaces and also prevent to generate additional constraints. 3.4. Nonlinear programming of kinematic limit analysis The kinematic limit analysis (11) provides a direct calculation of the limit load k t of structures. Based on the mathematical programming theory and the solution of the plastic dissipation power (25), the kinematic limit analysis of the anisotropic material can then be formulated as the following nonlinear programming problem

8  rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ffi R 1 > T 1 > _ p ÞT P 1 e_ p  12 ðe_ p ÞT P 1 Q dv ; > ð4 þ Q k ¼ min P Q Þ  ð e > V 2 > e_ p > < R T _ s:t: Ct t ud C ¼ 1; > > > _ in V; _ p ¼ HðuÞ > e > > : u_ ¼ 0 on Cu ; ð26Þ

where t is the reference load column vector of surface tractions. Finally, kinematic limit analysis of anisotropic materials is formulated as a nonlinear minimum optimization problem subject to equality constraints. An upper bound to the limit load multiplier k can be found by solving Eq. (26), with k  t denoting the limit load of the structure. 4. Limit analysis of microstructures The limit load of a macro-structure can be obtained by solving the nonlinear programming problem (26). In order to extend it for a micro-structure, e.g. a RVE from a composite or heterogeneous material, the homogenization theory is introduced into limit analysis in this paper. 4.1. Limit analysis based on homogenization theory According to the homogenization technique, all local variables, e.g. the stress, strain and displacement fields, are decomposed into two parts: the average macroscopic fields and the fluctuation microscopic fields. Then, the objective function of the nonlinear programming problem (26) can be re-written as  Z  rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   1 1 T T ð4 þ Q T P 1 Q Þ  ðe_ p Þ P 1 e_ p  ðe_ p Þ P 1 Q dv ; 2 V 2  Z  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 T T ¼ ð4 þ Q T P 1 Q Þ  ðð~e_ p þ E_ p Þ P 1 ð~e_ p þ E_ p ÞÞ  ð~e_ p þ E_ p Þ P 1 Q dv ; 2 V 2 Z  rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 1 ð4 þ Q T P 1 Q Þ  ð~e_ p ÞT P 1 ~e_ p þ 2ð~e_ p ÞT P 1 E_ p þ ðE_ p ÞT P 1 E_ p ¼ V 2  1 T  ð~e_ p þ E_ p Þ P 1 Q dv : ð27Þ 2

Meanwhile, the normalization condition in the nonlinear programming problem (26) can be re-expressed as

Z

_ C¼ t T ud

Ct

Z

rT e_ p dv ¼ V 0 ðRT E_ p Þ ¼ 1;

ð28Þ

V

where V0 is the volume of the RVE. Accordingly, the geometry condition and boundary condition are also imposed on the microscopic fluctuation variables. Then, the kinematic limit analysis based on the homogenization theory for a microstructure can be expressed as the following nonlinear programming problem: 2 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 8  ffi 3 T T T T > > R 6 12 ð4 þ Q P 1 Q Þ  ð~e_ p Þ P 1 ~e_ p þ 2ð~e_ p Þ P 1 E_ p þ ðE_ p Þ P 1 E_ p 7 > > > k ¼ min V 4 5dv ; > > ~e_ p ;E_ p > T > <  12 ð~e_ p þ E_ p Þ P 1 Q > s: t: > > > > > > > > :

V 0 ðRT E_ p Þ ¼ 1; ~_ Þ in V; ~e_ p ¼ Hðu ~_ is periodic on C; u

ð29Þ

Finally, the microscopic kinematic limit analysis for a microstructure is formulated as a nonlinear mathematical programming problem subject to equality constraints. An

579

H.X. Li / Mechanics of Materials 43 (2011) 574–585

upper bound to the limit load multiplier k can be obtained by solving Eq. (29), with kR denoting the limit load or macroscopic strength domain of the composite. 4.2. Finite element modelling In order to numerically implement the limit analysis for the RVE, the displacement-based finite element method is used. Based on the finite element technique, the microstructure RVE is discretized by finite elements S V e ðV ¼ Ne¼1 V e Þ. Then, the displacement velocities and strain rates can be interpolated in terms of unknown nodal displacement velocities as:

~_ e ðxÞ ¼ N e ðxÞ~d_ e ; u ~e_ ðxÞ ¼ B ðxÞ~d_ ; e

e

e

K i ¼ C Te ðK e Þi C e ;

ð37Þ

ð31Þ

Gi ¼ C Te ðGe Þi :

ð38Þ

Z  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð4 þ Q T P 1 Q Þ  ðð~e_ p ÞT P 1 ~e_ p þ 2ð~e_ p ÞT P 1 E_ p þ ðE_ p ÞT P 1 E_ p Þ V 2 1  ð~e_ p þ E_ p ÞT P 1 Q dv 2 N Z  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 1 _ _ _ ð4 þ Q T P 1 Q Þ  ððBe d~e ÞT P 1 Be d~e þ 2ðBe d~e ÞT P 1 E_ p þ ðE_ p ÞT P 1 E_ p Þ 2 e¼1 V e

 1  ðBe ~d_ e þ E_ P ÞT P 1 Q dv 2

¼

N X IG X

ðqe Þi jJji

e¼1 i¼1

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð4 þ Q T P 1 Q Þ  ð~d_ Te ðK e Þi ~d_ e þ 2~d_ Te ðGe Þi E_ p þ ðE_ p ÞT P 1 E_ P Þ 2

 1 1  ~d_ Te ðGe Þi Q  ðE_ p ÞT P 1 Q ; 2 2

ð32Þ

where (qe)i is the Gaussian integral weight at the i-th Gaussian integral point in the element e, jJji is the determinant of the Jacobian matrix at the i-th Gaussian integral point and IG is the number of Gaussian integral points in the finite element e, and

K e ¼ BTe P 1 Be ; Ge ¼

BTe P 1 :

where I denotes the set of all Gaussian integration points of the discretized RVE and

ð30Þ

_ where, with reference to the e-th finite element, ~ de is the nodal fluctuation displacement velocity, Ne(x) is the interpolation function matrix and Be(x) is the strain function matrix. By means of the Gaussian integration technique, the objective function in Eq. (29) can be discretized and expressed in terms of the nodal displacement velocity as follows

¼

Z  rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 1 ð4 þ Q T P 1 Q Þ  ð~e_ p ÞT P 1 ~e_ p þ 2ð~e_ p ÞT P 1 E_ p þ ðE_ p ÞT P 1 E_ p V 2  N X 1  ð~e_ p þ E_ p ÞT P 1 Q dv ¼ 2 e¼1  rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi IG   X 1 ð4 þ Q T P 1 Q Þ  d~_ Te ðK e Þi d~_ e þ 2d~_ Te ðGe Þi E_ p þ ðE_ p ÞT P 1 E_ p ðqe Þi jJji 2  X i¼1 1 1  d~_ Te ðGe Þi Q  ðE_ p ÞT P 1 Q ¼ 2rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 i2I    1 ð4 þ Q T P 1 Q Þ  ~d_ T K i ~d_ þ 2~d_ T Gi E_ p þ ðE_ p ÞT P 1 E_ p qi jJji 2  1 1 ð36Þ  ~d_ T Gi Q  ðE_ p ÞT P 1 Q ; 2 2

ð33Þ ð34Þ

Finally, the finite element model of the kinematic limit analysis for a microstructure is expressed as the following discretized nonlinear programming problem: 8 > > > > < k ¼ min > > > > :

1 V ~_ E_ p 0 d;

2 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi 3 T 1 1 ~d_ T K i ~d_ þ 2~d_ T Gi E_ p þ ðE_ p ÞT P 1 E_ p ð4 þ Q P Q Þ  7 2 qi jJji 6 4 5; i2I _  1 d~T Gi Q  1 ðE_ p ÞT P 1 Q

P

2

s:t:

2

RT E_ p ¼ 1:

ð39Þ It should be noted that the periodic feature of the nodal fluctuation displacement velocity ~ d_ on the boundary of the RVE must be enforced when the finite element technique is performed. Then the limit load multiplier k can be obtained by solving the above minimum optimization problem. 5. A direct iterative algorithm The kinematic limit analysis (39) is a nonlinear, minimum optimization problem subject to only one equality constraint. The equality constraint can be introduced into the programming problem by means of the Lagrangean method (Bazaraa et al., 2006; Nocedal and Wright, 2006). And the nonlinear objective function can be solved based on a Newton-type iterative algorithm. However, the nonlinear objective function in Eq. (39) is non-differentiable which may arise in another numerical difficulty. Therefore, a penalty factor is introduced in this paper to overcome the numerical difficulty from non-differentiability. First, based on the mathematical programming theory, the equality constraint of the normalization condition in Eq. (39) is introduced into the optimization problem by means of the Lagrangean method (Bazaraa et al., 2006; Nocedal and Wright, 2006). As a result, an unconstrained minimum optimization problem is obtained as follows

By introducing the transformation matrix of each element Ce, the nodal displacement velocity vector ~ d_ e for each element can be expressed by the global nodal displacement velocity vector ~ d_ for the structure as

2 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  3 _ _ _ T 1 ð4 þ Q T P 1 Q Þ  ~dT K i ~d þ 2~dT Gi E_ p þ ðE_ p Þ P 1 E_ p 7 1 X 2 _~ _ p F Pðd; E ;L Þ ¼ qi jJji 6 4 5 V 0 i2I _ T  12 ~dT Gi Q  12 ðE_ p Þ P 1 Q ð40Þ þ LF ð1  RT E_ p Þ;

_ ~d_ e ¼ C e  ~d:

where LF are the Lagrangean multiplier. Then, according to the Kuhn–Tucker stationarity condition (Himmelblau, 1972; Bazaraa et al., 2006), once apply@P ing @ P~_ ¼ 0; @@EP _ p ¼ 0; @LF ¼ 0 to Eq. (40), the following set of

ð35Þ

Then, the Eq. (32) can be expressed in terms of the global displacement velocity as follows:

@d

580

H.X. Li / Mechanics of Materials 43 (2011) 574–585

equations is obtained to solve the kinematic limit analysis problem (39):

" 8 P > > 1 1 > q jJj > V0 i i 2 > > i2I > < " P 1 1 > q jJj > > > V 0 i2I i i 2 > > > : T P R E_ ¼ 1:

# pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ T 1 P Q ÞðK i ~dþðE_ p ÞT GTi Þ pð4þQ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  12 Q T Gi T ¼ 0; _ ~d_ T G E_ p þðE_ p ÞT P 1 E_ ~d_ T K ~dþ2 i i

# pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ T 1 P Q ÞðP 1 E_ p þGTi ~d T Þ 1 1 pð4þQ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P Q  L F R ¼ 0;  2 _ _ _ T ~dT K ~dþ2~dT G E_ p þðE_ p Þ P 1 E_ i i

ð41a; b; cÞ It is difficult to directly solve the set of Eqs. (41) because it is nonlinear and non-differentiable as well. Then, an iteration technique based on the Newton method is developed to linearize the nonlinear equation set (41) and meanwhile the non-differentiable areas will be distinguished. In order to implement the iteration technique, the set of Eqs. (41) is re-expressed as follows 8   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  P ICP > 1 > > q jJji x2 ð4 þ Q T P 1 Q Þ K i ~d_ þ ðE_ p ÞT GTi  12 Q T GTi ¼ 0; i > V 0 > > < i2I   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  P 1 x ICP ð4 þ Q T P 1 Q Þ P 1 E _ p þ G T ~d_ T  1 P 1 Q  LF R ¼ 0; q jJj > i i i V 2 2 > > > 0 i2I > > : T _p R E ¼ 1; ð42a;b;cÞ

where x

xICP

ICP

is the coefficient parameter which is defined by

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 ~d_ T K i ~d_ þ 2~d_ T Gi E_ p þ ðE_ p ÞT P 1 E_ p ¼

ð43Þ

and the superscript ‘ICP’ indicates that xICP is an Iteration Control Parameter. By solving the linearized equation set (42), the variable ~ d_ and E_ p can be determined. Then the limit load multiplier can be calculated as 2 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  3 T 1 1 ~d_ T K i ~d_ þ 2~d_ T Gi E_ p þ ðE_ p ÞT P 1 E_ p ð4 þ Q P Q Þ  1 X 7 2 k¼ q jJj 6 4 5: V 0 i2I i i  1 d~_ T Gi Q  1 ðE_ p ÞT P 1 Q 2

step. The iteration starts from the hypothesis that the plastic strain rate is non-zero everywhere ((xICP)0 = 1.0), but from the second iterative step, the plastic strain rate will be updated and the non-plastic areas can be identified by qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi checking whether ~ d_ T K i ~ d_ þ 2~ d_ T Gi E_ p þ ðE_ p ÞT P 1 E_ p is equal to zero. All non-plastic areas will be found through a step-by-step technique. Then, from the second iterative step, the Gaussian integration point set I will be subdivided into two subsets: a subset (IR)k+1 where the rate of the plastic dissipation power is equal to zero and a subset (IP)k+1 where the plastic dissipation power occurs

I ¼ ðIR Þkþ1 [ ðIP Þkþ1 ; ð46aÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ðIR Þkþ1 ¼ i 2 I; ~d_ T K i ~d_ þ 2~d_ T Gi E_ p þ ðE_ p ÞT P 1 E_ p ¼ 0 ; ð46bÞ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðIP Þkþ1 ¼ i 2 I; ~d_ T K i ~d_ þ 2~d_ T Gi E_ p þ ðE_ p ÞT P 1 E_ p – 0 : ð46cÞ Once the non-plastic region (IR)k+1 is found, the objective function in Eq. (39) needs to be modified through introducing the non-differentiability by means of the penalty method. Finally, at each following iterative step, the equation set (42) will become 8  rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi     P > ICP _ > 1 > qi jJji x2 4 þ Q T P 1 Q K i ~d þ ðE_ p ÞT GTi  12 Q T GTi > V0 > > i2IP > > > > P > _ T T > > þb qi jJji ðK i ~d þ ðE_ p Þ Gi Þ ¼ 0; > > > < i2IR   qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  P ICP 1 qi jJji x2 ð4 þ Q T P 1 Q Þ P 1 E_ p þ GiT ~d_ T  12 P 1 Q > > V0 > > > i2IP >   > P > > > þb qi jJji P 1 E_ p þ GTi ~d_ T  LF R ¼ 0; > > > i2IR > > > : T _p R E ¼ 1;

where b is the penalty multiplier and its value varies from 106 to 1012 in practice. The iteration will be repeated until the following convergence conditions are satisfied:

2

ð44Þ

In order to trigger an iteration to linearize the equation set (42), an iteration seed is defined as:

ðxICP Þ0 ¼ 1:0;

ð45Þ

where the subscript ‘‘0’’ denotes that the variable is determined at step 0. Once, the iteration is triggered, the iteration factor xICP needs to be updated based on Eq. (43) at each following iterative step and the updated value of xICP will be used for calculating the equation set (42) at the next iterative step. However, it should be noted that when the value qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _~ dT K i ~ d_ þ 2~ d_ T Gi E_ p þ ðE_ p ÞT P 1 E_ p is equal to zero in some areas of the structure at an iterative step, the iterative factor xICP will become infinite and the objective function in Eq. (39) will be non-differentiable. The value qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ~ d_ þ 2~ d_ T Gi E_ p þ ðE_ p Þ T P 1 E_ p ¼ 0 corresponds to a d_ T K i ~ non-plastic deformation. In order to overcome the numerical difficult of non-differentiability, all non-plastic areas of the discretized RVE must be distinguished at each iterative

ð47a;b;cÞ

jkkþ1  kk j 6 g1 ; kkþ1

  ~_ _  dkþ1  ~dk  _ k~dkþ1 k

6 g2 ;

kE_ pkþ1  E_ kp k   6 g3 ; _p  Ekþ1 

ð48Þ

where k is the number of the iterative step, and g1, g2 and g3 are the computational error tolerances. Once the convergence conditions (48) are attained, the limit load multiplier at the current iterative step is the optimized solution. The above iterative process leads to the limit load multiplier k through a convergent sequence with monotonically decreasing series and a minimum optimal upper bound to the limit load multiplier will be obtained. 6. Applications The developed method is now applied to the limit analysis of fiber-reinforced composite materials at the microscopic scale through the study of RVEs. The effects of microstructures and anisotropy on macroscopic strength domain of heterogeneous and composite materials are numerically investigated.

581

H.X. Li / Mechanics of Materials 43 (2011) 574–585

Fiber-reinforced composite material is a typical heterogeneous material since reinforcing fibers embedded in a bonding matrix makes its mechanical properties rapidly varying from point to point. If the fibers are periodically embedded in the matrix, the average macroscopic properties of composites can be obtained through the study of RVE by means of the homogenization technique. Then, an equivalent macroscopic homogeneous material model can be proposed for further study of composite structures. In this section, by means of the developed method, microscopic limit analysis is implemented on unidirectional fiber-reinforced composites to get the macroscopic strength domain based on different microstructures and anisotropic constituents. The material of the matrix is assumed to be anisotropic and obeys the Tsai-Wu yield criterion. The plane stress model is adopted in the analysis and macroscopic loading is applied transverse to the fibers on the composites. The Tsai-Wu criterion in the plane stress condition can be expressed as FðrÞ ¼

A1

r0

r11 þ

A2

r0

r22 þ



B1

B2

B3

B4

r0

r0

r0

r0

r2 þ 2 11

r2 þ 2 22

r r þ 2 11 22

 1 ¼ 0;

r2 2 12



ð49Þ

where r0 is the reference yield stress. In order to investigate how anisotropy (the degree of anisotropy, the direction of anisotropy, and pressuredependence) affects the macroscopic strength domain of composites, three different anisotropic materials of the matrix are chosen in the numerical simulation, and their parameters used in the Tsai-Wu yield criterion (49) are

Table 1 The mechanical parameters of the matrix of a composite for the Tsai-Wu yield criterion. No. of the matrix

A1

A2

B1

B2

B3

B4

A B C

0 0 0.5

0 0 0.5

1 1 1

1 2 2

1 1.5 1.5

3 3 3

listed in Table 1. They represent three typical material behaviors: Matrix A is reduced to isotropic; Matrix B is orthotropic; Matrix C is generally anisotropic. The yield surfaces of these three materials in the plane (r11, r22) are plotted in Fig. 2. The material of the fiber is assumed to be isotropic and the von Mises criterion is used to describe its yield behavior. The yield strength of the fiber is chosen as rfs ¼ 20r0 in the numerical simulation. The fiber has a circular cross section and its volume fraction is Vf = 0.2. According to the arrangement of fibers in the matrix, both square and hexagonal RVEs chosen from the composite are investigated and they are discretized by eight-node quadrilateral finite elements, as shown in Fig. 3. The convergence tolerances chosen in the numerical iterative calculation are g1 = g2 = g3 = 103. By means of the developed numerical method, the macroscopic strength domain (R11, R22) of the composite are calculated and the numerical results are presented in Figs. 4–6. The numerical results in Figs. 4–6 clearly show that anisotropic constituents of a composite have a dramatic effect on the macroscopic average strength of the composite. Normally, due to the inhomogeneous feature of microstructures, even if all constituents of the composite are isotropic, the macroscopic strength of the composite still exhibits the orthotropic behavior. However, once there is a constituent of the composite obeys a generally anisotropic rule, e.g. the Tsai-Wu yield condition, the macroscopic strength domain of the composite will become significantly anisotropic and pressure-dependent. In order to investigate how reinforced fibers affect the macroscopic strength domain of composites, further numerical simulations are implemented for different fiber volume fractions. All material parameters for both the matrix and fiber are same with the above simulations and the numerical results are presented Figs. 7–10. From the numerical results, it can be concluded that in the plane stress condition, the macroscopic average strength of a composite is mainly determined by the weak phase. The strength and volume fraction of reinforced fibers only have

1.5 Matrix A

Sigma22/Sigma0

1 0.5 0 Matrix C

-0.5

Matrix B

-1 -1.5 -2 -2

-1.5

-1

-0.5 0 Sigma11/Sigma0

0.5

Fig. 2. The yield surfaces of the matrix materials.

1

1.5

582

H.X. Li / Mechanics of Materials 43 (2011) 574–585

Fig. 3. Finite element meshes of RVEs from a fibre-reinforced composite.

a slight reinforcing effect on the macroscopic strength domain of the composite. The weak phase determines the macroscopic average transverse strength of the composite and any microscopic material damage has a dramatic effect on the macroscopic strength of the composite.

7. Conclusions A general nonlinear programming approach is proposed to implement kinematic limit analysis for an anisotropic microstructure by means of introducing the homogenization theory into classical limit analysis. The macroscopic average strength domain of heterogeneous and composite

materials is calculated in a direct way. The homogenization theory provides an efficient bridge between the microscopic and macroscopic scales. Two well-separated scales can be well correlated. A nonlinear anisotropic material model is directly introduced into limit analysis without linearization. This reduces both computational error and constraint conditions. Meanwhile, a purely kinematic formulation is obtained to calculate the limit load multiplier, which does not involve any calculation of stress fields. The finite element model of limit analysis for a microstructure is finally formulated as a nonlinear programming problem subject to only one equality constraint. Therefore, the computational effort of the proposed approach is very modest. The developed method provides a powerful tool for the cal-

H.X. Li / Mechanics of Materials 43 (2011) 574–585

Fig. 4. The macroscopic strength domain of the composite (Vf=0.2, Matrix A).

Fig. 5. The macroscopic strength domain of the composite (Vf=0.2, Matrix B).

Fig. 6. The macroscopic strength domain of the composite (Vf=0.2, Matrix C).

culation of the macroscopic strength and the design of microscopic constituents of a composite. Numerical examples show that the anisotropic microstructure has a dramatic effect on the macroscopic average strength of a heterogeneous and composite material. A heterogeneous and composite material mostly exhibits anisotropic at the macroscopic scale, which results from the inhomogeneous feature of its microstructure. The inhomo-

583

Fig. 7. Effect of the volume fraction Vf of fibres on the macroscopic uniaxial strength R11 of the composite (Matrix A).

Fig. 8. Effect of the volume fraction Vf of fibres on the macroscopic uniaxial strength R11 of the composite (Matrix B).

Fig. 9. Effect of the volume fraction Vf of fibres on the macroscopic uniaxial tensile strength R11 of the composite (Matrix C).

geneous deformation of the microstructure makes the symmetry and pressure-dependence of an anisotropic composite much more complicated. The numerical examples also reveal that in the plane stress condition, the macroscopic strength domain of a composite is mainly determined by the weak phase. Therefore, a rational strength of the reinforced constituent is suggested based

584

H.X. Li / Mechanics of Materials 43 (2011) 574–585

Fig. 10. Effect of the volume fraction Vf of fibres on the macroscopic uniaxial compressive strength R11 of the composite (Matrix C).

on the microstructure of composites. This can provide a reference parameter for the microstructure design of composites. It should be emphasized that limit analysis is based on the condition of convexity and normality. This requires that the plastic flow of materials must be associated and the yield surface of materials must be convex. Meanwhile, convexity plays an extremely important role in nonlinear mathematical programming. For a minimum optimization problem subject to constraint conditions, if the feasible region is a convex set, any local minimum for the nonlinear programming is an optimal solution. Moreover, the necessary and sufficient condition which guarantees the local minimum obtained by the Kuhn–Tucker condition to be an optimal global solution, is that the objective function of the minimum optimization problem is convex. This requires that the mathematical expression of the yield function of materials must obey the condition of convexity, so that the fundamental theory of limit analysis is still provable from the extreme theorem and the optimization solution is the exact one. References Bazaraa, M.S., Sherali, H.D., Shetty, C.M., 2006. Nonlinear Programming: Theory and Algorithm. Wiley, New Jersey, USA. Buhan, P., Taliercio, A., 1991. A homogenization approach to the yield strength of composite materials. Eur. J. Mech. A/Solids 10, 129–154. Capsoni, A., Corradi, L., 1997. A finite element formulation of the rigidplastic limit analysis problem. Int. J. Numer. Methods Eng. 40, 2063– 2086. Capsoni, A., Corradi, L., Vena, P., 2001a. Limit analysis of orthotropic structures based on Hill’s yield condition. Int. J. Solids Struct. 38, 3945–3963. Capsoni, A., Corradi, L., Vena, P., 2001b. Limit analysis of anisotropic structures based on the kinematic theorem. Int. J. Plasticity 17, 1531– 1549. Chen, H.F., Liu, Y.H., Cen, Z.Z., Xu, B.Y., 1998. Numerical analysis of limit load and reference stress of defective pipelines under multi-loading systems. Int. J. Pres. Ves. Pip. 75, 105–114. Chen, M., Hachemi, A., Weichert, D., 2010. A non-conforming finite element for limit analysis of periodic composites. Proc. Appl. Math. Mech. 10, 405–406. Corradi, L., Vena, P., 2003. Limit analysis of orthotropic plates. Int. J. Plasticity 19, 1543–1566. Corradi, L., Luzzi, L., Vena, P., 2006. Finite element limit analysis of anisotropic structures. Comput. Methods Appl. Mech. Eng. 195, 5422–5436. Cristiansen, E., 1981. Computation of limit loads. Int. J. Numer. Method Eng. 17, 1547–1570.

Dallot, J., Sab, K., 2008a. Limit analysis of multi-layered plates, Part I: the homogenized Love–Kirchhoff model. J. Mech. Phys. Solids 56, 561– 580. Dallot, J., Sab, K., 2008b. Limit analysis of multi-layered plates. Part II: shear effects. J. Mech. Phys. Solids 56, 581–612. Drucker, D.C., Prager, W., Greenberg, H.J., 1952. Extended limit design theorems for continuous media. Quart. Appl. Math. 9, 381–389. Drucker, D.C., 1953. Limit analysis of two and three-dimensional soil mechanics problems. J. Mech. Phys. Solids 1, 217–226. Francescato, P., Pastor, J., 1997. Lower and upper numerical bounds to the off-axis strength of unidirectional fiber-reinforced composite by limit analysis methods. Eur. J. Mech. A/Solids 16, 213–234. Gao, Y., 1988. On the complementary bounding theorems for limit analysis. Int. J. Solids Struct. 24, 545–556. Hill, R., 1948. A variational principle of maximum plastic work in classical plasticity. Q. J. Mech. Appl. Math. 1, 18. Hill, R., 1951. On the state of stresses in a plastic-rigid body at the yield point. Phil. Mag. 42, 868–875. Himmelblau, D.M., 1972. Applied Nonlinear Programming. McGraw-Hill Book Company., New York. Koopman, D.C.A., Lance, R.H., 1965. On linear programming and plastic limit analysis. J. Mech. Phys. Solids 13, 77–87. Li, H.X., Liu, Y.H., Feng, X.Q., Cen, Z.Z., 2001. Micro/macromechanical plastic limit analyses of composite materials and structures. Acta Mech. Solida Sin. 14, 323–333. Li, H.X., Liu, Y.H., Feng, X.Q., Cen, Z.Z., 2003. Limit analysis of ductile composites based on homogenization theory. Proc. Roy. Soc. Lond. A 459, 659–675. Li, H.X., Yu, H.S., 2005. Kinematic limit analysis of frictional materials using nonlinear programming. Int. J. Solids Struct. 42, 4058–4076. Li, H.X., Yu, H.S., 2006a. Limit analysis of 2-D and 3-D structures based on an ellipsoid yield criterion. Acta Geotechnica 1, 179–193. Li, H.X., Yu, H.S., 2006b. Limit analysis of composite materials based on an ellipsoid yield criterion. Int. J. Plasticity 22, 1962–1987. Liu, Y.H., Cen, Z.Z., Xu, B.Y., 1995. A numerical method for plastic limit analysis of 3-D structures. Int. J. Solids Struct. 32, 1645–1658. Maier, G., Giacomini, S., Paterlini, F., 1979. Combined elastoplastic and limit analysis via restricted basis linear programming. Comput. Methods Appl. Mech. Eng. 19, 21–48. Makrodimopoulos, A., Martin, C.M., 2006. Lower bound limit analysis of cohesive-frictional materials using second-order cone programming. Inter. J. Numer. Methods Eng. 66, 604–634. Makrodimopoulos, A., Martin, C.M., 2007. Upper bound limit analysis using simplex strain elements and second-order cone programming. Int. J. Numer. Anal. Methods Geomech. 31, 835–865. Nocedal, J., Wright, S.J., 2006. Numerical Optimization. Springer, New York, USA. Pastor, J., Thai, T., Francescato, P., 2003. Interior point optimization and limit analysis: an application. Commun. Numer. Methods Eng. 19, 779–785. Pastor, F., Thore, P., et al., 2008. Convex optimization and limit analysis: application to Gurson and porous Drucker–Prager materials. Eng. Fract. Mech. 75, 1367–1383. Pastor, F., Loute, E., Pastor, J., 2009a. Limit analysis and convex programming: a decomposition approach of the kinematic mixed method. Int. J. Numer. Methods Eng. 78, 254–274. Pastor, F., Loute, E., Pastor, J., Trillat, M., 2009b. Mixed method and convex optimization for limit analysis of homogeneous Gurson materials: a kinematical approach. Eur. J. Mech. A/Solids 28, 25–35. Pellegrino, S., 1988. Efficient upper-bound limit analysis. Int. J. Mech. Sci. 30, 455–474. Ponter, A.R.S., Carter, K.F., 1997. Limit state solutions based upon linear elastic solutions with a spatially varying elastic modulus. Comput. Methods Appl. Mech. Eng. 140, 237–258. Ponter, A.R.S., Fuschi, P., Engelhardt, M., 2000. Limit analysis for a general class of yield conditions. Eur. J. Mech. A/Solids 19, 401–421. Suquet, P., 1987. Elements of homogenization for inelastic solid mechanics. In: Sanchez-Palencia, E., Zaoui, A. (Eds.), Homogenization Techniques for Composite Media, Lecture Notes in Physics, vol. 272. Springer, New York, pp. 193–278. Taliercio, A., 1992. Lower and upper bounds to the macroscopic strength domain of a fiber-reinforced composite material. Int. J. Plasticity 8, 741–762. Taliercio, A., Sagramoso, P., 1995. Uniaxial strength of polymeric-matrix fibrous composites predicted through a homogenization approach. Int. J. Solids Struct. 14, 2095–2123. Thai, T., Francescato, P., Pastor, J., 1998. Limit analysis of unidirectional porous media. Mech. Res. Commun. 25, 535–542.

H.X. Li / Mechanics of Materials 43 (2011) 574–585 Trillat, M., Pastor, J., 2005. Limit analysis and Gurson’s Model. Eur. J. Mech. A/Solids 24, 800–819. Turgeman, S., Pastor, J., 1982. Limit analysis: a linear formulation of the kinematic approach for axisymmetric mechanic problems. Int. J. Numer. Anal. Meth. Geomech. 6, 109–128. Weichert, D., Hachemi, A., Schwabe, F., 1999a. Application of shakedown analysis to the plastic design of composites. Archive Appl. Mech. 69, 623–633. Weichert, D., Hachemi, A., Schwabe, F., 1999b. Shakedown analysis composites. Mech. Res. Commun. 26, 309–318. Zhang, H., Liu, Y., Xu, B., 2004. Plastic limit analysis of periodic heterogeneous materials by a static approach. Key Eng. Mat., 739– 744.

585

Zhang, H., Liu, Y., Xu, B., 2009. Plastic limit analysis of ductile composite structures from micro- to macro-mechanical analysis. Acta Mech. Solida Sin. 22, 73–84. Zhang, P.X., Lu, M.W., Hwang, K., 1991. A mathematical programming algorithm for limit analysis. Acta Mech. Sin. 7, 267–274. Zhang, Y.G., Lu, M.W., 1994. Computational limit analysis of anisotropic axisymmetric shells. Int. J. Pres. ves. pip. 58, 283–287. Zhang, Y.G., Lu, M.W., 1995. An algorithm for plastic limit analysis. Comput. Methods Appl. Mech. Eng. 126, 333–341. Zouain, N., Herskovits, J., Borges, L.A., Feijóo, R.A., 1993. An iterative algorithm for limit analysis with nonlinear yield functions. Int. J. Solids Struct. 30, 1397–1417.