Thermal–mechanical cell model for unbalanced plain weave woven fabric composites

Composites: Part A 38 (2007) 1019–1037 www.elsevier.com/locate/compositesa

Thermal–mechanical cell model for unbalanced plain weave woven fabric composites James Lua

*

Global Engineering and Materials, Inc., 9 Goldfinch Terrace, East Lyme, CT 06333-1342, USA Received 13 January 2006; received in revised form 25 May 2006; accepted 13 June 2006

Abstract A high fidelity assessment of accumulative damage of woven fabric composite structures subjected to aggressive loadings is strongly reliant on the accurate characterization of the inherent multiscale microstructures and the underlying deformation phenomena. The stress and strain fields predicted at a global structural level are unable to determine the damage and failure mechanisms at the constituent level and the resulting stiffness degradation. To establish a mapping relation between the global and constituent response parameters, a new four-cell micromechanics model is developed for an unbalanced weave subjected to a thermal–mechanical loading. The developed four-cell micromechanics model not only bridges the material response from one length scale to another but also quantifies the composite thermal–mechanical properties at a given state of constituent damage. The thermal–mechanical mapping relations at different microstructural levels are derived based on the multicell homogenization, intercell compatibility conditions, and energy methods. Because of the high computational efficiency of the developed thermal–mechanical micromechanics model, it can be linked with a finite-element-based dynamic progressive failure model, where the response parameters at different microstructural levels can be extracted for each Gaussian point and at each time step. The accuracy and the dual function of the developed micromechanics model are demonstrated with its application to a balanced plain weave, an unbalanced plain weave, and failure mode simulation of a tensile coupon test.  2006 Elsevier Ltd. All rights reserved. Keywords: A. Fabrics/textiles; B. Physical properties; C. Micro-mechanics; D. Mechanical testing

1. Introduction Fiber-reinforced polymer (FRP) composites have been increasingly used in the current and next generation of naval ships (e.g., composite hull, composite topside, composite mast, composite sail, hybrid hull, etc.). The most compelling reasons to use composite materials in naval ship components are stealth, lower total ownership cost, and weight reduction. The composite materials have played key roles in reducing the magnetic, acoustic, hydrodynamic, radar, and thermal signatures and increasing pay-

*

Tel.: +1 860 398 5620; fax: +1 860 398 5621. E-mail address: [email protected]

1359-835X/$ - see front matter  2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compositesa.2006.06.023

load, top speed, and operation range. Preliminary design constraints for fiber-reinforced polymer (FRP) composite materials are mainly based on a test-driven certification procedure for which a large number of coupon tests have to be performed at various room/wet/hot conditions. A significant amount of resources will be invested in the material certification phase. The prohibitive nature of extensive component testing has led to the investigation of using a virtual testing tool to simulate the material behavior and structural response from small-scale and coupon-level test data. Woven fabric composites have unique tailoring capabilities for improving the ballistic resistance and high strength of naval ship structural components. However, the complex architecture of woven fabric composites makes response and failure prediction very difficult and computationally intensive. The limited understanding of the

1020

J. Lua / Composites: Part A 38 (2007) 1019–1037

presence of multiple failure modes, damage patterns, and their synergistic interactions toward the final rupture has resulted in a test-driven material and structural certification. Using a conservative ‘‘knock-down’’ factor based on limited tests adds unnecessary weight that will have a significant impact on the payload, top speed, and operation range of a ship. To reduce the number of full-scale tests required and enhance the level of design confidence, a reliable response and damage assessment tool has to be developed to capture the constituent failure mechanism and propagate the damage from one length scale to another. An accurate assessment of the accumulative damage to a composite ship structure subjected to an aggressive loading environment strongly relies on the accurate characterization of the inherent multiscale microstructures and the underlying deformation phenomena. The stress and strain fields predicted at a global structural level are unable to determine the damage and failure mechanisms at the constituent level and the resulting stiffness degradation. To capture the synergistic interaction of constituent failure modes that leads to material degradation and final rupture, the stress/strain response at each microstructural level has to be determined from given macro stress/strain response parameters. Given the constituent stress and strain and the associated mechanism-driven failure criteria, the material degradation at each constituent can be quantified. Finally, to propagate the damage from the constituent to the structural level, a multilevel homogenization has to be performed. Therefore, to implement a multiscale progressive failure analysis, a dual-function micromechanics model has to be developed to (1) decompose the response parameters from the global structural level to microstructural levels for the failure mode prediction and (2) homogenize damaged material properties at the global structural level based on the damaged material properties at microstructural levels for the damage prediction at the next loading stage. As both the decomposition and homogenization procedures have to be performed at each Gaussian point and at each time step during a finite-element-based dynamic progressive failure analysis, it is imperative to develop an efficient micromechanics model while retaining sufficient accuracy. A large body of work has been performed to develop micromechanics models of a woven fabric composite. Many of the existing analysis models can be classified as either a simplified structural representation model or a detailed finite-element-based micromechanics model. Several structural representation models have been developed such as the mosaic model [1], the fiber-undulation model [2], the bridging model [3], the binary model [4], the laminated strap model [5], and the cell model [6]. Because these simplified micromechanics models are unable to accurately capture the detailed microstructures and the associated stress gradient, various finite-element-based micromechanics models have been developed for a representative volume cell of a woven fabric composite [7–10]. The main focus of these micromechanics models is to quantify the effective

material properties of a balanced weave for given constituent properties (homogenization). In addition, the thermal residual stress resulting from the thermal mismatch of coefficients of thermal expansion (CTE) of the constituent material has been ignored. As shown by Key and Lua [11] recently, the temperature-induced thermal residual stress in each constituent of a woven fabric composite plays a detrimental role in failure initiation and damage progression. To accurately capture the geometric configuration of the woven fabric architecture, a B-spline representation of the yarn coupled with the meso-volume and sub-element partition has been developed by Bogdanovich and Pastore [12]. A system of basis spline functions is used to describe the discontinuity of stress field at the constituent interface. The heterogeneity induced spatial variability in mechanical properties has been studied for different types of textile composites. Recognized the presence of a hierarchy of structural and scale levels in textile composites, advanced modeling strategies have been developed recently and implemented in an integrated analysis and design tool [13–15] to establish a mapping sequence from fiber, yarn, textile, perform, to the final composite material system. A finite-element-based thermal–mechanical micromechanics model (TMAT) has been developed recently by Lua and O’Brien [16] based on a balanced weave configuration. Given the thermal–mechanical properties of each constituent along with their microstructural parameters, a finite-element model (FEM) associated with a commercial FEM solver (NASTRAN/LSDYNA3D) is generated to quantify the orthotropic elastic properties, coefficients of thermal expansions, and thermal conductivities. The accuracy of TMAT predictions has been demonstrated using composite material strain measurements from a digital image correlation system [17]. To assess the fire damage of a composite structure, TMAT has been used to generate the thermal–mechanical properties of a plain weave at various pyrolysis stages [18]. The finite-element-based micromechanics model can capture the detailed woven fabric architecture and incorporate the nonlinear constitutive behavior for each constituent. However, it is computationally intensive and not well suited for use in the dynamic progressive failure analysis where the finite-element-based micromechanics model has to be executed at each Gaussian point and at each loading stage. Using oversimplified micromechanics models cannot accurately capture the microstructures and the underlying deformation phenomena, which are crucial for the mechanism-driven failure prediction. Recognizing that the energy dissipates by various damage or fracturing processes at constituent levels, the use of constitutive stress and strain in a mechanism-driven failure criterion has become the recent subject of many research works [19–23]. To implement the constituent-based progressive failure analysis procedure, a dual-function micromechanics model has to be developed with an emphasis on high computational efficiency and sufficient accuracy.

J. Lua / Composites: Part A 38 (2007) 1019–1037

The main objective of this study is to develop an accurate and computationally efficient micromechanics model for (1) quantification of the thermal–mechanical composite material properties of an unbalanced weave at a given constituent damage state and (2) establishment of the mapping relation of the stress/strain response between the microstructural and macrostructural levels. To enhance the computational efficiency while retaining sufficient accuracy for response prediction at the constituent level, a dual-function micromechanics model based on the four-cell representation is developed to bridge the thermal–mechanical behavior of an unbalanced plain weave between the micro- and macro-levels. While the four-cell micromechanics model has been developed by Tabiei and Ivanov [22], it is only applicable for the mechanical loading because it does not address the temperature-dependent thermal residual stress. In addition, due to the use of antisymmetry in the model development, it cannot be directly applied for the unbalanced weave. To characterize the thermal–mechanical properties of a multicell system, an N-cell homogenization model is given in Section 2 based on the imposed iso-strain and iso-stress boundary conditions. Given the four-cell representation of the unit cell of an unbalanced weave shown in Section 3, the multicell homogenization procedure is then applied in Section 4 for each of these representative cells. With the established thermal–mechanical constitutive behavior of each cell, a mapping relation between the combined individual cell response and the unit cell response is established in Section 5 for solving unknown cell response variables from given unit cell response parameters. The numerical implementation of the dual-function micromechanics model is summarized in Section 6 followed by example applications in Section 7. 2. Multicell homogenization Consider a representative cell that consists of N subcells as shown in Fig. 1, where each subcell is assumed to be an orthotropic ply. The constitutive behavior of the kth sub-

1021

cell is described by an orthotropic thermal–mechanical model in its local material coordinate system by 8 0 9 3 2 k rxx > > C 11 C k12 C k13 0 0 0 > > > > > > 6 Ck Ck Ck > r0yy > > > 0 0 0 7 > > 7 6 12 22 23 > > > 0 > 7 6 k = < k k rzz 6 C 13 C 23 C 33 0 0 0 7 7 6 ¼ 6 0 > r0yz > 0 0 C k44 0 0 7 > > 7 6 > > > > 7 6 > k > > 5 4 0 r0zx > 0 0 0 C 0 > > 55 > > > ; : r0 > k 0 0 0 0 0 C 66 xy k 9 8 0 k exx  a1 DT > > > > > > > > > e0yy  ak2 DT > > > > > > > > = < e0  ak DT > zz 3  ðk ¼ 1; . . . ; NÞ ð1Þ 0 > > cyz > > > > > > > > > > c0zx > > > > > > ; : 0 c xy

k

0

where the prime ( ) stands for the response parameters associated with the local material coordinate system. In Eq. (1), C kij and aki denote the stiffness coefficients and the coefficients of thermal expansion of the kth subcell, respectively. To establish a mapping relation between the material response at the cell and its associated subcell levels, the stress and strain components defined in the global cell coordinate system shown in Fig. 1 are divided into two groups in terms of iso-strain and iso-stress. Based on the laminated orthotropic ply structure shown in Fig. 1, the strain of the kth subcell in the global cell coordinate system is assumed to be represented by fegk ¼ f fen gTk

fes gTk gk

T

ð2aÞ

cxy gTk czx gTk

ð2bÞ

fen gk ¼ f exx

eyy

fes gk ¼ f ezz

cyz

ð2cÞ

where the superscript T in Eqs. (2a)–(2c) stands for the transpose operator, and {en}k is the assumed iso-strain of the kth subcell. Similarly, the stress of the kth subcell in the global cell coordinate system is assumed to have the following form frgk ¼ f frn gTk

T

ð3aÞ

rzx gTk rxy gTk

ð3bÞ

T frs gk gk

frs gk ¼ f rzz

ryz

frn gk ¼ f rxx

ryy

ð3cÞ

where {rs}k defined in Eq. (3b) is the assumed iso-stress of the kth cell. Given the definition of the iso-strain (Eq. (2b)) and iso-stress (Eq. (3b)), the thermal–mechanical constitutive relation of the kth cell in the global cell coordinate system can be rewritten as      ðrn Þk ðen Þk  ðetn Þk ½C nn k ½C ns k ¼ ð4Þ ðrs Þk ½C sn k ½C ss k ðes Þk  ðets Þk Fig. 1. A multicell block consisting of N subcells.

