Global Sensitivity Analysis for the Elastic Properties of Unidirectional Carbon Fibre Reinforced Composites Based on Metamodels

Abstract

A fast and effective numerical method to predict mechanical properties of carbon fibre reinforced polymer (CFRP) composites, even elastic properties, is complicated due to the mismatch of mechanical properties among the constituents. Furthermore, it is not possible to completely characterise the influence of multiple parameters including mechanical and structural parameters on the bulk properties of CFRP by experiments. In this study, a three-phase finite-element model consisting of matrix, carbon fibre and interface was developed to predict the elastic mechanical behaviour of unidirectional CFRP. The elastic properties in terms of two Young's moduli, two Poisson's ratios and a shear modulus were calculated by means of a homogenisation method. High-accuracy Kriging surrogate models were constructed to fast-calculate the elastic responses for a large number of samples. Combining Kriging and high-dimensional model representation (HDMR) methods, a global sensitivity analysis was performed to study how the microscopic parameters influence the elastic responses to get a deeper understanding of elastic property-structure relationship. Eleven parameters, including mechanical and geometry properties of constituent phases, were chosen as inputs. Independent and cooperative effects of input parameters on the elastic properties of the studied composites were surveyed via first- and second-order sensitivity indices, respectively. An importance ranking of these parameters for each elastic response was derived directly by these indices. The procedure proposed in this work could serve as a theoretical guide for further design optimisation of CFRP.

1. Introduction

Carbon fibre-reinforced polymers (CFRP) are manufactured by compressing and curing polymeric matrix and carbon fibres under certain volume fractions. The general microstructure over the transverse cross-section of unidirectional CFRP is illustrated in Figure 1 ¹ . Fibre-reinforced composites have gained wide employments in industrial domains such as aerospace and automobile due to their outstanding specific strength, specific stiffness and excellent designability^2–4. Unidirectional CFRP is considered as a transversely isotropic material consisting of a matrix phase, a reinforcement phase and a transition region called the interphase between them, which presents high stiffness and strength along the fibre direction and was represented as the basic element in modelling laminates or fabrics. Evaluation of elastic properties is essential for exploring the applications of fibre-reinforced composites in industry.

Figure 1

Cross-section of unidirectional carbon fibre reinforced composites ¹

Many efforts have been made to predict the elastic properties of composites by linking microstructures together with mechanical parameters of constituents to the macro properties. Except for mechanical experiments which are of high cost, there exist two main kinds of quantitative methods to predict the mechanical properties of unidirectional fibre reinforced composites, namely the analytical method and the numerical method. Various micromechanical analytical models have been proposed, for instance, the Rule of Mixtures^5,6, Halpin-Tsai model ⁷ , Chamis model ⁸ and the Bridging model ⁹ . Analytical models could give direct mathematic formulations between elastic responses of composites and the parameters of constituents, which facilitate fast calculation of elastic constants. However, the mathematic expressions naturally limit in describing the details of stress and strain fields of microstructure. Besides, almost all analytical models were proposed on the basis that the composite is composed of only two phases: matrix and filler, which is in disagreement with the recognised fact that a transition zone of different properties exits between matrix and inclusions because of complex physical and chemical actions during manufacturing of composites ¹⁰ . Moreover, the assumption of isotropic property of fibre is incompatible with the property of carbon fibre.

Numerical methods, in particular, finite element method, were emerging as a widely used tool to study the mechanical behaviours of composites for its potential to overcome the disadvantages of analytical models^11–14. Within the framework of computational micromechanics, the transversely isotropic properties of carbon fibre can be modelled easily, and the microscopic stress and strain fields can be resolved, leading to the ability to better estimate not only the elastic properties but also damage and failure of CFRP ¹⁵ . Fibre-reinforced composites are heterogeneous on the microscale while they could be treated as a periodic array of a specific volume of material called Representative Volume Element (RVE). The homogenisation method could be applied to such materials to gain bulk properties by computing the representative unit when given the properties of constituents and microstructural parameters. Numerical modelling method has been recognised as a reliable tool, but the time consumed on the geometrical construction, model parameter definition and calculation represents a major drawback since the relation of results and input parameters of finite element simulation is not as explicit as mathematic formulations. In order to reduce the modelling and computational time and facilitate fast optimisation iterations when analysis of a composite material system is conducted, statistical methods such as surrogate models have been used, combined with the finite element simulations^16,17.

To better understand elastic responses of unidirectional CFRP, and quantify the influence of numerous input parameters, including microstructure parameters and mechanical properties of constituents, statistical methods such as sensitivity analysis are necessary. Existing research is limited to studying only the independent effects of a few microscopic parameters on the macromechanical properties of fibre-reinforced composites, mainly depending on analytical methods ¹⁸ , mechanical experiments^19,20, or finite element method via simply controlling a single variable^1,21. Various experiments have been developed to research the dependence of mechanical properties of unidirectional CFRP on the microscopic parameters. However, only partial parameters can be considered, and merely a trend of influence could be demonstrated by means of controlling variables during experiments. Meanwhile, there exist several analytical empirical formulae that could express the relationships between some microscopic parameters and elastic constants explicitly, but they cannot predict the elastic properties more precisely than numerical methods owing to several idealised hypotheses, limited parameters and inability to describe the detailed microscopic stress and strain fields. Neither experimental nor analytical methods can handle the elastic property-structure problem entirely in a quantitative way due to their inability to cover all parameters that may play definite roles in the mechanical behaviour of a unidirectional CFRP. Moreover, existing work on global sensitivity analysis for properties of composites relies on large amounts of finite element modelling and simulation ²² , along with relatively high computational costs. An efficient quantitative global sensitivity analysis to investigate both independent and cooperative effects of microscopic parameters on the elastic properties of unidirectional CFRP based on finite element modelling and high accuracy metamodels has not been reported.

