A new plate finite element using integral representation interpolation

Engineering Analysis with Boundary Elements 23 (1999) 549–554

A new plate finite element using integral representation interpolation G. Bezine*, B. Millet Laboratoire de Me´canique et Physique des Mate´riaux, UMR CNRS 6617-ENSMA, Site du Futuroscope, Chasseneuil du Poitou, BP109, 86960 Futuroscope Cedex, France Received 5 February 1998; accepted 18 February 1999

Abstract The purpose of this article is to propose a quadrilateral plate bending finite element with eight nodes. An integral representation of the second derivatives over the quadrilateral element is developed in terms of the deflection and the normal slopes. This representation allows us to obtain a displacement interpolation and consequently expresses the strain energy in terms of the usual degrees of freedom at the corner and the midside nodes, leading to the computation of an adequate stiffness matrix. Results prove the accuracy of the method and the efficiency of an interpolation on an element using integral representations. q 1999 Elsevier Science Ltd. All rights reserved. Keywords: Boundary integral equation; Kirchhoff’s plate bending; Finite element method; Stiffness matrix

1. Introduction It is well known that plate bending elements built with the classical Kirchhoff’s theory do not reliably preserve the continuity of the deflection w and of its two spatial derivatives at the interface of the neighbouring elements. Moreover, convergence criteria require that elements satisfy completeness conditions such as rigid body and constant curvature representations. By taking polynomial function for deflection interpolation, the choice of the deflection and of the two slopes as degrees of freedom (DOF) automatically verify the displacement continuity but does not respect the normal derivative continuity. To satisfy this requirement, refined triangular and rectangular elements have been developed by subdividing each element into several sub-triangles with different polynomial approximations [1,2]. In addition, it is necessary to take into account the second derivatives as DOF. The purpose of this article is to present the development of the stiffness matrix for a quadrilateral element comprising eight nodes and twenty DOF, which satisfies the continuity of displacement, and of normal and tangential derivatives at the interface between two elements. This element is obtained by introducing an integral representation as interpolation of the deflection and slopes over the domain. Thus, a linear piecewise approximation along the * Corresponding author. Tel.: 1 33-05-4949-8037; fax: 1 33-05-49498238. E-mail address: [email protected] (G. Bezine)

side of the quadrilateral element is sufficient to ensure the continuity of w and its derivatives.

2. Basic equations Consider a thin elastic plate of thickness h, defined using a cartesian coordinate system Oxyz describing domain V. Consider then its midsurface S, described by z ˆ 0, which has a boundary G in the xy-plane and a plate deflection w(x,y) due to the transverse load per unit area p. For convenience, we assume that no traction is present in the xy-plane. Using the classical assumptions of Kirchhoff’s plate theory [3], the deformed state of the plate can be entirely described by the deflection w and its two derivatives 2w/2x and 2w/2y on the midsurface, which are respectively the rotations 2 u y and u x. Hence, the displacements ux and uy are given by: ux ˆ 2z

2w ; 2x

uy ˆ 2z

2w : 2y

…1†

Consequently the expressions of the strains are:

1x ˆ 2z

22 w ; 2x2

22 w : 21xy ˆ 22z 2x2y

0955-7997/99/$ - see front matter q 1999 Elsevier Science Ltd. All rights reserved. PII: S0955-799 7(99)00009-0

1y ˆ 2z

22 w ; 2y2

…2†

We can then deduce the bending strain energy for a

550

G. Bezine, B. Millet / Engineering Analysis with Boundary Elements 23 (1999) 549–554

Kirchhoff plate: Z z2 {w 00 }T ‰CŠ{w 00 } dV; W ˆ 12

…3†

V

where w 00 is the curvature vector such as {w 00 }T ˆ {…22 w=2x2 †; …22 w=2y2 †; 2…22 w=2x2y†} and [C] by taking into consideration the stress-strain relationships for an elastic isotropic material in the following matrix: 2 3 1 n 0 6 7 E 6n 1 0 7 …4† ‰CŠ ˆ 6 7; 2 …1 2 n † 4 12n 5 0 0 2 where E and n are respectively the elastic modulus and Poisson’s ratio. Moreover, if we consider the equilibrium equations, it can be proved that, according to Kirchhoff’s plate theory, the transverse deflection satisfies the differential equation [3] on S: DDDw ˆ p…x; y†

…5†

in S;

where D is the Laplacian operator and D the flexural rigidity defined by Dˆ

Eh3 : 12…1 2 n2 †

