3 Elastoplastic Problems 3.1 Introduction There are numerous well-established techniques to calculate effective material characteristics for composite materials.In the case of composite components volume fractions only,one can use the closed form algebraic equations on upper and lower bounds or direct estimates for the effective material tensor components. Otherwise,the cell problems are formulated and solved using their Finite Element Method(FEM)or,alternatively,the Boundary Element Method(BEM)numerical implementations that enable direct computations of the effective characteristics. Recent advances in the area of computational methods in homogenisation of the nonlinear effective characterisation of heterogeneous materials and structures are reported in [4,85,86,107,112,136,250,325].In the same time,stochastic analysis is still being developed to estimate or to compute probabilistic moments of homogenised material tensors. Homogenisation of composite materials with elastoplastic constituents is presented below using the so-called Transformation Field Analysis (TFA) proposed by Dvorak and now applied to approximate the effective nonlinear behaviour of a three-component periodic composite.The self-consistent model and Mori-Tanaka theory,providing the estimation of the overall thermoelastic constants of composites on the basis of constituent properties and volume fractions, are partially incorporated in this model.Computational implementation of the method consists of the utilisation of the program ABAQUS to enable automatic homogenisation of n-component periodic composites in a general configuration of the components in the RVE.Numerical examples of the three-component periodic composite homogenisation make it possible to compare the nonlinear behaviour of a composite for its real and homogenised models in the case of the specific boundary problem defined for the cell.The next step in the development of this approach would be to determine the parameter sensitivity of the homogenised properties of the composite with respect to the material characteristics of the constituents as well as to some geometrical data defining the RVE.Statistical and stochastic simulation of probabilistic moments of the effective material tensors would be possible after such a sensitivity determination,taking into account the experimental knowledge of the statistical parameters of the composite constituents. 3.2 Homogenisation Method The periodic n-component composite in the plane orthogonal to the fibre direction is considered where perfectly bonded components are assumed to be
3 Elastoplastic Problems 3.1 Introduction There are numerous well-established techniques to calculate effective material characteristics for composite materials. In the case of composite components volume fractions only, one can use the closed form algebraic equations on upper and lower bounds or direct estimates for the effective material tensor components. Otherwise, the cell problems are formulated and solved using their Finite Element Method (FEM) or, alternatively, the Boundary Element Method (BEM) numerical implementations that enable direct computations of the effective characteristics. Recent advances in the area of computational methods in homogenisation of the nonlinear effective characterisation of heterogeneous materials and structures are reported in [4,85,86,107,112,136,250,325]. In the same time, stochastic analysis is still being developed to estimate or to compute probabilistic moments of homogenised material tensors. Homogenisation of composite materials with elastoplastic constituents is presented below using the so-called Transformation Field Analysis (TFA) proposed by Dvorak and now applied to approximate the effective nonlinear behaviour of a three-component periodic composite. The self-consistent model and Mori-Tanaka theory, providing the estimation of the overall thermoelastic constants of composites on the basis of constituent properties and volume fractions, are partially incorporated in this model. Computational implementation of the method consists of the utilisation of the program ABAQUS to enable automatic homogenisation of n-component periodic composites in a general configuration of the components in the RVE. Numerical examples of the three-component periodic composite homogenisation make it possible to compare the nonlinear behaviour of a composite for its real and homogenised models in the case of the specific boundary problem defined for the cell. The next step in the development of this approach would be to determine the parameter sensitivity of the homogenised properties of the composite with respect to the material characteristics of the constituents as well as to some geometrical data defining the RVE. Statistical and stochastic simulation of probabilistic moments of the effective material tensors would be possible after such a sensitivity determination, taking into account the experimental knowledge of the statistical parameters of the composite constituents. 3.2 Homogenisation Method The periodic n-component composite in the plane orthogonal to the fibre direction is considered where perfectly bonded components are assumed to be
164 Computational Mechanics of Composite Materials elastoplastic.Mechanical behaviour of the composite constituents is represented by time and temperature dependent constitutive relations under the assumption that for any time t the total strains and stresses can be decomposed as e(y,)=e(y,T)+e(y,) (3.1) o,(y,t)=o(y,t)+o,(y,t) (3.2) where e,denote elastic strain resulting from a given displacement boundary condition applied on the region and the eigenstrain in the same subregion, respectively:,stand for the elastic stress and eigenstress tensor components in The eigenstrain and eigenstress fields considered here as transformation field may be decomposed in the case of thermal and inelastic effects as y,t)=m,7(t)+n(y,T) (3.3) ,(y,t)=1a(t)+o(y,t) (3.4) where m denotes the thermal strain tensor,I,is the thermal stress tensor,while or(t)represents the linear thermal expansion coefficient.A procedure of determination of the effective thermal expansion coefficients for various composites has been described in [253.305.311].Since g is the inelastic strain and o"is the relaxation stress,(3.1)and (3.2)can be written as E,(y,t)=M,o,(y.t)+m,0(t)+gm(y.t) (3.5) o,(y.t)=C,E,(y,t)+1,0(t)+o(y.t) (3.6) where C.and M,are the elastic and compliance tensor components for the subregion ,Hence,it is possible to write the following relations between mMC and o M,=C- (3.7) m,=-M,1 (3.8) 1,=C,m, (3.9) Eimd(y.T)=-M,(y.t) (3.10) (y.)=-Cd(y.t) (3.11)
164 Computational Mechanics of Composite Materials elastoplastic. Mechanical behaviour of the composite constituents is represented by time and temperature dependent constitutive relations under the assumption that for any time τ the total strains and stresses can be decomposed as , ) , ) , ) r r r ε ( y τ = ε ( y τ + ε ( y τ el * (3.1) ( , ) ( , ) ( , ) * r r r σ y τ = σ y τ +σ y τ el (3.2) where el r ε , * r ε denote elastic strain resulting from a given displacement boundary condition applied on the region Ωr and the eigenstrain in the same subregion, respectively; el σ r , * σ r stand for the elastic stress and eigenstress tensor components in Ωr . The eigenstrain and eigenstress fields considered here as transformation field may be decomposed in the case of thermal and inelastic effects as ( , ) ( ) ( , ) * r ε y τ m α τ ε y τ inel = r T + r (3.3) ( , ) ( ) ( , ) * σ y τ l α τ σ y τ rel r = r T + r (3.4) where mr denotes the thermal strain tensor, rl is the thermal stress tensor, while αT(τ) represents the linear thermal expansion coefficient. A procedure of determination of the effective thermal expansion coefficients for various composites has been described in [253,305,311]. Since inel r ε is the inelastic strain and rel σ r is the relaxation stress, (3.1) and (3.2) can be written as ε (y,τ ) M σ (y,τ ) m θ ( ) ε (y,τ ) inel r r r r r = + t + (3.5) σ (y,τ ) C ε (y,τ ) l θ ( ) σ (y,τ ) rel r r r r r = + t + (3.6) where Cr and Mr are the elastic and compliance tensor components for the subregion Ωr . Hence, it is possible to write the following relations between inel r r r r r m ,l ,M ,C ,ε and rel σ r : −1 Mr = Cr (3.7) r r r m = −M l (3.8) r Crmr l = − (3.9) ε (y,τ ) M σ (y,τ ) rel r r inel r = − (3.10) σ (y,τ ) C ε (y,τ ) inel r r rel r = − (3.11)
Elastoplasticity problems 165 Mechanical and thermal elastic influence functions are given by the following relations: e,(,t)=A,(y)e(t)+D.(,ye(y,t) (3.12) o,(y,t)=B,(y)e(T)+F (y,yo;(y,t) (3.13) Matrices B,(y)and A,(y)in (3.12)and (3.13)denote stress and strain concentration factor tensors representing the volume averages of the corresponding functions over the periodicity cell,as is proposed in(3.14)to (3.17).To describe the overall homogenised response of volume the resulting strains and stresses are combined with their corresponding local components described by (3.3)to (3.6)as E(T)= e,n (3.14) =可。e.t0+ekn=e4a+ea oo)=回Jo,.ta2= (3.15) 司aic.t+dict-oo+oio Then,local elastic fields may be written as e4(t)= Aoem+aggaln (3.16) 1 o=可B,0ae+b,a,e1 (3.17) where a(y)and b(y)are the thermoelastic strain and stress concentration factors [86,94].The strain transformation field (y,t)defined in results in the displacements on the unconstrained part of surface o,while the transformation stress o(y,t)generates surface tractions on being constrained.The relation between the local and global transformation fields is proposed as e'(= le,)-Aeon= A,0e.a (3.18) oio=时a,-B,ao1n=向sain (3.19)
Elastoplasticity problems 165 Mechanical and thermal elastic influence functions are given by the following relations: ( , ) ( ) ( ) ( ) ( , ) * rs ε y τ A y ε τ D y, y ε y τ r r s = + ′ (3.12) ( , ) ( ) ( ) ( , ) ( , ) * rs σ y τ B y ε τ F y y σ y τ r r s = + ′ (3.13) Matrices ) B (y r and ) A (y r in (3.12) and (3.13) denote stress and strain concentration factor tensors representing the volume averages of the corresponding functions over the periodicity cell, as is proposed in (3.14) to (3.17). To describe the overall homogenised response of volume Ω , the resulting strains and stresses are combined with their corresponding local components described by (3.3) to (3.6) as ( , ) ( , )] ( ) ( ) 1 ( , ) 1 ( ) * * r el r ε τ ε τ ε τ ε τ ε τ ε τ + Ω = + Ω = Ω Ω = ∫ ∫ Ω Ω el r d d [ y y y (3.14) ( , ) ( , )] ( ) ( ) 1 ( , ) 1 ( ) * * r el r σ τ σ τ σ τ σ τ σ τ σ τ + Ω = + Ω Ω = Ω = ∫ ∫ Ω el V r d d [ y y y (3.15) Then, local elastic fields may be written as + Ω Ω = ∫ Ω r ( ) ( ) r ( ) T ( )] d 1 ( ) el ε τ [A y ε τ a y α τ (3.16) ∫ Ω + Ω Ω = r r T d el ( ) ( ) ( ) ( )] 1 σ (τ ) [B y σ τ b y α τ (3.17) where ar(y) and br(y) are the thermoelastic strain and stress concentration factors [86,94]. The strain transformation field ) ( , * ε y τ defined in Ω results in the displacements on the unconstrained part of surface ∂ Ω , while the transformation stress ) ( , * σ y τ generates surface tractions on Ω being constrained. The relation between the local and global transformation fields is proposed as ∫ ∫ Ω Ω Ω Ω − Ω = Ω = r r d r ( ) ( , ) d 1 ( , ) ( ) ( )] 1 ( ) * ε τ [ε y τ A y ε τ A y ε y τ * (3.18) ∫ ∫ Ω Ω Ω = Ω = Ω V T r d r ( ) ( , ) d 1 ( ) - ( ) ( )] 1 ( ) * r * σ τ [σ y B y σ τ B y σ y τ (3.19)
166 Computational Mechanics of Composite Materials The elastic local strain (y,t)and stress fields o(y,t)are found from a superposition of the elastic local fields (y,)and of(y,)with local eigenstrains E,(y,T)and eigenstresses o,(y,t),respectively;the same model in the context of global scale is introduced analogously.These two different scales are linked using the formulation for local strain and stress fields in the following form: e,(y')=A,(y')e+D.y,y')e(y') (3.20) o,(y)=B,(y)o+F.(y,y')o,(y') (3.21) D(y,y'),F(y,y')are transformation strain and stress influence functions, which enable us to relate the strain and stress tensor components on the macroscale defined by y and the microscale identified by y'.Solving the following boundary value problem on the RVE we get divo(y)= ∂o=0 ∂y (3.22) e,(y)=M,o,(y)+e(y) (3.23) 1 le,0an=ea (3.24) u(y)=ay+u"(y) (3.25) where the local uniform strain field E,is found using the matrices A,(y), D(y,y').Further,it is possible to determine the approximated expression of the averaged strain in the subvolume given as e,=A,e+∑D.e (3.26) Analogous to (3.26),the averaged stresses in the subregion 2,can be written in the form o,=B,0+∑Fn0i ,曰 (3.27) It is observed that F(y.y)and D.(y,y)are eigenstress and eigenstrain influence functions,that reflect the effect on the scale y resulting from a transformation on the scale y'under overall uniform applied stress or strain.The additional influence functions are determined in terms of averages and piecewise constant
166 Computational Mechanics of Composite Materials The elastic local strain ) ε (y,τ r and stress fields ) σ (y,τ r are found from a superposition of the elastic local fields ) ( , el r ε y τ and ) ( , el r σ y τ with local eigenstrains ) ( , * r ε y τ and eigenstresses ) ( , * r σ y τ , respectively; the same model in the context of global scale is introduced analogously. These two different scales are linked using the formulation for local strain and stress fields in the following form: ( ) ( ) ( , ) ( ) * y' A y' D y y' y' r r rs s ε = ε + ε (3.20) ( ) ( ) ( , ) ( ) * y B y F y y' y' σ r = r σ + rs σ s (3.21) D (y, y') rs , ) F (y, y' rs are transformation strain and stress influence functions, which enable us to relate the strain and stress tensor components on the macroscale defined by y and the microscale identified by y′. Solving the following boundary value problem on the RVE we get 0 ( ) ( ) = ∂ ∂ = y y y σ divσ (3.22) ( ) ( ) ( ) * r r r ε y = M σ y + ε y r (3.23) ∫ Ω Ω Ω = Ω ε d ε r ( ) 1 y (3.24) ( ) ( ) * u y = εy + u y (3.25) where the local uniform strain field r ε is found using the matrices ) A (y r , D (y, y') rs . Further, it is possible to determine the approximated expression of the averaged strain in the subvolume Ωr given as ∑ = = + N r s r r rs r , 1 * ε A ε D ε (3.26) Analogous to (3.26), the averaged stresses in the subregion Ωr can be written in the form ∑ = = + N r s r r rs r , 1 * σ B σ F σ (3.27) It is observed that Fr(y,y′) and Dr(y,y′) are eigenstress and eigenstrain influence functions, that reflect the effect on the scale y resulting from a transformation on the scale y’ under overall uniform applied stress or strain. The additional influence functions are determined in terms of averages and piecewise constant
Elastoplasticity problems 167 approximations in the introduced subregions inside the RVE.Therefore,under overall strain e(t)=0,the transformation concentration factor tensor Ds gives the strain induced in the subvolume by a unit uniform eigenstrain in ,Under overall stress o(t)=0,the concentration factor tensor Fs defines the stress in resulting from the unit eigenstrain located in ,It can be shown that these tensors can be determined in the case of multiphase medium as D.r=(I-A,)(C,-C)(6I-e,AT)C, (3.28) Fr =(I-B,)(M,-M)6I-c,B M, (3.29) (r,s=1,...,N,without summation over repeated indices) which for a two-component composite gives Dpa=(I-Ap)(Ca-CB)Ca (3.30) DpB=-(I-Ap)(Ca-Cp}'CB Fpa=(I-Bp)(Ma-MB)Ma (3.31) FpB=-(I-Bp)(Ma-Mg)'Mg forp=oβ This completes the description of the homogenisation method for a composite with elastoplastic coefficients by use of the Transformation Field Analysis (TFA).It should be underlined that,in comparison to the homogenisation model valid for the linear elastic range,the necessity of transformation matrix computations is crucial for the proposed nonlinear FEM analysis. 3.3 Finite Element Equations of Elastoplasticity The following boundary value problem is now considered [206,210]: △oHu=0;x∈2 (3.32) △6u=Cumn AEm;x∈2 (3.33) △Emm=2[△k.1+△.k+,k△1+△.k1+△.k△小；X∈2 (3.34) with the boundary conditions △ca乃=△g;xe02。,k=l,2,3 (3.35) △lk=Ait;xe∂2n,R=l,2,3 (3.36)
Elastoplasticity problems 167 approximations in the introduced subregions inside the RVE. Therefore, under overall strain ε(t)=0, the transformation concentration factor tensor Drs gives the strain induced in the subvolume Ωr by a unit uniform eigenstrain in Ωs . Under overall stress σ(t)=0, the concentration factor tensor Frs defines the stress in Ωr resulting from the unit eigenstrain located in Ωs . It can be shown that these tensors can be determined in the case of multiphase medium as ( )( ) ( ) s T Dsr = I − Ar Cr − C rsI − csAs C − δ 1 (3.28) ( )( ) ( ) s T Fsr = I − Br Mr − M rsI − csBs M − δ 1 (r,s=1,...,N, without summation over repeated indices) (3.29) which for a two-component composite gives ( )( ) α A Cα Cβ Cα D I −1 p = − p − ( )( ) β A Cα Cβ Cβ D I −1 p = − − p − (3.30) ( )( ) α B Mα Mβ Mα F I −1 p = − p − ( )( ) β B Mα Mβ Mβ F I −1 p = − − p − for p=α,β (3.31) This completes the description of the homogenisation method for a composite with elastoplastic coefficients by use of the Transformation Field Analysis (TFA). It should be underlined that, in comparison to the homogenisation model valid for the linear elastic range, the necessity of transformation matrix computations is crucial for the proposed nonlinear FEM analysis. 3.3 Finite Element Equations of Elastoplasticity The following boundary value problem is now considered [206,210]: ∆σ kl ,l = 0 ; x ∈Ω (3.32) kl Cklmn mn ∆σ = ∆ε ~ ; x ∈Ω (3.33) [ ] 2 , , , , , , , , 1 ∆ mn = ∆uk l + ∆ul k + ui k∆ui l + ∆ui kui l + ∆ui k∆ui l ε ; x ∈Ω (3.34) with the boundary conditions kl l k ∆σ n = ∆t ; σ x ∈∂ Ω , k = 1,2,3 (3.35) k k u ˆ u ˆ ∆ = ∆ ˆ ; ∈ Ωu x ∂ , 3 1,2, ˆ k = (3.36)
168 Computational Mechanics of Composite Materials This problem is solved for displacements u(x),strain Ex)and stress ox) tensor components fulfilling equilibrium equations(3.32)-(3.36).Let us note that the stress tensor increments Ao(x),Ad(x)denote here the first and second Piola-Kirchhoff tensors △oH=△Fm△dnm+fFm△Gnm+△FOm;x∈2 (3.37) where △Fm=△lkm;X∈2 (3.38) To get the solution,the following functional defined on the displacement increments as Au is introduced: JAue)=∫CumAE△emn+ouA:kAu:u)dQ-∫△ieAu:d旧D (3.39) Let us note that this methodology is common for homogeneous m aterials as well as heterogeneous media.In case of composites,the last equation can be decomposed into the integrals valid for particular constituents and their boundaries and interfaces,separately. Now,let us introduce the displacement increment function Au(x)being continuous and differentiable on and,consequently,including all geometrically continuous and coherent subsets (finite elements).,e=1,...,E discretising the entire It is not assumed that Au(x)is differentiable on the interelement surfaces and boundaries (for ef=1.....,ef).Next,let us consider the following approximation of Au(x)for xe: Ne Au:x)=paK)△ (3.40) =1 wherex)are the shape functions for node k.Au)represents the degrees of freedom (DOF)vector,while Ne is the total number of the DOF in this node. Considering above,the displacements and strains gradients are rewritten as follows: Ak()=p(x)△n (3.41) △EuK)=[+B]△”=略△ (3.42) △Eu(K)=B△u△n (3.43)
168 Computational Mechanics of Composite Materials This problem is solved for displacements ( ) x k u , strain ( ) x kl ε and stress ( ) x σ kl tensor components fulfilling equilibrium equations (3.32)-(3.36). Let us note that the stress tensor increments ( ) x ∆σ kl , ( ) x σ kl ~ ∆ denote here the first and second Piola-Kirchhoff tensors σ kl Fkm σ ml Fkm σ ml Fkmσ ml ~ ~ ~ ∆ = ∆ ∆ + ∆ + ∆ ; x ∈Ω (3.37) where ∆Fkm = ∆uk ,m ; x ∈Ω (3.38) To get the solution, the following functional defined on the displacement increments as k ∆u is introduced: ( ) ( ) ( ) ∫ ∫ Ω Ω ∆ = ∆ ∆ + ∆ ∆ Ω − ∆ ∆ Ω ∂ J uk Cklmn ε kl ε mn σ kl ui k ui l d t ˆ k ukd ∂ ~ 2 , , 1 2 1 (3.39) Let us note that this methodology is common for homogeneous m aterials as well as heterogeneous media. In case of composites, the last equation can be decomposed into the integrals valid for particular constituents and their boundaries and interfaces, separately. Now, let us introduce the displacement increment function ( ) x k ∆u being continuous and differentiable on Ω and, consequently, including all geometrically continuous and coherent subsets (finite elements) Ωe , e=1,...,E discretising the entire Ω. It is not assumed that ( ) x k ∆u is differentiable on the interelement surfaces and boundaries ∂Ωef (for e,f=1,...,E, e ≠ f ). Next, let us consider the following approximation of ( ) x k ∆u for x∈Ω : () () ( ) 1 N N uk k u e ζ ζ ∆ = ∑ϕζ ∆ = x x (3.40) where ϕζk ( ) x are the shape functions for node k, (N ) ∆uζ represents the degrees of freedom (DOF) vector, while Ne is the total number of the DOF in this node. Considering above, the displacements and strains gradients are rewritten as follows: () () ( ) , , N k l k l u uζ ζ ∆ x = ϕ x ∆ (3.41) ( ) (1) (2) ( ) ( ) [ ] N kl N kl kl kl B B u B uζ ζ ζ ζ ζ ∆ε x = + ∆ = ∆ (3.42) ( ) (N ) (N) kl Bkl uζ uξ ζξ ∆ε x = ∆ ∆ (3.43)
Elastoplasticity problems 169 and finally △eH(x)=△Eu&)+△Eu(x) (3.44) The following notation is applied (3.42)and (3.43): Bx)=p{(&) (3.45) B2s(区)=p&p(xg (3.46) B(x)=xi(x) (3.47) All these equations are substituted into the variational formulation of the problem(cf.(3.39)).There holds 3Chn△Eu△emn=Chmm△Eu+△Eu)△Em+△Enm） =Cklmn (AEgAEmn+△Eu△E,n+△Eu△Emm+AEgAEmn） (3.48) =号Chn⑧i△nB点An”+略AnB△△ +B5△△BA”+BfA△n△A） 6u△lu△4u=6up(K△upi(xAu (3.49) Next,the following notation is applied: kge=∫6up(ug9i(aa2 (3.50) 2 kg=∫CumB Bd2 (3.51) 42 k号=∫分Cmn®2+2S+2Ba)dn (3.52) 2 where 缨=kg+k+k号 (3.53) and for the second and third order terms
Elastoplasticity problems 169 and finally ( ) x ( ) x ( ) x kl kl kl ∆ε = ∆ε + ∆ε (3.44) The following notation is applied (3.42) and (3.43): () () x x ζ ζ Bkl ϕk ,l (1) = (3.45) () () () ( ) , , (2) N kl i k i l B uξ ζ ζ ξ x = ϕ x ϕ x (3.46) ( ) x () () x x ζξ ζ ξ Bkl 2 ϕi,k ϕi,l 1 = (3.47) All these equations are substituted into the variational formulation of the problem (cf. (3.39)). There holds ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 1 2 1 2 1 2 1 N N mn N N kl N mn N N kl N N mn N kl N mn N klmn kl klmn kl mn kl mn kl mn kl mn klmn kl mn klmn kl kl mn mn B u u B u B u u B u u C B u B u B u B u u C C C µ ν µν ζ ζ ζξ µ µ ζ ξ ζξ µ ν µν ζ ζ ξ ξ ζ ζ ε ε ε ε ε ε ε ε ε ε ε ε ε ε + ∆ ∆ ∆ + ∆ ∆ ∆ ∆ = ∆ ∆ + ∆ ∆ ∆ = ∆ ∆ + ∆ ∆ + ∆ ∆ + ∆ ∆ ∆ ∆ = ∆ + ∆ ∆ + ∆ (3.48) () () ( ) , ( ) 2 , 1 2 , , 1 ~ ~ N i l N kl i k i l kl i k u u u uξ ξ ζ ζ σ ∆ ∆ = σ ϕ x ∆ ϕ x ∆ (3.49) Next, the following notation is applied: () () ∫ Ω = Ω e k u d i l N kl i k e x x ξ ζ σ ζ ζξ σ ϕ ϕ , ( ) , ( ) ~ (3.50) ∫ Ω = Ω e k CklmnBkl Bmn d con e ζ ξ ζξ (1) (1) 2 ( ) 1 (3.51) ∫ ( ) Ω = + + Ω e k C B B B B B B d klmn kl mn kl mn kl mn u e ζ ξ ζ ξ ζ ξ ζξ (1) (2) (2) (1) (2) (2) 2 ( ) 1 (3.52) where e e con e u e k k k k (1) ( ) ( ) ( ) ζξ ζξ σ ζξ = ζξ + + (3.53) and for the second and third order terms
170 Computational Mechanics of Composite Materials (3.54) =J房Cm@号△△ArW+Bf△A△d k2=∫2Chn同△n△△r△）a (3.55) Ω Introducingk for i=1,2.3 to the functional (Au)in(3.39)and applying the transformation from the local to the global system by the use of the following formula,typical for the FEM implementation: △ng=asa△qa (3.56) it is obtained that J△qa)=Kw△qa△qB+K编△qa△qB△q7 (3.57) +tKs△qa△qB△q,△qs-△0.△ga The stationarity of the functional /(Ada)leads to the following algebraic equation: Kw△qB+KAqB△q,+Ks△qB△q,Aq6=△0 (3.58) being fulfilled for any configuration of The iterative numerical solution of this equation makes it possible,according to the specified boundary conditions,to compute the effective constitutive tensor components of the homogenised composite.It should be stressed that the first two components of the stiffness matrix are considered only in further numerical analysis(geometrical nonlinearity is omitted in the homogenisation process);a detailed description of the numerical integration and solution of (3.58)can be found in [12,72,271,276],for instance. 3.4 Numerical Analysis As was mentioned above,the main goal of the TFA approach is to compute the transformation matrices A,,Ds that are determined only once for the original geometry of the composite and assuming initially linear elastic characteristics of the constituents.There holds that (da,=C(de,)C(des+de)
170 Computational Mechanics of Composite Materials ( ) ∫ Ω = ∆ ∆ ∆ + ∆ ∆ ∆ Ω e C B u B u u B u u B u d k N mn N N kl N N mn N klmn kl e ( ) ( ) ( ) ( ) ( ) ( ) 2 3 (2) µ µ ζ ξ ζξ µ ν µν ζ ζ ζξ (3.54) ∫ ( ) Ω = ∆ ∆ ∆ ∆ Ω e k C B u u B u u d N N mn N N klmn kl (3)e ( ) ( ) ( ) ( ) 2 µ ν µν ζ ξ ζξ ζξ (3.55) Introducing (i) kζξ for i=1,2,3 to the functional ( ) k J ∆u in (3.39) and applying the transformation from the local to the global system by the use of the following formula, typical for the FEM implementation: uζ aξα qα N ∆ = ∆ ( ) (3.56) it is obtained that ( ) αβγδ α β γ δ α α α αβ α β αβγ α β γ K q q q q Q q J q K q q K q q q + ∆ ∆ ∆ ∆ − ∆ ∆ ∆ = ∆ ∆ + ∆ ∆ ∆ (3) 4 1 (2) 3 (1) 1 2 1 (3.57) The stationarity of the functional ( ) α J ∆q leads to the following algebraic equation: Kαβ∆qβ + Kαβγ∆qβ∆qγ + Kαβγδ∆qβ∆qγ∆qδ = ∆Qα (1) (2) (3) (3.58) being fulfilled for any configuration of Ω. The iterative numerical solution of this equation makes it possible, according to the specified boundary conditions, to compute the effective constitutive tensor components of the homogenised composite. It should be stressed that the first two components of the stiffness matrix are considered only in further numerical analysis (geometrical nonlinearity is omitted in the homogenisation process); a detailed description of the numerical integration and solution of (3.58) can be found in [12,72,271,276], for instance. 3.4 Numerical Analysis As was mentioned above, the main goal of the TFA approach is to compute the transformation matrices Ar, Drs that are determined only once for the original geometry of the composite and assuming initially linear elastic characteristics of the constituents. There holds that r r r d ( d d ) el r r Ω Ω Ω σ ε ε εinel r eff r r eff r = C d = C +
Elastoplasticity problems 171 a=Ae+2ga (3.59) ao）n,=ce）n,and）n=c.(ae）a (3.60) Cem =C +C (3.61) Further,using spatial averaging definitions,the averaged stress tensor components are calculated as follows: (e,)o=8 and (o)o=o (3.62) Hence,the effective elasticity tensor components C are derived for a given increment as do=C de (3.63) c-立cc+cb:。 (3.64) In the particular case of a two-component composite,the transformation and concentration matrices are obtained as,cf.(3.30)and(3.31) D=I-A)(C-C2)C (3.65) D21=(I-A2)C1-C2)PC (3.66) D2=-(I-A)C1-C2)-C3 (3.67) D2=-I-A2)C-C2)'C2 (3.68) C.C2 denote here the components corresponding to elastic properties,while A,,A,are mechanical concentration matrices.Finally,using (3.64)it is obtained that ci-立GC+c,Cg+sDe':o+e,De'o (3.69) +cD:d+Dd: The FEM aspects of TFA computational implementation are discussed in detail in Section 3.4 below.Further,it should be noticed that there were some approaches in the elastoplastic approach to composites where,analogously to the linear
Elastoplasticity problems 171 ∑ = = + N A D r,s 1 r s r d d d Ω Ω ε ε εinel r rs s (3.59) r r el r el r d d Ω Ω σ C ε el = r and r r inel r inel r d d Ω Ω σ C ε = r (3.60) r el C = Cr + C eff r (3.61) Further, using spatial averaging definitions, the averaged stress tensor components are calculated as follows: ε ε Ω r = and σ σ Ω r = (3.62) Hence, the effective elasticity tensor components eff C are derived for a given increment as σ = dε eff d C (3.63) ( ) inel r N r,s inel rs N r el r r eff C c C D ε : σ 1 1 r s 1 c − = = = ∑ + ∑ (3.64) In the particular case of a two-component composite, the transformation and concentration matrices are obtained as, cf. (3.30) and (3.31) 1 1 11 1 1 2 D (I A )(C C ) C− = − − (3.65) 1 1 21 2 1 2 D (I A )(C C ) C− = − − (3.66) 2 1 12 1 1 2 D (I A )(C C ) C− = − − − (3.67) 2 1 22 2 1 2 D (I A )(C C ) C− = − − − (3.68) 1 2 C ,C denote here the components corresponding to elastic properties, while 1 2 A ,A are mechanical concentration matrices. Finally, using (3.64) it is obtained that ∑ ∑ () () = − − = + + + N r eff el el inel inel inel inel c c 2 2 1 1 2 12 2 1 1 1 2 2 1 11 1 C C C c D ε :σ c D ε :σ ∑ () () − − + + inel inel inel inel 2 1 1 2 22 2 1 1 21 1 c D ε :σ c D ε :σ (3.69) The FEM aspects of TFA computational implementation are discussed in detail in Section 3.4 below. Further, it should be noticed that there were some approaches in the elastoplastic approach to composites where, analogously to the linear
172 Computational Mechanics of Composite Materials elasticity homogenisation method,the approximation of the effective yield limit stresses of a composite is proposed as a quite simple closed form function (3.70) or,in terms of the effective yield surface,in the following form: (o) aow-bzy (3.71) where m(y=)=3 uV(u)and Vis any estimate of the viscosity compliance tensor defined using the viscosities u and u2.A review of the most recent theories in this field can be found in [381],for instance. The main aim of computational experiment presented is to determine the global nonlinear homogenised constitutive law for two component composites with elastoplastic components;the FEM based program ABAQUS [1]is used in all computational procedures.However the method presented can be implemented in any nonlinear FEM plane strain/stress code such as [60],for instance.The numerical experiments are carried out in the microstructural(RVE)level,and that is why the global response of the composite is predicted starting from the behaviour of the periodicity cell.The numerical micromechanical model consists of a three-component periodicity cell with horizontal and vertical symmetry axes and dimensions 3.0 cm (horizontal)X 2.13 cm (vertical)(cf.Figure 3.1 and 3.2). The composite is made of epoxy matrix and metal reinforcement with material properties of the components collected in Table 1.The void embedded into the steel casting simulates a lack of any matrix in the periodicity cell.Some nonzero material data are introduced to avoid numerical singularities during the homogenisation problem solution. The 10-node biquadratic,quadrilateral hybrid linear pressure reduced integration plane strain finite elements with 4 integration Gaussian points are used to discretise the cell.Periodic boundary conditions are imposed to ensure periodic character of the entire structure behaviour.A suitable formulation of displacement boundary conditions has the following form: 4,=e(y(P2)-y(B) (3.72) where=represents the displacement function components,E is the global total strain imposed on the periodicity cell,while y(P)and y(P)denote coordinates of the points lying on the opposite sides of the RVE
172 Computational Mechanics of Composite Materials elasticity homogenisation method, the approximation of the effective yield limit stresses of a composite is proposed as a quite simple closed form function 1 2 σ = σ σ (eff) (3.70) or, in terms of the effective yield surface, in the following form: ( ) [ ] ( ) ⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ − ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Φ = + − ≥ 2 (2) 1 (2) (1) 2 1 0 max ( ) y y y y c c y σ m y σ σ σ σ σ (3.71) where ( ) 3 2 1 = = µ µ m y ( ) 2 1 2 µ V µ ,µ and V is any estimate of the viscosity compliance tensor defined using the viscosities µ1 and µ2. A review of the most recent theories in this field can be found in [381], for instance. The main aim of computational experiment presented is to determine the global nonlinear homogenised constitutive law for two component composites with elastoplastic components; the FEM based program ABAQUS [1] is used in all computational procedures. However the method presented can be implemented in any nonlinear FEM plane strain/stress code such as [60], for instance. The numerical experiments are carried out in the microstructural (RVE) level, and that is why the global response of the composite is predicted starting from the behaviour of the periodicity cell. The numerical micromechanical model consists of a three-component periodicity cell with horizontal and vertical symmetry axes and dimensions 3.0 cm (horizontal) × 2.13 cm (vertical) (cf. Figure 3.1 and 3.2). The composite is made of epoxy matrix and metal reinforcement with material properties of the components collected in Table 1. The void embedded into the steel casting simulates a lack of any matrix in the periodicity cell. Some nonzero material data are introduced to avoid numerical singularities during the homogenisation problem solution. The 10-node biquadratic, quadrilateral hybrid linear pressure reduced integration plane strain finite elements with 4 integration Gaussian points are used to discretise the cell. Periodic boundary conditions are imposed to ensure periodic character of the entire structure behaviour. A suitable formulation of displacement boundary conditions has the following form: ( ( ) ( )) 2 P1 u y P y i = ε ij − (3.72) where { } 1 2 ui = u ,u represents the displacement function components, ij ε is the global total strain imposed on the periodicity cell, while ) (P1 y and ) (P2 y denote coordinates of the points lying on the opposite sides of the RVE