# Application of the spectral stochastic finite element method for performance prediction of composite

Composite Structures 78 (2007) 447–456 www.elsevier.com/locate/compstruct

Application of the spectral stochastic ?nite element method for performance prediction of composite structures

M.F. Ngah *, A. Young

Applied Materials Group, Security and Dual Use Division, QinetiQ, Cody Technology Park, Farnborough GU14 0LX, United Kingdom Available online 22 December 2005

Abstract This paper demonstrates an application of the spectral stochastic ?nite element method (SSFEM) for predicting the performance of a composite structure with variable material constitutive properties. The SSFEM can provide an e?cient means of calculating statistical moments (e.g. mean and variance) and probability distributions associated with the probabilistic behaviour of a structure, and should therefore be suitable as an alternative to Monte Carlo Simulation for assessing structural reliability. A summary of the SSFEM theory is given, and it is shown how the method can be used to calculate the variability of strain and stress in a composite structure with variable material properties. An example is presented to demonstrate application of the method, and solutions are compared with results from other methods of analysis. ? 2005 Elsevier Ltd. All rights reserved.

` Keywords: Stochastic ?nite elements; Karhunen–Loeve expansion; Polynomial chaos; Composites

1. Introduction Variability in the performance of composite materials arises from the variability in constituent properties, ?bre distribution, structural geometry and loading conditions. Safety factors for loading and allowables for material strength have been used in structural design to account for the variability. Such a design procedure could result in a costly and unnecessary conservatism. Probabilistic design methodologies take into account the actual variability exhibited in materials property data, specimen geometry and applied loading. With the probabilistic approach, designs are quali?ed to have a given reliability value, in other words, to exhibit only a given acceptable level of failure probability. The USAF general speci?cation for aircraft structures, AFGS 87221A, speci?es the maximum acceptable frequency of structural failure leading to the loss of the aircraft as 1 · 10?7 occurrences per ?ight or 1 occurrence in 107 ?ights [1], which is further

*

Corresponding author. E-mail address: mfngah@qinetiq.com (M.F. Ngah).

reduced to 1 · 10?8 or 1 · 10?9 for the design of structural elements within the airframe. Current methods for modelling structural reliability are based on Monte Carlo Simulation (MCS) [2] schemes or approximate methods such as the First Order Reliability Methods (FORM) [3] or the Second Order Reliability Methods (SORM) [4]. MCS becomes increasingly expensive for high reliability calculations involving many random variables. The FORM/SORM methods on the other hand tend to su?er non-convergence and accuracy problems due to their approximate nature. Of immediate interest for structural reliability assessment is the spectral stochastic ?nite element method (SSFEM) [5–9] using random ?elds theory to model the spatial and stochastic variability within a composite structure. In the SSFEM, random ?elds representing the material constitutive properties and the response are expressed explicitly as series expansions in a set of normalised Gaussian random variables. The material properties are repre` sented as a Karhunen–Loeve Expansion (KLE) consistent with mean and autocovariance functions estimated from empirical measurement. A Polynomial Chaos Expansion

0263-8223/$ - see front matter ? 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruct.2005.11.009

448

M.F. Ngah, A. Young / Composite Structures 78 (2007) 447–456

(PCE) is used to represent the unknown stochastic response, and involves a set of polynomial basis functions of the random variables with coe?cients in the form of vectors of nodal parameters. The polynomial basis functions are de?ned to be mutually orthogonal with respect to expectation, and facilitate the subsequent numerical reduction of the problem. Using a Galerkin approach, the stochastic problem is converted into a large deterministic system of equations in the unknown nodal parameters from the PCE, in which the coe?cient matrix is blocksparse. Once the solution vectors of nodal parameters have been computed, the stochastic response in the form of the PCE can be used to generate statistical moments and probability distributions. The SSFEM is therefore suitable as an alternative to MCS for modelling stochastic structures and for assessing structural reliability. This paper describes how the SSFEM may be applied to composite structures, and is set out in four sections. Section 2 presents a summary of the SSFEM, and shows how it can be used to calculate the variability of strain and stress. In Section 3, an example is presented to demonstrate application and results of the method, and conclusions from the study are given in Section 4. 2. The spectral stochastic ?nite element method In the following, spatial position is denoted by the vector x, and the solid volume occupied by the composite structure is denoted by V. Random variability is denoted by a parameter x, and expectation of a random ?eld f(x, x) is denoted by hf(x, x)i. A ?nite element model of a structure involves a representation of the solid geometry (element mesh), applied loading (boundary conditions) and material properties. When the material properties have random variability, the matrix of constitutive coe?cients D, relating stress and strain r = D ? e, can be de?ned as a continuous random ?eld D(x, x) consistent with a given mean D?x? ? hD?x; x?i and autocovariance function CDD ?x; x0 ? ? h?D?x; x? ? D?x?? ?D?x0 ; x? ? D?x0 ??i ?1?

