Modelization and numerical approximation of piezoelectric thin shells

Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073 www.elsevier.com/locate/cma Modelization and numerical approximation of piezoelectric thin sh...

Download PDF

1001KB Sizes 4 Downloads 49 Views

Report

PDF Reader
Full Text

Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073 www.elsevier.com/locate/cma

Modelization and numerical approximation of piezoelectric thin shells Part II: Approximation by ﬁnite element methods and numerical experiments Michel Bernadou *, Christophe Haenel P^ ole Universitaire L eonard de Vinci, 92916 Paris la D efense Cedex, France INRIA, Domaine de Voluceau Rocquencourt, B.P. 105, 78153 Le Chesnay Cedex, France Received 11 April 2002; received in revised form 17 January 2003; accepted 27 February 2003

Abstract A two-dimensional modelization of piezoelectric thin shells has been developed in the ﬁrst part of this work. Three equivalent variational formulations have been considered: • an homogeneous one (with respect to the potential) the bilinear form of which is positive deﬁnite but not symmetric; • the associated nonhomogeneous one which is directly related to the natural boundary conditions on the potential; • a third one, whose bilinear form is now symmetric but no longer positive deﬁnite. In this second part of this work, the approximation of the second formulation by a conforming ﬁnite element method is analyzed. It takes into account the use of numerical integration techniques and gives criteria on the choice of suitable numerical schemes. Moreover, the drawbacks linked to the unsymmetrical character of the associated square matrix of the linear system are circumvent by using a condensation method which consists to eliminate the potential and then to only solve a symmetrical and elliptic linear system with respect to the displacement. This method is most eﬀective than a full implementation of the third variational formulation. Finally some numerical experiments on plate and thin cylindrical shells prove the eﬀectiveness of this method. 2003 Elsevier B.V. All rights reserved.

1. Introduction In the ﬁrst part of this work [1], we have constructed a two-dimensional modelization of piezoelectric thin shells and we have considered three equivalent variational formulations:

*

Corresponding author. Address: P^ ole Universitaire Leonard de Vinci, 92916 Paris la Defense Cedex, France. Tel.: +33-1-41167170; fax: +33-1-41167171. E-mail address: [email protected] (M. Bernadou). 0045-7825/$ - see front matter 2003 Elsevier B.V. All rights reserved. doi:10.1016/S0045-7825(03)00362-1

4046

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

ii(i) an homogeneous one with respect to the potential (see Eqs. (I,(136))) the bilinear form of which is positive deﬁnite but not symmetric; this formulation is convenient to prove an existence and uniqueness theorem; i(ii) the associated nonhomogeneous one (see Eq. (I,(148))) which is directly connected to the natural boundary conditions on the potential; (iii) a third one (see Eq. (I,(165))), the bilinear form of which is now symmetrical but no longer positive deﬁnite. This formulation could be interesting for a direct approximation by ﬁnite element methods since it leads to a linear system whose matrix is symmetrical and invertible, even if no positive deﬁnite. For the subsequent developments, we retain the variational formulation (ii) and we analyze its approximation by a conforming ﬁnite element method. The associated essential boundary conditions are easy to impose since they are zero for the displacement and constant for the potential along the appropriate parts of the boundary. Such formulation and approximation are convenient because they allow to adapt directly the approximation studies thoroughly developed in [2,3]. For the implementation, at the ﬁrst glance, the most attractive way seems to be to associate with the second discrete formulation of type (ii), its third associated formulation of type (iii); thus, we would obtain a linear system whose matrix is symmetrical and invertible. In fact, it is most eﬀective to use the second formulation (ii) and to make a condensation upon the potential. Thus, we keep as single unknown the approximate displacement and the matrix of the associated linear system is symmetrical, positive and deﬁnite. Finally, some numerical experiments on plate and cylindrical shells, with various boundary conditions, illustrate the eﬀectiveness of such an approximation. Thus, the contents of this second part of the paper can be summarized as follows. We start by giving a most convenient matrix form of the continuous second variational formulation (ii). Next, we approximate the solution ð~ u; ðu0 ; u1 ; u2 ÞÞ by a conforming ﬁnite element method which combines the Ganev and Dimitrov [4] and Argyris et al. [5] elements and we prove an a priori asymptotic error estimate theorem. For the solution of the ﬁnal matrix system, we propose a condensation method for the potential degrees of freedom in order to only retain the mechanical degrees of freedom. Next, we compare this method with the solution of the full system. Finally, some numerical experiments prove the eﬀectiveness of this method. 1.1. Notations We use the same notations than in Part I and references to formulas, theorems, ﬁgures of this Part I are made by (I,. . .).

2. A matrix formulation for the continuous problem For the analysis of the approximation by a conforming ﬁnite element method and for the implementation, it is convenient to consider a matrix formulation of the continuous problem. As explained in the introduction, we use here the two-dimensional nonhomogeneous unsymmetrical variational formulation (I,(148)). 2.1. The ﬁrst hand member Let us set t

U ¼ ½u1 u1;1 u1;2 u2 u2;1 u2;2 u3 u3;1 u3;2 u3;11 u3;12 u3;22 ;

t

U ¼ ½u0;1 u0;2 u1 u1;1 u1;2 u2 u2;1 u2;2 ð1Þ

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4047

and, in a similar way, t

V ¼ ½v1 v1;1 v1;2 v2 v2;1 v2;2 v3 v3;1 v3;2 v3;11 v3;12 v3;22 ;

t

W ¼ ½w0;1 w0;2 w1 w1;1 w1;2 w2 w2;1 w2;2 :

ð2Þ

Theorem 1. The bilinear form a1 ½; deﬁned by relations (I,(149,150)) can be rewritten Z V u; ðu0 ; u1 ; u2 ÞÞ; ð~ v; ðw0 ; w1 ; w2 ÞÞ ¼ ðt U t UÞ½AIJ dn1 dn2 ; a1 ½ð~ W x with

½AIJ ¼

½A

½B

ð3Þ

ð4Þ

2 M2020 ; t ½B ½C where matrices A, B and C are explicited by the subsequent relations (21). Proof. Let us remind that a1 ½ð~ u; ðu0 ; u1 ; u2 ÞÞ; ð~ v; ðw0 ; w1 ; w2 ÞÞ ¼

R x

9 > =

fN kl ckl ð~ vÞ þ M kl qkl ð~ vÞ

pﬃﬃﬃ S0k w0;k S1k w1;k S2k w2;k þ T1 w1 þ T2 w2 g a dn1 dn2 > ; 8ð~ v; ðw0 ; w1 ; w2 ÞÞ 2 V~KL ðxÞ VP ðxÞ;

ð5Þ

0

with

9 e2 abkl > > N kl ¼ e C cab ð~ uÞ þ p eakl u0;a þ u2;a þ p e3kl u1 ; > > 24 > > > > 3 > e > abkl akl 3kl > kl > C qab ð~ uÞ þ p e u1;a þ p e u2 ; M ¼ > > 12 > > > = 2 3 e e kab kb 3k k kab kb 3k k uÞ d u0;b þ u2;b d u1 ; S1 ¼ ð~ u Þ d u d u e q ; S0 ¼ e p e cab ð~ p ab 1;b 2 > 24 12 > > > > 3 2 > e 3e > kab kb 3k k > u c ð~ u Þ d u þ u e S2 ¼ d ; > p ab 0;b 2;b 1 > 24 40 > > > > 2 3 > > e e 3ab 3b 33 > 3ab 3b 33 ; uÞ d u0;b þ u2;b d u1 ; T2 ¼ ð~ u Þ d u d u e q T1 ¼ e p e cab ð~ p ab 1;b 2 24 12 ð6Þ

and 1 cab ð~ uÞ ¼ ðua;b þ ub;a Þ Cqab uq bab u3 ; 2

ð7Þ

qab ð~ uÞ ¼ ððbkajb bqa Ckqb bqb Ckqa Þuk þ bka uk;b þ bkb uk;a bka bkb u3 Ckab u3;k þ u3;ab Þ:

ð8Þ

Thus, if we set 1 1 1 1 2 2 1 1 2 2 1 2 2 1 2 1 t 1 Tab ¼ Cab ; da db ; ðda db þ da db Þ; Cab ; ðda db þ da db Þ; da db ; bab ; 0; 0; 0; 0; 0 2 2 and t

h 2 Tab ¼ b1ajb bka C1kb bkb C1ka ; d1b b1a þ d1a b1b ; d2b b1a þ d2a b1b ; b2ajb bka C2kb bkb C2ka ;

9 =

i d1b b2a þ d1a b2b ; d2b b2a þ d2a b2b ; bka bkb ; C1ab ; C2ab ; d1a d1b ; d1a d2b þ d2a d1b ; d2a d2b ; ;

ð9Þ

