Multiscale computational homogenization of woven composites from microscale to mesoscale using data-driven self-consistent clustering analysis

Composite Structures 220 (2019) 760–768

Contents lists available at ScienceDirect

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

Multiscale computational homogenization of woven composites from microscale to mesoscale using data-driven self-consistent clustering analysis

T

⁎

Xinxing Han, Chenghai Xu , Weihua Xie, Songhe Meng National Key Laboratory of Science and Technology for National Defence on Advanced Composites in Special Environments, Harbin Institute of Technology, Harbin 150001, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Multiscale computational homogenization Woven composites Self-consistent clustering analysis Reduced order model

Compared with the traditional phenomenological method, the multiscale simulation has significant advantages. This paper presents a methodology and computational homogenization framework to predict the macroscale behavior of woven composites based on fiber and matrix in microscale. The major challenge conducting multiscale analysis is the huge computational cost. To improve the efficiency, one of the reduced order models, which is called the data-driven self-consistent clustering analysis (SCA), is introduced and the multiscale framework is proposed by integrating two SCA solvers from different scales. The macroscale performance of 4-H satin weave carbon/carbon composites is investigated using the proposed framework. In order to reconstruct a real microstructure representative volume element (RVE), both microscale and mesoscale architectures are observed using scanning electron microscopy (SEM) and optical microscopy, and statistical geometry features are obtained. In addition, the SCA method is also verified by comparing the results with the finite element method (FEM). The uniaxial tension process is simulated using the multiscale approach, and strain/stress fields in both mesoscale and microscale can be captured simultaneously. Moreover, the uniaxial tensile experiments are also carried out to validate this framework, which shows high efficiency and great accuracy.

1. Introduction Woven composites are made of yarn and matrix. Given their outstanding mechanical performance, especially high stiffness and strength to weight ratio, they have been widely used in space, automotive and other industry fields [1,2]. These outstanding properties depend on the fibrous reinforcements and their mesoscale architecture. However, it is still a challenge to conduct efficient and accurate mechanical analysis, because of their complex microstructure and strong heterogeneous features. Traditional phenomenological models [3–5], which characterize the average behavior of the composites, cannot take the mechanical behavior of each component and its microstructure into account. In particular, if the constituents are characterized by the nonlinear response, the overall model will be more complex, which leads to burdensome experimental characterization and tedious calibration. In addition, these models often fail to capture the highly localized nonlinear material behavior. In contrast, the multiscale methods [6–9] provide the possibility to link the microstructure to the overall properties. The representative volume element (RVE) [10–14] is a popular multiscale approach, which is widely used to determine the performance of

⁎

homogenized equivalent material. For the woven composites, most RVE models in the literature [15–17] consider the yarn and matrix as homogenized material, and the finite element method (FEM) is usually used to conduct the computation of the RVE problem. These models cannot capture the local fields in the fiber microscale. Wang, L. et al. [18] bridged the mesoscale and microscale using the stress amplification factors (SAFs), where only several key points in the microscale are selected to obtain the SAFs. To take into account the microscale, the effective properties of the tow are predicted through the RVE model in microscale, this method is also widely used. [19] Patel, D. K. et al. [20] linked the microscale to mesoscale through a 6 × 6 transformation matrix. Huang, T. et al. [21] predicted the properties of fiber tows using the analytical approach. The SAFs, transformation matrix and analytical approach are very efficient, however, the microscale fields are calculated less precisely, as a result, it cannot capture the local effects very well. Computational homogenization (CH) framework [22,23] provides a basic idea to establish the relation between mesoscale and microscale in this paper. For every material point in mesoscale, a corresponding microscale RVE will be solved and then homogenized in response to this material point. The challenge to realize the computational

Corresponding author. E-mail address: [email protected] (C. Xu).

https://doi.org/10.1016/j.compstruct.2019.03.053 Received 9 January 2019; Received in revised form 11 March 2019; Accepted 18 March 2019 Available online 20 March 2019 0263-8223/ © 2019 Elsevier Ltd. All rights reserved.

Composite Structures 220 (2019) 760–768

X. Han, et al.

data science field is introduced. After that, the degrees of freedom are significantly reduced. Then the corresponding discrete incremental Lippmann-Schwinger integral equations are solved using NewtonRaphson’s iteration algorithm to ensure the convergence.