where RDD(x, x 0 ) is a scalar autocovariance function. The random ?eld D(x, x) associated with a given autocovari` ance function (2) may be represented as a Karhunen–Loeve Expansion (KLE) 1 X p????? ^ D?x; x? ? D?x? ? D nm ?x? km Km ?x? ?3?

m?1

de?ned in terms of a set of eigenvalues km and eigenfunctions Km(x) obtained from the spectral decomposition of the scalar autocovariance function 1 X RDD ?x; x0 ? ? km Km ?x?Km ?x0 ? ?4?

m?1

by solution of the integral eigenproblem [5] Z RDD ?x; x0 ?Km ?x? dV ?x? ? km Km ?x0 ?; V Z Km ?x?Kn ?x? dV ?x? ? dmn

V

?5?

In the KLE (3), the stochastic variability of D(x, x) is ex^ pressed in terms of a matrix D, that incorporates the degree of variability (i.e. the standard deviation) of each component of D(x, x), and a set of uncorrelated stochastic functions nm(x), m = 1, 2, . . . , which are regarded as independent normalised Gaussian random variables with hnmi = 0 and hnmnni = dmn. The sequence of km and Km(x), m = 1, 2, . . . , is ordered in decreasing magnitude of the eigenvalues, so that the in?nite series (3) can be truncated after NKL terms to provide a manageable approximation to the random ?eld D(x, x). When a truncated KLE for the material properties D(x, x) is used in a ?nite element scheme, there results a system of equations of the form

N KL X m?0

nm ?K m ?fug ? ff g ?

N KL X m?0

~ u nm ?K m ?f~g ?

N KL X m?0

nm ffm g

?6?

where denotes outer (tensor) product. The autocovariance function quanti?es the stochastic variability of a random ?eld about its mean and how that variability correlates at di?erent spatial locations. In general, each component of the constitutive material matrix would have its own mean and autocovariance functions, which could be estimated from statistical measurements on samples of the material. However, in order to simplify the description and application of the method, it is assumed that the constitutive coe?cient matrix is of the form D?x; x? ? ^ D?x? ? f ?x; x?D, where f(x, x) is a zero-mean scalar ran^ dom ?eld and D is a constant matrix. Then the autocovariance function (1) is of the form ^ ^ CDD ?x; x0 ? ? D DRDD ?x; x0 ? ?2?

where an extra parameter n0 = 1 is de?ned for notational convenience. The unknowns are contained in the nodal response vector {u} and the boundary conditions in the nodal force vectors {fm}, m = 0, 1, 2, . . . , NKL, which are composed of an applied traction/force vector {f} and speci?ed nodal response values f~g. The matrix terms are scheu matically of the form Z T ?K 0 ? ? fg0 ?x?g ? D?x? ? fg0 ?x?g dV ?x? ZV p????? ?7? T ^ ?K m ? ? fg0 ?x?g ? D ? fg0 ?x?g km Km ?x? dV ?x?;

V

m ? 1; 2; . . . ; N KL where {g 0 (x)} involves shape function derivatives associated with strain at x expressed in terms of the response at the nodal locations. The matrix [K0] is the ?nite element global sti?ness matrix obtained from the mean material properties D?x?, while the subsequent matrices [Km], m = 1, 2, . . . , NKL, are calculated analogously using the p????? ^ KLE terms D km Km ?x? in place of D?x?.

M.F. Ngah, A. Young / Composite Structures 78 (2007) 447–456

449

Using computer-generated pseudo-random numbers for the Gaussian variables n1 ; n2 ; . . . ; nN KL , the ?nite element equation system (6) can be used as a basis for MCS, from which realizations of the stochastic response u(x, x) may be computed and used to accumulate statistical moments. The SSFEM [5] provides an alternative solution methodology, in which the random ?eld u(x, x) is approximated as a Polynomial Chaos Expansion (PCE) up to polynomial order p of the form u?x; x? ?

Np X j?0

where the coe?cients cmjk ? hnm Wj ?n?Wk ?n?i ?14?

Wj ?n?vj ?x?

?8?

where Wj(n), j = 0, 1, 2, . . . , Np, are chaos polynomials in the NKL random variables n ? fn1 ; n2 ; . . . ; nN KL g, vj(x) are a set of spatial functions to be determined, and the number of terms in the PCE of order p in NKL random variables is ?p ? N KL ?! Np ? 1 ? ?p?!?N KL ?! ?9?

are integers and may be evaluated as described in Appendix A. Owing to the orthogonality properties of the chaos polynomials (10), the coe?cients cmjk are mostly zero. The set of equations (13) may be regarded as a single super-system in which the solution super-vector is composed of a concatenation of the nodal vectors {vj}, j = 0, 1, 2, . . . , Np, and the corresponding coe?cient super-matrix is block-sparse. The super-matrix may not be positive-de?nite, and so (13) is solved with the preconditioned MINRES algorithm [10] using the Choleski decomposition of [K0]. Once the solution vectors {vj} have been determined, the mean and autocovariance of the response may be obtained as u?x? ?

Np X j?0

hWj ?n?ivj ?x? ? v0 ?x?

Np X 2 h?Wj ?n?? ivj ?x? vj ?x0 ? j?1

?15? ?16?

Eq. (8) is essentially equivalent to a power series expansion of order p in the NKL random variables. However, expressing the response in terms of chaos polynomials is computationally advantageous in that it results in a sparse deterministic system of equations derived from the stochastic Eq. (6). Chaos polynomials are de?ned to be orthogonal with respect to expectation, i.e. hWj ?n?Wk ?n?i ? 0 if j 6? k ?10?

Cuu ?x; x0 ? ?

where vj(x) = {g(x)}T{vj} is the contribution to the response at location x from the jth nodal solution vector and h(Wj(n))2i can be evaluated as shown in Appendix A. The strain and stress at a location x within the structure are given by e?x; x? ?

Np X j?0

and for convenience are ordered such that the ?rst NKL + 1 polynomials are the order 0 term W0(n) = 1 and the order 1 terms Wj(n) = nj (j = 1, 2, . . . , NKL). For Gaussian random variables n, the sequence of chaos polynomials may be constructed as products of one-dimensional monic Hermite polynomials, as described in Appendix A. In a ?nite element scheme, the nodal response vector corresponding to the expansion (8) is fu?x?g ?

Np X j?0

Wj ?n?ej ?x?

?17?

r?x; x? ? D?x; x? ? e?x; x? ?

N KL N p XX m?0 j?0

nm Wj ?n?Dm ?x? ? ej ?x?

?18?

Wj ?n?fvj g

?11?

where each vj(x) = {g(x)}T{vj} is represented in terms of nodal values {vj} and ?nite element shape functions {g(x)}. Then, using the expansion (11) in the equation system (6) gives

N p N KL XX j?0 m?0

where ej(x) = {g 0 (x)}T{vj} are strain-type values at x derived from the jth nodal solution vector and the e?ective material coe?cients are given by D0 ?x? ? D?x? and p????? ^ Dm ?x? ? D km Km ?x? for m = 1, 2, . . . , NKL. De?ning rmj(x) = Dm(x) ? ej(x), the mean and autocovariance of strain and stress can be obtained using

Np X e?x? ? hWj ?n?iej ?x? ? e0 ?x? j?0 Np X j?1

?19? ?20? ?21?

nm Wj ?n??K m ?fvj g ?

N KL X m?0

nm ffm g

?12?

Cee ?x; x0 ? ? r?x? ?

h?Wj ?n?? iej ?x? ej ?x0 ?

N KL X m?0

2

The orthogonality of the chaos polynomials may be exploited in a Galerkin procedure to reduce (12) to a block-sparse system of deterministic equations in the (Np + 1) sets of nodal values {vj}, j = 0, 1, 2, . . . , Np. Multiplying (12) through by Wk(n) and taking expectation gives

N p N KL XX j?0 m?0

N KL N p XX m?0 j?0 0

hnm Wj ?n?irmj ?x? ?

rmm ?x?

Np N p N KL N KL XXXX Crr ?x; x ? ? ?X mnjk rmj ?x? rnk ?x0 ?? k?0 j?0 n?0 m?0

cmjk ?K m ?fvj g ?

N KL X m?0

? r?x? r?x0 ? dkm ffm g; k ? 0; 1; 2; . . . ; N p ?13?

?22?

The non-zero coe?cients Xmnjk = hnmnnWj(n)Wk(n)i can be identi?ed and calculated analogously to the coe?cients cmjk de?ned in (14) using a recurrence formula for

450

M.F. Ngah, A. Young / Composite Structures 78 (2007) 447–456

250mm

25mm

1

2

3

4

5

6

7

8

9

10

Fig. 1. 250 mm · 25 mm panel under uniaxial loading modelled using ten nine-noded ?nite elements.

0.035

0.035

0.03

FOSM

0.03

FOSM

0.025 SSFEM: p= 1

0.025 SSFEM: p= 1

cv(eps11)

cv(eps11)

0.02

0.02

0.015

SSFEM: p= 4

0.015

SSFEM: p= 4

0.01 MCS (10000)

0.01 MCS (10000)

0.005

0.005

0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x/L

x/L

cv(D11) = 0.04, NKL = 4

0.2 0.18 FOSM 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 MCS (10000) SSFEM: p= 4 SSFEM: p= 1 0.16 0.14 0.12 0.1 0.08 0.06 0.04 0.02 0 0 0.1 0.2 0.2 0.18

cv(D11) = 0.04, NKL = 8

FOSM

cv(eps11)

cv(eps11)

SSFEM: p= 1

SSFEM: p= 4

MCS (10000)

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x/L

x/L

cv(D11) = 0.20, NKL = 4

0.35 0.35

cv(D11) = 0.20, NKL = 8

0.3

FOSM

0.3

FOSM

0.25 SSFEM: p= 1

0.25 SSFEM: p= 1

cv(eps11)

cv(eps11)

0.2

0.2

0.15

SSFEM: p= 4

0.15

SSFEM: p= 4

0.1 MCS (10000)

0.1 MCS (10000)

0.05

0.05

0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x/L

x/L

cv(D11) = 0.32, NKL = 4

cv(D11) = 0.32, NKL = 8

Fig. 2. Strain variability cv(e11) plotted against distance x along the length of the panel (L = 250 mm).

M.F. Ngah, A. Young / Composite Structures 78 (2007) 447–456

451

one-dimensional Hermite polynomials, as described in Appendix A. The cumulative distribution function (CDF) and the probability density function (PDF) of the response or strain or stress may be estimated very e?ciently from Monte Carlo realizations of the polynomial expansions (11) or (17) or (18) respectively. 3. Example problem Application of the SSFEM is illustrated with an example in which a rectangular panel of T800/924 unidirectional carbon ?bre composite with a given variability in elastic properties is loaded under uniaxial tension. SSFEM solutions are compared with results from Monte Carlo Simulation (MCS) and the First Order Second Moment (FOSM) Perturbation Method summarised in Appendix B. The unidirectional rectangular panel, illustrated in Fig. 1, has dimensions 0 mm 6 x 6 250 mm by 0 mm 6 y 6 25 mm (nominal thickness 2 mm) and is considered as an orthotropic two-dimensional sheet under generalised plane stress. The nominal material properties in the longitudinal ?bre 1-direction (i.e. x) and the transverse 2direction (i.e. y) are de?ned by the stress–strain law 0 1 2 30 1 r11 D11 D12 e11 0 B C 6 7B C ?23? r22 A ? 4 D12 D22 0 5@ e22 A @ 0 0 D33 r12 c12 and are assumed to have mean values D11 ? 159:85 GPa, D22 ? 9:45 GPa, D12 ? 2:84 GPa, D33 ? 4:00 GPa and ^ ^ standard deviations D11 ? 0:04D11 , D22 ? 0:02D22 , ^ 12 ? 0:07D12 , D33 ? 0:11D33 . These parameters de?ne the ^ D ^ material matrices D and D in the KLE (3) for the random ?eld D(x, x), while the scalar autocovariance function in (4) and (5) is assumed to be in the form of an inverse exponential RDD ?x; x0 ? ? e?jx?x j=b ? e?jy?y j=b

0 0

?24?

with correlation length b = 125 mm. In general, the eigenvalues and eigenfunctions required in the KLE (3) would ` be obtained by numerical solution of the Karhunen–Loeve eigenproblem (5). However, for the idealised autocovariance function (24), the solution of (5) is available in closed-form [5]. The panel is modelled using ten quadrilateral ninenoded quadratic ?nite elements, and the applied loading is deterministic, with one end, x = 0, ?xed and the other end, x = 250 mm, subjected to a unit extension u1 = 1 mm. The variable chosen to characterise the solutions is the strain component e11(x, x) calculated at the centre of each ?nite element. Its mean value is found to be almost uniform along the length of the panel, and so the variability of e11 is illustrated in terms of the coe?cient of variation ??????????????????????? p cv(e11) de?ned as the ratio of the standard deviation C e11 e11 ?x; x? to the mean e11 ?x?. The stochastic models MCS, FOSM and SSFEM are all based on a system of ?nite element Eq. (6) in which the stochastic variability of the material sti?ness and of the response is expressed in terms of NKL random variables. In each case, the size of the numerical problem to be solved is directly related to the number of random variables used, and so NKL should not be larger than is necessary to provide a realistic representation of the required stochastic problem. In the SSFEM, the polynomial order p of the PCE is also a parameter that ideally should be kept as low as possible for computational e?ciency. In MCS, the computational e?ort is proportional to the number of simulations Nmcs. Therefore, results are presented below to demonstrate how the numerical solutions depend on the parameters NKL, p and Nmcs. The convergence of the solutions is found to depend on the stochastic variability of the material, and this is also considered by varying the compo^ nent D11 such that the coe?cient of variation ^ cv?D11 ? ? D11 =D11 is increased from its nominal value of 0.04 up to 0.32.