ð10Þ

4048

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

we can check that 1 1 uÞ ¼ t Tab U ¼ t U Tab ; cab ð~

ð11Þ

2 2 uÞ ¼ t Tab U ¼ t U Tab : qab ð~

ð12Þ

Next, we set t

R0a ¼ ½d1a ; d2a ; 0; 0; 0; 0; 0; 0;

t

R1a ¼ ½0; 0; 0; d1a ; d2a ; 0; 0; 0;

t

R2a ¼ ½0; 0; 0; 0; 0; 0; d1a ; d2a ;

t

R4 ¼ ½0; 0; 1; 0; 0; 0; 0; 0;

t

R5 ¼ ½0; 0; 0; 0; 0; 1; 0; 0;

so that we obtain u1 ¼ t R4 U ¼ t U R4 ; u2 ¼ t R5 U ¼ t U R5 ; u0;a ¼ t R0a U ¼ t U R0a ; u1;a ¼ t R1a U ¼ t U R1a ; u2;a ¼ t R2a U ¼ t U R2a : Then, relations (6) involve 2 N kl

M

kl

ð13Þ

3 1 Tab 6 7 ¼ e½t U t U4 5; e2 p eakl R0a þ R2a þ p e3kl R4 24

ð14Þ

" # abkl 2 e3 t t C Tab ¼ ½ U U ; 12 p eakl R1a þ p e3kl R5

ð15Þ

2 6 S0k ¼ e½t U t U4

C

abkl

1 ekab Tab p

3

7 5; e2 2 kb 0 3k 4 Rb þ Rb d R d 24

ð16Þ

2 3 kab 2 3 T e e p ab S1k ¼ ½t U t U4 kb 1 3k 5 5; 12 d Rb d R 2 3

S2k ¼

1 ekab Tab p

ð17Þ 3

e t t 6 7 ½ U U4 kb 0 3e2 2 5; 3k 4 24 R d R Rb þ d 40 b 2

6 T 1 ¼ e½t U t U4

1 e3ab Tab p

ð18Þ

3

7 5 e2 2 3b 0 Rb þ Rb d33 R4 d 24

ð19Þ

and 2 3 3ab 2 3 T e e p ab T 2 ¼ ½t U t U4 3b 1 33 5 5: 12 d Rb d R

ð20Þ

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4049

Then, the substitution of relations (7), (8), (13)–(20) into the expression (5) gives relations (3) and (4) with 9 pﬃﬃﬃ abkl 1 t 1 e2 2 t 2 > > ½A ¼ e aC Tab Tkl þ Tab Tkl ; > > 12 > > > > 2 2 > pﬃﬃﬃ akl 1 t 0 e2 2 e 2 t 1 e 2 t 5 > 3kl 1 t 4 Tkl Ra þ Ra þ Tkl Ra p e Tkl R þ Tkl R ½B ¼ e a p e ;> > > > 24 12 12 > > = 2 2 2 2 pﬃﬃﬃ ab 0 t 0 e 2 e 1 t 1 e 2 t 0 3e 2 ð21Þ Rb Ra Rb þ Rb þ Ra Rb þ Ra Rb þ ½C ¼ e a d > 24 12 24 40 > > > > e2 e2 > > > d3b R4 t R0b þ R2b þ R5 t R1b > > 24 12 > > > 2 2 2 > > e e e > 2 1 t 5 5 t 5 a3 0 t 4 33 4 t 4 ; R a þ Ra R þ R a R þ d R R þ R R d : 24 12 12 2.2. The second hand member Theorem 2. The linear form fM ðÞ deﬁned by relations (I,(110)) can be written Z Z t t fM ð~ vÞ ¼ F V dn1 dn2 þ G V dc; with

ð22Þ

cM 1

x

pﬃﬃﬃ F ¼ paﬃﬃﬃ½P1 Qa b1a ; 0; 0; P2 Qa b2a ; 0; 0; P3 ; Q1 ; Q2 ; 0; 0; 0; t G ¼ a½N 1 M a b1a ; 0; 0; N 2 M a b2a ; 0; 0; N 3 ; M 1 ; M 2 ; 0; 0; 0: t

ð23Þ

Proof. If suﬃces to substitute relation (2) into the expression (I,(110)), i.e., Z pﬃﬃﬃ vÞ ¼ fðPk Qa bka Þvk þ P3 v3 Qa v3;a g a dn1 dn2 fM ð~ xZ þ fðN k M a bka Þvk þ N 3 v3 M a v3;a g dc: cM 1

3. Approximation by a conforming ﬁnite element method 3.1. The ﬁnite element spaces Firstly, let us deﬁne suitable ﬁnite element subsets of the admissible sets V~KL ðxÞ and VPhd respectively deﬁned by (I,(114) and (93)), i.e., ov3 2 1 2 M ~ ¼ 0 on c0 ; VKL ðxÞ ¼ ~ v 2 ðH ðxÞÞ H ðxÞ; vi ¼ om ( VPhd ¼ ðH 1 ðxÞÞ3 \

ðh0 ; h1 ; h2 Þ such that h0 ¼ hd0 and h1 ¼ h2 ¼ 0 on cE0 ;

hþ hþ e2 d þ hd 1 d hd 1 h2 ¼ 1 and h1 ¼ 1 in x1 ; 8 2 e ) 2 2 e e e e h0 þ h1 þ h2 ¼ hþ h2 ¼ h d2 in x2 ; h0 h1 þ d3 in x3 : 2 2 8 8

h0 þ

4050

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

In this way, for simplicity, we assume that x is a polygonal domain which can be exactly covered by a regular family Th of triangulations in the following sense, i(i) there exists a constant r > 0 such that qhKK 6 r 8K 2 Th , where, for each triangle K, we denote hK ¼ diamðKÞ and qK ¼ supfdiamðSÞÞ; S is a ball contained in Kg; (ii) the quantity h ¼ maxK2Th hK ! 0. With each triangulation Th we associate the ﬁnite element spaces Xh1 and Xh2 constructed from the Ganev and Dimitrov [4] and Argyris et al. [5] ﬁnite elements (see Figs. 1 and 2). In order to take easily into account the boundary conditions and potential data, we assume in addition E ii(i) that parts cM 0 and c0 of boundary c ¼ ox are the exact unions of triangle sides; i(ii) that sub-domains x1 , x2 and x3 of x (see Fig. 3) are the exact unions of triangles K 2 Th ; (iii) that thickness eðn1 ; n2 Þ is constant upon x; this hypothesis is not restrictive since piezoelectric devices have generally constant thickness; (iv) that applied potentials ud are constant upon each connected component of cE0 , x1 f 2eg, x2 fþ 2eg and x3 f 2eg.

Then, we deﬁne the ﬁnite element space V~KLh ðxÞ ¼ Vh1 ðxÞ Vh1 ðxÞ Vh2 ðxÞ, with Vh1 ðxÞ ¼ fvh 2 Xh1 ; M vh ¼ 0 on cM 0 g and Vh2 ðxÞ ¼ fvh 2 Xh2 ; vh ¼ om vh ¼ 0 on c0 g.

Fig. 1. Ganev and DimitrovÕs triangle [4].

Fig. 2. Argyris et al.Õs triangle [5].

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4051

Fig. 3. A generic distribution of electrodes (displayed on the reference domain).

Likewise, for h ¼ 0 or h ¼ ud0 (we denote h the value of h upon C ), we deﬁne 3

VPhh ¼ fðw0h ; w1h ; w2h Þ 2 ðXh1 Þ such that w0h ¼ h and w1h ¼ w2h ¼ 0 on cE0 ; w0h þ

e2 h þ þ h hþ h w2h ¼ ¼ cst: and w1h ¼ ¼ cst: in x1 ; 8 2 e

e e2 w0h þ w1h þ w2h ¼ hþ ¼ cst: in x2 ; 2 8

9 > > > > > > =

> > > > > e e2 ; w0h w1h þ w2h ¼ h ¼ cst: in x3 g: > 2 8

ð24Þ

Due to previous assumptions (i) to (iv), we obtain the following inclusions: ~KLh ðxÞ V~KL ðxÞ; V

VPud h ðxÞ VPud ðxÞ;

VP0h ðxÞ VP0 ðxÞ;

which lead to a conforming ﬁnite element approximation of the problem (I,(148)). For completeness, let us mention how it is possible to realize the constraints which appear in deﬁnition (24). Since the set VPhh is constructed from the Ganev and DimitrovÕs ﬁnite element, we are led to impose the following constraints to the degrees of freedom of triangle K: • if K x1 (note that hþ and h are constant upon x1 ), for i ¼ 1; 2; 3: e2 e2 hþ þ h w2h ðai Þ ¼ w0h ðbi Þ þ w2h ðbi Þ ¼ ; 8 8 2 e2 e2 D w0h ðai Þ þ w2h ðai Þ ðaj ai Þ ¼ 0; j 6¼ i; D w0h ðbi Þ þ w2h ðbi Þ ðci bi Þ ¼ 0; 8 8 w0h ðai Þ þ