In this work, three-dimensional unit cell finite element models composed of matrix, carbon fibre and interphase were developed to predict the transversely isotopic elastic properties of unidirectional CFRP in terms of moduli and Poisson's ratios. High-accuracy Kriging models of elastic response were then constructed to facilitate rapid calculation of responses for further sensitivity analysis instead of relatively expensive finite-element simulations. To parameterise the relationship between input parameters and elastic properties of composite in a quantitative manner, a global sensitivity analysis based on HDMR method using the computational results of constructed Kriging models was finally performed to characterise the independent and cooperative effects of microscopic input parameters on the elastic properties of unidirectional CFRP.

2. Finite-Element Simulation

The finite-element models to predict the elastic properties should first be developed instead of numerous experiments. Once the accuracy is verified, the models could be used to predict elastic constants under different combinations of input parameters.

2.1 Finite-element Modelling

It has been widely recognised that the unit cell homogenisation method, whose main idea is to compute a homogeneous medium equivalent to the macro heterogeneous composites, is accurate and efficient to calculate the mechanical properties of fibre-reinforced composites under elastic conditions. The general steps are as follows: firstly, an average stress or average deformation is applied to a detailed model (unit cell) on a given length scale in which the material is heterogeneous; secondly, equilibrium microscale stress and strain fields are computed by numerical methods such as finite-element simulations under specific loading conditions; finally, macroscopic responses are gained via a spatial average of microscale stress and strain fields. The periodicity of microscale stress and strain fields is characterised as:

$$\sigma \left(X+n{\lambda }_i{e}_i\right)=\sigma \left(X\right)$$

(1)

$$\epsilon \left(X+n{\lambda }_i{e}_i\right)=\epsilon \left(X\right)$$

(2)

λ _i is the length of unit cell, e _i is the unit vector, n is a positive integer.

The microscale stress and strain tensors are governed by constitutive equations of constituents below:

$$\sigma \left(X\right)=C\left(X\right)\epsilon \left(X\right)$$

(3)

$$\epsilon \left(X\right)=S\left(X\right)\sigma \left(X\right)$$

(4)

in which C(x) is the stiffness matrix of constituents and S(x) is the compliance matrix. For the linear elasticity of composites, the macroscopic constitutive formulation is also purely elastic and can be expressed as:

$${\sigma }^{*}=C^{*}{\epsilon }^{*}$$

(5)

$${\epsilon }^{*}=S^{*}{\sigma }^{*}$$

(6)

C^* is the effective stiffness tensor and S^* is the effective compliance tensor. σ^* and ε^* are the effective stress and strain tensors, respectively, which are calculated by the volume averaging over the RVE and described as:

$${\sigma }^{*}=\frac{1}{V}{\displaystyle \underset{V}{\int }\sigma (X)dV}$$

(7)

$${\epsilon }^{*}=\frac{1}{V}{\displaystyle \underset{V}{\int }\epsilon (X)dV}$$

(8)

where V is the volume of the RVE.

The stress and strain tensors are calculated through applying appropriate boundary conditions under given macroscopic loads and then the elastic properties of composites can be calculated by using Equations (7) and (8).

For unidirectional CFRP, the fibres are often distributed randomly over the cross-section as shown in Figure 1 . Thus, a RVE which can represent the composite substantially would contain several fibres located randomly in the matrix when RVEs are computed to calculate the mechanical properties of composite. However, Trias et al. ²³ pointed that periodic models could be used in the case of computing elastic constants while random models should be used for failure and damage simulations. Figure 2 illustrates the idealised periodic microstructure array of unidirectional CFRP and a 3D unit cell consisted of one fibre with diameter D _f centred in a cubic matrix of length size L. Then the volume fraction of fibre is expressed as:

$$V_f=\frac{\pi }{4}{\left(\frac{D_f}{L}\right)}^2$$

(9)

Figure 2

Idealised microstructure of unidirectional carbon fibre reinforced composite and a 3D unit cell

In order to take into consideration the interface, which is recognised as a transition zone between fibre and matrix, a microstructure-based three-phase unit cell model composed of fibre, matrix and interface was developed. The epoxy resin matrix was presumed to be an isotropic elastic material and the carbon fibre was modelled as a transversely isotropic material. An eight-node cohesive element was adopted in the unit cell model to simulate the mechanical behaviour of the interface, as shown in Figure 3 .

Figure 3

Description of 8-node cohesive element

Deformation of a cohesive element includes relative normal motions and relative tangential displacements between top face and bottom face, which measure the normal traction due to opening or closing and transverse shear deformation of the interface, respectively. The mechanical behaviour of cohesive elements is governed by a traction-separation law of various forms, e.g. bilinear, linear-parabolic, exponential, trapezoidal ²⁴ . A typical bilinear traction-separation law was selected in this study, as illustrated in Figure 4 .

Figure 4

Bilinear traction-separation law for cohesive element

The nominal traction stress vector t is composed of three components: t _n , t _s , and t _t , which represent the normal and two shear tractions, respectively. The corresponding separations are denoted by δ _n , δ _s , and δ _t . The nominal strains can be defined as:

$${\epsilon }_n=\frac{{\delta }_n}{T_i},{\epsilon }_S=\frac{{\delta }_s}{T_i},{\epsilon }_t=\frac{{\delta }_t}{T_i}$$

(10)

in which T _i denotes the original thickness of the cohesive element. The elastic behaviour can be written as:

$$t=\left\{\begin{array}{c}{t}_n\\ t_s\\ t_t\end{array}\right\}=\left[\begin{array}{ccc}{K}_{nn}& K_{ns}& K_{nt}\\ K_{ns}& K_{ss}& K_{st}\\ K_{nt}& K_{st}& K_{tt}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_n\\ {\epsilon }_s\\ {\epsilon }_t\end{array}\right\}=K\epsilon$$

(11)

For simplicity of computation, an uncoupled relationship of normal and tangential stiffness coefficients was introduced and two coefficients in tangent directions were assumed to be equal: K _ss = K _tt . The properties of the interface may vary along its thickness ²⁵ . To simplify the models, an effective interface where elastic properties are constant through thickness was assumed, which is recognised to be reasonable when applied to fibre-reinforced composites ²⁶ . Thus, the properties of the interface using a cohesive element could be determined by three parameters, namely, K _in , K _it , and T _i , which represent the normal stiffness, the tangent stiffness and the thickness of the interface, respectively.