where ðetn Þk and ðets Þk are the stress-free thermal strains that are grouped based on the iso-strain definition.

1022

J. Lua / Composites: Part A 38 (2007) 1019–1037

By applying the same iso-strain and iso-stress grouping, the constitutive model of the homogenized representative cell shown in Fig. 1 in its global coordinate system can be written in terms of its volume average stress f rg and volume average strain feg as #   "   t n r ð rn Þ ðen Þ ½C nn  ½C ns  ð5Þ ¼ þ ts r ð rs Þ ðes Þ ½C sn  ½C ss    T T are the thermal-induced residual In Eq. (5), ½ rts  rtn  ½ stress components. Among these four 3 · 3 matrices, ½C nn  and ½C ss  are symmetrical matrices and ½C sn  ¼ ½C ns T . To represent f½C ij ; i; j ¼ n; sg as a function of the four 3 · 3 stiffness matrices of each subcell {[Cnn]k; [Cns]k; [Csn]k; [Css]k; k = 1, 2, . . . , N}, both the rule of mixtures and the compatibility conditions posed on the iso-stress and iso-strain are employed. Given the volume fraction of each subcell Vk (k = 1, 2, . . . , N), the volume average strain and stress can be determined by feg ¼

N X

V k fegk

ð6Þ

k¼1

f rg ¼

N X

V k frgk

ð7Þ

k¼1

where the Di (i = 1, 2, . . . , 4) matrices are defined by ½D1  ¼

N X

1

V k ð½C nn k  ½C ns k ½C ss k ½C sn k Þ

½D2  ¼

N X

1

V k ð½C ns k ½C ss k Þ

½D3  ¼

N X

1

V k ð½C ss k Þ

ð16Þ

V k ð½C ss 1 k ½C sn k Þ

ð17Þ

k¼1

½D4  ¼

N X k¼1

Finally, to establish the mapping relation among the cell stiffness matrices and the stiffness matrices of associated subcells, Eqs. (12) and (13) are used to cast the functional form of Eq. (5). Solving for ½ rs  from Eq. (13) and substituting the result into Eq. (12), the final cell constitutive relation can be written in the form of Eq. (5). The four 3 · 3 stiffness matrices shown in Eq. (5) along with the thermal stress components are defined here. 1

½C nn  ¼ ð½D1  þ ½D2 ½D3  ½D4 Þ ½C ns  ¼ ð½D2 ½D3  Þ

Applying the compatibility conditions associated with the iso-strain and iso-stress, we have

½C sn  ¼ ð½D3  ½D4 Þ

fen g ¼ fen gk ;

k ¼ 1; 2; . . . ; N

ð8Þ

½C ss  ¼ ½D3 

f rs g ¼ frs gk ;

k ¼ 1; 2; . . . ; N

ð9Þ

þ

ð11Þ

After determining the strain and stress response of each subcell, the unknown strain and stress at the global cell level can be determined by substituting Eqs. (10) and (11) into Eqs. (6) and (7), respectively. The resulting cell stress and strain are given by f rn g ¼ ½D1 fen g þ ½D2 f rs g 1

ð12Þ

V k ð½C ss k ½C sn k fetn gk þ fets gk Þ

ð13Þ

V k ð½C ns k ½C ss k ½C sn k  ½C nn k Þfetn gk

k¼1

rs g  ½D4 fen g fes g ¼ ½D3 f þ

N X k¼1

1

N X

N X

h i t t V k ð½C ss 1 ½C  Þfe g þ fe g sn n k s k k k

h i t V k ð½C ns k ½C ss 1 ½C   ½C  Þfe g sn nn n k k k k

1

 ½C nn k Þfetn gk

N X

ð21Þ

ð22Þ

k¼1

þ ½C ns k ½C ss k f rs g þ ð½C ns k ½C ss k ½C sn k

þ

1

The cell-level residual thermal stress components that result from the mismatch of coefficients of thermal expansion of individual subcells are determined by

ð10Þ

1 ½C ns k ½C ss k ½C sn k Þfen g 1

ð20Þ

k¼1

1

þ ½C ss k ½C sn k fetn gk þ fets gk frn gk ¼ ð½C nn k 

ð18Þ ð19Þ

1

f rtn g ¼ ½C ns 

rs g  ½C ss 1 en g fes gk ¼ ½C ss 1 k f k ½C sn k f

ð15Þ

k¼1

1

By substituting Eqs. (8) and (9) into Eq. (4) and solving for the unknown constituent strain {es}k and constituent stress {rn}k, we have

ð14Þ

k¼1

f rts g ¼ ½D3 1

N X

t t V k ð½C ss 1 k ½C sn k fen gk þ fes gk Þ

ð23Þ

k¼1

The cell homogenization model derived, based on the rule of mixture and the posed iso-strain and iso-stress conditions within the cell, can be used to (1) determine the cell effective thermal–mechanical properties from subcell constituent properties and (2) compute the subcell stress/strain for failure mode prediction from a given cell-level stress/ strain. Because of this dual functionality, the cell homogenization model can be used as a material analyzer and a stress/strain decomposition module. For given thermal– mechanical constituent properties of each subcell, the effective cell properties can be determined from Eq. (5) along with Eqs. (18)–(21). Given the cell stress ð rÞ/strain ðeÞ, the resulting decomposed stress/strain for each subcell can be calculated from Eqs. (8)–(11). The developed cell

J. Lua / Composites: Part A 38 (2007) 1019–1037

homogenization model is applied next to each building block of the four-cell representation of an unbalanced plain weave. 3. A four-cell representation of an unbalanced weave The commonly used naval composites are unbalanced plain weaves fabricated using vacuum assisted resin transfer molding (VARTM). The composite panels produced by Seemann Composites Inc. have the warp tows spaced at 5.1 mm (5 tows per inch) and fill tows spaced at 6.4 mm (4 tows per inch) as shown in Fig. 2. In a previous approach [16–18], an FEM-based micromechanics tool (TMAT) was used by replacing the unbalanced fabric with an equivalent balanced one in evaluating its thermal– mechanical properties (see Fig. 2). To account for the difference in X-tow and Y-tow spacing, a higher fiber volume fraction within a tow is used for the tows with closer spacing. Despite the versatility of TMAT in characterizing the tow path and tow cross-sectional area, the computational effort associated with a dynamic progressive failure analysis is very intensive which prevents its use for the failure prediction of a larger scale composite structure. To truly capture the unbalanced weave architecture while retaining sufficient accuracy, the cell homogenization model derived in Section 2 is used for the micromechanics characterization of the unbalanced weave. Given the complexity of the woven fabric architecture, the multicell homogenization model cannot be directly applied to the representative volume element (RVE) of the unbalanced weave. To remedy this difficulty, the RVE of the

1023

plain weave is further divided into four basic building blocks. The cell decomposition procedure is illustrated here using the TMAT micromechanics model of a balanced weave. As shown in Fig. 3, the RVE is first simplified by using quarter-symmetry. Next, the 1/4th RVE is further decomposed into the XY-, XR-, YR-, and RR-cells. The XY-cell consists of primarily the X- and Y-tow with a small volume fraction in resin (R). In XR- and YR-cells, the tow (X or Y) volume fraction and the resin volume fraction (R) is assumed the same. For the RR-cell, 100% of the material is resin (resin pocket). These four individual cells form the basic building blocks of the 1/4th unit cell model. To characterize each of these four basic building blocks by the cell homogenization model, two simplifications are used as follows: (1) the tow cross-sectional shape is assumed to be rectangular; and (2) the undulated X- and Y-tow is replaced with a straight line segment. For a typical woven fabric composites with a flat tow cross-sectional shape and low waviness ratio of its tow path, the resulting error from using these two assumptions will be small. The final simplified four-cell representation of the unbalanced weave is depicted in Fig. 4. Four microstructural parameters are used to characterize the weave architecture: (1) tow thickness (Htow); (2) tow volume fraction within the 1/4th unit cell (Vtow); (3) the X-tow spacing (Ly); and (4) the Y-tow spacing (Lx). The size of the 1/4th unit cell in dimensionless space is defined by Ly ¼ Ly =Ly ¼ 1:0, Lx ¼ Lx =Ly , Lf ¼ Lf =Ly , and H ¼ H =Ly . Using these dimensionless parameters, the length of the straight portion of either the X- or Y-tow (Lf ), the X-tow angle (hx) in the

Fig. 2. Comparison of balanced and unbalanced weave representation of a Seemann composite system.

1024

J. Lua / Composites: Part A 38 (2007) 1019–1037 V1 L7

Unit Cell

C1 G5

V1

1/4th Unit Cell

L7

Resin (R)

C1 G9

Z

X-tow (X)

Y X

Y-tow (Y)

Z Y X

V1 V1

V1

L7

L7

C1

L7 C1 G18

C1

G16

G17

Z

Z

Y X

Y X

Z Z

XY-subcell

XR-subcell

Y X

Y X

YR-subcell

RR-subcell

Fig. 3. Illustration of basic steps in cell decomposition.

z Resin

Resin

x

=Lf . As can be seen from Fig. 4, except for the RRV RVE f cell, each cell consists of two material phases. Using the multicell homogenization procedure given in Section 2, a mapping relation between each subcell and its associated constituents is established to characterize effective thermal–mechanical behavior.

x

X-To w

H=2Htow

Lx

Y-T ow

y

y

x

Ly

4. Thermal–mechanical homogenization of selected subcells RR y

XR Lf

XR-subcell, and the Y-tow angle (hy) in the YR-subcell can be determined by

To illustrate the homogenization procedure, the derivation of the thermal–mechanical constitutive model of the XY-cell and the XR-cell is given in detail below. The thermal–mechanical models of the YR- and RR-cell are summarized in Appendix A. Because the cell-coordinate system of both the XY- and RR-cell is selected to be the same as the global unit cell coordinate system shown in Fig. 4, the multicell homogenization is applied directly in the global unit cell coordinate system with no coordinate transformation. However, due to presence of the inclined tow in both the XR- and YR-cell, the multicell homogenization is employed in its local cell coordinate system followed by a coordinate transformation to the global cell coordinate system.

Lf ¼ 2Lx V tow =ð1 þ Lx Þ

4.1. XY-cell

XY

YR

Lx

Lf

Ly RR

XR Lx L f

YR

XY

1 Lf

Lf

Lf

Lx = Lx / L y

L y = L y / L y = 1.0

Fig. 4. Four-cell representation of an unbalanced plain weave.

1

ð24Þ

hx ¼ tan ½H =4ðLx  Lf Þ

ð25Þ

hy ¼ tan1 ½H =4ð1  Lf Þ

ð26Þ

Given the overall fiber volume fraction in RVE of V RVE , f the fiber volume fraction within a tow can be calculated by ¼ V tow f

2V RVE Lx f Lf ð1 þ Lx Þ

ð27Þ

For a balanced weave, we have Lx ¼ Lx =Ly ¼ 1:0. Using Eqs. (24) and (27), we have Lf ¼ V tow and V tow ¼ f

As shown in Fig. 4, the XY-cell is filled completely with the straight part of the X-tow and Y-tow with the same subcell volume fraction. Unlike the balanced weave, a distinct thermal–mechanical constitutive model is used to characterize both the X- and the Y-tow. For the XY-cell, it is assumed that the global unit cell coordinate system shown in Fig. 4 coincides with the XY-cell coordinate system. Given the orthotropic behavior of the X-tow and Y-tow defined in the material coordinate system, the

J. Lua / Composites: Part A 38 (2007) 1019–1037