R which is reduced to dWe ˆ h3 =12 Se {dw 00 }T ‰CŠ{w 00 } dS after integration on the thickness h. In a standard finite element formulation, the deflection over each element can be expressed in terms of DOF (the vector {qe} in which there are w and its derivatives at all the nodes of the element) by mean of interpolation functions Ni …x; y† such as w…x; y† ˆ ‰N…x; y†Š{qe }. Then, the strain energy variation is given by: Z ‰N 00 ŠT ‰CŠ‰N 00 Š dS{qe }; …10† dWe ˆ d{qe }T h3 =12 Se

which defines the elementary stiffness matrix [ke] by: Z ‰N 00 ŠT ‰CŠ‰N 00 Š dS: …11† ‰ke Š ˆ h3 =12 Se

Usually, the interpolation over elements is performed by polynomial fields. Nevertheless, Kirchhoff plate elements require C 1 continuities for the transversal deflection w. Therefore 2w/2x and 2w/2y as well as w must be continuous across the element interfaces. In general this continuity of slopes and deflection at the interfaces cannot be satisfied except for very complex elements. Hence, in the following, we will propose an interpolation of displacement and slopes by a boundary integral representation, which satisfies the continuity conditions.

4. Integral formulation 3. Finite element formulation A classic finite element formulation can be established from a variational formulation. The most popular ones are the virtual work principle and the minimum of potential energy. For instance, if we use the latter, we have to minimise the total potential energy P of the system with respect to the displacement field u so that: dP ˆ 0

;du

…6†

in V;

where du represents the displacement variation, and Z Z Z d1ij sij dV 2 dui fi dV 2 dui pi dS: dP ˆ V

V

S

…7†

The first term is the total strain energy variation and the last two represents the variation of the potential energy of loads per volume (fi) or area unit (pi). By dividing domain V into finite elements Ve we can split the integrals into a summation of integrals over all the element domains X dP e …8† dP ˆ e

with dWe ˆ

Z Ve

z2 {dw 00 }T ‰CŠ{w 00 } dV;

…9†