homogenization is that it is highly expensive as many RVEs need to be calculated for every iteration. The finite element method (FEM) is an accurate solving approach, but it will lead to very high computational burden and a longer computational time if used at both microscale and mesoscale level [24]. To reduce the computational cost, the analytical micromechanical method [25–28] is a very good option, but it will lose accuracy when facing nonlinear problems with complex microstructure due to their many assumptions. In addition, some numerical methods are also proposed. The fast Fourier transform (FFT)-based method [29,30] is more efficient than FEM, but there will be some convergence problems for high phase contrast and complex microstructure with the nonlinear process. The transformation field analysis (TFA) [31–33], the nonuniform transformation field analysis (NTFA) [34–37] and proper orthogonal decomposition (POD) [38–40] are reduced order models. TFA groups the material points according to the geometry features and assumes the strain field is uniform in each group, which will lose its predictive capability when the microstructure is complex, such as woven composites. The NTFA improves the TFA method and nonuniform strain field is considered for each group, but the strain field modes need to be predefined, and these modes are not always known in advance. POD is a mathematical method to solve the partial differential equations, where the solution modes also need to be assumed. Again, this information is also not always known for complex microstructure. The data-driven self-consistent clustering analysis (SCA) [41] is an effective and efficient method to solve the RVE problem with complex microstructure for irreversible processes such as inelastic deformation. This makes it particularly attractive for integration in the framework of computational homogenization simulation. This method uses a clustering algorithm to reduce the overall degrees of freedom of the RVE, which is called offline stage. Based on the reduced order RVE database generated in offline stage, the online stage is conducted by solving the discretized incremental Lippmann-Schwinger integral equations to obtain the stress and strain fields. This efficient method has been used for 2-dimensional (2D) two-phase composites and 3-dimensional (3D) hard inclusion materials considering nonlinear elastoplastic damage softening effect [42] and crystal plasticity analysis [43]. This shows that the SCA has good efficiency as well as high accuracy. Additionally, Bessa, M. A., et al. [44] used this method for data mining and uncertainty analysis. The overall objective of this work is to develop an efficient multiscale computational homogenization framework for woven composites based on the SCA reduced order model. Three scales are involved in this paper: macroscale, mesoscale and microscale. The macroscale properties are obtained through the homogenization of the mesoscale RVE, the response of every material point in the yarn is determined by the solution of corresponding microscale unidirectional (UD) RVE. So, the details of the methodology and framework are introduced in Section 2. In Section 3, the process of solving the unidirectional (UD) RVE in microscale and the efficiency are presented. The results are also verified by comparing with the FEM method in this section. In addition, the multiscale computational homogenization (CH) and experimental validation are fully illustrated in Section 4. Moreover, woven RVE is generated through the mesoscale observation and measurement, and the average response is calculated using the proposed framework. In order to validate the entire computation, the typical uniaxial tensile experiment is carried out and stress-strain curves are used as the comparison. Concluding remarks are provided in Section 5.

2.1.1. Offline stage: model reduction To cluster the material points with the similar mechanical response into the same cluster. Firstly, the similarity between two material points is measured using the strain concentration tensor Am (x) , which is defined as following form:

εm (x) = Am (x): εM

in

(1)

Ω

where εM is the elastic average strain under periodic boundary conditions imposed on the RVE, and εm is the elastic local strain at point x in the RVE domain Ω. After computing the strain concentration tensor Am (x) , the clustering method is used to cluster the material points. The main assumption is that the local strain field is uniform within each cluster. This assumption reduces the degrees of freedom for the following corresponding incremental Lippmann-Schwinger integral equations. 2.1.2. Online stage: solving the incremental Lippmann-Schwinger integral equations The equilibrium equation in the RVE can be rewritten as a continuous incremental Lippmann-Schwinger integral equation by introducing a homogeneous isotropic linear elastic reference material:

Δεm (x) +

∫Ω Φ0 (x,x′): [Δσm (x′) − C0: Δεm (x′)] d x′ − Δε0 = 0

(2)

Δε0

is the uniform far-field strain increment in the RVE domain where which controls the evolution of the local strain. Δεm (x) and Δσm (x) are increment variables of local strain and local stress respectively. Φ0 (x, x′) is the Green’s function which indicates the strain at x contributed by a concentrated external stress at x′ in the homogeneous reference material. The reference material is isotropic linear elastic and the stiffness tensor C0 depends on two independent Lame parameters λ0 and μ0 :

C0 = λ0I ⊗ I + μ0 II

(3)

where I is the second-order identity tensor, II is the symmetric fourthorder identity tensor and ⊗ means tensor product. Solving the Eq. (2) at every material point in RVE domain will be very computationally expensive. So, the continuous increment Lippmann-Schwinger integral equation is homogenized in each cluster to reduce the degrees of freedom: 1 c I | Ω| 1 c I | Ω|

∫ΩI ∫Ω

Φ0 (x,

∫ΩI Δεm (x) d (x)+

x′): [Δσm (x′) − C0 : Δεm (x′)] d x′d x − Δε0 = 0

(4)

cI

where is the volume fraction of I th cluster, |Ω| indicates the volume of the RVE domain Ω. With the assumption of the uniform strain field within each cluster, piecewise uniform distribution is used to express the strain and stress fields: K

K

Δεm (x) =

∑ χ I (x)ΔεmI , Δσm (x) = ∑ χ I (x)Δσ mI I=1

where K is the number of clusters in RVE, function in the I th cluster, defined as:

2. Methodology and framework

(5)

I=1

χI

(x) is the characteristic