constituent behavior of the X- and Y-tow in the XY-cell coordinate system is defined by 9 8 x 9 2 x 38 x exx  ax1 DT > rxx > C 11 C x12 C x12 0 > 0 0 > > > > > > > > > x > > x > x 7> 6Cx Cx Cx > > > > e r  a DT > > > > 0 0 0 yy yy 2 > > > 7> 6 12 22 23 > > > > > > > > 7 6 x x < rx = 6 C x C x C x < 0 0 0 7 ezz  a3 DT = zz 12 23 22 7 6 ¼ > > cxyz rx > 6 0 0 0 C x44 0 0 7 > > > 7> > 6 > > yz > > > > > > 7 6 x x x > > > > > > > > 5 4 r c 0 0 0 0 C 0 > > > zx > zx 55 > > > > > > > > ; : rx ; : x x cxy 0 0 0 0 0 C 55 xy 8 y 9 r > > > > > > xx > y > > > r > > yy > > > > > = < ry > zz

> ryyz > > > > > > y > > > > r > > zx > > > > ; : ry > xy

2

C y22 6Cy 6 12 6 y 6 C 23 6

¼6 6 0 6 4 0

0

C y12 C y11 C y12

C y23 C y12 C y22

0 0

0 0

0

0

0

0

0 0

0 0

C y55 0 0 C y44 0

0

ð28Þ 9 38 y y exx  a2 DT > 0 > > > > > > > y y 7 > > e  a DT > 0 7> 1 yy > > > > > 7> y y < 7 ezz  a3 DT = 0 7 > cyyz 0 7 > > 7> > > > 7> > > y > c 0 5> > > zx > > > > ; : y y cxy C 55 ð29Þ

where the superscripts x and y in Eqs. (28) and (29) denote the parameter associated with the X-tow and Y-tow, respectively, and the subscripts 1, 2, and 3 given in the stiffness matrix [C] and the coefficient of thermal expansion vector {a} represent the property variable defined along the tow (1) and transverse-to-the-tow (2, 3) direction. Using the iso-strain and iso-stress definitions given by Eqs. (2b) and (3b), Eqs. (28) and (29) can be written in the form of Eq. (4) where the stiffness matrices of the Xand Y-tow are given by " #   ½C nn y ½C ns y ½C nn x ½C ns x ½Cx ¼ ; ½Cy ¼ ð30Þ ½C sn y ½C ss y ½C sn x ½C ss x Substituting these stiffness matrices into Eqs. (14)–(17) and using the subcell volume fraction Vx = Vy = 0.5, the Di (i = 1, 2, . . . , 4) matrices for the XY-cell can be obtained. Finally, using Eqs. (18)–(21), the effective stiffness matrices of the XY-cell can be determined from its constituent stiffness matrices {[C]x, [C]y}. The resulting thermal–mechanical constitutive relation of the XY-cell can be expressed as 3 2 11 8 9 11 rxx > 0 0 0 ½C nn XY ½C 12 > nn XY ½C ns XY > > 7 6 12 > > > 21 6 ½C nn XY ½C 22 > > ryy > 0 0 0 7 > > nn XY ½C ns XY > > 7 6 > > > = > 7 6 11 11 < 7 6 ½C ns XY ½C 21 rzz  ½C  0 0 0 ns XY ss XY 7 ¼6 7 6 33 > rxy > > 0 0 ½C nn XY 0 0 7 > 6 0 > > > > 7 6 > > > 6 0 > ryz > > 0 0 0 ½C 22 0 7 > > 5 4 ss XY > > : ; 33 rzx XY 0 0 0 0 0 ½C ss XY 8 t 9 8 9 e r xx > xx > > > > > > > > > > > > > > t > > > > > e r yy > > > yy > > > > > > > > > > > = > t > < < ezz rzz =  þ ð31Þ > > cxy > 0 > > > > > > > > > > > > > > > > > >c > > > > > 0 > > > > > yz > > > > : > : ; ; 0 XY czx XY

1025

where the thermal-mismatch-induced residual stress comT ponents f rtxx rtyy rtzz g are determined from Eqs. (22) and (23). In Eq. (31), the stiffness coefficients are given by " # 1 ðC x12  C y23 Þ2 y 11 x C 11 þ C 22  ½C nn XY ¼ ð32aÞ 2 C y22 þ C x22   1 ðC x23  C y12 ÞðC y23  C x12 Þ y 12 x ð32bÞ ½C nn XY ¼ C 12 þ C 12 þ 2 C y22 þ C x22 " # 1 ðC y  C x23 Þ2 C x22 þ C y11  12y ½C 22 ð32cÞ nn XY ¼ 2 C 22 þ C x22 1 x y ½C 33 nn XY ¼ ½C 55 þ C 55  2 C x12 C y22 þ C y23 C x22 ½C 11  ¼ ns XY C y22 þ C x22 x C 23 C y22 þ C y12 C x22 ½C 21  ¼ ns XY C y22 þ C x22 2C x22 C y22 ½C 11  ¼ ss XY C y22 þ C x22 2C x44 C y55 ½C 22  ¼ ss XY C y55 þ C x44 2C x55 C y44 ½C 33  ¼ ss XY C y44 þ C x55

ð32dÞ ð32eÞ ð32fÞ ð32gÞ ð32hÞ ð32iÞ

In Eq. (31), the thermal residual stress components resulting from a temperature increase, DT, are calculated by  x y   C 12 C 22 þ C y23 C x22 C x12 x C y12 y C x23 x C y23 y y x y x x a1 þ y a1 þ x a2 þ y a2 þ a3 þ a3 C 22 þ C 22 C 22 C 22 C 22 C 22 ! ! 2 2 C x12 x C y23 y y x þ C 11  x a1 þ C 22  y a2 C 22 C 22    C y23 C y12 y Cx Cx y x þ C 12  ð33aÞ a1 þ C 12  12 x 23 ax2 y C C 22  x 22 y  y y y x  x x DT C 23 C 22 þ C 12 C 22 C 12 x C 12 y C 23 x C 23 y y x ¼ y x x a1 þ y a1 þ x a2 þ y a2 þ a3 þ a3 2 C 22 þ C 22 C 22 C 22 C 22 C 22 ! ! 2 2 C y12 y C x23 x y x þ C 11  y a1 þ C 22  x a2 C 22 C 22    x x C C Cy Cy þ C x12  12 x 23 ax1 þ C y12  23 y 12 ay2 ð33bÞ C C 22  x22 y  x  y y x C 22 C 22 C 12 x C 12 y C 23 x C 23 y ¼ DT a þ a þ a þ a þ ax3 þ ay3 C y22 þ C x22 C x22 1 C y22 1 C x22 2 C y22 2

½rtxx XY ¼ 

½rtyy XY

½rtzz XY

DT 2

ð33cÞ Thus, given the volume average strain {e}XY and the uniform temperature increase DT in the XY-cell, the resulting volume average stress in the XY-cell is computed from Eqs. (31)–(33). With the aid of the volume average stress and strain in the XY-cell [r, e]XY, the stress and strain in its constituents, X- and Y-tow {[r, e]X; [r, e]Y}, can be determined from Eqs. (8)–(11). Therefore, a thermal–mechanical mapping relation between the XY-cell and its constituents (X/Y tow) is established. 4.2. XR-cell Different from the XY-cell, the XR-cell consists of an inclined X-tow with the tow angle of hx and the isotropic

1026

J. Lua / Composites: Part A 38 (2007) 1019–1037

resin R (see Fig. 4). Given the tow and resin configuration of the XR cell shown in Fig. 4, the volume fraction of the X-tow (X) and the resin (R) within the XR-cell can be shown to be the same. Because of the orthotropic behavior of the X-tow associated with its local material coordinate system, two steps are used in deriving the thermal–mechanical behavior of the XR-cell. The multicell homogenization is first applied in the local material coordinate system associated with the XR-cell followed by the coordinate transformation into the unit-cell coordinate system (x–y–z) shown in Fig. 4. The local material coordinate system of the XR-cell (x 0 –y 0 –z 0 ) is obtained from the rotation of the x–y–z system about the positive y-axis with an angle hx. With respect to the local XR-cell coordinate system (x 0 –y 0 –z 0 ), the constitutive model of the X-tow is given by Eq. (28) while the constitutive model of the matrix is defined by 90 8 x 90 2 m 38 m exx  am1 DT > rxx > > C 11 C m12 C m12 0 0 0 > > > > > > > > > m > > x > m 7> 6Cm Cm Cm 0 > > > > e r  a DT > > > > 0 0 yy yy 1 > > > > 7 6 12 11 12 > > > > > 6 m > > x = 7> m m m m < < rzz 6 C 12 C 12 C 11 0 0 0 7 ezz  a1 DT = 7 6 ¼ > 6 > > rxyz > cmyz 0 0 C m44 0 0 7 > > > 7> 6 0 > > > > > > > > 7 6 m x > m > > > > > > > 5 4 r c 0 0 0 0 C 0 > > > > zx zx 44 > > > > > > > ; ; : rx > : m m c 0 0 0 0 0 C xy xy 44 ð34Þ where the superscript ( 0 ) denotes the quantity defined in the local coordinate system. By applying the same multicell homogenization procedure used for the XY-cell, the thermal–mechanical constitutive model of the XR-cell in its local coordinate system can be written as 2 11 3 8 90 11 rxx > 0 0 0 ½C nn XR ½C 12 > nn XR ½C ns XR > > 6 12 7 > > 21 > 6 ½C nn XR ½C 22 > > ryy > 0 0 0 7 > > nn XR ½C ns XR > > 6 7 > = > < 6 ½C 11  ½C 21  ½C 11  rzz 0 0 0 7 6 ns XR 7 ns XR ss XR ¼6 7 33 > > 6 0 7 r 0 0 ½C  0 0 xy > > nn XR > > 6 7 > > > > 6 22 > 0 0 0 ½C ss XR 0 7 > ryz > > 4 0 5 > > : ; rzx XR  0 0 0 0 0 ½C 33 ss XR 8 t 90 8 90 rxx > exx > > > > > > > > > > rt > > eyy > > > > > > > > > > > yy > > > > > > < rt > = = zz zz þ ð35Þ  > > > > > > > 0 > > cxy > > > > > > > > > > > > > 0 > c > > > > > > > : > : yz > ; ; 0 XR czx XR

where the coefficients of the stiffness matrix defined in Eq. (35) are given by ½C 11 nn XR ½C 12 nn XR

" # 2 1 ðC x12  C m12 Þ x m C 11 þ C 11  ¼ 2 C x22 þ C m11   1 ðC x  C m12 ÞðC m12  C x23 Þ ¼ C x12 þ C m12 þ 12 2 C m11 þ C x22

ð36aÞ ð36bÞ

½C 22 nn XR

" # 1 ðC m12  C x23 Þ2 x m C 22 þ C 11  ¼ 2 C x22 þ C m11

1 x m ½C 33 nn XR ¼ ½C 55 þ C 44  2 C x12 C m11 þ C m12 C x22 ½C 11 ns XR ¼ C x22 þ C m11 C x23 C m11 þ C m12 C x22 ½C 21 ns XR ¼ C m11 þ C x22 2C x22 C m11 ½C 11 ss XR ¼ C m11 þ C x22 2C x44 C m44 ½C 22 ss XR ¼ C m44 þ C x44 2C x55 C m44 ½C 33 ss XR ¼ C m44 þ C x55

ð36cÞ ð36dÞ ð36eÞ ð36fÞ ð36gÞ ð36hÞ ð36iÞ

In Eq. (35), the thermal residual stress components resulting from a temperature increase, DT, in the local coordinate system are calculated by

