Consistent thick shell element

Pergamon PII:

CONSISTENT

Compulers & Srrucrures Vol. 65, No 4, pp 531-549, 1997 Q 1997 Published by Elwner Science Ltd All nghts reserved Printed I” Great Britam SOO45-7949(96)004&t-2 0045-7949197 s17.00 + 0.00

THICK

SHELL

ELEMENT

Brad L. Kozieyt and Farooque A. Mirza Department of Civil Engineering, McMaster University, Hamilton, Ontario, Canada LSS 4L7 (Received 15 November 1995)

Abstract-Thick shell finite element with transverse shear deformation have required the use of reduced integration t,o provide improved results for thin plates and shells, due to the presence of spurious transverse shear strain modes. It has been found that the spurious transverse shear strain modes result from inconsistencies in the displacement fields used in the formulation of these elements. A new thick shell element has been formulated. By providing cubic polynomials for approximation of displacements, and quadratic polynomials for approximation of rotations a consistent formulation is ensured thereby eliminating :the spurious modes. Rotational degrees of freedom which vary quadratically through the

thickness of the element are included. This allows for a parabolic variation of the transverse shear strams and hence eliminates the need for use of the shear correction factor K as required by the Mindlin plate theory. These rotational degrees of freedom also provide cubic variations of the displacements through the thickness of the element. Thus, the normal to the middle surface is neither straight nor normal after shearing and bending, allowing for warping of the cross-section. Material non-linearities are also incorporated, along with the modified Newton-Raphson method for nonlinear analysis. Comparisons are made with the available elasticity solutions and those predicted by the eight and nine-node isoparametric shell elemenrs. The results indicate that the consistent thick shell element provides excellent predictions of the displacements, stresses and plastic zones for both thin and thick plates and shells. 0 1997 Published by Elsevier !3cience Ltd.

INTRODUCTION