w1h ðai Þ ¼ w1h ðbi Þ ¼

hþ h ; e

Dðw1h ðbi ÞÞðci bi Þ ¼ 0;

Dðw1h ðai ÞÞðaj ai Þ ¼ 0; j 6¼ i;

• if K x2 (note that hþ is constant upon x2 ), for i ¼ 1; 2; 3: e e2 w0h ðai Þ þ w1h ðai Þ þ w2h ðai Þ ¼ hþ ; 2 8

e e2 w0h ðbi Þ þ w1h ðbi Þ þ w2h ðbi Þ ¼ hþ ; 2 8

4052

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

e e2 D w0h ðai Þ þ w1h ðai Þ þ w2h ðai Þ ðaj ai Þ ¼ 0; 2 8 e e2 D w0h ðbi Þ þ w1h ðbi Þ þ w2h ðbi Þ ðci bi Þ ¼ 0; 2 8

j 6¼ i;

• if K x3 (note that h is constant upon x3 ), for i ¼ 1; 2; 3: e e2 e e2 w0h ðai Þ w1h ðai Þ þ w2h ðai Þ ¼ h ; w0h ðbi Þ w1h ðbi Þ þ w2h ðbi Þ ¼ h ; 2 2 8 8 e e2 D w0h ðai Þ w1h ðai Þ þ w2h ðai Þ ðaj ai Þ ¼ 0; j 6¼ i; 2 8 e e2 D w0h ðbi Þ w1h ðbi Þ þ w2h ðbi Þ ðci bi Þ ¼ 0: 2 8 For simplicity, let us start with the deﬁnition and properties of the associated discrete problems without any use of numerical integration. 3.2. The ﬁrst discrete problem The ﬁrst discrete problem is the conforming ﬁnite element approximation of the nonhomogeneous problem (I,(148)) associated with the ﬁnite element spaces deﬁned in Section 3.1. 9 Find ð~ u ; ðu ; u ; u ÞÞ 2 V~KLh ðxÞ VPud h ðxÞ such that > h 0h 1h 2h = ~KLh ðxÞ VP ðxÞ: > a1 ½ð~ u ; ðu ; u ; u ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ ¼ fM ð~ vh Þ 8ð~ vh ; ðw0h ; w1h ; w2h ÞÞ 2 V ; 0h h

0h 1h 2h

ð25Þ Then, we obtain the existence and uniqueness of a solution: Theorem 3. The ﬁrst discrete problem (25) has a unique solution ð~ u ; ðu0h ; u1h ; u2h ÞÞ. h

~KL ðxÞ Proof. Firstly, let us observe that V~KLh ðxÞ VP0h ðxÞ is a closed ﬁnite dimensional subspace of V VP0 ðxÞ. Thus, the proof of existence and uniqueness of a solution follows the same lines than the proof of ^ h 2 VPud h ðxÞ such that u ^ h jCE ¼ ud . h (I,Theorem 6): it suﬃces to introduce an extension function u 0

Moreover, we obtain the following a priori asymptotic error estimate theorem: Theorem 4. We assume that the solution ð~ u; ðu0 ; u1 ; u2 ÞÞ of the continuous problem (I,(148)) belongs to the 2 3 space ðH 5 ðxÞÞ H 6 ðxÞ ðH 5 ðxÞÞ , that AIJ 2 W 4;1 ðxÞ and that Pi 2 W 4;q ðxÞ, q P 2. Then, there exists a constant C > 0, independent of h, such that k~ u ~ u k~KL ðxÞ þ kðu0 ; u1 ; u2 Þ ðu0h ; u1h ; u2h ÞkVP V h

6 Ch4

8 < X 2 :

a¼1

2

2

kua k5;x þ ku3 k6;x þ

2 X i¼0

0

ðxÞ

!1=2 2

kui k5;x

þ

3 X i¼1

q

kPi k4;q;x

!1=q 9 = ;

;

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4053

where ð~ uh ; ðu0h ; u1h ; u2h ÞÞ is the solution of the discrete problem (25) and where k kV~KL ðxÞ and

kðu0 ; u1 ; u2 ÞkVP

ðxÞ 0

denote the usual cross-product norms on ðH 1 ðxÞÞ2 H 2 ðxÞ and on ðH 1 ðxÞÞ3 .

^ h jCE to associate with problem (25) an homogeneous one Proof. Once again, we use the extension function u 0 and then we conclude with the usual properties of conforming ﬁnite element approximations [2,3]. h 3.3. The second discrete problem It is generally impossible to exactly compute the integrals which appear in statement (25). Thus, we are led integration techniques, i.e., over each triangle K 2 Th , we introduce an approximation R to use numerical PL hðnÞ dn ’ ‘¼1 x‘;K hðb‘;K Þ, where x‘;K and b‘;K denote respectively the weights and the nodes of the K . In this way, we obtain numerical integration scheme while h denotes a continuous function deﬁned on K the statement of the second discrete problem: 9 Find ð~ uh ; ðu0h ; u1h ; u2h ÞÞ 2 V~KLh ðxÞ VPud h ðxÞ such that > > > = ð26Þ a1h ½ð~ uh ; ðu0h ; u1h ; u2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ ¼ fMh ð~ vh Þ > > > ; 8ð~ vh ; ðw0h ; w1h ; w2h ÞÞ 2 V~KLh ðxÞ VP0h ðxÞ; with, according to relations (3) and (22): a1h ½ð~ uh ; ðu0h ; u1h ; u2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ ¼

P

K2Th

PL

t ‘¼1 x‘;K ðUh

Uh Þðb‘;K Þ½AIJ ðb‘;K Þ

Vh Wh

!

9 > > = ðb‘;K Þ > > ;

ð27Þ

and fMh ð~ vh Þ ¼

L X X

x‘;K t F ðb‘;K ÞVh ðb‘;K Þ:

K2Th ‘¼1

For simplicity, we assume that we do not have any mechanical loading upon cM 1 ; if we have some, we have just to use an another numerical integration scheme upon the concerned sides of the triangles. Of course, to be able to compute the matrices ½AIJ ðb‘;K Þ and t F ðb‘;K Þ, we need to assume that the diﬀerent coeﬃcients are continuous; a convenient way for that (but not strictly mandatory) is to assume that ~S 2 ðC3 ðx ÞÞ3 . U Now, we can extend to the approximation of piezoelectric thin shells the results obtained in [2] for thin shells made of classical materials. The framework of this study is given by the following abstract error estimate theorem: Theorem 5 (Abstract error estimate). Assume that the bilinear form a1h ½; deﬁned by (27) is uniformly elliptic on the discrete ﬁnite element space V~KLh ðxÞ VP0h ðxÞ, e.g., there exist constants a > 0 and h1 > 0, independent of h, such that, for any h < h1 , we have ) aðk~ vh k2V~KLðxÞ þ kðw0h ; w1h ; w2h Þk2VP ðxÞ Þ 0 ð28Þ 6 a1h ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ 8ð~ vh ; ðw0h ; w1h ; w2h ÞÞ 2 V~KLh ðxÞ VP0h ðxÞ:

4054

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Then, there exist a constant C, independent of h, such that: 9 k~ u ~ uh kV~KL ðxÞ þ kðu0 ; u1 ; u2 Þ ðu0h ; u1h ; u2h ÞkVP ðxÞ > 0 > ( > > n > > > 6 C inf ð~vh ;ðw0h ;w1h ;w2h ÞÞ2V~KL ðxÞVP ðxÞ k~ u ~ vh kV~KL ðxÞ : þ kðu0 ; u1 ; u2 Þ ðw0h ; w1h ; w2h ÞkVP ðxÞ > > > h ud h 0 = o > þ supð~wh ;ðh0h ;h1h ;h2h ÞÞ2V~KL ðxÞVP ðxÞ Nh ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ wh ; ðh0h ; h1h ; h2h ÞÞ > > h 0h > ) > > > jfM ð~ wh Þ fMh ð~ wh Þj > > > þ sup~wh 2V~KL ðxÞ ; ; h k~ wk

ð29Þ

~KL ðxÞ h V

where ð~ u; ðu0 ; u1 ; u2 ÞÞ and ð~ uh ; ðu0h ; u1h ; u2h ÞÞ are the solutions of the continuous problem (I,(148)) and of the second discrete problem (26), and where vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ wh ; ðh0h ; h1h ; h2h ÞÞ Nh ½ð~ ¼