(  m   DT C x23 C m11 þ C m12 C x22 C x12 x C x23 x 2C 12 m x ¼ m x x a1 þ x a2 þ m þ 1 a1 þ a3 2 C 11 þ C 22 C 22 C 22 C 11 !  x x x2 C C C þ C x12  12 x 23 ax1 þ C x22  23 ax C 22 C x22 2 ! ) 2 2C m þ C m11 þ C m12  m12 am1 ð37aÞ C 11 (  m   DT C x12 C m11 þ C m12 C x22 C x12 x C x23 x 2C 12 m x ½rtyy 0XR ¼  a þ a þ þ 1 a þ a 1 3 2 C m11 þ C x22 C x22 1 C x22 2 C m11 !  2 Cx Cx Cx þ C x11  12 ax þ C x12  12 x 23 ax2 C x22 1 C 22 ! ) m2 2C þ C m11 þ C m12  m12 am1 ð37bÞ C 11  m  x m  x  x C 22 C 11 C 12 x C 23 x 2C 12 ½rtzz 0XR ¼ DT a þ a þ þ 1 am1 þ ax3 C m11 þ C x22 C x22 1 C x22 2 C m11

½rtxx 0XR

ð37cÞ After homogenizing the XR-cell in its local coordinate system, the thermal–mechanical model of the XR-cell in the global unit cell coordinate system (x–y–z) can be determined by a coordinate transformation. The thermal– mechanical constitutive model of the XR-cell in the global unit cell coordinate system can be expressed as 8 9 2 ½C 11 XR ½C 12 XR ½C 13 XR 0 0 rxx > > > > > > 6 ½C 12  ½C 22  ½C 23  >r > > > 0 0 yy > > 6 XR XR XR > > > 6 = 6 ½C 13 XR ½C 23 XR ½C 33 XR 0 0 zz ¼6 6 > > 0 0 0 ½C  ½C 44 XR 45 XR > 6 > rxy > > > > 6 > > > ryz > > 4 0 0 0 ½C 45 XR ½C 55 XR > > : > ; ½C 16 XR ½C 26 XR ½C 36 XR 0 0 rzx XR 8 9 8 t 9 rxx > exx > > > > > > > > > > > > > > > > > > rtyy > e yy > > > > > > > > > > > = zz zz  þ > > 0 > cxy > > > > > > > > > > > > > > > > > > > > > 0 c > > > yz > ; > > > > : : ; t rzx XR czx XR

3 ½C 16 XR ½C 26 XR 7 7 7 ½C 36 XR 7 7 0 7 7 7 0 5 ½C 66 XR

ð38Þ

J. Lua / Composites: Part A 38 (2007) 1019–1037

where the stress [r]XR and strain [e]XR tensor in Eq. (38) are T T computed by ½rXR ¼ ½T 1 ½r0 XR and ½eXR ¼ ½T 2 ½e0 XR , respectively, and [T1] and [T2] are the transformation matrices. Given Cx = cos(hx), Sx = sin(hx), and CSx = cos(hx)sin(hx), the coefficients of the stiffness matrix defined in Eq. (38) can be calculated by 4 11 ½C 11 XR ¼ C 4x ½C 11 nn XR þ S x ½C ss XR   33 þ 2CS 2x ½C 11 ns XR þ 2½C ss XR

½C 12 XR ½C 13 XR ½C 16 XR ½C 22 XR ½C 23 XR ½C 26 XR ½C 33 XR ½C 36 XR ½C 44 XR ½C 45 XR ½C 55 XR ½C 66 XR

2 21 ¼ C 2x ½C 12 nn XR þ S x ½C ns XR  11 2 ¼ ðC 4x þ S 4x Þ½C 11 ns XR þ CS x ½C nn XR  33 þ½C 11 ss XR  4½C ss XR  11 2 11 ¼ CS x C 2x ð½C 11 ns XR  ½C nn XR Þ þ S x ð½C ss XR  2 2 33 ½C 11 ns XR Þ þ 2ðC x  S x Þ½C ss XR ¼ ½C 22 nn XR 2 2 21 ¼ S x ½C 12 nn XR þ C x ½C ns XR  21  ¼ CS x ½C ns XR  ½C 12 nn XR 4 11 ¼ S 4x ½C 11 nn XR þ C x ½C ss XR

 33 þ 2CS 2x ½C 11 ns XR þ 2½C ss XR 



11 11 2 ¼ CS x S 2x ½C 11 ns XR  ½C nn XR þ C x ½C ss XR   2 2 33 ½C 11 ns XR  2ðC x  S x Þ½C ss XR 2 22 ¼ C 2x ½C 33 nn XR þ S x ½C ss XR

22  ¼ CS x ½C ss XR  ½C 33 nn XR 2 33 ¼ C 2x ½C 22 ss XR þ S x ½C nn XR   11 11 ¼ CS 2x ½C 11 nn XR þ ½C ss XR  2½C ns XR

2 þ C 2x  S 2x ½C 33 ss XR

ð39aÞ ð39bÞ ð39cÞ ð39dÞ ð39eÞ ð39fÞ ð39gÞ ð39hÞ ð39iÞ ð39jÞ ð39kÞ ð39lÞ ð39mÞ

In addition, the thermal-mismatch-induced residual stress components in Eq. (35) are given by

1027

To determine the strain response of each cell from the unit-cell response parameters, a parallel–series system is employed to establish the physical linkage of these four cells in forming the unit-cell model shown in Fig. 4. The same approach has been used by Tabiei and Yi [24] for a balanced weave to achieve the assumed stress and strain compatibility conditions between the unit cell and its associated 4 subcells. Using the analogy associated with a parallel system, it is assumed that (1) the volume average transverse normal strain (ezz) and in-plane shear strain (cxy) of each cell are the same as the corresponding volume average strain component of the unit cell; and (2) the volume average transverse shear stresses ryz and rzx of each cell are the same as the corresponding volume average stress of the unit cell. Given these two assumptions, we have ½ezz XY ¼ ½ezz XR ¼ ½ezz YR ¼ ½ezz RR ¼ ezz ½cxy XY ¼ ½cxy XR ¼ ½cxy YR ¼ ½cxy RR ¼ cxy

ð42Þ ð43Þ

yz ½ryz XY ¼ ½ryz XR ¼ ½ryz YR ¼ ½ryz RR ¼ r zx ½rzx XY ¼ ½rzx XR ¼ ½rzx YR ¼ ½rzx RR ¼ r

ð44Þ ð45Þ

For a given volume average strain of the unit cell, e ¼ f exx eyy ezz cxy cyz czx g, the total number of unknown strains in these four cells is 24. Using Eqs. (42) and (43), the number of the remaining unknown strains is 16. Adding the unknowns in Eqs. (44) and (45), we can conclude that the total number of unknowns is 18 given yz and r zx . by {exx, eyy, ezz, cxy, cyz, czx}IC=XX,XR,YR,RR, and r To solve these 18 unknowns, the 18 governing equations are summarized here. Using Eqs. (44), (45), (31), (38), (A.4) and (A.7), the first 8 equations can be written as yz ½C 22 ss XY fcyz gXY ¼ r

ð46Þ

½C 33 ss XY fczx gXY

zx ¼r

ð47Þ

yz ½C 45 XR fcxy g þ ½C 55 XR fcyz gXR ¼ r

ð48Þ

½C 16 XR fexx gXR þ ½C 26 XR feyy gXR þ ½C 36 XR fezz g zx þ ½C 66 XR fczx gXR þ ½rtzx XR ¼ r

ð49Þ

½C 15 YR fexx gYR þ ½C 25 YR feyy gYR þ ½C 35 YR fezz g yz þ ½C 55 YR fcyz gYR þ ½rtyz YR ¼ r

ð50Þ

5. Determination of four-cell response from unit-cell response parameters

zx ½C 46 YR fcxy g þ ½C 66 YR fczx gYR ¼ r yz ½C 44 RR fcyz gRR ¼ r

ð51Þ ð52Þ

Based on Eqs. (31), (38), (A.4) and (A.7) for the XY-, XR-, YR- and RR-cells respectively, the thermal–mechanical model of these four cells defined in the global unit-cell coordinate system can be summarized as

zx ½C 44 RR fczx gRR ¼ r

ð53Þ

0

0

½rtxx XR ¼ C 2x ½rtxx XR þ S 2x ½rtzz XR 0

½rtyy XR ¼ ½rtyy XR 0

ð40bÞ 0

½rtzz XR ¼ S 2x ½rtxx XR þ C 2x ½rtzz XR

0  0 0 ½rtzx XR ¼ CS x ½rtzz XR  ½rtxx XR

½rIC ¼ ½CIC ½eIC þ ½rt IC ;

ð40aÞ

IC ¼ XY ; XR; YR; RR

ð40cÞ ð40dÞ

Using the stress and strain compatibility conditions associated with the parallel–series system of the four-cell model (see Fig. 4), the remaining equations are given by

ð41Þ

where [r]IC is the column vector of the stress; [e]IC is strain; [rt]IC is thermal stress of the IC-th cell; and [C]IC is the stiffness matrix of the IC-th cell. To compute the cell stress using Eq. (41), the cell strain has to be determined first from a given volume average strain ½e and temperature increase (DT) of the unit cell.

½rxx YR ¼ ½rxx RR ½rxx XY ¼ ½rxx XR

ð54Þ ð55Þ

½ryy XY ¼ ½ryy YR ½ryy XR ¼ ½ryy RR

ð56Þ ð57Þ

Lf fexx gXY þ ðLx  Lf Þfexx gXR ¼ exx Lx

ð58Þ

1028

J. Lua / Composites: Part A 38 (2007) 1019–1037

Lf fexx gYR þ ðLx  Lf Þfexx gRR ¼ exx Lx

ð59Þ

Lf feyy gXY þ ð1  Lf Þfeyy gYR ¼ eyy

ð60Þ

Lf feyy gXR þ ð1  Lf Þfeyy gRR ¼ eyy

ð61Þ

Lf ð1  Lf Þfcyz gYR þ Lf ðLx  Lf Þfcyz gXR þ ð1  Lf ÞðLx  Lf Þfcyz gRR þ L2f fcyz gXY ¼ cyz Lx

ð62Þ

Lf ð1  Lf Þfczx gYR þ Lf ðLx  Lf Þfczx gXR þ ð1  Lf ÞðLx  Lf Þfczx gRR þ L2f fczx gXY ¼ czx Lx

ð63Þ

Using Eqs. (31), (38), (A.4) and (A.7), an explicit form of Eqs. (54)–(57) can be derived in terms of strain and thermal stress components of the relevant cells. The Lf and Lx parameters used in Eqs. (58)–(63) are defined in Section 3. After a few manipulations, these 18 linear system equations can be reduced to a system of 8 equations in terms of 8 independent variables, namely, {exx}XR, {exx}YR, {eyy}XY, {eyy}RR, {cyz}YR, {czx}YR, {cyz}XR, and {czx}XR. Given feg and DT of the unit cell, the six strain components of each cell and the transverse shear stress of the unit cell, zx g can be obtained from solving the linear system f ryz ; r equations. The cell stress components are then determined from the cell strain using Eqs. (31), (38), (A.4) and (A.7). 6. Numerical implementation of dual-function micromechanics model A micromechanics analysis tool, CELLMAT, was developed based on the four-cell representation of an unbalanced plain weave (UPW) unit-cell model. For a given volume average strain ðeÞ and the temperature increase (DT) of the UPW, dual analysis modules have been implemented in CELLMAT for (1) computing the effective thermal–mechanical properties of the UPW from its constituent properties (woven fabric analyzer) and (2) determining the constituent stress and strain in each cell from given e and DT (stress/strain decomposition module). The determined constituent stress and strain can be used in a mechanism-driven failure criterion for a rational stiffness degradation prediction during a progressive failure analysis. The numerical procedures associated with these two solution modules in CELLMAT are summarized in Sections 6.1 and 6.2. 6.1. Woven fabric analyzer (WFA) As shown in Section 5, both the volume average stress and strain of these four cells along with the transverse zx g of UPW can be determined for a shear stress f ryz ; r given feg and DT of UPW. With the aid of cell stress ([r]IC, IC = XY, XR, YR, RR), the remaining stress compoyy ; r zz ; r xy g can be determined from the nents of UPW f rxx ; r cell compatibility condition as shown here. 1h xx ¼ Lf frxx gXY þ ð1  Lf Þfrxx gYR r 2 i þ Lf frxx gXR þ ð1  Lf Þfrxx gRR ð64Þ