Figure 5 illustrates the established three-phase cubic unit cell finite element model for unidirectional CFRP. The unit cell was meshed by 8-node hexahedral elements. The element types of matrix and fibre were selected as C3D8I, and the interface as COH3D8 using the ABAQUS preprocessor.

Figure 5

Three-phase finite element unit cell model

In order to fulfil the periodical characteristics of composites, the unified periodic boundary conditions ²⁷ defined in terms of displacements, which satisfies not only the boundary displacement periodicity but also boundary traction periodicity when using displacement-based finite element simulation, should be applied to the unit cell model. The unified periodic boundary conditions are expressed as follows:

$$\overrightarrow{u}(0,X_2,X_3)-\overrightarrow{u}(L_1,X_2,X_3)={\overrightarrow{U}}_1$$

(12)

$$\overrightarrow{u}(X_1,0,X_3)-\overrightarrow{u}(X_1,L_2,X_3)={\overrightarrow{U}}_2$$

(13)

$$\overrightarrow{u}(X_1,X_2,0)-\overrightarrow{u}(X_1,X_2,L_3)={\overrightarrow{U}}_3$$

(14)

${\overrightarrow{U}}_1$ , ${\overrightarrow{U}}_2$ and ${\overrightarrow{U}}_3$ denote displacement vectors of the opposite faces, L₁, L ₂ and L ₃ are the lengths of the unit cell along three orthogonal directions, respectively.

The unified periodic boundary conditions above were realised by applying displacement boundary conditions to the unit cell model. For the considerations of symmetry characteristic and reducing computing time, one-eighth of the cubic unit cell model was used when calculating the Young's modulus and Poisson's ratio, as shown in Figure 6a . The boundary conditions applied to the finite-element model are as follows: (i) Symmetric boundary conditions were applied to the three symmetric planes; (ii) One plane on the tensile direction was used as the loading plane on which uniform normal displacement was applied to all the nodes of the plane. Among the three planes, the plane parallel to the fibre cross-section was chosen as the loading plane when calculating longitudinal properties, while the plane perpendicular to the fibre cross-section was selected to calculate transverse responses; (iii) As for the other two planes, nodes were coupled in the normal direction to assure the periodicity.

Figure 6

Models for different loading conditions: (a) normal loading (b) shear loading

As for calculating shear modulus, half of the unit cell model was modelled as shown in Figure 6b . The boundary conditions were applied as: (i) Uniform tangent displacement along fibre direction when calculating in-plane shear modulus is applied to the shear plane that is vertical to fibre cross section; (ii) To facilitate the periodic boundary conditions, the mesh in opposite boundary planes was nearly the same. A constraint equation was imposed on each pair of displacement components at the two corresponding nodes of the pair opposite boundary planes with the same in-plane coordinates; (iii) Rigid motions of the nodes on the plane opposite to the shear plane on the direction of shear loading were removed, and the nodes on the lines vertical to shear direction on these two opposite planes were applied to fix the rigid body degree in the direction perpendicular to the shear plane.

2.2 Results of Numerical Simulations

Under the unified periodic boundary conditions, the macroscopic strain in Equation (8) can be calculated as:

$${\epsilon }_{ii}^{*}=\frac{{\overrightarrow{U}}_i}{L_i},{\epsilon }_{ij}^{*}=\frac{{\overrightarrow{U}}_i}{L_j},\left(i\ne j\right)$$

(15)

the subscript i denotes the direction of displacement loading. The macroscopic stress in Equation (7) is calculated as:

$${\sigma }_{ij}^{*}=\frac{{(P_i)}_j}{S_j}$$

(16)

S _j is the area of jth boundary surface on which the displacement is applied, (P _i ) _j is the traction force on the loading surface. Thus the overall elastic responses of the composite in terms of Young's modulus, Poisson's ratio and shear modulus can be calculated as:

$$E_{ii}=\frac{{\sigma }_{ii}^{*}}{{\epsilon }_{ii}^{*}},{\nu }_{ij}=-\frac{{\epsilon }_{jj}^{*}}{{\epsilon }_{ii}^{*}},G_{ij}=\frac{{\sigma }_{ij}^{*}}{{\epsilon }_{ij}^{*}}$$

(17)

The elastic responses of the unit cell models subject to tensile and shear loads were simulated by the FE code ABAQUS/Standard. The macroscopic stress and strain were calculated through Equations (16) and (15) after recording the required traction force during simulation. Hence, the elastic constants of the studied composite were obtained, based on Equation (17). The deformation mode and stress distribution showed reasonable patterns under each load conditions, which in turn ensured the mesh convergence of the finite element models.

Based on the developed finite element models and the material parameters of constituents as listed in Table 1 ²⁸ and Table 2 ²⁹ , in which the matrix was isotropic and the carbon fibre was modelled as transversely isotropic, elastic responses were calculated using finite element simulating results. The results, as illustrated in Table 3 , demonstrated that the unit cell finite element simulation method was of high accuracy to predict five individual elastic properties of unidirectional CFRP, namely longitudinal Young's modulus E₁₁, major Poisson's ratio v₁₂, transverse Young's modulus E₂₂, transverse Poisson's ratio v₂₃ and in-plane shear modulus G₁₂.

Table 1

Mechanical properties of epoxy resin matrix²⁸

Mechanical properties	BSL914C epoxy
Modulus, E m (GPa)	4.0
Poisson's ratio, v m	0.35

Table 2

Mechanical properties of carbon fibre²⁹

Mechanical properties	Carbon fibre (T300)
Longitudinal modulus, E f 11 (GPa)	230
Transverse modulus, E f 22 (GPa)	18.6
Major Poisson's ratio, v f 12	0.255
In-plane shear modulus, G f 12 (GPa)	20.5
Transverse shear modulus, G f 23 (GPa)	7

Table 3

Comparison of finite element predicted results and experimental data

Mechanical properties	Experiment results28	Predicted results	Relative error (%)
Longitudinal modulus, E11 (GPa)	138	139.4	1.01
Major Poisson's ratio, v12	0.28	0.281	0.36
Transverse modulus, E22 (GPa)	11	10.9	0.91
Transverse Poisson's ratio, v23	0.4	0.402	0.50
In-plane shear modulus, G12 (GPa)	5.5	5.4	1.82