ja1 ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ wh ; ðh0h ; h1h ; h2h ÞÞ a1h ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ wh ; ðh0h ; h1h ; h2h ÞÞj : k~ wh kV~KL ðxÞ þ kðh0h ; h1h ; h2h ÞkVP ðxÞ 0

Proof. Firstly, let us quote that the V~KLh ðxÞ VP0h ðxÞ-ellipticity (28) of a1h ½; insures the existence and uniqueness property for the solution ð~ uh ; ðu0h ; u1h ; u2h ÞÞ of problem (26). For any ðw0h ; w1h ; w2h Þ 2 VPud h ðxÞ, the diﬀerence ðu0h ; u1h ; u2h Þ ðw0h ; w1h ; w2h Þ 2 VP0h ðxÞ. Moreover, relation (28) involves: 2

2

aðk~ uh ~ vh kV~KL ðxÞ þ kðu0h ; u1h ; u2h Þ ðw0h ; w1h ; w2h ÞkVP

0

ðxÞ Þ

6 a1h ½ð~ uh ; ðu0h ; u1h ; u2h ÞÞ ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ uh ; ðu0h ; u1h ; u2h ÞÞ ð~ vh ; ðw0h ; w1h ; w2h ÞÞ: This inequality can be developed as follows: 2

2

aðk~ uh ~ vh kV~KL ðxÞ þ kðu0h ; u1h ; u2h Þ ðw0h ; w1h ; w2h ÞkVP

0

ðxÞ Þ

6 a1 ½ð~ u; ðu0 ; u1 ; u2 ÞÞ ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ uh ; ðu0h ; u1h ; u2h ÞÞ ð~ vh ; ðw0h ; w1h ; w2h ÞÞ þ a1 ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ uh ; ðu0h ; u1h ; u2h ÞÞ ð~ vh ; ðw0h ; w1h ; w2h ÞÞ a1h ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ uh ; ðu0h ; u1h ; u2h ÞÞ ð~ vh ; ðw0h ; w1h ; w2h ÞÞ a1 ½ð~ u; ðu0 ; u1 ; u2 ÞÞ; ð~ uh ; ðu0h ; u1h ; u2h ÞÞ ð~ vh ; ðw0h ; w1h ; w2h ÞÞ þ a1h ½ð~ uh ; ðu0h ; u1h ; u2h ÞÞ; ð~ uh ; ðu0h ; u1h ; u2h ÞÞ ð~ vh ; ðw0h ; w1h ; w2h ÞÞ: u; ðu0 ; u1 ; u2 ÞÞ and ð~ uh ; ðu0h ; u1h , u2h ÞÞ are the Since the bilinear form a1 ½; is continuous and since ð~ respective solutions of the problems (I,(148)) and (26), there exists a constant M > 0 such that 2

2

vh kV~KL ðxÞ þ kðu0h ; u1h ; u2h Þ ðw0h ; w1h ; w2h ÞkVP aðk~ uh ~

0

ðxÞ Þ

6 Mðk~ u ~ vh kV~KL ðxÞ þ kðu0 ; u1 ; u2 Þ ðw0h ; w1h ; w2h ÞkVP þ kðu0h ; u1h ; u2h Þ ðw0h ; w1h ; w2h ÞkVP

0

ðxÞ Þ

0

uh ðxÞ Þðk~

~ vh kV~KL ðxÞ

þ ja1 ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ uh ; ðu0h ; u1h ; u2h ÞÞ

ð~ vh ; ðw0h ; w1h ; w2h ÞÞ a1h ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ uh ; ðu0h ; u1h ; u2h ÞÞ ð~ vh ; ðw0h ; w1h ; w2h ÞÞj þ jfM ð~ uh ~ vh Þ fMh ð~ uh ~ vh Þj:

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4055

Hence, by using the triangular inequality k~ u ~ uh kV~KL ðxÞ þ kðu0 ; u1 ; u2 Þ ðu0h ; u1h ; u2h ÞkVP

0

ðxÞ

6 k~ u ~ vh kV~KL ðxÞ þ kðu0 ; u1 ; u2 Þ ðw0h ; w1h ; w2h ÞkVP þ kðw0h ; w1h ; w2h Þ ðu0h ; u1h ; u2h ÞkVP

0

0

ðxÞ

þ k~ vh ~ uh kV~KL ðxÞ

ðxÞ

and, by taking the minimum with respect to ð~ vh ; ðw0h ; w1h ; w2h ÞÞ 2 V~KLh ðxÞ VPud h ðxÞ, we obtain the inequality (29). h Now we give suﬃcient conditions to insure the property of uniform ellipticity for a1 ½; . Þ, 1 6 I, Theorem 6. Let us assume that the coeﬃcients which appear in relation (4) satisfy AIJ 2 W 1;1 ðx J 6 20 and that the numerical integration scheme satisﬁes: ðiÞ the integration nodes b‘;K 2 K

8K 2 Th ; ‘ ¼ 1; . . . ; L;

ð30Þ

ðiiÞ the scheme is exact for polynomials of degree 6; i:e:; R numericalPintegration L P ðxÞ dx x P ðb 8P 2 P6 ðKÞ 8K 2 Th : ‘;K ‘;K Þ ¼ 0 ‘¼1 K

ð31Þ

Then, the bilinear form a1h ½; is uniformly elliptic on V~KLh ðxÞ VP0h ðxÞ, e.g., the inequality (28) is satisﬁed. Proof. Since V~KLh ðxÞ VP0h ðxÞ V~KL ðxÞ VP0 ðxÞ, we obtain for any ð~ vh ; ðw0h ; w1h , w2h ÞÞ 2 V~KLh ðxÞ VP0h ðxÞ: a1h ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ ¼ a1 ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ þ a1h ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ: a1 ½ð~ Moreover, the bilinear form a1 ½; is V~KL ðxÞ VP0 ðxÞ-elliptic: there exists a constant b > 0, independent of h, such that (I,(146)): 2

2

a1 ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ P bðk~ vh kV~KL ðxÞ þ kðw0h ; w1h ; w2h ÞkVP

0

ðxÞ Þ

ð32Þ

for any ð~ vh ; ðw0h ; w1h ; w2h ÞÞ 2 V~KLh ðxÞ VP0h ðxÞ. Then, a proof entirely similar to that of [2, Theorem 1.3.4] shows that, under the assumptions (30) and (31) on the numerical integration scheme, there exists a constant C > 0, independent of h, such that 9 > ja1 ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ a1h ½ð~ vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞj = ð33Þ 2 2 6 Chðk~ vh kV~KL ðxÞ þ kðw0h ; w1h ; w2h ÞkVP ðxÞ Þ > ; 0 ~KLh ðxÞ VP ðxÞ. for any ð~ vh ; ðw0h ; w1h ; w2h ÞÞ 2 V 0h By combining the inequalities (32) and (33), we ﬁnally obtain for any ð~ vh ; ðw0h ; w1h ; w2h ÞÞ 2 V~KLh ðxÞ V~P0h ðxÞ 2

2

vh ; ðw0h ; w1h ; w2h ÞÞ; ð~ vh ; ðw0h ; w1h ; w2h ÞÞ P ðb ChÞðk~ vh kV~KL ðxÞ þ kðw0h ; w1h ; w2h ÞkVP a1h ½ð~ b . so that we obtain (28) with, for instance, a ¼ b2 and h1 ¼ 2C

0

ðxÞ Þ;

4056

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Fig. 4. A numerical integration scheme exact for polynomials of degree 6 [6,7].

In Fig. 4, we recall a numerical integration scheme which satisﬁes property (31).

h

Theorem 7 (A priori asymptotic error estimate). Let us assume that 2

• the solution ð~ u; ðu0 ; u1 ; u2 ÞÞ of the continuous problem (I,(148)) belongs to the space ðH 5 ðxÞÞ 3 6 5 H ðxÞ ðH ðxÞÞ , AIJ 2 W 4;1 ðxÞ and Pi 2 W 4;q ðxÞ, q P 2; • the integration numerical scheme satisﬁes the hypotheses (30) and (31). Then, there exist constants C > 0 and h1 > 0, independent of h, such that, for any h < h1 , 9 k~ u ~ uh kV~KL ðxÞ þ kðu0 ; u1 ; u2 Þ ðu0h ; u1h ; u2h ÞkVP ðxÞ > = 0 1=2 1=q ð34Þ P P P 2 2 2 2 2 3 i q 6 Ch4 þ ;> ; a¼1 kua k5;x þ ku3 k6;x þ i¼0 kui k5;x i¼1 kP k4;q;x where ð~ uh ; ðu0h ; u1h ; u2h ÞÞ denotes the solution of problem (26).