1 h Lf fryy gXY þ ðLx  Lf Þfryy gXR 2Lx i þ Lf fryy gYR þ ðLx  Lf Þfryy gRR 1 h zz ¼ Lf ð1  Lf Þfrzz gYR þ Lf ðLx  Lf Þfrzz gXR r Lx i þ ð1  Lf ÞðLx  Lf Þfrzz gRR þ L2f frzz gXY 1 h xy ¼ Lf ð1  Lf Þfrxy gYR þ Lf ðLx  Lf Þfrxy gXR r Lx i þ ð1  Lf ÞðLx  Lf Þfrxy gRR þ L2f frxy gXY yy ¼ r

ð65Þ

ð66Þ

ð67Þ

Thus, a mapping relation between the calculated volume  ¼ f yy ; r zz ; r xy ; r yz ; r zx g and a given pair average stress r rxx ; r of fe; DT g is established for the UPW unit cell. This mapping relation paves the way to compute the effective thermal–mechanical properties of the UPW from its constituent and microstructural parameters. To be consistent with a user-defined material model associated with a commercial finite element solver, a prescribed strain is selected as our load case. Since the behavior of UPW can be characterized as an effective orthotropic material, nine stiffness coefficients (C11, C12, C13, C22, C23, C33, C44, C55, C66) and three coefficients of thermal expansion (a1, a2, a3) need to be determined using nine mechanical load cases and one thermal load case. The nine stiffness coefficients are determined first by applying the energy equivalence of the UPW unit cell subi jected to nine independent strain fields, feg , i = 1, 2, . . . , 9. i Given the ith volume average strain feg ¼ feixx ; eiyy ; eizz ; i i i cxy ; cyz ; czx g and DT = 0, the resulting volume average stress iyy ; r izz ; r ixy ; r iyz ; r izx g of the UPW can be deterf rgi ¼ f rixx ; r mined from the developed mapping relation. The resulting energy equivalence states 2

2

2

C 11eixx þ C 22eiyy þ C 33eizz þ 2C 12ei11ei22 þ 2C 13eixxeizz 2

2

2

i iT þ 2C 23eiyyeizz þ C 44 ci12 þ C 66 ci13 þ C 55 ci23 ¼ feg f rg

ð68Þ The left side of Eq. (68) is the strain energy computed from a smeared orthotropic constitutive model for the UPW while the right side represents the strain energy determined from the four-cell homogenization model of the UPW. To determine the nine stiffness coefficients (Cij) in Eq. (68), nine independent load cases are used with the corresponding strain field defined by 3 8 19 2 exx 0 0 0 0 0 e > > > > > 6 > > > > 0 0 0 7 e2 > 7 6 0 eyy 0 > > > > 7 6 > > > 3> 6 > > 0 ezz 0 0 0 7 > 7 > 6 0 > e > > > > 7 4> > 6 > 0 0 0 c12 0 0 7 > = 6 < e > 7 6 9 ð69Þ ½e ¼ e5 ¼ 6 0 0 0 0 c23 0 7 7 6 > 6 > > 7 6> > > 0 0 0 0 0 c31 7 > > 6 > e > > > 7 6 > > 7 > > 6 > 6 exx eyy 0 > e > 0 0 0 7 > 7 > > > 7 6 > > 8> > >   0 e 0 0 0 e  5 4 e yy zz > > > ; : 9> exx 0 ezz 0 0 0 e

J. Lua / Composites: Part A 38 (2007) 1019–1037

After solving Eq. (68) for the stiffness matrix [C] using the nine load cases given in Eq. (69), the nine engineering elastic constants (E1, E2, E3, m21, m31, m32, G12, G23, G31) can be computed from the compliance matrix [S], which is determined from the inversion of [C]. Given the determined [C], the three equivalent coefficients of thermal expansion {a1, a2, a3} are calculated from the thermal– mechanical model of the UPW by 2 38 9 8 t 9 C 11 C 12 C 13 > = > = < a1 > < rxx > 6 7 rtyy ¼ ð70Þ C C C a 4 12 22 23 5 2 > ; > ; : > : t > rzz C 13 C 23 C 33 a3 where rtxx , rtyy , and rtzz are the components of the residual thermal stress of UPW associated with a unit temperature increase (DT = 1.0) subjected to a zero-displacement constraint ðexx ¼ eyy ¼ ezz ¼ 0:0Þ in the UPW unit cell. 6.2. Stress/strain decomposition module The stress/strain decomposition module in CELLMAT is developed based on the multicell homogenization at the unit cell (UPW) and the four-cell level. The mapping relation between the UPW cell response and the response associated with each of its four cells (XY, XR, YR, RR) is established in Section 5. After determining the cell stress/strain response parameters from fe; DT g of UPW (level 1 stress/strain decomposition), the mapping relation between the cell response to the response in each of its constituents is established with the level 2 stress/strain decomposition. To illustrate this two-level decomposition process, the XY-cell of UPW is considered here. After solving the linear system equations given in Section 5, the cell strain components {e}XY are determined from a given fe; DT g of UPW. Using Eq. (31), the cell stress {r}XY can be obtained. Given {e}XY and {r}XY and the definition of iso-strain and isostress (Eqs. (2b) and (3b)), the following stress/strain components in X-tow and Y-tow of the XY-cell are given by f rzz f exx

ryz eyy

rzx gX ¼ f rzz ¼ f rzz

ryz

rzx gY

ryz

rzx gXY ¼ frs gXY

cxy gX ¼ f exx

eyy

cxy gY

¼ f exx

eyy

cxy gXY ¼ fen gXY

ð71Þ

Using Eqs. (10) and (11), the remaining unknown constituent strain fes gX ;Y ¼ f ezz cyz czx gX ;Y and constituent stress frn gX ;Y ¼ f rxx ryy rxy gX ;Y in X- and Y-tow are computed by 1

1029

where k = X, Y and all the C-matrices {[Cij]k} are computed from Eqs. (18)–(21). The stress-free thermal strains T in Eqs. (72) and (73) are determined by f etn ets gk ¼ T ½ak DT. By applying this two-level decomposition for all four cells, the stress and strain can be quantified for each constituent within each cell. To retain physical identity of an individual phase and be consistent with the mechanism-driven failure criteria, three distinct constituent phases are selected as the multiphase representation of the plain weave woven fabric composite. They are: (1) X-tow; (2) Y-tow; and (3) resin (matrix). To compute the average r/e an Xtow, a volume average is performed for r/e in the X-tow of the XY- and XR-cells. Similarly, the r/e in a Y-tow is computed from the volume average of r/e in the Y-tow of the XY- and YR-cells. Finally, the stress and strain in the resin is determined by performing the volume average of r/e of R in the XR-, YR-, and RR-cells. The determined volume average stress and strain in X-tow, Y-tow, and resin can be used in the mechanism-driven failure criteria to predict the failure mode and the resulting stiffness degradation during a progressive damage analysis. Thus, the two-level stress/strain decomposition module in CELLMAT provides additional windows to explore the multiscale damage and the synergistic interaction toward the final rupture of a structural element. 6.3. Remarks The four-cell representation of an unbalanced plain weave unit cell provides a hierarchical micromechanics framework for quantification of constituent damage induced thermal–mechanical properties and stress/strain response at multiscale levels. The WFA in CELLMAT serves as a material virtual testing tool to determine the effective thermal–mechanical properties at the composite level for given microstructural and constituent material parameters. The two-level stress/strain decomposition in CELLMAT establishes a link among the response parameters at different microstructural levels. In comparison with the FEM-based micromechanics model (TMAT), CELLMAT requires only eight linear system equations for solving a mechanical load case and three linear system equations for a thermal load case. Because of the high computational efficiency, CELLMAT can be used as a numerical material model associated with a finite-element-based failure prediction where the degraded thermal–mechanical properties have to be determined at each Gaussian point and at each time step.

1

fes gk ¼ ½C ss k frs gXY  ½C ss k ½C sn k fen gXY t t þ ½C ss 1 k ½C sn k fen gk þ fes gk

1 frn gk ¼ ½C nn k  ½C ns k ½C ss k ½C sn k fen gXY

7. Example applications ð72Þ

1

þ ½C ns k ½C ss k frs gXY

t þ ½C ns k ½C ss 1 k ½C sn k  ½C nn k fen gk

ð73Þ

To demonstrate the accuracy and numerical efficiency of the dual function of the CELLMAT solver, three numerical examples are given here. The application of CELLMAT for thermal–mechanical material quantifications of balanced and unbalanced plain weave woven fabric laminates

1030

J. Lua / Composites: Part A 38 (2007) 1019–1037

is summarized in Sections 7.1 and 7.2, respectively, along with comparisons to available test and analysis results. To illustrate the constituent failure induced stiffness degradation, the CELLMAT solver is applied to a tensile specimen using one element model where the constituent stress components are used in the maximum stress criteria. Both the failure sequence and damage mode are predicted from the constituent-based failure model.

constituent properties of the isotropic matrix (epoxy) and the transversely isotropic tow are summarized here. 1. Epoxy matrix: E ¼ 3:5 GPa; m ¼ 0:35; 2. E-glass/epoxy tow Exx ¼ 47:77 GPa; Eyy ¼ Ezz ¼ 18:02 GPa; Gxy ¼ Gxz ¼ 5:494 GPa; Gyz ¼ 3:877 GPa;

7.1. Balanced E-glass/epoxy plain weave

mxy ¼ 0:314;

The balanced plain weave used in Ref. [24] consists of 35% E-glass fiber with the fiber volume fraction within the tow of 65%. Given these two volume fractions, the tow volume fraction (Vtow) within the unit cell is 0.5385, which is calculated by the ratio of the overall fiber volume fraction (35%) to the fiber volume fraction in the tow (65%). Using Eq. (24) with Lx of 1.0 for the balanced weave, we have Lf ¼ V tow ¼ 0:5385. Assuming the same tow cross-sectional geometry as shown in Fig. 2, the spacing between either the X-tow or the Y-tow is 4.1275 mm to achieve Vtow of 0.5385. Based on the SEM examination of the tow architecture (see Fig. 2), the width and height of the tow is 5.461 mm and 0.4826 mm, respectively. Using these woven fabric microstructural parameters, the tow angle computed from Eqs. (25) and (26) is hx = hy = 7.22. While the orthotropic thermal–mechanical tow properties can be generated in CELLMAT using its triangular packing micromechanics model, we assumed that the mechanical properties of a tow are the same as those defined in Refs. [24,25]. The use of the same constituent properties can assist us in judging the accuracy of the CELLMAT tool by comparing its predictions directly with the results shown in Refs. [24,25]. As the coefficients of thermal expansion (CTE) of the tow (atow) are not given in either Ref. [24] or Ref. [25], they are generated from the CELLMATs triangular packing micromechanics model using (1) the CTE of the E-glass fiber of 5.4 le/K; (2) the CTE of the epoxy resin of 65 le/K; and (3) the fiber volume fraction within the tow of 0.65. The thermal mechanical

a

am ¼ 65 le=K

myz ¼ 0:249;

mzx ¼ 0:118;