Numerous attempts had been made at formulating a shell element for the analysis of general shell structures. The nine-node isoparametric shell element, developed by Ahmed et al. [I], has been a successful attempt and has enjoyed much use since its introduction. It includes transverse shear deformation through the use of the Reissner [2] and Mindlin [3] theories, which assume the transverse shearing strains to be constant through the thickness of the shell. It is well known that the element performs well in the analysis of moderately thick plates and shells. However, as the shell thickness decreases the performance of the element has been found to deteriorate. This behaviour is commonly referred to as shear and/or membrane locking [4] and leads to poor predictions of both the displacements and stresse.s. Many suggestions have been made to obtain improved results when thin plates or shells are analysed. These include selective or reduced integration [S-7], heterosis formulation [8], and application of stabilization matrices [9-l 11. Other schemes, which have also been employed, but with less success, include the assumed strain method [12, 131 and the addition of nonconforming displacement modes [14]. While some of these methods perform ,well, there is as yet no method that is generally accepted because of the inherent t To whom all correspondence should be addressed.

limitations of each method. For example, the application of selective or reduced integration leads to the development of zero energy modes, which for certain boundary conditions may lead to singularity of the global stiffness matrix, or to a solution obscured by spurious kinematic modes. In general, these methods are usually complicated, difficult to implement, and cannot be applied with complete confidence to general situations, especially for analyses beyond the elastic range. Recently, several papers have been published on the formulation of plate and shell elements that include transverse shear deformation and are based on new higher order plate theories. A Co finite element formulation of the higher order theory, due to Lo et al. [15], has been presented by Kant et al. [16]. In their formulation the displacements are expressed using a Taylor series expansion in terms of the thickness co-ordinate. For the inplane displacements the series is truncated after the third degree term, while the transverse displacement series is truncated after the second degree term. Thus, the formulation incorporates a quadratic variation of the transverse shearing strains and a linear variation of the transverse normal strain. As a result the transverse normal stress is included in the formulation. The resulting element employs the nine-noded Lagrange quadrilateral with six degrees of freedom per node. Although good results have been reported,

531

532

B. L. Koziey and F. A. Mirza

the element is not free of spurious shear strain modes. It was also pointed out that further investigation was required to determine the behaviour of the element for those boundary conditions known to lead to the development of zero energy modes. Another finite element formulation based on a new higher order plate theory is that given by Reddy [17], in which the displacement field given by Levinson [ 181 and Murthy [19] is corrected and applied to laminated composite plates. In this case the inplane displacements are expanded as cubic functions of the thickness co-ordinate while the transverse deflection is independent of the thickness co-ordinate. Thus, the transverse normal stress is assumed to be negligible. The four additional higher order functions used in the definition of the inplane displacements are eliminated and expressed in terms of the lower order functions by satisfying the zero transverse shear stress conditions on the top and bottom surfaces of the plate. The resulting formulation contains second order derivatives of the transverse displacement in the strain energy expression and consequently the formulation requires the use of C’ continuous shape functions. In addition, Ren and Hinton [20] pointed out that for thin plates the transverse shear strains do not tend towards zero at the mid-plane as required. They overcame this difficulty by modifying the displacement field given by Reddy. However, like the previous formulation they also required the use of C’ continuous shape functions. Based on the preceding considerations an ‘ideal’ shell element should (i) not exhibit shear and or membrane locking; (ii) not contain any spurious zero energy modes; (iii) satisfy the usual isotropy and convergence requirements; (iv) not be based on any numerically adjusted factors; (v) be capable of predicting accurately displacements and stresses in both thin and thick plates and shells; and (vi) be easy to implement and use. Much of the effort in recent years has been directed towards the development of such an element. However, at this time there appears to be no element available which satisfies all of the above criteria. A new Co shell element, which includes transverse shear deformation, is formulated and is referred to as the consistent thick shell element. The order of the polynomials used for approximations of displacements and rotations within the new element are chosen to ensure that the element is free of the spurious transverse shear strain modes that plagued other shell elements in the past. In addition, the adopted polynomials are also complete, since a triangular element is employed instead of a rectangular one. New rotational degrees of freedom which vary quadratically through the thickness of the element are included. These new rotations allow for a parabolic representation of the transverse shear strains through the thickness of the element. This eliminates the need for use of a shear correction factor in the formulation, which is usually required

when the Mindlin [3] plate theory is used. In addition, the quadratically varying rotations produce a cubic variation of the tangential displacements and, hence, the ‘normal’ to the middle surface does not remain straight or normal after shearing and bending. The stress normal to the middle surface is assumed to be negligible. For material non-linearity an incremental elasto-plastic constitutive relation is used which employs the Huber-von Mises yield criterion, an associated flow rule and isotropic hardening. The nonlinear finite element analysis is performed using a modified Newton-Raphson method. Special attention is given to the efficiency of the numerical implementation. Numerical results for a variety of problems are presented for performance evaluation of the consistent thick shell element. These include elastic, and elasto-plastic analyses of thin and thick plates and shells, and comparisons are carried out with results predicted by other shear deformable shell elements, such as the eight and nine-node isoparametric shell elements. CONSISTENCY OF DISPLACEMENT FIELD

The inconsistencies in the formulation of the nine-node isoparametric shell element by Ahmed et al. [ 11,and other similar elements developed in the past, is demonstrated by considering the elemental strain energy u’ which employs the Mindlin plate theory [3] u’=

l-v a!, + B: + 2va,$,, + --2-

+ 4

1

(a,, + p.V)2

b-~,~- a)2 + (w,, - /I)‘] dA.

(1)

In the expression above D = flexural rigidity of the plate = Eh3/(12(1 - v)), E = modulus of elasticity, v = Poisson’s ratio, G = shear modulus, K = shear correction factor, w = transverse displacement, a = rotation due to flexure about y-axis, /3 = rotation due to flexure about x-axis, h = thickness and A = planar area of the plate. Also note that the subscripts ,I and ,, represent partial derivatives with respect to the x and y coordinates, respectively. Only the bending action is required to demonstrate the inconsistency in the formulation of the above mentioned elements. The first term, i.e. D[. . .] in eqn (1) represents strain energy due to bending only, and the second term i.e. (G~)/K[. . .] is the contribution to the elemental strain energy due to transverse shearing. The Mindlin assumption of the normal to the middle surface remaining straight, but not normal, after bending and shearing is employed. The factor K accounts for the known parabolic shear strain

Consistent thick shell element distribution across the depth of the plate. Hence, W, and w,. are the total rotations about the y- and x-axes, respectively. In the past the same degree of polynomials had been employed for approximation of all three variables, w, CI .and /?, e.g. the eight and ninenode isoparametric shell elements. To demonstrate the inconsistency in the formulation of such elements and to assert the consistency of the displacement field, consider the following quadratic approximations x(x, Y) = a~ + a2x + a3y + (14x2+ aSxy + a6y2 = i a,xm,y”, (2) ,=I

KG Y) = b, + b2X + bsy + tpax2+ bsxy + &y2 = f: b,Xmyn, (3) t-1

4-G Y) = 6

Cl +

c2x

+

c,y +

&X2

+

csxy

+

csy2 =

Substitution of these approximations energy expression in eqn (1) yields

1

,=I

c,x"'*y"'.

(4)

into the strain

(m,a,xm,- ‘~q)~+ (n,b,~“‘y”-l)~

+ 2v(m,a,.P- 'ynf)(n,a,xm,ynt1)

533

+1-v -j-(n,a,xmynf-l

+ ?

+

[(m,c,x”-lynl -

m,b,xT-

'Y%)~

a,xTy"))2

+ (n,C,Xm,y”, - ’ - b,~“‘y”,)~] dA.

r Active Degrees Of Fdom

yv v

(5)

The discrete equations obtained via the first variation of Ue with respect to a, and b,, which involve the bending and shear strain energy terms, present no problems when the coefftcients a, and b, are solved for in terms of c,. However, the discrete equations obtained via the first variation with respect to c,, which involve only the shear strain energy of the plate, are problematic. The terms w,, and w,! are linear in x and y, whereas both c( and /I are quadratic. Therefore, the total rotation terms w,, and w,~lack the matching quadratic terms in c( and B. More specifically the quadratic modes associated with coefficients a4, a5 and a6 for c(, and b4, bs and b6 for fi are not matched and, hence, lead to the so called spurious shear strain modes unless they can be determined uniquely from flexural considerations, either via the rotational equilibrium or from the boundary conditions. To the authors’ best knowledge, one such case of unique determination is that of cylindrical bending in which the coefficients are equal to zero. This discrepancy is observed regardless of the aspect ratio (typical length to thickness). The degree of inconsistency also depends upon the term (Gh)/rc in eqn (1) which acts as a Lagrange multiplier for the constraints, i.e. w,, = a and w,,=b for a thin plate. As such, for small values of (G~)/K the

2\

./

1

0

‘i

‘i

wi

3

Pi

4i

Vi

Fig. 1. Consistent thick shell element co-ordinate systems and nodal degrees of freedom.

B. L. Koziey and F. A. Mirza

534

2

t

//’

Y X

Fig. 2. Nodal unit vectors and rotational degrees of freedom.

constraint assumes a less active role thus leading to larger errors. It is obvious from the discussion above that if a cubic were used for the transverse displacement w, the above mentioned inconsistency would disappear. Therefore, the use of a polynomial one degree higher for approximation of w than that used for CIand /I will always ensure a consistent formulation. This is the basis for the formulation of the consistent thick shell element in this paper. It should be noted that recognition of the unmatching modes is easy for complete polynomials which occur in the case of triangular elements. The use of incomplete polynomials for displacements and rotations for rectangular elements, e.g. the nine-node isoparametric shell element, makes it extremely difficult to identify such lack of matching. This might very well have been the reason for the delayed recognition of the problem with the eight and nine-node isoparametric shell elements.

CONSISTENT

THICK SHELL ELEMENT

FORMULATION

Co-ordinate systems and geometry

Co-ordinate systems used in the formulation of the consistent thick shell element are shown in Fig. 1, and are defined in the following manner. (1) Global Cartesian co-ordinate system x, y and z and the corresponding global displacements u, v and w. (2) Curvilinear co-ordinate system (r, s, t), r and s are tangent to the surface defined by t = constant, and t is not necessarily normal to the r-s tangent plane. (3) Local Cartesian co-ordinate system x’, y’ and z’ to define local strains and stresses, z’ is normal to the surface defined by t = constant and x’ and y’ are tangent to the surface. For geometric distortion of the element a curvilinear transformation is used in terms of the parent r and

Fig. 3. Cross-sectional deformation in the L-x’ plane due to rotations a, and 4,.

Consistent thick shell element

535

Displacement field

0.8 0

100

200

I

1

I

300

400

500

Length/Thickness

Fig, 4. Normalized center plate deflection vs width to thickness ratio for a clamped square plate subjected to a uniformly distributed load. s co-ordinates, and a linear transformation is used in terms of the t co-ordinate. The location of any point within the element in the global co-ordinate system is determined by the co-ordinates of the nodal points (x,, y,, z,) and vector VX, at each node as

To construct cubic approximations for global displacements u, v and W, the displacement degrees of freedom at the corner nodes, one-third side nodes and center element node are used. Quadratic approximations of rotations c(, /S, C#J and cp, defined below, are achieved using the corresponding rotational degrees of freedom at the corner and mid-side nodes as shown in Fig. 2. Rotations c( and 4 are about the local y/-axis and rotations B and cp are about the local x/-axis as shown in Fig. 2. Thus, the number of degrees of freedom per node is different, as shown in Fig. 1. Assuming that each node has seven degrees of freedom (u,, v,, w,, GI,,B!, $,, cp,) then there are 91 degrees of freedom per element. However, only 54 of these are active. It should be noted that rotations LXand /I are constant through the depth of the element while rotations 4 and cp vary quadratically. Thus, doand /I’ provide a linear variation of displacements u, v and w along t while C$and cp lead to a cubic variation of U, v and w as depicted in Fig. 3. The global displacements (u, v, W) can now be written in terms of the element degrees of freedom as

(6)

where quadratic interpolation functions N, are given in Appendix A and vector V,,, which connects the upper and lower points of the shell at the ith-node, is defined as (7)

(8) where N, and fl# are the quadratic and cubic shape functions, respectively, and are given in Appendix A. Matrix [E] = [V,,, - V,,] where unit vectors V,, and V,, are directed along the x’- and y’-axes, respectively, and form an orthogonal basis with unit vector V,, at

L = 120.0

q(x) = llpmw

H = 60.0 unit width

s, = looo.o E=

10’

v = 0.25

L

I

Fig. 5. Dimensions, loading and material properties of simply supported semi-infinite thick plate.

536

B. L. Koziey

and F. A. Mirza

.----O-Node

Fig. 6. Semi-infinite

thick plate normal

the ith-node. The procedure employed for construction of this basis is described in Appendix B. Shape functions MI and Mz approximate the displacement field through the depth due to the constant and quadratic rotations, respectively, and are given by

where h, is the shell thickness at the ith-node. The derivation of the through thickness displacement field due to quadratic rotations Q, and cp is given in Appendix C. Note that the consistent thick shell element is subparametric since a, v and u’ are approximated cubically while the geometric transformations are performed using quadratic expressions. Strain-displacement and stress-strain relationships

For small deflections the strains (6’) along the local axes are given by

stress distributions

where u’, v’ and w’ are displacements directed along the local x’, y’ and z’ axes, respectively, and [L] is the linear operator matrix. The strain in the z’ direction is ignored for obvious reasons. Transformation matrix [O] is established using the Jacobian matrix as described in Appendix B. It contains the direction cosines between the local and the global axes at the point where the displacements are to be transformed and is given as

x

z

1

(11)

Using the transformation matrix [8] and the expression for the global displacements, eqn (8), the displacements along the local axes in terms of the

FF

0

0

0

va

0

111 W’

=

Y

1, ml n, I2 m2 n2 z’ 1, m3 n3 L

x’ [13]= y’

a

u’ = [L] v’

at x = L/4.

a

v

a

22 0

a 32

U

O a

O ax;

a &Fa

3-7

v

WI ii

W

(10)

Consistent thick shell element

z

537

0.5

H

0.4 0.3 0.2 0.1 0.0 .-0.1 .-0.2 .-0.3 w-o.4 , .-0.5 Fig. 7. Semi-infinite thick plate shear stress distributions at x = 0.

thickness = 0.25 v = 0.0 E = 4.32 x 10’ weight = 90 per unit area Fig. 8. Cylindrical shell roof dimensions, loading and material properties.

Table 1. Cylindrical shell roof (Ref. Defln = 3.610) CONSHL Mesh 1x1 2x2 4x4 6x6 8x8

Elements, 2 8 32 -

Defln A 2.28660 3.10276 3.51078 -

Norm. 0.633 0.859 0.973 -

9 Node Isoparametric d.o.f. 46 162 604 -

Elements 4 16 36 64

Defln A 0.96239 2.90531 3.38604 3.48602

Norm.

d.o.f

0.267 0.805 0.938 0.966

96 352 768 1344

538

B. L. Koziey

and F. A. Mirza

R = 25 in. L = 50 in. t = 0.25 in. P = 10 lbs. E=3OOOOksi v = 0.3

Fig. 9. Pinched

element

cylindrical

shell dimensions,

loading

and material

properties.

the required differentiation yields the desired relationship between the local strains and the element degrees of freedom as

degrees of freedom are given by

(6’) =

PI(d)

(13)

where {d} is the vector of element degrees of freedom matrix [B’] is given in Appendix D. The following form of the generalized Hooke’s Law has been employed

(12) and the strain+lisplacement where

= [D’][{t’} - {EL}1+ {06} (14)

=

(-l,h,

-

mlm2,

-

nln2,)

mlml,

+

wz,,)

(-l,l,,

-

m2m2#

-

n2n21)

m3mh

+

n3nlz)

(-IA

-

m3m21

-

n,n,,)

(1111, + rnlml, + ntn,,) (12lh

+

WI, + Substituting

eqn (12) into eqn (10) and carrying

Table 2. Pinched

cylindrical

I

out

where (66) and {CA}are the initial local strain and . . local stress vectors, respectively.

shell, diaphragmed

ends (Ref. Defln = 0.218987) 9 Node Isoparametric

CONSHL Mesh 2x2 4x4 6x6 8x8 10 x 12 x 14 x 16 x 18 x 20 x 22 x 24 x 26 x

Elements

10 12 14 16 18 20 22 24 26

8 32 72 128 200 288 392 512 648 -

Defln

Norm.

d.o.f.

Elements

Defln

Norm.

d.o.f

0.02785 0.11367 0.16921 0.19308 0.20439 0.21031 0.21366 0.21570 0.21703 -

0.127 0.519 0.773 0.882 0.933 0.960 0.976 0.985 0.991 -

152 578 1284 2293 3536 5082 6908 9014 1 1,400

4 16 36 64 100 144 196 256 324 400 484 576 676

0.01089 0.03530 0.07965 0.12345 0.15596 0.17676 0.18976 0.19813 0.20376 0.20769 0.21054 0.21266 0.21428

0.050 0.161 0.364 0.564 0.712 0.807 0.867 0.905 0.930 0.948 0.961 0.971 0.978

86 330 734 1298 2022 2906 3950 5154 6518 8042 9726 11,570 13,574

Consistent thick shell element

539

Table 3. Pinched cylindrical shell, free ends (Ref. Defln 5.4236) CONSHL

9 Node Isoparametric

Mesh

Element 3

Defln

Norm.

d.o.f.

Elements

Defln

Norm.

d.o.f

2x2 4x4 6x6 8x8 10 x 12 x 14 x 16 x 18 x 20 x

8 32 72 128

3.17303 4.88251 5.07955 5.12217 -

0.585 0.900 0.937 0.944 -

165 607 1329 2331 -

4 16 36 64 100 144 196 256 324 400

0.02909 0.33731 1.32131 2.67958 3.74367 4.37334 4.71702 4.90646 5.01516 5.08053

0.005 0.062 0.244 0.494 0.690 0.806 0.870 0.905 0.925 0.937

95 351 767 1343 2079 2975 4031 5247 6623 8159

10 12 14 16 18 20

-

-

The elasticity matrix [D’], for IJ;,;,= 0, is of the following form: 1

P’l=j$ L

v

V

1

0

02

0 0

0 0

0

0

0

0

0 0

1-vo

0

0 oy OF0

1

1 (‘5)

The above equation is integrated numerically in terms of the non-dimensional co-ordinates (r, s, t) using the Gaussian Quadrature scheme with seven integration points in the r--s plane and five integration points along t. Before integrating eqn (16) the partial derivatives with respect to the local co-ordinates (x’,y’, z’) appearing in [B’] must be transformed to derivatives in terms of the non-dimensional co-ordinates (r, s, t). This is accomplished using the following Jacobian matrix [J’]

where E and v are the elastic modulus and Poisson’s ratio, respectively. It should be noted that the shear correction factor K is not required in [D’] because of the parabolic approximation of the transverse shear strains through the thickness. Furthermore, at this stage we can incorporate any general constitutive equation, e.g. for any orthotropic or anisotropic material. Stiffness matrix

The element stiffness matrix [k] is derived in the usual manner as

[k] =

[B’]T[D’][B’]dJ’. sA

a Z a

ad ayf af --ar ar ar ax< ayf az’ as-dsBsZ a ax’ af aZf at star-z

1;lI(i

8X’ axI ax ax’ ay aX( az -=zar+3jG+azar ar

=/&+m!2+naZ ’ ar

A-,k-

v ?j

0.8

$ Q 0.4

$ z

0.4

t p 0.2

E z 0.2

0

3000

8000

9000

12000 15000

Degrees of Freedom Fig. 10. Normalized deflection under applied load P vs number of degrees of freedom for pinched cylindrical shell with diaphragmed ends.

’ ar

’ ar. (‘8)

_&--A

P’

.P = 0.8

$ 0.8

(17)

The components of [J’] are found using eqn (6) and the matrix of direction cosines [8]. For example the first term of the matrix is given by

(16)

.g 0.8

-a ax' -a ayf. a 32

0

2500

5000

7500

10000

Degrees of Freedom Fig. 11. Normalized deflection under applied load P vs number of degrees of freedom for pinched cylindrical shell with free ends.

B. L. Koziey and F. A. Mirza

540

R=

10 in.

t = 0.2 in. P = 100 lbs. E = 30000 ksi v = 0.3

Fig. 12. Thin pinched spherical shell dimensions, loading and material properties.

components of [J’] are determined in a similar manner. Finally, the differential volume d V can be written in terms of the curvilinear co-ordinates as dV= ]]J]]drdsdt (19) The remaining

where )I.I 11is the determinant

of the Jacobian matrix

[Jl. It should be noted that an eigenvalue test of the element stiffness matrix [k] was performed. The results verified the existence of the required rigid body modes and the absence of any spurious zero energy deformation modes. In addition, the validity of the formulation was further verified through constant strain and constant curvature patch tests.

Yielding occurs when the yield function F, eqn (20), is equal to zero i.e. when the effective stress d equals the uniaxial yield stress by. Beyond the proportional limit the following incremental stress-strain relationship is employed {da’} = [D&]{dc’}

(23)

where [D&] is the instantaneous elasto-plastic constitutive matrix in the local co-ordinate system. For an associated flow rule and isotropic hardening it takes the following form [21]

[D&l= [D‘I - [D‘I

Elasto-plastic constitutive relationship

The onset of plastic deformation is determined by the yield criterion. The yield function F corresponding to the Huber-von Mises yield criterion and a general state of stress is given by F=&-

ay=d-oy=O

(20)

where cry and 6 are the uniaxial yield stress and the effective stress, respectively. The second stress invariant J2 can be written in terms of the deviatoric stress tensor s, as J2 = f s,s,

SV= s, - &a,,,

1 a,, = - u,, . 3

(21)

(22)

where H’ is the plastic modulus of the slope of the uniaxial stress-strain curve beyond the yield limit. For use of a bilinear stress-strain curve for elasto-plastic material behaviour the plastic modulus is given by H’ = ET/(1 - ET/E)

(25)

where E and ET are the elastic and tangent moduli, respectively. After substituting Jz, eqn (21), into the yield function F, eqn (20), the vector {c~F/%} can be

Consistent thick shell element obtained as

has converged to the desired tolerance which occurs when

(1

aFT 3 = GS” aa

@$#
which, with a,.,, = 0 for plate and shell applications, takes the following form cry,,,- U,Y/2

{}I 1 aF 1 a0 :=d

6,,,, - bYY/2 37,.,,

.

(27)

37,y

37,y

IMPLEMENTATION

Non-linear finite er’ement analysis procedure

Using the principle of virtual displacements the discretized equations of equilibrium can be written as

s

[B’]‘{a’);dV Y

- {R},, = 0

(28)

where {R} is the applied load vector. In an incremental analysis the total load {R} acting on the structure is applied in increments. Because of the impending non-1inr:arity a certain unbalance of forces will exist between the external and the internal forces. This unbalance is minimized to a tolerable magnitude by employing an iterative procedure. For the mth iteration during the nth load increment eqn (28) can be written as

{@}:: = 5 {B’}T{Aa’}:dV - {AR}* )I

(29)

where {@}; is the unbalanced force vector called the residual load vector. In this study a modified Newton-Raphson iterative method is adopted. The displacement increment corresponding to the residual vector {@I:: is calculated as {Ad}::’ = [&(a’):]-‘{@}::

n

(31)

which compares the mth increment in the internal energy (i.e. the amount of work done by the out of balance loads on the displacement increments) to the initial internal energy increment, Bathe [22]. A typical value for TOL is between 1 and 5%. The residual forces are added to the load vector for the following load increments in order to prevent excessive drifting. Reduced strain-displacement matrix

Having obtained all of the terms required, the elastic-plastic constitutive matrix [O&l in eqn (24) can now be computed. NUMERICAL

541

As mentioned before, for the same number of degrees of freedom per node there are ninety-one degrees of freedom per element. However, the number of active degrees of freedom per element is only 54. When calculating the element stiffness matrix through multiplication of [B’lTID][B’] the presence of the inactive degrees of freedom results in many unnecessary multiplications with zero. This leads to a significant increase in the computer time required for solution of a problem and also can effect accuracy. However, by eliminating the columns of matrix [B’] corresponding to the inactive degrees of freedom unnecessary multiplications can be avoided. Re-writing eqn (16) in terms of the column reduced strain-displacement matrix [&‘I gives

[k,] =

[B:lTID’][B:ld I’ sY

(32)

where [BJ and [k,] are the column reduced strain-displacement and element stiffness matrix, respectively. Obviously, the reduced element stiffness matrix is 54 x 54 instead of 9 1 x 91. This approach requires minimal modification of the program code. Prior to evaluation of eqn (32) matrix [&‘I is formed using the appropriate destination vectors. After multiplication and integration, matrix [k,] is then expanded using the destination vectors to its full dimension. Implementation in this manner means that existing code need not be modified to allow for nodes which have different numbers of degrees of freedom. APPLICATIONS

(30)

where [&(a’):] is the tangential stiffness matrix evaluated at the beginning of each load increment. Using the incremental displacements the corresponding strain increment {AC’): and stress increment {AC’}; are computed. The displacement, strain and stress increments are added to the levels at the end of the previous iteration. Equilibrium is checked using eqn (28) and the iterations proceed until the solution

Convergence to thin plate solution

A clamp>d square plate subjected to a uniformly distributed load has been analysed. Due to double symmetry only one quarter of the plate has been modelled using 8 (2 x 2) consistent thick shell elements. This problem has been used in the literature as a test for shear locking. For this purpose the ratio of the plate width to thickness is increased and the predicted solutions are checked for convergence to

B. L. Koziey and F. A. Mirza

542

+

01 0

g

O-Node

I

I

1

I

I

1000

2000

3000

4000

5000

Number of Degrees of Freedom Fig. 13. Normalized vertical deflection at the boss vs number of degrees of freedom for thin pinched spherical shell.

the theoretical solution. In the present study plates having width to thickness ratios ranging from 10 to 500 have been analysed. The predicted deflections of the center of the plate have been normalized with respect to the theoretical solutions given by Timoshenko [23] and plotted in Fig. 4. Also plotted in the figure are results obtained using 16 eight-node isoparametric shell elements and 16 nine-node isoparametric shell elements with full 3 x 3 quadrature in the r-s plane. As the length to thickness ratio increases the transverse deflection predicted by the consistent thick shell element model converges to the thin plate solution. However, the eight- and nine-node isoparametric shell element models predict displacements which are significantly less than the reference solution. This is due to the well known shear locking phenomena. The provision of a consistent formulation eliminates the spurious transverse shear strain modes which plaque the eight- and nine-node shell elements. The analyses of the plate have also been carried out using the full strain-displacement matrix (FSDM) version and the reduced strain-displacement matrix (RSDM) version of the consistent thick shell element. For the FSDM implementation the computer time required for analysis was 28 s and for the RSDM version the computer time required was 18 s. Thus, by employing the reduced strain-displacement matrix a 36% reduction in the computer time required for problem solution can be achieved.

shown in Fig. 5. Half of the plate has been modelled using eight consistent thick shell elements by making use of symmetry. The plate segment has also been analysed using 4 nine-node isoparametric shell elements. The predicted mid-span deflection and stress distributions are compared with those determined using an elasticity solution by Little [24]. The predicted normal stress and transverse shear stress distributions through the plate thickness are plotted in Figs 6 and 7, respectively. Also plotted are the stress distributions given by the elasticity solution. As can be observed the normal stress and transverse shear stress distributions predicted by the consistent thick shell model are in excellent agreement with the elasticity solution unlike the stress distributions predicted by the nine-node isoparametric shell element. Furthermore, the variation of the transverse shear stress along the length of the plate was found to be free of spurious transverse shear modes. This, however, was not the case for the nine-node isoparametric shell element. The mid-span deflection predicted by the consistent thick shell

1.2

0.8 0.6

Semi-infinite thick plate

0.4

An analysis of a simply supported semi-infinite thick plate has been carried out to illustrate the ability of the consistent thick shell element to predict the normal stress and the transverse shear stress distributions through the depth of the plate. The material properties, dimensions, and loading are as

0.2 0 0

0.02

0.04

0.06

Fig. 14. Applied load vs center plate displacement for elasto-plastic analysis of thick clamped circular plate.

Consistent thick shell element

543

converges much more rapidly than the nine-node shell element. Approximately twice as many degrees of freedom are required by the nine-node shell element model to obtain results similar to those predicted by the consistent shell element. Because shear and membrane locking have been eliminated from the consistent shell element by ensuring consistency in the formulation it demonstrates good accuracy even for coarse meshes. Pinched cylindrical shell 0

1

2

3

4

5

The thin cylindrical shell depicted in Fig. 9, which is loaded by two centrally located and diametrically opposed concentrated forces is analysed. Two type of boundary conditions are considered; (i) the ends are covered with rigid diaphragms which allow displacement in the axial direction and rotation about the tangent to the shell boundary, and (ii) the ends are completely unrestrained. Using the double symmetry only one-eighth of the cylinder has been modelled using consistent shell elements and nine-node isoparametric shell elements. The displacement under the point force for both boundary conditions predicted by the finite element models are compared to the reference solutions. For the rigid diaphragm case the exact displacement is taken as 0.2189866 which has been computed by Lindberg et al. [25] based on Flugge’s equations. For short and thin cylinders with free ends a good approximation of the displacement, 6 under the force is given by [26]

6

Center Plate Displacement Fig. 15. Applied load vs center plate displacement for elasto-plastic analysis of thin clamped circular plate.

model and the nine-node isoparametric shell model are normalized with respect to the deflection given by the elasticity solution. The deflection predicted by the consistent thick shell element model is 1.09, which is considerably better than 1.20 as predicted by the nine-node isoparametric shell element. Cylindrical shell mof In this example the cylindrical panel shown in Fig. 8 is analysed. The shell is loaded by its own uniform dead weight and is supported by diaphragms at the ends, but is free along the sides. Using symmetry only one-quarter of the panel has been modeled using the consistent thick shell element. For the purpose of comparison the shell quadrant has also been analysed using the nine-node isoparametric shell element. The vertical deflection at point A predicted by the finite element models is compared to 3.610 in Ref. [7], which has been obtained from a deep shell theory solution. The results of the analysis using progressively finer meshes are given in Table 1. It should be noted that the nine-node isoparametric shell element used full 3 x 3 quadrature in the r-s plane. A comparison of the number of degrees of freedom shows that the consistent thick shell element 2 0.5 r

6 = 0.0745PR’ 2DL

(33)

where R is the radius, L the half-length, I the thickness, D is the flexural rigidity and P the magnitude of the applied force. Using the geometry and material properties for the cylinder given in Fig. 9 the reference displacement is found to be equal to 5.4236. This value has been found to be in good agreement with experiments [I 11. The results of the analysis of the cylinder with diaphragmed ends and

I,/

‘i

7rz (W

-

0.1

“““‘I’_ 0.0

0.5

0.4

0.3

0.2

0.i : #’

-----

Axisymmetric CONSHL

#’

__-- _* , -.__

OG

-0.1

- -0.2 - -0.3 - -0.4 -0.5

Fig. 16. normal stress and shear stress distributions at clamped edge of plate for q = 0.625 ksi.

544

B. L. Koziey and F. A. Mirza

Trz (ksi) II 1.0

/’

t

-0.4

-

I 0.8

I,

0.6

.- - - - Axisymmetric CONSHL

-0.5 -

Fig. 17. Normal stress and shear stress distributions at clamped edge of plate for q = 1.02 ksi.

Thin pinched spherical shell

the cylinder with free ends are given in Tables 2 and 3, respectively, for progressively finer finite element meshes. These results are also shown graphically in Figs 10 and 11, respectively, and show the variation of the normalized deflection with the number of degrees of freedom. The cylinder with free end diaphragms presents a more severe state of bending and membrane actions than the cylinder with free ends. This is demonstrated by the consistent shell element results. The results show that in order to obtain a similar level of accuracy in the predictions for the different boundary conditions requires that about twice as many degrees of freedom be used for the diaphragmed end case as compared to the case with the free ends. This however is not the case for the nine-node isoparametric shell element. The results for the nine-node shell element show similar convergence for both cases i.e. diaphragmed and free. This is because the nine-node shell element suffers from shear and membrane locking regardless of the boundary conditions. The consistent shell element which is free of spurious modes converges more rapidly to the reference solution for the cylinder with free ends as expected.

The spherical shell shown in Fig. 12 is subjected to concentrated loading at each pole. The load is transmitted to the shell through an infinitely rigid boss which has a half-angle /I0 = lo” as shown in the figure. The shell is assumed to be rigidly attached to each boss, and hence only vertical motion is allowed along the boundary. Using the symmetry about the equator and the poles only one-eighth of the shell has been modelled using consistent shell elements and nine-node isoparametric shell elements. Lindberg et al. [25] obtained an exact solution for the vertical displacement at the boss through integration of Timoshenko’s equations [23]. For the geometry and material properties given in Fig. 12 the exact deflection is equal to 3.5849. The shell was analysed using progressively finer meshes and the predicted vertical deflection at the boss, normalized with respect to the exact solution, is plotted in Fig. 13. In this example the elements are doubly curved unlike the previous examples. Again the vertical displacement predicted by the consistent shell element model converges quickly to the reference solution unlike the

z 0.5 To.4

0.4 H

0.3

0.3

0.2

0.2

7rz (ksi)

I

’

0.5

I 0.4

*

I

1

0.3 - -0.2

.---- Axisymmetric

-

CONSHL

Fig. 18. Normal stress and shear stress distributions at R = 450in for q = 0.625 ksi.

Consistent thick shell element

Axisymmetric CONSHL

Fi,g. 19. Normal stress and shear stress distributions at R = 450 in for q = 1.02 ksi.

nine-node isoparametric shell model predictions which are plagued by the presence of the spurious transverse shear strain modes. Elasto-plastic ana(ysis of a thin and a thick circular plate

Elastic-perfectly plastic analyses of thin and thick clamped circular plates subjected to uniformly distributed loads are performed using 85 consistent thick shell elements to model one-quarter of the plate. The radius, elastic modulus, Poisson’s ratio, and yield stress for both plates are 500 in, 30,000 ksi, 0.3, and 1.0 ksi, respectively. The ratio of the radius to thickness for the thin plate is 250, and for the thick plate 2.5. The predicted central plate deflections, normal stresses 6,, , and transverse shear stresses t,; are compared to the results obtained using the cubic displacement field ten-node axisymmetric triangular elasto-plastic element. Six axisymmetric elasto-plastic elements through the thickness 20 and in the radial direction have been employed for a total of 240 elements.

1.2 5 .-

I

I

Plots of the applied loads vs the central plate deflections for the thick and thin plates are shown in Figs 14 and 15, respectively. The responses predicted by the consistent shell element model are in very good agreement with those given by the axisymmetric model. The distributions of the normal stress, u,, and shear stress, T,: through the depth of the thick plate at its edge (i.e. R = 500 in.) are plotted in Figs 16 and 17 for load levels of q = 0.625 ksi and q = 1.02 ksi, respectively. For q = 0.625 ksi it is observed that the consistent thick shell element model predicts a normal stress distribution which is very close to that given by the axisymmetric elements. Even after a significant amount of yielding of the cross-section (i.e. at q = 1.02 ksi) has occurred the consistent thick shell element model predicts a normal stress distribution which is very close to that given by the axisymmetric model. For the shear stress distributions however, it is observed that there are significant differences between the predictions from the two models at both load levels. It is believed that if the meshes used in

I

I

I

1.0

0

200

400 Number

600 of

800

1000

1200

D.0.F

Fig. 20. Convergence of central plate deflection for a clamped thin square plate subjected to a u.d.1. CAS 65/4--c

B. L. Koziey and F. A. Mirza

546

I

I

I

I

Z I

I

I

I

50

100

150

200

0.0

0

Number

of

250

D.0.F

Fig. 21. Convergence of central plate deflection for a clamped thick circular plate subjected to a u.d.1.

both models were refined near the edge of the plate, the predicted distributions would be more accurate. It must also be pointed out that the above stress distributions are right at the fixed supports and, hence, one of a very complicated nature. It is also possible that the two distributions, especially with some plastic regions having developed, can be inherently different because of the two different modelling procedures involved. The variations of the normal and shear stresses through the thickness of the thick plate away from the edge at R = 450 in are plotted in Figs 18 and 19 for loads of q = 0.625 ksi and q = 1.02 ksi, respectively. For both load levels the normal and shear stress distributions predicted by the consistent thick shell element model are observed to be in excellent agreement with the those given by the axisymmetric model. Consistent thick shell element convergence Consider an element similar to the consistent thick shell element, but which employs quadratic and linear approximations of the displacements and rotations, respectively, instead of cubic and quadratic approximations. In this example the consistent thick shell element is referred to as the cubic-quadratic shell element and its lower order counterpart will be referred to as the quadratic-linear shell element. The purpose of this example is to compare the convergence and accuracies of the refined and crude consistent thick shell elements to each other. A clamped thin square plate subjected to a uniformly distributed load is analysed using the cubic-quadratic and quadratic-linear consistent shell elements. The elastic modulus, Poisson’s ratio and length to thickness ratio of the plate are 30,000 ksi, 0.3 and 200, respectively. Due to symmetry only one quarter of the plate has been analysed. The central plate deflection predicted by the finite element models is normalized with respect to the solution given by

Timoshenko [23] and plotted against the number of degrees of freedom in Fig. 20 along with the results obtained from the nine-node isoparametric shell element. For each finite element model, subsequent analyses have been carried out after uniform refinement of the previous mesh. The cubic-quadratic model converges more rapidly to the reference solution than the nine-node isoparametric shell element and the quadratic-linear model. It should be noted that the refined and the crude consistent shell elements lead to convergent results. More importantly, it demonstrates the consistency of formulation using two different polynomials for approximations. Furthermore, the faster convergence of the cubicquadratic consistent shell element than the nine-node isoparametric shell element is also free of any spurious shear stress modes. In addition, stiffer results are expected beyond the length to thickness ratio of 200 using the nine-node isoparametric shell element. A clamped thick circular plate subjected to a uniformly distributed load has also been analysed using the cubic-quadratic and quadratic-linear shell elements. The radius, thickness, elastic modulus and Poisson’s ratio for the plate are 5 in, 2 in, 30,000 ksi and 0.3, respectively. The central plate deflection predicted by the finite element models is normalized with respect to the solution given by Young [27] and plotted against the number of degrees of freedom in Fig. 21 along with the results from the nine-node isoparametric shell element. Again the accuracy and rate of convergence of the cubic-quadratic shell element are excellent. CONCLUSIONS

Shell finite elements, which employ Mindlin plate theory to account for transverse shear deformation, experience an inconsistent formulation if polynomials of the same degree are used to approximate both the

Consistent thick shell element displacements and the rotations (e.g. nine-node isoparametric shell element). A consistent formulation is ensured for such elements by using a polynomial one degree higher for approximation of the displacement than that used for approximation of the rotation. The numerical results presented clearly illustrate that the provision of a consistent formulation eliminates both of the spurious transverse shear strain modes, and the overly stiff response predicted in the case of thin plates and shells. Thus, application of the reduced integration technique is not required, especially in the inelastic range. A new consist,ent thick shell element has been formulated. In the formulation of the element two new rotations, which vary quadratically through the thickness, have been introduced. This provides a quadratic variation of the transverse shear strains, and a cubic variation of the displacements through the thickness. Inclusion of the quadratic transverse shear strain variation eliminates the need for use of the shear correction factor K, and also accounts for warping of the cross section. Based on the results presented, it is concluded that introduction of the quadratically varying rotations allows for an excellent representation of the transverse shear stresses and normal stresses through the thickness of the element for thick plates and shells. Various numerical examples, carried out and presented, demonstrate the s,uperior performance of the newly developed consistent thick shell element including fast convergence to the known solutions and degeneration frorn thick to thin plates and shells without producing overly stiff results.

REFERENCES

I. Ahmed, S., Irons, B. M. and Zienkiewiz, 0. C., Analysis of thick and thin shell structures by curved finite element. Int. *I. numer. Methods Engng, 1970, 2, 419451. 2. Reissner, E., The effect of transverse shear deformation on the bending of elastic plates. J. appl. Mech., 1945, 12, 69-76. 3. Mindlin, R. D., Influence of rotary inertia and shear in flexural rotations of isotropic elastic plates. J. apd. __ Mech., 1951, 18#, 1031-1038. 4. Stolarski, H. and Belvtschko. T.. Shear and membrane locking in curved Ci eleme&s.’ Comput. Meth. appl. Mech. Engng.,

541

shell elements. Inr. J. numer. Methods Engng, 1989, 28, 385-414. IO. Briassoulis, D., The zero-energy modes problem of the nine

node

lagrangian

degenerated

shell element.

Compur. Srruct., 1988, 30, 1389-1402.

11. Vu-Quoc, L. and Mora, J. A., A class of simple and efficient degenerated shell elements-analysis of global spurious-mode filtering. Comp. Merhoris appl. Mech. Engng, 1989, 74, 117-175. 12. Huang, H. C. and Hinton, E., A new nine node degenerated shell element with enhanced membrane and shear interpolation. Inl. J. numer. Methods Engng, 1986, 22, 73-92. 13. Hinton, E. and Huang, H. C., A family of quadrilateral Mindlin plate elements with substitute shear strain fields. Comput. Struct., 1986, 23, 409431. 14. Choi, C. K. and Kim, S. H., Coupled use of reduced integration and nonconforming modes in quadratic Mindlin element. Inr. numer. Merhodr Engng, 1989,28, 1909-1928. 15. Lo, K. H., Christensen,

R. M. and Wu, E. M., A high-order theory of plate deformation. J. appl. Mech.

ASME, 1977, A4, 663-676. 16. Kant, T., Owen, D. R. J. and Zienkiewicz, 0. C., A refined higher-order Co plate bending element. Comp. Strucr., 1982, 15, 177-183.

17. Reddy, J. N., A simple higher-order theory for laminated composite plates. J. Appl. Mech., 1984, 51, 745-752. 18. Levinson, M., An accurate simple theory of the statics and dynamics of elastic plates. Mech Res. Commun., 1980. 7, 343-350. 19. Murthy, M. V. V., An improved

deformation

transverse shear theory for laminated anisotropic plates.

NASA technical paper, 1903, 1981. 20. Ren, J. G. and Hinton, E., The fnute element analysis

of homogeneous and laminated composite plates using a simple higher order theory. Communs appl. numer. Meth., 1986, 2, 217-228. 21. Zienkiewicz, 0. C., The Finrre Element Method, 3rd edn.

McGraw-Hill, London, 1977. 22. Bathe, K. J., Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ, 1982. 23. Timoshenko, S., Theory of Plates and Shells. McGrawHill, New York, 1959. 24. Little, R. W., Elasticity. Prentice-Hall, Englewood Cliffs, NJ, 1973. 25. Lindberg, G. M., Olson, M. D. and Cowper, G. R., New developments in the finite element analysis of shells. In Quar. Bul. Div. Mech. Engng and Nat. Areo. Esrab. National Research Council of Canada, Division of Mechanical Engineering, Vol. 4, pp. l-38, 1969. 26. Lukasiewicz, S., Local Loads in Plates and Shells. Sijthoff & Noordhoff, Alphen an den Rijn, 1979. 27. Young, W. C., Roark’s Formulas for Stress and Strain, 6th edn. McGraw-Hill, New York, 1989.

1983, 41, 275-290.

5. Zienkiewicz, 0. C., Taylor, R. L. and Too, J. M., Reduced integration technique in general analysis of plates and shells. Znt. J. numer. Methods Engng, 1971, 3, 275-290.

6. Pawsey, S. F. and Clough, R. L., Improved numerical integration of ihick shell elements. ht. J. number. Methods Engng. 1971, 3, 575-586.

7. Parish, H., A critical survey of the 9-node degenerated shell element with special emphasis on thin shell application and reduced integration. Comput. Meth. appl. Mech. Engng, 1979, 20, 323-350.

8. Hughes, T. J. R. and Cohen, M., The ‘Heterosis’ finite element for plate bending. Camp. Struct.,

1978, 9,

445450.

9. Belytschko, T., Wong, B. L. and Stolarski, H., Assumed strain stabilization procedure for the 9-node Lagrangian

APPENDIX A

Interpolarion jiuncfions N, and N, Quadratic interpolation functions N, N, = L,(2L, - 1)

Nj = b(2L,

- 1)

N, = L1(2L2 - 1) N4 = 0

Ns = 4L, L2

Ng = 0

Nt = 0

Ns = 4LzL,

Nq = 0

NW,= 0

N,, = 4L, L,

NII = 0

N,, = 0

B. L. Koziey and F. A. Mirza

548 Cubic interpolation functions

R,

APPENDIX

Derivation

m, = f L,(3L, - 1)(3L, - 2)

ml = ; L2(3L*- 1)(3L2 - 2)

I& = ; L,(3L, - 1)(3L, - 2)

9 N4 = z L, L2(3L, - 1)

iirS= 0

Ng = z9 L,L2(3Lz - 1)

N7 = -9 L*L,(3L2 - 1) 2

& = 0

C

of interpolation function

M2

Let rotation I$ vary quadratically through the thickness and be approximated by 4(r) =

a + 65 + cl’

(A5)

where 5 is the dimensional co-ordinate in the thickness direction and a, b and c are unknown constants. The displacement U due to 4 is calculated as L’(5) = r$(T)dT.

N,o = -9 L,L1(3L, 2

Ns = -9 L2L3(3Ln- 1) 2

- 1)

Substituting the expression for 4 m eqn (AS), into the above equation and integrating gives

9 Nu = Z L,L,(3L, - 1)

NII =0

(A@

U(t) = a{ + bi;’ + g + d 2 3

(A7)

N,, = 27L,LlL, L,, LZ and L, designate the area coordinates of the triangular

where d is an additional unknown constant. The following boundary conditions are imposed

parent element. APPENDIX

Construction

u h =o (2)

B

u -4 =o ( 2)

U(0) =o

4(O)= 4, (A8)

of orthogonal basis

The Jacobian matrix can be written as

(Al)

where 4, is the allowed rotational degree of freedom. Substituting the above boundary conditions into the approximations for 4 and U i.e. eqns (A5) and (A7). respectively, and subsequently solvmg for the unknowns yields a=&

b=O

c=-$

d=O.

where vectors R and S are tangent to the surface defined by t = constant. A vector V, normal to this surface is found as Vj = R x S.

(A2)

(A9)

Backsubstitutlon of the above into the approximations for 4 and U yield

The remaining vectors VZand VI of the orthogonal basis are given by VZ= Vj x R

V, = Vz x V,.

(A3)

Normalizmg VI, V? and V, gives the set of unit vectors V,, 0, and V, from which the transformation matrix of direction cosines is constructed as

Converting the dimensional co-ordinate 5 to dimensionless co-ordinate t using 5 = ht/2 enables eqns (AIO) to be re-written in terms of t as u(t) = ?(I

f&r) = (1 - 3tJ)f$,

- P)f#I,= M$,.

(All)

Interpolation function Ml defines the displacement variation through the thickness due to rotation 4.

APPENDIX

Consistent

I

D

shell element strain-displacement

,a

'ai

@

*af

I ,

m’ax’

I

m2m

aN,

matrix

I ,

ma

I

n2w

[B’]

a& afiT

afit

I I I

Consistent thick shell element I I I

I )

I

1 C;‘(M2~+~N,)

aN,

+C;’ )

Mzu+-gN,

ahh )

(

1 C~‘~2~+~N,)+C:‘~2~+~N.)

1

aN#

+C;’

)

Mzdy’+‘N,

(

aM2

aY

)

1 c:‘~2~+~N,)+C:‘(M2~+~N.)

I C? (

aN, ahf2 M2ayl+vN,

+C?

M2$+$$N, (

)

I i I

(

)

I I