0.35

FOSM

0.3

SSFEM: p = 1

0.25

cv (eps11)

SSFEM: p = 2

0.2

SSFEM: p = 3

0.15

SSFEM: p = 4

0.1

MCS (10000)

0.05

MCS (50000)

0 0 0.04 0.08 0.12 0.16 0.2 0.24 0.28 0.32

cv (D11)

Fig. 3. Strain variability cv(e11) at the centre of element 10.

452

M.F. Ngah, A. Young / Composite Structures 78 (2007) 447–456

3.1. Convergence of solutions Solutions obtained from the SSFEM with PCE orders p = 1 and p = 4, from MCS with Nmcs = 10 000 simulations and from the FOSM Perturbation Method are compared in Fig. 2. The normalised standard deviation cv(e11) calculated at the centres of the ten elements is plotted against distance x along the panel, and results are given for both NKL = 4 and NKL = 8 random variables, with material var4500 4000 3500 3000

iability cv(D11) = 0.04, 0.20 and 0.32. The SSFEM with p = 1 and the FOSM give almost identical results that agree with the MCS only for low material variability (i.e. cv(D11) < 0.1). In contrast, it is found that the SSFEM with p = 4 gives reasonable agreement with MCS over the entire range of cv(D11) considered. Using more random variables in the KLE (3) provides a better approximation to the random ?eld D(x, x) consistent with the exponential autocovariance (24), and the solutions using NKL = 8 are