The integral formulations of plate bending problems by the BIEM were first obtained on the Green’s identity [4] but, in the aim of taking into account the boundary’s discontinuities, several researchers have established the generalised Rayleigh–Green identity for the bending of Kirchhoff isotropic plates with arbitrary boundaries (for instance, Be´zine [5,6] and Stern [7]). However, the first formulation with Green’s identity [8,9] is easier to compute and does not introduce the additional unknowns Mnt(w). Consequently, to build a finite element, it is more effective to use this formulation. The integral representations of the bending analysis governed by the Eq. (1) are summarised thereafter. First, an integral representation of plate deflection inside the plate domain S is obtained [9] Z p 1 Z v dS 2 bw…P† ˆ D G S D   2Dv 2w 2v 2Dw w 2 Dv 1 Dw 2 v  ds …12† 2n 2n 2n 2n with b ˆ 1 if P [ S; b ˆ 0.5 if P [ G and for P located on a smooth part of the boundary G and b ˆ u=2p [7] for any corner point P of the boundary G for which the tangential discontinuity is u . A second equation is necessary, then the derivative is formed in a fixed direction z at the point P. If this direction

G. Bezine, B. Millet / Engineering Analysis with Boundary Elements 23 (1999) 549–554

551

the two directions, due to a sufficiently regular expression for w. Therefore, by taking two arbitrary but fixed directions z1 and z2 in Oxy, we can evaluate the expression of 2 2w/ 2z12z2 for a point P considered strictly inside S. This leads to the following integral representation:

Fig. 1. Orientation of tangent and normal vectors at the points P and Q of the boundary.

z is taken as n0, the following relation is deduced: Z p 2v 2w 1 Z b …P† ˆ dS 2 2n0 D G S D 2n0 " 2 # 2 Dv 2Dv 2w 22 v 2v 2Dw w2 Dw 2 1 ds × 2n0 2n 2n0 2n 2n0 2n 2n0 2n …13†

Z p 22 v 22 w …P† ˆ dS 2z1 2z2 S D 2z1 2z2 " 3 1 Z 2 Dv 22 Dv 2w 23 v w2 Dw 1 2 D G 2z1 2z2 2n 2z1 2z2 2n 2z1 2z2 2n # 22 v 2Dw ds: …14† 2 2z1 2z2 2n By taking successively: z1 ; x and z2 ; x; z1 ; x and z2 ; y; z1 ; y and z2 ; y it is possible to compute …22 w=2x2 †; …22 w=2x2y† and 22 w=2y2 . No singularity appears since r ± 0 and P is strictly inside domain S. The boundary integral Eqs. (12)–(14) are sufficient to express the strains for any plate flexure problem.

with b ˆ 0.5 for a point P located on a smooth part of boundary G and b ˆ u=2p for any corner point, where v ˆ …1=8p†r 2 log r with r ˆ iPQi; P and Q represents respectively a fixed point of G and a moving point of G. This function v(P,Q), called the fundamental solution, satisfies DDv ˆ d (P,Q) for which d is the Dirac function with its origin at point P. n is the outward normal at a point Q of G and n0 is the outward unit normal at a point P of G (Fig. 1). Du is the Laplacian associated with the deflection field u and p is the load per unit area. As a result, it is possible to express the deflection at any point as a function of the nodal deflections, slopes, Laplacian and Laplacian’s normal derivatives. Since the DOF are the deflections and the slopes, consequently it is necessary to eliminate the other unknowns. Some boundary integrals contains strongly singular integrands of order r 21 or r 22 that can be regularized with the first two terms of a Taylor expansion of the displacements about the source point [10]. Then the integrals are considered in a Hadamard principal-value sense [11].

Consider a quadrilateral element with eight nodes, which are at the four corners and at the four middles of the sides (Fig. 2). This finite element has three DOF (w, 2w/2x and 2w/2y) per corner node and two DOF per interface node (w and 2w/ 2n). Hence this element present twenty DOF. The interpolation of w, 2w/2x and 2w/2y are provided by the integral representations (12)–(14) with p ˆ 0 on each element. But as it is well known, in the boundary integral equation method, an interpolation of all the quantities w, 2w/2x, 2w/ 2y, Dw, 2Dw/2x and 2Dw/2y on the boundaries is necessary. This implies that over every segment of the boundary element, the deflection w, the normal derivative 2w/2n, the Laplacian Dw and its normal derivative 2Dw/2n are supposed to be varying linearly between their nodal values. Numerical integration of the rigidity matrix will be computed with k Gauss integration points inside the elements.

5. Evaluation of deflection’s second derivatives

7. Matrix formulation of boundary integral equation

Considering the integral representation (12) of w for P inside domain S, it is possible to differentiate (12) twice in

Using the approach proposed in [12], a matrix formulation can be performed. The element boundary G is divided into four straight segments (Fig. 2). At the middle points Qi, (nodal points Q2, Q4, Q6 and Q8) of which we define the nodal values of deflection w(Qi), normal slope …2w=2n†…Qi †, Laplacian Dw(Qi) and Laplacian normal derivative…2Dw=2n†…Qi †. Similarly, at each corner Qi (nodal points Q1, Q3, Q5 and Q7), we define the nodal values of deflection w(Qi), Laplacian Dw(Qi), as well the nodal values of the two normal slopes …2w=2n†…Qi † and the two

Fig. 2. Element description.

6. Finite element description

552

G. Bezine, B. Millet / Engineering Analysis with Boundary Elements 23 (1999) 549–554

Fig. 3. Deflection along the axis of symmetry.

Laplacian normal derivatives …2Dw=2n†…Qi † since in a corner there are two normal. Assumes that …2Dw=2n†, Dw; 2w=2n and w vary linearly between the two nodes on each straight segment, Eq. (12) written at each node Qi can be discretised in the following way:     2Dw 2w 1 ‰BG Š{Dw} 1 ‰CG Š 0:5{w} ˆ ‰AG Š 2n 2n 1 ‰DG Š{w};

(

2Dw 2Dw 2Dw 2Dw {Dw} ˆ Dw1 ; ; ; Dw2 ; ; Dw3 ; ; 2x 1 2y 1 2n 2 2x 3 2Dw 2Dw 2Dw 2Dw  ; Dw4 ; ; Dw5 ; ; ; 2y 3 2n 4 2x 5 2y 5 T

) 2Dw 2Dw 2Dw 2Dw Dw6; ; Dw7 ; ; ; Dw8 ; 2n 6 2x 7 2y 7 2n 8

…15†

where [AG] and [CG] are 8 × 12 matrices, [BG] and [DG] are 8 × 8 matrices. By a similar treatment Eq. (13) yields:       2w 2Dw 2w ˆ ‰A 0G Š 1 ‰B 0G Š{Dw} 1 ‰C 0G Š 0:5 2n 2n 2n 1 ‰D 0G Š{w};

and

…16†

where [A 0 G] and [C 0 G] are 12 × 12 matrices, [B 0 G] and [D 0 G] are 12 × 8 matrices. The coefficients of matrices [AG], [BG], [CG], [DG] and [A 0 G], [B 0 G], [C 0 G], [D 0 G] result from curvilinear integrals respectively of Eqs. (12) and (13). The vectors {2Dw=2n}; {Dw}; {2w=2n}, and {w} are the column vectors of quantities defined on the 8 nodal points of G. The subscript G indicates that the matrices are obtained for points P on the boundary G. By introducing the two vectors {w} and {Dw} such as: ( 2w 2w 2w 2w 2w T {w} ˆ w1 ; ; ;w ; ;w ; ; ;w ; 2x 1 2y 1 2 2n 2 3 2x 3 2y 3 4

2w 2w 2w 2w 2w 2w  ;w ; ; ;w ;w ; ; ; 2n 4 5 2x 5 2y 5 6; 2n 6 7 2x 7 2y 7 ) 2w  w8 ; ; 2n 8

the following system can be finally obtained: ‰GG Š{w} 1 ‰JG Š{Dw} ˆ {0};

…17†

where [GG] is a 20 × 20 matrix deduced from [CG], [DG], [C 0 G] and [D 0 G]; [JG] is a 20 × 20 matrix deduced from [AG], [BG], [A 0 G] and [B 0 G]; {w} is the vector whose 20 components are the 20 DOF of the finite element. The vector {w} can be identified as the vector {qe}. Therefore, we have [GG]{qe} 1 [JG]{Dw} ˆ {0}. Likewise Eq. (14) written inside domain S gives (with p ˆ 0) {22 w} ˆ ‰GS Š{w} 1 ‰JS Š{Dw}:

…18†

Thus, {22 w} ˆ ‰GS Š{qe } 1 ‰JS Š{Dw} where [GS] and [JS] are 3k × 20 rectangular matrices; {2 2w} is the vector whose 3k components are the three curvatures …22 w=2x2 †; …22 w=2x2y† and …22 w=2y2 † at the k Gauss points P inside S. Subscript S indicates the matrices are obtained for points P inside domain S. As …2Dw=2n† and Dw are not DOF in finite element formulation, a very convenient method consists in eliminating them in Eq. (17) by using Eq. (16) in the following form: {Dw} ˆ 2‰JG Š21 ‰GG Š{qe }:

…19†

Substituting Eqs. (19) into (17) leads to the resulting

G. Bezine, B. Millet / Engineering Analysis with Boundary Elements 23 (1999) 549–554

553

Table 1 Deflection at the load point for a square cantilever plate Mesh

2 × 2 elements l

3 × 3 elements l

4 × 4 elements l

Abaqus results (S8R) BIE finite elements results Difference

0.3905 0.4278 9%

0.3843 0.3973 3.4%

0.3901 0.3965 1.6%

linear system {22 w} ˆ ‰GS Š{qe } 2 ‰JS Š‰JG Š21 ‰GG Š{qe };

…20†

which can be condensed in the following form: {22 w} ˆ ‰HŠ{qe };

…21†

where matrix [H] of dimension 3k × 20 is given by [H] ˆ [GS] 2 [JS][JG] 21[GG]. Consequently it is possible to compute all the curvatures …22 w=2x2 †, …22 w=2x2y† , …22 w=2y2 † in function of the DOF vector {qe}. It follows that with Eq. (11) the rigidity matrix can be computed as: ‰ke Š ˆ h3 =12

k X

‰HŠT ‰CŠ‰HŠ;

jˆ1

where k is the number of Gaussian integration points in order to perform the numerical integration. It is clear that the rigidity matrix [ke] is a symmetric matrix, since it results from the product [H] T[C][H], where [C] is a symmetric matrix. 8. Numerical results Two examples are treated. For all the numerical results, the value of Poisson’s ratio is 0.3 and the plate is supposed to be a square plate of side length a. Moreover the origin of the axis Oxy is located at the centre of the plate. The first illustrative example consists of a square clamped plate with an uniform load. This square plate is divided into 64(8 × 8) elements. The value of E, n and h are such that the bending rigidity D is equal to the unit. Fig. 3 compares our numerical results for the displacement along a central transversal line with the analytical solution given by Lekhnitskii [13]. The second example treats the problem of a cantilever plate with a concentrated load applied at the middle of the free side opposite to the clamped side. The results computed with the BIE finite elements are presented below: the normalised deflections l at the loaded point …l ˆ w × 10110 =Da2 † are compared with Abaqus [14] results obtained with the same meshes for S8R eight node elements. The domain of the square plate is divided into 4(2 × 2) or 9(3 × 3) or 16(4 × 4) finite elements. Table 1 shows that the comparison with the existing solutions is quite satisfactory. The difference is always inferior

to 9% between our BIE finite element and Abaqus finite element results. Considering with the finite element method boundary conditions concern w, 2w\2x and 2w\2y, and considering the choice of the Green identity as reciprocity theorem, it is assured that no difficulties appears at the corners whatever boundary conditions.

9. Conclusion The boundary integral equation method proposed in this article is an effective numerical technique to obtain interpolation functions for displacement in Kirchhoff’s plate bending problems. The element computed by this method was successfully evaluated for clamped and cantilever plates. The results prove the accuracy of these eight nodes quadrilateral elements.

Appendix Expressions of the kernels of the integral representation of w(P): vˆ

1 2 r log r; 8p

2v 1 ˆ …2 log r 1 1†r cos…g 2 a†; 2n 8p Dv ˆ

1 …log r 1 1†; 2p

2Dv 1 ˆ cos…g 2 a†: 2n 2pr Expressions of the kernels of the integral representation of 2w/2n: 2v 1 …2 log r 1 1†r cos…g 2 u†; ˆ2 2n0 8p 22 v 1 ˆ2 ‰…2 log r 1 1† cos…u 2 a† 1 2 cos…g 2n0 2n 8p 2 a† cos…g 2 u†Š;

554

G. Bezine, B. Millet / Engineering Analysis with Boundary Elements 23 (1999) 549–554

2Dv 1 cos…g 2 u†; ˆ2 2n0 2pr

23 Dv 1 ˆ 2 3 cos…3g 2 a†; 2 2y 2n pr

22 Dv 1 ˆ cos…a 1 u 2 2g†: 2n0 2n 2pr2

23 Dv 1 ˆ sin…3g 2 a†; 2x2y2n pr 3

Expressions of the kernels of the integral representation of …22 v=2x2 †, …22 v=2x2y†, …22 v=2y2 †: 22 v 1 …2 log r 1 1 1 2 cos2 g†; ˆ 2 8p 2x 22 v 1 …2 log r 1 1 1 2 sin2 g†; ˆ 2 8p 2y 22 v 1 ˆ sin g cos g; 2x2y 4p 23 v 1 ‰cos a cos g 1 sin g sin…2g 2 a†Š; ˆ 2 4pr 2x 2n 23 v 1 ‰sin a sin g 1 cos g cos…2g 2 a†Š; ˆ 2 4pr 2y 2n 23 v 1 ˆ ‰sin a cosg 1 sin g cos…2g 2 a†Š; 2x2y2n 4pr 22 Dv 1 ˆ2 cos 2g; 2x2 2pr2 22 Dv 1 ˆ cos 2g; 2y2 2pr2 22 Dv 1 ˆ2 sin 2g; 2x2y 2pr2 23 Dv 1 cos…3g 2 a†; ˆ 2x2 2n pr3

References [1] Clough RW, Tocher JL. Finite element stiffness matrices for analysis of plate bending, Matrix methods in structural mechanics, AFFDL TR-66-80, November 1965, p. 515–545. [2] Fraeijs De Veubeke B. A conforming finite element for plate bending. Int. J. Solids Struct. 1968;4(1):95–108. [3] Timoshenko S, Woinowsky-Krieger S. The´orie des plaques et Coques, Librairie Polytechnique Ch. Beranger, Paris, 1961. [4] Bergman S, Schiffer M. Kernel functions and elliptic differential equation in mathematical physics. New York: Academic Press, 1953. [5] Be´zine G, Gamby D. A new integral equation formulation for plate bending problems. In: Brebbia CA, editor. Recent advances in boundary element methods, London: Pentech Press, 1978. pp. 327–342. [6] Be´zine G. Boundary integral formulation for plate flexure with arbitrary boundary conditions. Mech. Res. Comm. 1978;5(4):197–206. [7] Stern M. A general integral boundary formulation for the numerical solution of plate bending problems. Int. J. Solids Struct. 1979;15:769–782. [8] Maiti M, Chakrabarty SK. Integral equation solutions for simply supported polygonal plates. Int. J. Engng. Sci. 1994;12:793–806. [9] Hansen EB. Numerical solution of integro-differential and singular integral equations for plate bending problems. J. Elasticity 1976;6:39–56. [10] Hartmann F. Static analysis of plates. In: Beskos DE, editor. Boundary element analysis of plates and shells, Berlin: Springer, 1991. pp. 1. [11] Hadamard J. Lectures on Cauchy’s problem in linear partial differential equations. New Haven: Yale University Press, 1923. [12] Be´zine G. A boundary integral equation method for plate flexure with conditions inside the domain. Int. J. Numer. Methods Eng. 1981;17:1647–1657. [13] Lekhnitskii SG. Anisotropic plates. London: Gordon and Breach, 1968. [14] Abaqus, Finite element program developed by Hibbitt, Karlson and Sorenson, Inc., USA.