! Proof. The above assumptions allow us to apply the abstract error estimate (29). By denoting ph u the ~KLh ðxÞ-interpolent of ~ V u and ph ðu0 ; u1 ; u2 Þ the VPud h ðxÞ-interpolent of ðu0 ; u1 ; u2 Þ, estimate (29) gives 9 k~ u ~ uh kV~KL ðxÞ þ kðu0 ; u1 ; u2 Þ ðu0h ; u1h ; u2h ÞkVP ðxÞ > 0 > > ( > > > > ! > > 6 C k~ u ph ukV~KL ðxÞ þ kðu0 ; u1 ; u2 Þ ph ðu0 ; u1 ; u2 ÞkVP ðxÞ > > 0 = ð35Þ ! þ supð~wh ;ðh0h ;h1h ;h2h ÞÞ2V~KL ðxÞVP ðxÞ Nh ½ðph u; ph ðu0 ; u1 ; u2 ÞÞ; ð~ wh ; ðh0h ; h1h ; h2h ÞÞ > > > h 0h > ) > > > > jfM ð~ wh Þ fMh ð~ wh Þj > > þ sup~wh 2V~KL ðxÞ : > ; h k~ wh kV~KL ðxÞ From [3, Theorem 3.2.1], we obtain 9

! > k~ u ph u kV~KL ðxÞ þ kðu0 ; u1 ; u2 Þ ph ðu0 ; u1 ; u2 ÞkVP ðxÞ = 0 1=2 P 1=q P P 2 2 2 2 2 3 i q ; 6 Ch4 þ :> a¼1 kua k5;x þ ku3 k6;x þ i¼0 kui k5;x i¼1 kP k4;q;x

ð36Þ

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Next, a proof entirely similar to that of [2, Theorems 1.3.3 and 1.3.5] gives the estimates 9 ! þ supð~wh ;ðh0h ;h1h ;h2h ÞÞ2V~KL ðxÞVP ðxÞ Nh ½ðph u; ph ðu0 ; u1 ; u2 ÞÞ; ð~ wh ; ðh0h ; h1h ; h2h ÞÞ = h 0h P P 1=2 P2 20 2 2 2 2 6 Ch4 ; ; I;J ¼1 kAIJ k4;1;x a¼1 kua k5;x þ ku3 k6;x þ i¼0 kui k5;x

sup ~KL ðxÞ ~ wh 2V h

jfM ð~ wh Þ fMh ð~ wh Þj 6 Ch4 k~ wh kV~KL ðxÞ

3 X

4057

ð37Þ

!1=q q kPi k4;q;x

:

ð38Þ

i¼1

Thus, we obtain estimate (34) by combining estimates (35)–(38). h

4. Solution by condensation 4.1. The linear system into consideration The implementation of the second discrete problem (26) leads, after assemblage and before taking into account any boundary conditions, to a linear system of the type: 2 3 0 1 0 1 0 0 Kuu Kuu Kuu Kuu U F 6 t 0 7 00 0 00 7B 0C 6 Kuu BF0 C Kuu Kuu Kuu U 6 t 7B C ¼ B C ; ð39Þ 0 7@ 6 Kuu t K 0 A @0A Kuu Kuu uu 4 5 U 0 00 t 0 00 t Kuu t Kuu Kuu Kuu U0 0 where superscripts (0 ) and (00 ) denote the degrees of freedom associated with the boundary conditions which appear in the space V~KL ðxÞ VPud h ðxÞ and the sub-matrices of the linear system which are concerned by these degrees of freedom. But these boundary conditions can be expressed as fU 0 g ¼ f0g : clamped conditions along cM 0 ;

fU0 g ¼ fud g : applied potential along cE0 :

ð40Þ

The substitution of conditions (40) into the system (39) leads to the ﬁnal system associated with problem (26): ! 0 fud g F ½Kuu Kuu Kuu U : ð41Þ ¼ 0 ½Kuu fud g t Kuu Kuu U Remark 1. Instead to discretize problem (I,(148)), we could, in the same way, discretize problem (I,(165)). The result is quite the same, except that the sign of the third and fourth lines of (39) is changed and also the sign of the second line of (41). It is easy to check that the condensation process proposed in Section 4.2 leads exactly to the same linear system. In other words, the discretization of problems (I,(148)) and (I,(165)) leads, after discretization, to the same discrete problem. 4.2. The condensation method According to Theorem 6, the matrix of the linear system (41) is positive deﬁnite so that the submatrix Kuu is itself positive deﬁnite. Thus, matrix Kuu is invertible and relation (41) implies 1

0 fUg ¼ ½Kuu ðt ½Kuu fU g ½Kuu fud gÞ:

ð42Þ

4058

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Substituting this expression of fUg into the ﬁrst line of system (41), we ﬁnally obtain that fU g is solution of the system e uu fU g ¼ f Fe g; ½K

ð43Þ

e uu and f Fe g are deﬁned by where ½ K e uu ¼ ½Kuu þ ½Kuu ½Kuu 1 ðt ½Kuu Þ; ½K

0 0 f Fe g ¼ fF g þ ð½Kuu ½Kuu 1 ½Kuu ½Kuu Þfud g:

ð44Þ

e uu is Theorem 8. The second discrete problem (26) is equivalent to solve the linear system (43) whose matrix ½ K symmetric positive deﬁnite. From fU g we can determine the potential fUg by using relation (42). Proof. The equivalence between the second discrete problem (26) and the solution of the linear system (43) results from the above considerations. e uu is a direct consequence of its deﬁnition. The argument which is used Next, the symmetry of matrix ½ K to prove that Kuu is positive deﬁnite is still available for Kuu which is itself symmetric positive deﬁnite. Since 0 ½Kuu 1 is likewise symmetric positive deﬁnite, the last term of (44), i.e., ½Kuu ½Kuu 1 ðt ½Kuu Þ is symmetric and e uu is symmetric positive deﬁnite. h positive. Thus, ½ K 4.2.1. Solution of the linear system (43) e uu shows that a direct solution is not really appropriate. Thus, to solve the The observation of matrix ½ K e uu , we use a conjugate gradient method with presystem (43) combined with the expression (44) of ½ K conditioning. The preconditioning matrix is the matrix ½Kuu factorized in the incomplete Cholesky form. For the solution of the systems involved in our examples, the number of iterations is three or four. A comparison between • the solution of the complete system (without condensation) by the Crout method; • the solution of the condensed system (43) by a conjugate gradient method, with or without preconditioning, is reported in Fig. 5: it includes some curves giving the dependence on the necessary memory allocation and curves giving the dependence on the computing time, both with respect to the number of mechanical de-

Fig. 5. Solution of the second discrete problem. Comparison between diﬀerent solution methods: CS ¼ complete system; Cond. S ¼ condensed system; SCG ¼ simple conjugate gradient; PCG ¼ preconditioned conjugate gradient.

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4059

grees of freedom. The conclusion is that the conjugate gradient with preconditioning is the most eﬀective method with respect to memory allocation and computing time as well. 4.2.2. Implementation The implementation is done by using MODULEF code [8,9]. For a brief description of the new diﬀerent modules the reader is referred to [10, Section 2.6].

5. Some numerical experiments In this paragraph, we analyze some numerical experiments obtained on plate and cylindrical shells by using the 2D-modelization developed in Part I and the second discrete problem deﬁned by relation (26). These examples illustrate the actuator eﬀect since we apply to these piezoelectric thin shells some potentials and we observe the associated deformations. In practice, let us remind that to apply a potential to such a shell, we need to metallize some parts of the shell which serve as electrodes. In our examples, these electrodes are sticked on some parts or on the whole of the upper and lower surfaces of the plate or of the shell. The potential is constant upon each electrode. Since these electrodes are very thin and light, we neglect their mechanical counterparts. 5.1. Some piezoelectric thin plate problems Let us consider a rectangular plate of length L, width ‘ and thickness e, the middle surface of which is presented in Fig. 6. We assume that this plate is made of a piezoelectric material, polarized through the thickness, and whose mechanical, piezoelectric and dielectric coeﬃcients (I,(13,14)) are given by relations (45)–(47) and Table 1.

Fig. 6. The middle surface of the plate into consideration.

Table 1 Mechanical, piezoelectric and dielectric coeﬃcients of the PZT material C m2

GPa

F m1

C11

C12

C13

C33

C44

C66

p e31

p e33

p e15

d11

d33

126

79.5

84.1

117

23

23.25

)6.5

23.3

17.0

1.503

1.3

4060

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

2

6C 6 6 6 6 6 ijk‘ ½C ¼ 6 6 6 6 6 4

1111

C C 2222

e 111 6p 6 e 211 4p 311 pe

2

11

1133

C C 2233 C 3333

1123

C C 2223 C 3323 C 2323

122 pe 222 pe 322 pe

12

d 6 d ½d ij ¼ 6 d 22 4 sym:

133 pe 233 pe 333 pe

123 pe 223 pe 323 pe

3

