Extended multiscale finite element method for small-deflection analysis of thin composite plates with aperiodic microstructure characteristics

Composite Structures 160 (2017) 422–434

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Extended multiscale finite element method for small-deflection analysis of thin composite plates with aperiodic microstructure characteristics Mingfa Ren, Jie Cong, Bo Wang ⇑, Xinglin Guo Department of Engineering Mechanics, State Key Laboratory of Structural Analyses for Industrial Equipment, Dalian University of Technology, Dalian 116024, China

a r t i c l e

i n f o

Article history: Received 21 July 2016 Revised 27 September 2016 Accepted 18 October 2016 Available online 19 October 2016 Keywords: Extended multiscale finite element method Thin composite plate Coupling effect Aperiodic microstructure characteristics

a b s t r a c t Based on the theoretical framework of the extended multiscale finite element method, an efficient multiscale finite element method is developed for small-deflection analysis of thin composite plates with aperiodic microstructure characteristics. First of all, the decoupling displacement boundary conditions for deflections and rotations are reconstructed based on their coupled displacement modes of thin composite plates. Then, the multiscale base functions can be constructed numerically under the boundary conditions through the micro scale computations. Moreover, the coupling effects of composite laminates, especially asymmetric laminates are considered by introducing the additional coupling terms among translations and rotations into the multiscale base functions. Finally, the macroscopic and microscopic response fields are obtained through the macro scale and downscale computation respectively on the basis of multiscale base functions. Numerical examples demonstrate that the developed method possesses high computing accuracy and efficiency compared with the conventional finite element method. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Fiber-reinforced composite material has been widely used in the field of aerospace due to its high strength-weight, high stiffness-weight and outstanding designability characteristic. It is worth noting that the analysis of composite structures is usually a multiscale problem as illustrated in Fig. 1. The relevant levels of composite structures from low to high are constituent level, lamina or laminate level and structural level [1]. Consider that the properties and responses of composite structures at higher levels are influenced by those at lower levels. In order to utilize and design composite structures better, it is important to implement the multiscale analysis. As shown in Fig. 1, the microscopic characters of composite structures at different levels could be various and aperiodic, such as the distributions of constituents, microscopic holes and hierarchical stiffeners [2–6]. Since the conventional finite element method (FEM) needs mesh refinement, which leads to a huge number of elements and low computing efficiency. Therefore, to find a multiscale numerical computation method with high computing efficiency and required precision has become a research hotspot in recent years. The multiscale method in the field of fiber-reinforced composite material is mainly the micromechanics method [7,8]. It aims at ⇑ Corresponding author. E-mail address: [email protected] (B. Wang). http://dx.doi.org/10.1016/j.compstruct.2016.10.073 0263-8223/Ó 2016 Elsevier Ltd. All rights reserved.

establishing the relationship between the macroscopic mechanical properties and microstructure characteristics. At present, the micromechanics methods which have been widely applied mainly include the asymptotic homogenization method, the representative volume element method, the multiscale finite element method, etc. The asymptotic homogenization method has strict mathematical basis, which was first proposed by Bensoussan [9] for the heterogeneous material with periodic structural characteristics. Based on the asymptotic homogenization theory, Guedes [10] proposed the computational homogenization method. In this method, the FEM is applied to obtain the homogenized macroscopic material parameter and the microscopic valuable information quantities [11]. The computational homogenization method has been applied to solve the nonlinear problems of the composites [12–15]. For the heterogeneous material with random microscopic characteristics, the homogenization theory generally needs to select a relatively large representative volume element to contain its microscopic characteristics due to the restrictions of periodicity assumption. That leads to a loss of computing efficiency. The representative volume element method is a kind of the homogenization method based on numerical values. This method can obtain the equivalent mechanical properties of heterogeneous material as well as the macroscopic and microscopic stress and strain fields by choosing a representative volume element and making analysis after setting specific boundary conditions [16].

M. Ren et al. / Composite Structures 160 (2017) 422–434

423

Fig. 1. Illustration of the relevant levels for fiber-reinforced composite structures.

The boundary conditions frequently used in the representative volume element method include Dirichlet boundary condition, Neumann boundary condition and periodic boundary condition. Hazanov [17] compared the effect of different boundary conditions on the accuracy of the representative volume element method and showed that the macroscopic elastic modulus of heterogeneous material tend to upper bound and lower bound respectively under Dirichlet and Neumann boundary conditions. Because the selection of boundary conditions of the representative volume element method is demanding, there is no guarantee that the boundary conditions imposed forcibly consist with actual deformations. Therefore, the predicted results may be distorted. In addition, the selection of representative volume elements depends entirely on the user’s experience. Thus, the size and size effect of representative volume elements should be predetermined [18,19]. The basic idea of the multiscale finite element method (MsFEM) can date back to the research work of Babuska [20,21]. Afterwards, Hou et al. [22,23] further developed and officially named it as the multiscale finite element method. This method achieves multiscale analysis by constructing the multiscale base functions which reflect the microstructure characteristics. The construction of the multiscale base functions is different from the conventional finite element base functions. They are obtained by numerically solving the displacement field of macroscopic elements, so the material properties within the macroscopic elements do not have to be uniform. Since the problems can be solved at the macro scale directly on the basis of the multiscale base functions, the computing efficiency is improved. Compared with other multiscale methods, the MsFEM is not restricted by the periodicity assumption and can execute the downscale computation easily. The MsFEM has been used for the flow numerical simulation [24,25]. Moreover, several similar multiscale computational methods have been developed, including the variational multiscale method [26,27], the Voronoi cell method [28,29] the multiscale finite volume method [30,31], generalized multiscale finite element method [32–34], etc. However, for the vector field (the stress and strain fields) problems of composite materials, the methods above are still not applicable. The extended multiscale finite element method (EMsFEM) is developed from the MsFEM by Zhang et al. [35]. In order to consider the coupling effects among different directions for the multi-dimensional vector field problems, the additional coupling terms of the multiscale base functions are introduced. The research [35] shows that the EMsFEM has the advantages of high computing efficiency and accuracy. It has been successfully applied in the research of the vector field problems of heterogeneous porous media and periodic truss material [36–39]. At present, less literature has discussed about the applications of the EMsFEM for the vector filed problems of thin composite plates. The difficulty lies in the following aspects. First, the anisotropy and layup configurations of composite plates increase the complexity of multiscale analysis. Second, the additional terms

