Multiscale homogenization for nearly periodic structures

Composite Structures 153 (2016) 345–355

Contents lists available at ScienceDirect

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

Multiscale homogenization for nearly periodic structures Akinori Yoshimura a,b,⇑, Anthony M. Waas b, Yoshiyasu Hirano a a b

Aeronautical Technology Directorate, Japan Aerospace Exploration Agency, 6-13-1 Osawa, Mitaka, Tokyo 181-0015, Japan William E. Boeing Department of Aeronautics and Astronautics, University of Washington, Guggenheim Hall, Seattle, WA 98195, USA

a r t i c l e

i n f o

Article history: Received 2 February 2016 Accepted 1 June 2016 Available online 7 June 2016 Keywords: Finite element analysis Homogenization method Carbon fiber reinforced plastic Effect of defect Fiber waviness Perturbation method

a b s t r a c t In this paper, we propose a new numerical multiscale homogenization method for nearly periodic structures, where perfect periodicity is disturbed due to unintended geometrical imperfections. Geometrical imperfection is introduced as an input to the calculation, and the calculation is always done using the idealized perfect geometry. This has a distinct advantage in numerical calculations where the same mesh can be used for a series of imperfect geometries. The method is implemented using commercially available finite element code ABAQUS. The proposed method is validated by comparing the results to those of a conventional calculation. The results demonstrate that the proposed method can accurately capture the stress distribution when the amplitude of the imperfection is small. Finally, the method is applied to the analysis of the effect of fiber waviness on the mechanical properties of a low cost CFRP. The simulation results reveal that the fiber waviness has little effect on the macroscopic stiffness, but significantly degrades the uniaxial tensile strength. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Materials with microstructure require multiscale analysis methods in order to capture effects at different length scales. This is significant in a wide variety of engineering materials, including structural foams and carbon fiber reinforced plastic (CFRP) laminates. In many instances, the microstructure geometry contains symmetries and repeat units, which allows simplifying the resulting analysis through the choice of a representative repeat unit cell. This has lead to a vast literature that deals with the mechanics of periodic structures and the resulting homogenization analysis methods were pioneered by a number of researchers, including, Babuška [1], Bensoussan [2], Sanchez-Palencia [3], Bakhvalov and Panasenko [4], etc. In many structural materials with microstructure, it is not unusual to find internal microscopic imperfections including tow distortion, voids, or resin pockets. Low-cost CFRP manufacturing techniques such as vacuum-assisted resin transfer molding (VaRTM) [5–8] or Out-of-Autoclave (OoA) prepreg [9] have been studied to apply these techniques to large aerospace structures [10–14]. Compared to prepreg-autoclave process, composites produced using low-cost manufacturing techniques tend to contain more imperfections in the microstructure. ⇑ Corresponding author at: Aeronautical Technology Directorate, Japan Aerospace Exploration Agency, 6-13-1 Osawa, Mitaka, Tokyo 181-0015, Japan. E-mail address: [email protected] (A. Yoshimura). http://dx.doi.org/10.1016/j.compstruct.2016.06.002 0263-8223/Ó 2016 Elsevier Ltd. All rights reserved.

These microscopic imperfections significantly affect the macroscopic material properties, such as stiffness and strength. A number of researchers have conducted studies on the compressive strength of polymer composites [15–21]. These studies have revealed that initial fiber misalignment has significant effect on the failure process and compressive strength of composites, which includes kink banding and fiber–matrix splitting instabilities. Moreover, Cox et al. [22] and Mouritz and Cox [23] revealed that fiber misalignment can have significant effect on the tensile failure process and strength of 3D woven CFRP and stitched CFRP. For VaRTM specimens, Hirano et al. [12] compared the macroscopic material properties between VaRTM and conventional prepreg-autoclave CFRP specimens. Compressive strengths of VaRTM specimens were significantly lower than those of the conventional CFRP laminates. For OoA CFRP, Agius and Fox [24] investigated the relationship between void content and fracture toughness, and revealed that mode-II interlaminar fracture toughness was significantly increased by reducing void content. As described above, microscopic geometrical imperfections in the composites significantly affect macroscopic material properties. Therefore, observation and evaluation techniques of microscopic imperfections and their effects on macroscopic material properties are very important for characterizing low cost CFRP as they are used in large aerospace structures. Past studies have reported results on the observation of microscopic imperfection by utilizing micro-focus X-ray computed tomography (CT). Schell et al. [25] produced 3D images of glass fiber reinforced plastic woven composites from 3D X-ray CT images.

346

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355

Djukic et al. [26–28] applied contrast enhancement methods to CFRP and measured tow architecture of the woven materials. Iida et al. [29] also applied a similar method to measure the distribution of fiber orientation in short-fiber CFRP. Centea and Hubert [30] acquired X-ray CT images for several out-of-autoclave CFRP samples which were partially processed, and investigated the porosity in the material during out-of-autoclave curing process. Scott et al. [31] employed micro-focus X-ray CT together with synchrotron radiation X-ray CT, and observed microscopic voids and damage processes of Aluminum-CFRP-GFRP structure. Yoshimura et al. [32] applied image processing technique to the micro-focus X-ray CT images and measured microscopic fiber orientation and its deviation in CFRP. In summary, detailed observation of the microscopic imperfection in composite materials are possible because of advances in micro-focus X-ray CT devices in recent years, together with the innovative methods developed by researchers. It is well recognized that microstructure drives the performance of a material. Thus, once the state of the microstructure is identified, analysis methods need to be developed to understand and predict the effect of the microstructure on performance. An effective method is finite element analysis using a three-dimensional finite element model constructed from 3D micro-focus X-ray CT images. Naouar et al. [33] proposed a method to produce a meso-scale finite element model directly from micro-focus X-ray CT image and conducted nonlinear simulations. Czabaj et al. [34] proposed image processing algorithm for detecting the positions of the fiber in the CT image, and applied it to small CFRP sample. They then constructed microscale (fiber and matrix) finite element model which contains realistic fiber orientation distortion. These methods can introduce realistic microscopic geometrical imperfection into the finite element model. However, in order to produce the finite element mesh directly from 3D X-ray CT image, many nodes and elements are necessary so as to trace the details of the microscopic imperfection. This can lead to large computational costs. Additionally, the mesh generation algorithm might produce distorted elements because of the CT data. In order to avoid such distorted elements, the meshing operation becomes a time-consuming process. When considering the structural simulation of a large aerospace structure, the imperfection pattern is different in each part of the structure. Therefore, it is necessary to conduct multiple calculations which include many imperfection patterns. However, it is very difficult to conduct multiple analyses based on these methods because of difficulty in meshing and the resulting large degrees of freedom (DoF). In this paper, we propose a method to overcome these issues by developing a calculation method in which a single idealized mesh is used, and each imperfection pattern is introduced to the model as input. The theoretical developments that are needed to facilitate this method are presented, based on the classical homogenization method for periodic structures. A new numerical simulation method is developed for a plate structure, in which unintended initial geometrical imperfections are considered. The method is implemented using the finite element (FE) method and a regular perturbation method. The geometrical imperfection is introduced as input for the calculation, and the calculation is always done using an idealized regular mesh. For large aerospace structural simulations, multiscale analysis is necessary because the characteristic length of the structure is much larger than the microstructure of the material. A number of researches have reported multiscale simulation methods for composite structures [35–40], which are relevant to the present study. The paper is organized as follows: In Section 2, formulation of new method is introduced. In Section 3 the proposed method is validated by comparing the results of the proposed method to the conventional calculation method. In Section 4, the proposed method is applied to the analysis of the effect of the imperfection