2 d11 d 7 7 4 d 23 5 ¼ 0 0 d 33

13

131 pe 231 pe 331 pe

3

C 7 2 C11 7 2212 7 C C 7 6 6 12 3312 7 6 C 7 6 C13 7¼6 2312 7 0 C 7 6 4 0 3112 7 C 7 5 0 1212 C

1131

C C 2231 C 3331 C 2331 C 3131

sym:

2 ½p ekij ¼

1122

1112

112 pe 212 pe 312 pe

3

2 0 7 7¼4 0 5 p e31

C12 C11 C13 0 0 0

0 0 e p 31

C13 C13 C33 0 0 0

0 0 e p 33

0 0 0 C44 0 0

0 e p 15 0 0

0 0 0 0 C44 0

p e15

0 0

3 0 0 7 7 0 7 7; 0 7 7 0 5 C66

ð45Þ

3 0 0 5; 0

ð46Þ

3 0 0 5: d33

0 d11 0

ð47Þ

Remark 2. Deﬁnitions (45)–(47) take into account the symmetrical properties of the tensors C ijk‘ , p e kij and d ij , e.g., C ijk‘ ¼ C jik‘ ¼ C k‘ij , p ekij ¼ p e kji and d ij ¼ d ji . Thus, in the most general case, the 81 coeﬃcients C ijk‘ kij are reduced to 21, the 27 coeﬃcients e are reduced to 18 and the 9 coeﬃcients d ij are reduced to 6. To obtain the two-dimensional formulation, we have made the hypothesis of approximatively plane stresses during the deformation. This assumption has allowed us to introduce a new set of coeﬃcients given by relations (I,(56)), i.e.,

C abkl ¼ C abkl C abj3 C i3kl C ji ;

dij ¼ d ij þ C k‘p e ik3 p e j‘3

ekab p

¼ p e kab C abj3 C jip eki3 ;

1

with ðC ij Þi;j¼1;2;3 ¼ ðC k3‘3 Þk;‘¼1;2;3 :

For the type of material into consideration the characterizations of which are given by the relations (45)– (47), we obtain 1 0 1 0 0 C B C44 C B 1 C B ðC ji Þ ¼ B 0 0 C; C B C44 @ 1 A sym: C33 so that 1 1 0 0 ðC13 Þ2 ðC13 Þ2 C11 C12 0 C B C 1111 C 1122 C 1112 C B C33 C33 C B C B C B 2 ½C abkl ¼ B ð48Þ 2222 2212 C ¼ B C; ðC13 Þ C B B C C C 0 C A @ 11 A @ C33 sym: C 1212 sym: C66

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

1

0

B p e 111 B B kab ½p e ¼ B p e 211 B @ 311 e p 0 "

#

dij

e 122 p e 222 p e 322 p

kab

1 0 0C C; A 0

ð49Þ

1 0 d11 þ

ðp e15 Þ C44

0 2

d33 þ

C C C C C: 0 C C 2A ðp e33 Þ

ð50Þ

C33

and dij are independent of n so that we also have, according to (I,(68))

ekab p

0 0 C13 p e33 p e31 C33

0

ðp e15 Þ2 d þ B 11 C44 d13 C B B C B C d23 C ¼ B B A B 33 @ d sym:

¼ C abkl ;

0 0 C C B C B 0 C¼@ C13 C A p e33 p e31 C33

1

These values of C abkl , p e abkl

e 112 p e 212 p e 312 p

B d11 d12 B ¼B d22 B @ sym:

C

4061

3

dij ¼ dij :

¼ p e kab ;

ð51Þ

~S : ðn1 ; n2 Þ 2 x 7! U ~S ðn1 ; n2 Þ 2 S is reduced to the identity. Thus, Moreover, for this plate, the mapping U aab ¼ dab ;

bab ¼ 0;

Ckab ¼ 0:

ð52Þ

Then, the substitution of relations (48)–(52) into relations (6)–(8) gives:

9 e3 abkl > 3kl > N ¼ e½C cab ð~ uÞ þ p e u1 ; M ¼ ½C qab ð~ uÞ þ p e u2 ; > > > 12 > = 2 3 3 2 e e e 3e kb k kb k kb k S0 ¼ ed u2;b ; u0;b þ u2;b ; S1 ¼ d u1;b ; S2 ¼ d u0;b þ > 24 12 24 40 > > > 3 > e > 3ab 33 3ab 33 ; uÞ d u1 ; T2 ¼ ½p e qab ð~ uÞ d u2 ; T1 ¼ e½p e cab ð~ 12

ð53Þ

cab ð~ uÞ ¼ 12 ðua;b þ ub;a Þ;

ð54Þ

kl

abkl

3kl

kl

with qab ð~ uÞ ¼ u3;ab :

From now on, by using these data, we consider four examples of piezoelectric plates associated with four diﬀerent sets of loadings and boundary conditions. Example 1. Clamped plate with applied potential on the whole upper and lower surfaces. 1 L The plate described in Fig. 6 is supposed to be clamped along the side n2 ¼ 0, i.e., cM 0 ¼ ðn ; 0Þ; 2 6 1 L þ n 6 2g. The upper and lower faces of the plate are entirely metallized so that we can apply a potential V on the upper face and 0 on the lower face. According to the partition of CE0 discussed in (I,Section 3.2, in particular Fig. 3), we have CE0 1 ¼ x f 2eg and x1 ¼ x. These electrical boundary conditions involve the following relations upon the three terms of the development of the potential through the thickness (see (I,(92)): u0 þ

e2 Vþ u2 ¼ 8 2

and

u1 ¼

Vþ e

in x1 ¼ x:

2

Thus, u0;a ¼ e8 u2;a and u1;a ¼ 0 in x, so that for this example, we obtain from (53): S0k ¼

e3 kb d u2;b ; 12

S1k ¼ 0;

S2k ¼

e5 kb d u2;b : 480

ð55Þ

4062

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Of course, according to (I,(94)) the data (55) are essential constraints for the deﬁnition of space VP0 ðxÞ so that e2 ð56Þ w0 þ w2 ¼ 0 and w1 ¼ 0 in x 8ðw0 ; w1 ; w2 Þ 2 VP0 ðxÞ: 8 2 As a consequence, w0;a ¼ e8 w2;a and w1;a ¼ 0, so that the two-dimensional variational equation (I,(151)) can be splitted in two independent equations with unknowns ua for the ﬁrst one (note that u1 is þ equal to Ve according to (55)) and, u3 and u2 for the second one (note that we do not apply any mechanical loading and that, according to (56), we have used w1 ¼ 0): ) 1 M Find ua 2 H 1 ðxÞ; ua ¼ 0 on cM 0 such that 8va 2 H ðxÞ; va ¼ 0 on c0 ; R R ð57Þ pﬃﬃﬃ pﬃﬃﬃ abkl eC cab ð~ uÞckl ð~ vÞ a dn1 dn2 ¼ x p e3kl V þ ckl ð~ vÞ a dn1 dn2 : x 9 1 Find u3 2 H 2 ðxÞ; u3 ¼ @m u3 ¼ 0 on cM 0 and u2 2 H ðxÞ such that > > Z 3 > > > e e3 abkl 3kl 3ab 33 > ðC qab ð~ uÞ þ p e u2 Þ qkl ð~ vÞ ðp e qab ð~ uÞ d u2 Þw2 > = 12 12 x ð58Þ pﬃﬃﬃ 1 2 > e5 kb > > þ a dn dn ¼ 0 d u2;b w2;k > > 120 > > ; 2 M 1 8v3 2 H ðxÞ; v3 ¼ @m v3 ¼ 0 on c0 ; 8w2 2 H ðxÞ: Let us observe that the unique solution of Eq. (58) is given by u3 0 and u2 0 in x, while Eq. (57) is a pure membrane problem. 5.1.1. Numerical solution In fact, as explained in paragraphs 3 and 4, the implementation use Ganev and DimitrovÕs ﬁnite element for the approximation of u1 , u2 , u0 , u1 , u2 and Argyris et al.Õs ﬁnite element for the approximation of u3 . In this framework, the data are taken according to Table 2. For the numerical experiments, we have used the following data: ‘ ¼ 5 cm, L ¼ p 2:5 cm, e ¼ 0:1 cm and V þ ¼ 100 V. Of course, we have found back u3 ¼ 0, u2 ¼ 0 in x. The associated plane deformations of the plate are represented in Fig. 7 with an ampliﬁcation factor by 500 (the maximum of the displacement is close to 10 lm); for clarity, the boundary of the undeformed plate is displayed in dotted line. When V þ ¼ 100 V, the plate takes extension while, for V þ ¼ 100 V, we obtain retraction. Example 2. Same data than in Example 1 but the plate is no longer clamped. The clamped condition in Example 1 introduces unsymmetric deformations. Here we do not impose these clamped conditions. In order to kill the associated rigid body motions, we ﬁx the center of the plate, i.e., Table 2 Constraints on some degrees of freedom u