within the multiscale base functions are only employed to consider the coupling effects among different directions in the present EMsFEM. But for the fiber-reinforced composite laminates, especially asymmetric laminates, the other additional coupling terms need to be introduced to consider the coupling effects among translations and rotations. Last but not least, it is necessary to reconstruct the decoupling displacement boundary conditions for deflections and rotations when constructing the multiscale base functions of thin composite plates because the displacement modes of deflections and rotations are both nonlinear and coupled with each other based on the thin plate theory. Based on the theoretical framework of the EMsFEM, a novel multiscale finite element method suitable for small-deflection analysis of thin composite plates is developed in this paper. First, the macro scale and micro scale finite element formulations of composite plates with the anisotropy and layup configurations are derived. Then, the decoupling displacement boundary conditions for deflections and rotations are reconstructed by assigning a unit value and zero on the macroscopic nodes and interpolating on the microscopic nodes according to their coupled displacement modes of thin composite plates. Further, the multiscale base functions are constructed numerically under the boundary conditions through the micro scale computations. Consider that the coupling effects of composite laminates, especially asymmetric laminates are strong. The other additional coupling terms among translations and rotations are introduced into the multiscale base functions. Finally, the macroscopic and microscopic response fields are obtained through the macro scale and downscale computation respectively on the basis of multiscale base functions. The numerical examples show that the developed method possesses high computing accuracy for the thin composite plates in consideration of layup configurations and aperiodic microstructure characteristics. In addition, compared with the conventional finite element method, the computing efficiency is improved significantly. 2. Methodology 2.1. Basic principle The developed method consists of three steps: the micro scale, macro scale and downscale computations. The micro scale computation is employed to construct the multiscale base functions, and then the macroscopic displacement fields can be obtained at the macro scale. The downscale computation is employed to obtain the microscopic response fields, including the strain and stress fields. The characteristics of the EMsFEM for small-deflection analysis of thin composite plates are as follows. First, the anisotropy and layup configurations of composite plates are considered at the micro scale computation and associated with the macroscopic computation through the multiscale base functions. Second, the

424

M. Ren et al. / Composite Structures 160 (2017) 422–434

additional coupling terms of the multiscale base functions are introduced to consider the coupling effects among translations and rotations of composite laminates. Third, the decoupling displacement boundary conditions for deflections and rotations are reconstructed according to their coupled displacement modes of thin composite plates. 2.2. Macro scale and micro scale finite element formulations The small-deflection theory for thin composite plates is assumed: A normal perpendicular to the middle surface before deformation is still perpendicular to the middle surface after deformation and the length of the normal remains constant. Since it is difficult to guarantee that the stress components satisfy the stress boundary conditions accurately. Then, only the internal forces obtained by integrating stress components along thickness direction are required to satisfy the boundary conditions. The equilibrium equations and boundary conditions for the smalldeflection problems of thin composite plates are given as:

8
> < KF þ f ¼ 0 in X nF ¼ T
on Cr > :
a¼a on Cu

ð1Þ

where K is the differential operator, F is the internal force vector of unit width,
f is the force vector of unit area, n is the unit normal
is the displacement matrix, T
is the force boundary condition, and a boundary condition. At the micro scale (i.e., one macroscopic element E), a composite plate is discretized into the microscopic finite elements (e) as shown in Fig. 2. The displacement field a within a microscopic element can be expressed as:

a ¼ Ne ae

By constructing a transformation matrix GE between the microscopic nodes within a microscopic element and the microscopic nodes within all the microscopic elements, the displacement vector ae can be expressed as

ae ¼ GE aE

where aE is the displacement vector of microscopic nodes within all the microscopic elements. With Eqs. (3) and (4), the micro scale finite element formulation is derived as Eq. (5) based on the minimum potential energy principle.

K E aE ¼ PE

    XZ X1 Z A B 0 B ae dxdy  aTe B0T aTe NTe
f dxdy 2 Xe B D Xe e e ( ) X Z
 aTe N Te TdS

PE ¼

e

ð5Þ

where

   X T R A B 0 B dxdy GE GE Xe B0T B D e o X T nR R T

PE ¼ GE Xe Ne f dxdyþ Ce NTe TdS KE ¼

ð6Þ

r

e

Eqs. (5) and (6) are the micro scale finite element formulations, which can be used to construct the multiscale base functions. At the macro scale, a composite plate is discretized into the macroscopic finite elements. It is assumed that the displacement vector uE of macroscopic nodes within a macroscopic element satisfy

ae ¼ NE uE

ð7Þ

where N E is the multiscale base function matrix of the microscopic element. With Eqs. (3) and (7), the total potential energy P for all macroscopic elements can be expressed as:

P¼

XX1 Z E

ð2Þ

where N e is the base function of a microscopic element and ae is the displacement vector of microscopic nodes within a microscopic element. According to Eqs. (1) and (2), the total potential energy PE for all microscopic elements can be expressed as:

