Composites: Part B 60 (2014) 378–391 Contents lists available at ScienceDirect Composites: Part B journal homepage: www.elsevier.com/locate/compositesb Thermomechanical properties and stress analysis of 3-D textile composites by asymptotic expansion homogenization method Muhammad Ridlo Erdata Nasution a,⇑, Naoyuki Watanabe a, Atsushi Kondo a,b, Arief Yudhanto a a b Department of Aerospace Engineering, Tokyo Metropolitan University, 6-6 Asahigaoka, Hino-shi, Tokyo 191-0065, Japan e-Xtream Engineering, MSC Software Company, 1-23-7 Nishishinjuku, Shinjuku-ku, Tokyo 160-0023, Japan a r t i c l e i n f o Article history: Received 11 September 2013 Received in revised form 2 November 2013 Accepted 22 December 2013 Available online 3 January 2014 Keywords: A. Fabrics/textile B. Thermomechanical C. Computational modelling a b s t r a c t This paper formulates asymptotic expansion homogenization method for the analysis of 3-D composites whereby thermomechanical effects are considered. The equivalent thermomechanical properties, viz. elastic constants and coefficient of thermal expansion, of 3-D orthogonal interlock composites are obtained based on a unit-cell derived from the optical microscopy observation. The Young’s modulus and Poisson’s ratio obtained from the analysis are compared with those obtained by experiments. The results show a good agreement especially for the Young’s modulus. Localization analysis based on the present formulation is also presented to assess the stresses within the constituents of 3-D composites (fiber tows and resin region) due to mechanical loading and temperature difference. In order to facilitate the localization analysis, the unit-cell of 3-D composites is first subjected to temperature difference in order to obtain thermal residual stresses. Subsequently, a beam model utilizing homogenized thermomechanical properties of 3-D composites is built to obtain stresses due to mechanical loading. Two loading cases are considered: beam under uniaxial tension (Case 1) and beam under bending load (Case 2). The thermal residual stresses and stresses due to mechanical loading in each case are then superposed to obtain total stresses. The total stresses are then used to understand the complete response of composite constituents under mechanical and thermal loadings. Ó 2013 Elsevier Ltd. All rights reserved. 1. Introduction Composite materials have been widely used in aerospace structures due to their excellent properties as compared to the metallic materials, e.g. high strength-to-weight ratio, good corrosion resistance. The recent development of composite materials aims to obtain better characteristics, such as enhanced damage tolerance and delamination resistance. However, the complexity and heterogeneity of composites often come with costly computational efforts. A modeling approach whereby all composite constituents (fibers and matrix) are explicitly built renders cumbersome meshing techniques, and requires high computational resources. An excellent technique to cope with this problem is to idealize the composite as a unit-cell representing the whole constituents and the corresponding responses, while the equivalent properties and other detailed stresses can be obtained. Representative unit-cell approach has been commonly used in the analysis of 3-D textile composites. Several studies employing this approach in order to predict the failure strength and ballistic impact damage behavior ⇑ Corresponding author. Tel./fax: +81 42 585 8619. E-mail address: nasution-ridlo@ed.tmu.ac.jp (M.R.E. Nasution). 1359-8368/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compositesb.2013.12.038 of 3-D orthogonal woven composites were performed by Tan et al. [1] and Sun et al. [2], respectively. Dealing with the so-called unit-cell, asymptotic expansion homogenization (AEH) method is often employed. This method assumes that the unit-cell can be repeated in specific directions to represent periodicity. AEH regards the composite structure as a homogeneous macrostructure, which is formed by a large number of periodic and heterogeneous microstructure. AEH incorporates asymptotic expansion analysis to couple between microscopic and macroscopic behavior of a system [3]. Generally, AEH is often formulated to obtain thermomechanical properties of heterogeneous structures. Guedes and Kikuchi [4] derived a rigorous formulation of AEH method, and applied the formulation to analyze sandwich honeycomb plate. AEH is also used to analyze textile composites [5–8] and metal-matrix composites [9]. Francfort [10] developed asymptotic expansion technique for the case of linear thermoelasticity. The incorporation of thermal effects into the heterogeneous unit-cell was performed by several authors to obtain mechanical properties as well as the coefficient of thermal expansions (CTE) [11–13]. In order to investigate the response of a unit-cell, localization analysis is usually performed in the framework of AEH. The localization analysis can be seen as a reverse scheme of the 379 M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 homogenization method that aims to obtain the local stresses of a unit-cell in a microscopic scale due to the external loads applied onto the homogeneous macrostructure [5,14]. However, in the studies described in Refs. [5,14], thermal residual stresses were not considered. In fact, thermal residual stresses may affect the damage behavior of composites [15,16]. In this paper, a detailed mathematical and finite element formulation of AEH method is proposed to study the thermomechanical response of 3-D textile composite. The textile composite under consideration is 3-D orthogonal interlocked composites. The numerical analysis includes the calculation of equivalent thermomechanical properties, and the localization analysis whereby the critical regions in the constituents of 3-D orthogonal interlock composites are identified based on the maximum stresses. Two load cases, namely uniaxial tension and bending, are studied to determine the critical loading condition by which the damage may initiate in 3-D orthogonal interlock composites. The method developed in this paper may potentially aid the material selection processes, damage analysis of composites, design of textile composites and other engineering steps in various industries (aerospace, automotive, marine, electronic devices [17], medical [18], etc.). Apparently, the composite material selected for present analysis, i.e. 3-D orthogonal interlock composites, has a limited commercial application due to high cost pertaining to the manufacturing processes. This specific type of 3-D composites has been applied only for a few specific areas where 2-D laminates and metallic materials cannot satisfy the desired properties [19]. Nevertheless, the type of 3-D composites selected herein is deemed appropriate to test the capability of homogenization method because of their microstructural complexities. and u0 denote actual displacement and macroscopic displacement. Similar to actual stress, actual displacement is also periodic. The macroscopic displacement is changing linearly from A to A0 . 2.2. Periodic vector function From the viewpoint of microscopic and macroscopic scales, microstructural variables within unit-cell vary. Accordingly, periodic vector function g which is a function of macroscopic coordinate x and microscopic coordinate y is introduced g e ðxÞ ¼ gðx; yÞ; e gðx; yÞ ¼ gðx; y þ YÞ ð1Þ ð2Þ where Y is the dimension of the unit-cell. Derivatives of g with respect to the macroscopic coordinate x can be obtained according to the chain rule: @ @g 1 @g ½gðx; yÞ ¼ þ @xi @xi e @yi ð3Þ The following expressions are used to describe the limit of integration of Y-periodicity vector function as e approaches zero from the positive side e!0 Z 2. Asymptotic expansion homogenization method The homogenization method is generally employed in a multiscale analysis, which involves at least two spatial scales, namely microscopic and macroscopic scales. In this regard, homogenization method correlates both scales by considering a periodic and heterogeneous microstructure in the microscopic scale as a homogeneous macrostructure in macroscopic scale. Viewed from macroscopic scale, the periodicity of the microstructure is assumed very small. Two aforementioned spatial scales can be seen in Fig. 1, whereby Fig. 1(a) is a unit-cell in microscopic scale, whilst Figs. 1(b) and (c) are heterogeneous and homogenous macrostructure in macroscopic scale. In Fig. 1(a), it is noted that the unit-cell is heterogeneous and periodic. The coordinate system used in microscopic scale is y = (y1, y2, y3), while that used in macroscopic scale is x = (x1, x2, x3). Both scales are correlated by a very small positive number e = x/y. When e approaches zero, the heterogeneous macrostructure (Fig. 1(b)) can be regarded as a homogeneous macrostructure (Fig. 1(c)). The periodicity of a very small microstructure in the macroscopic scale necessitates the use of asymptotic expansion analysis. Asymptotic expansion analysis has a powerful capability to transfer the properties and behavior of the microstructure into the macroscopic representation. Consider a beam subjected to uniaxial traction force shown in Fig. 2(a). Four unit-cells (UC), namely UC1, UC2, UC3 and UC4, can be identified as a subset of the beam. Each unit-cell is assumed to have two black bricks and two white bricks. The distance between UC1 and UC4 is represented by A–A0 . Figs. 2(b) and (c) show the responses of the unit-cell in terms of stress (r11) and displacement (u1). In Fig. 2(b), r0 and rH denote actual stress and homogenized stress, respectively. The actual stress is periodic, while homogenized stress is constant. In Fig. 2(c), ue x The vector function considers the periodicity on the microscopic coordinate y, which is called Y-periodic. Homogenization method is implemented by assuming that the unit-cell is periodic, viz. infinitely repeated unit-cell in three directions (x1, x2 and x3). Due to Y-periodicity, Eq. (1) can be written as follows limþ 2.1. General concept where y ¼ limþ e e!0 Xe Z gðx; yÞdX ! Se 1 jYj Z Z 1 jYj Z Z gðx; yÞdX ! X gðx; yÞdYdX ð4Þ gðx; yÞdSdX ð5Þ U S X where |Y| denotes unit-cell volume. 2.3. Principle of virtual works Structural problems in equilibrium condition can be effectively solved by using the principle of virtual work. The principle specifies that the equilibrium of a system is satisfied if and only if the integration of the product of the forces and virtual displacement on every point of the system is equal to zero [20]. Eq. (6) mathematically expresses the principle of virtual work by involving the equilibrium equation of force.  Z  @ rij þ fi v i dX ¼ 0 @xj X ð6Þ where rij denotes the stresses acting on the domain X, fi is body force, vi is virtual displacement, and xj is macroscopic (global) coordinate system. The application of integration by parts and divergence theorem to Eq. (6) yields Z X rij @v i dX ¼ @xj Z fi v i dX þ X Z t i v i dC ð7Þ Ct where ti is traction force, which is the product of stresses and normal vector on the boundary Ct. In this study, the constitutive equation is linear, and the deformation is small. Stress–strain relationship in linear thermomechanical problem can be expressed as follows t rij ¼ C ijkl ðetot kl  ekl Þ ð8Þ t where Cijkl is elastic constants tensor, etot kl is total strain, and ekl denotes thermal strain. The mechanical strain–displacement 380 M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 Fig. 1. (a) Heterogeneous unit-cell, (b) heterogeneous macrostructure and (c) homogeneous macrostructure. Fig. 2. (a) Beam subjected to traction force, (b) stress response of the unit-cell and (c) displacement response of the unit-cell. relationship and thermal strain are described by Eqs. (9) and (10), respectively. em kl ¼   1 @uk @ul þ 2 @xl @xk etkl ¼ akl DT ð9Þ ð10Þ where u is the actual displacement, akl is the coefficient of thermal expansion and DT is the temperature difference. By substituting Eqs. (9) and (10) into Eq. (8), and considering the symmetry of elastic tensor, stresses can be re-expressed as follows rij ¼ C ijkl   @uk  akl DT @xl ð11Þ Consider again the heterogeneous unit-cell and macrostructure as shown in Figs. 1(a) and (b), respectively. The macrostructure consists of a large amount of unit-cells, where the unit-cell is composed of ‘‘solid part’’ ¥ and ‘‘hollow part’’ h. The unit-cell is subjected to traction p on surface S in microscopic field Y. The heterogeneous microstructure in Fig. 1(b) is subjected to body forces f and tractions t. Displacement is prescribed on the boundary Cd. The principle of virtual works is then employed to solve the thermomechanical problem by substituting Eq. (11) into Eq. (7) and 381 M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 by incorporating Se pei v i dS into the equation, which yields a governing equation as follows R Z Xe C eijkl  e  Z Z Z @uk @v i  aekl DT dX ¼ fie v i dX þ t i v i dC þ pei v i dS @xl @xj Xe Ct Se ð12Þ Superscript e in Eq. (12) indicates that the variable is a function of total region (e.g. including microstructure). 2.4. Formulation of AEH method le The actual displacement asymptotic expansion series: k in Eq. (12) can be represented by uek ¼ u0k þ eu1k þ e2 u2k þ      @ðu0k þ eu1k þ e2 u2k þ   Þ @v i C eijkl  aekl DT dX @xl @xj Xe Z Z Z ¼ fie v i dX þ t i v i dC þ pei v i dS Z Se Ct e2 : Order of 1 ¼ Z e : 1 Z e2 1 e Xe Z C eijkl e Xe C ijkl @u0k @ v i dX ¼ 0 @yl @yj " ð14Þ #  0  @u0k @ v i @uk @u1k @v i e dX þ þ  akl DT @yl @xj @xl @yl @yj ð16Þ " #   1  @u0k @u1k @uk @u2k @ v i @v i e dX Order of e : C ijkl þ  akl DT þ þ @xl @yl @xj @xl @yl @yj Xe Z Z ¼ fie v i dX þ t i v i dC ð17Þ Z 0 Xe " #  0  @u0k @ v i @uk @u1k @v i C ijkl þ þ  akl DT dYdX @yl @xj @xl @yl @yj X U Z Z 1 ¼ p v i dSdX jYj X S i 1 jYj Z Z By choosing v = v(x) and by considering the expression (20), following equation is obtained pi v i ðxÞdS ¼ 0 ð22Þ e Ct By choosing v = v(y) and considering Eq. (22), following equation is obtained 1 jYj " #  @u0k ðxÞ @u1k @ v i ðyÞ dYdX ¼ 0 C ijkl þ  akl DT @yj @yl @xl U Z Z X u1k ðx; yÞ ¼ vpq k ðyÞ @u0p ðxÞ  wk ðyÞ @xq ð24Þ where v and w denote the characteristic displacement vectors for elastic and thermal problem, respectively, which are used to separately solve microscopic problem and macroscopic problem. By substituting Eq. (24) into Eq. (23), following equation is obtained Z Z @ vpq @u0p ðxÞ @u0k ðxÞ @ v i ðyÞ 1 dYdX  C ijkl k @yj jYj X U @yl @xq @xl X U   Z Z @ v i ðyÞ 1 @wk @ v i ðyÞ  dYdX  C ijkl þ akl DT dYdX @yj jYj X U @yj @yl 1 jYj Z Z C ijkl ¼0 ð25Þ 2.5.1. Order of e2 By multiplying Eq. (15) with e2, and by taking the limit as e ? 0+, (i.e. using expression (4)), following equation is obtained @u0 @ v i C ijkl k dYdX ¼ 0 @yl @yj U Z Z X ð23Þ where u0k represents macroscopic displacement, which is a function of macroscopic coordinate x, while u1k represents microscopic displacement obtained by involving the characteristic displacement vectors. Microscopic displacement u1k can be written as follows 2.5. Solving the hierarchical equations 1 jYj ð21Þ S ð15Þ pei v i dS Se ð20Þ 2.5.2. Order of e1 By multiplying Eq. (16) with e1, and by taking the limit as e ? 0+ (i.e. using expression (4) and (5)), following equation is obtained Z It is important to note that the thermal strain has been incorporated into Eq. (14). By taking the derivatives of actual and virtual displacement, i.e. by using Eq. (3), and by assuming that the limit of all integral equations in Eq. (14) exist when e ? 0+, Eq. (14) can be satisfied if the terms of each order of e is zero. Therefore, Eq. (14) is arranged into three hierarchical equations based on the order of e as follows Order of u0k ¼ u0k ðxÞ ð13Þ By substituting Eq. (13) into Eq. (12), an expanded form of the governing equation can be obtained Xe The third and fourth terms of square brackets in Eq. (19) cancel each other due to the Y-periodicity, i.e. in each pair of corresponding surfaces Ya and Yb. Based on mathematical treatment in Ref. [4], the remaining equation implies that ð18Þ Eq. (25) can be expressed in a compact form as follows 1 jYj Z Z " X U ! !# @ vkl @wp @u0k ðxÞ p C ijkl  C ijpq  C ijpq þ apq DT @yq @yq @xl @ v i ðyÞ dYdX ¼ 0 @yj ð26Þ Virtual displacement v is arbitrary, and it can be a function of macroscopic coordinate x (v = v(x)) or microscopic coordinate y (v = v(y)). If v = v(x), Eq. (18) can be automatically satisfied. If v = v(y), Eq. (18) can be expanded by applying integration by parts and Gauss’ divergence theorem as follows Eq. (26) can be decoupled into two equations representing the elastic and thermal problems to separately calculate the solution of two characteristic displacement vectors (v and w). Two microscopic equilibrium equations can be respectively written as follows 3   R @u0k @u0k @ Z  C v ðyÞdY þ C n v ðyÞdS i j i ijkl ijkl U @yj S 1 @yl @yl 7 6 5dX ¼ 0 4 R R @u0k @u0k jYj X þ Y a C ijkl @y nj v i ðyÞdY a þ Y b C ijkl @y nj v i ðyÞdY b Elastic problem : 2 R l Z C ijkl  C ijpq U l ð19Þ Thermal problem : Z U C ijpq ! @ vkl @ v i ðyÞ p dY ¼ 0 @yq @yj @wp þ apq DT @yq ! @ v i ðyÞ dY ¼ 0 @yj ð27Þ ð28Þ 382 M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 2.5.3. Order of e0 By taking the limit of Eq. (17) as e ? 0+ (i.e. using expression (4)), following equation is obtained Z f2gT fCgkl dY n  1 jYj Z _ f2gT ½CfwgdY n þ Z X " #   1  @u0k @u1k @uk @u2k @ v i @v i  C ijkl þ  akl DT þ þ dYdX @xl @yl @xj @xl @yl @yj U Z Z Z 1 ¼ fi v i dYdX þ t i v i dC jYj X U Ct ð29Þ By selecting v = v(x), macroscopic equilibrium equations can be obtained.  0   Z @uk ðxÞ @u1k 1 @ v i ðxÞ C ijkl þ  akl DT dY dX @xj @yl @xl X jYj U  Z  Z Z 1 ¼ fi dY v i ðxÞdX þ t i v i ðxÞdC X jYj U Ct Z  ð30Þ Substitution of Eq. (24) into Eq. (30) yields  @u0k ðxÞ @ v i ðxÞ dY dX @xj @xl X U ! Z @ vpq @u0p ðxÞ 1 @ v i ðxÞ  C ijkl k dY dX @yl @xq @xj X jYj U    Z  Z 1 @wk @ v i ðxÞ  C ijkl þ akl DT dY dX jYj @xj @y X U l  Z Z  Z 1 ¼ fi dY v i ðxÞdX þ t i v i ðxÞdC X jYj U Ct 1 jYj Z Z  Z C ijkl X C 0ijkl ðxÞ @xl @ v i ðxÞ dX ¼ @xj ð31Þ sij ðxÞ X sij ðxÞ ¼ rtij ðxÞ ¼ bi ðxÞ ¼ Z C ijkl  C ijpq U 1 jYj Z C ijpq 1 jYj Z C ijpq apq DTdY 1 jYj U ! @ vkl p dY @yq @wp dY @yq ð32Þ fi dY where kl fv_ g ¼ ½Bfvgkl ð39Þ _ ¼ ½Bfwg fwg ð40Þ f2g ¼ ½Bfv g ð41Þ kl where {v} denotes the elastic corrector of mode ‘kl’, {w} is the thermal corrector, {v} is the virtual displacement vector, [C] is the elastic constants matrix, {C}kl is the column ‘kl’ of matrix [C], and [B] is the strain shape function matrix obtained from the derivatives of shape function matrix [N]. ¥n indicates the n-th element of domain ¥. Eqs. (37) and (38) can be re-arranged by removing the virtual displacement vectors Z ½BT ½C½BdY n fvgkl ¼ Z Z ½BT ½C½BdY n fwg ¼  ½BT fCgkl dY n ð42Þ Un Z ½BT ½CfaDTgdY n ð43Þ Un It is important to note that the elastic correctors have a symmetric property, e.g. vkl = vlk. In other words, six independent modes of v11, v22, v33, v12, v23 and v13 have to be calculated independently. vkl ð0; y2 ; y3 Þ ¼ vkl ðl1 ; y2 ; y3 Þ vkl ðy1 ; 0; y3 Þ ¼ vkl ðy1 ; l2 ; y3 Þ vkl ðy1 ; y2 ; 0Þ ¼ vkl ðy1 ; y2 ; l3 Þ ð44Þ ð33Þ where l1, l2 and l3 are the dimension of unit-cell in each axis. To ensure that the unique solution of periodic boundary condition is attained on each node of the opposite faces in the finite element model, following equation is applied: ð34Þ fvgkl ¼ fvgkl  1 jYj Z fvgkl dY ð45Þ U ð35Þ The periodic boundary condition is represented by elastic corrector v. Eqs. (44) and (45) are also valid for thermal corrector w by replacing {v}kl by {w}. ð36Þ 3.3. Homogenized thermomechanical properties U Z ð38Þ Periodic boundary conditions are applied on the surfaces of the unit-cell. Periodic boundary conditions in in-plane and out-ofplane directions of the unit-cell can be written as follows where 1 jYj f2gT ½CfaDTgdY n ¼ 0 3.2. Periodic boundary conditions Z Ct C 0ijkl ðxÞ ¼ ð37Þ Un Un @ v i ðxÞ dX þ rtij ðxÞ @xj X Z @ v i ðxÞ  dX þ bi ðxÞv i ðxÞdX @xj X Z þ ti ðxÞv i ðxÞdC Z Z Un The compact form of Eq. (31) can be written as follows @u0k ðxÞ kl f2gT ½Cfv_ g dY n ¼ 0 Un Un Z Z Un Z U 3. Finite element formulation 3.1. Characteristic displacement vectors Characteristic displacement vectors of v and w (also known as correctors) already given in Eqs. (27) and (28) can be obtained by finite element method. Eqs. (27) and (28) re-expressed in the matrices and vectors reflecting their element parameters as follows Homogenized elastic constants calculated using Eq. (33) are described in the ‘kl’ column of the matrix [C] by using following equation: fC 0 gkl ¼ X 1 Z fCgkl  ½C½Bfvgkl dY jYj U n n ð46Þ It is necessary to note that due to symmetry of stresses and strain energy densities [21], [C] has the following properties: C ijkl ¼ C jikl ¼ C klij ¼ C ijlk ð47Þ M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 Homogenized coefficient of thermal expansion is calculated by following equation: a0kl ¼ " ! # Z S0klij 1 @wp C ijpq þ apq DT dY DT jYj Un @yq ð48Þ Eq. (48) can be re-expressed in the matrix and vector forms as follows fa0 g ¼ ½S0  X 1 DT n jYj Z 383 Substitution of Eqs. (33) and (48) into Eq. (56) yields the same result with the right hand side of Eq. (57). The relationship between rHij and r0ij is thus expressed as follows rHij ¼ 1 jYj Z r0ij dY ð58Þ U 4. Numerical models, assumptions and analysis sequence ½Cð½Bfwg þ fagDTÞdY ð49Þ 4.1. Numerical models Un 0 where [S ] is homogenized compliance tensors, i.e. inverse of homogenized elastic constants. 3.4. Macroscopic displacement Macroscopic displacement is obtained by solving Eq. (32) as follows Z ½BT ½C 0 ½BdXn fu0 g ¼ Xn Z ½BT fsgdXn þ Xn þ Z ½BT frt gdXn Z ½N T fbgdXn þ ½NT ftgdC Xn Z Xn ð50Þ Ct 3.5. Calculation of stresses Stresses on each point of domain can be obtained by following equation reij ¼ C ijkl  e  @uk  aekl DT @xl ð51Þ By substituting the actual displacement uek stated in Eq. (13) and taking its derivatives by using Eq. (3), following equation is obtained reij ¼ C ijkl  1   0   @uk @u1k @uk @u2k þ e2 ð  Þ þ  aekl DT þ e þ @xl @yl @xl @yl ð52Þ Eq. (52) can be re-expressed in a compact form as follows e rij ¼ r0ij þ er1ij þ e2 ð  Þ ð53Þ where r0ij is the first approximation of stresses. Substitution of Eq. (24) into Eq. (52) yields reij  r0ij ðx; yÞ ¼ ! @ vkl @wp @u0k p C ijkl  C ijpq  C ijpq  C ijpq apq DT @yq @xl @yq ð54Þ In finite element method viewpoint, the first approximation of stress is fr0 g ¼ ½½Ckl  ½C½B½vkl ½Bfu0 g  ½C½Bfwg  ½CfagDT ð55Þ Subscript (and superscript) ‘kl’ indicates that the matrices are arranged based on column ‘kl’. Homogenized stresses rHij of the unit-cell is given as follows rHij ¼ C 0ijkl  0  @uk  a0kl DT @xl Relationship between rHij and averaging operators as follows 1 jYj Z U r 0 ij dY ð56Þ r0ij can be established by using ! Z @ vkl @u0 1 p C ijkl  C ijpq dY k  C ijpq @yq @xl jYj U U Z @wp 1  dY  C ijpq apq DTdY jYj U @yq 1 ¼ jYj Z Fig. 3(a) shows the schematic of 3-D orthogonal interlock composites, and the side-view observed from optical microscopy is shown in Fig. 3(b). A unit-cell of 3-D orthogonal interlocked composites was identified, and the idealized geometry is depicted in Fig. 3(c). The unit-cell consists of x-tow, y-tow, z-tow and resin region. In Fig. 3, the subscripts 1, 2 and 3 of x and y coordinate systems denote the fiber direction of x-, y- and z-tows, respectively. This designation also corresponds to the coordinate of macroscopic beam model in Figs. 4 and 5. The fiber tows are assumed to consist of fiber and matrix phases arranged in hexagonal pattern. For x- and y- tows, the fiber is T300, while the resin is epoxy EP828. The fiber in z-tow is T-900, whilst the resin is also EP828. Resin region contains EP828. The thermomechanical properties of T-300 [1,22], T-900 [22,23] (assumed to have similar properties with T-800) and EP828 [24,25] are given in Table 1. Thermomechanical properties of fiber tows are obtained from the homogenization analysis of the hexagonal model. For hexagonal model, the unit-cell consists of 2664 elements and 8773 nodes, whilst that of 3-D orthogonal interlock composite consists of 12,800 elements and 58,097 nodes. 4.2. Assumptions As mentioned in Section 2.3, the deformation assumed for the numerical models is considerably small, and geometric nonlinearity is not included. Nevertheless, present method can potentially be extended to facilitate nonlinear, large deformation case. One technique can be employed to further improve present method, e.g. modification of constitutive model involving the second Piola–Kirchhoff stress tensor and Cauchy-Green strain tensor. Such technique was implemented in the analysis of plain woven fabrics [26]. Another technique was carried out in the analysis of rubbersbraided fabric composite hose described in Ref. [27]. In the study in Ref. [27], orthotropic braided fabric is firstly homogenized in the framework of linear elasticity. Afterwards, the homogenized rubbers-braided fabric composite hose is obtained, and nonlinear large deformation analysis is carried out by using the obtained homogenized properties. This technique could be a good alternative particularly when the stress response in the heterogeneous microstructure is not the variable of interest. It is also important to note that the analysis in this paper implements linear elastic assumption to all of the unit-cell constituents. Viscoelasticity, which exhibits time-dependent strain of the resin, is not included. However, in the extended application, consideration of viscoelastic which may create matrix creep can be performed by modifying its constitutive equation [28]. The modification described in Ref. [28] was performed by involving Findley’s and Norton’s equations, and showed a good agreement with the experimental results. 4.3. Analysis sequence ð57Þ Fig. 4 shows the analysis sequence using two models: unit-cell and beam models. The analysis is carried out using an in-house 384 M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 Fig. 3. 3-D orthogonal interlocked composite: (a) schematic, (b) side-view observed from optical microscopy and (c) unit-cell model. Fig. 4. Analysis sequence. code developed in Fortran 90. The unit-cell is used to perform homogenization and localization analyses. In the homogenization analysis, characteristic displacements (or correctors) and homoge- nized thermomechanical properties are obtained while the localization analysis acquires the local stresses within the unit-cell model. The beam model is a macroscopic model used to calculate M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 385 Fig. 5. Region of localization analysis: (a) beam model, (b) five elements of interest and (c) four representative lines. Table 1 Thermomechanical properties of carbon fibers and epoxy resin. x-tow and resin; Line C intersects resin and y-tow; Line D intersects z-tow. Properties Carbon fiber T-300 Carbon fiber T-900 Epoxy EP828 EL (GPa) ET (GPa) GLT (GPa) GTT (GPa) 220 13.8 11.35 5.5 0.20 0.25 4.10  107 5.60  106 294 6.5 18.2 6.5 0.32 0.41 5.60  107 5.60  106 3.4 3.4 1.26 1.26 0.35 0.35 6.45  105 6.45  105 mLT mTT aL (/°C) aT (/°C) the macroscopic displacement (u0). The beam model utilizes homogenized thermomechanical properties obtained from homogenization analysis of the unit-cell model of 3-D orthogonal interlock composites. After the properties have been defined, the beam is then subjected to temperature difference DT of 130 °C. Subsequently, the cross-sectional area of the beam tip is subjected to mechanical loading depending on the case. In Case 1, the beam is subjected to uniaxial tensile loading, while in Case 2, the beam is subjected to bending load. Localization analysis is divided into three steps: (i) calculation of thermal residual stresses, (ii) calculation of stresses due to mechanical loading (tension or bending loads), (iii) calculation of total stresses, which are the superposition of thermal residual stresses and stresses due to mechanical loading. For the purpose of localization analysis, the beam model consisting of 625 elements and 3396 nodes is built (Fig. 5(a)). Five elements in the center of the beam model are of interest, namely Elements 63, 188, 313, 438 and 563 (Fig. 5(b)). Each element represents the whole volume of the unit-cell model. In localization analysis, the calculation of stress requires the macroscopic displacement of each node in the unit-cell. The displacement is obtained by interpolating the macroscopic displacements of the selected elements in the beam model. In order to study the stress distribution in the unit-cell of 3-D orthogonal interlock composites, the stresses are extracted from four representative lines intersecting the constituents of the unitcell, namely Line A, Line B, Line C and Line D as shown in Fig. 5(c). Line A intersects x-tow and y-tow; Line B intersects 5. Results and discussion 5.1. Homogenization analysis 5.1.1. Modes of elastic and thermal correctors The homogenization analysis produces elastic and thermal correctors. Fig. 6 shows the six modes of elastic correctors, which consist of Mode 11, Mode 22, Mode 33, Mode 23, Mode 31 and Mode 12. The correctors are shown in terms of normalized magnitude with reference to the maximum value of each mode. Figs. 6(a)– (f) show that the deformation of opposite faces of each mode is the same in terms of magnitude and direction. Because temperature difference is included, present analysis also produces thermal corrector. Fig. 7 shows the mode of thermal corrector of 3-D orthogonal interlock composites. The typical modes of elastic correctors are similar with the ones displayed in Refs. [5,9]. This shows that the correctors obtained from present results are reasonable albeit different materials understudy. 5.1.2. Homogenized thermomechanical properties The homogenization analysis produces homogenized thermomechanical properties (Young’s moduli (E1, E2, E3), shear moduli (G12, G31, G23), Poisson’s ratios (m23, m31, m12) and coefficients of thermal expansion (a1, a2, a3)). These properties are acquired after the calculation of elastic and thermal correctors. Because the properties are sensitive to the number of fibers within fiber tows, the effect of fiber volume fraction of each tow (Vft) on the properties is discussed in this section. It is important to note that Vft is the fractional volume of fibers within a fiber tow, while Vf is defined as fractional volume of fibers within a whole composite. In Fig. 8, the effect of Vft of x-, y-, and z-tows on the homogenized properties of 3-D orthogonal interlock composites is shown. Some properties have the same value (E1 = E2, G31 = G23, a1 = a2) due to in-plane symmetry of the unit-cell. It is important to note that parametric studies are performed to study the effect of changing the Vft of three tows, while the total volume of each tow is maintained. Figs. 8(a)–(d) show that the increase of Vft of x- and 386 M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 Fig. 6. Elastic correctors {v}: (a) mode 11, (b) mode 22, (c) mode 33, (d) mode 23, (e) mode 31 and (f) mode 12. subjected to low rate compression loading. Three specimens, made by Resin Transfer Molding (RTM) technique, were tested. Each specimen was uniaxially compressed under constant displacement rate of 0.5 mm/min. In the experiment, contribution of fiber volume fractions (Vf) from the fibers within x-, y- and z-tows are 24.3%, 24.7% and 0.5%, respectively. Hence, the corresponding total Vf of the specimen was 49.5%. Numerical results are obtained by varying the Vf of x- and y-tows with following values: 22.5%, 24.75% and 27%. Vf of z-tow is maintained at 0.5%. A considerably good agreement between experiment and numerical analysis is shown in Fig. 9(a). The average value of E11 obtained from the experiment is 59.13 GPa, while that of numerical analysis is 58.53 GPa. The value of 58.53 GPa is calculated from the linear interpolation of homogenization results when Vf of unit-cell is 49.5%. The difference of only 1% indicates that the analysis yields a good result. For the Poisson’s ratio, shown in Fig. 9(b), the difference of 28% between experimental and analysis results is compared, despite the fact that the numerical result almost intersects the lower cap of the experimental error bar. Several factors may contribute to the difference of E11 and v12 between the experimental and analysis results. In the unit-cell modeling, idealization of complex structure of 3-D orthogonal interlock composites is inevitable: (i) the cross-section of fiber tow is rectangular; (ii) fiber undulation is not modeled; (iii) a horizontal portion of z-tow binding the fiber tows is not included; (iv) selvage yarn is not modeled; (v) the unit-cell model of 3-D composites may not capture the variation of thickness or width of fiber tows and other phases. In terms of experiments, particularly in the data reduction of Poisson’s ratio, the strain gages used to calculate the strain were attached on the surface of specimen, which is usually a resin-rich layer. In that case, the placement of strain gages may produce larger deformation, and may not represent the average strain ratio of the whole structure [22]. With regard to the numerical analysis, homogenized Poisson’s ratio v12 is sensitive to the particular utilized fiber properties of x- and y-tows. Further analysis is performed to study the effect of fiber properties on homogenized properties of 3-D composites. Authors found that EL, ET and vLT of the x- and y-tows’ fiber have a significant influence on the homogenized v12. Figs. 10(b) and (c) show that the increase of ET and vLT will increase the homogenized v12. Contrarily, the increase of EL decreases the homogenized in-plane Poisson’s ratio as shown in Fig. 10(a). However, the influence of GLT, GTT and vTT to the homogenized v12 is negligible, even though the fiber portion of x- and ytows within the unit-cell are significant. 5.2. Localization analysis Fig. 7. Thermal corrector {w}. y-tows tends to increase the elastic moduli (i.e. Young’s and shear moduli) while Poisson’s ratio and CTE results shown in Figs. 8(e)– (i) are decreased. The effect of Vft of z-tow is also studied. Fig. 8 shows that E1, E2, G12, G23, G31 and m31 are insensitive to the Vft of z-tow. The increase of Vft of z-tow increases E3, m12, a1 and a2, and reduces m23 and a3. In general, it is found that the effect of Vft of z-tow is significant for E3 and m23, as compared to the effect of Vft of x- and y- tows on those properties. 5.1.3. Comparison between simulation and experiment In Fig. 9, E11 and v12 obtained from homogenization analysis are compared with the experimental results. The experimental results are obtained from Ref. [29]. The experiment was a compression test of composites based on ASTM-D695 (Standard Test Method for Compressive Properties of Rigid Plastics), which covers the determination of mechanical properties of reinforced plastics 5.2.1. Thermal residual stresses In localization analysis, the unit-cell model of 3-D composites is subjected to DT of 130 °C. Six surfaces of the unit-cell model are fully constrained during the application of DT. Due to such constraint, the total strain of homogeneous macroscopic model is zero so that Eq. (54) can be re-expressed as r0ij ðx; yÞ ¼ C ijpq @wp  C ijpq apq DT @yq ð59Þ Therefore, the thermal residual stresses depend only on the microscopic variables. Consequently, the distribution of the thermal residual stresses is the same for the selected elements in the composite beam (elements 63, 188, 313, 438 and 563). In localization analysis, stresses calculation is carried out by using a unit-cell model utilizing homogenized thermomechanical properties as shown in Table 2. Six stresses in the unit-cell were obtained. However, r11 is selected here to represent the dominant stresses. Fig. 11 shows the distribution of thermal residual stress M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 387 Fig. 8. Homogenized thermo-mechanical properties. Fig. 9. Numerical and experimental results (Vf of z-tow = 0.5%): (a) E1, (b) m12. r11 along the Lines A, B, C and D. The vertical axis denotes thickness coordinate in the microscopic scale (y3). In Figs. 11(a)–(d), the first approximation of stresses (r0) and the homogenized stress of unit-cell (rH) are shown. The obtained stresses (r0) indicate that, in general, thermal stresses of 3-D orthogonal interlock composites are influenced by CTE. The region having relatively larger CTE in particular direction than its adjacent region of the same line tends to exhibit higher stresses. For example, a1 of x-tow is 5.03  107/°C while that of resin is 6.45  105/°C. Correspondingly, Fig. 11(b) shows that average r11 of x-tow and resin along Line B are 31.2 MPa and 80.2 MPa, respectively. A different result is seen when comparing Fig. 11(c) and Fig. 11(d) whereby ztow exhibits lower stresses than y-tow despite the fact that a1 of ztow is higher than that of y-tow. Such behavior exists because E1 of y-tow is higher than that of compared to z-tow while the difference of their CTE is insignificant. Figs. 11(b)–(d) show the curved stress pattern along thickness direction. The stress pattern is caused by the complex constituent of the structure with highly different 388 M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 Fig. 10. Effect of x- and y- tows’ fiber properties to the homogenized Poisson’s ratio m12: (a) EL (ELref = 220 GPa), (b) ET (ETref = 13.8 GPa) and (c) mLT (mLTref = 0.20). Table 2 Tows and resin properties. Properties X-tow T-300 (Vft = 55%) Y-tow T-300 (Vft = 55%) Z-tow T-900 (Vft = 50%) Resin region EP828 E1 (GPa) E2 (GPa) E3 (GPa) G12 (GPa) G31 (GPa) G23 (GPa) 122.55 7.13 7.13 3.25 3.25 2.53 0.263 0.015 0.414 5.03  107 3.94  105 3.94  105 7.13 122.55 7.13 3.25 2.53 3.25 0.015 0.414 0.263 3.94  105 5.03  107 3.94  105 4.96 4.96 148.70 2.45 3.21 3.21 0.476 0.335 0.011 4.50  105 4.50  105 2.00  107 3.40 3.40 3.40 1.26 1.26 1.26 0.35 0.35 0.35 6.45  105 6.45  105 6.45  105 m12 m31 m23 a1 (/°C) a2 (/°C) a3 (/°C) Fig. 11. Thermal residual stress r11 (DT = 130 °C): (a) Line A, (b) Line B, (c) Line C and (d) Line D. material properties (a and E) which create non-uniform strain distribution. From Fig. 11, it is clearly shown that resin region experiences the highest tensile stress r11. The maximum stress of resin region is 85.8 MPa. 5.2.2. Stresses due to mechanical loading 5.2.2.1. Case 1: Uniaxial tensile loading. Fig. 12 shows the distribution of r11 due to uniaxial tensile loading. The vertical axis in Fig. 12 denotes the normalized thickness coordinate in macro- M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 scopic scale (x3/h), where h is the thickness of the macroscopic model (h = 2.98 mm). Fig. 12 shows the distribution of r11 in 3-D orthogonal interlock composites obtained by the application of F = 250 N. The first approximation of stresses (r0) and the homogenized stress of unit-cell (rH), as well as the finite element and analytical results are included to identify the difference of localization result as compared to the results of the homogeneous macroscopic model. The r0 and rH are calculated five times by involving the macroscopic displacement (u0) of element 63 to 563 along the thickness direction. The r0 are taken along Lines A, B, C and D. The analytical stress can be simply calculated by dividing the tensile force F with the cross-sectional area. In Figs. 12(a)–(d), FEM, analytical, and rH yield the same value of r11 so that those three results coincide. Figs. 12(a)–(d) also show that stresses in each of five selected elements, which consists of two constituents except for Fig. 12(d), exhibits a repeating pattern through the thickness direction. It is also shown that the region having larger Young’s moduli in x-direction (E1) experiences higher stress. This phenomenon is clearly shown in Figs. 12(a) and (b), where x-tow (E1 = 122.55 GPa) experiences a relatively higher r11 than y-tow (E1 = 7.13 GPa) and resin (E1 = 3.4 GPa). Similar pattern is shown in Figs. 12(c) and (d) whereby y-tow experiences a slightly higher stress than z-tow and resin region. In general, among all constituents, x-tow exhibits the highest r11. The maximum r11 in x-tow obtained by the localization analysis is 11.82 MPa. As described in Section 5.2.1 (thermal residual stress analysis), the curved stress pattern along thickness direction is also clearly found particularly in Fig. 12(b) and (d). This behavior is due to non-uniform strain distribution predominantly caused by the different value of E1 between the adjacent constituents. 5.2.2.2. Case 2: Bending load. Fig. 13 shows the distribution of r11 along four different lines of 3-D composite unit-cell model due to bending load. Fig. 13 shows the stress distribution obtained by application of V = 250 N onto the beam model. In Fig. 13, FEM 389 and analytical results show the normal stress distribution of homogeneous beam due to bending load. The analytical stress is calculated by the following formulation: r11 = VLz/I; where V denotes the bending load, L is the horizontal distance between the center and the tip of beam, z is the vertical distance between point of analysis and neutral axis of the beam, and I is the cross-sectional area moment of inertia. As the consequence of the involvement of u0 calculated from homogeneous beam model, tension and compression stresses exist with the transition in the neutral axis of homogeneous beam (x3/h = 0.5). It is important to note that rH is the homogenized stress, which belongs to each the unit-cell. Fig. 13 shows that localization analysis obtains the typical stress distribution under bending load. Figs. 13(a) and (b) show that xtow experiences higher stress r11 than y-tow and resin due to its higher modulus. Figs. 13(a) and (b) also elucidate that the stresses, both tensile and compressive stresses, are dominantly experienced by x-tow as compared to other constituents. The maximum tensile stress of x-tow is 867.42 MPa while the maximum compressive stress is 687.33 MPa. The relatively low stresses are experienced by y-tow, z-tow and resin region as shown in Figs. 13(c) and (d). 5.2.3. Remarks on stresses under thermal and mechanical loading The stresses are investigated by calculating the normalized stress (ns) with respect to the strength of corresponding direction (see Table 3). Subscripts ‘L’ and ‘T’ in Table 3 denote longitudinal and transverse directions, respectively, whilst ‘t’ and ‘c’ denote tension and compression. In this analysis, the transverse strength of both in-plane and out-of-plane directions are assumed to be equal. The other assumption is that the strength of z-tow has the same values with those of x- and y-tows. The maximum values of ns of each constituent of 3-D orthogonal interlock composites are included in Table 4. The values are based on the stresses along Lines A, B, C and D. Table 4 shows that thermal residual stresses increase the criticality to experience damage initiation of every constituent of 3-D composite even Fig. 12. r11 due to tension loading F = 250 N: (a) Line A, (b) Line B, (c) Line C and (d) Line D. 390 M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391 Fig. 13. r11 Due to bending load V = 250 N: (a) Line A, (b) Line B, (c) Line C and (d) Line D. mechanical loading. In the case of total stress under uniaxial tensile loading, the maximum ns of x- and y-tows are obtained in zdirection (due to r33). Moreover, different direction of maximum ns is also found in the results of the total stresses under bending load. In this case, although r11 in x-tow is very high, the maximum ns is caused by r33 instead due to the low value of FTt. Table 3 Tow and resin strength properties (L and T directions correspond to those of x-tow) [30]. Component FLt (MPa) FLc (MPa) FTt (MPa) FTc (MPa) SLT (MPa) STT (MPa) Tow (T-300, Vft = 60%) Resin region (EP828) 1480 89 1490 – 80 – 280 – 93 – 93 – 6. Conclusions A formulation of homogenization method has been presented in the framework of asymptotic analysis in linear elasticity by including the thermomechanical effect. The formulation is used to obtain homogenized thermomechanical properties of 3-D orthogonal interlock composites, and to perform localization analysis whereby detailed stresses due to mechanical loading and thermal residual stresses are attained. Several conclusions can be summarized from the analysis: when the mechanical loading has not been applied. It is characterized by the high value of maximum ns (nearly 1). The maximum ns due to thermal loading is obtained from the in-plane stresses (r11, r22) in z-tow and resin, while maximum ns of x- and y-tows are obtained on out-of-plane direction. In the case of only mechanical loading, the maximum ns in every constituent of 3-D composite is due to r11. The results also show that, with the same amount of applied load, bending load is more critical than uniaxial tension loading. The effects of thermal residual stresses are also investigated by calculating the total stresses. The total stresses are obtained by the superposition of thermal residual stresses (only DT) and stresses under mechanical loading. The incorporation of thermal residual stresses also creates different behavior of stresses as compared to the results of only  The homogenization method developed in this paper is able to obtain both elastic and thermal correctors of unit-cell of 3-D orthogonal interlock composites. For the elastic correctors, the typical modes obtained from present analysis are similar with those obtained from several references, albeit different materials understudied. Table 4 Maximum values of ns. Loading DT Tension Tension + DT Bending Bending + DT X-tow Y-tow Resin Z-tow Max. ns Stress Max. ns Stress Max. ns Stress Max. ns Stress 0.96 0.01 0.96 0.59 1.10 r33 r11 r33 r11 r33 0.96 0.01 0.96 0.51 1.37 r33 r11 r33 r11 r11 0.84 0.01 0.85 0.47 1.29 r11, r22 r11 r11 r11 r11 0.96 0.01 0.97 0.44 1.31 r11, r22 r11 r11 r11 r11a M.R.E. Nasution et al. / Composites: Part B 60 (2014) 378–391  Homogenized thermomechanical properties of 3-D orthogonal interlock composites can be efficiently obtained using present method. The effects of fiber content within fiber tows (Vft) of 3-D composites, specifically of x-, y- and z-tows, on the thermomechanical properties are also investigated. The investigation shows that Vft of x- and y-tows play an important role in the calculation of equivalent thermomechanical properties, either of in-plane and out-of-plane properties. It is also found that the effect of Vft of z-tow is significant for E3 and m23 of 3-D orthogonal interlock composites.  Tensile modulus (E1) and Poisson’s ratio (m12) obtained from present homogenization method are compared with those obtained from experiments. The difference of around 1% is obtained when comparing the modulus. When Poisson’s ratio is compared between analysis and experiment, larger difference is found. Several factors are considered to be the reason behind the difference between analysis and experiment in terms of Poisson’s ratio: idealization in the modeling (i.e. simple unitcell model), larger deformation due to the placement of strain gages, and the sensitivity of Poisson’s ratio to the particular fiber properties (i.e. x- and y-tows).  Localization analysis shows that the stresses generated in the constituents of 3-D orthogonal interlock composites are greatly dependent on the loading conditions (purely mechanical, purely thermal or combined mechanical-thermal loadings). It is also found that the thermal residual stresses are highly influenced by the CTE of composite constituents.  Localization analysis also suggests that incorporating thermal effects in the analysis of 3-D composite is imperative because the resulting thermal residual stresses may change the damage mechanisms, i.e. damage initiation and the ensuing damage could be accelerated with the presence of thermal loading. Acknowledgement Authors would like to extend their gratitude to the Tokyo Metropolitan Government for the financial support under the project of Asian Network of Major Cities 21 (ANMC-21). References [1] Tan P, Tong L, Steven GP. Behavior of 3D orthogonal woven CFRP composites Part II FEA and analytical modeling approaches. Compos A Appl Sci Manufact 2000;31:273–81. [2] Sun B, Liu Y, Gu B. A unit cell approach of finite element calculation of ballistic impact damage of 3-D orthogonal woven composite. Compos B Eng 2009;40:552–60. [3] Bensoussan A, Lion JL, Papanicolaou G. Asymptotic analysis for periodic structures. Amsterdam: North-Holland Pub. Co.; 1978. [4] Guedes JM, Kikuchi N. Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Comput Methods Appl Mech Eng 1990;83:143–98. [5] Chung PW, Tamma KK, Namburu RR. Asymptotic expansion homogenization for heterogeneous media: computational issues and applications. Compos A Appl Sci Manufact 2001;32:1291–301. 391 [6] Matsuda T, Nimiya Y, Ohno N, Tokuda M. Elastic-viscoplastic behavior of plainwoven GFRP laminates: homogenization using a reduced domain of analysis. Compos Struct 2007;79:493–500. [7] Angioni SL, Meo M, Foreman A. A comparison of homogenization methods for 2-D woven composites. Compos B Eng 2011;42:181–9. [8] Peng X, Cao J. A dual homogenization and finite element approach for material characterization of textile composites. Compos B Eng 2002;33:45–56. [9] Oliveira JA, Pinho-da-Cruz J, Teixeira-Dias F. Asymptotic homogenisation in linear elasticity. Part II: Finite element procedures and multiscale applications. Comput Mater Sci 2009;45:1081–96. [10] Francfort GA. Homogenization and linear thermoelasticity. SIAM J Math Anal 1983;14(4):696–708. [11] Shabana YM, Noda N. Numerical evaluation of the thermomechanical effective properties of a functionally graded material using the homogenization method. Int J Solids Struct 2008;45:3494–506. [12] Dasgupta A, Agarwal RK, Bhandarkar SM. Three-dimensional modeling of woven-fabric composites for effective thermo-mechanical properties and thermal properties. Compos Sci Technol 1996;56:209–23. [13] Vel SS, Goupee AJ. Mutiscale thermoelastic analysis of random heterogeneous materials. Part I: Microstructure characterization and homogenization of material properties. Comput Mater Sci 2010;48:22–38. [14] Pinho-da-Cruz J, Oliveira JA, Teixeira-Dias F. Asymptotic homogenisation in linear elasticity. Part I: Mathematical formulation and finite element modelling. Comput Mater Sci 2009;45:1073–80. [15] Hobbiebrunken T, Fiedler B, Hojo M, Ochiai S, Schulte K. Microscopic yielding of CF/epoxy composites and the effect on the formation of thermal residual stresses. Compos Sci Technol 2005;65:1626–35. [16] Yang L, Yan Y, Ma J, Liu B. Effects of inter-fiber spacing and thermal residual stress on transverse failure of fiber-reinforced polymer-matrix composites. Comput Mater Sci 2013;68:255–62. [17] Yao L, Jiang M, Zhou D, Xu F, Zhao D, Zhang W, et al. Fabrication and characterization of microstrip array antennas integrated in the three dimensional orthogonal woven composite. Compos B Eng 2011;42:885–90. [18] Luo H, Xiong G, Yang Z, Raman SR, Li Q, Ma C, et al. Preparation of threedimensional braided carbon fiber-reinforced PEEK composites for potential load-bearing bone fixations. Part I. Mechanical properties and cytocompatibility. J Mech Behav Biomed Mater 2014;29:103–13. [19] Mouritz AP, Bannister MK, Falzon PJ, Leong KH. Review of applications for advanced three-dimensional fibre textile composites. Compos A Appl Sci Manufact 1999;30:1445–61. [20] Darrigol O. On the necessary truth of the laws of classical mechanics. Stud Hist Philos Mod Phys 2007;38:757–800. [21] Gibson RF. Principle of composite material mechanics. 3rd ed. Florida: CRC Press; 2012. [22] Yudhanto A. Mechanical characteristics and damage mechanisms of stitched carbon/epoxy composites under static and fatigue loads. PhD Thesis. Tokyo Metropolitan University; 2013. [23] Trinquecoste M, Carlier JL, Derre A, Delhaes P, Chadeyron P. High temperature thermal and mechanical properties of high tensile carbon single filaments. Carbon 1996;34(7):923–9. [24] Yudhanto A, Iwahori Y, Watanabe N, Hoshi H. Open hole fatigue characteristics and damage growth of stitched plain weave carbon/epoxy laminates. Int J Fatigue 2012;43:12–22. [25] Mibayashi H. Damage Development Analyses of 3D-CFRP by Using Node Separation Method. Master Thesis. Tokyo Metropolitan University; 2003. [26] Peng X, Guo Z, Du T, Yu WR. A simple anisotropic hyperelastic constitutive model for textile fabrics with application to forming simulation. Compos B Eng 2013;52:275–81. [27] Cho JR, Jee YB, Kim WJ, Han SR, Lee SB. Homogenization of braided fabric composite for reliable large deformation analysis of reinforced rubber hose. Compos B Eng 2013;53:112–20. [28] Pang F, Wang CH. Activation theory for creep of woven composites. Compos B Eng 1999;30:613–20. [29] Watanabe N, Mibayashi H. Experimental & analytical approach for compressive characteristic of 3-D composite. In: Proceedings of 42nd AIAA/ ASME/ASCE/ASC structures, structural dynamics, and materials conference. Seattle; April, 2001. Paper ID A01-25242. [30] Tanaka R. Damage development analysis for woven composite materials. Master Thesis. Tokyo Metropolitan University; 2007.