u;n1

u;n2

u;n1 n1

u;n1 n2

u;n2 n2

u;m

0 0 0

Free Free 0

· · 0

· · 0

· · Free

Free Free 0

u;n1

u;n2

u;m

0

0

0

0

0

0

2

On n ¼ 0 (clamped condition) u1 0 u2 0 0 u3 u In x (given potential) e2 Vþ u0 þ u2 8 2 Vþ u1 e

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4063

Fig. 7. Deformation of the clamped piezoelectric plate (ampliﬁed by a factor 500). (a) V þ ¼ 100 V and (b) V þ ¼ 100 V.

‘ ‘ ‘ u1 0; ¼ u2 0; ¼ u3 0; ¼0 2 2 2 and we prevent any rotations by imposing ‘ ‘ ‘ ðu1;n2 u2;n1 Þ 0; ¼ u3;n1 0; ¼ u3;n2 0; ¼ 0: 2 2 2

ð59Þ

ð60Þ

In this case, we ﬁnd back u3 0, u2 0; the associated deformations are symmetric and displayed in Fig. 8 with an ampliﬁcation factor of 500. Here again, we observe a dilatation for V þ ¼ 100 V and a retraction for V þ ¼ 100 V. Example 3. Same than Example 1 but only a part of the upper and lower faces of the plate is symmetrically anodised. Let us consider Example 1 but instead to metallize the whole upper and lower faces of the plate, we only metallize a part which is symmetrically distributed with respect to the middle surface (see Fig. 9). To ﬁx the ideas, let us say that the electrodes are located on x1 f 2eg with x1 ¼ L20 ; L20 ½0; ‘½; 0 < L0 < L. As above, we apply the following potentials n eo n eo u ¼ 0 on x1 and u ¼ V þ on x1 þ : ð61Þ 2 2 According to the partition of CE0 discussed in (I,Section 3.2 and Fig. 3), we obtain the following partition of x (see Fig. 10):

Fig. 8. Deformation of the piezoelectric plate without any clamped condition. (a) V þ ¼ 100 V and (b) V þ ¼ 100 V.

4064

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Fig. 9. The metallized parts of the plate.

Fig. 10. The reference domain and the 3D-plate.