ð4Þ



2

e

Xe

ðNE uE ÞT B0T

XXZ E

Xe

e



A B B D



 B0 NE uE dxdy

 XX(Z ðNE uE ÞT NTe
f dxdy  E

e

ð3Þ where A, B and D are the tension stiffness matrix, tension-bending coupling stiffness matrix and bending stiffness matrix for the composite plates respectively. B0 is the strain–displacement matrix.

Cer

ð8Þ By constructing a transformation matrix G between the macroscopic nodes within a macroscopic element and the macroscopic nodes within all the macroscopic elements, the displacement vector uE can be expressed as

uE ¼ Gu

Cer

)
ðNE uE ÞT NTe TdS

ð9Þ

where u is the displacement vector of macroscopic nodes within all the macroscopic elements. With Eqs. (8) and (9), the macro scale finite element formulation is derived as Eq. (10) based on the minimum potential energy principle:

Ku ¼ P

ð10Þ

where

( )   X T X TR A B 0 G NE Xe B0T B dxdyNE G B D e E ( ) X T X TR R
P¼ G NE Xe NTe
f dxdyþN TE Ce NTe TdS K¼

E

Fig. 2. Microscopic finite element model: (a) microscopic finite elements; (b) layup configuration.

e

ð11Þ

r

Eqs. (10) and (11) are the macro scale finite element formulations, which can be used to obtain the macroscopic displacement fields. Once the macroscopic displacement fields have been obtained, the microscopic displacement fields can be obtained by Eq. (7) and the microscopic stress and strain fields can be obtained easily by the downscale computation as follows:

425

M. Ren et al. / Composite Structures 160 (2017) 422–434

ee ¼ B0 ae ¼ B0 NE uE rke ¼ Q k ðe0e þ zk je Þ

ð12Þ

where Q k is the stiffness matrix of the kth layer, and

ee ¼ f e0e je g ¼ f e0x e0y c0xy jx jy jxy g rke ¼f rkx rky skxy g It can be seen that the multiscale base functions are employed to establish the relationship between the macro scale and the micro scale computation and they are the keys to the extended multiscale finite element method.

8 4 4 4 4 4 X X X X X > > > Nuui uEi þ Nuv i v Ei þ Nuwi wEi þ Nuxi hExi þ Nuyi hEyi >u ¼ > > > i¼1 i¼1 i¼1 i¼1 i¼1 > > > 4 4 4 4 4 > X X X X X > > > v ¼ Nv ui uEi þ Nvv i v Ei þ Nv wi wEi þ Nv xi hExi þ Nv yi hEyi > > > > i¼1 i¼1 i¼1 i¼1 i¼1 > > > < 4 4 4 4 4 X X X X X w¼ Nwui uEi þ Nwv i v Ei þ Nwwi wEi þ N wxi hExi þ Nwyi hEyi > > i¼1 i¼1 i¼1 i¼1 i¼1 > > > 4 4 4 4 4 > X X X X X > > E E E E > hx ¼ Nxui ui þ N xv i v i þ Nxwi wi þ Nxxi hxi þ Nxyi hEyi > > > > i¼1 i¼1 i¼1 i¼1 i¼1 > > > 4 4 4 4 4 > X X X X X > > E E E E > h ¼ N u þ N v þ N w þ N h þ Nyyi hEyi > y yui y v i ywi yxi xi i i i : i¼1

i¼1

2.3. Multiscale base functions

KF ¼ 0
E NE ¼ N

in X on Cu

ð13Þ


E is the displacement boundary condition. where N The common displacement boundary conditions include the linear displacement boundary conditions and the oscillatory displacement boundary conditions [35]. For the thin composite plates, the linear boundary conditions can be used for the in-plane displacements. While for deflections and rotations, new boundary conditions need to construct. The reason is that the displacement modes of deflections and rotations are generally nonlinear and coupled with each other. Taking a rectangular element as an example, the displacement modes of deflections and rotations are expressed as shown in Eq. (14).

8 4   X > > > w¼ Newwi wi þ Newxi hxi þ Newyi hyi > > > > i¼1 > > > < 4   @w X hx ¼ ¼ Nexwi wi þ Nexxi hxi þ Nexyi hyi > @y > i¼1 > > > > 4   > @w X > > > Neywi wi þ Neyxi hxi þ Neyyi hyi : hy ¼  @x ¼ i¼1

i¼1

i¼1

ð15Þ

The multiscale base functions are the displacement functions of the macroscopic elements which satisfy the equilibrium equations and displacement boundary conditions as shown in formula (13):



i¼1

ð14Þ

where N e is the base function of the microscopic elements. To construct the displacement boundary conditions for deflections and rotations, the displacement modes need to decouple as follows: when the multiscale base functions for deflections are constructed, the deflections of macroscopic nodes on the boundary satisfy the property of Kronecker delta, the rotations of macroscopic nodes on the boundary are zero and the deflections and rotations of microscopic nodes on the boundary are obtained by Eq. (14); Similarly, when the multiscale base functions for rotations are constructed, the rotations of macroscopic nodes on the boundary satisfy the property of Kronecker delta, the deflections of macroscopic nodes on the boundary are zero and the deflections and rotations of microscopic nodes on the boundary are obtained by Eq. (14). Then, the decoupling displacement boundary conditions for deflections and rotations can be constructed. With the decoupling displacement boundary conditions for deflections and rotations and the linear displacement boundary conditions for the in-plane displacements, the multiscale base functions can be obtained by solving Eqs. (5) and (13). Consider that the coupling effects are obvious as a result of layup configurations and coupled displacement modes between deflections and rotations for the thin composite plates. The displacements of a microscopic node within a macroscopic element can be expressed as