on the properties of VaRTM CFRP laminate. Conclusions from the present study are reported in the final section. 2. Formulation 2.1. Problem setting Fig. 1 shows the schematic of the problem considered in this study. Consider a macroscopic plate X that includes an in-plane periodic microscopic inhomogeneity. It is assumed that deformation in X can be described by Kirchhoff–Love’s plate theory. Let Y denote the representative volume (RUC; Representative unit cell) of the in-plane periodic inhomogenity. Geometrical imperfection is considered only in the RUC. Orthogonal coordinate system O
X 1 X 2 X 3 is used in order to indicate reference configuration of domain X. On the other hand, another orthogonal coordinate system O0
Y 1 Y 2 Y 3 is used for reference configuration of the RUC. Note that basis vectors of these two coordinate systems are respectively parallel. Scales of X 1 and X 2 are much larger than those of Y 1 and Y 2 . The scale of thickness direction X 3 is the same as that of Y 3 . In the macroscopic domain, the current configuration xi and reference configuration X i are related by,

xi ¼ X i þ uGi ðX Þ

ð1Þ

uGi

where denotes macroscopic displacement, which depends on macroscopic position. ^ Li ðY Þ is In the microscopic domain Y, geometrical imperfection eu considered. Scalar parameter e is used as perturbation parameter in the latter part of the paper. In order to consider geometrical imperfection, we used three configurations: reference configuration, imperfect configuration and current configuration (see Fig. 2). In the finite element analysis, reference configuration corresponds to regular mesh, and imperfect configuration corresponds to a distorted mesh. As described in Section 1, in order to reduce meshing cost, formulation using reference configuration is preferred. Positions in reference configuration, imperfect configura i and tion, and current configuration are expressed by vector Y i ; Y yi , respectively. These vectors are related as,

^ Li ðY Þ Y i ¼ Y i þ eu

ð2Þ

yi ¼ Y i þ uLi ðY Þ

ð3Þ

uLi

where denotes microscopic displacement, which depends on microscopic position. Because of the in-plane periodicity of Y, both ^ Li and uLi has in-plane periodicity. u 2.2. Strain and stress tensor Infinitesimal mechanical strain tensor is calculated as,

eij

@uGj 1 @uGi ¼ þ 2 @X j @X i

!

 1 0   ^ Lj @ uLj
eu ^ Li 1 @@ uLi
eu A þ þ 2 @ Y i @ Y j

Fig. 1. Schematic showing the problem considered.

ð4Þ

347

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355

Fig. 2. Three configurations used in this study.

It should be noted that the mechanical strain is caused by the deformation between current configuration and imperfect configuration. Thus microscopic displacement must be differentiated with respect  Because partial differentiation with to imperfect configuration Y. respect to reference configuration is preferred, chain rule is considered as,

    ^ Lk ^ Lk @u @ @ @Y k @ @ Y k
eu @ ¼ ¼ ¼ d
e ik @Y k @ Y i @ Y i @ Y i @Y k @ Y i @Y k ^ Lk @ @ @u
e ¼ @Y i @Y k @ Y i

ð5Þ

ð6Þ

Note that only the 0th and 1st order terms are considered in the rest of this paper because terms higher than 2nd order have a negligible effect. By substituting Eq. (6) into Eq. (4), we get

eij

@uGj 1 @uGi ¼ þ 2 @X j @X i

!

L 1 @uLi @uj þ þ 2 @Y j @Y i ! ^ Lk @uLi @ u ^ Lk @uLj 1 @u þ
e 2 @Y j @Y k @Y i @Y k

!

^ Lj ^ Li @ u 1 @u þ
e 2 @Y j @Y i

!

ð7Þ

uL0 i

ð8Þ

where and denote 0th and 1st order components of the solution, respectively. By substituting Eq. (8) into Eq. (7), the following perturbation series for the strain tensor is obtained.

eij  e þ e G ij

L0 ij

  C ^L þ e eL1 ij
eij
eij

eGab ¼ E0ab þ Y 3 K ab eGi3 ¼ 0

ð12Þ

Eab

@uG0 1 @uG0 b a ¼ þ 2 @X b @X a

K ab

@ 2 uG0 3 ¼
@X a @X b

0

! ð13Þ

Here, E0ab corresponds to macroscopic membrane strain and K ab corresponds to macroscopic curvature. Note that Y 3 instead of X 3 is used because scale of X 3 is the same as that of Y 3 . Therefore Eq. (9) is written as, L1 C ^L eab ¼ E0ab þ Y 3 K ab þ eL0 ab þ e eab
eab
eab



L1 C ^L ei3 ¼ eL0 i3 þ e ei3
ei3
ei3

uL1 i

ð11Þ

where Greek subscripts are used for in-plane quantities, taking, values 1 and 2. Superscript 0 in Eq. (11) denotes mid-surface entities. Therefore, macroscopic strain eGij is computed as,



Here we assume that the solution (microscopic displacement uLi ) is expressed by the following perturbation series. L1 2 L2 L0 L1 uLi ¼ uL0 i þ eui þ e ui þ     ui þ eui

@uG0 3 ðX 1 ; X 2 Þ @X a

uG3 ðX 1 ; X 2 ; X 3 Þ ¼ uG0 3 ðX 1 ; X 2 Þ

where,

By substituting Eq. (5) into itself recursively, we acquire

^ Lk ^L ^L ^ Lk @ @ @ @u @ @ @u 2 @ @ ul @ uk ¼
e þ e
   
e @Y k @Y i @Y l @Y k @Y i @Y i @Y k @Y i @ Y i @Y i