ax ¼ 7:135 le=K; ay ¼ az ¼ 28:894 le=K The thermal–mechanical properties predicted from CELLMAT are summarized in Table 1. To demonstrate the validity and accuracy of the CELLMAT prediction, the prediction from the FEM-based micromechanical model [25] and Tabiei’s cell model [24] are also given in Table 1. In addition, the TMAT tool [16] is also applied for the present balanced plain weave to demonstrate the computational efficiency associated with CELLMAT. In comparison with CELLMAT where only four cells are used in the microstructural representation of the plain weave, the FEM-based micromechanics model generated from TMAT has 6805 nodes and 7888 solid elements. Both the deformed shape and the transverse shear stress distribuTable 1 Comparison of CELLMAT predictions of a balanced weave with TMAT and other published data Material properties

CELLMAT (four-cell model)

TMAT (FEM-model)

Ref. [21] (FEMmodel)

Tabiei’s four-cell model [20]

Exx/Eyy (GPa) Ezz (GPa) Gxy (GPa) Gyz/Gxz (GPa) mxy myz/mxz axx/ayy (le/K) azz (le/K)

17.935 9.820 3.542 2.323 0.172 0.328 22.352 50.868

18.552 9.403 3.009 2.372 0.184 0.334 21.06 54.690

18.634 8.346 3.190 2.422 0.175 0.372 – –

17.853 9.788 3.529 2.497 0.172 0.332 – –

b

V1 L7 C1 G4

1.644 1.536 1.428

V1 L7 C1

1.522 1.422

1.32

1.322

1.212

1.222

1.104

1.122

0.996

1.022

0.888

0.922

0.78

0.821

0.672

0.721

0.564

0.621

0.456

0.521

0.348

0.421 0.24

Z

0.321 0.132

X

Y

Output Set: MSC/NASTRAN Case4 Deformed(18.3): Total Translation Contour: SolidYZ Shear Stress

Z

0.221

0.0237

0.121 -0.0843

Y X Output Set: MSC/NASTRAN Case4 Deformed(18.3): Total Translation Contour: Solid YZ Shear Stress

0.0204 -0.0797

Fig. 5. TMAT prediction of the transverse shear stress (ryz) distribution subjected to a transverse shear load case: (a) in fabric; (b) in unit cell.

J. Lua / Composites: Part A 38 (2007) 1019–1037

tion (ryz) resulting from the transverse shear load case are shown in Fig. 5. As shown in Table 1, the mechanical properties predicted from CELLMAT are in good agreement with both the FEM-based tool and Tabiei’s four-cell model. Different from Tabiei’s four-cell model, CELLMAT also computes the orthotropic coefficient of thermal expansions (CTEs) at the composite level for a given CTE at the constituent level. With this additional capability, both the processingand environmental-induced thermal residual stress can be included during the progressive failure analysis. As can be seen from Table 1, the CTEs predicted from CELLMAT are in good agreement with the results from the finite-element-based micromechanics solver, TMAT. The small discrepancy in thermal mechanical property prediction from using the CELLMAT and TMAT tool is mainly attributed to the geometric approximation of woven fabric architecture and the assumed compatibility conditions used in CELLMAT. 7.2. Unbalanced E-glass/vinyl ester plain weave An unbalanced plain weave fabricated by Seemann Composites, Inc., is selected here (see Fig. 2). The thermal–mechanical properties of its constituents (E-glass fiber, vinyl ester resin) are summarized in Table 2. The X-tows are spaced 6.4 mm (Ly = 6.4 mm), while the Y-tows are spaced 5.1 mm (Lx = 5.1 mm). The fiber volume fraction of the unbalanced E-glass/vinyl ester plain weave is 51%, which is much higher than the balanced weave example in Section 7.1. Fiber volume fraction within the tow and tow volume fraction within the unit cell are determined from the microscopic examination (see Fig. 2) and are given by 75% and 68%, respectively. Based on the SEM picture shown in Fig. 2, the width of X/Y tows is 5.461 mm, and the tow height is 0.483 mm. By substituting these geometrical parameters in Eqs. (24)–(26), the Lf parameter and the tow angles (hx, hy) are given by 0.605, 11.092, and 5.446, respectively. Numerous tensile and in-plane shear coupon tests have been performed by Lopez-Anido et al. [26] for the unbalanced plain weave based on ASTM D3039 and ASTM D4255. As reported by Fayad [27], the in-plane elastic properties along with the statistical dispersion are shown in Table 3. The thermal–mechanical properties predicted from CELLMAT are summarized in Table 4. To demon-

Table 2 Summary of constituent properties of Seemann composite (plain weave) Thermal/mechanical properties

E-glass (fiber)

Vinyl ester (matrix)

Elastic modulus E (GPa) Poisson’s ratio m Coefficient of thermal expansion j (le/K)

72.45 0.225 5.4

3.25 0.325 65

1031

Table 3 Summary of test data on in-plane elastic properties along with their statistical parameters In-plane elastic properties

Mean

Coefficient of variation (%)

Standard deviation

Test method

Exx (GPa) Eyy (GPa) m12 G12 (GPa)

22.61 26.57 0.136 4.857

1.64 1.11 8.78 3.5

0.371 0.295 0.012 0.170

ASTM ASTM ASTM ASTM

D3039 D3039 D3039 D4255

Table 4 Comparison of CELLMAT results with TMAT prediction and test data Thermal–mechanical properties

CELLMAT prediction

TMAT prediction

UMaine test data [22]

Exx (GPa) Eyy (GPa) Ezz (GPa) Gxy (GPa) Gyz (GPa) Gzx (GPa) mxy myz mxz axx (le/K) ayy (le/K) azz (le/K)

22.185 26.160 11.909 5.041 2.879 2.877 0.132 0.309 0.309 16.423 13.806 35.766

24.821 24.821 10.811 4.495 3.371 3.371 0.138 0.316 0.316 13.986 13.986 40.917

22.61 (ASTM D3039) 26.57 (ASTM D3039) – 4.857 (ASTM D4255) – – 0.136 – – – – –

strate the validity and accuracy of the CELLMAT prediction, the prediction from the FEM-based micromechanics model (TMAT) and the mean value shown in Table 3 are also given in Table 4. Because of the use of the balanced weave model in TMAT, the elastic constants along the Xand Y-tow are identical in the TMAT prediction. The inplane material properties measured from UMaine’s test are in excellent agreement with the CELLMAT predictions. The marginal difference between the CELLMAT and TMAT predictions are mainly due to the balanced weave assumption used in TMAT and the simplified geometric and structural model used in CELLMAT. The Halpin-Tsai equation [28] has been considered the most popular empirical model in characterizing the effective properties of a composite system by taking into account the geometry of its constituents. To assess the validity and accuracy of using the Halpin-Tsai model for quantification of mechanical properties of a plain weave woven fabric composite, the current E-glass/vinyl ester plain weave system is considered here. Because of the interlacing of X- and Y-tows (see Fig. 2), the Halpin-Tsai model cannot be directly applied to the plain weave ply. To remedy this difficulty, the plain weave ply is divided into two sub-plies where each ply consists of tows in one direction. With this unidirectional ply representation of a plain weave, the Halpin-Tsai model is applied first to determine the ply elastic constants followed by the lamination theory for quantifying the smeared in-plane elastic constants of the plain weave.

1032

J. Lua / Composites: Part A 38 (2007) 1019–1037

In these two sub-plies, it is assumed that the X sub-ply has the tow direction along the x-axis, while the Y subply has the tow along the y-axis. The overall fiber volume fraction in the X- and Y-sub-ply is 51%, and the fiber and matrix elastic properties are the same as those given in Table 2. Based on the Halpin-Tsai model for the continuous fiber configuration, the conventional rule of mixtures is used for the axial modulus E11 and the Poisson’s ration m12. The remaining elastic properties of the sub-ply are determined by pc =pm ¼ ð1 þ fgV f Þ=ð1  gV f Þ

ð74Þ

where g ¼ ðpf =pm  1Þ=ðpf =pm þ fÞ

ð75Þ

As shown in Eq. (74), a new empirical parameter f is introduced which has no scientific basis or relation to any material or geometry property. In Eqs. (74) and (75), the subscripts c, m, and f denote the parameters associated with the composite, matrix, and fiber, respectively. Eqs. (74) and (75) are used to determine the ply properties of G12, E22, G23, and m23, where the empirical parameter f is calculated for each of these properties based on the continuous fiber model given by Halpin [28]. The resulting subply material properties calculated from the rule of mixture and Eqs. (74) and (75) are given as E11 ¼ 38:542 GPa;

E22 ¼ E33 ¼ 11:228 GPa;

G12 ¼ G13 ¼ 3:440 GPa;

G23 ¼ 2:746 GPa; m12 ¼ 0:274 ð76Þ

Because m23 as determined from Halpin’s model (m12 = 1  E2/G23) is a negative number, it is assumed that m23 = m12. Using the ply-properties given by Eq. (76) and the classical lamination theory, the in-plane smeared properties of the plain weave can be calculated and given by Exx ¼ Eyy ¼ 25:052 GPa; Gxy ¼ 3:439 GPa

mxy ¼ 0:1236; ð77Þ

Comparing the above model predictions with the data shown in Table 4, both the elastic modulus (Exx/Eyy) and

the Poisson’s ratio (mxy) are in good agreement with the TMAT prediction. However, the matrix dominant shear modulus Gxy has about 30% error as compared with the experimental data. While a new curve-fit constant f12 can be re-derived for use in Eqs. (74) and (75) to reduce the model prediction error in Gxy, this process would not be feasible during a progressive failure analysis as the test data of Gxy at any arbitrary constituent damage state are not available. Therefore, for a woven fabric composite, a reliable micromechanics model that captures the dominant feature of the woven fabric architecture has to be used to achieve sufficient accuracy and computational efficiency. To demonstrate the stress/strain decomposition capability in CELLMAT, consider the same unbalanced plain weave subjected to an in-plane shear strain c12 of 0.002. A two-level decomposition is performed to establish (1) the mapping relation of the global composite response to the four-cell response and (2) the mapping relation of each four-cell response to the response in each of its constituents. Table 5 lists the constituent stress components of each cell defined in Fig. 4. These stress components are given in the local material coordinate system for each constituent with x1 along the tow direction. For the given c12 of 0.002, the effective composite shear stress r12 can be calculated by r12 = G12 * c12. Using G12 of 5.04 GPa given in Table 4, the average in-plane shear stress at the composite level is 10.08 MPa. As can be seen from Table 5, the volume average shear stress in either the X-tow or the Y-tow is about 13 MPa, which is 30% higher than the corresponding average composite stress. In addition, the shear stress in the resin (matrix) constituent is about 2.5 MPa, which is significantly lower than the tow shear stress. Given the in-plane shear failure strength of the tow and in-situ resin of 60.8 MPa and 34.4 MPa, respectively [29], the shear failure will initiate in a tow first. The constituent damage induced stiffness degradation and the rotation of X and Y tows will result in a nonlinear in-plane shear response under subsequent shear straining. Despite the pure in-plane shear defined at the composite level, small axial and transverse shear stresses are introduced in each constituent as

Table 5 Summary of local stress components in each constituent of a cell after the decomposition Subcell/constituent

Constituent stress in local material coordinate system (MPa) r11

r22

r33

s12

s23

s31

XY – subcell X-tow Y-tow

0.124 0.102

0.0451 0.0509

0.0244 0.0244

13.722 13.722

0.144 0.153

0.153 0.144

XR – subcell X-tow Resin

0.295 0.0607

0.0138 0.0247

0.0703 0.0222

13.104 2.560

1.361 0.887

0.122 0.124

YR – subcell Y-tow Resin

0.203 0.0123

0.00589 0.00286

0.0353 0.00978

13.592 2.476

0.608 0.134

0.134 0.376

0.00544

0.00281

2.453

RR – subcell Resin

0.00320

0.1442

0.153

J. Lua / Composites: Part A 38 (2007) 1019–1037