where N is the multiscale base function of the macroscopic element. Take an additional coupling terms Nuv i as an example. It means the displacement field u when unit displacement v is applied on the macroscopic node i. Eq. (15) can be expressed in a unified form as Eq. (7) for all microscopic nodes within a microscopic element. Take N w1 ¼ fN uw1 ; N v w1 ; N ww1 ; N xw1 ; N yw1 gT as an example to show the construction process of the multiscale base functions. N w1 means the displacement fields u, v , w, hx and hy within a macroscopic element when a unit deflection is applied on the macroscopic node 1. In other words, the construction of the multiscale base functions N w 1 corresponds to the same decoupling displacement boundary conditions. The decoupling displacement boundary conditions of macroscopic nodes on the macroscopic element boundary can be expressed as follows:

8
j ¼0 N > > > uw1 i >
> > < Nv w1 ji ¼ 0
Nww1 ji ¼ d1m > > >
xw1 j ¼ 0 > N > i > :
Nyw1 ji ¼ 0

ði ¼ 1; 2; 3; 4Þ

ð16Þ

The decoupling displacement boundary conditions of microscopic nodes on the macroscopic element boundary can be expressed as follows:

8 4  X  > >
uw1 ¼
uw1  > N Neui N > > i > > i¼1 > > > 4  > X  > >
v w1 
v w1 ¼ > Nev i N N > > i > > i¼1 > > > < 4  X e 
ww1 
ww1 ¼ Nwwi N N i > > i¼1 > > > 4 >   X e  > >
xw1 
> ¼ Nxwi N N > > i > xw1 > i¼1 > > > 4  > X >  >

yw1  > Neywi N > : Nyw1 ¼ i

ð17Þ

i¼1

Eqs. (16) and (17) are the decoupling displacement boundary conditions corresponding to the multiscale base functions N w 1. Under these decoupling displacement boundary conditions, the multiscale base functions N w 1 can be constructed by solving Eq. (13). The other multiscale base functions can be constructed in a similar way. 2.4. Oversampling technique For the thin composite plates, the decoupling displacement boundary conditions will impose restrictions on the deformations

426

M. Ren et al. / Composite Structures 160 (2017) 422–434

of macroscopic elements. By adopting the oversampling technique, the relaxed decoupling displacement boundary conditions can be obtained which are more suitable to simulate the boundary deformations of macroscopic elements, especially those with aperiodic microstructure characteristics. Consider a lager domain E0 that covers the original macroscopic element E as shown in Fig. 3. By applying the decoupling displacement boundary conditions to the oversampling element, the temporary base functions N 0 can be constructed. Then, the base functions N can be constructed by the linear combination of N 0 with the coefficient matrix C. Since the coupling effects are strong for composite laminates, especially asymmetric laminates, the coupling terms within the equations between N and N 0 should not be neglected. Taking the base function N iww as an example, it can be expressed as follows

Niww ¼

4 4 4 4 X X X X v 0 0 0 0 C wu Cw C ww C wx ij N jwu þ ij N jww þ ij N jwhx ij N jwv þ j¼1

þ

j¼1

j¼1

4 X

0 C wy ij N jwhy

j¼1

ð18Þ

j¼1

where C wu , C wv , C wx and C wy are the coupling terms. It is important to emphasize that the coupling terms can’t be omitted. It is not enough to obtain the coefficient matrix C by the condition N iww jj ¼ dij alone. In order to obtain the coefficient matrix C, it is necessary to construct all the other base functions to build the simultaneous equations as follows.

N ¼ CN0

ð19Þ

where

2

Nuu 6N 6 vu 6 N¼6 6 Nwu 6 4 Nxu Nyu 2 uu C 6 vu 6C 6 wu C¼6 6C 6 xu 4C

N uv

Nuw

Nux

Nvv

Nv w

Nv x

N wv

Nww

Nwx

N xv N yv

Nxw Nyw

Nxx Nyx

C uv C

vv

C wv C

yu

xv

yv

C uw C

vw

C ux C

vx

C ww

C wx

xw

xx

C

C

Nuy

3

Nv y 7 7 7 Nwy 7 7; 7 Nxy 5 Nyy 3 C uy 7 Cvy 7 7 wy 7 C 7; 7 C xy 5

yw

C C C C yx C yy 3 2 0 0 0 Nuu Nuv N uw N0ux N0uy 7 6 0 6 N v u N0vv N0v w N0v x N0v y 7 7 6 6 0 0 0 0 0 7 N 0 ¼ 6 Nwu Nwv Nww Nwx Nwy 7 7 6 0 0 0 0 0 7 6N 4 xu Nxv N xw Nxx Nxy 5 N0yu N0yv N0yw N0yx N0yy

Fig. 3. Illustration of the oversampling technique.

which is composed of 10  10 macroscopic elements (see Fig. 4a). The length of plate is 100 mm and the thickness is 1 mm with layup [0/90/0]. Each layer possesses the same thickness. The lamina constants are given: E11 ¼ 25E22 , G12 ¼ 0:5E22 , m12 ¼0:25. Three kinds of microscopic FE models are considered, where the numbers of microscopic elements within a macroscopic element are 1  1, 5  5 and 10  10 (see Fig. 4b). The three multiscale numerical models are designated as EMsFEM_1  1, EMsFEM_5  5 and EMsFEM_10  10 respectively. The nondimensionalized center deflection and stress can be expressed as:

8 3 2 4 >
> < x¼xðE2 h =a q0 Þ10 2 2 r
xx ¼rxx ðh =a q0 Þ > > :
ryy ¼ryy ðh2 =a2 q0 Þ