1.2

MCS

1

SSFEM: p=1 SSFEM: p=4

0.8

PDF

CDF

2500 2000 1500 1000

0.6

MCS

0.4

SSFEM: p=1 SSFEM: p=4

0.2

500 0 0.0035

0 0.0035

0.004

0.0045

0.004

0.0045

eps11

eps11

PDF with cv(D11) = 0.04

900 800 700 600 1.2

CDF with cv(D11) = 0.04

MCS

1

SSFEM: p=2 SSFEM: p=4

0.8

CDF

PDF

500 400 300 200

0.6

MCS

0.4

SSFEM: p=2

0.2

100 0 0.001 0 0.001

SSFEM: p=4

0.004 0.007

0.004

0.007

eps11

eps11

PDF with cv(D11) = 0.20

600 1.2

CDF with cv(D11) = 0.20

MCS

500 1

SSFEM: p=2

400

SSFEM: p=4

0.8

300

CDF

PDF

0.6

MCS

200 0.4

SSFEM: p=2 SSFEM: p=4

100

0.2

0 0 0.004 0.008

0 0 0.004 0.008

eps11

eps11

PDF with cv(D11) = 0.32

CDF with cv(D11) = 0.32

Fig. 4. PDF and CDF for the strain e11 at the centre of element 10 obtained from MCS (100 000) and from SSFEM solutions.