x ¼ x0 [ x1 with x0 ¼

i

i i L L h L L0 h L0 L h 0 0 ; ; 0; ‘½ ; x1 ¼ ; 0; ‘½ [ 0; ‘½: 2 2 2 2 2 2 2

þ

þ

Here again, we obtain, by similarity with (55), u0 þ e8 u2 ¼ V2 and u1 ¼ Ve in x1 . Then, the continuous two-dimensional problem (I,(151)) becomes (compare with Eqs. (57) and (58)): 9 1 Find ua 2 H 1 ðxÞ with ua ¼ 0 on cM > 0 ; and find u1 2 H ðx0 Þ such that > > R > pﬃﬃﬃ 1 2 abkl > > eC c ð~ u Þc ð~ v Þ a dn dn > ab kl x > > Z = 2 p ﬃﬃ ﬃ e 1 2 3kl 3ab 33 kb ð62Þ a dn dn þ e p e u1 ckl ð~ vÞ ðp e cab ð~ uÞ d u1 Þw1 þ d u1;b w1;k > 12 x0 > > R > pﬃﬃﬃ > > ¼ x1 V þ p e3kl ckl ð~ vÞ a dn1 dn2 > > > ; 1 M 1 8va 2 H ðxÞ; va ¼ 0 on c0 and 8w1 2 H ðx0 Þ: 9 Find u3 2 H 2 ðxÞ with u3 ¼ @m u3 ¼ 0 on cM > 0 ; > > 2 þ > > e V > 1 > in x1 such that and find u0 and u2 2 H ðxÞ with u0 þ u2 ¼ > > 8 2 = R p ﬃﬃ ﬃ 1 2 kl k k ðM qkl ð~ vÞ S0 w0;k S2 w2;k þ T2 w2 Þ a dn dn ¼ 0 x > > > e2 > 2 2 1 > 8v3 2 H ðxÞ; v3 ¼ @m v3 ¼ 0 and 8ðw0 ; w2 Þ 2 ðH ðxÞÞ ; w0 þ w2 ¼ 0 in x1 ; > > > 8 > ; where the expressions of M kl ; S0k ; S2k and T2 are given by relations ð53Þ:

ð63Þ

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4065

Once more, these equations are independent and Eq (63) involves u3 0;

u2 0 and

u0

Vþ 2

in x:

ð64Þ

For the implementation, we have considered the general equations (62) and (63). The clamped condition along the side n2 ¼ 0 and the conditions (61) on the potential lead to ﬁx some degrees of freedom in the same manner than in Table 2; of course, for the potential degrees of freedom, these restrictions are only applied in x1 . The numerical experiments give back the results (64) and the associated deformations of the plate are represented in Fig. 11 with an ampliﬁcation factor of 500. Example 4. Same than Example 3 but the plate is no longer clamped. The transition between Examples 3 and 4 is entirely similar to that between Examples 1 and 2. That means only the clamped condition is relaxed and, in order to kill the rigid body motions, we apply conditions (59) and (60). The associated deformations are represented in Fig. 12 with an ampliﬁcation factor of 500. Example 5. Same than Example 1 but the position of the electrodes on the upper and lower faces is no longer symmetric. Now, the deformation includes some bending eﬀects. The new data are illustrated in Fig. 13: the upper face is the same than in Example 3 while the lower face is completely metallized. With the partition of CE0

Fig. 11. Deformation of the piezoelectric plate shown in Fig. 10. (a) V þ ¼ 100 V and (b) V þ ¼ 100 V.

Fig. 12. Deformation of the piezoelectric plate of Fig. 10 without clamped boundary condition. (a) V þ ¼ 100 V and (b) V þ ¼ 100 V.

4066

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Fig. 13. The reference domain, the 3D plate and its metallized parts.

introduced in (I,Section 3.2 and Fig. 3), we obtain x ¼ x1 [ x3 , with x3 ¼ ð L2 ; L20 ½0; ‘½Þ [ ð L20 ; L2 ½ 0; ‘½Þ and x1 ¼ L20 ; L20 ½ 0; ‘½. The potentials are applied as follows: u ¼ 0 on x f 2eg and u ¼ V þ on x1 fþ 2eg, so that we obtain the following relations for the development of the potential through the thickness: u0 þ

e2 Vþ u2 ¼ 8 2

and

u1 ¼

Vþ e e2 in x1 ; u0 u1 þ u2 ¼ 0 in x3 ¼ x x1 : 2 e 8

ð65Þ

The last constraint (65) prevents the splitting in two independent variational formulations like in (62) and (63). Thus, the continuous two-dimensional problem (I,(151)) becomes: ~KL ðxÞ and ðu0 ; u1 ; u2 Þ 2 VPu ðxÞ such that Find ~ u2V d Z pﬃﬃﬃ fN kl ckl ð~ vÞ þ M kl qkl ð~ vÞ S0k w0;k S1k w1;k S2k w2;k þ T1 w1 þ T2 w2 g a dn1 dn2 ¼ 0 x

8ð~ v; ðw0 ; w1 ; w2 ÞÞ 2 V~KL ðxÞ VP0 ðxÞ; where the diﬀerent terms are deﬁned by relations (53) and (54) and where V~KL ðxÞ ¼ f~ v 2 ðH 1 ðxÞÞ2 H 2 ðxÞ; vi ¼ om v3 ¼ 0 on cM 0 g; e2 1 VPud ¼ ðH 1 ðxÞÞ3 \ ðh0 ; h1 ; h2 Þ such that h0 þ h2 ¼ V þ and 2 8 1 þ e e2 h1 ¼ V in x1 ; and h0 h1 þ h2 ¼ 0 in x3 : e 2 8 In particular, u3 and u2 are no longer zero while u0 6 V þ =2. All these remarks are clearly observed in the numerical experiments. In Figs. 14 and 15, we present for V þ ¼ 100 V, the isovalues of u3 and the values of u0 . It is worth to note that the values of u3 are really small compared with those of u1 and u2 . Thus, we qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ PNOE 2 PNOE 2 2 obtain k¼1 ðu3 ðkÞÞ = k¼1 ðua ðkÞÞ 6 10 , for a ¼ 1; 2, where NOE is the number of nodes of the mesh. On the other hand, if we denote by ~ u1 the solution of Example 3 and by ~ u2 the present solution, we observe that the membrane displacement in both examples are really close, i.e., ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ sP NOE 2 1 2 k¼1 ðua ðkÞ ua ðkÞÞ 6 102 : PNOE 1 2 ðu ðkÞÞ a k¼1

ð66Þ

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4067

Fig. 14. Isovalues of u3 in lm.

Fig. 15. Values of u0 in kV.

5.2. Some piezoelectric circular cylindrical thin shell problems Now, let us consider a similar study than in Section 5.1 but with a half-cylinder of length ‘, radius R and thickness e. Its geometrical deﬁnition appears in Fig. 16. We assume that this thin shell is constituted by the same piezoelectric material than the plate analyzed in Section 5.1. Since the system of curvilinear coordinates which describes this half-cylinder is orthogonal, the abkl tensors C , p ekkl and dij have the same form than in (48)–(51), but now, due to the curvature of the shell the tensorial components cab ð~ uÞ and qab ð~ uÞ are replaced by (combine (I,(152)) with the geometrical pa~S given in Fig. 16): rameters derived from the mapping U 1 1 cab ð~ uÞ ¼ ðua;b þ ub;a Þ þ u3 d1a d1b ; 2 R

qab ð~ uÞ ¼ u3;ab þ

1 1 u3 d1a d1b þ ðd1a u1;b þ d1b u1;a Þ: R2 R

Here again, we consider diﬀerent loadings and boundary conditions for this thin shell.

4068

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Fig. 16. The circular cylindrical thin shell: ~ US ðn1 ; n2 Þ ¼ ðR sinðn1 =RÞ; n2 ; R cosðn1 =RÞÞ.

Example 6. A clamped circular cylinder with an applied potential on the whole upper and lower surfaces. Let us assume that the half-cylinder displayed in Fig. 16 is • clamped along the side n2 ¼ 0; • submitted to a potential V þ ¼ 100 V on its upper face and a potential V ¼ 0 on its lower face. Both potentials are applied through electrodes which cover the entire upper and lower faces. Then, we can check that the associated deformations of this shell are described by a system similar to (57) and (58) but now these equations are no longer independent since the curvature relates ua and u3 . For the numerical experiments, we take ‘ ¼ 5 cm;

R ¼ 2:5 cm;

e ¼ 0:1 cm; V þ ¼ 100 V;

V ¼0 V

and the clamped conditions are speciﬁed according to Table 2. The associated deformation of the cylinder are represented in Fig. 17 with an ampliﬁcation vector of 500. As mentioned before, all the unknowns are coupled and, in particular, u3 is no longer zero. In Fig. 18, we present the value of u3 (in lm) at points ðn1 ¼ 0; n2 ¼ ‘Þ and ðn1 ¼ 0; 625p cm; n2 ¼ ‘Þ as a function of the curvature radius R. Concerning the potential, its quadratic development involves now a term u2 which is no longer zero and a term u0 which is no longer constant: thus the potential is no longer linear through the thickness. In Fig. 19, we present the isovalues of u0 (when R ¼ 2:5 cm) and, in Fig. 20, we present the

Fig. 17. Deformation of the half clamped cylinder. (a) V þ ¼ 100 V, V ¼ 0 V and (b) V þ ¼ 100 V, V ¼ 0 V.

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Fig. 18. Variation of u3 as a function of R measured at points ð0; ‘Þ and ð0; 625p cm; ‘Þ for V þ ¼ 100 V.

Fig. 19. Isovalues (in kV) of u0 for the half-cylinder for V þ ¼ 100 V.

Fig. 20. Maximum and minimum values of u0 (in kV) as functions of R for V þ ¼ 100 V.

4069

4070

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Fig. 21. Relative errors as functions of R.

maximum and minimum values of u0 for various values of R: it is worth to note that these values converge to the plate value when R is increasing. Let us make some comments about the relative errors done on the displacement if we replace the quadratic development of the potential u by a linear one. Of course, if we assume a linear potential through the thickness, it is completely determined by its values on the upper and lower faces. Then, it acts like an applied load to the shell and the only remaining unknown is the mechanical displacement ~ u. By denoting ~ ulin the solution associated with a linear development of u, and ~ u the two-dimensional solution associated with the quadratic development of u through the thickness, we deﬁne the relative errors Erri ¼ ﬃ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ PNOE 2 PNOE 2 , i ¼ 1; 2; 3, where NOE means the number of nodes of the mesh. ðu ðkÞ u ðkÞÞ = ðu ðkÞÞ i lin;i i k¼1 k¼1 These relative errors are presented in Fig. 21 for various values of R. Of course, when R is increasing, ua converges to the plate solution for which the potential u takes a linear development: thus, we can check that Erra converges to 0. Example 7. Same than Example 6 except that the shell is no longer clamped. The statement is entirely similar; the only diﬀerence is the replacement of clamped conditions by a set of conditions which prevents the rigid body motions, i.e., conditions of type (59) and (60). Corresponding results are presented in Fig. 22 with an ampliﬁcation factor of 500. Example 8. Same than Example 6 except that the electrodes are symmetrically localised only on a part of the upper and lower faces of the shell. These new data are displayed in Fig. 23. The associated results are presented in Fig. 24 with an ampliﬁcation factor of 500. Example 9. Same than Example 8 but now the cylinder is no longer clamped. The associated results are presented in Fig. 25. Example 10. Same than Example 8 but now the potential is applied on a part of the upper surface and on the whole of the lower surface.

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4071

Fig. 22. Deformation of the half unclamped cylinder. (a) V þ ¼ 100 V, V ¼ 0 V and (b) V þ ¼ 100 V, V ¼ 0 V.

Fig. 23. The electrodes are partly and symmetrically distributed on the upper and lower faces.

Fig. 24. Deformations of the clamped half-cylinder when the potential is applied only on symmetrical parts of the upper and lower faces. (a) V þ ¼ 100 V, V ¼ 0 V and (b) V þ ¼ 100 V, V ¼ 0 V.

4072

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

Fig. 25. Deformations of the unclamped half-cylinder when the potential is applied only on symmetrical parts of the upper and lower faces. (a) V þ ¼ 100 V and (b) V þ ¼ 100 V.

Fig. 26. Relative diﬀerences di as a function of R.

If we denote by ~ u1 (resp. ~ u2 ) the qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬃ discrete solutions of Example 8 (resp. 10), we deﬁne di ¼ PNOE 1 P 2 NOE 2 1 k¼1 ðui ðkÞ ui ðkÞÞ = k¼1 ðui ðkÞÞ, where NOE means the number of nodes of the mesh. Then, Fig. 26 shows the evolution of the diﬀerences di as functions of the curvature radius R of the cylinder. Like for the plate examples, the behaviour of the half-cylinder depends on the location of the electrodes. u2 of Examples 8 and 10 are of the same order than those The diﬀerences between the solutions ~ u1 and ~ obtained for the plates when R is increasing, i.e., of order 102 according to relation (66) and Fig. 26.

References [1] M. Bernadou, C. Haenel, Modelization and numerical approximation of piezoelectric thin shells. Part I: The continuous problem, Comput. Methods Appl. Mech. Engrg. 192 (37–38) (2003) 4003–4043. [2] M. Bernadou, Finite Element Methods for Thin Shell Problems, John Wiley, Chichester, 1996. [3] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [4] H.G. Ganev, Tch.T. Dimitrov, Calculation of arch dams as a shell using IBM-370 Computer and curved ﬁnite elements, in: Theory of Shells, North-Holland, Amsterdam, 1980, pp. 691–696.

M. Bernadou, C. Haenel / Comput. Methods Appl. Mech. Engrg. 192 (2003) 4045–4073

4073

[5] J.H. Argyris, I. Fried, D.W. Scharpf, The TUBA family of plate elements for the matrix displacement method, Aero. J. Royal, Aeronaut. Soc. 72 (1968) 701–709. [6] J.N. Lyness, D. Jespersen, Moderate degree symmetric quadrature rules for the triangle, J. Inst. Maths. Applicat. 15 (1975) 19–32. [7] D.A. Dunavant, High degree eﬃcient symmetrical gaussian quadrature rules for the triangle, Int. J. Numer. Methods Engrg. 21 (1985) 1129–1148. [8] M. Bernadou, P.L. George, A. Hassim, P. Joly, P. Laug, B. Muller, A. Perronnet, E. Saltel, D. Steer, G. Vanderbork, M. Vidrascu, MODULEF: Une Bibliotheque Modulaire dÕElements Finis, Editions INRIA, Rocquencourt, Deuxieme edition, 1988. [9] MODULEF. Available from . [10] C. Haenel, Modelisation, Analyse et Simulation Numerique de Coques Piezoelectriques, These de lÕUniversite Pierre et Marie Curie, Paris, France, 2000.

Modelization and numerical approximation of piezoelectric thin shells

Modelization and numerical approximation of piezoelectric thin shells

Recommend Documents