uGa ðX 1 ; X 2 ; X 3 Þ ¼ uG0 a ðX 1 ; X 2 Þ
X 3



ð14Þ



It is assumed that materials considered are linear elastic. Then, the Cauchy stress tensor rij is calculated as,









L1 C ^L rij ¼ cijkl ekl ¼ cijab E0ab þ Y 3 K ab þ cijkl eL0 kl þ e ekl
ekl
ekl

 ð15Þ

where cijkl denotes the elasticity tensor.

ð9Þ 2.3. Principle of virtual work

where G ij

e

L0 ij

e

@uGj 1 @uGi ¼ þ 2 @X j @X i

!

@uL0 1 @uL0 j i ¼ þ 2 @Y j @Y i @uL1 j

1 @uL1 i þ 2 @Y j @Y i ! ^ Lj ^ Li @ u 1 @u ^eLij ¼ þ 2 @Y j @Y i

eL1 ij ¼

C ij

e

The principle of virtual work when applied over a representative unit cell (RUC) in equilibrium, results in,

Z

!

Y

!

^ Lk @uL0 ^ L @uL0 @u 1 @u j i ¼ þ k 2 @Y j @Y k @Y i @Y k

ð10Þ

!

As stated in Section 2.1, macroscopic deformation in X can be described by Kirchhoff–Love’s plate theory.

deLij rij dY ¼ 0

ð16Þ

where deLij denotes microscopic virtual strain. It should be noted that virtual work should be computed over the volume in the imperfect configuration. Because of the in-plane periodicity of the RUC, the right-hand-side of the Eq. (16) is equal to zero. As mentioned in Section 1, calculation with respect to reference configuration is preferred. The volume of integration in Eq. (16) can be transformed as,

Z

Y

deLij rij det jJjdY ¼ 0

ð17Þ

where det jJj denotes determinant of the Jacobian matrix, which is calculated as,

348

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355

@ Y

1

@Y 1

 det jJj ¼ det @@YY 2

1

@ Y

3





 @Y 2

@Y 3



@ Y 3

@ Y 1 @Y 2

@ Y 1 @Y 3

@ Y 2 @Y 2 @ Y 3 @Y 2

@Y 1

ð18Þ

@Y 3

2.4. Multiscale analysis

By considering the chain rule (Eq. (6)), Eq. (18) becomes,



^L ^L ^L



@u @u @u

1 þ e @Y1 e @Y1 e @Y13 1 2





@ u^L L L ^ ^ @ u @ u

2 2 det jJj ¼

e @Y2 1 þ e @Y 2 e @Y 3

1



L L L

@ u^3 ^ ^

@u @u 3 3

e @Y e 1 þ e @Y 2 @Y 3 1  L  L L ^1 @ u ^ ^ @u @u 1þe þ 2þ 3 @Y 1 @Y 2 @Y 3

The macroscopic stress resultant N ab and macroscopic moment resultant Mab are calculated as,

ð19Þ

deLij ¼

1 2

@ Y j

þ

@duLj

! ð20Þ

@ Y i

where duLi denotes microscopic virtual displacement. By substituting Eq. (6) into Eq. (20), virtual strain is written in the reference configuration as, C deLij ¼ deL0 ij
edeij

ð21Þ

L 1 @duLi @duj ¼ þ 2 @Y j @Y i

!

^ Lk @duLi 1 @u ¼ þ 2 @Y j @Y k @Y i @Y k

de

!

ð22Þ

By substituting Eqs. (15), (19), (21) into Eq. (17), results in,

Z Z   0 L0 L1 deL0 deL0 deL0 ij c ijab Eab þ Y 3 K ab dY þ ij c ijkl ekl dY þ e ij c ijkl ekl dY Y Y Y Z Z C ^L deL0
e deL0 ij c ijkl ekl dY
e ij cijkl ekl dY

Z

Y

Y

Z

Z   @u ^ Lp ^ Lp @u L0 þ e de E0ab þ Y 3 K ab dY þ e deL0 dY ij c ijkl ekl @Y p @Y p Y ZY Z  
e deCij cijab E0ab þ Y 3 K ab dY
e deCij cijkl eL0 kl dY ¼ 0 Y

ð23Þ

Y

Eq. (23) must be valid for any e including 0. Eq. (23) can be separated into 0th order and 1st order equations.

Y

Z   0 L0 deL0 c E þ Y K deL0 dY þ 3 ab ab ij ijab ij c ijkl ekl dY ¼ 0

ð24Þ

1 N0ab ¼  S

Z

Z L1 deL0 ij c ijkl ekl dY ¼

Y

Y

^L deL0 ij c ijkl ekl dY þ

C deL0 ij c ijkl ekl dY




Y

Y

ð25Þ

Y

The 0th order Eq. (24) is exactly the same equation as the equation that results from the classical homogenization method, and it corresponds to a perfect structure. The 0th order equation can be solved independently from the 1st order equation. By solving L0 the 0th order Eq. (24), uL0 i is first obtained. By substituting ui into right-hand-side terms of the 1st order Eq. (25), we then acquire uL1 . By using Eq. (8), the overall displacement vector is finally calculated.

ð29Þ

Y

1 M0ab ¼  S

Z Y

Z   Y 3 cijab E0ab þ Y 3 K ab dY þ Y 3 cijkl eL0 dY kl Y