M.F. Ngah, A. Young / Composite Structures 78 (2007) 447–456

453

noticeably smoother that those using NKL = 4. However, there is not a great di?erence in the results obtained using NKL = 4 and 8, and it is found that good estimates of the standard deviation and mean of the strain or stress can be obtained using a low order KLE.

Convergence of solutions from the SSFEM with increasing polynomial order p (= 1, 2, 3 and 4) and from MCS with number of simulations Nmcs (= 10 000 and 50 000) is shown in Fig. 3, where values of cv(e11) calculated at the centre of element 10 (see Fig. 1) are plotted against material

1 0.99999 0.99998 0.99997 0.99996

MCS

SSFEM: p=2 SSFEM: p=3 SSFEM: p=4 SSFEM: p=6 SSFEM: p=8

0.00445 0.0045 0.00455

CDF

0.99995 0.99994 0.99993 0.99992 0.99991 0.9999 0.0044

(a)

1 0.99998 0.99996 0.99994 0.99992

eps11

MCS

SSFEM: p=2 SSFEM: p=3 SSFEM: p=4 SSFEM: p=6 SSFEM: p=8

CDF

0.9999 0.99988 0.99986 0.99984 0.99982 0.9998 0.006

0.007

0.008

0.009

0.01

(b)

1

eps11

MCS

0.995

SSFEM: p=2

CDF

0.99

SSFEM: p=3 SSFEM: p=4 SSFEM: p=6

0.985

0.98 0.004

0.009

0.014

0.019

(c)

eps11

Fig. 5. Upper tail of the CDF for the strain e11 at the centre of element 10 with material variability: (a) cv(D11) = 0.04, (b) cv(D11) = 0.20 and (c) cv(D11) = 0.32.

454

M.F. Ngah, A. Young / Composite Structures 78 (2007) 447–456

variability cv(D11). It can be seen that the MCS results differ noticeably from the FOSM results for cv(D11) > 0.16, but are estimated fairly well by the SSFEM with p P 3. However, there are discrepancies apparent for the highest variability cv(D11) = 0.32 that do not diminish as the polynomial order is increased, and it is found that the SSFEM solution diverges as p is increased further; e.g. using p = 6 gives cv(e11) = 0.58. Furthermore, there are indications that the MCS solutions for cv(D11) = 0.32 are not converging with increasing number of simulations Nmcs. 3.2. Estimation of response probability distributions For low material variability, the FOSM Perturbation Method is capable of providing accurate estimates of the statistical moments of the response with minimal computational e?ort. In such cases, the SSFEM with polynomial order p = 1 gives similar results. However, an advantage of using SSFEM is that the probability distributions associated with the response can also be obtained. Once the solution (11) has been computed, Monte Carlo realizations of the response (8), strain (17) or stress (18) may be used to estimate the cumulative distribution function (CDF) and the probability density function (PDF). Examples are shown in Fig. 4, where the CDF and PDF for the strain e11 at the centre of element 10 obtained from Nmcs = 100 000 MCS solutions of the stochastic equation system (6) are compared with those from SSFEM solutions using 100 000 realizations of the strain (17). The results indicate that SSFEM with polynomial order as low as p = 1 or 2 can give reasonable overall approximations to the probability distributions.