ð20Þ

The results obtained by the method developed in this paper and the Navier method [40] are presented in Table 1 respectively. The results obtained by the Navier method are employed as reference solutions. The relative error calculation formula is defined as follows:

d¼

jREMsFEM  RNav uer j  100% jRNav uer j

ð21Þ

where REMsFEM are the results obtained by the developed method, RNav ier are the results obtained by the Navier method. By comparing the results of the two methods, we can see that: (1) The results obtained by the developed method are in good agreement with the analytical solutions. It shows that the developed method has high computing accuracy.

All the terms of N, C and N 0 are the n-order matrices, n is the node number of a macroscopic element. According to the property of Kronecker delta, the coefficient matrix C can be determined by the condition N ¼ I. I is a unit matrix. Then, the relaxed decoupling displacement boundary conditions can be obtained which can be used to construct the final multiscale base functions by the method mentioned in Section 2.3. 3. Validation In order to validate the computing accuracy of the developed method, two groups of examples are considered in this section and the results are compared with those in the literatures [40,41]. Example 3.1 Consider a thin square composite plate with all edges simply supported under the uniform distributed load q0 ,

Fig. 4. Multiscale finite element model of a thin square composite plate: (a) macroscopic finite element model; (b) three kinds of microscopic finite element models.

427

M. Ren et al. / Composite Structures 160 (2017) 422–434 Table 1 The nondimensionalized center deflections and stresses. Variables


x r
xx r
yy

The Navier method [40]

The developed method

Results

11

0.6660 0.8075 0.0306

55

10  10

Results

Errors (%)

Results

Errors (%)

Results

Errors (%)

0.6778 0.8199 0.0307

1.77 1.54 0.33

0.6778 0.8254 0.0307

1.77 2.22 0.33

0.6778 0.8261 0.0306

1.77 2.30 0

(2) With the refinement of the microscopic FE models, the results obtained by the developed method tend to be stable. It means the method has a good convergence. (3) For the incompatible rectangular elements of thin composite laminates, the convergence may be non-monotonic. In other words, the results obtained by the developed method may be the upper or lower bounds of the reference solutions. Example 3.2 Consider a beam-like composite laminate that is composed of 20  4 macroscopic elements (see Fig. 5a). The microscopic FE model is considered as Fig. 5b, in which the number of microscopic elements within a macroscopic element is 20  20 and the size of microscopic elements is similar to the examples in the literature [41]. This multiscale numerical model is designated as EMsFEM_20  20. The laminate is 40 mm thickness with layup [0/90/0/90]s. Each layer is 5 mm thickness. The laminate is clamped at x = 0 and subjected to 55.56 N/mm uniform distributed linear load in x-direction and 55.56 N/mm uniform distributed linear load in z-direction at x = 1000 mm. Note that the uniform distributed linear load is equivalent to the point force mentioned in the literature [41]. The lamina constants are given: E11 ¼ 110:50 GPa, E22 ¼ 13:64 GPa, G12 ¼ 3:26 GPa, m12 ¼ 0:329. The macroscopic displacements and microscopic stresses obtained by the EMsFEM are compared with the results calculated by Mechanics of Structure Genome (MSG) and 3D finite element analysis (3D-FEA) in the literature [41]. As Figs. 6 and 7 show that the macroscopic displacements u and w along x-direction on the central line (y = 90 mm) of the middle plane (z = 0) and microscopic stresses components rxx and ryy along z-direction on the central line (y = 90 mm) of the middle cross section (x = 500 mm) obtained by the EMsFEM agree well with those calculated by the above methods. Besides, the computing time for the EMsFEM is 103 s including the time taken by micro scale, macro scale and downscale computations, which is far less than the computing time for 3D-FEA. 4. Numerical examples In this section, two groups of examples are designed to discuss the effects of aperiodic microstructures and layup configurations on the accuracy of the developed method. In order to verify the

Fig. 6. Displacement of the beam-like composite laminate.

accuracy of the developed method, the results obtained by the conventional FEM on the fine model (FEM_F) are employed as reference solutions. Example 4.1 Consider a thin composite plate which is composed of 15  3 macroscopic elements (see Fig. 8a). For the macroscopic element E, three kinds of microscopic holes with different radius are designed (see Fig. 8b–d). According to whether the oversampling technique is used to construct the relaxed boundary conditions for the macroscopic element E, the multiscale finite element models are designated as EMsFEM and EMsFEM_O. The oversampling element E0 is shown in Fig. 8a. The composite plate is 2 mm thickness with cross-ply layup [02/902]2. Each layer is 0.25 mm thickness. The left side is clamped and the right side is subjected to a uniform distributed transverse linear load F = 103 N/mm. Three kinds of microscopic finite element models are shown as Fig. 9. For the other macroscopic elements without holes, the microscopic finite element models are the same as microscopic finite element model 1 (see Fig. 9a). Take composite as the anisotropic material. The properties of lamina are shown in Table 2 [42]: The microscopic nodal displacements of middle surface along the positive direction of X axis (Y = 15 mm) are shown in Figs. 10–12:

Fig. 5. Multiscale finite element model of a beam-like composite laminate: (a) macroscopic finite element model; (b) microscopic finite element model.

428

M. Ren et al. / Composite Structures 160 (2017) 422–434 Table 2 Material properties of T300/BSL914 epoxy [42]. Variable

Value

Longitudinal modulus E11 Transverse modulus E22 In-plane shearing modulus G12 Poisson ratio m12

138 GPa 11 GPa 5.5 GPa 0.28

Fig. 10. Microscopic nodal displacements u and w of the composite plate with a hole (microscopic model 1).