{

I χ I (x) = 1 x ∈ Ω 0 otherwise

2.1. Data-driven self-consistent clustering analysis

(6)

which has the averaging relationship as following:

∫Ω χ I (x)[∙] d x ≡ ∫Ω

The data-driven self-consistent clustering analysis (SCA) is an efficient RVE solving method, which groups the material points with the similar mechanical response into the same cluster using cluster analysis technology in offline stage. In this stage, the cluster analysis from the

I

[∙] d x

(7)

where [∙] indicates any quantity to be averaged in the I th cluster domain. 761

Composite Structures 220 (2019) 760–768

X. Han, et al.

aspects. The first one is the RVE clusters results using data compression algorithm-cluster analysis. For the composites, the reinforcement phase and the matrix are grouped separately. In this paper, the k-means clustering algorithm [45] is used to cluster the microstructure. The second one is the reference material stiffness matrix C0 in Lippmann-Schwinger equations, which is set to be approximately the ¯ even after nonlinear desame as the macroscopic tangent modulus C formation occurs. Based on micromechanics, the macroscopic tangent ¯ can be expressed directly as: modulus C

Substituting Eq. (5) into Eq. (4) and using the relationship (7), the reduced order discrete incremental Lippmann-Schwinger integral equation for Ith cluster can be written as: K I Δεm +

∑ D IJ : [Δσ mJ − C0: ΔεmJ] − Δε0 = 0 J =1

(8)

where D IJ is the fourth-order interaction tensor which depends on the microstructure of the RVE. It can be expressed as:

1 c I |Ω|

D IJ =

∫Ω ∫Ω χ I (x) χ J (x′) Φ0 (x, x′) d x′d x

k

(9)

¯ = C

For the whole RVE domain, a discrete incremental LippmannSchwinger integral equation set is organized as: K ∑J = 1 K ∑J = 1

1 ⎧ Δεm

J [Δσ m

I=1

is the tangent modulus of the material in the I th cluster where which is an output of the local constitutive model. AI is the average strain concentration matrix of the I th cluster and, for the 3-dimensional case, AI is a 6 by 6 matrix with 36 components, which can be obtained under 6 orthogonal loading conditions respectively in the offline stage. ¯ is always anisotropic for woven composites The macroscopic tangent C in the global system. Projection method proposed in [43] is used to find the closest isotropic material stiffness matrix C0 based on known ani¯: sotropic material stiffness matrix C

(10)

The Newton-Raphson’s iteration algorithm is adopted to solve Eq. (10) in the online stage, which can improve the accuracy and convergence, especially for nonlinear material behavior. For simplicity, the reference material is chosen as isotropic linear elastic, so that the Green’s function can be explicitly expressed in Fourier space as following [30]:

λ0 + μ0 1 ̂1 0 2 Φ̂ (ξ) = Φ (ξ) + 0 0 Φ̂ (ξ) μ (λ + 2μ0 ) 4μ0

̂1

¯ ) J + 1 (K: :C ¯ )K C0 = (J: :C 5

Iijkl = (11)

̂2

1 = 2 (δik ξ j ξl + δil ξ j ξk + δjl ξi ξk + δjk ξi ξl ) |ξ|

2

̂ (ξ ) = − Φijkl

Data-driven self-consistent clustering analysis is a numerical method to solve the RVE problem. As for the woven composites, multiple material scales are involved, including the fiber/matrix microscale, the yarn architecture mesoscale and the specimen macroscale. Computational homogenization provides an idea to build the bridge among these scales. The response in macroscale is calculated by homogenizing the woven RVE in mesoscale, and each material point within the yarn domain is connected to a corresponding unidirectional (UD) RVE in microscale. So, the behavior of these material points is determined through averaging the UD RVE in microscale, as shown in the following Fig. 1. The performance of the woven composites in macroscale can be predicted only based on the material behavior of fiber and matrix in microscale. While, the real microstructure of the composites is very complex, to simplify the model, only the main features are considered. The interface and defects, including voids and microcracks, are not considered in this paper. Finally, the algorithm of SCA × SCA multiscale simulation is given in the following Box I.

(12)

(13)

So, the part of the integration in Eq. (9) can be easily calculated based on the Fourier transformation:

∫Ω χ I (x′) Φ0 (x, x′) d x′ = F−1 (χ Î (ξ) Φ̂0 (ξ))

(14)

Substituting Eq. (11) into Eq. (14) and using the interaction tensor expression (9), the D IJ can be rewritten as the sum of two parts:

D IJ =

λ0 + μ0 1 IJ D1 + 0 0 D2IJ 4μ0 μ (λ + 2μ0 )

where

D1IJ

D1IJ =

1 c I |Ω|

∫Ω χ I (x)[F−1 (χ Î (ξ) Φ̂1 (ξ))] d x

(16)

D2IJ =

1 c I |Ω|

∫Ω χ I (x)[F−1 (χ Î (ξ) Φ̂2 (ξ))] d x

(17)

D1IJ

and

D2IJ

(15)

have the following form:

Box I. Algorithm of SCA×SCA multiscale simulation 1. Database generation and model reduction in the offline stage. 1.1. Build the woven RVE and UD RVE geometric model using the scanning electron microscope image data. 1.2. Reduce the order of woven RVE and UD RVE based on cluster analysis. 1.3. Calculate the interaction tensor D IJ of woven RVE and UD RVE respectively. 1.4. Save the above two reduced order RVEs and interaction tensors as the database. 2. SCA × SCA multiscale simulation in the online stage. 2.1. Apply the macro-strain to the woven RVE 2.2. Start Newton-Raphson iteration and give the meso-strain increment. 2.3. If the meso-strain increment is from matrix material, use the matrix material law to obtain the meso-stress increment. 2.4. If the meso-strain increment is from yarn material, go to 2.4.1. 2.4.1. Apply meso-strain increment to the UD RVE. 2.4.2. Start Newton-Raphson iteration and give the micro-strain increment. 2.4.3. If the micro-strain increment is from matrix material, use the matrix material law to obtain the micro-stress increment,

D2IJ

and do not depend on the material properties, they can be As calculated in the offline stage. Only the coefficients before them need to be updated if the Lame parameters λ0 and μ0 changed. Meanwhile, the macroscopic boundary conditions are also discretized and the discrete form of the macro-strain constraints can be written as: K

∑ c I ΔεmI = ΔεM or

Δε0 = ΔεM

I=1

(18)

The discrete form of the macro-stress constraints become: k

∑ c I Δσ mI = ΔσM I=1

1 1 (δik δjl + δil δjk ), Jijkl = δij δkl 2 3

2.2. SCA × SCA multiscale simulation framework

ξi ξ j ξk ξl |ξ|4

(21)

where K = I − J , I and J have components as

where Φ (ξ) and Φ (ξ) can be written as: 1 ̂ (ξ ) Φijkl

(20)

I Calg

J ] C0 : Δεm

+ − − Δε0 = 0 D1J : ⎪ 2 J J J 2 0 ⎪ Δεm + D : [Δσ m − C : Δεm] − Δε0 = 0 ⎨⋮ ⎪ K J J ⎪ Δεm + ∑KJ = 1 D KJ : [Δσ m ] − Δε0 = 0 − C0 : Δεm ⎩

∑ c I CaI lg ·AI

(19)

As shown in [41,42], the solution of the Eq. (10) depends on two 762

Composite Structures 220 (2019) 760–768

X. Han, et al.

Fig. 1. SCA × SCA multiscale simulation framework.

(b)

(a)

(c)

Fig. 2. Microstructure in yarn and UD RVE: (a), (b) micrographs of the yarn cross section, and (c) UD RVE model. 2.4.4. If the micro-strain increment is from fiber material, use the fiber material law to decide the micro-stress increment. 2.4.5. Check convergence, if it is satisfied, calculate the meso-stress increment by averaging the micro-stress increment, if not, go to 2.4.2 and loop. 2.5. Check convergence, if it is satisfied, calculate the macro-stress increment by averaging the meso-stress increment, if not, go to 2.2 and loop.

3. The application of SCA for yarn 3.1. UD RVE generation and model reduction The 4-H satin weave carbon/carbon composites are selected to illustrate the application of the SCA method and multiscale framework. For every material point in yarn, a corresponding UD RVE model is built to determine the behavior of this point. To generate the real UD RVE model, the statistical microstructure information is obtained from the 763

Composite Structures 220 (2019) 760–768

X. Han, et al.

(b)

(a)

0.25

FEM U-4-4 U-8-4 U-16-4

(a)

Stress(GPa)

0.20

22

23

Tension

0.15

0.10

Shear 0.05

0.00 0.000

Fig. 3. Clustering results using k-means: (a) matrix, and (b) fiber.

0.005

0.010

0.015

0.020

Strain

0.25

FEM U-4-4 U-4-8 U-4-16

(b)

Stress(GPa)

0.20

22

Tension

0.15

23

0.10

Shear

Fig. 4. The local frame of the fiber. 0.05

Table 1 Material properties of the constituents [46,47]. E1 (GPa) T300 carbon fiber Carbon matrix

362.7 10.6

E2 (GPa) 9.9 –

v12 0.41 0.15

v 23 0.45 –

0.00 0.000

G12 (GPa)

11 Am12 Am13 Am14 Am15 Am16 ⎤ ⎡ εM ⎤ 22 ⎥ Am22 Am23 Am24 Am25 Am26 ⎥ ⎢ εM ⎥⎢ ⎥ 33 ⎥ Am32 Am33 Am34 Am35 Am36 ⎥ ⎢ εM 42 43 44 45 46 ⎥ ⎢ 23 ⎥ Am Am Am Am Am ⎥ ⎢ εM ⎥ 13 ⎥ Am52 Am53 Am54 Am55 Am56 ⎥ ⎢ εM ⎥⎢ ⎥ 12 62 63 64 65 66 Am Am Am Am Am ⎥ ⎣ εM ⎥ ⎦ ⎦⎢

0.010

0.015

0.020

Strain

19.9 –

Fig. 5. Stress-strain curves of FEM and SCA under tension and shear: (a) change the cluster number in the matrix, and (b) change the cluster number in fiber.

cross section micrograph of the yarn by scanning electron microscopy (SEM), as shown in Fig. 2. The key parameters, such as the volume fraction and the diameter of the fiber, are identified using the image recognition algorithm. As a result, the fiber volume fraction in the yarn is 88.12% and the average diameter of the fiber is 7 μm. Since the fiber is well organized, hexagonal packing RVE is adopted to model the distribution of the fiber. Since the periodic boundary conditions will be applied to the UD RVE, only one layer voxel elements are built along the fiber direction, which is enough to include the whole features in this direction. Finally, the UD RVE is generated with 101 voxel elements in width and 176 voxel elements in length, and the volume fraction is 88.14%. In the 3D case, the strain concentration tensor has 36 components, which can be expressed in the following form, 11 11 ⎡ εm ⎤ ⎡ Am ⎢ ε 22 ⎥ ⎢ A 21 ⎢ m⎥ ⎢ m ⎢ εm33 ⎥ ⎢ Am31 ⎢ 23 ⎥ = ⎢ 41 ⎢ εm ⎥ ⎢ Am ⎢ ε 13 ⎥ ⎢ A 51 ⎢ m⎥ ⎢ m 12 61 ⎢ ⎣ Am ⎦ ⎢ ⎣ εm ⎥

0.005

Table 2 Computational cost and efficiency comparison. Elements/Clusters

FEM U-4-4 U-8-4 U-16-4 U-4-8 U-4-16 Max. Speed Up

17,776 8 12 20 12 20 2222

Equations

54,162 48 72 120 72 120 1128.4

Computational Time/s Uniaxial Tension

Pure Shear

482.7 0.54 0.61 0.92 0.59 0.99 893.9

474.6 0.55 0.56 0.71 0.61 0.73 862.9

elements will be grouped into different clusters using the k-means cluster analysis. As a result, the degrees of freedom are significantly reduced. The process of clustering is carried out in the fiber and matrix domains separately. Fig. 3 shows the clustering results, and the number of clusters is set to be 4 in both matrix and fiber domains. Since the cluster analysis is based on the strain concentration tensor, the voxel elements with the most similar mechanical response will be grouped into the same cluster. This process corresponds to the uniform strain field assumption within every single cluster of the SCA method. In addition, the shape of the cluster is irregular, which gives it the capability to handle the complex microstructure.

(22)

Six orthogonal loading conditions are used to compute the strain concentration tensor Am . The microscale strain field within the UD RVE domain is calculated using the Finite element method (FEM) based on the eight-node brick with reduced integration (C3D8R) element. After the strain concentration tensor is obtained at every voxel element, these 764

Composite Structures 220 (2019) 760–768

X. Han, et al.

(a)

(b)

(d)

(c)

Fig. 6. Woven RVE architecture: (a)(b) microstructure observation using optical microscopy, and (c)(d) woven RVE reconstruction. ν

Parameters Yarn dimensions

width 2a¯

Yarn spacing

thickness 2b¯ d¯1 d¯2

RVE

length w¯ height h¯

Mean value (μm)

Standard dev (μm)

771.8 70.8

69.5 6.2

800

32.25

29.2

1.4

3200 160

129.0 12.4

ν

1 ⎡ E −E −E ⎢ ν 1 ν ε ⎡ 11 ⎤ ⎢− E E − E ε ⎢ 22 ν ν ⎢ ⎥ 1 ⎢ ε33 ⎥ ⎢− E − E E = ⎢ γ ⎢ 23 ⎥ 0 0 ⎢γ ⎥ ⎢ 0 ⎢ 13 ⎥ ⎢ γ 0 0 0 ⎣ 12 ⎦ ⎢ ⎢ 0 0 ⎢ 0 ⎣

Table 3 Woven RVE architecture parameters.

0

0

0

0

0

0

2 + 2ν E

0

0

2 + 2ν E

0

0

0 ⎤ ⎥ 0 ⎥ σ11 ⎥ ⎡ σ22 ⎤ 0 ⎥ ⎢σ ⎥ 33 ⎥⎢ ⎥ 0 ⎥ ⎢ σ23 ⎥ ⎥ ⎢ σ13 ⎥ ⎥ 0 ⎥⎢ ⎣ σ12 ⎦ ⎥ 2 + 2ν ⎥ E ⎦

The material law of the fiber is orientation dependent, so, the local frame is built in Fig. 4. In the 2–3 plane, the properties of the fiber are assumed to be isotropic, so, transversely isotropic material law is adopted for the fiber, which has the following form:

3.2. Microscale material model Since carbon/carbon composites have brittle properties, the elastic material law is considered before the damage occurs. For the matrix, isotropic elastic is used and the stress-strain relationship can be expressed in the following familiar form:

(b)

(a)

(23)

(c)

Fig. 7. Cluster results for woven RVE: (a) Cluster results for the matrix, (b) Cluster results for yarn, and (c) Cluster results for a single yarn. 765

Composite Structures 220 (2019) 760–768

X. Han, et al.

(a)

(c)

(b)

Fig. 8. Uniaxial tension experiment: (a)(b) uniaxial tension specimen, and (c) experimental process.

(a)

0.10

Stress(GPa)

0.08

Table 4 Uniaxial tensile modulus comparison.

Experiment W-8-128 W-16-128 W-16-256

Method

Average uniaxial tensile modulus (GPa)

Error

Experiment Multiscale simulation

98.7 98.4

– 0.3%

0.06 ν

0.04

0.02

0.00 0.0000

0.0002

0.0004

0.0006

0.0008

ν

1 21 21 ⎡ E1 − E2 − E2 ⎢ ν12 ν 1 − E23 ⎢− ε E2 2 ⎡ 11 ⎤ ⎢ E1 ε 1 ⎢ 22 ⎥ ⎢− ν12 − ν23 E2 E2 ⎢ ε33 ⎥ ⎢ E1 ⎢ γ23 ⎥ = ⎢ 0 0 0 ⎢γ ⎥ ⎢ ⎢ 13 ⎥ ⎢ 0 0 ⎣ γ12 ⎦ ⎢ 0 ⎢ 0 0 0 ⎢ ⎣

0.0010

Strain

(b)

0

0

0

0

0

0

2 + 2ν 23 E2

0

0

1 G12

0

0

0 ⎤ ⎥ 0 ⎥ σ11 ⎥ ⎡σ ⎤ 0 ⎥ ⎢ σ22 ⎥ ⎥ ⎢ 33 ⎥ 0 ⎥ ⎢ σ23 ⎥ ⎥ ⎢ σ13 ⎥ σ12 ⎥ 0 ⎥⎢ ⎥⎣ ⎦ 1 ⎥ G12 ⎥ ⎦

(24)

The material parameters of the fiber and matrix are listed in Table 1. Carbon fiber utilized in fabrication undergo property changes due to the high temperature heating cycles during the carbon/carbon composites process. As a result, the modules of the T300 carbon fiber are increased. 3.3. Results and discussion The comparison of stress and strain curves are plotted in Fig. 5, and we can conclude the results are in good agreement with each other. The U-8-4 indicates that there are 8 clusters in the matrix domain and 4 clusters in the fiber domain. With the number of clusters increasing in the matrix and in the fiber respectively, the results are converged to FEM both in tension and shear, which can demonstrate stability and accuracy of the SCA method. The computational cost and efficiency comparison are presented in Table 2. All the tasks are run on Intel(R) Xeon(R) processor with 24 cores and 64 GB memory, and FEM and SCA are implemented using commercial software ABAQUS and MATLAB. As for the SCA, online running time is listed since only this stage needs to be integrated into computational homogenization framework, while the FEM method needs the full analysis. It can be concluded that the FEM method is more expensive than SCA method. The FEM requires 1128 times

Fig. 9. Multiscale computational results: (a) Macroscale stress curves, and (b) Stress fields in mesoscale and microscale.

766

Composite Structures 220 (2019) 760–768

X. Han, et al.

Table 5 Efficiency comparison for the different multiscale framework. Method

Mesoscale cost

Meso step

Microscale cost

Micro step

Total cost

Time (s)

FEM × FEM SCA × SCA

1,638,400 8 + 128 16 + 128 16 + 256

10 10 10 10

17,776 4+8 4+8 4+8

1 1 1 1

291.24 billion 16,320 17,280 32,640

– 10.46 11.05 20.95

stress-strain curves with different numbers of clusters and experiments are plotted in Fig. 9(a). It can be concluded that the computational results are in good agreement with experimental data. Moreover, the computational results converge when using different numbers of the cluster in matrix and yarn, which demonstrates the stability and accuracy of the SCA × SCA multiscale framework. Fig. 9(b) shows stress fields along the loading direction in both mesoscale and microscale. Every cluster in mesoscale corresponds to a UD RVE in microscale, and only three clusters are selected to plot the stress fields in microscale. In the yarn along the loading direction, the stress values in the fiber are bigger than those in the matrix. So, the fiber has a greater impact on the behavior of yarn. While, in the yarn which is vertical to the loading direction, the stress values in the matrix and the fiber are closed. In addition, the uniaxial tensile modulus from experiments and simulations are calculated in Table 4, which also shows good agreement. Efficiency comparison is shown in Table 5. The cost is characterized using the number of elements or clusters. So, the whole computational cost is 291.24 billion elements using FEM × FEM, which is hard to handle even using high performance computing (HPC) clusters. However, the computational cost is significantly reduced using SCA × SCA multiscale framework, which makes it efficient to handle this task.

equations than SCA method, which means the SCA saves much more memory. For the computational time, the SCA method is more than 800 times faster than the FEM in obtaining the RVE response under different loading conditions. Due to its significant efficiency, lots of RVEs in different loading cases can be solved in a very short time, which is essential to conduct the multiscale computational homogenization analysis. 4. Multiscale computational homogenization simulation 4.1. Woven RVE and model reduction The architecture of 4-H satin weave is characterized by optical microscopy, as shown in Fig. 6(a) and (b). In order to reconstruct woven RVE, the cross section of the yarn is assumed to be an ellipse, and the curved part of the centerline in the yarn is modeled using the sine function. These two features are depicted using Eq (25). Major involved architecture parameters are shown in Fig. 6(c) and (d), which can be identified using an image recognition algorithm. As a result, the mean values and standard deviations for each parameter are listed in Table 3. In addition, the interface and micro defects are not considered in this paper. Based on these parameters, 3D woven voxel RVE is generated in Fig. 6(c) and (d), which has 320 voxel elements in both width and length, 16 voxel elements in height.

y= x2 a¯ 2

d¯2 2

5. Conclusion Multiscale simulation is a powerful method to calculate the behavior of woven composites, which can consider microstructures in different levels. However, the major challenge to carry out this simulation is extremely expensive. To reduce the computational cost, the reduced order model, data-driven self-consistent clustering analysis (SCA), is introduced, which shows the great efficiency and high accuracy. Based on data-driven self-consistent clustering analysis (SCA), multiscale computational homogenization framework is proposed, which is called SCA × SCA. Using this framework, the macroscale behavior of 4-H satin weave carbon/carbon composites is predicted only using the properties of fiber and matrix, the computational results are in good agreement with experimental results. In addition, stress and strain fields can be captured both in microscale and mesoscale at the same time, this information provides reliable data to improve the analysis and design of the woven composites.

sin(4d¯1·x )

+

y2 b¯2

=1

(25)

For the woven composites, cluster analysis for the matrix is simpler as the material law is isotropic, so it only needs to be carried out once using k-means based on the strain concentration tensor Am . While, cluster analysis for the yarn is more complex, the orientation of the yarn is defined as tangent of the centerline in each yarn. So, it will be changed along the curved part of the centerline. Since one single cluster has only one material orientation, the clustering process should be conducted twice. The first time, a single yarn will be clustered based on the material orientation using k-means. The points with the closest material orientation will be grouped into the same cluster, then the material points in one of the orientation clusters will be clustered again according to the strain concentration tensor Am . After these two steps, the material points with the closest orientation and the most similar mechanical response will be grouped into one cluster, as shown in Fig. 7.

Acknowledgments The authors warmly acknowledge the financial support of the China Scholarship Council to enable this work, and the support of National Natural Science Foundation of China (Grant Nos. 11472092, 11672088 and 11502058), the National Basic Research Program of China (973 Program; Grant No. 2015CB655200) and Science & Technology on Reliability & Environmental Engineering Laboratory. The first author wants to thank Prof. Wing Kim Liu for studying in his research group in the Northwestern. In addition, we thank Erin Anderson for checking the writing of this paper in the Northwestern Graduate Writing Center.

4.2. Experimental test of the woven specimen Uniaxial tensile specimens are cut from a single part of the composite to make the yarn in one direction parallel to the loading axis. The geometry and dimensions of the specimen are shown in Fig. 8(a) and (b). In the gripping area of the specimen, thin aluminum pieces are stuck on the surface to protect the composites. The strain is measured using strain gauges and the loading rate is 1 mm/min.

References 4.3. Simulation results and experimental validation [1] Brandt J, Drechsler K, Arendts FJ. Mechanical performance of composites based on various three-dimensional woven-fibre preforms. Compos Sci Technol

Multiscale simulation results are presented in Fig. 9. The macroscale 767

Composite Structures 220 (2019) 760–768

X. Han, et al.

homogenization: trends and challenges. J Comput Appl Math 2010;234:2175–82. [24] Ozdemir I, Brekelmans WAM, Geers MGD. FE2 computational homogenization for the thermo-mechanical analysis of heterogeneous solids. Comput Methods Appl Mech Eng 2008;198:602–13. [25] Hashin Z. Citation classic – the elastic-moduli of heterogeneous materials. Curr Contents/Eng Technol Appl Sci 1988:14. [26] Hashin Z. Citation classic – the elastic-moduli of heterogeneous materials – the elastic-moduli of fiber-reinforced materials. Curr Contents/Phys Chem Earth Sci 1988:14. [27] Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc R Soc London Ser a-Math Phys Sci 1957;241:376–96. [28] Mura T. Micromechanics of Defects in Solids. Springer Science & Business Media; 2013. [29] Michel JC, Moulinec H, Suquet P. Effective properties of composite materials with periodic microstructure: a computational approach. Comput Methods Appl Mech Eng 1999;172:109–43. [30] Moulinec H, Suquet P. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput Methods Appl Mech Eng 1998;157:69–94. [31] Dvorak GJ, Baheieldin YA, Wafa AM. Implementation of the transformation field analysis for inelastic composite-materials. Comput Mech 1994;14:201–28. [32] Dvorak GJ, Baheieldin YA, Wafa AM. The modeling of inelastic composite-materials with the transformation field analysis. Modell Simul Mater Sci Eng 1994;2:571–86. [33] Dvorak GJ. Transformation field analysis of inelastic composite-materials. Proc R Soc London Ser a-Math Phys Eng Sci 1992;437:311–27. [34] Largenton R, Michel JC, Suquet P. Extension of the Nonuniform Transformation Field Analysis to linear viscoelastic composites in the presence of aging and swelling. Mech Mater 2014;73:76–100. [35] Michel JC, Suquet P. Nonuniform transformation field analysis. Int J Solids Struct 2003;40:6937–55. [36] Roussette S, Michel JC, Suquet P. Nonuniform transformation field analysis of elastic-viscoplastic composites. Compos Sci Technol 2009;69:22–7. [37] Fritzen F, Bohlke T. Nonuniform transformation field analysis of materials with morphological anisotropy. Compos Sci Technol 2011;71:433–42. [38] Willcox K, Peraire J. Balanced model reduction via the proper orthogonal decomposition. Aiaa J 2002;40:2323–30. [39] Schmid PJ. Dynamic mode decomposition of numerical and experimental data. J Fluid Mech 2010;656:5–28. [40] Berkooz G, Holmes P, Lumley JL. The proper orthogonal decomposition in the analysis of turbulent flows. Annu Rev Fluid Mech 1993;25:539–75. [41] Liu ZL, Bessa MA, Liu WK. Self-consistent clustering analysis: an efficient multiscale scheme for inelastic heterogeneous materials. Comput Methods Appl Mech Eng 2016;306:319–41. [42] Liu ZL, Fleming M, Liu WK. Microstructural material database for self-consistent clustering analysis of elastoplastic strain softening materials. Comput Methods Appl Mech Eng 2018;330:547–77. [43] Liu Z, Kafka OL, Yu C, Liu WK. Data-driven self-consistent clustering analysis of heterogeneous materials with crystal plasticity. Adv Comput Plast: Springer 2018:221–42. [44] Bessa MA, Bostanabad R, Liu Z, Hu A, Apley DW, Brinson C, et al. A framework for data-driven analysis of materials under uncertainty: countering the curse of dimensionality. Comput Methods Appl Mech Eng 2017;320:633–67. [45] Kanungo T, Mount DM, Netanyahu NS, Piatko CD, Silverman R, Wu AY. An efficient k-means clustering algorithm: analysis and implementation. IEEE Trans Pattern Anal Mach Intell 2002;24:881–92. [46] Buckley JD, Edie DD. Carbon-carbon Materials and Composites. William Andrew; 1993. [47] Bacarreza O, Abe D, Aliabadi M, Kopula Ragavan N. Micromechanical modeling of advanced composites. J Multiscale Modell 2012;4:1250005.

1996;56:381–6. [2] Mouritz AP, Bannister MK, Falzon PJ, Leong KH. Review of applications for advanced three-dimensional fibre textile composites. Compos Part a-Appl Sci Manuf 1999;30:1445–61. [3] Ryou H, Chung K, Yu WR. Constitutive modeling of woven composites considering asymmetric/anisotropic, rate dependent, and nonlinear behavior. Compos Part aAppl Sci Manuf 2007;38:2500–10. [4] Bohm R, Gude M, Hufenbach W. A phenomenologically based damage model for textile composites with crimped reinforcement. Compos Sci Technol 2010;70:81–7. [5] Skordos AA, Aceves CM, Sutcliffe MF. A simplified rate dependent model of forming and wrinkling of pre-impregnated woven composites. Compos Part a-Appl Sci Manuf 2007;38:1318–30. [6] Aboudi J, Arnold SM, Bednarcyk BA. Micromechanics of Composite Materials: A Generalized Multiscale Analysis Approach. Micromechanics of Composite Materials: A Generalized Multiscale Analysis Approach. 2013. p. 1–984. [7] Hou TY, Wu XH. A multiscale finite element method for elliptic problems in composite materials and porous media. J Comput Phys 1997;134:169–89. [8] Kanoute P, Boso DP, Chaboche JL, Schrefler BA. Multiscale methods for composites: a review. Arch Comput Methods Eng 2009;16:31–75. [9] Fish J. Multiscale modeling and simulation of composite materials and structures. Multiscale Methods Comput Mech: Prog Accomplishments 2011;55:215–31. [10] Kanit T, Forest S, Galliet I, Mounoury V, Jeulin D. Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int J Solids Struct 2003;40:3647–79. [11] Liu YJ, Chen XL. Evaluations of the effective material properties of carbon nanotube-based composites using a nanoscale representative volume element. Mech Mater 2003;35:69–81. [12] Sun CT, Vaidya RS. Prediction of composite properties, from a representative volume element. Compos Sci Technol 1996;56:171–9. [13] Chen XL, Liu YJ. Square representative volume elements for evaluating the effective material properties of carbon nanotube-based composites. Comput Mater Sci 2004;29:1–11. [14] Xia ZH, Zhang YF, Ellyin F. A unified periodical boundary conditions for representative volume elements of composites and applications. Int J Solids Struct 2003;40:1907–21. [15] Doitrand A, Fagiano C, Chiaruttini V, Leroy FH, Mavel A, Hirsekorn M. Experimental characterization and numerical modeling of damage at the mesoscopic scale of woven polymer matrix composites under quasi-static tensile loading. Compos Sci Technol 2015;119:1–11. [16] Daggumati S, Van Paepegem W, Degrieck J, Xu J, Lomov SV, Verpoest I. Local damage in a 5-harness satin weave composite under static tension: Part II – Meso-FE modelling. Compos Sci Technol 2010;70:1934–41. [17] Dai S, Cunningham PR. Multi-scale damage modelling of 3D woven composites under uni-axial tension. Compos Struct 2016;142:298–312. [18] Wang L, Wu JY, Chen CY, Zheng CX, Li B, Joshi SC, et al. Progressive failure analysis of 2D woven composites at the meso-micro scale. Compos Struct 2017;178:395–405. [19] Wang L, Zhao B, Wu J, Chen C, Zhou K. Experimental and numerical investigation on mechanical behaviors of woven fabric composites under off-axial loading. Int J Mech Sci 2018;141:157–67. [20] Patel DK, Waas AM, Yen CF. Direct numerical simulation of 3D woven textile composites subjected to tensile loading: an experimentally validated multiscale approach. Composites Part B-Eng 2018;152:102–15. [21] Huang T, Gong YH. A multiscale analysis for predicting the elastic properties of 3D woven composites containing void defects. Compos Struct 2018;185:401–10. [22] Miehe C, Schroder J, Becker M. Computational homogenization analysis in finite elasticity: material and structural instabilities on the micro- and macro-scales of periodic composites and their interaction. Comput Methods Appl Mech Eng 2002;191:4971–5005. [23] Geers MGD, Kouznetsova VG, Brekelmans WAM. Multi-scale computational

768