shown in Table 5 due to the material heterogeneity. These constituent stresses will play a key role in driving the mechanism-based constituent failure mode and the resulting stiffness degradation at the composite level. Because of its high computational efficiency and versatility, CELLMAT is well suited for integration with a progressive failure analysis module of a large composite structure. 7.3. Failure mode and stiffness degradation prediction with a one-element model

Table 6 Summary of strength parameters of tow (X/Y) and resin (R) YT(ZT)

YC(ZC)

TXY

TXZ(TYZ)

Tow (unidirectional composite) X/Y constituent 817.5 759.7

45.26

144.30

60.8

48.52

Resin pocket (in-situ property) R constituent 37.1 118.1

37.1

118.1

34.4

34.4

Strengths (MPa)

ri ¼

A typical tensile stress–strain curve obtained from a tensile specimen in Ref. [26] is shown in Fig. 6. A linear curvefit is performed to determine the elastic modulus at the initial virgin and the final damage state. The elastic modulus determined from the curve-fit of experimental data is 26.35 MPa and 20.21 MPa at undamaged and final damage state, respectively. The goodness of fit measured by (R2) is also given in Fig. 6 with R2 of 1 for the perfect fit. Conventional progressive analyses are unable to quantify the correlation between the observed stiffness degradation in the composite and the degree of damage in each of its constituents. To bridge the global composite response and the damage evolution in each constituent, a progressive failure analysis is performed on the tensile specimen using CELLMAT coupled with the maximum constituent stress failure criteria. For a given load increment, a stress–strain decomposition is performed using CELLMAT. Based on the volume average stress in the X/Y tow and resin pocket, the maximum stress criteria are applied. As shown in Table 6, the strength parameters given by Mayers [29] are used for the tow and resin constituent. For a damaged constituent, a continuous damage function [30] is used to compute the stiffness degradation in each damaged constituent. Given the maximum ith stress (ri) and it associated failure strength (Si) of a constituent, a damage threshold based on the maximum stress criterion for the ith failure mode can be computed by

1033

XC

XT

 ri Si

ð78Þ

Given the damage thresholds ri, associated with the ith failure mode, the corresponding damage variable xi can be determined by 1

m

xi ¼ 1  emeð1ri Þ

ðno sum on index iÞ

where m is a strain rate softening constant. Using Eq. (79), the constituent induced stiffness reduction can be characterized by Ci = (1  xi)Ci0, where Ci0 is the constituent stiffness at its virgin state. After computing the degraded stiffness for each constituent, the CELLMAT solver is used again to determine the effective composite stiffness from the damaged constituent properties. By applying a uniaxial tensile load in the Y-tow direction, the constituent stress components at some critical load stages are summarized in Table 7. The bold numbers in Table 7 indicate the stress associated with the initial constituent failure. As shown in Table 7, resin matrix cracking initiates at e1 = 7.4069e03, which is near the transition point of a change in slope of the stress and strain curve (see Fig. 6). At e1 = 1.3565e02, the X-tow experiences transverse tensile failure, where the composite stiffness reduces from 26.127 GPa to 23.81 GPa. The elastic modulus at e1 = 1.6565e02 is 20.155 GPa, which is consistent with the final modulus shown in Fig. 6 before the brittle rupture. At r1 = 461.29 MPa, the ultimate longitudinal tow strength in Y-tow has been reached and catastrophic failure follows. The predicted ultimate tensile strength of 461.29 MPa is in

500 450 400

Final Chord Modulus Range

Initial Chord Modulus Range

Stress (MPa)

350 300

y = 20.218x + 42.524

250

R2 = 0.9999

200 150 100 50 0 0.000000

y = 26.351x + 1.1639 R2 = 0.9998 0.005000

0.010000

ð79Þ

0.015000

0.020000

Strain (mm/mm) Fig. 6. A tensile stress–strain curve based on an ASTM D3039 test.

0.025000

1034

J. Lua / Composites: Part A 38 (2007) 1019–1037

Table 7 Summary of local stress components in each constituent of a cell at different failure stages Subcell/constituent

Constituent stress in local material coordinate system (MPa) s12 (MPa)

s23 (MPa)

s31 (MPa)

At initial matrix failure: e1 = 7.4609e03, r1 = 195.04 MPa, E = 26.127 GPa Y-tow 375.26 12.035 0.6286 X-tow 28.544 106.01 5.4664 Resin 37.107 11.338 7.2850

r11 (MPa)

r22 (MPa)

r33 (MPa)

0.08409 0.8033 0.1079

0.5526 6.4869 0.050

3.5931 0.220 1.8063

At e1 = 1.0565e02, r1 = 273.744 MPa, E = 25.321 GPa Y-tow 525.921 16.927 X-tow 39.202 146.99 Resin 48.901 15.318

0.01242 1.2951 0.15986

0.7457 9.3623 0.06024

4.7361 0.2953 2.3846

At initial Y-tow transverse failure: e1 = 1.3565e02, r1 = 345.27 MPa, E = 23.81 GPa Y-tow 656.27 20.434 0.7909 X-tow 48.836 179.59 9.7825 Resin 55.375 17.882 12.025

0.07358 2.200 0.21426

0.9095 12.179 0.06703

5.3652 0.3486 2.6835

At e1 = 1.6565e02, r1 = 405.8125 MPa, E = 20.155 GPa Y-tow 744.386 21.534 X-tow 56.885 196.368 Resin 58.111 18.998

1.8169 13.083 12.941

0.22494 4.2627 0.25352

1.0633 14.229 0.07089

5.5077 0.3921 2.7523

At final X-tow fiber failure: e1 = 1.9565e02, r1 = 461.29 MPa, E = 18.49 GPa Y-tow 816.776 22.063 4.9506 X-tow 63.997 209.27 16.552 Resin 60.123 19.833 13.652

0.40772 6.7181 0.28520

1.2061 15.932 0.07367

5.5572 0.4345 2.7921

1.1140 7.4581 10.090

Table 8 Comparison of composite tensile strain predictions at four load levels

good agreement with the ultimate strength (450 MPa) shown in Fig. 6. To further demonstrate the validity of the model prediction, the predicted composite tensile strains (e1) at these five stress levels are compared with the corresponding experimental data given in Fig. 6. As can be seen from Table 8, the strain predicted from the constituent-based progressive failure model agrees well with the test data.

and unbalanced weave with available numerical results and test data. Despite the use of the simple maximum stress criteria, both the failure events and the associated stiffness degradation agree well with the experimental observation from a tensile coupon test. Thus, the integration of the dual-function micromechanics model with a constituentbased failure criterion provides an additional window to reveal the failure mode and failure sequence during a progressive failure analysis. The developed four-cell micromechanics model will be integrated with a general purpose laminator to form a stress/strain translator among different length scales (fiber/matrix/tow/ply/laminated plate). The physical mapping established among these different length scales will provide a rational way to explore an optimal material and processing procedure for a given design requirement of a composite ship structure.

8. Conclusions

Acknowledgements

A dual-function micromechanics model is developed for an unbalanced plain weave subjected to a thermal–mechanical load. The developed micromechanics model can not only characterize the effective thermal–mechanical properties of the unbalanced weave for given constituent damage, but can also compute the stress and strain at each constituent. The calculated constituent stress and strain can be used in a mechanism-driven failure criterion to predict the failure mode, failure sequence, and the synergistic interaction that leads to global stiffness degradation and the final rupture. Accuracy and computational efficiency are demonstrated by comparing the predictions of a balanced

The support of Dr. Luise Couchman and Dr. Roshdy Barsoum at the Office of Naval Research is gratefully acknowledged.

Stress level (MPa)

Model prediction

UMaine test data

Relative error (%)

195.04 273.74 345.27 405.81 450 (failure stress)

7.4609e03 0.010565 0.013565 0.016565 0.019565

8.0769e03 0.01115 0.01481 0.01788 0.0202

7.6 5.2 8.4 7.3 3.1

Appendix A A.1. Homogenization of YR-cell As shown in Fig. 4, the YR-cell is filled completely with the inclined part of the Y-tow of orientation angle of hy and the isotropic resin (R). The volume fraction of the resin pocket and the Y-tow in the YR-cell can be shown to be the

J. Lua / Composites: Part A 38 (2007) 1019–1037

same. Similar to the XR-cell described in Section 4.2, the local cell coordinate system of the YR-cell is obtained from rotation through a hy angle about the positive x-axis. Given the local cell coordinate system, the orthotropic thermal– mechanical model (Eq. (29)) is used for the Y-tow, and the isotropic model (Eq. (34)) is applied for the resin of the YR-cell. The resulting thermal–mechanical model of the YR-cell in its local coordinate system is given by 3 2 11 8 90 11 rxx > 0 0 0 ½C nn YR ½C 12 > nn YR ½C ns YR > > > 6 ½C 12  ½C 22  ½C 21  > > ryy > 0 0 0 7 > > 7 6 nn YR nn YR ns YR > = 7 6 11 21 11 0 0 0 7 6 ½C ns YR ½C ns YR ½C ss YR zz ¼6 7 33 > rxy > 6 0 0 0 ½C nn YR 0 0 7 > > > > 7 6 > > 22 > > 4 0 0 5 0 0 0 ½C ss YR > ; : ryz > 33 rzx YR 0 0 0 0 0 ½C ss YR 8 t 90 8 90 exx > rxx > > > > > > > > > > t > > > > e yy > > > > > > > > ryy > = < < t = ezz rzz  þ ðA:1Þ > > > 0 > > cxy > > > > > > > > > > > > > > > > > ; ; : 0 > : cyz > 0 YR czx YR where the coefficients of the stiffness matrix defined in Eq. (A.1) are given by " # y m 2 1 ðC  C Þ 12 C m11 þ C y22  23y ½C 11 ðA:2aÞ nn YR ¼ 2 C 22 þ C m11

   ðC m12  C y23 Þ C y12  C m12 1 y m C ½C 12  ¼ þ C þ ðA:2bÞ 12 nn YR 12 2 C m11 þ C y22 " # 2 1 ðC y  C m12 Þ C y11 þ C m11  12y ½C 22 ðA:2cÞ m nn YR ¼ 2 C 22 þ C 11 1 y m ðA:2dÞ ½C 33 nn YR ¼ ½C 55 þ C 44  2 y y m m C 23 C 11 þ C 12 C 22 ½C 11 ðA:2eÞ ns YR ¼ C y22 þ C m11 y y m m C 12 C 11 þ C 12 C 22 ½C 21 ðA:2fÞ ns YR ¼ C m11 þ C y22 2C y22 C m11 ½C 11 ðA:2gÞ ss XR ¼ C m11 þ C y22 y m 2C 55 C 44 ½C 22 ðA:2hÞ ss YR ¼ C m44 þ C y55 2C y44 C m44 ½C 33 ðA:2iÞ ss YR ¼ C m44 þ C y44 In Eq. (A.1), the thermal residual stress components resulting from a temperature increase, DT, in the local coordinate system are calculated by (  DT C y23 C m11 þ C m12 C y22 C y12 y C y23 y t 0 ½rxx YR ¼  a þ a 2 C m11 þ C y22 C y22 1 C y22 2  m   2C 12 C y12 C y23 y y y m þ þ 1 a1 þ a3 þ C 12  a1 C m11 C y22 ! ! ) 2 2 C y23 y 2C m12 m y m m þ C 22  y a2 þ C 11 þ C 12  m a1 C 22 C 11 ðA:3aÞ

0 ½rtyy YR

½rtzz 0YR

1035

(  DT C y12 C m11 þ C m12 C y22 C y12 y C y23 y ¼ a þ a 2 C m11 þ C y22 C y22 1 C x22 2 !  m  2 2C 12 C y12 y y y m þ 1 a1 þ a3 þ C 11  y a1 þ C m11 C 22 ! )  2 y y C 12 C 23 y 2C m12 m y m m þ C 12  a2 þ C 11 þ C 12  m a1 C y22 C 11 

 y C y22 C m11 C 12 y C y23 y ¼ DT a þ a y m C 11 þ C 22 C y22 1 C y22 2  m  2C 12 þ þ 1 am1 þ ay3 C m11