3. Construction of Kriging Metamodels

Metamodelling technique has been widely used to construct surrogate models for the purpose of parameterising the relationship between inputs and outputs and avoiding adding more relatively expensive finite element simulations. Once the approximate mathematic relationship based on sufficient experimental or computational simulation data was established, metamodels could be used to predict the outputs at given untried inputs immediately. Then the predictions could be adopted for further use of computation or optimisation, and so on, as long as the predicted accuracy was verified to be high enough.

3.1 Kriging Metamodels

Kriging metamodel is an unbiased estimate model within which sample points are interpolated to predict the values of unknown points and estimate the trend of stochastic process used in metamodelling by Gaussian random function. Kriging metamodel was adopted as the surrogate model for its high global approximate accuracy. The basic formulation of Kriging model is composed of a global model of polynomial regression and a stochastic term:

$$y(x)=f^T(x)\beta +Z(x)$$

(18)

y(x) is the unknown function of interest, f ^T (x)β is the known approximate function as a linear regression of the data, within which f (x) is the basis functions and (3 is the corresponding coefficient. Z(x) is a stochastic process to represent the departure from the linear model with mean 0, variance σ ² , and nonzero covariance. The covariance matrix of Z(x) is expressed as:

$$Cov[Z(x^i),Z(x^j)]={\sigma }^{2}R;i,j=1,2,\dots ,n_s$$

(19)

where n _s is the number of sample points, R stands for the correlation matrix. Gaussian function is chosen as the correlation function, then the general form of R can be expressed as:

$$R(i,j)=R(x^i,x^j)={\Pi }_{k=1}^{n_d}\exp \left(-{\theta }_k|x_k^i-X_k^j{|}^{p_k}\right)$$

(20)

R(x ⁱ , x ^j ) (i, j=1,2…,n _s ) represents the correlation function between sample x ⁱ and x ^j which is related to the distance between the two points; n _d denotes the dimension number of the design space; θ _k ≥ 0 and 0 ≤ p _k ≤ 2 are two parameters that can be determined by the Maximum Likelihood Estimates (MLE) method:

$$\begin{array}{c}\text{Maximise}:-[n_{s}In\left({\widehat{\sigma }}^2\right)+In|R|]/2\\ \text{Subject to}:\{\begin{array}{c}{\theta }_k\ge 0,\\ 0\le p_k\le 2,\end{array}k=1,2,\dots ,n_d\end{array}$$

(21)

where ${\widehat{\sigma }}^2$ is a function of θ _k and p _k :

$${\widehat{\sigma }}^2={\left(Y-\widehat{\beta }F\right)}^T{R}^{-1}\left(Y-\widehat{\beta }F\right)/n_s$$

(22)

$$\widehat{\beta }={\left(F^T{R}^{-1}F\right)}^{-1}\left(F^T{R}^{-1}Y\right)$$

(23)

where F is a unit matrix. Then the predicted value of any unknown point x is obtained as:

$$\widehat{Y}\left(x\right)=f^t\left(x\right)\beta +r^T\left(x\right)R^{-1}\left(Y-F^T\beta \right)$$

(24)

3.2 Metamodels Construction of Elastic Properties

In this work, eleven parameters, including geometric and mechanical parameters, were selected as input parameters as listed in Table 4 , where V _f is the volume fraction of fibre, and E _m and v _m denote the Young's modulus and Poisson's ratio of the matrix, respectively. E ^f ₁₁, E ^f ₂₂ and G ^f ₁₂ represent longitudinal, transverse Young's modulus and shear modulus of carbon fibre, respectively. v ^f ₁₂ and v ^f ₂₃ are major and transverse Poisson's ratio of carbon fibre. K _in , K _it and T _i denote the normal, tangent interfacial stiffness and thickness of interface, respectively.

Table 4

Ranges of input parameters

Input parameter	Lower bound	Upper bound	Unit
V f	20	65	%
E m	1	7	GPa
vm	0.1	0.46	-
E f 11	100	360	GPa
E f 22	5	25	GPa
G f 12	5	25	GPa
v f 12	0.1	0.46	-
v f 23	0.01	0.46	-
K in	1	1000	GPa•mm−1
K it	1	1000	GPa•mm−1
T i	0.05	0.4	μm

Sufficient data of different set of input parameters and corresponding target responses are needed to construct the metamodels, and several additional data should be used to evaluate the accuracy of the constructed models. In this work, 150 training sample points and 30 test points are generated by Optimal Latin Hypercube Design (OLHD) ³⁰ . Following the lower bounds and upper bounds of input parameters, the OLHD sampling method could generate sample points uniformly filled in the entire design space most. Corresponding finite element model for each set of input parameters was established following the modelling method developed in section 2, five individual elastic responses of each sample were calculated subsequently based on simulation results operated on parallel computing facilities. A Gaussian function was used as the correlation function and the surrogate models of five outputs were constructed through the Kriging metamodelling technique. To verify the accuracy of these metamodels, coefficient of determination (R ² ) criterion as expressed below was selected:

$$R^2=\frac{{\displaystyle \sum _{i=1}^{n_{test}}{\left(y_i-{\widehat{y}}_i\right)}^2}}{{\displaystyle \sum _{i=1}^{n_{test}}{\left(y_i-{\bar{y}}_i\right)}^2}}$$

(25)

where y _i is the calculated value of origin model, ŷ _i is the predicted value through metamodels, ȳ _i represents the mean value of simulated results of origin model, n _test is the number of test samples.

The values of R ² are listed as Table 5 . The Kriging metamodels of five outputs are of specifically high accuracy, which may result from the sufficient number of training samples and relatively weak nonlinear of elastic properties. Thus, the surrogate models are available for further fast computing procedures.

Table 5

Coefficient of determination (R²) of Kriging metamodels

Response	E 11	v 12	E 22	v 23	G12
R 2	0.9999	0.9991	0.9975	0.9879	0.9994

4. Global Sensitivity Analysis using Kriging and Rs-Hdmr

To better understand the bulk properties of unidirectional CFRP and fully characterise the influence of various microscopic parameters, a global sensitivity analysis based on RS-HDMR is performed. High dimensional model representation (HDMR) ³¹ is a set of quantitative model assessment method for expressing the input-output relationships of a complex system with a high dimensional input space. HDMR allows the model responses to be captured by the first few lower-order terms as higher-order correlated behaviour of input parameters is expected to be weak for most well-defined systems. HDMR reveals the hierarchy of correlations among input variables for each component function represents a unique contribution of the input parameters to dependently or cooperatively influence the output response. The main purpose of sensitivity analysis is to investigate either isolated or combined effects of inputs on the outputs of model, thereby the main contributors were identified. HDMR presents a model replacement that can be conveniently employed within sensitivity analysis, while relatively large numbers of samples are needed to achieve higher accuracy. In this work, the HDMR method combined with the Kriging models established above were employed to perform the global sensitivity analysis for the elastic responses of unidirectional CFRP.

4.1 Theoretical Framework of HDMR

As the impact of multiple input variables on the output can be independent and cooperative, the mapping between the input variables x₁, …, x _n and the output variables f(x) = f(x₁, …, x _n ) in the domain Rⁿ can be expressed in the following form³¹:

$$\begin{array}{c}f\left(x\right)=f_0+{\displaystyle \sum _{i=1}^n{f}_i\left(X_i\right)+{\displaystyle \sum _{1\le i<j\le n}{f}_{ij}\left(X_i,X_j\right)+\dots }}\\ +f_{12\cdots n}\left(X_1,X_2,\cdots ,X_n\right)\end{array}$$

(26)

where f₀ is a constant that denotes the zeroth-order effect in the domain. The first-order term f _i (x _i ) gives the effect of the input variable x _i acting independently on the output f(x). The function f _ij (x _i , x _j ) is a second-order term describing the joint effects of the variables x _i and x _j upon the output f(x). Higher-order terms represent the cooperative effects of increasing number of variables acting together on f(x). If there is no interaction between the input variables, only zeroth-order component f₀ and first-order term f _i (x _i ) will exist in HDMR expansions. These component functions are thus hierarchically tailored to f(x) over the entire domain Rⁿ.

Experience has shown that even a low order HDMR expansion provides already a sufficient description of f(x) in real problems ³¹ . This motivates the usage of the second-order expansion in this work:

$$f\left(x\right)=f_0+{\displaystyle \sum _{i=1}^n{f}_i\left(X_i\right)+{\displaystyle \sum _{1\le i<j\le n}{f}_{ij}\left(X_i,X_j\right)}}$$

(27)

Such a function approximation facilitates not only data reduction but also the decrease of computational complexity.

There are two commonly used HDMR expansions ³¹ : Cut-HDMR, which depends on the value of the value of f(x) at a specific reference point, and RS-HDMR, which depends on the average value of f(x) over the whole domain. The Cut-HDMR expansions can be constructed when ordered sampling of f(x) is possible, while the RS-HDMR component functions can be determined with only one set of random samples. In this work, RS-HDMR is applied to the global sensitivity analysis.

To construct RS-HDMR expansion, a set of random sampling points over the entire domain Rⁿ is used to determine the component functions via averaging process. The input variable x _i is rescaled to the range [0,1] for all i, then the output function is defined in the unit hypercube K ⁿ ={(x₁,…,x _n ),i=1,…,n}. The component functions of RS-HDMR are expressed as follows:

$$f_0={\displaystyle \underset{k^n}{\int }f\left(x\right)dx}$$

(28)

$$f_i\left(x_i\right)={\displaystyle \underset{k^{n-1}}{\int }f\left(x\right)d{x}^i-f_0}$$

(29)

$$f_{ij}\left(X_i,X_j\right)={\displaystyle \underset{k^{n-2}}{\int }f\left(x\right)d{x}^{ij}-f_i\left(X_i\right)-f_j\left(X_j\right)-f_0}$$

(30)

where dx ⁱ denotes the product dx₁dx₂…dx _n without dx _i , dx ^ij stands for the same product without dx _i and dx _i . The constant term f₀ can be approximated by averaging the value of f(x) for all x^(s) = (x₁^(s),x₂^(s),…,x _n ^(s))(s = 1,2,…, N):

$$f_0={\displaystyle \underset{k^n}{\int }f\left(x\right)dx\approx \frac{1}{N}{\displaystyle \sum _{s=1}^{N}f\left(x^{\left(S\right)}\right)}}$$

(31)

The higher-order component functions can be evaluated via direct Monte Carlo integrals which is prohibitively expensive due to large number of full model runs. To reduce the sampling effort and facilitate the calculation of sensitivity indices in the next section, the component functions can be approximated by analytical basis functions (orthonormal polynomials, for instance) ³² :

$$f_i\left(X_i\right)\approx {\displaystyle \sum _{r=1}^k{\alpha }_r^i{\phi }_r\left(x_i\right)}$$

(32)

$$f_{ij}(X_i,X_j)\approx {\displaystyle \sum _{p=1}^I{\displaystyle \sum _{q=1}^{I^{\prime }}{\beta }_{pq}^{ij}{\phi }_p(X_i){\phi }_q(X_j)}}$$

(33)

where k,l,l’ denote the orders of the polynomial expansion, φ _r (x _i ), φ _p (x _i ) and φ _q (x _i ) are the orthonormal basis functions, α ⁱ _r and β ^ij _pq are coefficients to be determined by a minimisation and integration ³² due to the orthonormal property of basis functions, which can be calculated as:

$${\alpha }_r^i\approx \frac{1}{N}{\displaystyle \sum _{S=1}^{N}f\left(x^{\left(S\right)}\right){\phi }_r\left(X_I^{\left(S\right)}\right)}$$

(34)

$${\beta }_{pq}^{ij}\approx \frac{1}{N}{\displaystyle \sum _{S=1}^{N}f\left(x^{\left(S\right)}\right){\phi }_p\left(X_I^{\left(S\right)}\right)}{\phi }_q\left(X_j^{\left(S\right)}\right)$$

(35)

4.2 Construction of RS-HDMR for Unidirectional CFRP

The same eleven parameters as listed in Table 4 were selected as input parameters. Five individual elastic responses of unidirectional CFRP in terms of longitudinal Young's modulus E₁₁, major Poisson's ratio v₁₂, transverse Young's modulus E₂₂, transverse Poisson's ratio v₂₃ and in-plane shear modulus G₁₂ were chosen as outputs.

As for the sampling process, five sets of training samples with different sizes N = (200, 300, 400, 500, 600) and a set of test samples with size N=100 were generated using OLHD ³⁰ to get the samples that uniformly filled the entire design space to the best. All the input variables were rescaled by mathematical transformations such that x _i ∈[0,1]. The outputs of each sample point were calculated by high-accuracy Kriging metamodels developed in Section 3.2 rapidly, which dramatically spare the computational efforts compared to direct FE simulations to facilitate specifically high efficiency of sensitivity analysis. RS-HDMR expansions up to second-order of each output were adopted. A set of orthonormal polynomials ³² shown as follows were chosen to approximate the component functions:

$${\phi }_1\left(X\right)=\sqrt{3}(2x-1)$$

(36)

$${\phi }_2\left(X\right)=6\sqrt{5}(x^2-x+\frac{1}{6})$$

(37)

$${\phi }_3\left(X\right)=20\sqrt{7}(x^2-\frac{3}{2}{x}^2+\frac{3}{5}\times -\frac{1}{20})$$

(38)

The best polynomial order that approximating each component function is determined by an optimisation algorithm developed by T. Ziehn et al ³³ based on a least-squares method. Sobol's method ³⁴ is widely used in global sensitivity analysis, and RS-HDMR is convenient to imbed into sensitivity analysis as it is conceptually the same as Sobol's method. The total variance D is given as:

$$D={\displaystyle \underset{K^n}{\int }{f}^2\left(x\right)dx-f_0^2}$$

(39)

The partial variances are given by:

$$D_i={\displaystyle {\int }_0^1{f}_i^2\left(X_i\right)d{x}_i}$$

(40)

$$D_{ij}={\displaystyle {\int }_0^1{\displaystyle {\int }_0^1{f}_{ij}^2\left(X_i,X_j\right)d{x}_{i}d{x}_j}}$$

(41)

Considering the orthonormal property of approximated component functions, the partial variances can be calculated as follows ³⁵ :

$$D_i\approx {\displaystyle \sum _{r=1}^{k_i}{\left({\alpha }_r^i\right)}^2}$$

(42)

$$D_{ij}\approx {\displaystyle \sum }_{p=1}^{l_i}{\displaystyle \sum }_{q=1}^{{l^{\prime }}_j}{({\beta }_{pq}^{ij})}^2$$

(43)

Finally, the global sensitivity indices can be defined as:

$$S_{i_1,\cdots ,i_S}=\frac{D_{i_1,\cdots ,i_s}}{D}$$

(44)

The sensitivity indices can then be used for importance ranking of input variables on each output.

5 Results and Discussion

The overall accuracy of constructed RS-HDMR expansions was also tested by the coefficient of determination (R ² ) as shown in Equation (25). Values of R ² for each elastic output response with different sample sizes are listed in Table 6 . The second-order RS-HDMR models present a higher accuracy than first-order models as they include both independent and cooperative effects of input parameters. Generally, with an increase of sample size, the values of R ² are on the rising trend and converge to 1. The relatively small set of samples of size N = 600 is sufficient to accurately estimate the elastic responses using second-order RS-HDMR expansions. The accuracy of the second-order RS-HDMR expansion with the sample size N = 600 can also be illustrated intuitively by the scatter plots as shown in Figure 7 , which demonstrates the relationship between calculated Kriging metamodel values and the predicted values via RS-HDMR. Meanwhile it can be assumed that there are no crucial higher-order effects of inputs due to the sufficiently high accuracy.

Figure 7

Scatter plots for the calculated and predicted elastic constants. (a) longitudinal Young's modulus. (b) major Poisson's ratio. (c) transverse Young's modulus. (d) transverse Poisson's ratio. (e) in-plane shear modulus

Table 6

Coefficient of determination (R²) of RS-HDMR expansions under different sample sizes

Sample size (N)		200	300	400	500	600
R2 (E11) (%)	1st-order	89.96	93.95	94.21	94.31	94.43
	2nd-order	93.04	94.89	99.23	99.81	99.87
R2 (v12) (%)	1st-order	91.89	93.09	93.36	93.46	93.61
	2nd-order	92.94	95.72	96.36	98.84	99.02
R2 (E22) (%)	1st-order	88.21	90.27	89.03	90.53	91.29
	2nd-order	92.34	96.52	96.12	96.33	97.37
R2 (v23) (%)	1st-order	88.62	89.69	90.82	91.36	91.88
	2nd-order	95.34	96.38	97.73	97.67	98.51
R2 (G12) (%)	1st-order	93.12	93.41	93.50	93.63	94.09
	2nd-order	94.63	95.92	97.70	97.92	98.48

The calculated global sensitivity indices for five outputs with sample size N=600 are illustrated in Tables 7–11 which give the first five first- and second-order sensitivity indices, namely S _i and S _ij , and their rankings, in which S _i and S _ij represent the degree of individual and coupled influence of input parameters on the output elastic properties, respectively, the higher the index, the greater the effect on the outputs.

Table 7

First- and second-order sensitivity indices for longitudinal Young's modulus

Output	E 11
Ranking	Input	S i	Input	S ij
1	E f 11	0.49786	V f & E f 11	0.04410
2	V f	0.41257	G f 12 & v f 12	1.41E-05
3	E m	0.00116	E f 22 & T i	1.23E-05
4	K it	0.00011	E f 22 & K it	7.53E-06
5	E f 22	0.00011	V f & v f 23	5.06E-06
ΣS i (ΣS ij )		0.9121		0.0441
ΣS i + ΣS ij	0.9562

Table 8

First- and second-order sensitivity indices for major Poisson's ratio

Output	v 12
Ranking	Input	S i	Input	S ij
1	v f 12	0.47029	V f & v m	0.03539
2	m	0.44825	V f & v f 12	0.03344
3	K in	0.00064	E m & v f 12	0.00044
4	K it	0.00061	E f 22 & v f 12	0.00020
5	E f 22	0.00050	E f 22 & K in	0.00006
ΣS i (ΣS ij )		0.9202		0.0696
ΣS i + ΣS ij	0.9898

Table 9

First- and second-order sensitivity indices for transverse Young's modulus

Output	E 22
Ranking	Input	S i	Input	S ij
1	E m	0.49206	V f & E f 22	0.03874
2	V f	0.17893	E m & E f 22	0.03080
3	E f 22	0.17377	V f & v m	0.00100
4	v m	0.03947	G f 12 & K in	0.00014
5	v f 12	0.00069	v f 23 & T i	0.00012
ΣS i (ΣS ij )		0.8855		0.071
ΣS i + ΣS ij	0.9565

Table 10

First- and second-order sensitivity indices for transverse Poisson's ratio

Output	v 23
Ranking	Input	S i	Input	S ij
1	v m	0.78816	V f & v m	0.03152
2	v f 23	0.08468	V f & v f 23	0.01438
3	V f	0.05590	E m & v f 23	0.00213
4	K in	0.00117	E f 22 & v f 23	0.00212
5	E f 22	0.00116	E m & v m	0.00204
ΣS i (ΣS ij )		0.9322		0.0533
ΣS i + ΣS ij	0.9855

Table 11

First- and second-order sensitivity indices for in-plane shear modulus

Output	G 12
Ranking	Input	S i	Input	S ij
1	E m	0.49143	V f & G f 12	0.01938
2	V f	0.33418	V f & E m	0.01499
3	G f 12	0.05247	E m & G f 12	0.01429
4	vm	0.01402	E m & v m	0.00117
5	v f 12	0.00032	K it & T i	0.00005
ΣS i (ΣS ij )		0.8928		0.0501
ΣS i + ΣS ij	0.9429

The sensitivity indices were ranked to show which input parameter or combination of two input parameters contribute most to the overall variance, acting as the criteria that determine the relative importance of each parameter. According to the overall values of first- and second-order sensitivity indices for each elastic response of composite, all five elastic properties are affected by input parameters acting both independently and interactively while with different rankings. Almost the independent influences of input parameters account for most of the variance; meanwhile, the interactions of parameters also play a relatively important role. For each elastic property, all first- and second-order sensitivity indices add up to nearly 0.95 or above, indicating that the majority of the total variance has been covered, which further proves the accuracy and good approximation of second-order RS-HDMR for predicting the elastic properties of the studied composites.

The optimal orders of orthonormal polynomials used to approximate the first five first-order RS-HDMR component functions that influence corresponding output the most are listed in Table 12 . The highest order of polynomials adopted was 5. Only the orders of polynomials approximating first-order component function of V _f and T _i with regard to the output E₁₁, and G ^f ₁₂ with respect to G₁₂ are equal to 1, which means the first-order responses to input V _f and T _i are linear for the output E₁₁, and the same as input G ^f ₁₂ for output G₁₂. All the other first-order component functions are represented by second- or higher-order polynomials, which indicates that the first-order response to each set of parameters is nonlinear for their corresponding elastic outputs to some degree. Furthermore, the component functions of RS-HDMR expansions could be plotted to illustrate the relationship between one certain input parameter with the interested elastic output. The first-order component functions of five elastic properties for the two or three most important input parameters are plotted as Figure 8–12 respectively. The curves demonstrate the degree of linearity between a certain input parameter and a specific elastic constant, the slope of which can reflect the rate of change of elastic response along with the variation of input parameter to a certain extent.

Figure 8

First-order component functions of longitudinal Young's modulus for the two most important parameters: (a) longitudinal Young's modulus of carbon fibre. (b) volume fraction of carbon fibre

Figure 9

First-order component functions of major Poisson's ratio for the two most important parameters: (a) major Poisson's ratio of carbon fibre. (b) Poisson's ratio of matrix

Figure 10

First-order component functions of transverse Young's modulus for the three most important parameters: (a) Young's modulus of matrix. (b) volume fraction of carbon fibre. (c) transverse Young's modulus of carbon fibre

Figure 11

First-order component functions of transverse Poisson's ratio for the three most important parameters: (a) Poisson's ratio of matrix. (b) transverse Poisson's ratio of carbon fibre. (c) volume fraction of carbon fibre

Figure 12

First-order component functions of in-plane shear modulus for the three most important parameters: (a) Young's modulus of matrix. (b) volume fraction of carbon fibre. (c) in-plane shear modulus of carbon fibre

Table 12

Order of polynomials used for the approximation of first-order RS-HDMR component functions of elastic models using a sample size N=600

Output	Polynomial orders of first-order component functions
	f1(V f )	f2(E m )	f3(v m )	f4(E f 11)	f5(E f 22)	f6(G f 12)	f7(v f 12)	f8(v f 23)	f9(K in )	f10(K it )	f11(T i )
E 11	1	5		5			3				1
v 12			3		5		3		4	5
E 22	5	5	2		3		3				5
v 23	3		4		2			4	4		3
G 12	3	2	4			1	5				3

As shown in Table 7 , for the longitudinal Young's modulus of composite E₁₁, longitudinal Young's modulus of carbon fibre E ^f ₁₁ and its volume fraction V _f as the most important parameters, the individual contributions of E ^f ₁₁ and V _f as well as their joint contribution to the overall variance of E₁₁ are 49.8%, 41.3%, and 4.4%, respectively. The properties of matrix and other microstructure parameters only have very slight influence on E₁₁. In Figure 8 , a pure linear relationship between E ^f ₁₁ and E₁₁ can be seen; meanwhile, the relationship between V _f and E₁₁ is almost linear.

For major Poisson's ratio v₁₂, major Poisson's ratio of carbon fibre v ^f ₁₂ and Poisson's ratio of matrix v _m contribute most to the variance of v₁₂ by the close percentages of 47% and 44.8% respectively as listed in Table 8 . The volume fraction V _f does not affect v₁₂ independently, but takes effect with these two most important individual parameters v ^f ₁₂ and v _m cooperatively with the second-order sensitivity indices 3.5% and 3.3% respectively. Other than these three parameters, the rest nearly have no influence on v₁₂. v₁₂ changes nearly linearly with v ^f ₁₂ as shown in Figure 9a , while the relationship between v₁₂ and v _m is slightly nonlinear, as illustrated in Figure 9b .

As for the transverse Young's modulus E₂₂ shown in Table 9 , the first four significant input parameters affecting E₂₂ independently are ordered as: Young's modulus of matrix E _m , volume fraction V _f , transverse Young's modulus of carbon fibre E ^f ₂₂, Poisson's ratio of matrix v _m . E _m influences E₂₂ more obviously strongly than the other parameters which contributes 49.2% to the overall variance of E₂₂. Among the first three important parameters E _m , V _f , and E ^f ₂₂, cooperative effects of E ^f ₂₂ together with V _f and E _m rank in front of second-order indices. The relationships between E₂₂ and E _m , E ^f ₂₂ are nonlinear with monotonically-decreasing slopes, as shown in Figure 10a, 10c; In Figure 10b , a highly nonlinear relationship between V_f and E₂₂ is presented.

As shown in Table 10 , the transverse Poisson's ratio of composite v₂₃ heavily depends on the Poisson's ratio of matrix v _m which contributes specifically high percentage to the overall variance of v₂₃ with 78.8%. v₂₃ is further influenced by the transverse Poisson's ratio of carbon fibre v ^f ₂₃ and volume fraction V _f . Besides, the coupled parameter pairs: V _f with v _m and V _f with v ^f ₂₃, take effect cooperatively on v₂₃. v₂₃ changes nonlinearly with v _m in a monotonic increasing slope as shown in Figure 11a ; Obvious nonlinear relationships between v₂₃ and v ^f ₂₃, V _f are presented in Figures 11b and 11c.

As illustrated in Table 11 , for the in-plane shear modulus of composite G₁₂, Young's modulus of matrix E _m and volume fraction V _f are the most important parameters, with contributions of 49.1% and 33.4% to the overall variance respectively. The in-plane shear modulus of carbon fibre G ^f ₁₂ and Poisson's ratio of matrix v _m rank in second place influencing G₁₂. Cooperative effects exist in the first three important individual parameters are listed in Table 11 . As shown in Figure 12a , G₁₂ changes slightly nonlinearly with E _m in a monotonic decreasing slope; In contrast, with V _f , G₁₂ changes with a monotonically-increasing slope, as illustrated in Figure 12b . Meanwhile, Figure 12c shows a pure linear relationship between G₁₂ and G ^f ₁₂.

The sensitivity indices and rankings can also be explained by physical reality to some extent. Firstly, V _f nearly ranked as the main affecting factor for all elastic properties corresponds to the fact that the constitution of composite decides the main properties, which is directly reflected in V _f . Secondly, as carbon fibre presents much higher modulus and lower major Poisson's ratio than the epoxy matrix, carbon fibre is stronger to resist deformation; meanwhile, the transverse deformation of carbon fibre is extremely tiny under transverse load. For longitudinal loading conditions, the longitudinal modulus of carbon fibre leads the whole longitudinal deformation; meanwhile, Poisson's ratio is the combined action of Poisson's ratios of fibre and matrix. For transverse and in-plane shear loading conditions, the matrix is in the first stage to bear the deformation caused by transverse tensile and in-plane shear, so that the modulus and Poisson's ratio of the matrix seem to have more influence on the transverse and in-plane shear properties of the composites.

6. Conclusion

In this work, a global sensitivity analysis for the elastic constants of unidirectional CFRP based on finite element simulations and metamodels was performed to characterise the independent and cooperative effects of various microscopic parameters on the elastic properties and quantify the importance rankings of these parameters. A reliable three-phase unit cell finite element model in which the interface was modelled by cohesive elements governed by a bilinear traction-separation law was simulated to obtain the elastic properties of unidirectional CFRP by homogenisation method. By comparison with the experimental data, the modelling method was verified to be valid. Some other microstructure features of the studied composite such as non-periodic arrangement of carbon fibres should be taken into consideration for a finer finite-element model to gain accurate predicted results. Although this issue is beyond the scope of this paper, it is an important issue especially for predicting damage and failure, which deserves further research.

Kriging metamodels were constructed to establish a parameterised relationship between input microscopic parameters and elastic constants of composites by combining Optimal Latin Hypercube Design (OLHD) of the experiment method and finite element simulations. The Kriging metamodels present specifically high accuracy that are reliable for further directly predicting the elastic properties at any given input parameters, which greatly reduces the computational cost in the next RS-HDMR-based sensitivity analysis. RS-HDMR expansions were then efficiently constructed, based on rapid calculation of numerous sample points using Kriging metamodels. RS-HDMR expansions up to second-order were verified to be sufficient to represent the studied elastic constants of composites. Global sensitivity indices were computed following Sobol's method to quantify the influence of eleven input parameters on the elastic properties, the ranked indices clearly identify the most important parameters for each elastic property.

From the perspective of sensitivity indices of input parameters on the elastic outputs, the elastic properties of unidirectional CFRPs strongly rely on the elastic properties of carbon fibre and matrix, while slightly affected by the elastic properties of interface. As an important manufacturing parameter, carbon fibre volume fraction V_f has an obvious independent influence on five elastic properties, except for the major Poisson's ratio v₁₂. Moreover, V_f participates in the cooperative effects for all the individual elastic responses of composites. In the relationship between eleven microscopic input parameters and macroscopic elastic properties of studied composites, there exist some identical parameters that have obvious influence on different outputs, namely the coupled variables, which could be considered into further research on the influence of input uncertainties on the macroscopic properties.

The procedure proposed in this work is very efficient to calculate the global sensitivity indices without the need for large amounts of highly time-consuming simulations, which provides a straightforward method to study the property-structure relations, proving to be a useful tool to providing information for further design optimisation of CFRP.