Figs. 10–12 illustrate that, for the three kinds of microscopic models, both the displacements u and w obtained by the EMsFEM and EMsFEM_O agree well with those of the FEM_F. It indicates that the method developed in this paper is of high computing accuracy. Consider that there are slight errors near the holes, so the relative errors within the macroscopic element E are analyzed. The relative error calculation formula is defined as follows: Fig. 7. Stress distribution along the thickness of the beam-like composite laminate: (a) rx ; (b) ry .

d¼

jREMsFEM  RFEM F j  100% jRFEM F j

ð22Þ

Fig. 8. Macroscopic finite element model with three kinds of microscopic models: (a) macroscopic finite element model; (b) microscopic model 1; (c) microscopic model 2; (d) microscopic model 3.

Fig. 9. Microscopic finite element models: (a) microscopic finite element model 1; (b) microscopic finite element model 2; (c) microscopic finite element model 3.

M. Ren et al. / Composite Structures 160 (2017) 422–434

Fig. 11. Microscopic nodal displacements u and w of the composite plate with a hole (microscopic model 2).

Fig. 12. Microscopic nodal displacements u and w of the composite plate with a hole (microscopic model 3).

where REMsFEM is the result obtained by the developed method, RFEM F is the result obtained by the FEM. According to formula (22), the error distributions of displacement w within the macroscopic element E are shown in Figs. 13–15. It is known from Figs. 13–15 that the maximum errors of macroscopic displacement w obtained by the EMsFEM and EMsFEM_O are all below 0.04%. It means that no matter the oversampling technique is used or not, the developed method can simulate the deformations accurately.

429

According to formula (22), the error distributions of microscopic stress rx on the first layer within the macroscopic element E are shown in Figs. 16–18. Figs 16–18 indicate that microscopic stresses obtained by EMsFEM_O are more accurate than those obtained by EMsFEM for all models. It means that the relaxed boundary conditions obtained by the oversampling technique can simulate the boundary deformations of the macroscopic elements with holes more reasonably. Besides, the reason that the maximum errors increase with the radius of the hole increasing is because the deformations of the macroscopic elements will be more complicated and the displacement boundary conditions will impose more restrictions on the deformations of macroscopic elements. Thus, the developed method is more applicable to the composite structures with multiscale features obviously. Example 4.2 Consider a thin rectangular composite plate with quasi-isotropic layup [0/±45/90]s. Each layer is 0.25 mm thickness. The geometry information of the plate are the same as example 4.2 in sec. 4. The left side is clamped and the right side is subjected to a combined load of tension and bending. Consider that the tension rigidity is higher than the bending rigidity, so the uniform distributed tension linear load is F1 = 10 N/mm and the uniform distributed transverse linear load is F2 = 103N/mm. The composite plate is composed of 15  3 macroscopic elements and there are three microscopic holes distributed in the macroscopic elements as shown in Fig. 19. The microscopic finite element models for the macroscopic elements E1 , E2 and E3 are all the same as microscopic finite element model 2 in Fig. 9b and the others are all the same as microscopic finite element model 1 in Fig. 9a. The properties are shown in Table 2. The microscopic nodal displacements of middle surface along the positive direction of X axis (Y = 15 mm) are shown in Fig. 20. Fig. 20 illustrates that both the displacements u and w obtained by the EMsFEM and EMsFEM_O agree well with those of the FEM_F. It indicates that the developed method can simulate the deformations of composite plates with multi-holes at the macro scale accurately. According to formula (22), the error distributions of microscopic stress rx on the first layer within the macroscopic elements E1 , E2 and E3 are shown in Figs. 21–23. Figs 21–23 indicate that the microscopic stresses obtained by the EMsFEM_O are also more accurate than those obtained by the EMsFEM for all the macroscopic elements with holes and the maximum errors are all nearly the same. Although the maximum errors are higher than those of Example 1, the overall accuracy can still meet the requirements for the EMsFEM_O.

Fig. 13. The error distributions of displacement w within the macroscopic element E (microscopic model 1): (a) EMsFEM; (b) EMsFEM_O.

430

M. Ren et al. / Composite Structures 160 (2017) 422–434

Fig. 14. The error distributions of displacement w within the macroscopic element E (microscopic model 2): (a) EMsFEM; (b) EMsFEM_O.

Fig. 15. The error distributions of displacement w within the macroscopic element E (microscopic model 3): (a) EMsFEM; b EMsFEM_O.

Fig. 16. The error distributions of microscopic stress

rx within the macroscopic element E (microscopic model 1): (a) EMsFEM; (b) EMsFEM_O.

5. Computer memory and computing time In addition to the high computing accuracy, the computing efficiency of the method developed in this paper is also significantly higher than the conventional finite element method. In our method, the structures can be solved at the macro scale rather than micro scale on the basis of the multiscale base functions which are employed to bring the microstructure characteristics into the macroscopic stiff matrices and force vectors. Thus, the computer

memory is reduced and the computing efficiency is improved significantly. In order to analyze the computing efficiency of the developed method, the computer memory and computing time are estimated roughly and compared with the conventional finite element method as follows. Since the computer memory mainly depends on the storage of the stiff matrices, in the case of the full storage of the stiff matrices, the computer memory needed in our method is equivalent to:

M. Ren et al. / Composite Structures 160 (2017) 422–434

Fig. 17. The error distributions of microscopic stress

rx within the macroscopic element E (microscopic model 2): (a) EMsFEM; (b) EMsFEM_O.

Fig. 18. The error distributions of microscopic stress

rx within the macroscopic element E (microscopic model 3): (a) EMsFEM; (b) EMsFEM_O.

431

Fig. 19. Macroscopic finite element model with three holes.