(Z

Z   @u ^ Lp ^ Lp @u Y 3 cijab E0ab þ Y 3 K ab dY þ Y 3 cijkl eL0 dY kl @Y p @Y p Y Y Z Z   ^ L3 cijab E0ab þ Y 3 K ab dY
u ^ L3 cijkl eL0
u kl dY

1 M1ab ¼  S

Y

Z

Z Y

Z   cijab E0ab þ Y 3 K ab dY þ cijkl eL0 dY kl

Z   @u ^ Lp ^Lp @u cijab E0ab þ Y 3 K ab dY þ cijkl eL0 dY kl @Y p @Y p Y Y Z   C ^L þ cijkl eL1 kl
ekl
ekl dY

1 N1ab ¼  S

Y 3 cijkl Y

Z   @u ^ Lp ^Lp @u 0 L0 deL0 dY
deL0 dY ij c ijab Eab þ Y 3 K ab ij c ijkl ekl @Y p @Y p Y ZY Z   þ deCij cijab E0ab þ Y 3 K ab dY þ deCij cijkl eL0 kl dY Z

ð28Þ

(Z

þ

Z Y

where

Y

and

ð27Þ

By substituting (15) into Eq. (27), the following is obtained;

and

L0 ij c ijab

Z

Z Z Z ^L @u 1 1 1 Nab ¼  rab det jJjdY ¼  rab dY þ e  rab p dY @Y p S Y S Y S Y Z   1 L ^ 3 rab det jJjdY Y 3
eu M ab ¼  S Y Z Z Z ^ Lp @u 1 1 1 ^ L3 rab dY þ e Y 3 rab dY
e  Y r dY ¼ u 3 a b S Y @Y p S Y S Y

M ab ¼ M0ab þ eM1ab

^ Lk @duLj @u

C ij

ð26Þ

Nab ¼ N0ab þ eN1ab

where

deL0 ij

Z 1 N ab ¼  rab dY S Y Z 1 Y 3 rab dY M ab ¼  S Y

1
Y  2 plane. Eq. (26) where  S denotes area of the RUC projected to Y is re-written by using reference configuration as,

Virtual strain in Eq. (17) can be calculated as,

@duLi

In Appendix A, the proposed method is applied to a simple onedimensional problem in order to illustrate an analytical closedform solution without resorting to numerical implementation.



 C

^L eL1 kl
ekl
ekl dY



Y

ð30Þ

In the right-hand-side of the second equation of Eq. (29), first and second terms represent the effect of volume change, third term represents first order displacement. In the right-hand-side of the second equation of (30), first and second terms capture the effect of volume change, third and forth terms represent the effect of the displacement of the moment arms and the fifth term corresponds to the first order displacement. Fig. 3 shows overall flowchart of the multiscale simulation. The simulation was conducted based on the numerical material testing concept proposed by Terada et al. [39–41]. First of all, macroscopic stress responses were calculated for six basic deformation patterns. In this paper, Y 1 direction tension, Y 2 direction tension, Y 1 Y 2 direction shear, Y 1 direction bending, Y 2 direction bending and Y 1 Y 2 direction torsion were chosen for basic patterns. For example in Y 1 -tension calculation,

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355

349

solved. This procedure is called localization, while the former is referred to as homogenization. 3. Validation 3.1. Method In this section, the simulation method proposed in the previous section is validated. Fig. 4 shows procedure of the validation method. In this study, a regular FE mesh and a distorted FE mesh were prepared. 1% Y 1 direction uniaxial tensile macroscopic strain was applied for both models. In the distorted mesh model, geometrical imperfection was explicitly included. In other words, each ^ Li . For the regular node was displaced according to imperfection eu mesh model, the effect of imperfection is captured by the method described in the previous section. In the distorted mesh model, only one calculation was conducted and the obtained results were checked for convergence. Microscopic stresses calculated by the proposed method were compared to the results obtained with the distorted mesh calculations. Of interest is a quasi-isotropic 16-ply CFRP laminate model. Such a model is representative of many low-cost CFRP configurations currently in use. Overviews of the regular mesh model and distorted mesh model are illustrated in Figs. 5 and 6. The lamina stacking sequence is, [45/0/
45/90]2s. Dimensions of the model were 3.0 mm  3.0 mm  3.09 mm. Thickness of each ply tp was

Fig. 3. Flowchart of multiscale calculation.

8 0 9 8 9 E11 > > 1> > > > > > > 0 > > > > > > > > > > > > 0> E > > > 22 > > > > > > 0 = > > < <0> = E12 ¼ > 0> > K 11 > > > > > > > > > > > > > > > > > > > > > K 22 > > > >0> > > > > > : ; : ; 0 K 12

ð31Þ

are substituted to E0ab and K ab . Then by using Eqs. (24), (25) and (28)–(30), macroscopic stress resultant and moment resultant were calculated. By using the calculated resultants and moments for 6 patterns, macroscopic homogenized plate stiffness

8 9 2 N11 > A11 > > > > > 6A > > > > N 22 > > 6 21 > > > > 6 < N12 = 6 A61 ¼6 6 > > M 11 > > 6 C 11 > > > > 6 > > > > 4 C 21 M 22 > > > > : ; M 12 C 61

A12 A22

A16 A26

B11 B21

B12 B22

A62

A66

B61

B62

C 12

C 16

D11

D12

C 22

C 26

D21

D22

C 62

C 66

D61

D62

38 0 9 E11 > B16 > > > > > > 0 > > > > B26 7 > E22 > > 7> > > 7> < B66 7 E0 = 7 12 D16 7 > K 11 > > 7> > > > 7> > > > 5 D26 > > > K 22 > > > > : ; D66 K 12

ð32Þ

can be calculated. By using the calculated plate stiffness, the macroscopic problem can be solved. After that assigning the calculated E0ab and K ab for Eqs. (24), (25) and (15), microscopic stress distribution can be

Fig. 4. Procedure of validation calculation.

Fig. 5. Overview of the regular mesh model.

350

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355

Table 1 shows the material properties of the lamina. The material properties were decided based on the material test results for T800S/#3900-2B (Toray Industries, Inc.) measured by Japan Aerospace Exploration Agency (JAXA) [42]. For calculation, commercial general-purpose finite element simulation code, ABAQUS 6.14-2 was employed. From 1st to 4th terms in the right-hand-side of Eq. (17) were implemented by using user subroutine SIGINI. On the other hand, fifth and sixth terms in right-hand-side were implemented by using user subroutine UEL. The model consists of 47,089 nodes. 43,200 isoparametric brick solid elements (C3D8) were used for calculation. Each ply was divided by three elements in the thickness direction. In Y 1 and Y 2 direction, lengths of each element were 0.1 mm. 3.2. Results and discussion

Fig. 6. Overview of the distorted mesh model.

Table 1 Material properties of lamina. Longitudinal Young’s modulus E11 (GPa) Transverse Young’s modulus E22 (GPa) Through-the-thickness Young’s modulus E33 (GPa) In-plane Poisson’s ratio m12 Transverse Poisson’s ratio m23 Out-of-plane Poisson’s ratio m13 In-plane shear modulus G12 (GPa) Transverse shear modulus G23 (GPa) Out-of-plane shear modulus G13 (GPa)

153 8.0 8.0 0.34 0.45 0.34 4.03 3.4 4.03

0.193125 mm. In this validation, a sinusoidal shape imperfection was considered.

 L L L ^2 ; u ^3 ¼ ^1 ; u u

   tp 2p Y 1 0; 0; sin 4 3

ð33Þ

Fig. 7. Contour plot of

Figs. 7 and 8 plot contours of the Y 1 direction tensile stress r11 of the regular mesh model and the distorted mesh model. Displacement fields agree very well among the two models, and almost similar stress distributions are calculated. Tensile stress concentration occurred near the peaks of sinusoidal imperfection. Maximum stress of the regular mesh model was 2107.74 MPa, and minimum stress was 63.21 MPa. For the distorted mesh model, the maximum stress was 1947.12 MPa and minimum stress was 65.10 MPa. Differences between the two models were 9.4% for maximum stress and 4.5% for minimum stress. Table 2 summarizes the maximum and minimum values of all stress components. For all components, almost similar values were calculated for the regular mesh model and the distorted mesh model. Stress distributions were also calculated for sinusoidal imperfections of various amplitudes. Perturbation parameter was varied ^ Li was used as Eq. (33). Fig. 9 as e ¼ 0:25; 0:5; 1:5; 2:0, and the same u plots the relationship between amplitude of the imperfection and maximum, and minimum values of r11 . Because we considered only 0th and 1st order terms in perturbation method, maximum and minimum values of stress linearly changed for regular mesh model. On the other hand, maximum and minimum values changed nonlinearly for the distorted mesh model. We considered only 0th and 1st order terms in perturbation method. Differences between regular mesh model and distorted

r11 of regular mesh model.

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355

Fig. 8. Contour plot of

351

r11 of distorted mesh model.

Table 2 Maximum and minimum values of stress tensor. Regular mesh

r11 (MPa) r22 (MPa) r33 (MPa) s12 (MPa) s13 (MPa) s23 (MPa)

Distorted mesh

Maximum

Minimum

Maximum

Minimum

2107.74 526.54 29.70 446.90 75.15 37.40

63.21
2.16
29.70
492.46
75.15
37.40

1947.12 498.62 27.30 435.15 69.11 36.03

65.10
2.01
26.88
465.36
69.11
36.03

mesh model were caused by terms higher than 2nd order. By considering these higher order terms, accuracy of the result is improved. However, calculation of the higher order terms need computational time. Therefore in this method, there is a trade-off between accuracy and computational cost on the one hand, and accuracy and ease of making a mesh on the other. In many microstructures, lamina undulations (such as in textile composites, see [43]) result in areas of microstructure where the mesh can get distorted due to unintended fiber tow waviness. In those case, the present method becomes advantageous and cost effective because all computations are performed on the same perfect mesh. As demonstrated here, the proposed perturbation method can lead to very accurate stress distributions when the amplitude of the imperfection is small. As the amplitudes become larger higher order corrections as discussed must be included to improve accuracy. 4. Effect of imperfection on the material properties of VaRTM CFRP 4.1. Experiments In this section, the method proposed in the Section 2 is applied to the analysis of the effect of fiber waviness on the material properties of CFRP laminates. Tensile tests were conducted for two types of CFRP laminates: CFRP manufactured by conventional prepreg-autoclave process and CFRP manufactured by VaRTM process. Conventional prepreg-autoclave CFRP was made from T800S/ #3900-2B (Toray Industries, Inc.) prepreg. On the other hand,

Fig. 9. Relationship between amplitude of imperfection and accuracy of the result.

uni-directional (UD) T800SC non-crimp fabric and modified low viscosity epoxy resin system XNR6809-XNH6809 (Nagase Elecs., Co. Ltd.) were used for VaRTM CFRP.

352

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355

Fig. 11. Comparison of calculated macroscopic stiffness matrix (A matrix).

Fig. 10. Micrograph of the VaRTM CFRP specimen.

VaRTM CFRP was manufactured as follows. Dry non-crimp fabrics were lay-upped on tool plate and covered by vacuum bag and vacuum-tight seal tape. Uncured resin was then drawn into the preform under vacuum. Curing of VaRTM CFRP was done under the atmospheric ambient pressure at 120 °C, whereas curing of prepreg-autoclave CFRP was done under high pressure provided by autoclave at 180 °C. Fig. 10 shows the micrograph of the section of the VaRTM CFRP. No voids are observed in the material, but fiber waviness was induced because of the low pressure during curing. Fiber volume fractions of prepreg-autoclave specimen and VaRTM specimen were 55.0% and 56.5%, respectively. Stacking sequences of both laminates were [0]6. Both specimens were cut by diamond saw into rectangle shape. Specimen dimensions of prepreg-autoclave specimens were 254 mm in length, 10 mm in width, and 1.10 mm in thickness. Those of VaRTM specimens were 254 mm in length, 15 mm in width, and 1.19 mm in thickness. Uniaxial tensile tests were conducted for both types of CFRP. Test method was based on SACMA SRM-4R [44]. For both materials, crosshead speed was 1.0 mm/min. Strain was measured by strain gages attached on both sides of each specimen. Stress was calculated from applied load measured by load cell and dimensions of the specimen. 6 specimens were tested for each type of CFRP. Measured Young’s moduli EL , Poisson’s ratios mLT , tensile strengths rmax and failure strains emax are summarized in Table 3. Young’s moduli and Poisson’s ratios were almost identical for both laminates. In contrast, strength and failure strain were significantly lower in the VaRTM laminates. Typically, 0° tensile strength is dominated by the failure of the carbon fiber, and fiber volume fractions were almost the same between prepreg-autoclave CFRP and VaRTM CFRP. Therefore, the lower strength of the VaRTM CFRP is caused by fiber waviness. 4.2. Analysis and discussion In order to calculate the effect of the waviness on the stiffness and strength of the unidirectional CFRP, numerical analysis was conducted. Calculations were based on the method proposed in

Section 2. Numerical simulation was conducted by ABAQUS/Standard 6.14-2. Material properties used were the same as those used in Section 3 (see Table 1). Dimensions of the RUC were determined from the specimen dimensions. According to Fig. 10, the wavelength of one fiber waviness is about 2.2 mm. Averaged thickness of specimens was 1.19 mm. Therefore dimensions of numerical RUC were 2.2 mm  2.2 mm in Y 1  Y 2 directions, and 1.19 mm in Y 3 direction. The model was divided into 0.05 mm  0.05 mm  0.0661 mm brick solid elements (C3D8). In this study, shape of the imperfection is assumed to be sinusoidal, and it is assumed to be uniform in the thickness direction. 0.052 mm wave height was assumed according to measurement in Fig. 10. Then, imperfection ^ i can be written as, u

 L L L ^2 ; u ^3 ¼ ^1 ; u u

   0:052 2pY 1 sin 0; 0; 2 2:2

ð34Þ

Additionally, the perturbation parameter e is set to 1. Stiffness matrices of perfect and imperfect RUCs were calculated using Eqs. (26)–(32). Calculated membrane stiffness matrices (A matrices) of perfect and imperfect models were compared in Fig. 11. Except for the trivial numerical errors (highlighted by green), both stiffness matrices were identical. Although imperfection causes stress concentration, stress is integrated over the RUC volume in the calculation of the stiffness as seen in Eqs. (29) and (30). Therefore high-stress portion and low-stress portion tend to cancel each other, and macroscopic stiffness is insensitive to the fiber waviness. Microscopic stress distribution of the imperfect model was calculated under uniaxial tensile load (e11 ¼ 0:16%). Calculated tensile stress in Y 1 direction is plotted in Fig. 12. Imperfection causes high-stress portion and low-stress portion on the top and bottom surfaces. Tensile stress of high-stress portion reached 3390 MPa, meanwhile macroscopic stress of this RUC was 2448 MPa. Therefore fiber tensile failure may occur at the regions of high-stress because the tensile stress is larger than the 0 tensile strength of prepreg-autoclave CFRP. This result implies that tensile strength is sensitive to the fiber waviness. These calculation results agreed qualitatively to the experimental results. Therefore, we can conclude that the fiber waviness induced in the manufacturing process did not affect the elastic moduli, but reduced the strength of the VaRTM CFRP.

Table 3 Summary of the results of static tensile tests.

Longitudinal Young’s modulus EL (GPa) In-plane Poisson’s ratio mLT 0 tensile strength rmax (MPa) 0 failure strain emax (%)

T800S/#3900-2B (Prepreg-autoclave)

T800SC/XNR6809-XNH6809 (VaRTM)

152.67 0.34 3100 1.80

151.43 0.35 2636 1.60

353

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355

Fig. 12. Tensile stress

r11 distribution of imperfect model. Displacement is magnified by the factor of 5.0.

It should be noted that the ratio of the strength of prepregautoclave CFRP to that of the VaRTM CFRP (1.176) is smaller than the stress concentration factor in the imperfect model simulation (1.385). Two facts cause this difference. Prepreg-autoclave CFRP also contains imperfections; these are smaller than that of VaRTM CFRP, in reality. Moreover, in this simulation, imperfection is assumed to be uniform in the thickness direction. This assumption causes larger stress concentration than real imperfection, the latter is not uniform in the thickness direction, as was assumed.

scopic uniform tensile strain eG is applied to the bar. The representative unit cell (RUC) consists of two elements, whose Young’s moduli are E1 and E2 , respectively. Lengths of the elements are pl and ð1
pÞl, where 0 6 p 6 1. In this example, the bar has an unintended initial geometrical imperfection in the lengths of elements, which disrupts the periodicity. Imperfect displacement  L L L ^1 ; u ^2 ; u ^ 3 ¼ f0; ql; 0g is considered. u Strains of element 1 (eL1 ) and of element 2 (eL2 ) are calculated as,

(

5. Conclusion In this paper, a novel multiscale simulation method for a laminated plate structure, in which geometrical imperfection is considered has been presented. The method is based on the finite element method, perturbation method, and homogenization method. Geometrical imperfection is introduced as input for the calculation, and the calculation is always done using an idealized regular mesh. The method was implemented using the general purpose finite element code ABAQUS by using user subroutines. The proposed method was validated by comparing calculated stress distribution to the results of conventional calculation results. It was demonstrated that the proposed method can calculate accurate stress distribution when the amplitude of the imperfection is small. When the imperfection amplitude is large, differences in outcome are not negligible. The differences are caused by neglecting higher terms of the perturbation method. By considering these higher order terms, the accuracy of the results can be improved. The proposed method was applied to calculate the effect of the fiber waviness on the macroscopic tensile properties of unidirectional VaRTM CFRP. Calculation results revealed that the fiber waviness hardly affects the macroscopic stiffness of the CFRP, because the effect of the waviness is cancelled in the volume average calculation. In contrast, fiber waviness significantly increases the stress concentration, which causes strength reduction in the VaRTM CFRP laminate.

L 1 L 2

e e

)

" ¼


pl1

1 pl

0

1
ð1
p Þl

8 9 8 L9 #> uL > > < 1= < u1 > = L ¼ ½ B  u uL2 2 1 > > > : L> ; ð1
pÞl : L ; u3 u3 0

Therefore, the 0th order principle of virtual work Eq. (24) can be computed as,

( ½BT

2

)

G

AplE1 e

Að1
pÞlE2 eG

AE1 pl

6 AE1 6
4 pl 0


AEpl1 AE1 pl

þ ½BT

AplE1

0

0

Að1
pÞlE2



8 L0 9 8 9 > < u1 > = > <0> = ½B uL0 ¼ 0 2 > : L0 > ; > : > ; 0 u3

38 9 8 9 L0 AE1 eG > > < u1 > < = > = 7 AE2 7 L0 G
ð1
p ¼ u A ð
E þ E Þ e 1 2 Þl 5 2 > > > > : L0 ; : ; AE2 u3
AE2 eG 0

AE2 þ ð1
p Þl

AE2
ð1
p Þl

ð1
pÞl

Appendix A. One dimensional problem To illustrate the proposed method, a closed-form analytical solution to a one-dimensional problem is presented. Fig. A.1 shows a schematic of the example problem. A one-dimensional, bi-material, periodic bar of sectional area A is considered. Macro-

ðA:1Þ

Fig. A.1. 1 dimensional example problem.

ðA:2Þ

354

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355

L0 L0 An additional condition, uL0 has 1 ¼ u3 ¼ 0 is assigned so that u periodicity and that the rigid displacement is equal to zero. Eq. (A.2) is then solved as,

9 8 L0 9 8 0 > > u > > > < 1 = < E
E G > = 2 1 L0 e ¼ E1 þ E2 u2 pl ð1
pÞl > > > > : L0 ; > > : ; u3 0

ðA:3Þ

The 1st order principle of virtual work Eq. (25) is then computed using the 0th order equation results Eq. (A.3). Strain caused by imperfection ^e and coupling strain eC can be calculated as,



^e1



^e2 (

eC1 eC2

8 9 ( q ) > <0> = p ¼ ½B ql ¼ q >
1
p : > ; 0

)

" ¼


pq

8 9 8 L0 9 #> uL0 > < 1 = h i> < u1 > = 0 C L0 ¼ B u uL0 2 2 q 1 >
ð1
pÞ ð1
pÞl > : L0 > ; : L0 > ; u3 u3

q 1 p pl

1 pl

q 1 ð1
pÞ ð1
pÞl

0

ðA:4Þ

ðA:5Þ The First order principle of virtual work Eq. (25) is then calculated as,

" ½BT

AplE1

0

0

Að1
pÞlE2

" ¼ ½B

T

(

AplE1

0

0

Að1
pÞlE2

)

q p q
1
p

"
½BT

C

þ B

" þ ½B

T

AplE1 qp 0

" iT AplE 1 0

8 L0 9 < u1 > = i> BC uL0 2 > Að1
pÞlE2 : L0 > ; u3 #( ) 0 eG

eG

Að1
pÞlE2




1 E1

p E1


E12

þ

1
p E2

C G 2 qAe

ðA:9Þ

Therefore macroscopic Young’s modulus EG ðp; qÞ is,

ðA:6Þ

EG ðp; qÞ ¼

1 p E1

þ


 ð1
pÞ

1 E1

p E1

E2


E12

þ 1
p E2

2 q

ðA:10Þ

Fig. A.2 plots macroscopic Young’s modulus calculated by the proposed method. In this chart, E1 ¼ 1000 (MPa) and E2 ¼ 5000 (MPa) were used. For this problem, the macroscopic Young’s modulus E A can be analytically calculated in closed-form as,

eG

8 L0 9 > < u1 > = ½B uL0 2 > : L0 > ; u3

#

0

1

Þ þ ð1
p E2 E1

8 9 # > uL0 < 1 > = 0 ½B uL0 2 q >
Að1
pÞlE2 ð1
pÞ : L0 > ; u3 #( ) 0 eG Að1
pÞlE2

ðA:8Þ

1

r1 ¼ r2 ¼ B @p

E A ðp; qÞ ¼ pþq E1

where 1st and 2nd terms on the right-hand-side mean direct effect of the imperfection and the effect of the coupling strain, respectively. 3rd and 4th terms correspond to the effect of the volume change caused by imperfection. 5th and 6th terms show the coupling effect between virtual displacement and imperfection. Eq. (A.6) is simplified as,

2

8 9 0  > 8 L1 9 > > >  > > > > E E > > q
21 þ 2 2 ðE1
E2 Þ < u1 > = > < = p l ð1
pÞ l G L1 e ql þ   2 ¼ u2 E1 E > > þ 2 : L1 > ; > > > pl ð1
pÞl > > > > u3 > > : ; 0

By using 0th and 1st order solutions, stresses in elements 1 and 2 are calculated from Eq. (15) as,

# h

0

0

L1 Similar to the 0th equation, additional condition uL1 1 ¼ u3 ¼ 0 is applied. By substituting Eq. (A.3), Eq. (A.7) is then solved as,

0

q
Að1
pÞlE2 ð1
pÞ

" iT AplE 1 C þ B 0 h

AplE1

AplE1 qp 0

"
½BT h

8 L1 9 > < u1 > = ½B uL1 2 > : L1 > ; u3 #

#

Fig. A.2. Macroscopic Young’s modulus of 1 dimensional example problem calculated by the proposed method.

3

AE1 8 L1 9
AEpl1 0 pl < u1 > = 7> 6 6
AE1 AE1 þ AE2
AE2 7 uL1 7 6 pl pl ð1
pÞl ð1
pÞl 2 5> 4 : L1 > ; AE2 AE2 u3 0
ð1
p Þl ð1
pÞl 8 9   L0 > >
AE1 qp
AE1 plq2 uL0 > > 2
u1 > > > > < = n o  q q q q L0 L0 ¼ AE1 p þ AE2 ð1
pÞ þ
AE1 p2 l þ AE2 ð1
pÞ2 l u1
u2 > > > >  L0  > > > > q q L0 : ;
AE2 ð1
pÞ
AE2 ð1
p u
u 1 2 Þ2 l

1 Þ þ ð1
p
q E2

ðA:11Þ

It is immediately seen that the macroscopic Young’s modulus calculated by the proposed method (A.10) is equal to the 1st order Taylor approximation of the analytical solution around the point q ¼ 0.

EG ðp; qÞ ¼ E A ðp; 0Þ þ

@E A ðp; 0Þq @q

ðA:12Þ

The above discussion demonstrates that an approximate solution can be acquired by using the proposed method, and it corresponds to the 1st order Taylor approximation of the exact solution around the solution for the perfect structure. Therefore, the proposed method provides a good approximation when the imperfection amplitude is small. References

ðA:7Þ

[1] Babuška I. Homogenization approach in engineering. In: Glowinski R, Lions JL, editors. Computing Methods in Applied Sciences and Engineering. Lecture Notes in Economics and Mathematical Systems, vol. 134. Berlin Heidelberg: Springer; 1976. p. 137–53. http://dx.doi.org/10.1007/978-3-64285972-4_8.

A. Yoshimura et al. / Composite Structures 153 (2016) 345–355 [2] Bensoussan A, Lions J-L, Papanicolaou G. Asymptotic Analysis for Periodic Structures. AMS Chelsea Publishing; 1978. [3] Sanchez-Palencia E. Non-Homogeneous Media and Vibration Theory. Lecture Notes in Physics, vol. 127. Berlin Heidelberg: Springer-Verlag; 1980. http://dx. doi.org/10.1007/3-540-10000-8. [4] Bakhvalov NS, Panasenko G. Homogenisation: Averaging Processes in Periodic Media. Mathematics and its Applications, volume 36. Netherlands: Springer; 1989. http://dx.doi.org/10.1007/978-94-009-2247-1. [5] Yokozeki T, Kobayashi Y, Aoki T, Yoshida D, Hirata T. VaRTM process of composites using porous mold. Adv Compos Mater 2013;22:99–107. [6] Matsuzaki R, Seto D, Todoroki A, Mizutani Y. In situ void content measurements during resin transfer molding. Adv Compos Mater 2013;22:239–54. [7] Yokozeki T, Yamazaki W, Kobayashi Y. Investigation into property control of VaRTM composites by resin infusion process. Adv Compos Mater 2015;24:495–507. [8] Yokozeki T, Iwamoto K. Effects of core machining configuration on the debonding toughness of foam core sandwich panels. Adv Compos Mater 2016;25:45–58. http://dx.doi.org/10.1080/09243046.2014.958302. [9] Centea T, Grunenfelder L, Nutt S. A review of out-of-autoclave prepregs – material properties, process phenomena, and manufacturing considerations. Composites Part A 2015;70:132–54. [10] Sugimoto S, Aoki Y, Hirano Y, Nagao Y. A study of quality assurance of VaRTM composite wing structure. In: 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. American Institute of Aeronautics and Astronautics; 2007. http://dx.doi.org/10.2514/6.2007-2338. [11] Hirano Y, Mizutani T, Iwahori Y, Nagao Y. An investigation on spring-in behavior of Va-RTM composite wing structure. In: Proceedings of 16th International Conference on Composite Materials, Kyoto, Japan. [12] Hirano Y, Aoki Y, Iwahori Y, Sugimoto S, Nagao Y. Development and mechanical properties verification of VaRTM composite full-scale wing demonstrator. In: Proceedings of SAMPE 2008. Long Beach, CA, USA: Society for the Advancement of Material and Process Engineering; 2008. [13] Lane SA, Higgins J, Biskner A, Sanford G, Springer C, Berg J. Out-of-autoclave composite fairing design, fabrication, and test. ASME J Manuf Sci Eng 2011;133:031020-1–031020-11. [14] Aoki Y, Hirano Y, Sugimoto S, Iwahori Y, Nagao Y, Ohnuki T. Durability and damage tolerance evaluation of VaRTM composite wing structure. In: ICAF 2011 Structural Integrity: Influence of Efficiency and Green Imperatives. p. 561–72. http://dx.doi.org/10.1007/978-94-007-1664-3_45. [15] Budiansky B, Fleck N. Compressive failure of fibre composites. J Mech Phys Solids 1993;41:183–211. [16] Kyriakides S, Arseculeratne R, Perry E, Liechti K. On the compressive failure of fiber reinforced composites. Int J Solids Struct 1995;32:689–738. [17] Schultheisz CR, Waas AM. Compressive failure of composites, part I: testing and micromechanical theories. Prog Aerosp Sci 1996;32:1–42. [18] Waas AM, Schultheisz CR. Compressive failure of composites, part II: experimental studies. Prog Aerosp Sci 1996;32:43–78. [19] Lee S, Waas A. Compressive response and failure of fiber reinforced unidirectional composites. Int J Fract 1999;100:275–306. [20] Prabhakar P, Waas AM. Interaction between kinking and splitting in the compressive failure of unidirectional fiber reinforced laminated composites. Compos Struct 2013;98:85–92. [21] Prabhakar P, Waas AM. Micromechanical modeling to determine the compressive strength and failure mode interaction of multidirectional laminates. Composites Part A 2013;50:11–21. [22] Cox BN, Dadkhah MS, Morris W. On the tensile failure of 3D woven composites. Composites Part A 1996;27:447–58. [23] Mouritz A, Cox B. A mechanistic approach to the properties of stitched laminates. Composites Part A 2000;31:1–27.

355

[24] Agius S, Fox B. Rapidly cured out-of-autoclave laminates: understanding and controlling the effect of voids on laminate fracture toughness. Composites Part A 2015;73:186–94. [25] Schell J, Renggli M, van Lenthe G, Müller R, Ermanni P. Micro-computed tomography determination of glass fibre reinforced polymer meso-structure. Compos Sci Technol 2006;66:2016–22. [26] Djukic LP, Herszberg I, Walsh WR, Schoeppner GA, Prusty BG. Contrast enhancement in visualisation of woven composite architecture using a microCT scanner. Part 2: tow and preform coatings. Composites Part A 2009;40:1870–9. [27] Djukic LP, Herszberg I, Walsh WR, Schoeppner GA, Prusty BG, Kelly DW. Contrast enhancement in visualisation of woven composite tow architecture using a microCT scanner. Part 1: fabric coating and resin additives. Composites Part A 2009;40:553–65. [28] Djukic LP, Pearce GM, Herszberg I, Bannister MK, Mollenhauer DH. Contrast enhancement of microCT scans to aid 3D modelling of carbon fibre fabric composites. Appl Compos Mater 2013;20:1215–30. [29] Iida A, Okamoto K, Ito S, Imai Y, Shimamoto D, Hotta Y. Observation in widely various the CFRP containing metal-plated carbon fiber as a tracer. In: Proceedings of JCCM-5, Kyoto, Japan (in-Japanese). [30] Centea T, Hubert P. Measuring the impregnation of an out-of-autoclave prepreg by micro-CT. Compos Sci Technol 2011;71:593–9. [31] Scott A, Sinclair I, Spearing S, Mavrogordato M, Hepples W. Influence of voids on damage mechanisms in carbon/epoxy composites determined via high resolution computed tomography. Compos Sci Technol 2014;90:147–53. [32] Yoshimura A, Hosoya R, Koyanagi J, Ogasawara T. X-ray computed tomography used to measure fiber orientation in CFRP laminates. Adv Compos Mater 2015. in press. [33] Naouar N, Vidal-Sallé E, Schneider J, Maire E, Boisse P. Meso-scale FE analyses of textile composite reinforcement deformation based on X-ray computed tomography. Compos Struct 2014;116:165–76. [34] Czabaj MW, Riccio ML, Whitacre WW. Numerical reconstruction of graphite/ epoxy composite microstructure based on sub-micron resolution X-ray computed tomography. Compos Sci Technol 2014;105:174–82. [35] Takano N, Zako M, Kubo F, Kimura K. Microstructure-based stress analysis and evaluation for porous ceramics by homogenization method with digital imagebased modeling. Int J Solids Struct 2003;40:1225–42. [36] Matsuda T, Okumura D, Ohno N, Kawai M. Three-dimensional microscopic interlaminar analysis of cross-ply laminates based on a homogenization theory. Int J Solids Struct 2007;44:8274–84. [37] Yoshimura A, Okabe T. Damage growth analysis in particle-reinforced composite using cohesive element. Adv Compos Mater 2011;20:569–83. [38] Lee C-Y, Yu W. Homogenization and dimensional reduction of composite plates with in-plane heterogeneity. Int J Solids Struct 2011;48:1474–84. [39] Terada K, Kato J, Hirayama N, Inugai T, Yamamoto K. A method of two-scale analysis with micro-macro decoupling scheme: application to hyperelastic composite materials. Comput Mech 2013;52:1199–219. [40] Terada K, Hirayama N, Yamamoto K, Kato J, Kyoya T, Matsubara S, Arakawa Y, Ueno Y, Miyanaga N. Applicability of micro-macro decoupling scheme to twoscale analysis of fiber-reinforced plastics. Adv Compos Mater 2014;23:421–50. [41] Terada K, Hirayama N, Yamamoto K, Matsubara S. Numerical plate testing for linear multi-scale analysis of composite plates. Trans Jpn Soc Comput Eng Sci 2015;2015 (in-Japanese). [42] Japan Aeroscape Exploration Agency (JAXA). In: JAXA Advanced Composites Database System: JAXA-ACDB. URL: . [43] Song S, Waas AM, Shahwan KW, Faruque O, Xiao X. Compression response of 2d braided textile composites: single cell and multiple cell micromechanics based strength predictions. J Compos Mater 2008;42:2461–82. [44] SACMA SRM-4R. Tensile Properties of Oriented Fiber-Resin Composites. Suppliers of Advanced Composite Materials Association (SACMA); 1994.