ðA:3bÞ

ðA:3cÞ

Given the thermal–mechanical constitutive model in the local YR-cell coordinate system (Eq. (A.1)), the constitutive behavior in the global unit-cell coordinate system can be obtained by a rotation through a hy angle about the negative x-axis. The thermal–mechanical constitutive model of the YR-cell in the global unit-cell coordinate system can be expressed as 9 8 2 ½C 11 YR ½C 12 YR ½C 13 YR rxx > > > > > > > > ryy > 6 ½C 12 YR ½C 22 YR ½C 23 YR > > > 6 = < 6 ½C 13 YR ½C 23 YR ½C 33 YR rzz ¼6 6 0 0 0 > > rxy > > 6 > > > > 4 ½C 15  > > ½C 25 YR ½C 35 YR r > > yz ; YR : 0 0 0 rzx YR 8 9 8 t 9 rxx > exx > > > > > > > > > > > > > rt > > > > eyy > > > > > = < > < yy t = ezz rzz þ  0 > > > > > > cxy > > > > > > > > > t > > > > > r c > > > yz ; : ; : yz > czx YR 0 YR

0 0 0 ½C 44 YR 0 ½C 46 YR

½C 15 YR ½C 25 YR ½C 35 YR 0 ½C 55 YR 0

0 0 0

3

7 7 7 7 ½C 46 YR 7 7 0 5 ½C 66 YR

ðA:4Þ

Given Cy = cos(hy), Sy = sin(hy), and CSy = cos(hy)sin(hy), the coefficients of the stiffness matrix defined in Eq. (A.4) can be calculated by ½C 11 YR ¼ ½C 11 nn YR 2 2 11 ½C 12 YR ¼ C y ½C 12 nn YR þ S y ½C ns YR  þ C 2y ½C 11 ½C 13 YR ¼ S 2y ½C 12 ns YR nn YR12  ½C 15 YR ¼ CS x ½C nn YR  ½C 11 ns YR 4 11 ½C 22 YR ¼ C 4y ½C 22 nn YR þ S y ½C ss YR   22 þ 2CS 2x ½C 21 ns YR þ 2½C ss YR  2 22 ½C 23 YR ¼ ðC 4y þ S 4y Þ½C 21 ns YR þ CS y ½C nn YR  þ ½C 11   4½C 22 ss ss YR n YR 2 21 ½C 25 YR ¼ CS y S y ð½C ns YR  ½C 11 ss YR Þ

ðA:5aÞ ðA:5bÞ ðA:5cÞ ðA:5dÞ ðA:5eÞ ðA:5fÞ

21 2 2 22 þ C 2y ð½C 22 nn YR  ½C ns YR Þ  2ðC y  S y Þ½C ss YR

o ðA:5gÞ

4 11 ½C 33 YR ¼ S 4y ½C 22 nn YR þ C y ½C ss YR

 þ 2CS 2x ½C 21  þ 2½C 22 ss YR n ns YR  11 ½C 35 YR ¼ CS x C 2y ½C 21 ns YR  ½C ss YR

ðA:5hÞ



o

 21 2 2 22 þ 2 C   ½C   S  ½C þ S 2y ½C 22 nn YR ns YR y y ss YR ðA:5iÞ

1036

J. Lua / Composites: Part A 38 (2007) 1019–1037

2 33 ½C 44 YR ¼ C 2y ½C 33 nn YR þ S y ½C ss YR

33  ½C 46 YR ¼ CS y ½C nn YR  ½C 33 ss YR   11 21 ½C 55 YR ¼ CS 2y ½C 22 nn YR þ ½C ss YR  2½C ns YR

þ C 2y  S 2y ½C 22 ss YR

ðA:5jÞ ðA:5kÞ

ðA:5lÞ

2 33 ½C 66 YR ¼ C 2y ½C 33 ss YR þ S y ½C nn YR

ðA:5mÞ

In addition, the thermal mismatch induced residual stress components in Eq. (A.4) are given by 0

½rtxx YR ¼ ½rtxx YR

ðA:6aÞ

t0

t0

ðA:6bÞ

t0

½rtzz YR ¼ S 2y ½ryy YR þ C 2y ½rzz YR

ðA:6cÞ

0 ½rtyz YR

ðA:6dÞ

½rtyy YR ¼ C 2y ½ryy YR þ S 2y ½rzz YR t0

¼

0 CS y ð½rtyy YR



0 ½rtzz YR Þ

A.2. Homogenization of RR-cell The RR-cell shown in Fig. 4 is filled with the resin alone. Given the isotropic nature of the resin constituent, the thermal–mechanical model of the RR-cell in the global unit-cell coordinate system can be written as 8 9 rxx > > > > > > > > > > > r > > yy > > > > > > > > < rzz > =

2

> rxy > > > > > > > > > > > > > > > r yz > > > > > : > ; rzx RR

½C 11 RR ½C 12 RR ½C 12 RR

0

0

6 6 ½C 12 RR ½C 11 RR ½C 12 RR 0 0 6 6 6 ½C 12 RR ½C 12 RR ½C 11 RR 0 0 6 ¼6 6 0 0 0 ½C 44 RR 0 6 6 6 0 0 0 0 ½C 44 RR 4



0 8 9 exx > > > > > > > > > > > > e > yy > > > > > > > > < ezz > = > cxy > > > > > > > > > > > > > > c yz > > > > > > : > ; czx RR

þ

0 0 8 t 9 rxx > > > > > > > > > t > > > r > yy > > > > > > > > < rt > = zz > 0 > > > > > > > > > > > > > > > 0 > > > > > : > ; 0 RR

0

0

0 0 0 0 0

3 7 7 7 7 7 7 7 7 7 7 7 5

½C 44 RR

ðA:7Þ

Given the elastic modulus, E, and Poisson’s ratio, m, of the resin material, the coefficients of the stiffness matrix in Eq. (A.7) can be expressed as Eð1  mÞ ð1 þ mÞð1  2mÞ Em ¼ ð1 þ mÞð1  2mÞ E ¼ 2ð1 þ mÞ

½C 11 RR ¼

ðA:8aÞ

½C 12 RR

ðA:8bÞ

½C 44 RR

ðA:8cÞ

The thermal stress components resulting from a uniform temperature increase (DT) in the RR-cell are given by ½rtxx RR ¼ ½rtyy RR ¼ ½rtzz RR ¼ f½C m11 RR þ 2½C m12 RR gam1 DT ðA:9Þ

References [1] Ishikawa T. Antisymmetric elastic properties of composite plates of satin weave cloth. Fiber Sci Technol 1981;15:127–45. [2] Ishikawa T, Chou TW. One dimensional micromechanics analysis of woven hybrid composites. AIAA J 1983;21:1714–21. [3] Ishikawa T, Chou TW. Stiffness and strength behavior of fabric composites. J Mater Sci 1982;17:3211–20. [4] Cox BN, Carter WC, Fleck NA. A binary model of textile composites. I Formulation. Acta Metall Mater 1994;42:3463–79. [5] Raju IS, Wang JT. Classical laminate theory models for woven-fabric composites. J Compos Technol Res 1994;16:289–303. [6] Tabiei A, Jiang Y. Woven fabric composite material model with material nonlinearity for nonlinear finite element simulation. Int J Solids Struct 1999;36(18):2757–71. [7] Whitcomb J, Woo K, Gundapaneni S. Macro finite element analysis of textile composites. J Compos Mater 1994;27:607–18. [8] Tan P, Tong L, Stever GP. Micromechanics models for the elastic constants and failure strengths. Compos Struct 1999;47:797–804. [9] Carvelli V, Poggi C. A homogenization procedure for the numerical analysis of woven fabric composites. Composites: Part A 2001;32: 1425–32. [10] Choi J, Tamma KK. Woven fabric composites – Part I: Predictions of homogenized elastic properties and micromechanical damage analysis. Int J Numer Meth Eng 2001;50:2285–98. [11] Key C, Lua J. Constituent based analysis of composite materials subjected to fire conditions. Composites: Part A 2006;37(7):1024–39. [12] Bogdanovich AE, Pastore CM. Material-smart analysis of textilereinforced structures. Compos Sci Technol 1996;56:291–309. [13] Lomov SV, Gusakov AV, Huysmans G, Prodromou A, Verpoest I. Textile geometry preprocessor for meso-mechanical models of woven composites. Compos Sci Technol 2000;60:2083–95. [14] Lomov SV, Huysmans G, Luo Y, Parnas RS, Prodromou A, Verpoest I, et al. Textile composites: modeling strategies. Composites: Part A 2001;32:1379–94. [15] Verpoest I, Lomov SV. Virtual textile composites software WiseTex: integration with micro-mechanical, permeability and structural analysis. Compos Sci Technol 2005;65:2563–74. [16] Lua J, O’Brien J. Fire simulation for woven fabric composites with temperature and mass dependent thermal–mechanical properties. In: Proceeding of Composites in Fire 3, September 9–10, 2003, Newcastle upon Tyne, UK. [17] Lua J, Key CK, O’Brien J, Lopez-Anido RA, El-Chiti FW, Dagher H, et al. Validation of a woven fabric analyzer using composite material strain measurements from a digital image correlation system. In: Proceedings of the joint American society for composites/ American society for testing and materials committee D30. October 17–20, 2004, Atlanta, GA. [18] Lua J, O’Brien J, Key C. A temperature and mass dependent thermal model for fire response prediction of marine composites. Composites: Part A 2006;37(7):1005–41. [19] Lua J, Key C, Schumacher SC, Hansen A. Rate dependent multicontinuum progressive failure analysis of woven fabric composite structures under dynamic impact. Shock Vib 2004;11(2): 103–17. [20] Bahei-El_Din YA, Rajendran AM, Zikry MA. A micromechanical model for damage progression in woven composite systems. Int J Solids Struct 2004;41:2307–30. [21] Zako M, Uetsuji Y, Kurashiki T. Finite element analysis of damaged woven fabric composite materials. Compos Sci Technol 2003;63: 507–16. [22] Tabiei A, Ivanov I. Materially and geometrically non-linear woven composite micro-mechanical model with failure for finite element simulations. Int J Non-Linear Mech 2004;39:174–88. [23] Key CT, Six RW, Hansen AC. A three-constituent multicontinuum theory for woven fabric composite materials. Compos Sci Technol 2003;63:1857–64.

J. Lua / Composites: Part A 38 (2007) 1019–1037 [24] Tabiei A, Yi W. Comparative study of predictive methods for woven fabric composite elastic properties. Compos Struct 2002;58: 149–64. [25] Chung PW, Tamma KK. Woven fabric composites-developments in engineering bounds, homogenization and applications. Int J Numer Meth Eng 1999;45:1757–90. [26] Lopez-Anido R, El-Chiti F, Muszynski L, Dagher H, Thompson L, Hess, P. Composite material testing using a 3-D digital image correlation system. Composites 2004, American Composites Manufacturers Association, Tampa, Florida.

1037

[27] Fayad GN. Probabilistic finite element analysis of marine grade composites. Master Thesis 2005; Department of Civil and Mechanical Engineering, The University of Maine. [28] Halpin JC. Primer on Composite Materials Analysis. second ed., Revised. Technomic Publishing Company, Inc.; 1992. [29] Mayers JS. Micromechanics based failure analysis of composite structural laminates. Technical Report of Naval Surface Warfare Center Carderock Division 1999, NSWCCD-65-TR-1999/15. [30] Matzenmiller A, Lublinear J, Taylor RL. A constitutive model for anisotropic damage in fiber-composites. Mech Mater 1995;20:125–52.