U ¼ ðn  Nmacro Þ2 þ ðn  Nmicro Þ2  Mmicro

ð23Þ

where U is the computer memory, Nmacro is the total number of macroscopic nodes, N micro is the number of microscopic nodes in a macroscopic element, Mmicro is the total number of macroscopic elements and n is the number of degrees of freedom on a single node. If the oversampling technique is adopted, the computer memory needed in our method is equivalent to:

U ¼ ðn  Nmacro Þ2 þðn  Nmicro Þ2  Mmicro þ ðn  NO Þ2  MO

ð24Þ

where M O is the number of oversampling elements and N O is the number of microscopic nodes in a oversampling element. The memory needed in the FEM is equivalent to:

U ¼ ðn  NÞ2 Fig. 20. Microscopic nodal displacements u and w of the composite plate with three holes.

where N is the total number of nodes.

ð25Þ

432

M. Ren et al. / Composite Structures 160 (2017) 422–434

Fig. 21. The error distributions of microscopic stress

rx within the macroscopic element E1 : (a) EMsFEM; (b) EMsFEM_O.

Fig. 22. The error distributions of microscopic stress

rx within the macroscopic element E2 : (a) EMsFEM; (b) EMsFEM_O.

Fig. 23. The error distributions of microscopic stress

rx within the macroscopic element E3 : (a) EMsFEM; (b) EMsFEM_O.

Take the care of Example 4.2 in Section 4. The composite plate is composed of 15  3 macroscopic elements at the macro scale and 40  40 microscopic elements at the micro scale. According to formulas (23)-(25), the computer memory needed in the FEM_F, EMsFEM and EMsFEM_O is shown in Fig. 24. Fig. 24 indicates that the computer memory needed in our method is far less than the FEM. The computer memory needed in the EMsFEM_O is slightly more than the EMsFEM owing to the construction of the relaxed boundary conditions.

The computing time can be roughly estimated on the basis of the same software and hardware environments. Take the care of Example 4.2 in Section 4. The computing time for different models is listed in Table 3. Note that for the same macroscopic elements, the multiscale base functions are constructed only once. Table 3 shows that the computing time needed for EMsFEM and EMsFEM_O is reduced dramatically compared with the FEM_F under the same software and hardware environments. It demonstrates that the developed method possesses high computing efficiency.

M. Ren et al. / Composite Structures 160 (2017) 422–434

433

References

Fig. 24. The comparison of the computer memory needed in different numerical models.

Table 3 Computing time needed for the different numerical models. Model Time (s)

EMsFEM 583

EMsFEM_O 1329

FEM_F 14626

Particularly, the micro scale computation employed to construct the multiscale base functions can execute independently. Then, the micro scale computation can be extended for parallel computing to further reduce computational cost. 6. Conclusion In this paper, a novel multiscale finite element method suitable for small-deflection analysis of thin composite plates with aperiodic macrostructure features is developed on the basis of the theoretical framework of the extended multiscale finite element method. Specifically: (1) the decoupling displacement boundary conditions and the relaxed decoupling displacement boundary conditions for deflections and rotations are reconstructed which can simulate the boundary deformations of the macroscopic elements accurately. Thus, the multiscale base functions employed to bring the microstructure characteristics into the macroscopic deformations are constructed accurately (2) by introducing the additional coupling terms among translations and rotations into the multiscale bases functions, the coupling effects of composite laminates, especially asymmetric composite laminates can be considered. The numerical examples illustrate that the developed method has the characteristic of high computing accuracy for thin composite plates in consideration of layup configurations and aperiodic microstructure characteristics, especially under the relaxed decoupling displacement boundary conditions. Meanwhile, the computer memory is reduced and the computing efficiency is improved significantly compared with the conventional finite element method. Consider that the microscopic response fields can be obtained accurately and efficiently. The developed method has great potential for multiscale damage analysis and strength prediction for composite structures. That will be carried out in our future work. Acknowledgements This work was supported by the National Natural Science Foundation of China (No. 11372058), the National Basic Research Program of China (Nos. 2014CB046506, 2014CB049000).