However, in order to capture the characteristics of the tails of the distributions (as required for assessing structural reliability), a higher polynomial order may be required, particularly for high material variability. This is illustrated in Fig. 5a–c, which show the upper tail of the CDF for the cases cv(D11) = 0.04, 0.20 and 0.32, respectively, obtained from MCS with Nmcs = 1 000 000 and from SSFEM with 1 000 000 realizations of the strain solution (17). For cv(D11) = 0.04, a good estimate of the CDF can be obtained using SSFEM with polynomial order as low as p = 2. For cv(D11) = 0.20, the SSFEM distributions converge towards the MCS distribution as the polynomial order p is increased, and the results using p P 6 are very accurate for CDF(e11) up to 0.9999. However, the SSFEM distributions for cv(D11) = 0.32 approximate the MCS distribution only for p 6 4, and solutions appear to diverge as the polynomial order is increased further (no results were obtained for p = 8, since iterative solution of the SSFEM equation system (13) did not converge). 3.3. Comparison of computer run times Tables 1 and 2 illustrate the relative accuracy of results and the computer run times (in seconds on a 3 GHz PC) obtained for the example test case using each method with NKL = 4 random variables. Accuracy of cv(e11) is assessed by comparison with results using MCS with Nmcs = 100 000 simulations. The FOSM Perturbation Method is the fastest, and can give good estimates of the mean and standard deviation for low material variability cv(D11) < 0.1. The SSFEM with polynomial order p = 2 is accurate for a

Table 1 Relative di?erence in the values cv(e11) calculated at the centre of element 10 for each method compared with MCS (100 000) using NKL = 4 cv(D11) 0.04 0.08 0.12 0.16 0.20 0.24 0.28 0.32 FOSM (%) SSFEM (%) p=1 ?0.3 ?0.8 ?2.1 ?3.9 ?6.4 ?9.9 ?17.6 ?26.4 ?0.3 ?0.8 ?2 ?3.8 ?6.2 ?9.6 ?17.3 ?26 p=2 ?0.1 0.1 ?0.1 ?0.4 ?1.1 ?2.3 ?8 ?14.9 p=3 ?0.1 0.1 0 ?0.2 ?0.4 ?0.9 ?5.4 ?10.4 p=4 ?0.1 0.1 0 ?0.2 ?0.3 ?0.5 ?4.2 ?6.9 p=5 ?0.1 0.1 0 ?0.2 ?0.3 ?0.4 ?3.3 9 p=6 ?0.1 0.1 0 ?0.2 ?0.3 ?0.3 ?1 92.9 MCS (%) (10 000) ?0.4 ?0.6 0.4 0.7 ?0.1 1.1 ?4.7 ?3.8 (50 000) 0.2 ?0.3 ?0.1 ?0.2 0.2 0 ?3.2 5.7

Table 2 Computer run times (seconds) for each method using NKL = 4 cv(D11) 0.04 0.08 0.12 0.16 0.20 0.24 0.28 0.32 FOSM SSFEM p=1 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 1.2 1.2 1.3 1.2 1.2 1.2 1.2 1.2 p=2 1.5 1.4 1.5 1.5 1.5 1.5 1.5 1.6 p=3 2.2 2.2 2.3 2.5 2.6 2.7 2.9 3.2 p=4 4.7 4.9 5.3 6.0 6.6 7.5 8.5 10.4 p=5 13.4 14.6 16.0 19.1 21.8 25.2 32.5 75.2 p=6 44.1 47.5 56.6 65.0 77.7 97.9 178.1 1062.2 MCS (10 000) 22.1 22.8 24.7 26.5 28.1 29.6 31.3 33.2 (50 000) 108.1 111.3 121.4 129.7 137.9 145.7 153.8 162.5

M.F. Ngah, A. Young / Composite Structures 78 (2007) 447–456

455

