- Email: [email protected]
Stochastic finite element analysis of shells with combined random material and geometric properties
Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160 www.elsevier.com/locate/cma
Stochastic finite element analysis of shells with combined random material and geometric properties G. Stefanou, M. Papadrakakis
*
Institute of Structural Analysis and Seismic Research, National Technical University, Athens 15780, Greece Received 12 June 2003; received in revised form 24 September 2003; accepted 3 October 2003
Abstract In this work, an extended stochastic formulation of the triangular composite facet shell element TRIC is presented for the case of combined uncertain material (YoungÕs modulus, PoissonÕs ratio) and geometric (thickness) properties. These properties are assumed to be described by uncorrelated two-dimensional homogeneous stochastic fields. The stochastic finite element analysis of shell structures is performed using the spectral representation method for the description of the random fields in conjunction with Monte Carlo simulation (MCS) for the computation of the response variability. Useful conclusions regarding the influence of each one of the structural parameters on the response variability are derived from the numerical tests examined. 2003 Elsevier B.V. All rights reserved. Keywords: Shell finite element; Stochastic analysis; Response variability; Spectral representation
1. Introduction The analysis of shells presents a challenge, since their formulation may become cumbersome and their behavior can be unpredictable with regard to geometry or support conditions. For this reason, shells have been considered as the ‘‘prima donna of structures’’ in the sense that their performance depends very much on how they are designed and how they are treated. The extreme sensitivity of thin shells to imperfections in material, geometry and boundary conditions can be described reasonably well, if the assumption of a random fluctuation of these imperfections is induced in the analysis, in other words if a stochastic finite element analysis is performed. The need for a robust, accurate and computationally efficient shell element becomes even greater for the computationally expensive stochastic finite element analysis of large realistic structures. The TRIC shell element [4] is a simple but sophisticated triangular, shear-deformable facet shell element suitable for the analysis of thin and moderately thick isotropic as well as composite plate and shell structures. Its formulation is based on the natural mode finite element method, which has a number of
*
Corresponding author. Tel.: +301-772-1694; fax: +301-772-1693. E-mail addresses: [email protected] (G. Stefanou), [email protected] (M. Papadrakakis).
0045-7825/$ - see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.10.001
140
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
computational advantages compared to the conventional isoparametric finite element formulations. The treatment of the element kinematics (inclusion of the transverse shear deformations in its formulation based on a first-order shear-deformable beam theory) eliminates locking phenomena in a physical manner. Several investigations have confirmed the rigorous theoretical basis of the element, while numerical examples have verified its locking free properties as well as its computational efficiency in linear and nonlinear structural applications [1,2,4]. The most important advantage of the TRIC element from the stochastic analystÕs point of view is the fact that there is no need to perform numerical integration for the computation of its deterministic stiffness matrix, which is carried out in closed form. As a result of this special feature, the TRIC element provides an ideal basis for the formulation of a computationally efficient stochastic stiffness matrix and for the use of the element in large-scale finite element stochastic computations. The majority of work done in the context of stochastic finite elements has been based on the assumption that only one material or geometric property is described by a stochastic field, e.g. [3,5,9,10,17,21]. For example, the assumption is often made that YoungÕs modulus is randomly varying over space, while PoissonÕs ratio is a deterministic constant. However, many of the physical mechanisms that lead to random variations of YoungÕs modulus also lead to random variations in other material properties such as PoissonÕs ratio. In this work, the stochastic stiffness matrix of the TRIC shell element is computed taking into account random variations of YoungÕs modulus and PoissonÕs ratio as well as the thickness of the shell. These random parameters are assumed to be uncorrelated. They are described by independent two-dimensional univariate (2D-1V) homogeneous Gaussian stochastic fields, which are represented via the spectral representation method. The stochastic stiffness matrix of the shell element is derived in the form of local average [15,20] and weighted integral [11,15] methods leading to a cost-effective matrix, which depends on a minimum number of random variables representing the stochastic field. Under the assumption of a prespecified power spectral density function of the stochastic fields, it is possible to compute the response (displacement and stress) variability of the shell structure using the direct MCS technique. It is also possible to investigate the influence of the variation of each random parameter as well as the effect of particular stochastic field characteristics, such as the correlation scale and structure, on the response variability.
2. Stochastic stiffness matrices The multi-layered triangular shell element TRIC has six Cartesian degrees of freedom per node. As it is stated in Section 1, the element formulation is based on the natural mode finite element method. A schematic representation of the element appears in Fig. 1. Its natural stiffness is based only on deformations and not on associated rigid-body motions. The element TRIC comprises 18 nodal displacements which lead to
1 side γ
2
z ψ
y ϕ θ ¸x
side α
side β
3
Fig. 1. The multi-layered triangular shell element TRIC.
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
141
12 straining modes (18 Cartesian freedoms––6 rigid-body freedoms): three axial straining modes, three symmetric bending modes, three anti-symmetric bending + shear modes and three in-plane rotations. In order to generate these pure deformational modes, which are called natural modes, a projection of the nodal displacements and rotations as well as the corresponding forces and moments on the edges of the triangle takes place. The 12 · 12 stiffness matrix kN corresponding to the natural modes is denoted as the natural stiffness matrix of the element. 2.1. Random variation of Young’s modulus and Poisson’s ratio 2.1.1. Weighted integral method As the elements of the stiffness matrix are nonlinear functions of PoissonÕs ratio m, the analysis of a shell structure involving both stochastic YoungÕs modulus and PoissonÕs ratio requires the introduction of additional approximations. However, an alternative weighted integral formulation of the stochastic stiffness matrix is possible, considering LameÕs constants k and l as the two uncertain material properties since all elements of the stiffness matrix are linear functions of k and l. For the TRIC shell element and in the case of isotropic material, the elasticity matrix j12 in the material coordinate system can be expressed as 2 3 1 m 0 E 4 j12 ¼ m 1 0 5: ð1Þ 1 m2 0 0 1m The LameÕs constants are defined as mE E : ; l¼G¼ k¼ 1 m2 2ð1 þ mÞ After some algebra, Eq. (1) becomes 2 3 2 3 2 3 k þ 2l k 0 k k 0 2l 0 0 ð1Þ ð2Þ j12 ¼ 4 k k þ 2l 0 5 ¼ 4 k k 0 5 þ 4 0 2l 0 5 ¼ j12 þ j12 : 0 0 2l 0 0 0 0 0 2l |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} j12 ðkÞ
ð2Þ
ð3Þ
j12 ðlÞ
The two material properties are assumed to vary randomly along the element surface X according to kðx; yÞ ¼ k0 ½1 þ fk ðx; yÞ;
ð4aÞ
lðx; yÞ ¼ l0 ½1 þ fl ðx; yÞ; ð4bÞ where k0 , l0 are the mean values of LameÕs constants and fk ðx; yÞ, fl ðx; yÞ are 2D-1V zero-mean homogeneous stochastic fields corresponding to the spatial variation of k and l respectively. Substituting Eqs. (4a) and (4b) into the expression of the elasticity matrix j12 we have ð1Þ
ð2Þ
j12 ¼ ½j12 0 ½1 þ fk ðx; yÞ þ ½j12 0 ½1 þ fl ðx; yÞ
ð5Þ
and the constitutive matrix jct takes now the following form ð1Þ
ð2Þ
jct ¼ jct0 þ jct0 fk ðx; yÞ þ jct0 fl ðx; yÞ or
2
jaa jct ¼ 4 jab jac
jab jbb jbc
3 jac jbc 5; jcc
where subscript zero denotes the corresponding mean value.
ð6Þ
142
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
In the following, the TRIC element stiffness terms are derived from their definitions and Eq. (6). They are finally grouped in a 12 · 12 matrix which constitutes the stochastic natural stiffness matrix of the TRIC element 2 3 0 0 kqc 6 ð66Þ 7 6 0 0 7 kqh 7: kN ¼ 6 ð7Þ 6 7 ð33Þ ð1212Þ 4 5 0 0 kaz ð33Þ
For an in-depth study of the specific features of the TRIC shell element and for any details regarding the symbols used, the reader can refer to the Appendix A of this work and to publications [1,4]. Axial and symmetric bending stiffness terms Definition: kqc ¼
Z
BTtc jct Btc dV :
ð8Þ
V
ð66Þ
Stochastic form: ð1Þ
ð2Þ
kqc ¼ kqc0 þ ak kqc0 þ al kqc0 ;
ð9Þ
ð66Þ
1 ak ¼ X
Z
fk ðx; yÞ dX;
X
1 al ¼ X
Z
fl ðx; yÞ dX:
ð10Þ
X
Anti-symmetric bending and shear stiffness terms Definition (bending): Z b kqh ¼ BTth jct Bth dV :
ð11Þ
V
ð33Þ
Stochastic form: b b ¼ kqh0 þ kqh ð33Þ
6 X
Xi DkXi þ
i¼1
6 X
ð12Þ
Yi DkYi :
i¼1
Xi , Yi are the weighted integrals corresponding to stochastic fields fk ðx; yÞ and fl ðx; yÞ respectively. Xi have the following expressions Z Z X1 ¼ f2a fk ðfa ; fb ; fc Þ dX; X2 ¼ f2b fk ðfa ; fb ; fc Þ dX; X X Z Z X3 ¼ f2c fk ðfa ; fb ; fc Þ dX; X4 ¼ fa fb fk ðfa ; fb ; fc Þ dX; ð13Þ XZ ZX X5 ¼ fb fc fk ðfa ; fb ; fc Þ dX; X6 ¼ fa fc fk ðfa ; fb ; fc Þ dX; X
X
while Yi are identical to Xi but they contain the stochastic field fl instead of fk . fa , fb , fc are the area coordinates of the triangular shell element [3]. Definition (shear): Z s BTsh Xs Bsh dV : kqh ¼ ð33Þ
V
ð14Þ
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
143
Stochastic form: s kqh ð33Þ
s ¼ ð1 þ al Þkqh0 :
ð15Þ
Finally these terms are coupled in the following way b 1 s 1 1 kqh ¼ f½kqh þ ½kqh g :
ð16Þ
ð33Þ
It is worth noting that the stochastic anti-symmetric shear terms are influenced by the average value al only. This is in accordance with the nature of these terms, since al describes the variation of the shear modulus within a finite element. Stiffness terms due to in-plane rotations (azimuth stiffness) Definition:
2
( kz ¼ X max
3 0:5 0:5 5; 1
0:5 1 0:5
1 kaz ¼ kz 4 0:5 ð33Þ 0:5 1 l2a
Z
h=2
1 z jaa dz; 2 lb h=2 2
ð17Þ Z
h=2
1 z jbb dz; 2 lc h=2 2
Z
)
h=2 2
z jcc dz :
ð18Þ
h=2
Stochastic form: ð1Þ
ð2Þ
kz ¼ kz0 þ ak kz0 þ al kz0 :
ð19Þ
In (18), h is the shell thickness. 2.1.2. Local average method As it can be seen from the previous section, the axial and symmetric bending, the anti-symmetric shear and the azimuth stiffness terms have an a priori (inherent) local average form. The expression of the antisymmetric bending stiffness terms in the context of the local average method is the following bð1Þ
bð2Þ
b b kqh ¼ kqh0 þ ak kqh0 þ al kqh0 :
ð20Þ
ð33Þ
From Eq. (12), it can be seen that the number of random variables per finite element in the weighted integral formulation of the anti-symmetric bending stiffness terms is twelve compared to the two random variables per element derived in the local average approach. This gives an indication of the computational overhead associated with the weighted integral formulation. 2.2. Random variation of thickness The consideration of random variations in shell thickness complicates significantly the application of the weighted integral method [7]. Since the thickness appears into the integrals expressing the stiffness matrix of any type of shell element, the application of the weighted integral approach for the derivation of the stochastic stiffness matrix is not straightforward. In this case, the use of the local average method is recommended. Moreover, in a recent study by Argyris et al. [3], it was shown that for cylindrical shells and random variation of YoungÕs modulus, the local average method is superior to the weighted integral method in terms of simplicity and computational efficiency, while the results obtained with the use of the local average were very close to those of the weighted integral method.
144
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
Using local averages, the randomly varying shell thickness within a finite element is expressed as hðx; yÞ ¼ h0 ð1 þ ah Þ; where 1 ah ¼ X
Z
ð21Þ
fh ðx; yÞ dX
ð22Þ
X
and fh ðx; yÞ is the stochastic field describing the variation of the shell thickness. Taking into account the definitions of the various stiffness terms as stated in Eqs. (8), (11), (14) and (17), the following stochastic terms are derived • Axial and symmetric bending stiffness terms kij ¼ ðkij Þ0 ð1 þ ah Þ;
ð23aÞ
zkij ðzkij Þ0 ð1 þ 2ah Þ;
ð23bÞ
z2 kij ðz2 kij Þ0 ð1 þ 3ah Þ: • Anti-symmetric bending and shear stiffness terms
ð23cÞ
b b kqh kqh0 ð1 þ 3ah Þ;
ð24Þ
ð33Þ s s kqh ¼ kqh0 ð1 þ ah Þ:
ð25Þ
ð33Þ
• Stiffness terms due to in-plane rotations (azimuth stiffness) kaz ð33Þ
kaz0 ð1 þ 3ah Þ:
ð26Þ
For the sake of clarity in our presentation, we indicate in the following lines how Eqs. (23a) and (23b) were derived. All other expressions may be produced in a similar way. Z h=2 N N X ð21Þ X k 0 0 kij ¼ jij dz ¼ jkij ðzk zkþ1 Þ ¼ jij ðzk z0kþ1 Þ ð1 þ ah Þ ¼ ðkij Þ0 ð1 þ ah Þ; ð23a Þ h=2
k¼1
k¼1
|fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} ðkij Þ0
zkij ¼
Z
h=2
h=2
zjij dz ¼
N N 1X ð21Þ 1 X k jkij ðz2k z2kþ1 Þ ¼ j ½ðz0 Þ2 ðz0kþ1 Þ2 ð1 þ ah Þ2 2 k¼1 2 k¼1 ij k |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ðzkij Þ0
¼ ðzkij Þ0 ð1 þ 2ah þ
a2h Þ
ðzkij Þ0 ð1 þ 2ah Þ:
0
ð23b Þ
In the above equations, subscript or superscript k indicates the layer number and subscripts i, j are equal to a, b, c indicating the edges of the triangle. It can be seen that some stiffness terms are approximate since their higher order stochastic parts have been neglected. 2.3. Random variation of Young’s modulus and thickness In this section, it is considered that YoungÕs modulus and the thickness of the shell are varying simultaneously within each finite element according to
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
Eðx; yÞ ¼ E0 ð1 þ aE Þ; hðx; yÞ ¼ h0 ð1 þ ah Þ; with
Z 1 aE ¼ fE ðx; yÞ dX; X X Z 1 fh ðx; yÞ dX; ah ¼ X X
145
ð27Þ
ð28Þ
being the local averages describing the random fluctuation of the two parameters. Following a similar procedure to that of the previous section the stochastic stiffness terms have finally the following forms: • Axial and symmetric bending stiffness terms kij ¼ ðkij Þ0 ð1 þ aE þ ah þ aE ah Þ;
ð29aÞ
zkij ðzkij Þ0 ð1 þ aE þ 2ah þ 2aE ah Þ;
ð29bÞ
z2 kij ðz2 kij Þ0 ð1 þ aE þ 3ah þ 3aE ah Þ: • Anti-symmetric bending and shear stiffness terms
ð29cÞ
b kqh ð33Þ
b kqh0 ð1 þ aE þ 3ah þ 3aE ah Þ;
ð30Þ
s kqh ð33Þ
s ¼ kqh0 ð1 þ aE þ ah þ aE ah Þ:
ð31Þ
• Stiffness terms due to in-plane rotations (azimuth stiffness) kaz ð33Þ
kaz0 ð1 þ aE þ 3ah þ 3aE ah Þ:
ð32Þ
2.4. Discussion of the statistical independence assumption between the stochastic fields In the preceding sections, the cross-correlation which could eventually exist between the random parameters considered, has been neglected. The same assumption has been made in the past by several researchers who have considered randomness in more than one material or geometric property [6,12,14,19,22]. In fact, although it is likely that a certain degree of cross-correlation might exist, especially between the YoungÕs modulus and the PoissonÕs ratio of a structure (due to the similarity of the physical mechanisms that lead to random variations of material properties), there is no available evidence regarding this issue in the literature. It is also reasonable to assume that these two material properties are not correlated with a geometric property, such as the shell thickness. The assumption made in the present work is in accordance with the results presented in a recent paper by Graham and Deodatis [11]. From the numerical tests examined in that paper, it becomes evident that the cross-correlation effect does not affect significantly the response variability when static problems are considered. This remark is not valid when a (dynamic) random eigenvalue problem is treated. In this case, the results show that the eigenvalue variability and its upper bounds are highly influenced by the degree of cross-correlation assumed between YoungÕs modulus and mass density. The substantial errors which can be introduced in random vibration problems by neglecting the cross-correlations, have been also investigated by Elishakoff in [8]. It is therefore evident that the consideration of the statistical dependence existing
146
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
between structural properties becomes crucial for the case of dynamic problems. Since in the present paper, only static shell problems are examined, the statistical independence assumption will not practically affect the accuracy of the results.
3. Computation of stress variability The main scope of this work is to calculate the response statistics of shell structures with multiple random material and geometric properties subjected to static loads. In the framework of the Monte Carlo simulation (MCS) method used here, the computation of both displacement and stress variability is straightforward even for large finite element models. A few comments regarding the calculation of the TRIC shell element stresses are following. In the natural mode method, the Cartesian strains have been replaced by the total natural strains 2 3 cta ct ¼ 4 ctb 5: ð33Þ ctc These strains are measured directly parallel to the triangleÕs edges. To the total natural axial stains ct correspond component natural stresses rc grouped in the vector 2 3 rca ð34Þ rc ¼ 4 rcb 5: rcc The natural direct stresses are for each layer ðkÞ computed from the constitutive relation rkc ¼ jkct ckt
ð35Þ
and transformed to local Cartesian stresses via r0k ¼ Brkc ;
ð36Þ
pffiffiffi T where r0 ¼ rx ry 2rxy is the local Cartesian stress vector as it is defined for the TRIC shell element. Matrix B depends only on the geometry of each element and is given in the Appendix A. In the section of numerical tests, the variability of the components of the local Cartesian stress vector (36) is examined. In addition, the variability of the equivalent Von Mises stress is calculated and compared to that of each component separately. The expression of this equivalent stress for the case of plane stress problems is used in this work. It is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi seq ¼ r2x þ r2y rx ry þ 3r2xy : ð37Þ For the analysis of composite structures the normal stress rz is also of interest since it may contribute to damage initiation. However, the material of the structures considered as test examples is supposed to be isotropic. This is the main reason for which only the variability of components rx , ry and rxy is finally examined. From Eqs. (35) and (36), it is also deduced that the variability of the stress components is based on the stochastic variation of both jct and ct , since stresses are calculated by a post-processing procedure. However, for the sake of simplicity and computational efficiency, only the contribution of ct in the random variation of the stresses of Eq. (35) is taken into account in the present work.
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
147
4. Representation of the stochastic field The generation of sample functions of the Gaussian zero-mean homogeneous stochastic fields which describe the random parameters of the structure, is done using the series of cosines formula of the spectral representation method. This method is well suited in the context of Monte Carlo simulation (MCS) technique used for calculating the response variability of stochastic structural systems [18]. The fast Fourier transform version of the method, used in many cases in order to reduce the computational burden of the simulation, is not exploited here because values of the field at non-uniformly spaced points are needed. For a 2D-1V stochastic field and for a specific simulation ðiÞ, we have [16] f ðiÞ ðx1 ; x2 Þ ¼
1 1 N 2 1 X pffiffiffi NX ð1ÞðiÞ ð2ÞðiÞ ð2Þ 2 ½Að1Þ n1 n2 cosðj1n1 x1 þ j2n2 x2 þ /n1 n2 Þ þ An1 n2 cosðj1n1 x1 j2n2 x2 þ /n1 n2 Þ;
ð38Þ
n1 ¼0 n2 ¼0
where /ðjÞðiÞ n1 n2 , j ¼ 1; 2 represent the realization for the ðiÞ simulation of the independent random phase angles uniformly distributed in the range ½0; 2p. Anð1Þ , Anð2Þ have the following expressions 1 n2 1 n2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð39aÞ ¼ 2Sff ðj1n1 ; j2n2 ÞDj1 Dj2 ; Anð1Þ 1 n2 Anð2Þ ¼ 1 n2
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2Sff ðj1n1 ; j2n2 ÞDj1 Dj2 ;
ð39bÞ
where j1n1 ¼ n1 Dj1 ; Dj1 ¼
j1u ; N1
j2n2 ¼ n2 Dj2 ;
ð40Þ
j2u ; N2
ð41Þ
Dj2 ¼
n1 ¼ 0; 1; . . . ; N1 1;
n2 ¼ 0; 1; . . . ; N2 1;
ð42Þ
Nj , j ¼ 1, 2, represent the number of intervals in which the wave number axes are subdivided and jju , j ¼ 1, 2, are the upper cut-off wave numbers which define the active region of the power spectrum Sff of the stochastic field. The last means that the power spectral density function Sff ðj1 ; j2 Þ is assumed to be zero outside the region defined by j1u 6 j1 6 j1u
and
j2u 6 j2 6 j2u :
ð43Þ
The two-sided power spectral density functions Sff used in the numerical tests section of this study, are given by " 2 2 # b1 j1 b2 j2 2 b1 b2 Sff ðj1 ; j2 Þ ¼ rf exp ; ð44aÞ 4p 2 2 r2f b1 b2 ðb21 j21 þ b22 j22 Þ exp½ðb21 j21 þ b22 j22 Þ; ð44bÞ p where rf denotes the standard deviation of the stochastic field and b1 , b2 denote the parameters that influence the shape of the spectrum and hence the scale of correlation [17]. It is reminded here that the power spectrum reflects the correlation structure of a random field in the wave number domain. Through random number generation and using Eq. (38), a large number of sample functions of the stochastic field is produced. These sample functions are generated in accordance with the pre-assigned second-order statistics of the simulated stochastic field (mean and variance). Sff ðj1 ; j2 Þ ¼
148
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
Before proceeding to the numerical tests section, a brief comment is made here regarding the choice of the distribution function (Gaussian) and of the power spectrum of the stochastic fields used in this paper. Strictly speaking, the choice of both functions is substantiated by available experimental data. If such data is not available, as in the case of the two shell structures considered as test examples, it is reasonable to resort to the aforementioned assumptions in order to perform a parametric investigation. Methods for the generation of sample functions of Gaussian stochastic fields are well established in the literature (e.g. spectral representation used in this paragraph) in contrast to the non-Gaussian case. However, with the choice of a Gaussian model for variables bounded by physical constraints (e.g. elastic constants), there is a non-zero probability that a violation of these constraints might occur in an actual MC simulation. The violation of the physical constraints is due to large negative realizations of the stochastic fields leading subsequently to possible non-positive definite stiffness matrices. This issue was handled in this paper by discarding the samples corresponding to a non-positive definite stiffness matrix without introducing a severe bias in the simulation. This is the case when the input stochastic fields are characterized by small coefficients of variation ðrf < 0:25Þ and thus, a small number of large negative outcomes is produced. For the two specific expressions of the power spectrum (and the underlying autocorrelation function), the square exponential type is adopted. This follows from the observation by Li and Der Kiureghian [13] who have shown that the error in the discretization (e.g. by local average method) of a homogeneous Gaussian stochastic field is negligible even for small correlation lengths when a square exponential type of correlation is adopted.
5. Numerical tests Two shell structures of different shape, denoted as Tests 1 and 2, are used in this section in order to demonstrate the proposed methodologies. Test 1 is the Scordelis-Lo shell roof, a cylindric shell subjected to gravity loading. The geometrical and material properties of the structure are presented in Fig. 2. It is assumed that the two longitudinal edges are free and the two circular edges are supported by diaphragms. The structure is discretized with 2 · 28 · 28 ¼ 1568 shell elements resulting in 4872 degrees of freedom (d.o.f.). Test 2 is a simply supported hyperboloid shell of revolution. The meridional line of the shell is described by the following hyperbola equation R1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi RðyÞ ¼ ð45Þ c2 þ y 2 ; c where c is a real constant and R1 the radius at the top of the shell. The structure is loaded by a single concentrated load of 10 kN acting in the x-direction at l ¼ 2=3L, with L ¼ 20 m the total height of the shell. The top and bottom radii are 4.8 and 8.0 m, respectively. Geometry and material data are shown in Fig. 3. A finite element mesh consisting of 1024 triangular elements and resulting in 3306 d.o.f. is used for the discretization of the shell. The spatial variability in YoungÕs modulus and thickness of the shells is described by two uncorrelated 2D-1V zero-mean homogeneous stochastic fields with coefficient of variation equal to 0.1, while the applied loads are considered to be deterministic. The same assumption is made for the stochastic fields describing the random variation of YoungÕs modulus and PoissonÕs ratio. The power spectral density function Sff characterizing the two stochastic fields in both cases is assumed to have the form of Eq. (44a) or (44b). It should be mentioned that the results appearing in Figs. 4, 5, 7 and 8 are based on the power spectrum of Eq. (44a), while the power spectrum of Eq. (44b) is used only for comparison purposes in Fig. 6. In both examples, the coefficient of variation (COV) of displacement wc at the characteristic node C of the structure, is calculated using the local average and, when feasible, the weighted integral method in conjunction with MCS. The COV of a response quantity, e.g. displacement ui , is defined as
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
149
h rigid diaphragm
q D rigid diaphragm
C A
L B
z y R x
θ
R = 7600 mm L = 15200 mm h = 76 mm θ = 40˚ 2 E = 21000 N/mm ν = 0.3
Fig. 2. Material and geometric data of the Scordelis-Lo shell roof.
Fig. 3. Material and geometric data of the hyperboloid shell.
COVðui Þ ¼
rðui Þ ; Eðui Þ
ð46Þ
where r and E are the standard deviation and mean value of ui , respectively. The response statistics have been computed using a sample size of 500. This number of samples can represent sufficiently well the preassigned first two moments of the stochastic field and statistical convergence of the results is achieved within this number of simulations. The sensitivity of structural response with regard to the scale of correlation of the stochastic fields, quantified via the correlation length parameter b is first examined. For this purpose, several sets of sample functions are generated using Eq. (38) each one for a different value of parameter b. These values are related to the dimensions of the structures along the x- and y-axes. For the Scordelis-Lo shell one value b ¼ b1 ¼ b2 is considered to be representative of both dimensions. For the hyperboloid shell, where the difference
150
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160 0.3 E h E,h
0.2
COV of w at node C
COV of w at node C
0.25
0.15 0.1
0.25 E
0.2
h E,h
0.15 0.1
0.05
0.05 0
0 0.065
0.65
(a)
1.3
2.6
26
65
0.2
130
2
(b)
Correlation length parameter b
4
8
80
200
Correlation length parameter b1
Fig. 4. (a) Scordelis-Lo shell: COV of vertical deflection wc as a function of correlation length parameter b (local average method). (b) Hyperboloid shell: COV of displacement wc as a function of correlation length parameter b1 (local average method).
0.1
COV of w at node C
COV of w at node C
0.15
0.1
0.05
0.05
0 0.065
0.65
1.3
2.6
26
65
0
130
0.2
Correlation length parameter b
(a)
E
E,v (Local average)
E,v (Weighted integral)
2
4
8
80
200
Correlation length parameter b1
(b)
E
E (Weighted integral)
E,v (Local average)
E,v (Weighted integral)
Fig. 5. (a) Scordelis-Lo shell: comparison between local average and weighted integral methods––COV of vertical deflection wc as a function of correlation length parameter b for the case of combined YoungÕs modulus and PoissonÕs ratio variation. (b) Hyperboloid shell: comparison between local average and weighted integral methods––COV of displacement wc as a function of correlation length parameter b1 for the case of combined YoungÕs modulus and PoissonÕs ratio variation.
between the two dimensions is greater, two values b1 6¼ b2 are used, each one related to the corresponding dimension. Four different cases of randomly varying parameters are considered: (i) stochastic YoungÕs modulus, (ii) stochastic thickness, (iii) combined YoungÕs modulus and thickness, and (iv) combined YoungÕs modulus and PoissonÕs ratio. The local average method is used in the first three cases for reasons explained in Section 2.2, while in the fourth case, a comparison between the local average and weighted integral methods is performed. The results concerning the displacement variability are depicted in Fig. 4a and b for the Scordelis-Lo and the hyperboloid shell cases, respectively. From these plots, it can be seen that, in both shell geometries, the displacement variability shows similar trends, starting from small values for small correlation lengths corresponding to white noise stochastic fields, up to large values for large correlation lengths. Especially, when YoungÕs modulus is the only random parameter of the problem, the COV of the selected displacement tends to the COV of YoungÕs modulus for large values of parameter b, as expected. It is also evident that randomness in shell thickness has significant effect on the response variability, compared to the effect of random variation of YoungÕs modulus [7]. When both quantities are supposed to vary randomly, the response COV tends to values that are 2–2.5 times greater than the corresponding input COV.
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160 0.25
0.2
COV of w at node C
COV of w at node C
0.25
151
h h(P.s.2)
0.15 0.1 0.05
0.2 E,h E,h(P.s.2)
0.15 0.1 0.05 0
0 0.065
(a)
0.65
1.3
2.6
26
65
0.065
130
(b)
Correlation length parameter b
0.65
1.3
2.6
26
65
130
Correlation length parameter b
0.3 0.35 0.3
COV of w at node C
COV of w at node C
0.25 0.2 0.15 h h(P.s.2)
0.1 0.05
0.2 0.15 0.1 0.05 0
0 0.2
(c)
E,h E,h(P.s.2)
0.25
2
4
8
80
Correlation length parameter b1
200
0.2
(d)
2
4
8
80
200
Correlation length parameter b1
Fig. 6. Scordelis-Lo shell: comparison between the two power spectra––COV of vertical deflection wc as a function of correlation length parameter b for the case of (a) thickness variation, (b) combined YoungÕs modulus and thickness variation. Hyperboloid shell: comparison between the two power spectra––COV of displacement wc as a function of correlation length parameter b1 for the case of (c) thickness variation, (d) combined YoungÕs modulus and thickness variation.
The very small effect of the PoissonÕs ratio stochastic variation is evidenced by Fig. 5a and b. The displacement variability for the case of randomness in both LameÕs constants is approximately the same or even slightly smaller compared to the case of YoungÕs modulus stochasticity. Similar conclusions have been reported in [22] for plane stress and plate problems. A comparison of the local average and weighted integral methods is performed next. For the ScordelisLo shell roof, these differences are negligible as it can be seen in Fig. 5a, while the discrepancy becomes more pronounced in the hyperboloid shell example (Fig. 5b). This can be explained by the fact that a relatively coarse mesh (compared to that used for the Scordelis-Lo shell) is used for the hyperboloid shell. This behavior confirms the observation that a relatively fine mesh is required by the local average method in order to produce accurate results. Furthermore, the effect of type of correlation, expressed by the mathematical form of the power spectrum, is investigated using the expressions (44a) and (44b). Two cases of random variation are considered for each structure: (i) stochastic thickness, (ii) stochastic YoungÕs modulus and thickness. The results are depicted in Fig. 6a–d. It is observed that, in both shell cases, some differences exist between the COVs calculated using the two spectra. For the hyperboloid shell example these differences are substantial in the whole range of b examined, while for the Scordelis-Lo shell the discrepancy becomes more pronounced as the correlation length parameter increases. This observation can be justified by the fact that the variance of a random response quantity (in this case the selected displacement) is mathematically related to the power spectral density function of the stochastic field describing the random input quantities (in this case the YoungÕs modulus and the thickness) [11]. It is therefore reasonable that spectra, having their spectral power concentrated over different wave number regions, produce slightly different response variations.
152
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
In Tables 1 and 2, a complete image of the response statistics of the two shells under consideration is presented, using the local average approach. The mean value, standard deviation and COV of displacement wc at the characteristic node C of the structures are calculated for the following cases: (i) in a ‘‘mixed’’ way, i.e. when the two random parameters are introduced simultaneously in the model and (ii) approximately, using the SRSS concept. SRSS denotes the square root of the sum of squares of the two outputs by two random parameters, E and h, considered separately. The assumption of independence between the stochastic fields describing the two properties is well reflected by the first four entries of the sixth column of the tables, where the ratio of mean values is approximately equal to the square root of 2. It is also evident that the SRSS formula can offer an accurate prediction of the results produced by the combined variation case [6]. Finally, it can be seen once again, that randomness in shell thickness has a significant effect on the displacement variability, compared to the effect of random variation in YoungÕs modulus. Up to this moment, the displacement variability of the two shell structures has been considered. In the context of this study, similar investigations have been performed with regard to the stress variability. The stress components examined are those described in Section 3. The influence of the correlation length parameter b on the stress variability is first tested and the corresponding plots of stresses inside some specific elements are depicted in Figs. 7 and 8. By comparing these plots with those of the displacement variability, different trends of variation can be observed among several stress components. It can be also seen that the ‘‘distance’’ separating the results produced by random thickness variation and by combined YoungÕs modulus and thickness variation, is greater for the stress variability case than for the case of displacement variability. In order to further examine the variability of stress components within several finite elements, Tables 3 and 4 are presented. The results contained in these tables correspond to a specific value of correlation length parameter b for each structure (b ¼ 65:0 for the Scordelis-Lo shell and b1 ¼ 200:0 for the hyperboloid shell). It is clear that the fluctuation is substantially different not only between several stress components but even between the same stress component within different elements. This variation is more pronounced for the case of random thickness fluctuation which indicates how important may be the calculation of probability- and spectral-distribution-free upper bounds of the response (displacement and stress) variability.
Table 1 Scordelis-Lo shell: displacement variability as a function of correlation length parameter b (local average method and power spectrum of Eq. (44a)) Items
Correlation length parameter b
Random parameters YoungÕs modulus E
Thickness h
E and h (mixed)
SRSS/ mixed
Mean value
0.065 0.65 2.6 65
41.8757 41.9306 42.0061 41.8467
42.1692 42.3711 42.5307 43.212
42.2713 42.5249 42.7531 43.2767
1.4059 1.4018 1.3982 1.39
Standard deviation
0.065 0.65 2.6 65
0.2046 0.4701 1.6144 3.9851
0.2168 0.8473 3.5007 8.4818
0.2904 0.9834 3.8937 9.3016
1.0265 0.9853 0.9901 1.0075
COV
0.065 0.65 2.6 65
0.00488 0.01121 0.03843 0.09523
0.00514 0.01999 0.08231 0.19628
0.00687 0.02312 0.09107 0.21493
1.0317 0.9913 0.9975 1.015
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
153
Table 2 Hyperboloid shell: displacement variability as a function of correlation length parameter b1 (local average method and power spectrum of Eq. (44a)) Items
Correlation length parameter b1
Random parameters YoungÕs modulus E
Thickness h
E and h (mixed)
Mean value
0.2 2 8 80
4.786 4.7949 4.8214 4.7804
4.8337 4.8798 4.9828 5.0057
4.8441 4.8974 5.0305 5.0089
1.4042 1.3969 1.3783 1.3819
Standard deviation
0.2 2 8 80
0.0546 0.1443 0.2833 0.4486
0.1475 0.4036 0.8697 1.2018
0.1565 0.4189 0.9357 1.2655
1.005 1.0232 0.9775 1.0137
COV
0.2 2 8 80
0.01141 0.03009 0.05877 0.09384
0.03052 0.0827 0.17453 0.24009
0.03232 0.08554 0.18601 0.25266
1.0081 1.0288 0.9901 1.0203
0.25
SRSS/mixed
0.25 σy(E)
σx(h)
0.2
COV of σy at point A
COV of σx at point A
σx(E) σx(E,h)
0.15 0.1 0.05
σy(h)
0.2
σy(E,h)
0.15 0.1 0.05 0
0 0.065
0.65
(a)
1.3
2.6
26
65
0.065
130
(b)
Correlation length parameter b
2.6
26
65
130
65
130
τeq(E)
σxy(h)
0.4
σxy(E,h)
0.3 0.2 0.1
COV of τeq at point A
σxy(E)
COV of σxy at point A
1.3
Correlation length parameter b
0.25
0.5
τeq(h)
0.2
τeq(E,h)
0.15 0.1 0.05
0
0 0.065
(c)
0.65
0.65
1.3
2.6
26
Correlation length parameter b
65
130
0.065
(d)
0.65
1.3
2.6
26
Correlation length parameter b
Fig. 7. (a) Scordelis-Lo shell: COV of stress rx at point A as a function of correlation length parameter b (local average method). (b) Scordelis-Lo shell: COV of stress ry at point A as a function of correlation length parameter b (local average method). (c) Scordelis-Lo shell: COV of stress rxy at point A as a function of correlation length parameter b (local average method). (d) Scordelis-Lo shell: COV of equivalent stress seq at point A as a function of correlation length parameter b (local average method).
154
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160 0.2
σy(E)
σx(E)
0.2
COV of σy at point A
COV of σx at point A
0.25
σx(h) σx(E,h)
0.15 0.1 0.05 0
(a)
2
4
8
80
0.1
0.05
0.2
200
(b)
Correlation length parameter b1 0.2
2
4
8
80
200
Correlation length parameter b1 0.2
σxy(E)
τeq(E)
σxy(h)
0.15
COV of τeq at point A
COV of σxy at point A
σy(E,h)
0 0.2
σxy(E,h)
0.1
0.05
τeq(h)
0.15
τeq(E,h)
0.1
0.05
0
0 0.2
(c)
σy(h)
0.15
2
4
8
80
Correlation length parameter b1
0.2
200
(d)
2
4
8
80
200
Correlation length parameter b1
Fig. 8. (a) Hyperboloid shell: COV of stress rx at point A as a function of correlation length parameter b1 (local average method). (b) Hyperboloid shell: COV of stress ry at point A as a function of correlation length parameter b1 (local average method). (c) Hyperboloid shell: COV of stress rxy at point A as a function of correlation length parameter b1 (local average method). (d) Hyperboloid shell: COV of equivalent stress seq at point A as a function of correlation length parameter b1 (local average method).
6. Summary and conclusions In this work, a stochastic finite element analysis of shell structures with combined uncertain material and geometric properties is performed. For this purpose, a stochastic formulation of the triangular composite facet shell element TRIC is derived assuming random variation of the YoungÕs modulus, the PoissonÕs ratio and the thickness of the shell. As a result of the proposed formulation and the special features of the element, the stochastic stiffness matrix of TRIC depends finally on a minimum number of random variables representing the stochastic field. This fact leads to a cost-effective stiffness matrix, which is very important in the case of a time consuming stochastic analysis of real world structures. The spatial fluctuation of the mechanical and geometric properties is described by uncorrelated stochastic fields. Under the assumption of a pre-specified power spectral density function for these stochastic fields, it is possible to compute the response variability of a realistic shell structure using the direct MCS technique. Two shells of cylindrical and hyperboloid shape are tested. The sensitivity of response statistics with regard to the scale of correlation of the stochastic fields, quantified via the correlation length parameter b is first examined. For both shell geometries, the displacement variability shows similar trends, starting from small values for small correlation lengths, up to large values for large correlation lengths. When the YoungÕs modulus is the only random parameter of the problem, the COV of the selected displacement tends to the COV of YoungÕs modulus for large values of parameter b, as expected. It is also evident that random variation in the shell thickness has significant effect on the displacement variability compared to the effect of
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
155
Table 3 Scordelis-Lo shell: stress variability for b ¼ 65:0 (local average method and power spectrum of Eq. (44a)) Items
Stress component, point, layer
Random parameters YoungÕs modulus E
Thickness h
E and h (mixed)
SRSS/mixed
Mean (absolute) value
rx , B, 2 rxy , B, 2 ry , B, 1 ry , A, 2 rx , D, 1 ry , D, 2 rx , C, 1
2.0766 0.2837 15.9683 1.2376 0.4926 7.7248 1.3981
2.1221 0.2393 16.2278 1.2724 0.5081 7.9053 1.3835
2.124 0.2399 16.2437 1.2739 0.5086 7.9134 1.3864
1.3979 1.5471 1.4016 1.3934 1.3914 1.3967 1.4187
Standard deviation
rx , B, 2 rxy , B, 2 ry , B, 1 ry , A, 2 rx , D, 1 ry , D, 2 rx , C, 1
0.1952 0.0277 1.5013 0.1169 0.0465 0.7301 0.1316
0.3251 0.1585 2.1852 0.2589 0.0963 1.2837 0.1429
0.3784 0.1562 2.659 0.2843 0.1063 1.4752 0.2008
1.0021 1.0301 0.9971 0.9992 1.006 1.0011 0.9675
COV
rx , B, 2 rxy , B, 2 ry , B, 1 ry , A, 2 rx , D, 1 ry , D, 2 rx , C, 1
0.094 0.0975 0.094 0.0945 0.0945 0.0945 0.0941
0.1532 0.6623 0.1347 0.2035 0.1895 0.1624 0.1033
0.1782 0.6511 0.1637 0.2231 0.2089 0.1864 0.1448
1.0086 1.0282 1.0034 1.0057 1.0137 1.008 0.965
Table 4 Hyperboloid shell: stress variability for b1 ¼ 200:0 (local average method and power spectrum of Eq. (44a)) Items
Stress component, point, layer
Random parameters YoungÕs modulus E
Thickness h
E and h (mixed)
SRSS/mixed
Mean (absolute) value
rx , D, 2 rxy , D, 2 ry , D, 1 rx , B, 1 ry , B, 2 rxy , B, 2 ry , A, 2
0.0363 0.0272 0.0485 0.1299 0.0278 0.0261 0.2806
0.0361 0.0275 0.0513 0.1298 0.0278 0.0253 0.2857
0.0362 0.0276 0.0514 0.1309 0.0279 0.0254 0.2862
1.4142 1.4014 1.3735 1.4029 1.4091 1.4311 1.3992
Standard deviation
rx , D, 2 rxy , D, 2 ry , D, 1 rx , B, 1 ry , B, 2 rxy , B, 2 ry , A, 2
0.0034 0.0026 0.0046 0.0123 0.0026 0.0025 0.0267
0.0032 0.0037 0.0142 0.0116 0.0031 0.0027 0.0403
0.0048 0.0045 0.0148 0.0166 0.0042 0.0037 0.0485
0.9727 1.0049 1.0085 1.0185 0.9633 0.9945 0.9967
COV
rx , D, 2 rxy , D, 2 ry , D, 1 rx , B, 1 ry , B, 2 rxy , B, 2 ry , A, 2
0.0937 0.0955 0.0948 0.0947 0.0935 0.0958 0.0952
0.0893 0.1345 0.2768 0.0894 0.1115 0.1067 0.1411
0.1326 0.1631 0.2879 0.1268 0.1505 0.1457 0.1695
0.9762 1.0114 1.0163 1.0271 0.9669 0.9842 1.0042
156
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
random YoungÕs modulus. When both quantities are supposed to vary randomly, the response COV tends to values that are 2–2.5 times greater than the input COV. In the case of combined fluctuation of YoungÕs modulus and PoissonÕs ratio, the displacement variability curve indicates the small effect of the PoissonÕs ratio on the results. Slight differences are found between the results given by the local average method and those of the weighted integral method in the case of Scordelis-Lo shell. These differences become somewhat more pronounced in the hyperboloid shell case. In addition, the effect of type of correlation, expressed in this work by the mathematical form of the power spectrum, is investigated. Due to the relation existing between the power spectral density function of the stochastic field describing the random input quantities and the variance of a response quantity, slight differences are once again observed. Finally, investigations similar to the aforementioned were made with regard to the stress variability. The main observation in this case is that the stress fluctuation is substantially different not only between several stress components but even between the same stress component within different elements. This variation takes its maximum value when the YoungÕs modulus and the shell thickness are simultaneously varying. Acknowledgement This work was partially supported by the research project ‘‘Thalis’’ of the National Technical University of Athens. Appendix A. Stochastic stiffness matrix components In Eq. (8) of the present work, where the definition for the axial and symmetric bending stiffness terms of the TRIC shell element is given, Btc is a matrix of the form 3 2 z 0 0 1 0 0 7 6 la 7 6 60 1 0 0 z 0 7 ðA:1Þ Btc ¼ 6 7; lb 7 6 ð36Þ 4 z5 0 0 1 0 0 lc where la , lb , lc are the lengths of the triangleÕs edges. In Eq. (9), the mean value of matrix kqc is defined as 2 zkaa zkab zkac 3 kab kac kaa 6 la lb lc 7 7 6 7 6 zk zk zk ab bb bc 7 6 kbb kbc 6 la lb lc 7 7 6 7 6 zkac zkbc zkcc 7 6 7 6 kcc 6 la lb lc 7 7 6 ; kqc0 ¼ X 6 z2 kaa z2 kab z2 kac 7 7 6 symm: ð66Þ 7 6 l2a la lb la lc 7 6 7 6 2 6 z kbb z2 kbc 7 7 6 6 lb lc 7 l2b 7 6 7 6 4 z2 kcc 5 l2c
ðA:2Þ
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
157
where kij ¼
Z
h=2
jij dz ¼
h=2
zkij ¼
Z
h=2
zjij dz ¼
Z
jkij ðzk zkþ1 Þ;
i; j ¼ a; b; c;
ðA:3Þ
k¼1
h=2
z2 kij ¼
N X
N 1X jk ðz2 z2kþ1 Þ; 2 k¼1 ij k
h=2
z2 jij dz ¼
h=2
N 1X jk ðz3 z3kþ1 Þ; 3 k¼1 ij k
i; j ¼ a; b; c;
i; j ¼ a; b; c;
ðA:4Þ
ðA:5Þ
X is the area of the triangle, h is the element thickness, N the number of layers and zk the coordinate of the layer k as shown in Fig. 9. ð1Þ ð2Þ Also in Eq. (9), matrices kqc0 , kqc0 are identical in form with matrix kqc0 but they depend only on k and l, respectively. b In the definition of the anti-symmetric bending stiffness terms kqh (Eq. (11)) Bth ¼ zl2 Ah l; l is a diagonal 2 la l¼40 0
ðA:6Þ
matrix, which contains the lengths of the triangleÕs edges 3 0 0 lb 0 5 0 lc
and matrix Ah is a function of the triangular or natural coordinates fa , fb , fc [3] 2 3 fa fa 3ðfc fb Þ 5: 3ðfa fc Þ fb Ah ¼ 4 fb fc fc 3ðfb fa Þ b The mean value of matrix kqh is defined as 2 b11 3 b12 b13 kqh kqh kqh 6 7 b b22 b23 7 kqh kqh ¼6 kqh0 4 5: b33 symm: kqh
ðA:7Þ
ðA:8Þ
ðA:9Þ
b The entries of matrix kqh0 depend on the geometric characteristics of the TRIC element (area and lengths of the triangleÕs edges) and on the stiffness quantities defined in (A.5) [3].
Fig. 9. Coordinates z of the layers.
158
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
All matrices Dk belonging to the fluctuating part of the stiffness terms given by Eq. (12) are deterministic and have the following forms 2 DkX1
6 6 6 ¼6 6 4
0
l2b 2 ð1Þ zk l4a aa
0 þ
3
0
6 2 ð1Þ 9 2 ð1Þ zk þ zk l2a ab l2b bb
symm:
lb lc 2 ð1Þ 3lc 2 ð1Þ 3lb 2 ð1Þ 9 2 ð1Þ 7 z kaa 2 z kab 2 z kac zk 7 4 lb lc bc 7 la la l b la lc 7; 7 l2c 2 ð1Þ 6 2 ð1Þ 9 2 ð1Þ 5 z k þ z k þ z k l2c cc l4a aa l2a ac ðA:10Þ
2
DkX2
9 2 ð1Þ 6 2 ð1Þ l2a 2 ð1Þ 6 l2 z kaa þ l2 z kab þ l4 z kbb b b 6 a ¼6 6 4 symm: 2
DkX3
9 2 ð1Þ 6 2 ð1Þ l2a 2 ð1Þ 6 l2 z kaa þ l2 z kac þ l4 z kcc 6 a c c 6 ¼6 6 4 symm: 2
DkX4
60 6 6 6 ¼6 6 6 4 symm:
10 2 ð1Þ zk la lb ab 0
0 0
3 3lc 2 ð1Þ la lc 2 ð1Þ 9 2 ð1Þ 3la 2 ð1Þ z k z k z k z k bb la lc ac la l2b ab l4b l2b lc bc 7 7 7; 0 7 2 5 lc 2 ð1Þ 6 2 ð1Þ 9 2 ð1Þ z k þ z k þ z k bb bc cc 4 2 2 lc lb lb
9 2 ð1Þ 3la 2 ð1Þ 3lb 2 ð1Þ la lb 2 ð1Þ zk zk z k 4 z kcc la lb ab lb l2c bc la l2c ac lc 9 2 ð1Þ 6 2 ð1Þ l2b 2 ð1Þ z k þ z k þ 4 z kcc lc l2b bb l2c bc
2 DkX6
6 6 6 ¼6 6 4
0
symm:
07 7 7 7; 07 5
ðA:12Þ
0
2
DkX5
3
3
lc 2 ð1Þ 9 2 ð1Þ 3la 2 ð1Þ zk þ z k þ 2 z kbc 7 la l2b ab la lc ac lb lc 7 7 lc 2 ð1Þ 3lb 2 ð1Þ 9 2 ð1Þ 7 z kab þ 2 z kac þ z kbc 7; 2 lb lc la lb l a lc 7 2l2c 2 ð1Þ 6 2 ð1Þ 6 2 ð1Þ 18 2 ð1Þ 7 5 2 2 z kab 2 z kac 2 z kbc 2 z kcc la lc lb la lb
18 2 ð1Þ 6 2 ð1Þ 6 2 ð1Þ 2l2a 2 ð1Þ 6 l2 z kaa l2 z kab l2 z kac l2 l2 z kbc 6 a c b b c 6 9 2 ð1Þ 3lb 2 ð1Þ la 2 ð1Þ 6 zk þ z k þ 2 z kbc ¼6 la lb ab la l2c ac lb lc 6 6 4 3lc 2 ð1Þ la 2 ð1Þ 9 2 ð1Þ z kab þ 2 z kbc þ zk 2 la lc ac la lb lb lc
ðA:11Þ
0
10 2 ð1Þ zk lb lc bc
lb 2 ð1Þ 9 2 ð1Þ 3la 2 ð1Þ zk þ zk þ zk la lb ab lb l2c bc la l2c ac 2l2b 6 18 ð1Þ 6 ð1Þ ð1Þ ð1Þ 2 z2 kab 2 2 z2 kac 2 z2 kbb 2 z2 kbc la lc la lc lb
3 symm: 7 7 7 7 7; 7 7 5 0
ðA:13Þ
ðA:14Þ
3 10 2 ð1Þ z kac 7 la lc 7 7 3lc 2 ð1Þ lb 2 ð1Þ 9 2 ð1Þ 7: 7 z k þ z k þ z k lb lc bc 5 l2a lb ab l2a lc ac 0
ðA:15Þ Matrices DkYi are identical to DkXi , with superscript (2) in the place of (1). Superscript (1) means dependence on k, while superscript (2) means dependence on l.
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
159
The mean value of the stochastic anti-symmetric shear stiffness terms given by Eq. (14) is contained in a matrix of the form 2
s kqh ¼
1 4 12
K111 þ K322
K112 11 K2 þ K122
symm:
3 K312 5: K212 11 K3 þ K222
ðA:16Þ
The entries of this matrix are functions of the geometry and of the mean shear modulus of the shell [3]. Finally matrix B, which appears in Eq. (36) and relates the local Cartesian stresses to the direct natural stresses, has the following expression 2
c2ax0 6 c2 B ¼ 4 bx0 c2cx0
s2ax0 s2bx0 s2cx0
3 2cax0 sax0 2cbx0 sbx0 7 5; 2ccx0 scx0
ðA:17Þ
where ax0 , bx0 , cx0 are the angles that the triangleÕs edges a, b, c form with the local x0 axis respectively and s, c are the sines and cosines of these angles.
References [1] J. Argyris, M. Papadrakakis, C. Apostolopoulou, S. Koutsourelakis, The TRIC shell element: theoretical and numerical investigation, Comput. Methods Appl. Mech. Engrg. 182 (2000) 217–245. [2] J. Argyris, M. Papadrakakis, Z.S. Mouroutis, Nonlinear dynamic analysis of shells with the triangular element TRIC, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3005–3038. [3] J. Argyris, M. Papadrakakis, G. Stefanou, Stochastic finite element analysis of shells, Comput. Methods Appl. Mech. Engrg. 191 (2002) 4781–4804. [4] J. Argyris, L. Tenek, L. Oloffson, TRIC: a simple but sophisticated 3-node triangular element based on 6 rigid-body and 12 straining modes for fast computational simulations of arbitrary isotropic and laminated composite shells, Comput. Methods Appl. Mech. Engrg. 145 (1997) 11–85. [5] B. Bhattacharyya, S. Chakraborty, NE-MCS technique for stochastic structural response sensitivity, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5631–5645. [6] C.K. Choi, H.C. Noh, Stochastic analysis of shape imperfections in RC cooling tower shells, J. Struct. Engrg. (ASCE) 126 (2000) 417–423. [7] C.K. Choi, H.C. Noh, Stochastic finite element analysis of plate structures by weighted integral method, Struct. Engrg. Mech. 4 (1996) 703–715. [8] I. Elishakoff, On the role of cross-correlations in random vibrations of shells, J. Sound Vibr. 50 (1977) 239–252. [9] G. Falsone, N. Impollonia, A new approach for the stochastic analysis of finite element modeled structures with uncertain parameters, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5067–5085. [10] R.G. Ghanem, R.M. Kruger, Numerical solution of spectral stochastic finite element systems, Comput. Methods Appl. Mech. Engrg. 129 (1996) 289–303. [11] L.L. Graham, G. Deodatis, Response and eigenvalue analysis of stochastic finite element systems with multiple correlated material and geometric properties, Probabilist. Engrg. Mech. 16 (2001) 11–29. [12] L.L. Graham, E.F. Siragy, Stochastic finite element analysis for elastic buckling of stiffened panels, J. Engrg. Mech. (ASCE) 127 (2001) 91–97. [13] C.C. Li, A. Der Kiureghian, Optimal discretization of random fields, J. Engrg. Mech. (ASCE) 119 (1993) 1136–1154. [14] P.L. Liu, A. Der Kiureghian, Finite element reliability of geometrically nonlinear uncertain structures, J. Engrg. Mech. (ASCE) 117 (1991) 1806–1825. [15] H.G. Matthies, C.E. Brenner, C.G. Bucher, C. Guedes Soares, Uncertainties in probabilistic numerical analysis of structures and solids––stochastic finite elements, Struct. Safety 19 (1997) 283–336. [16] M. Shinozuka, G. Deodatis, Simulation of multi-dimensional Gaussian stochastic fields by spectral representation, Appl. Mech. Rev. (ASME) 49 (1996) 29–53.
160
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
[17] C.A. Schenk, G.I. Schu€eller, Buckling analysis of cylindrical shells with random geometric imperfections, Int. J. Non-Linear Mech. 38 (2003) 1119–1132. [18] P.D. Spanos, B.A. Zeldin, Monte Carlo treatment of random fields: a broad perspective, Appl. Mech. Rev. (ASME) 51 (1998) 219–237. [19] B. Van den Nieuwenhof, J.P. Coyette, Modal approaches for the stochastic finite element analysis of structures with material and geometric uncertainties, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3705–3729. [20] E. Vanmarcke, Random Fields: Analysis and Synthesis, The MIT Press, Cambridge, MA, 1983. [21] J. Zhang, B. Ellingwood, Orthogonal series expansions of random fields in reliability analysis, J. Engrg. Mech. (ASCE) 120 (1994) 2660–2677. [22] W.Q. Zhu, Y.J. Ren, W.Q. Wu, Stochastic FEM based on local averages of random vector fields, J. Engrg. Mech. (ASCE) 118 (1992) 496–511.
Stochastic finite element analysis of shells with combined random material and geometric properties G. Stefanou, M. Papadrakakis
*
Institute of Structural Analysis and Seismic Research, National Technical University, Athens 15780, Greece Received 12 June 2003; received in revised form 24 September 2003; accepted 3 October 2003
Abstract In this work, an extended stochastic formulation of the triangular composite facet shell element TRIC is presented for the case of combined uncertain material (YoungÕs modulus, PoissonÕs ratio) and geometric (thickness) properties. These properties are assumed to be described by uncorrelated two-dimensional homogeneous stochastic fields. The stochastic finite element analysis of shell structures is performed using the spectral representation method for the description of the random fields in conjunction with Monte Carlo simulation (MCS) for the computation of the response variability. Useful conclusions regarding the influence of each one of the structural parameters on the response variability are derived from the numerical tests examined. 2003 Elsevier B.V. All rights reserved. Keywords: Shell finite element; Stochastic analysis; Response variability; Spectral representation
1. Introduction The analysis of shells presents a challenge, since their formulation may become cumbersome and their behavior can be unpredictable with regard to geometry or support conditions. For this reason, shells have been considered as the ‘‘prima donna of structures’’ in the sense that their performance depends very much on how they are designed and how they are treated. The extreme sensitivity of thin shells to imperfections in material, geometry and boundary conditions can be described reasonably well, if the assumption of a random fluctuation of these imperfections is induced in the analysis, in other words if a stochastic finite element analysis is performed. The need for a robust, accurate and computationally efficient shell element becomes even greater for the computationally expensive stochastic finite element analysis of large realistic structures. The TRIC shell element [4] is a simple but sophisticated triangular, shear-deformable facet shell element suitable for the analysis of thin and moderately thick isotropic as well as composite plate and shell structures. Its formulation is based on the natural mode finite element method, which has a number of
*
Corresponding author. Tel.: +301-772-1694; fax: +301-772-1693. E-mail addresses: [email protected] (G. Stefanou), [email protected] (M. Papadrakakis).
0045-7825/$ - see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.10.001
140
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
computational advantages compared to the conventional isoparametric finite element formulations. The treatment of the element kinematics (inclusion of the transverse shear deformations in its formulation based on a first-order shear-deformable beam theory) eliminates locking phenomena in a physical manner. Several investigations have confirmed the rigorous theoretical basis of the element, while numerical examples have verified its locking free properties as well as its computational efficiency in linear and nonlinear structural applications [1,2,4]. The most important advantage of the TRIC element from the stochastic analystÕs point of view is the fact that there is no need to perform numerical integration for the computation of its deterministic stiffness matrix, which is carried out in closed form. As a result of this special feature, the TRIC element provides an ideal basis for the formulation of a computationally efficient stochastic stiffness matrix and for the use of the element in large-scale finite element stochastic computations. The majority of work done in the context of stochastic finite elements has been based on the assumption that only one material or geometric property is described by a stochastic field, e.g. [3,5,9,10,17,21]. For example, the assumption is often made that YoungÕs modulus is randomly varying over space, while PoissonÕs ratio is a deterministic constant. However, many of the physical mechanisms that lead to random variations of YoungÕs modulus also lead to random variations in other material properties such as PoissonÕs ratio. In this work, the stochastic stiffness matrix of the TRIC shell element is computed taking into account random variations of YoungÕs modulus and PoissonÕs ratio as well as the thickness of the shell. These random parameters are assumed to be uncorrelated. They are described by independent two-dimensional univariate (2D-1V) homogeneous Gaussian stochastic fields, which are represented via the spectral representation method. The stochastic stiffness matrix of the shell element is derived in the form of local average [15,20] and weighted integral [11,15] methods leading to a cost-effective matrix, which depends on a minimum number of random variables representing the stochastic field. Under the assumption of a prespecified power spectral density function of the stochastic fields, it is possible to compute the response (displacement and stress) variability of the shell structure using the direct MCS technique. It is also possible to investigate the influence of the variation of each random parameter as well as the effect of particular stochastic field characteristics, such as the correlation scale and structure, on the response variability.
2. Stochastic stiffness matrices The multi-layered triangular shell element TRIC has six Cartesian degrees of freedom per node. As it is stated in Section 1, the element formulation is based on the natural mode finite element method. A schematic representation of the element appears in Fig. 1. Its natural stiffness is based only on deformations and not on associated rigid-body motions. The element TRIC comprises 18 nodal displacements which lead to
1 side γ
2
z ψ
y ϕ θ ¸x
side α
side β
3
Fig. 1. The multi-layered triangular shell element TRIC.
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
141
12 straining modes (18 Cartesian freedoms––6 rigid-body freedoms): three axial straining modes, three symmetric bending modes, three anti-symmetric bending + shear modes and three in-plane rotations. In order to generate these pure deformational modes, which are called natural modes, a projection of the nodal displacements and rotations as well as the corresponding forces and moments on the edges of the triangle takes place. The 12 · 12 stiffness matrix kN corresponding to the natural modes is denoted as the natural stiffness matrix of the element. 2.1. Random variation of Young’s modulus and Poisson’s ratio 2.1.1. Weighted integral method As the elements of the stiffness matrix are nonlinear functions of PoissonÕs ratio m, the analysis of a shell structure involving both stochastic YoungÕs modulus and PoissonÕs ratio requires the introduction of additional approximations. However, an alternative weighted integral formulation of the stochastic stiffness matrix is possible, considering LameÕs constants k and l as the two uncertain material properties since all elements of the stiffness matrix are linear functions of k and l. For the TRIC shell element and in the case of isotropic material, the elasticity matrix j12 in the material coordinate system can be expressed as 2 3 1 m 0 E 4 j12 ¼ m 1 0 5: ð1Þ 1 m2 0 0 1m The LameÕs constants are defined as mE E : ; l¼G¼ k¼ 1 m2 2ð1 þ mÞ After some algebra, Eq. (1) becomes 2 3 2 3 2 3 k þ 2l k 0 k k 0 2l 0 0 ð1Þ ð2Þ j12 ¼ 4 k k þ 2l 0 5 ¼ 4 k k 0 5 þ 4 0 2l 0 5 ¼ j12 þ j12 : 0 0 2l 0 0 0 0 0 2l |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} j12 ðkÞ
ð2Þ
ð3Þ
j12 ðlÞ
The two material properties are assumed to vary randomly along the element surface X according to kðx; yÞ ¼ k0 ½1 þ fk ðx; yÞ;
ð4aÞ
lðx; yÞ ¼ l0 ½1 þ fl ðx; yÞ; ð4bÞ where k0 , l0 are the mean values of LameÕs constants and fk ðx; yÞ, fl ðx; yÞ are 2D-1V zero-mean homogeneous stochastic fields corresponding to the spatial variation of k and l respectively. Substituting Eqs. (4a) and (4b) into the expression of the elasticity matrix j12 we have ð1Þ
ð2Þ
j12 ¼ ½j12 0 ½1 þ fk ðx; yÞ þ ½j12 0 ½1 þ fl ðx; yÞ
ð5Þ
and the constitutive matrix jct takes now the following form ð1Þ
ð2Þ
jct ¼ jct0 þ jct0 fk ðx; yÞ þ jct0 fl ðx; yÞ or
2
jaa jct ¼ 4 jab jac
jab jbb jbc
3 jac jbc 5; jcc
where subscript zero denotes the corresponding mean value.
ð6Þ
142
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
In the following, the TRIC element stiffness terms are derived from their definitions and Eq. (6). They are finally grouped in a 12 · 12 matrix which constitutes the stochastic natural stiffness matrix of the TRIC element 2 3 0 0 kqc 6 ð66Þ 7 6 0 0 7 kqh 7: kN ¼ 6 ð7Þ 6 7 ð33Þ ð1212Þ 4 5 0 0 kaz ð33Þ
For an in-depth study of the specific features of the TRIC shell element and for any details regarding the symbols used, the reader can refer to the Appendix A of this work and to publications [1,4]. Axial and symmetric bending stiffness terms Definition: kqc ¼
Z
BTtc jct Btc dV :
ð8Þ
V
ð66Þ
Stochastic form: ð1Þ
ð2Þ
kqc ¼ kqc0 þ ak kqc0 þ al kqc0 ;
ð9Þ
ð66Þ
1 ak ¼ X
Z
fk ðx; yÞ dX;
X
1 al ¼ X
Z
fl ðx; yÞ dX:
ð10Þ
X
Anti-symmetric bending and shear stiffness terms Definition (bending): Z b kqh ¼ BTth jct Bth dV :
ð11Þ
V
ð33Þ
Stochastic form: b b ¼ kqh0 þ kqh ð33Þ
6 X
Xi DkXi þ
i¼1
6 X
ð12Þ
Yi DkYi :
i¼1
Xi , Yi are the weighted integrals corresponding to stochastic fields fk ðx; yÞ and fl ðx; yÞ respectively. Xi have the following expressions Z Z X1 ¼ f2a fk ðfa ; fb ; fc Þ dX; X2 ¼ f2b fk ðfa ; fb ; fc Þ dX; X X Z Z X3 ¼ f2c fk ðfa ; fb ; fc Þ dX; X4 ¼ fa fb fk ðfa ; fb ; fc Þ dX; ð13Þ XZ ZX X5 ¼ fb fc fk ðfa ; fb ; fc Þ dX; X6 ¼ fa fc fk ðfa ; fb ; fc Þ dX; X
X
while Yi are identical to Xi but they contain the stochastic field fl instead of fk . fa , fb , fc are the area coordinates of the triangular shell element [3]. Definition (shear): Z s BTsh Xs Bsh dV : kqh ¼ ð33Þ
V
ð14Þ
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
143
Stochastic form: s kqh ð33Þ
s ¼ ð1 þ al Þkqh0 :
ð15Þ
Finally these terms are coupled in the following way b 1 s 1 1 kqh ¼ f½kqh þ ½kqh g :
ð16Þ
ð33Þ
It is worth noting that the stochastic anti-symmetric shear terms are influenced by the average value al only. This is in accordance with the nature of these terms, since al describes the variation of the shear modulus within a finite element. Stiffness terms due to in-plane rotations (azimuth stiffness) Definition:
2
( kz ¼ X max
3 0:5 0:5 5; 1
0:5 1 0:5
1 kaz ¼ kz 4 0:5 ð33Þ 0:5 1 l2a
Z
h=2
1 z jaa dz; 2 lb h=2 2
ð17Þ Z
h=2
1 z jbb dz; 2 lc h=2 2
Z
)
h=2 2
z jcc dz :
ð18Þ
h=2
Stochastic form: ð1Þ
ð2Þ
kz ¼ kz0 þ ak kz0 þ al kz0 :
ð19Þ
In (18), h is the shell thickness. 2.1.2. Local average method As it can be seen from the previous section, the axial and symmetric bending, the anti-symmetric shear and the azimuth stiffness terms have an a priori (inherent) local average form. The expression of the antisymmetric bending stiffness terms in the context of the local average method is the following bð1Þ
bð2Þ
b b kqh ¼ kqh0 þ ak kqh0 þ al kqh0 :
ð20Þ
ð33Þ
From Eq. (12), it can be seen that the number of random variables per finite element in the weighted integral formulation of the anti-symmetric bending stiffness terms is twelve compared to the two random variables per element derived in the local average approach. This gives an indication of the computational overhead associated with the weighted integral formulation. 2.2. Random variation of thickness The consideration of random variations in shell thickness complicates significantly the application of the weighted integral method [7]. Since the thickness appears into the integrals expressing the stiffness matrix of any type of shell element, the application of the weighted integral approach for the derivation of the stochastic stiffness matrix is not straightforward. In this case, the use of the local average method is recommended. Moreover, in a recent study by Argyris et al. [3], it was shown that for cylindrical shells and random variation of YoungÕs modulus, the local average method is superior to the weighted integral method in terms of simplicity and computational efficiency, while the results obtained with the use of the local average were very close to those of the weighted integral method.
144
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
Using local averages, the randomly varying shell thickness within a finite element is expressed as hðx; yÞ ¼ h0 ð1 þ ah Þ; where 1 ah ¼ X
Z
ð21Þ
fh ðx; yÞ dX
ð22Þ
X
and fh ðx; yÞ is the stochastic field describing the variation of the shell thickness. Taking into account the definitions of the various stiffness terms as stated in Eqs. (8), (11), (14) and (17), the following stochastic terms are derived • Axial and symmetric bending stiffness terms kij ¼ ðkij Þ0 ð1 þ ah Þ;
ð23aÞ
zkij ðzkij Þ0 ð1 þ 2ah Þ;
ð23bÞ
z2 kij ðz2 kij Þ0 ð1 þ 3ah Þ: • Anti-symmetric bending and shear stiffness terms
ð23cÞ
b b kqh kqh0 ð1 þ 3ah Þ;
ð24Þ
ð33Þ s s kqh ¼ kqh0 ð1 þ ah Þ:
ð25Þ
ð33Þ
• Stiffness terms due to in-plane rotations (azimuth stiffness) kaz ð33Þ
kaz0 ð1 þ 3ah Þ:
ð26Þ
For the sake of clarity in our presentation, we indicate in the following lines how Eqs. (23a) and (23b) were derived. All other expressions may be produced in a similar way. Z h=2 N N X ð21Þ X k 0 0 kij ¼ jij dz ¼ jkij ðzk zkþ1 Þ ¼ jij ðzk z0kþ1 Þ ð1 þ ah Þ ¼ ðkij Þ0 ð1 þ ah Þ; ð23a Þ h=2
k¼1
k¼1
|fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} ðkij Þ0
zkij ¼
Z
h=2
h=2
zjij dz ¼
N N 1X ð21Þ 1 X k jkij ðz2k z2kþ1 Þ ¼ j ½ðz0 Þ2 ðz0kþ1 Þ2 ð1 þ ah Þ2 2 k¼1 2 k¼1 ij k |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ðzkij Þ0
¼ ðzkij Þ0 ð1 þ 2ah þ
a2h Þ
ðzkij Þ0 ð1 þ 2ah Þ:
0
ð23b Þ
In the above equations, subscript or superscript k indicates the layer number and subscripts i, j are equal to a, b, c indicating the edges of the triangle. It can be seen that some stiffness terms are approximate since their higher order stochastic parts have been neglected. 2.3. Random variation of Young’s modulus and thickness In this section, it is considered that YoungÕs modulus and the thickness of the shell are varying simultaneously within each finite element according to
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
Eðx; yÞ ¼ E0 ð1 þ aE Þ; hðx; yÞ ¼ h0 ð1 þ ah Þ; with
Z 1 aE ¼ fE ðx; yÞ dX; X X Z 1 fh ðx; yÞ dX; ah ¼ X X
145
ð27Þ
ð28Þ
being the local averages describing the random fluctuation of the two parameters. Following a similar procedure to that of the previous section the stochastic stiffness terms have finally the following forms: • Axial and symmetric bending stiffness terms kij ¼ ðkij Þ0 ð1 þ aE þ ah þ aE ah Þ;
ð29aÞ
zkij ðzkij Þ0 ð1 þ aE þ 2ah þ 2aE ah Þ;
ð29bÞ
z2 kij ðz2 kij Þ0 ð1 þ aE þ 3ah þ 3aE ah Þ: • Anti-symmetric bending and shear stiffness terms
ð29cÞ
b kqh ð33Þ
b kqh0 ð1 þ aE þ 3ah þ 3aE ah Þ;
ð30Þ
s kqh ð33Þ
s ¼ kqh0 ð1 þ aE þ ah þ aE ah Þ:
ð31Þ
• Stiffness terms due to in-plane rotations (azimuth stiffness) kaz ð33Þ
kaz0 ð1 þ aE þ 3ah þ 3aE ah Þ:
ð32Þ
2.4. Discussion of the statistical independence assumption between the stochastic fields In the preceding sections, the cross-correlation which could eventually exist between the random parameters considered, has been neglected. The same assumption has been made in the past by several researchers who have considered randomness in more than one material or geometric property [6,12,14,19,22]. In fact, although it is likely that a certain degree of cross-correlation might exist, especially between the YoungÕs modulus and the PoissonÕs ratio of a structure (due to the similarity of the physical mechanisms that lead to random variations of material properties), there is no available evidence regarding this issue in the literature. It is also reasonable to assume that these two material properties are not correlated with a geometric property, such as the shell thickness. The assumption made in the present work is in accordance with the results presented in a recent paper by Graham and Deodatis [11]. From the numerical tests examined in that paper, it becomes evident that the cross-correlation effect does not affect significantly the response variability when static problems are considered. This remark is not valid when a (dynamic) random eigenvalue problem is treated. In this case, the results show that the eigenvalue variability and its upper bounds are highly influenced by the degree of cross-correlation assumed between YoungÕs modulus and mass density. The substantial errors which can be introduced in random vibration problems by neglecting the cross-correlations, have been also investigated by Elishakoff in [8]. It is therefore evident that the consideration of the statistical dependence existing
146
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
between structural properties becomes crucial for the case of dynamic problems. Since in the present paper, only static shell problems are examined, the statistical independence assumption will not practically affect the accuracy of the results.
3. Computation of stress variability The main scope of this work is to calculate the response statistics of shell structures with multiple random material and geometric properties subjected to static loads. In the framework of the Monte Carlo simulation (MCS) method used here, the computation of both displacement and stress variability is straightforward even for large finite element models. A few comments regarding the calculation of the TRIC shell element stresses are following. In the natural mode method, the Cartesian strains have been replaced by the total natural strains 2 3 cta ct ¼ 4 ctb 5: ð33Þ ctc These strains are measured directly parallel to the triangleÕs edges. To the total natural axial stains ct correspond component natural stresses rc grouped in the vector 2 3 rca ð34Þ rc ¼ 4 rcb 5: rcc The natural direct stresses are for each layer ðkÞ computed from the constitutive relation rkc ¼ jkct ckt
ð35Þ
and transformed to local Cartesian stresses via r0k ¼ Brkc ;
ð36Þ
pffiffiffi T where r0 ¼ rx ry 2rxy is the local Cartesian stress vector as it is defined for the TRIC shell element. Matrix B depends only on the geometry of each element and is given in the Appendix A. In the section of numerical tests, the variability of the components of the local Cartesian stress vector (36) is examined. In addition, the variability of the equivalent Von Mises stress is calculated and compared to that of each component separately. The expression of this equivalent stress for the case of plane stress problems is used in this work. It is defined as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi seq ¼ r2x þ r2y rx ry þ 3r2xy : ð37Þ For the analysis of composite structures the normal stress rz is also of interest since it may contribute to damage initiation. However, the material of the structures considered as test examples is supposed to be isotropic. This is the main reason for which only the variability of components rx , ry and rxy is finally examined. From Eqs. (35) and (36), it is also deduced that the variability of the stress components is based on the stochastic variation of both jct and ct , since stresses are calculated by a post-processing procedure. However, for the sake of simplicity and computational efficiency, only the contribution of ct in the random variation of the stresses of Eq. (35) is taken into account in the present work.
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
147
4. Representation of the stochastic field The generation of sample functions of the Gaussian zero-mean homogeneous stochastic fields which describe the random parameters of the structure, is done using the series of cosines formula of the spectral representation method. This method is well suited in the context of Monte Carlo simulation (MCS) technique used for calculating the response variability of stochastic structural systems [18]. The fast Fourier transform version of the method, used in many cases in order to reduce the computational burden of the simulation, is not exploited here because values of the field at non-uniformly spaced points are needed. For a 2D-1V stochastic field and for a specific simulation ðiÞ, we have [16] f ðiÞ ðx1 ; x2 Þ ¼
1 1 N 2 1 X pffiffiffi NX ð1ÞðiÞ ð2ÞðiÞ ð2Þ 2 ½Að1Þ n1 n2 cosðj1n1 x1 þ j2n2 x2 þ /n1 n2 Þ þ An1 n2 cosðj1n1 x1 j2n2 x2 þ /n1 n2 Þ;
ð38Þ
n1 ¼0 n2 ¼0
where /ðjÞðiÞ n1 n2 , j ¼ 1; 2 represent the realization for the ðiÞ simulation of the independent random phase angles uniformly distributed in the range ½0; 2p. Anð1Þ , Anð2Þ have the following expressions 1 n2 1 n2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð39aÞ ¼ 2Sff ðj1n1 ; j2n2 ÞDj1 Dj2 ; Anð1Þ 1 n2 Anð2Þ ¼ 1 n2
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2Sff ðj1n1 ; j2n2 ÞDj1 Dj2 ;
ð39bÞ
where j1n1 ¼ n1 Dj1 ; Dj1 ¼
j1u ; N1
j2n2 ¼ n2 Dj2 ;
ð40Þ
j2u ; N2
ð41Þ
Dj2 ¼
n1 ¼ 0; 1; . . . ; N1 1;
n2 ¼ 0; 1; . . . ; N2 1;
ð42Þ
Nj , j ¼ 1, 2, represent the number of intervals in which the wave number axes are subdivided and jju , j ¼ 1, 2, are the upper cut-off wave numbers which define the active region of the power spectrum Sff of the stochastic field. The last means that the power spectral density function Sff ðj1 ; j2 Þ is assumed to be zero outside the region defined by j1u 6 j1 6 j1u
and
j2u 6 j2 6 j2u :
ð43Þ
The two-sided power spectral density functions Sff used in the numerical tests section of this study, are given by " 2 2 # b1 j1 b2 j2 2 b1 b2 Sff ðj1 ; j2 Þ ¼ rf exp ; ð44aÞ 4p 2 2 r2f b1 b2 ðb21 j21 þ b22 j22 Þ exp½ðb21 j21 þ b22 j22 Þ; ð44bÞ p where rf denotes the standard deviation of the stochastic field and b1 , b2 denote the parameters that influence the shape of the spectrum and hence the scale of correlation [17]. It is reminded here that the power spectrum reflects the correlation structure of a random field in the wave number domain. Through random number generation and using Eq. (38), a large number of sample functions of the stochastic field is produced. These sample functions are generated in accordance with the pre-assigned second-order statistics of the simulated stochastic field (mean and variance). Sff ðj1 ; j2 Þ ¼
148
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
Before proceeding to the numerical tests section, a brief comment is made here regarding the choice of the distribution function (Gaussian) and of the power spectrum of the stochastic fields used in this paper. Strictly speaking, the choice of both functions is substantiated by available experimental data. If such data is not available, as in the case of the two shell structures considered as test examples, it is reasonable to resort to the aforementioned assumptions in order to perform a parametric investigation. Methods for the generation of sample functions of Gaussian stochastic fields are well established in the literature (e.g. spectral representation used in this paragraph) in contrast to the non-Gaussian case. However, with the choice of a Gaussian model for variables bounded by physical constraints (e.g. elastic constants), there is a non-zero probability that a violation of these constraints might occur in an actual MC simulation. The violation of the physical constraints is due to large negative realizations of the stochastic fields leading subsequently to possible non-positive definite stiffness matrices. This issue was handled in this paper by discarding the samples corresponding to a non-positive definite stiffness matrix without introducing a severe bias in the simulation. This is the case when the input stochastic fields are characterized by small coefficients of variation ðrf < 0:25Þ and thus, a small number of large negative outcomes is produced. For the two specific expressions of the power spectrum (and the underlying autocorrelation function), the square exponential type is adopted. This follows from the observation by Li and Der Kiureghian [13] who have shown that the error in the discretization (e.g. by local average method) of a homogeneous Gaussian stochastic field is negligible even for small correlation lengths when a square exponential type of correlation is adopted.
5. Numerical tests Two shell structures of different shape, denoted as Tests 1 and 2, are used in this section in order to demonstrate the proposed methodologies. Test 1 is the Scordelis-Lo shell roof, a cylindric shell subjected to gravity loading. The geometrical and material properties of the structure are presented in Fig. 2. It is assumed that the two longitudinal edges are free and the two circular edges are supported by diaphragms. The structure is discretized with 2 · 28 · 28 ¼ 1568 shell elements resulting in 4872 degrees of freedom (d.o.f.). Test 2 is a simply supported hyperboloid shell of revolution. The meridional line of the shell is described by the following hyperbola equation R1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi RðyÞ ¼ ð45Þ c2 þ y 2 ; c where c is a real constant and R1 the radius at the top of the shell. The structure is loaded by a single concentrated load of 10 kN acting in the x-direction at l ¼ 2=3L, with L ¼ 20 m the total height of the shell. The top and bottom radii are 4.8 and 8.0 m, respectively. Geometry and material data are shown in Fig. 3. A finite element mesh consisting of 1024 triangular elements and resulting in 3306 d.o.f. is used for the discretization of the shell. The spatial variability in YoungÕs modulus and thickness of the shells is described by two uncorrelated 2D-1V zero-mean homogeneous stochastic fields with coefficient of variation equal to 0.1, while the applied loads are considered to be deterministic. The same assumption is made for the stochastic fields describing the random variation of YoungÕs modulus and PoissonÕs ratio. The power spectral density function Sff characterizing the two stochastic fields in both cases is assumed to have the form of Eq. (44a) or (44b). It should be mentioned that the results appearing in Figs. 4, 5, 7 and 8 are based on the power spectrum of Eq. (44a), while the power spectrum of Eq. (44b) is used only for comparison purposes in Fig. 6. In both examples, the coefficient of variation (COV) of displacement wc at the characteristic node C of the structure, is calculated using the local average and, when feasible, the weighted integral method in conjunction with MCS. The COV of a response quantity, e.g. displacement ui , is defined as
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
149
h rigid diaphragm
q D rigid diaphragm
C A
L B
z y R x
θ
R = 7600 mm L = 15200 mm h = 76 mm θ = 40˚ 2 E = 21000 N/mm ν = 0.3
Fig. 2. Material and geometric data of the Scordelis-Lo shell roof.
Fig. 3. Material and geometric data of the hyperboloid shell.
COVðui Þ ¼
rðui Þ ; Eðui Þ
ð46Þ
where r and E are the standard deviation and mean value of ui , respectively. The response statistics have been computed using a sample size of 500. This number of samples can represent sufficiently well the preassigned first two moments of the stochastic field and statistical convergence of the results is achieved within this number of simulations. The sensitivity of structural response with regard to the scale of correlation of the stochastic fields, quantified via the correlation length parameter b is first examined. For this purpose, several sets of sample functions are generated using Eq. (38) each one for a different value of parameter b. These values are related to the dimensions of the structures along the x- and y-axes. For the Scordelis-Lo shell one value b ¼ b1 ¼ b2 is considered to be representative of both dimensions. For the hyperboloid shell, where the difference
150
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160 0.3 E h E,h
0.2
COV of w at node C
COV of w at node C
0.25
0.15 0.1
0.25 E
0.2
h E,h
0.15 0.1
0.05
0.05 0
0 0.065
0.65
(a)
1.3
2.6
26
65
0.2
130
2
(b)
Correlation length parameter b
4
8
80
200
Correlation length parameter b1
Fig. 4. (a) Scordelis-Lo shell: COV of vertical deflection wc as a function of correlation length parameter b (local average method). (b) Hyperboloid shell: COV of displacement wc as a function of correlation length parameter b1 (local average method).
0.1
COV of w at node C
COV of w at node C
0.15
0.1
0.05
0.05
0 0.065
0.65
1.3
2.6
26
65
0
130
0.2
Correlation length parameter b
(a)
E
E,v (Local average)
E,v (Weighted integral)
2
4
8
80
200
Correlation length parameter b1
(b)
E
E (Weighted integral)
E,v (Local average)
E,v (Weighted integral)
Fig. 5. (a) Scordelis-Lo shell: comparison between local average and weighted integral methods––COV of vertical deflection wc as a function of correlation length parameter b for the case of combined YoungÕs modulus and PoissonÕs ratio variation. (b) Hyperboloid shell: comparison between local average and weighted integral methods––COV of displacement wc as a function of correlation length parameter b1 for the case of combined YoungÕs modulus and PoissonÕs ratio variation.
between the two dimensions is greater, two values b1 6¼ b2 are used, each one related to the corresponding dimension. Four different cases of randomly varying parameters are considered: (i) stochastic YoungÕs modulus, (ii) stochastic thickness, (iii) combined YoungÕs modulus and thickness, and (iv) combined YoungÕs modulus and PoissonÕs ratio. The local average method is used in the first three cases for reasons explained in Section 2.2, while in the fourth case, a comparison between the local average and weighted integral methods is performed. The results concerning the displacement variability are depicted in Fig. 4a and b for the Scordelis-Lo and the hyperboloid shell cases, respectively. From these plots, it can be seen that, in both shell geometries, the displacement variability shows similar trends, starting from small values for small correlation lengths corresponding to white noise stochastic fields, up to large values for large correlation lengths. Especially, when YoungÕs modulus is the only random parameter of the problem, the COV of the selected displacement tends to the COV of YoungÕs modulus for large values of parameter b, as expected. It is also evident that randomness in shell thickness has significant effect on the response variability, compared to the effect of random variation of YoungÕs modulus [7]. When both quantities are supposed to vary randomly, the response COV tends to values that are 2–2.5 times greater than the corresponding input COV.
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160 0.25
0.2
COV of w at node C
COV of w at node C
0.25
151
h h(P.s.2)
0.15 0.1 0.05
0.2 E,h E,h(P.s.2)
0.15 0.1 0.05 0
0 0.065
(a)
0.65
1.3
2.6
26
65
0.065
130
(b)
Correlation length parameter b
0.65
1.3
2.6
26
65
130
Correlation length parameter b
0.3 0.35 0.3
COV of w at node C
COV of w at node C
0.25 0.2 0.15 h h(P.s.2)
0.1 0.05
0.2 0.15 0.1 0.05 0
0 0.2
(c)
E,h E,h(P.s.2)
0.25
2
4
8
80
Correlation length parameter b1
200
0.2
(d)
2
4
8
80
200
Correlation length parameter b1
Fig. 6. Scordelis-Lo shell: comparison between the two power spectra––COV of vertical deflection wc as a function of correlation length parameter b for the case of (a) thickness variation, (b) combined YoungÕs modulus and thickness variation. Hyperboloid shell: comparison between the two power spectra––COV of displacement wc as a function of correlation length parameter b1 for the case of (c) thickness variation, (d) combined YoungÕs modulus and thickness variation.
The very small effect of the PoissonÕs ratio stochastic variation is evidenced by Fig. 5a and b. The displacement variability for the case of randomness in both LameÕs constants is approximately the same or even slightly smaller compared to the case of YoungÕs modulus stochasticity. Similar conclusions have been reported in [22] for plane stress and plate problems. A comparison of the local average and weighted integral methods is performed next. For the ScordelisLo shell roof, these differences are negligible as it can be seen in Fig. 5a, while the discrepancy becomes more pronounced in the hyperboloid shell example (Fig. 5b). This can be explained by the fact that a relatively coarse mesh (compared to that used for the Scordelis-Lo shell) is used for the hyperboloid shell. This behavior confirms the observation that a relatively fine mesh is required by the local average method in order to produce accurate results. Furthermore, the effect of type of correlation, expressed by the mathematical form of the power spectrum, is investigated using the expressions (44a) and (44b). Two cases of random variation are considered for each structure: (i) stochastic thickness, (ii) stochastic YoungÕs modulus and thickness. The results are depicted in Fig. 6a–d. It is observed that, in both shell cases, some differences exist between the COVs calculated using the two spectra. For the hyperboloid shell example these differences are substantial in the whole range of b examined, while for the Scordelis-Lo shell the discrepancy becomes more pronounced as the correlation length parameter increases. This observation can be justified by the fact that the variance of a random response quantity (in this case the selected displacement) is mathematically related to the power spectral density function of the stochastic field describing the random input quantities (in this case the YoungÕs modulus and the thickness) [11]. It is therefore reasonable that spectra, having their spectral power concentrated over different wave number regions, produce slightly different response variations.
152
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
In Tables 1 and 2, a complete image of the response statistics of the two shells under consideration is presented, using the local average approach. The mean value, standard deviation and COV of displacement wc at the characteristic node C of the structures are calculated for the following cases: (i) in a ‘‘mixed’’ way, i.e. when the two random parameters are introduced simultaneously in the model and (ii) approximately, using the SRSS concept. SRSS denotes the square root of the sum of squares of the two outputs by two random parameters, E and h, considered separately. The assumption of independence between the stochastic fields describing the two properties is well reflected by the first four entries of the sixth column of the tables, where the ratio of mean values is approximately equal to the square root of 2. It is also evident that the SRSS formula can offer an accurate prediction of the results produced by the combined variation case [6]. Finally, it can be seen once again, that randomness in shell thickness has a significant effect on the displacement variability, compared to the effect of random variation in YoungÕs modulus. Up to this moment, the displacement variability of the two shell structures has been considered. In the context of this study, similar investigations have been performed with regard to the stress variability. The stress components examined are those described in Section 3. The influence of the correlation length parameter b on the stress variability is first tested and the corresponding plots of stresses inside some specific elements are depicted in Figs. 7 and 8. By comparing these plots with those of the displacement variability, different trends of variation can be observed among several stress components. It can be also seen that the ‘‘distance’’ separating the results produced by random thickness variation and by combined YoungÕs modulus and thickness variation, is greater for the stress variability case than for the case of displacement variability. In order to further examine the variability of stress components within several finite elements, Tables 3 and 4 are presented. The results contained in these tables correspond to a specific value of correlation length parameter b for each structure (b ¼ 65:0 for the Scordelis-Lo shell and b1 ¼ 200:0 for the hyperboloid shell). It is clear that the fluctuation is substantially different not only between several stress components but even between the same stress component within different elements. This variation is more pronounced for the case of random thickness fluctuation which indicates how important may be the calculation of probability- and spectral-distribution-free upper bounds of the response (displacement and stress) variability.
Table 1 Scordelis-Lo shell: displacement variability as a function of correlation length parameter b (local average method and power spectrum of Eq. (44a)) Items
Correlation length parameter b
Random parameters YoungÕs modulus E
Thickness h
E and h (mixed)
SRSS/ mixed
Mean value
0.065 0.65 2.6 65
41.8757 41.9306 42.0061 41.8467
42.1692 42.3711 42.5307 43.212
42.2713 42.5249 42.7531 43.2767
1.4059 1.4018 1.3982 1.39
Standard deviation
0.065 0.65 2.6 65
0.2046 0.4701 1.6144 3.9851
0.2168 0.8473 3.5007 8.4818
0.2904 0.9834 3.8937 9.3016
1.0265 0.9853 0.9901 1.0075
COV
0.065 0.65 2.6 65
0.00488 0.01121 0.03843 0.09523
0.00514 0.01999 0.08231 0.19628
0.00687 0.02312 0.09107 0.21493
1.0317 0.9913 0.9975 1.015
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
153
Table 2 Hyperboloid shell: displacement variability as a function of correlation length parameter b1 (local average method and power spectrum of Eq. (44a)) Items
Correlation length parameter b1
Random parameters YoungÕs modulus E
Thickness h
E and h (mixed)
Mean value
0.2 2 8 80
4.786 4.7949 4.8214 4.7804
4.8337 4.8798 4.9828 5.0057
4.8441 4.8974 5.0305 5.0089
1.4042 1.3969 1.3783 1.3819
Standard deviation
0.2 2 8 80
0.0546 0.1443 0.2833 0.4486
0.1475 0.4036 0.8697 1.2018
0.1565 0.4189 0.9357 1.2655
1.005 1.0232 0.9775 1.0137
COV
0.2 2 8 80
0.01141 0.03009 0.05877 0.09384
0.03052 0.0827 0.17453 0.24009
0.03232 0.08554 0.18601 0.25266
1.0081 1.0288 0.9901 1.0203
0.25
SRSS/mixed
0.25 σy(E)
σx(h)
0.2
COV of σy at point A
COV of σx at point A
σx(E) σx(E,h)
0.15 0.1 0.05
σy(h)
0.2
σy(E,h)
0.15 0.1 0.05 0
0 0.065
0.65
(a)
1.3
2.6
26
65
0.065
130
(b)
Correlation length parameter b
2.6
26
65
130
65
130
τeq(E)
σxy(h)
0.4
σxy(E,h)
0.3 0.2 0.1
COV of τeq at point A
σxy(E)
COV of σxy at point A
1.3
Correlation length parameter b
0.25
0.5
τeq(h)
0.2
τeq(E,h)
0.15 0.1 0.05
0
0 0.065
(c)
0.65
0.65
1.3
2.6
26
Correlation length parameter b
65
130
0.065
(d)
0.65
1.3
2.6
26
Correlation length parameter b
Fig. 7. (a) Scordelis-Lo shell: COV of stress rx at point A as a function of correlation length parameter b (local average method). (b) Scordelis-Lo shell: COV of stress ry at point A as a function of correlation length parameter b (local average method). (c) Scordelis-Lo shell: COV of stress rxy at point A as a function of correlation length parameter b (local average method). (d) Scordelis-Lo shell: COV of equivalent stress seq at point A as a function of correlation length parameter b (local average method).
154
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160 0.2
σy(E)
σx(E)
0.2
COV of σy at point A
COV of σx at point A
0.25
σx(h) σx(E,h)
0.15 0.1 0.05 0
(a)
2
4
8
80
0.1
0.05
0.2
200
(b)
Correlation length parameter b1 0.2
2
4
8
80
200
Correlation length parameter b1 0.2
σxy(E)
τeq(E)
σxy(h)
0.15
COV of τeq at point A
COV of σxy at point A
σy(E,h)
0 0.2
σxy(E,h)
0.1
0.05
τeq(h)
0.15
τeq(E,h)
0.1
0.05
0
0 0.2
(c)
σy(h)
0.15
2
4
8
80
Correlation length parameter b1
0.2
200
(d)
2
4
8
80
200
Correlation length parameter b1
Fig. 8. (a) Hyperboloid shell: COV of stress rx at point A as a function of correlation length parameter b1 (local average method). (b) Hyperboloid shell: COV of stress ry at point A as a function of correlation length parameter b1 (local average method). (c) Hyperboloid shell: COV of stress rxy at point A as a function of correlation length parameter b1 (local average method). (d) Hyperboloid shell: COV of equivalent stress seq at point A as a function of correlation length parameter b1 (local average method).
6. Summary and conclusions In this work, a stochastic finite element analysis of shell structures with combined uncertain material and geometric properties is performed. For this purpose, a stochastic formulation of the triangular composite facet shell element TRIC is derived assuming random variation of the YoungÕs modulus, the PoissonÕs ratio and the thickness of the shell. As a result of the proposed formulation and the special features of the element, the stochastic stiffness matrix of TRIC depends finally on a minimum number of random variables representing the stochastic field. This fact leads to a cost-effective stiffness matrix, which is very important in the case of a time consuming stochastic analysis of real world structures. The spatial fluctuation of the mechanical and geometric properties is described by uncorrelated stochastic fields. Under the assumption of a pre-specified power spectral density function for these stochastic fields, it is possible to compute the response variability of a realistic shell structure using the direct MCS technique. Two shells of cylindrical and hyperboloid shape are tested. The sensitivity of response statistics with regard to the scale of correlation of the stochastic fields, quantified via the correlation length parameter b is first examined. For both shell geometries, the displacement variability shows similar trends, starting from small values for small correlation lengths, up to large values for large correlation lengths. When the YoungÕs modulus is the only random parameter of the problem, the COV of the selected displacement tends to the COV of YoungÕs modulus for large values of parameter b, as expected. It is also evident that random variation in the shell thickness has significant effect on the displacement variability compared to the effect of
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
155
Table 3 Scordelis-Lo shell: stress variability for b ¼ 65:0 (local average method and power spectrum of Eq. (44a)) Items
Stress component, point, layer
Random parameters YoungÕs modulus E
Thickness h
E and h (mixed)
SRSS/mixed
Mean (absolute) value
rx , B, 2 rxy , B, 2 ry , B, 1 ry , A, 2 rx , D, 1 ry , D, 2 rx , C, 1
2.0766 0.2837 15.9683 1.2376 0.4926 7.7248 1.3981
2.1221 0.2393 16.2278 1.2724 0.5081 7.9053 1.3835
2.124 0.2399 16.2437 1.2739 0.5086 7.9134 1.3864
1.3979 1.5471 1.4016 1.3934 1.3914 1.3967 1.4187
Standard deviation
rx , B, 2 rxy , B, 2 ry , B, 1 ry , A, 2 rx , D, 1 ry , D, 2 rx , C, 1
0.1952 0.0277 1.5013 0.1169 0.0465 0.7301 0.1316
0.3251 0.1585 2.1852 0.2589 0.0963 1.2837 0.1429
0.3784 0.1562 2.659 0.2843 0.1063 1.4752 0.2008
1.0021 1.0301 0.9971 0.9992 1.006 1.0011 0.9675
COV
rx , B, 2 rxy , B, 2 ry , B, 1 ry , A, 2 rx , D, 1 ry , D, 2 rx , C, 1
0.094 0.0975 0.094 0.0945 0.0945 0.0945 0.0941
0.1532 0.6623 0.1347 0.2035 0.1895 0.1624 0.1033
0.1782 0.6511 0.1637 0.2231 0.2089 0.1864 0.1448
1.0086 1.0282 1.0034 1.0057 1.0137 1.008 0.965
Table 4 Hyperboloid shell: stress variability for b1 ¼ 200:0 (local average method and power spectrum of Eq. (44a)) Items
Stress component, point, layer
Random parameters YoungÕs modulus E
Thickness h
E and h (mixed)
SRSS/mixed
Mean (absolute) value
rx , D, 2 rxy , D, 2 ry , D, 1 rx , B, 1 ry , B, 2 rxy , B, 2 ry , A, 2
0.0363 0.0272 0.0485 0.1299 0.0278 0.0261 0.2806
0.0361 0.0275 0.0513 0.1298 0.0278 0.0253 0.2857
0.0362 0.0276 0.0514 0.1309 0.0279 0.0254 0.2862
1.4142 1.4014 1.3735 1.4029 1.4091 1.4311 1.3992
Standard deviation
rx , D, 2 rxy , D, 2 ry , D, 1 rx , B, 1 ry , B, 2 rxy , B, 2 ry , A, 2
0.0034 0.0026 0.0046 0.0123 0.0026 0.0025 0.0267
0.0032 0.0037 0.0142 0.0116 0.0031 0.0027 0.0403
0.0048 0.0045 0.0148 0.0166 0.0042 0.0037 0.0485
0.9727 1.0049 1.0085 1.0185 0.9633 0.9945 0.9967
COV
rx , D, 2 rxy , D, 2 ry , D, 1 rx , B, 1 ry , B, 2 rxy , B, 2 ry , A, 2
0.0937 0.0955 0.0948 0.0947 0.0935 0.0958 0.0952
0.0893 0.1345 0.2768 0.0894 0.1115 0.1067 0.1411
0.1326 0.1631 0.2879 0.1268 0.1505 0.1457 0.1695
0.9762 1.0114 1.0163 1.0271 0.9669 0.9842 1.0042
156
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
random YoungÕs modulus. When both quantities are supposed to vary randomly, the response COV tends to values that are 2–2.5 times greater than the input COV. In the case of combined fluctuation of YoungÕs modulus and PoissonÕs ratio, the displacement variability curve indicates the small effect of the PoissonÕs ratio on the results. Slight differences are found between the results given by the local average method and those of the weighted integral method in the case of Scordelis-Lo shell. These differences become somewhat more pronounced in the hyperboloid shell case. In addition, the effect of type of correlation, expressed in this work by the mathematical form of the power spectrum, is investigated. Due to the relation existing between the power spectral density function of the stochastic field describing the random input quantities and the variance of a response quantity, slight differences are once again observed. Finally, investigations similar to the aforementioned were made with regard to the stress variability. The main observation in this case is that the stress fluctuation is substantially different not only between several stress components but even between the same stress component within different elements. This variation takes its maximum value when the YoungÕs modulus and the shell thickness are simultaneously varying. Acknowledgement This work was partially supported by the research project ‘‘Thalis’’ of the National Technical University of Athens. Appendix A. Stochastic stiffness matrix components In Eq. (8) of the present work, where the definition for the axial and symmetric bending stiffness terms of the TRIC shell element is given, Btc is a matrix of the form 3 2 z 0 0 1 0 0 7 6 la 7 6 60 1 0 0 z 0 7 ðA:1Þ Btc ¼ 6 7; lb 7 6 ð36Þ 4 z5 0 0 1 0 0 lc where la , lb , lc are the lengths of the triangleÕs edges. In Eq. (9), the mean value of matrix kqc is defined as 2 zkaa zkab zkac 3 kab kac kaa 6 la lb lc 7 7 6 7 6 zk zk zk ab bb bc 7 6 kbb kbc 6 la lb lc 7 7 6 7 6 zkac zkbc zkcc 7 6 7 6 kcc 6 la lb lc 7 7 6 ; kqc0 ¼ X 6 z2 kaa z2 kab z2 kac 7 7 6 symm: ð66Þ 7 6 l2a la lb la lc 7 6 7 6 2 6 z kbb z2 kbc 7 7 6 6 lb lc 7 l2b 7 6 7 6 4 z2 kcc 5 l2c
ðA:2Þ
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
157
where kij ¼
Z
h=2
jij dz ¼
h=2
zkij ¼
Z
h=2
zjij dz ¼
Z
jkij ðzk zkþ1 Þ;
i; j ¼ a; b; c;
ðA:3Þ
k¼1
h=2
z2 kij ¼
N X
N 1X jk ðz2 z2kþ1 Þ; 2 k¼1 ij k
h=2
z2 jij dz ¼
h=2
N 1X jk ðz3 z3kþ1 Þ; 3 k¼1 ij k
i; j ¼ a; b; c;
i; j ¼ a; b; c;
ðA:4Þ
ðA:5Þ
X is the area of the triangle, h is the element thickness, N the number of layers and zk the coordinate of the layer k as shown in Fig. 9. ð1Þ ð2Þ Also in Eq. (9), matrices kqc0 , kqc0 are identical in form with matrix kqc0 but they depend only on k and l, respectively. b In the definition of the anti-symmetric bending stiffness terms kqh (Eq. (11)) Bth ¼ zl2 Ah l; l is a diagonal 2 la l¼40 0
ðA:6Þ
matrix, which contains the lengths of the triangleÕs edges 3 0 0 lb 0 5 0 lc
and matrix Ah is a function of the triangular or natural coordinates fa , fb , fc [3] 2 3 fa fa 3ðfc fb Þ 5: 3ðfa fc Þ fb Ah ¼ 4 fb fc fc 3ðfb fa Þ b The mean value of matrix kqh is defined as 2 b11 3 b12 b13 kqh kqh kqh 6 7 b b22 b23 7 kqh kqh ¼6 kqh0 4 5: b33 symm: kqh
ðA:7Þ
ðA:8Þ
ðA:9Þ
b The entries of matrix kqh0 depend on the geometric characteristics of the TRIC element (area and lengths of the triangleÕs edges) and on the stiffness quantities defined in (A.5) [3].
Fig. 9. Coordinates z of the layers.
158
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
All matrices Dk belonging to the fluctuating part of the stiffness terms given by Eq. (12) are deterministic and have the following forms 2 DkX1
6 6 6 ¼6 6 4
0
l2b 2 ð1Þ zk l4a aa
0 þ
3
0
6 2 ð1Þ 9 2 ð1Þ zk þ zk l2a ab l2b bb
symm:
lb lc 2 ð1Þ 3lc 2 ð1Þ 3lb 2 ð1Þ 9 2 ð1Þ 7 z kaa 2 z kab 2 z kac zk 7 4 lb lc bc 7 la la l b la lc 7; 7 l2c 2 ð1Þ 6 2 ð1Þ 9 2 ð1Þ 5 z k þ z k þ z k l2c cc l4a aa l2a ac ðA:10Þ
2
DkX2
9 2 ð1Þ 6 2 ð1Þ l2a 2 ð1Þ 6 l2 z kaa þ l2 z kab þ l4 z kbb b b 6 a ¼6 6 4 symm: 2
DkX3
9 2 ð1Þ 6 2 ð1Þ l2a 2 ð1Þ 6 l2 z kaa þ l2 z kac þ l4 z kcc 6 a c c 6 ¼6 6 4 symm: 2
DkX4
60 6 6 6 ¼6 6 6 4 symm:
10 2 ð1Þ zk la lb ab 0
0 0
3 3lc 2 ð1Þ la lc 2 ð1Þ 9 2 ð1Þ 3la 2 ð1Þ z k z k z k z k bb la lc ac la l2b ab l4b l2b lc bc 7 7 7; 0 7 2 5 lc 2 ð1Þ 6 2 ð1Þ 9 2 ð1Þ z k þ z k þ z k bb bc cc 4 2 2 lc lb lb
9 2 ð1Þ 3la 2 ð1Þ 3lb 2 ð1Þ la lb 2 ð1Þ zk zk z k 4 z kcc la lb ab lb l2c bc la l2c ac lc 9 2 ð1Þ 6 2 ð1Þ l2b 2 ð1Þ z k þ z k þ 4 z kcc lc l2b bb l2c bc
2 DkX6
6 6 6 ¼6 6 4
0
symm:
07 7 7 7; 07 5
ðA:12Þ
0
2
DkX5
3
3
lc 2 ð1Þ 9 2 ð1Þ 3la 2 ð1Þ zk þ z k þ 2 z kbc 7 la l2b ab la lc ac lb lc 7 7 lc 2 ð1Þ 3lb 2 ð1Þ 9 2 ð1Þ 7 z kab þ 2 z kac þ z kbc 7; 2 lb lc la lb l a lc 7 2l2c 2 ð1Þ 6 2 ð1Þ 6 2 ð1Þ 18 2 ð1Þ 7 5 2 2 z kab 2 z kac 2 z kbc 2 z kcc la lc lb la lb
18 2 ð1Þ 6 2 ð1Þ 6 2 ð1Þ 2l2a 2 ð1Þ 6 l2 z kaa l2 z kab l2 z kac l2 l2 z kbc 6 a c b b c 6 9 2 ð1Þ 3lb 2 ð1Þ la 2 ð1Þ 6 zk þ z k þ 2 z kbc ¼6 la lb ab la l2c ac lb lc 6 6 4 3lc 2 ð1Þ la 2 ð1Þ 9 2 ð1Þ z kab þ 2 z kbc þ zk 2 la lc ac la lb lb lc
ðA:11Þ
0
10 2 ð1Þ zk lb lc bc
lb 2 ð1Þ 9 2 ð1Þ 3la 2 ð1Þ zk þ zk þ zk la lb ab lb l2c bc la l2c ac 2l2b 6 18 ð1Þ 6 ð1Þ ð1Þ ð1Þ 2 z2 kab 2 2 z2 kac 2 z2 kbb 2 z2 kbc la lc la lc lb
3 symm: 7 7 7 7 7; 7 7 5 0
ðA:13Þ
ðA:14Þ
3 10 2 ð1Þ z kac 7 la lc 7 7 3lc 2 ð1Þ lb 2 ð1Þ 9 2 ð1Þ 7: 7 z k þ z k þ z k lb lc bc 5 l2a lb ab l2a lc ac 0
ðA:15Þ Matrices DkYi are identical to DkXi , with superscript (2) in the place of (1). Superscript (1) means dependence on k, while superscript (2) means dependence on l.
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
159
The mean value of the stochastic anti-symmetric shear stiffness terms given by Eq. (14) is contained in a matrix of the form 2
s kqh ¼
1 4 12
K111 þ K322
K112 11 K2 þ K122
symm:
3 K312 5: K212 11 K3 þ K222
ðA:16Þ
The entries of this matrix are functions of the geometry and of the mean shear modulus of the shell [3]. Finally matrix B, which appears in Eq. (36) and relates the local Cartesian stresses to the direct natural stresses, has the following expression 2
c2ax0 6 c2 B ¼ 4 bx0 c2cx0
s2ax0 s2bx0 s2cx0
3 2cax0 sax0 2cbx0 sbx0 7 5; 2ccx0 scx0
ðA:17Þ
where ax0 , bx0 , cx0 are the angles that the triangleÕs edges a, b, c form with the local x0 axis respectively and s, c are the sines and cosines of these angles.
References [1] J. Argyris, M. Papadrakakis, C. Apostolopoulou, S. Koutsourelakis, The TRIC shell element: theoretical and numerical investigation, Comput. Methods Appl. Mech. Engrg. 182 (2000) 217–245. [2] J. Argyris, M. Papadrakakis, Z.S. Mouroutis, Nonlinear dynamic analysis of shells with the triangular element TRIC, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3005–3038. [3] J. Argyris, M. Papadrakakis, G. Stefanou, Stochastic finite element analysis of shells, Comput. Methods Appl. Mech. Engrg. 191 (2002) 4781–4804. [4] J. Argyris, L. Tenek, L. Oloffson, TRIC: a simple but sophisticated 3-node triangular element based on 6 rigid-body and 12 straining modes for fast computational simulations of arbitrary isotropic and laminated composite shells, Comput. Methods Appl. Mech. Engrg. 145 (1997) 11–85. [5] B. Bhattacharyya, S. Chakraborty, NE-MCS technique for stochastic structural response sensitivity, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5631–5645. [6] C.K. Choi, H.C. Noh, Stochastic analysis of shape imperfections in RC cooling tower shells, J. Struct. Engrg. (ASCE) 126 (2000) 417–423. [7] C.K. Choi, H.C. Noh, Stochastic finite element analysis of plate structures by weighted integral method, Struct. Engrg. Mech. 4 (1996) 703–715. [8] I. Elishakoff, On the role of cross-correlations in random vibrations of shells, J. Sound Vibr. 50 (1977) 239–252. [9] G. Falsone, N. Impollonia, A new approach for the stochastic analysis of finite element modeled structures with uncertain parameters, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5067–5085. [10] R.G. Ghanem, R.M. Kruger, Numerical solution of spectral stochastic finite element systems, Comput. Methods Appl. Mech. Engrg. 129 (1996) 289–303. [11] L.L. Graham, G. Deodatis, Response and eigenvalue analysis of stochastic finite element systems with multiple correlated material and geometric properties, Probabilist. Engrg. Mech. 16 (2001) 11–29. [12] L.L. Graham, E.F. Siragy, Stochastic finite element analysis for elastic buckling of stiffened panels, J. Engrg. Mech. (ASCE) 127 (2001) 91–97. [13] C.C. Li, A. Der Kiureghian, Optimal discretization of random fields, J. Engrg. Mech. (ASCE) 119 (1993) 1136–1154. [14] P.L. Liu, A. Der Kiureghian, Finite element reliability of geometrically nonlinear uncertain structures, J. Engrg. Mech. (ASCE) 117 (1991) 1806–1825. [15] H.G. Matthies, C.E. Brenner, C.G. Bucher, C. Guedes Soares, Uncertainties in probabilistic numerical analysis of structures and solids––stochastic finite elements, Struct. Safety 19 (1997) 283–336. [16] M. Shinozuka, G. Deodatis, Simulation of multi-dimensional Gaussian stochastic fields by spectral representation, Appl. Mech. Rev. (ASME) 49 (1996) 29–53.
160
G. Stefanou, M. Papadrakakis / Comput. Methods Appl. Mech. Engrg. 193 (2004) 139–160
[17] C.A. Schenk, G.I. Schu€eller, Buckling analysis of cylindrical shells with random geometric imperfections, Int. J. Non-Linear Mech. 38 (2003) 1119–1132. [18] P.D. Spanos, B.A. Zeldin, Monte Carlo treatment of random fields: a broad perspective, Appl. Mech. Rev. (ASME) 51 (1998) 219–237. [19] B. Van den Nieuwenhof, J.P. Coyette, Modal approaches for the stochastic finite element analysis of structures with material and geometric uncertainties, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3705–3729. [20] E. Vanmarcke, Random Fields: Analysis and Synthesis, The MIT Press, Cambridge, MA, 1983. [21] J. Zhang, B. Ellingwood, Orthogonal series expansions of random fields in reliability analysis, J. Engrg. Mech. (ASCE) 120 (1994) 2660–2677. [22] W.Q. Zhu, Y.J. Ren, W.Q. Wu, Stochastic FEM based on local averages of random vector fields, J. Engrg. Mech. (ASCE) 118 (1992) 496–511.