[1] Aboudi J, Arnold SM, Bednarcyk BA. Micromechanics of composite materials: a generalized multiscale analysis approach; 2013. [2] Huang C, Ren MF, Li T, Chang X, Cong J, Lei YJ. Trans-scale modelling framework for failure analysis of cryogenic composite tanks. Compos Part BEng 2015;85:41–9. [3] Chen JF, Morozov EV, Shankar K. Progressive failure analysis of perforated aluminium/CFRP fibre metal laminates using a combined elastoplastic damage model and including delamination effects. Compos Struct 2014;114(1):64–79. [4] Bogdanor MJ, Oskay C, Clay SB. Multiscale modeling of failure in composites under model parameter uncertainty. Comput Mech 2015;56(3):1–16. [5] Wang B, Tian K, Hao P, et al. Hybrid analysis and optimization of hierarchical stiffened plates based on asymptotic homogenization method. Compos Struct 2015;132(11):136–47. [6] Wang B, Tian K, Hao P, et al. Numerical-based smeared stiffener method for global buckling analysis of grid-stiffened composite cylindrical shells. Compos Struct 2016. [7] Karkkainen RL, Sankar BV. A direct micromechanics method for analysis of failure initiation of plain weave textile composites. Compos Sci Technol 2006;66(1):137–50. [8] Sun H, Di S, Zhang N, et al. Micromechanics of composite materials using multivariable finite element method and homogenization theory. Int J Solids Struct 2001;38(17):3007–20. [9] Bensoussan A, Lions JL, Papanicolaou G. Asymptotic analysis for periodic structures. North-Holland Pub. Co.; 1978. [10] Guedes J, Kikuchi N. Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Comput Method Appl Mech 1990;83(2):143–98. [11] Terada K, Kikuchi N. Nonlinear homogenization method for practical applications. Distrib Comput 1995;17(2):171–89. [12] Fish J, Shek K, Pandheeradi M, et al. Computational plasticity for composite structures based on mathematical homogenization: theory and practice. Comput Methods Appl Mech Eng 1997;148(1):53–73. [13] Feyel F, Chaboche JL. FE2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials. Comput Methods Appl Mech Eng 2000;183(3–4):309–30. [14] Lee K, Moorthy S, Ghosh S. Multiple scale computational model for damage in composite materials. Comput Methods Appl Mech Eng 1999;172(1– 4):175–201. [15] Jacob F, Yu Q, Kamlun S. Computational damage mechanics for composite materials based on mathematical homogenization. Int J Numer Methods Eng 1999;45(11):1657–79. [16] Ma J, Zhang J, Li L, et al. Random homogenization analysis for heterogeneous materials with full randomness and correlation in microstructure based on finite element method and Monte-carlo method. Comput Mech 2014;54 (6):1395–414. [17] Hazanov S, Huet C. Order relationships for boundary conditions effect in heterogeneous bodies smaller than the representative volume. J Mech Phys Solids 1994;42:1995–2011. [18] Kanit T, Forest S, Galliet I, et al. Determination of the representative volume for random composites: statistical and numerical approach. Int J Solids Struct 2003;40(13–14):3647–79. [19] Khisaeva ZF, Ostoja-Starzewski M. On the size of RVE in finite elasticity of random composites. J Elasticity 2006;85(2):153–73. [20] Babuška I, Osborn JE. Generalized finite element methods: their performance and their relation to mixed methods. Siam J Numer Anal 1983;20(3):510–36. [21] Babuska I, Caloz G, Osborn JE. Special finite-element methods for a class of 2nd-order elliptic problems with rough coefficients. Siam J Numer Anal 1994;31(4):945–81. [22] Hou TY, Wu XH. A multiscale finite element method for elliptic problems in composite materials and porous media. J Comput Phys 1997;134(1):169–89. [23] Hou TY, Wu XH, Cai Z. Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients. Math Comput 1999;68 (227):913–43. [24] Efendiev YR, Durlofsky LJ. Accurate subgrid models for two-phase flow in heterogeneous reservoirs. Spe J 2003;9(2):219–26. [25] Hou TY. Multiscale modelling and computation of fluid flow. Int J Numer Meth Fluid 2005;47(47):707–19. [26] Hughes TJR, Feijóo GR, Mazzei L, et al. The variational multiscale method-a paradigm for computational mechanics. Comput Methods Appl Mech Eng 1998;166(1–2):3–24. [27] Rudraraju S, Salvi A, Garikipati K, et al. Predictions of crack propagation using a variational multiscale approach and its application to fracture in laminated fiber reinforced composites. Compos Struct 2012;94(11):3336–46. [28] Ghosh S, Lee K, Raghavan P. A multi-level computational model for multi-scale damage analysis in composite and porous materials. Int J Solids Struct 2001;38 (14):2335–85. [29] Raghavan P, Ghosh S. Concurrent multi-scale analysis of elastic composites by a multi-level computational model. Comput Method Appl Mech 2004;193 (6):497–538. [30] Jenny P, Lee SH, Tchelepi HA. Multi-scale finite-volume method for elliptic problems in subsurface flow simulation. J Comput Phys 2003;187(1):47–67. [31] Hajibeygi H, Jenny P. Multiscale finite-volume method for parabolic problems arising from compressible multiphase flow in porous media. J Comput Phys 2009;228(14):5129–47.

434

M. Ren et al. / Composite Structures 160 (2017) 422–434

[32] Hajibeygi H, Karvounis D, Jenny P. A hierarchical fracture model for the iterative multiscale finite volume method. J Comput Phys 2011;230 (24):8729–43. [33] Calo VM, Efendiev Y, Galvis J, et al. Multiscale empirical interpolation for solving nonlinear PDEs. J Comput Phys 2014;278(1):204–20. [34] Chung E, Efendiev Y, Hou TY. Adaptive multiscale model reduction with generalized multiscale finite element methods. J Comput Phys 2016;320:69–95. [35] Zhang HW, Wu JK, Fu ZD. Extended multiscale finite element method for mechanical analysis of periodic lattice truss materials. Int J Multiscale Comput 2010;8. 597-61. [36] Zhang HW, Wu JK, Lv J. A new multiscale computational method for elastoplastic analysis of heterogeneous materials. Comput Mech 2012;49:149–69. [37] Triantafyllou SP, Chatzi EN. A hysteretic multiscale formulation for nonlinear dynamic analysis of composite materials. Comput Mech 2014;54:763–87.

[38] Zhang HW, Liu Y, Zhang S, Tao J, Wu JK, Chen BS. Extended multiscale finite element method: its basis and applications for mechanical analysis of heterogeneous materials. Comput Mech 2014;53:659–85. [39] Liu H, Sun X, Xu Y, et al. A hierarchical multilevel finite element method for mechanical analyses of periodical composite structures. Compos Struct 2015;131:115–27. [40] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. Florida: CRC; 2004. [41] Liu X, Yu W. A novel approach to analyze beam-like composite structures using mechanics of structure genome. Adv Eng Softw 2016;100:238–51. [42] Ozsoy MA, Uzuntarla M, Ozer M. Current status and future trend of researches on the strength of fiber-reinforced composite-a summary of the results from a ‘‘failure olympics”. Adv Math 2007;23:1–4.