wider range of material variability cv(D11) < 0.2, while polynomial order p = 3 gives good results for material variability cv(D11) < 0.28. For these cases, the SSFEM is faster than MCS with 10 000 simulations. However, none of the methods perform well for high material variability cv(D11) P 0.28; there are signi?cant di?erences in the MCS results using 50 000 and 100 000 simulations, and the SSFEM solutions for cv(D11) = 0.32 are unstable and become highly inaccurate with polynomial order p = 6. The divergence of the stochastic response (using any method of analysis) when Gaussian random variables are used to represent the material properties is a known problem [8], and arises since the assumed random ?eld D(x, x) allows the possibility of zero (or nearly zero) structural sti?ness. In MCS, realizations that are obviously non-physical may be identi?ed (e.g. if D11 6 0 at any of the ?nite element quadrature points) and omitted from the simulation, but they are implicitly present within the SSFEM solutions, which are found to become unstable when the material variability is high. 4. Conclusions The SSFEM has been applied to a test problem where a rectangular unidirectional composite panel with variable material properties is loaded in tension, and solutions have been compared with results from MCS and from the FOSM Perturbation Method. Although the FOSM is very e?cient for low material variability (standard deviation up to about 10% of mean), the SSFEM is applicable over a wider range of material variability (standard deviation up to about 24% of mean). Furthermore, the SSFEM is capable of estimating probability distributions for the stochastic behaviour of the composite more e?ciently than MCS. This makes the SSFEM a potentially useful technique for application to structural reliability assessment. In the basic formulation of the SSFEM used above, the KLE for the material properties is a linear superposition of Gaussian random variables, so that the structural sti?ness is assumed to be a Gaussian random ?eld. Whereas this assumption is theoretically convenient, it is unrealistic in that it allows the possibility of zero or negative structural sti?ness, and it is found that stochastic solutions fail to converge when material variability is high [8]. This problem can be avoided by representing the material properties by a suitable non-Gaussian random ?eld that is consistent with observed material variability [6,9]. A reformulation of the SSFEM involving non-Gaussian material properties will provide a more realistic probabilistic modelling method for assessing reliability of composite structures. Acknowledgements The authors would like to acknowledge funding from UK MoD under EUT contract EUT/04/31/D03.

Appendix A. Chaos polynomials The set of chaos polynomials Wk(n), k = 0, 1, 2, . . . , Np, up to order p in NKL Gaussian random variables n ? fn1 ; n2 ; . . . ; nN KL g are constructed as products of onedimensional monic Hermite polynomials hr(n), r = 0, 1, 2, . . . , that are orthogonal with respect to expecta2 tion (with Gaussian weight function e?n =2 ) Z ?1 2 hr ?n?hs ?n?e?n =2 dn ? 0 if r 6? s hhr ?n?hs ?n?i ?

?1

?A:1? The Hermite polynomials are generated using the recurrence formula h0 ?n? ? 1; h1 ?n? ? n; hr?1 ?n? ? nhr ?n? ? rhr?1 ?n? ?A:2? from which the following expectations may be evaluated hhr ?n?i ? dr;0 ; hhr ?n?hs ?n?i ? r!dr;s ?A:3? ?A:4?

hnhr ?n?hs ?n?i ? ?r ? 1?!dr;s?1 ? r!dr;s?1 hn2 hr ?n?hs ?n?i ? ?r ? 2?!dr;s?2 ? r!?2r ? 1?dr;s ? r!dr;s?2

?A:5? Each NKL-dimensional chaos polynomial term Wk(n) is of the form Wk ?n? ?

N KL Y m?1

hik ?nm ? m

?A:6?

where the set of non-negative integer indices fik ; ik ; . . . ; ik KL g is unique to each index k, and the entire 1 2 N set k = 0, 1, 2, . . . , Np covers all combinations such that PN KL k 0 6 m?1 ?im ? 6 p. The non-zero NKL-dimensional coe?cients cmjk that arise in the SSFEM matrix system (13) Np can be identi?ed from the index lists fik ; ik ; . . . ; ik KL gk?0 1 2 N and evaluated using (A.3) and (A.4). For each index k, there are at most (2NKL + 1) non-zero values cmjk ? hnm hi?j? ?nm ?hi?k? ?nm ?i ?

m m

N KL Y r?1 ?r6?m?

hhi?j? ?nr ?hi?k? ?nr ?i

r r

?A:7?

with indices j that can be identi?ed as follows. For m = 0 (with n0 = 1), the only non-zero cmjk is the one for which each j-index is equal to the corresponding k-index, i.e. ij ? ik for r = 1, 2, . . . , NKL. Then j = k, and (A.3) gives r r c0kk ? h?Wk ?n?? i ?

2 N KL Y ?ik ?! r r?1

?A:8?

For each m 5 0, there are at most two non-zero cmjk associated with each index k: one where the index j is determined by ij ? ik ? 1 (provided that ik 6 p ? 1) and m m m ij ? ik ?r 6? m?, for which r r

456

N KL Y ?ij ?! r r?1

M.F. Ngah, A. Young / Composite Structures 78 (2007) 447–456

cmjk ?

?A:9?

into (6), ignoring all terms higher than second order in the random variables and isolating the coe?cients of the linearly independent random terms gives ?K 0 ?fz0 g ? ff0 g ?K 0 ?fzm g ? ffm g ? ?K m ?fz0 g; m ? 1; 2; . . . ; N KL N KL N KL X X ?K 0 ? fzmm g ? ?2 ?K m ?fzm g

m?1 m?1

and one where the index j is determined by ij ? ik ? 1 m m (provided that ik P 1) and ij ? ik (r 5 m), for which m r r cmjk ?

N KL Y ?ik ?! r r?1

?B:5? ?B:6? ?B:7?

?A:10?

The coe?cients Xmnjk required for the autocovariance of stress in (22), i.e. X mmjk ? h?nm ? hijm ?nm ?hik ?nm ?i ? m and X mnjk ? hnm hijm ?nm ?hik ?nm ?i ? hnn hijn ?nn ?hik ?nn ?i m n ?

N KL Y r?1 ?r6?m;n? 2 N KL Y

r?1 ?r6?m?

hhijr ?nr ?hik ?nr ?i r

?A:11?

hhijr ?nr ?hik ?nr ?i r

?A:12?

for m 5 n, can be evaluated analogously to the procedure above for cmjk, using the formulae (A.3)–(A.5). Appendix B. FOSM perturbation method The First Order Second Moment Perturbation Method provides a rapid means of estimating statistical moments of the solution of the stochastic system (6). If the response random ?eld is approximated as a second order polynomial in a set of independent Gaussian random variables n1 ; n2 ; . . . ; nN KL u?x; x? ? z0 ?x? ?

N KL X m?1

The set of (NKL + 2) nodal vectors {z0}, {zm} PN KL (m = 1, 2, . . . , NKL) and m?1 fzmm g required to evaluate the mean (B.2) and autocovariance (B.3) of the response can be recovered by solving (B.5), (B.6) and (B.7) in turn using the Choleski decomposition of the mean sti?ness matrix PN KL [K0]. The response functions z0(x), zm(x), and m?1 zmm ?x? and the corresponding strains e0(x), em(x) PN KL and m?1 emm ?x? may be determined for any required point x from the nodal solution vectors and the ?nite element shape functions in the usual way. The mean and autocovariance of the strain and stress are then evaluated using N KL 1X e?x? ? e0 ?x? ? emm ?x? ?B:8? 2 m?1 Cee ?x; x0 ? ?

N KL X m?1

em ?x? em ?x0 ?

N KL X m?1

?B:9? ?B:10?

r?x? ? D0 ?x? ? e?x? ? Crr ?x; x0 ? ?

Dm ?x?.em ?x?

N KL X ?D0 ?x?.em ?x? ? Dm ?x? ? e0 ?x?? m?1

?D0 ?x? ? em ?x0 ? ? Dm ?x? ? e0 ?x0 ?? References

?B:11?

nm zm ?x? ?

N KL N KL 1XX nm nn zmn ?x? 2 n?1 m?1

?B:1? with hnmi = 0 and hnmnni = dmn, then the mean and autocovariance of the response are given by u?x? ? z0 ?x? ? ? z0 ?x? ? Cuu ?x; x0 ? ? ?

N KL N KL 1XX hn n izmn ?x? 2 n?1 m?1 m n N KL 1X zmm ?x? 2 m?1

?B:2?

N KL N KL XX n?1 m?1 N KL X m?1

hnm nn izm ?x? zn ?x0 ? ?B:3?

zm ?x? zm ?x0 ?

Substituting the nodal values corresponding to the expansion (B.1) fug ? fz0 g ?

N KL X m?1

nm fzm g ?

N KL N KL 1XX n n fzmn g 2 n?1 m?1 m n

?B:4?

[1] Long MW, Narciso JD. Probabilistic design methodology for composite aircraft structures. DOT/FAA/AR-99/2, June 1999. [2] Mooney CZ. Monte Carlo simulation. Sage Publications Inc.; 1997. [3] Hasofer AM, Lind NC. Exact and invariant second moment code format. J Eng Mech Division, ASCE 1974;100:111–21. [4] Fiessler B, Neumann HJ, Rackwitz R. Quadratic limit states in structural reliability. J Eng Mech, ASCE 1979;105:661–76. [5] Ghanem RG, Spanos PD. Stochastic ?nite elements—a spectral approach. revised ed. New York: Dover Publications Inc.; 2002. [6] Ghanem R. Ingredients for a general purpose stochastic ?nite elements implementation. Comput Methods Appl Mech Eng 1999;168:19–34. [7] Keese A. A review of recent developments in the numerical solution of stochastic partial di?erential equations (stochastic ?nite elements). Institute of Scienti?c Computing, Technical University of Braunschweig, Brunswick, Germany, Informatikbericht Nr: 2003-06, 2003. Available from: http://opus.tu-bs.de/opus/volltexte/2003/504/. [8] Schevenels M, Lombaert G, Degrande G. Application of the stochastic ?nite element method for Gaussian and non-Gaussian systems. In: Proc ISMA2004, 2004, p. 3299–314. Available from: http://www.kuleuven.ac.be/bwm/papers/scheip04a.pdf. [9] Matthies HG, Keese A. Galerkin methods for linear and nonlinear elliptic stochastic partial di?erential equations. Comput Methods Appl Mech Eng 2005;194:1295–331. [10] Golub GH, Van Loan CF. Matrix computations. Baltimore: Johns Hopkins University Press; 1996.