Error estimation through the constitutive relation for Reissner–Mindlin plate bending finite elements

Computers and Structures 73 (1999) 615±627

www.elsevier.com/locate/compstruc

Error estimation through the constitutive relation for Reissner±Mindlin plate bending ®nite elements P. Boisse a, b,*, S. Perrin c, G. Cognal a, K. Hadjeb a a

Laboratoire de ModeÂlisation et MeÂcanique des Structures, URA CNRS 1776, ENSAM/Universite Paris 6/E.N.S. Cachan, Paris, France b Ecole SupeÂrieure de l'Energie et des MateÂriaux, Orleans, France c CADOE S.A., Villeurbanne, France Received 8 January 1997; accepted 3 October 1998

Abstract This paper deals with the so-called `error estimation through the constitutive relation' for Reissner±Mindlin plate bending ®nite element analysis, which is based on the calculation of statically admissible stress resultant ®elds. Further an energy norm through the constitutive relation is introduced to quantify the di€erence between the calculated statically admissible stress resultant ®elds and the stress resultants given by the ®nite element analysis. The statically admissible stress resultant ®elds are obtained by local calculations at the node or at the element level. The procedures involve an optimization for estimating the statically admissible quantities as close as possible to the ®nite element result, and thus lead to a more accurate error estimation, considering the Pragger±Synge theorem. The set of examples shows the good agreement between the proposed indicator and the error obtained from the analytical solution. It is concluded that plate bending analysis is a ®eld of application which is well ®tted for error estimation through the constitutive relation. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Error through the constitutive relation; Reissner±Mindlin; Plate bending; Statically admissible stress resultants

1. Introduction Knowing the quality of a ®nite element analysis is very important for the engineer. For this reason, error estimators are nowadays becoming an essential requirement for structural calculations. In addition to the knowledge of the quality of the analysis, they

* Corresponding author. Tel.: +33-144-24-62-85; fax: +33144-24-64-68. E-mail address: [email protected] (P. Boisse)

allow automatic remeshing and constitute a main part of the so-called `adaptive meshing' technique. Several methods or families of methods have been proposed in order to estimate the error of a given ®nite element analysis. Most of them give a posteriori values of the estimated error. A family of indicators is based on the element equilibrium residual calculation and is dealt with in the work of Babuska and Rhieinboldt [1,2] and Kelly et al. [3]. The most commonly used error estimators have been proposed by Zienkiewicz and Zhu [4,5], based on recovery methods, whereas Refs. [6±9]

0045-7949/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 9 8 ) 0 0 2 2 4 - 7

616

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

Nomenclature uh E …uh † s sh uà s^ C rO, rE e1, e2, e3 x1, x2, x3 np a, b k T, M, Q T M Q b

D Ds Sm @St @Su p3  M q e3 i DP(C ) G bGE, bGE '

nGE

®nite element displacement deformation related to uh …EE…uh † ˆ T 1 2 …grad u ‡ grad u†† stress tensor in the case of plates, s ˆ sab ea eb ‡ sa3 …ea e3 ‡ e3 ea † stress tensor solution of the ®nite element analysis kinematically admissible (KA) displacement ®eld derived from uh statically admissible (SA) stress elastic Hooke's tensor error through the constitutive equation on domain O, on element E natural orthogonal basis in plate analyses, e3 is normal to the plate coordinates in e1, e2, e3 unit vector in plane e1, e2 indices taking the value 1 or 2 indices taking the value 1±5 stress resultants: membrane forces, bending moments, shear force membrane tensor„ such that T ˆ ‡h=2 Tab ea eb ; Tab ˆ ÿh=2 sab dx 3 bending moment tensor „ ‡h=2 such that M=Mabea eb; Mab ˆ ÿh=2 x 3 sab dx 3 shear „force vector such that Q=Qaea; ‡h=2 Qa ˆ ÿh=2 sa3 dx 3 plate bending Hooke operator plate transverse shear Hooke operator mid-surface of a plate part of the boundary of Sm with prescribed loads part of the boundary of Sm with prescribed displacements shear load component prescribed on Sm prescribed moment density on @Sm prescribed shear load density on @Sm nodal component at node i, of a scalar density C along line D, i DP…C † ˆ „ i D N C dG segment, edge of an element or common to two adjacent elements orientation functions of two adjacent elements E and E' with common edge G. bGE=1 and bGE '=ÿ1 if E' has a greater index than E in the global numbering exterior normal to elements E, along edge G, in plane e1, e2 of the plate

nG lG Ni N B s im

MGm

QGm

^G M Q^ G C^ kG $

j Z ‰tk Š [skl] [F kE] gE gi Gint

normal to segment G, in plane e1, e2 such that nG=bGEnGE=bGE 'nGE ' length of edge G linear interpolation function (scalar) related to node i displacement interpolation matrix of an element strain interpolation matrix of an element mean stress classically calculated by ®nite element post-processors at node, i, p being the number of element Ej connected to the node i, s im ˆ …1=p†Sjˆ1,p s h …Ej † mean bending moment density obtained from the ®nite element computation along the edge G linking up nodes i and „ ‡h=2 j, and such that MGm ˆ ÿh=2 x 3 …N i s im ‡ N j s mj †
nG dx 3 mean shear load density obtained from the ®nite element computation along the edge G linking up nodes „ ‡h=2 i and j, and such that QGm ˆ … ÿh=2 …N i s im ‡ N j s mj †
nG dx 3 †
e3 `in equilibrium' bending moment density on G `in equilibrium' shear load density on G ^ G , Q^ G , M ^ G ˆ C^ 1 e1 ‡ component of M G C^ 2G e2 , Q^ G ˆ C^ 3G e3 normal to the plate virtual displacement of a point of the mid-surface of the plate satisfying $ ˆ 0 on @Su virtual rotation of the normal to the mid-surface of the plate satisfying j ˆ 0 on @Su virtual displacement vector satisfying Z ˆ 0 on @Su and such that Z ˆ$e3 ‡ x 3 e3 ^ j column matrix of the components of the vector t column matrix of the components of the second-order tensor s column matrix of the components of exterior„ nodal load„ and moments ‰F kE Š ˆ E NT ‰ fk Š dV ÿ @ Em NT ‰tk Š dS set of edges G of element E set of edges G connected to node i interface between two sub-triangles

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

MÃj, MÃk MÃPSijk

nodal value of MÃ at node i for the subtriangles j and k connected by Gint moments added at point Pij to the linear ®eld de®ned from values at node i, j and G

employ a superconvergent patch recovery technique. The concept of `error through the constitutive relation' was introduced very early by LadeveÁze [10]. It is based on the construction of strictly admissible static and kinematics ®elds from the ®nite element results. When considering them together, their lack of veri®cation of the constitutive relation de®nes the error which is calculated as a strain energy. One of the main interests of the method is, as shown in Ref. [11] from the Prager± Synge theorem, that the error through the constitutive relation provides an upper bound of the error referred to the exact solution. The method also has a strong mechanical base as the exact solution of the problem is strictly admissible. In particular, the lack of veri®cation of the boundary conditions is taken into account. The error through the constitutive relation method has been adapted to a large class of two-dimensional (2D) and three-dimensional (3D) linear and non-linear analyses. The important works are cited here, for example, for 2D and 3D linear analyses [11±16], elastoplasticity [17±19,31], viscoplasticity [21], vibrations [17,20] and large strains [21]. In practical ®nite elements analyses, a large class of structures presents one of the dimensions as being smaller than the others. They can be analysed using the plate or shell theories. Finite elements based on the kinematical assumptions of these theories are numerous and they are ecient and very often used in practice. The objective of this paper is to present an application of the error estimator through the constitutive relation method to Mindlin plate bending static analyses [22] and also to bring out the ecacy of the method in this ®eld of applications. In the case of a C 0 linear triangular plate element [23,24], the construction of the admissible resultant stress ®elds is presented. The ®rst step gives, for each element, a set of load and moment densities along the edges of the elements, which ensure the equilibrium of each element. These densities are the sum of densities derived from the ®nite element solution and of densities equilibrated on each element side. In the second step, the statically admissible stress resultants are calculated within each element from the previous loads and moments densities. These two steps allow an optimization of the static ®elds in order to get a more accurate indicator. The calculations are local. Hence, the approach is ecient,

MÃSPkiG

617

moments added at point PiG to the linear ®eld de®ned from values at node i, j and G

from the numerical point of view, for the large number of degrees of freedom (d.o.f.). Examples are given for the case of rectangular, circular and square bending plates subjected to normal pressure with di€erent boundary conditions. They show the good agreement between the indicator proposed for plate bending ®nite element analysis and the error obtained from the analytical approach using the plate bending theory. 2. Error through the constitutive relation for Mindlin plate bending analysis 2.1. De®nition of the error The study is restricted to geometric and material linearities. The stress ®eld s over the plate corresponds to stress resultant T, Q, M, as de®ned in the notation section at the beginning of the paper. We consider plate bending assuming the membrane stress resultant T is equal to zero. Hence the stress ®eld is de®ned by the bending moment tensor M and the transverse shear stress vector Q. QÃ, Mà are statically admissible stress resultants if they satisfy the local equilibrium à ‡ p3 ˆ 0 div M à ˆ0 à ÿQ div Q

…1†

and the boundary conditions Ã
np ˆ q Q

 Ã
np ˆ M M

…2†

Within a displacement approach, uh is kinematically admissible and uÃ, the admissible ®eld, used further is such that uà ˆ uh

…3†

The error on the constitutive relation rO on the domain O is de®ned from the lack of constitutive law veri®cation of the kinematically admissible displacement ®eld uà and statically admissible stress ®eld sà [11,13] as rO ˆ

ks^ ÿ C:EE…uà †kO ks^ ‡ C:EE…uà †kO

…4†

618

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

with … ks s kO ˆ

O

s :Cÿ1 :s s dV

1=2

…5†

In order to de®ne a map of the error over the whole structure, the contribution of each element to this error is considered as rE ˆ

ks^ ÿ C:EE …uà †kE ks^ ‡ C:EE …uà †kO

with r2O ˆ

X

r2E

…6†

elements

It has been shown in [11] that the error through the constitutive relation provides an upper bound of the error referring to the exact solution. This property is based on the Prager±Synge theorem and is of major importance for the practical use of the error estimation. Using the bending moment tensor and shear stress vector, the error through the constitutive equation can be written in the following form, for the case of Mindlin plate bending analyses: … rO ˆ …O O

admissible stress resultant ®elds as ecient as possible, i.e. which will lead to an error through the constitutive relation as close as possible to the error referring to the exact solution. Accounting for the fact that the error through the constitutive relation Eq. (7) is an upper bound of the `exact' error, the error will be minimized by choosing SA ®elds as close as possible to the stress ®eld results coming from the ®nite element computation. This is done at each step where the choice is possible. The calculation of the admissible stress resultant ®eld will be made as an optimization problem, where the SA ®elds are calculated as admissible but also rendering the error as small as possible. Finally, the method of calculation will be simple and economical from the numerical cost point of view. In particular, the computations will be made locally, i.e. they will concern an element, a node or an element edge. In this way, the computational e€ort will be proportional to the number of d.o.f. for the given problem

à ÿ Qh †:Dsÿ1 :…Q à ÿ Qh † dS à ÿ Mh †:Dbÿ1 :…M à ÿ Mh † ‡ …Q …M à ‡ Qh †:Ds :…Q à ‡ Qh † dS à ‡ Mh †:Db :…M à ‡ Mh † ‡ …Q …M ÿ1

ÿ1

2.2. Basic feature of the statically admissible stress resultant construction The calculation of the error through the constitutive relation, Eq. (7), requires the determination of statically admissible stress resultants QÃ, MÃ. They have to strictly satisfy the equilibrium equations (1) over the whole domain O and the load boundary conditions, Eq. (2), on the boundary @Sm. As it has been used in the case of 2D or 3D problems [11], the determination of statically admissible stress resultants QÃ, MÃ on the domain O is replaced by the calculation of QÃE, MÃE ®elds satisfying both the equilibrium on each element E as well as the boundary conditions with load and bending moment densities determined on each element edge. These densities must lead to systems on equilibrium on each element and ensure that the reciprocal actions between two neighbouring elements E and E ' are opposite. Consequently, the set of stress ®elds QÃE, MÃE over the total domain O is statically admissible. Statically admissible stress resultant ®elds are not unique. Therefore, it is the main aspect of the calculation which will be made in order to ®nd statically

1=2 1=2

…7†

under consideration (in contrast with the global calculations, the size of which increases in a drastic manner with the number of degrees of freedom) and will allow for the error estimate in case of structures with a large number of d.o.f for which adaptivity is a very important feature. 2.3. C0 plate bending ®nite elements The error calculation is made using three node C 0 plate ®nite elements with linear interpolations. The triangular shell ®nite element proposed in [23±25] will be used in the examples. Other plate elements with linear interpolation can be considered such as Belytshcko three-node plate element proposed in [26]. These elements use di€erent approaches in order to avoid shear locking phenomena (a mixed interpolation technique in [23±25] and a stress projection in [26]). This has no consequence in the present work if the interpolation matrix used here takes into account these modi®cations. Furthermore, the accuracy of plate analyses performed using any ®nite element can be determined using the method described in this paper and based on the use of linear interpolation functions. The

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

de®nition of the error Eq. (9) needs the ®nite element results Qh, Mh, which are obtained from a given plate element and QÃ, MÃ, which are derived from using the method described in this paper, and use a linear interpolation.

619

…

… E

E …Z Z†:s^ a dV ˆ … ‡

@E

… E

f
Z dV ‡

@ Em

p3 e3
Z dS …9†

^ Ga † dG $bGE …Q^ Ga
e3 † ‡ j
…bGE M

Considering this equality in the case of the interpolated ®elds, Eq. (9) becomes 3. Load and moment densities on the element edges

…

3.1. The two stages of the calculation

E

^ G and shear A set of bending moments densities M load densities Q^ G are calculated on each element edge G representing the action of the element E on the neighbouring element E'. G separates elements E and E', and the number of E' is greater than those of E in the global numbering. They satisfy the following properties: ^ G , Q^ G ensures the equilibrium on . The set of all M each element taking into account the possible exterior load within the element. ^ G and Q^ G are equal to the prescribed exterior . M loads on the part of the frontier @St. ^ G and Q^ G are as close as possible to the results of . M the ®nite element analysis. ^ G , Q^ G is restricted to those . The determination of M with a linear evolution on each edge G and such that ^ GˆM ^ Ga ‡ M ^ Gb M

Q^ G ˆ Q^ Ga ‡ Q^ Gb

…8†

^ Ga , Q^ Ga are bending moment and shear load densities M ensuring the equilibrium of each element. They satisfy the local boundary conditions and lead to a statically admissible stress ®eld s^ a within each element, the interior nodal loads of which are equal to those of the ®nite element stress solution s h . ^ Gb , Q^ Gb are bending moment and shear load denM sities such as the resultant action of each edge G is ^ G , Q^ G as equal to zero and they lead to densities M close as possible to the ®nite element solution. The ^ Gb and Q^ Gb appears as an correction made by M ^ improvement of MGa and Q^ Ga , preserving the equili^ Ga , M ^ Gb , Q^ Ga , brium of each element. Furthermore, M Q^ Gb variations are chosen to be linear along the edge G. 3.2. Calculation of the ®rst bending moment and shear ^ Ga and Q^ Ga load densities M The principle of virtual works applied to the plate ^ Ga and bGE Q^ Ga is veri®ed element E subjected to bGE M by accounting for the statically admissible nature of the related stress ®eld s^ a as, 8$ and j equal to zero on @Su,

…

BT ‰…s^ a †kl Š dV ÿ … ÿ

@E

E

NT ‰ fk Š dV ÿ

… @ Em

NT ‰tk Š dS …10†

N ‰bGE C^ kGa Š dG ˆ 0 T

By de®nition, the stress ®elds s^ a and s h have equal interior nodal e€orts … E

BT ‰…s^ a †kl Š dV ˆ

… E

BT ‰…sh †kl Š dV

…11†

consequently i ^k @ E P…bGE C Ga †

… ˆ

T

E

B ‰…sh †kl Š dV

k

ÿF kEext

…12†

Eq. (12) gives i@EP(bGECà kGa) from the ®nite element ^ Ga , calculation. The nodal components iGP(Cà kGa) of M Q^ Ga are calculated from those nodal components i Ãk @EP(bGEC Ga). i k à GP(C Ga) are prescribed on the frontier @St. The study is limited to linear exterior loads and moments over each element edge, which is usually sucient in ®nite element simulations. Considering a node i, the nodal components i Ãk @EP(bGEC Ga) are related to the components on each face i ^k @ E P…bGE C Ga †

ˆ

X G2gE \gi

… bGE

G

N i C^ kGa dG

…13†

At each node i, the following system is obtained: 8E containing node i X G2gE \gi

bGE Gi P…C^ kGa † ˆ@ EiP…bGE C^ kGa †

…14†

The solution of system Eq. (14) is unique if the component Cà kGa is imposed on one or two edges of gi. When the node i is interior, the system is undetermined, allowing the best choice with respect to an eciency criterion. In that case, a supplementary equation is obtained from the minimization of the distance between the unknown nodal components and of MGm, QGm, i.e.

620

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

X G2GE \Gi

‰iG P…C^ kGa † ÿGi P…C kGm †Š2 ˆ minimum

…15†

When deriving this equation with respect to a component iGP(Cà kGa), a determined system is obtained, the resolution of which gives the nodal load components i Ãk GP(C Ga). ^ Ga , Q^ Ga being The components Cà kG (k $ (1,3)) of M linear along the edge G linking up i and j, the knowledge of nodal load components iGP(Cà kGa) and jGP(Cà kGa) à kj determines Cà kGa. Denoting Cà ki Ga and C Ga as the values k à of C Ga at nodes i and j, j ^ kj C^ kGa ˆ N i C^ ki Ga ‡ N C Ga

…16†

and i ^k G P…C Ga †

… ˆ

G

i j ^ kj …N i †2 C^ ki Ga ‡ N N C Ga dG

…17†

Taking into account the fact that the interpolation functions are linear, Eq. (17) written for node i and j leads to 2 i ^k …C^ ki ‰2 P…C Ga † ÿGj P…C^ kGa †Š Ga † ˆ lG G

…18†

^ Ga are in equilibrium The densities bGE Q^ Ga and bGE M ^ Gb will allow on element E. The addition of Q^ Gb , M densities to be obtained closer to QGm, MGm, while conserving the equilibrium on each element.

3.3. Calculation of the complementary bending moment ^ Gb and Q^ Gb and shear load densities M ^ Gb , Q^ Gb are calculated in order to The densities M ^ G , Q^ G as close as possible to render the variations of M those of MGm, QGm, while preserving the equilibrium ^ Ga , Q^ Ga . properties of M The load and moment eGb and m 0Gb resulting from ^ Gb on an element edge G are the densities Q^ Gb and M considered as: … … ^ Gb eGb ˆ Q^ Gb dG m 0Gb ˆ M G

…

dG ‡

G

G

…19†

eGb ˆ e3Gb e3 e3Gb ˆ

G

… G

m02 Gb ˆ

C^ 1Gb dG ‡ … G

… G

C^ 2Gb dG ‡

x 2 C^ 3Gb dG …22†

… G

ÿ x 1 C^ 3Gb dG

The density components Cà 1Gb, Cà 2Gb and Cà 1Gb are calcu02 lated in order to render e 3Gb, m 01 Gb, m Gb, equal to zero 1 ^ G, à and to lead to three components C G, Cà 2G, Cà 3G of M ^QG , with slopes equal to those of the components Cà 1Gm, Cà 2Gm, Cà 3Gm of MGm, QGm, respectively. This leads to six conditions, that determine Cà 1Gb, Cà 2Gb, Cà 3Gb, accounting for their linear nature: C^ 1Gb ˆ C^ 1Gm ÿ C^ 1Ga 1 ÿ lG

… G

…C 1Gm ÿ C^ 1Ga ‡ x 2 C^ 3Gb † dG

…23†

C^ 2Gb ˆ C^ 2Gm ÿ C^ 2Ga ÿ

1 lG

… G

…C 2Gm ÿ C^ 2Ga ÿ x 1 C^ 3Gb † dG

1 C^ 3Gb ˆ C^ 3Gm ÿ C^ 3Ga ÿ lG

… G

…C 3Gm ÿ C^ 3Ga † dG

…24†

…25†

Remark. If the previous method to achieve densities closer to the ®nite element results by adding densities with a resultant equal to zero on each edge G is used in the case of a membrane (or 2D) analyses, it can be shown that the nature of the unknowns and equations does not permit the slope of the admissible membrane components to be rendered equal to those of ®nite element results. In this case, only a minimization of the distance of the slopes is possible [27]. Hence, this shows that the method is well suited to plate bending analyses. Remark. The proposed method, used in order to render the densities on the edges closer to the ®nite element results by rendering the slope of the load and moment densities in equilibrium equal to those of the ®nite element computation is not unique. It is chosen for eciency.

OM ^ Q^ Gb dG 4. Calculation of statically admissible stress resultants within each element

eGb and m 0Gb are prescribed to be equal to zero.

…

m01 Gb ˆ

02 m 0Gb ˆ m01 Gb e1 ‡ mGb e2

C^ 3Gb dG

…20† …21†

The knowledge of shear load bGE Q^ G and moment ^ G on the frontier @E ensuring the equilidensities bGE M brium of each element E and satisfying the boundary conditions on @Et, permits calculation of stress resul-

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

tants QÃE statically admissible (SA) within each element. The calculations of bending moment MÃE tensor and shear load QÃE, being related by Eq. (16), are simultaneous and result from the resolution of a single system. Di€erent strategies are possible for the determination of these ®elds within the element. The proposed method aims to be a simple while assuring the tensor Mà and the vector Qà to be strictly SA, i.e. satisfying local equilibrium Eq. (1) at any point and boundary conditions prescribed by the bending ^ G and shear force bGE Q^ G along the exmoments bGE M terior faces of the element. In the following section, all the calculations concern an element E. Consequently the subscript E will be omitted for MÃE and QÃE denoted MÃ, Qà for simplicity. The study is restricted to a constant exterior normal density p3 over each element. 4.1. Computation of MÃ, Qà at nodes Nodal values of Mà and Qà only depend on the boundary conditions. At a node i, the boundary conditions on the two edges GE and G E' of the element E linked to the node i ^G Ã
nGE ˆ bG M M E

^ G0 Ã
nG 0 E ˆ bG 0 M M E

…26†

^ G , Q^ G have do not have any solution, if the ®elds M not been calculated under special constraints. Consequently, the construction of a regular ®eld MÃ over the triangular plate element is not possible. As proposed in the case of plane stress analyses [16], a division of Watwood [28] in three sub-triangles Sk is used (Fig. 1) based on the centre G of the element. The use of new internal faces adds constraints since the two ®elds MÃ and QÃ have to satisfy the continuity condition of the bending moments (vector) and shear force (scalar). Consequently an additional set of equations corresponding to the continuity on the internal face connected with each node i is prescribed. ^G Ã j
nGE ˆ bG M M E

…27†

621

Fig. 1. Division of the element. (a) Watwood division of the element. (b) Sub-triangle Sk of the element.

^ G0 Ã k
nG 0 E ˆ bG 0 M M E

…28†

à j
nGint ˆ M à k
nGint M

…29†

à j
nG ˆ bG Q^ G Q E E

…30†

à k
nG 0 ˆ bG 0 Q^ G 0 Q E E

…31†

à j
nG ˆ Q à k
nG Q int int

…32†

The (6  6) system Eqs. (27)±(29) determines MÃj and MÃk. The system Eqs. (30)±(32) of three equations and four unknowns admits in®nite QÃj, QÃk. This will allow an optimization presented further. Remark. A regular ®eld Qà could be found only be Eqs. (30) and (31) so that the sub-triangulation may appear unavailing for the shear force. In fact, we choose to keep this method in order to take advantage of the possible optimization rendering the computer SA stress ®eld more ecient.

4.2. Stress resultants interpolation ^ G and Q^ G leads to interpolation of Qà with the linear functions and Equilibrium Eq. (1) and the linear nature of M Mà with quadratic functions over each sub-triangle Sk (linking up to node i, j and G, Fig. 1) j ^ j G à G à S ˆ Ni Q Ãi Q Sk Sk ‡ N Sk Q Sk ‡ N Sk Q Sk k

…33†

622

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

à Sk ˆ N i M à i ‡ Nj M à j ‡ NG M à G ‡ 4N i N j M à Pij ‡ 4N i N G M à PiG ‡ 4N i N G M à pjG M Sk Sk Sk Sk Sk Sk Sk Sk Sk Sk Sk Sk Sk Sk Sk

…34†

Accounting for the linearity of bending moments MG along all the faces of each sub-triangle MÃPijSk is a priori taken equal to zero. 4.3. Equilibrium equations Taking into account the above interpolations of Mà and QÃ, equilibrium Eq. (1a) lead to (k $ [1,2,3]) @ N jSk @ N jSk @ N iSk i @NG @ N iSk i @NG j j Sk ^ G Sk ^ G …Q^ Sk †1 ‡ …Q^ Sk †1 ‡ … Q Sk † 1 ‡ …Q^ Sk †2 ‡ …Q^ Sk †2 ‡ …QSk †2 ‡ p3 ˆ 0 @x @x @x @y @y @y

…35†

In the same way, the equilibrium Eq. (1b) leads to (k $ [1,2,3])   @ N jSk @ N iSk @NG @ N iSk G @NG i j Sk ^ G Sk i ^ ^ ^ PiG †11 …MSk †11 ‡ …MSk †11 ‡ …MSk †11 ‡ 4 N Sk ‡ N Sk …M Sk @x @x @x @x @x ! @ N jSk G @ N jSk @NG @ N iSk @NG Sk Sk ^ G ^ PjG †11 ‡ ^ i †12 ‡ ^ j †12 ‡ N Sk ‡ N jSk …M … M …M …MSk †12 ‡4 Sk Sk Sk @x @x @y @y @y ! ! @ N jSk G @ N iSk G @NG @NG PiG Sk S j i k ^ PjG †12 ÿ N i …Q^ i †1 ^ †12 ‡ 4 ‡4 N Sk ‡ N Sk N Sk ‡ N Sk …M …M Sk Sk Sk Sk @y @y @y @y j

…36†

G

^ ÿ N jSk …Q^ Sk †1 ÿ N G Sk …QSk †1 ˆ 0   @ N jSk @ N iSk @NG @ N iSk G @NG i j Sk ^ G Sk i ^ ^ ^ PiG †12 …MSk †12 † ‡ …MSk †12 ‡ …MSk †12 ‡ 4 N Sk ‡ N Sk …M Sk @x @x @x @x @x ! @ N jSk G @ N jSk @NG @ N iSk @NG Sk Sk ^ G ^ PjG †12 ‡ ^ i †22 ‡ ^ j †22 ‡ N Sk ‡ N jSk …M … M …M …MSk †22 ‡4 Sk Sk Sk @x @x @y @y @y ! ! G @ N jSk G @ N iSk G @NG PiG Sk j @ N Sk i ^ ^ PjG †22 ÿ N i …Q^ i †2 ‡4 N Sk ‡ N Sk …MSk †22 ‡ 4 N Sk ‡ N Sk …M Sk Sk Sk @y @y @y @y j

…37†

G

^ ÿ N jSk …Q^ Sk †2 ÿ N G Sk …QSk †2 ˆ 0

4.4. Continuity conditions The continuity conditions at the centre imposed for the shear force are ÃG Q S1


nGG3 ˆ

ÃG Q S2


nGG3

à G
nGG3 Ã G
nGG3 ˆ M M S1 S2

…41†

à G
nGG2 ˆ M à G
nGG2 M S1 S3

…42†

à G
nGG1 ˆ M à G
nGG1 M S2 S3

…43†

…38†

à G
nG ˆ Q à G
nG Q G2 G2 S1 S3

…39†

à G
nG ˆ Q à G
nG Q G1 G1 S2 S3

…40†

The continuity condition for bending moments at the centre G is ensured by

Eqs. (41)±(43) impose the two conditions à G à G MÃG The continuity must also be S1=MS2=MS3. respected at middle interior points à P3G
nGG3 ˆ M à P3G
nGG3 M S1 S2

…44†

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

à P2G
nGG2 ˆ M à P2G
nGG2 M S1 S3

…45†

@R ˆ0 ^ P1G †ij @ …M Sk

à P1G
nGG1 ˆ M à P1G
nGG1 M S2 S3

…46†

The above equations (equilibrium equations, boundary conditions for QÃ, continuity conditions within the element for Mà and QÃ) lead to a system with 39 equations and 39 unknowns. Due to the fact that the elementary external loading is already in equilibrium, the rank of this system is only 36. An optimization method described below is used to complete the system. 4.5. Stress resultant ®eld optimization The di€erence between the number of unknowns and the rank of the equations determining MÃ, Qà computation, leads to the introduction of three additional independent equations which allow a relative optimization. The error through the constitutive law being an upper bound of the exact error the system is solved in order to ®nd the stress resultant ®elds, corresponding to a minimal error on constitutive relation. No arbitrary constraint is introduced. Consequently, the stress resultants MÃ, Qà over an element are found as the solution of à † under the 36 independent à ,Q Minimize rO E…M conditions prescribed by Eqs: …35†±…46†

…47†

The minimization under constraints is obtained by a Lagrangian multiplier method. The unknown vector is now composed of two ®elds Mà and Qà and 36 Lagrangian multipliers. The following new function is introduced: à ,li † ˆ rO …M à †‡ à ,Q à ,Q R…M E

36 X li fi

…48†

iˆ1

where fi for i $ [1,. . . ,36] represents the scalar à TlTŠ à TQ equations given by (35)±(46). The unknown ‰M is the solution of the following problem @R ˆ0 à @M

@R ˆ0 à @Q

@R ˆ 0 8i 2 ‰1, . . . ,36Š @ li

G

^ †ij @ …M

ˆ 0 8i,j 2 ‰1,2Š  ‰1,2Š

…51†

8i,j,k,l 2 ‰1,2Š  ‰1,2Š  ‰1,2,3Š  ‰1,2Š @R

ˆ 0 8i,k 2 ‰1,2Š  ‰1,2,3Š

…52†

@R ˆ 0 8i,j,k 2 ‰1,2Š  ‰1,2Š  ‰1,2,3Š i @ …Q^ Sk †j

…53†

@R ˆ0 @ li

…54†

G

@ …Q^ Sk †i

8i 2 ‰1, . . . ,36Š

The resolution of this (78  78) system with a rank equal to 78 admits a single solution which renders the error through the constitutive relation minimal and leads to an optimal stress resultant ®eld over the element. Such a resolution for all the elements allows the SA stress ®eld construction over all the structure. The di€erent integration domains (global and elementary) will be used to de®ne either a local contribution or a global error. 4.5.1. Remarks . The great diversity among the SA computed stress resultant ®elds allows computation of some of them rather far from the exact solution, although keeping the SA property. Consequently, the optimization is of primary importance. . When considering plate equations, the optimization of the stress resultants within the element is natural. This is not the case for membrane (or 2D) cases where the classical approach leads to a unique stress ®eld [17]. . The equations presented in this section may seem somewhat complicated. In practice, the obtained system is rather particular (some equation are uncoupled). The main aspect that must be underlined again is that the calculations are made for each element, i.e. they are local which is very important in case analyses with large numbers of d.o.f. In these cases, the error calculation is not costly in comparison with the ®nite element computation.

…49†

All the partial derivatives of R have to be equal to zero. Eq. (9) can be written @R

623

…50†

5. Numerical examples In these examples, the shell element with linear interpolation presented in [23,24] is used for the ®nite element analyses. This element has been developed for

624

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

Fig. 2. Plate strip under constant pressure. (a) One clamped edge. (b) Simply supported. (c) Two clamped edges. (d) Type of meshes used.

non-linear analyses and especially for sheet forming simulations [29]. The analyses made with this element here, i.e. plate bending under geometrical linearity assumption, are not the best domain of application of this element, which uses linear interpolation functions. However, the goal of this section is not to test the eciency of the element but to show the quality of the error estimation based on the approach described above. In the same manner, it is also clear that the meshes used (regular size over the plate) are not opti-

Fig. 3. Error through the constitutive relation in comparison with the `exact error for a plate strip under constant pressure'. ECR: error through the constitutive relation. EXE: exact error.

mal, especially when the meshes are ®ne. The question here is to know whether the error estimation is in good agreement with an `exact error' obtained from an analytical solution. This solution is obtained from the Kirchho€ theory, although the analyses (and the ®nite elements) are based on the Reissner±Mindlin theory. The thicknesses of the plates considered in the di€erent examples are small enough in order to make these comparisons. In order to get an `exact error', the analytical solution is used in the error in energy norm, Eq. (7). In the di€erent examples, the error through the constitutive relation is compared to this `exact error', and the e€ective index, de®ned as the ratio (error

Fig. 4. E€ectivity index in cases of plate strips under constant pressure.

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

625

Fig. 5. Circular plate with a hole under constant pressure. (a) Simply supported exterior edge. (b) Clamped exterior edge. (c) Type of mesh used.

through the constitutive relation)/(`exact' error), is given for di€erent mesh re®nements. 5.1. Plate strip with di€erent boundary conditions Three di€erent boundary conditions are prescribed to a rectangular plate strip subjected to a constant normal pressure (Fig. 2). The meshes used are regular and cross diagonal. Fig. 3 shows, for di€erent mesh re®nements, that the estimated error through the constitutive relation is in very good agreement with the `exact' error pertaining to the three cases. The e€ectivity index (Fig. 4) is smaller than 1.1 and

becomes better as the mesh is re®ned. The property of the upper bound of the error through the constitutive relation in comparison to the `exact' error is veri®ed, and the estimated error is always larger than the `exact' error. 5.2. Circular plate with a hole The second analysis concerns a circular plate with a hole subjected to a normal pressure with two di€erent boundary conditions. The exterior edge is simply supported in the ®rst case and clamped in the second case (Fig. 5). The meshes used are shown in Fig. 5(c).

Fig. 6. Error on the constitutive relation in comparison with the `exact error' for a circular plate under constant pressure.

626

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

Fig. 7. E€ectivity index in the cases of a circular plate under constant pressure.

Fig. 10. E€ectivity index in the cases of a square plate under constant pressure.

Fig. 6 shows the correct agreement between the estimated error through the constitutive relation and the `exact' error concerning two cases. The e€ectivity index (Fig. 7) is good and smaller than 1.2 for the realistic meshes. The estimated error through the constitutive relation is always larger than the `exact' error. 5.3. Square plates

Fig. 8. Square plate under constant pressure. (a) Simply supported. (b) Clamped edges. (c) Type of meshes used.

The third analysis concerns a square plate with simply supported edges or clamped edges, and subjected to a constant pressure (Fig. 8). The meshes used are regular and cross diagonal (Fig. 8(c)). Fig. 9 shows the agreement between the estimated error through the constitutive relation and the `exact' error for the two cases. The e€ectivity index (Fig. 10) reveals a good behaviour of the error predictor (e€ectivity index <1.2 for realistic meshes). The estimated error through the constitutive relation is always larger than the `exact' error.

6. Conclusions and future prospects

Fig. 9. Error through the constitutive relation in comparison with the `exact error for a square plate under constant pressure'.

An error estimator based on the so-called `error through the constitutive relation' has been developed for Reissner±Mindlin plate bending analyses. The tests presented show the eciency of this estimator, and plate bending is a good ®eld of application for this method. The e€ectivity indices are satisfactory (<1.2) even for coarse meshes. It becomes very good and tends to 1 when the mesh is re®ned. The examples presented here con®rm in cases of plate bending analyses that the error through the constitutive relation is an upper bound of the exact error and the e€ectivity indices are always superior or equal to 1. This is very important from the practical engineering point of view. When an error on the constitutive law is estimated, the user knows that the result of his computation cannot be more di€erent from the exact solution. The plate

P. Boisse et al. / Computers and Structures 73 (1999) 615±627

bending analysis is a good ®eld of application of the error through the constitutive equation method because the nature of the di€erent equations is such that an optimization can be made within the statically admissible unknowns at di€erent steps of the calculation. This is the case while calculating bending moments and shear load densities on the element edges ensuring the element equilibrium and also when obtaining the stress resultant ®elds within the element. The applications of the proposed method to plate vibrations [30], Ahmad shell elements, and further to sheet metal forming simulations are now being studied [27]. References [1] Babuska I, Rheinboldt WC. A posteriori error estimates for the ®nite element method. Int J Numer Meth Engng 1978;12:1597±659. [2] Babuska I, Rheinboldt WC. Error estimates for adaptive ®nite element computations. SIAM J Numer Anal 1978;15(4):736±54. [3] Kelly DW, Cago JP, Zienkiewicz OC, Babuska I. A posteriori error analysis and adaptative processes in the ®nite element method. Part 1: error analysis. Int J Numer Meth Engng 1983;19:1593±619. [4] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptative procedure for practical engineering analysis. Int J Numer Meth Engng 1987;24:337±57. [5] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique. Part 2: Error estimates and adaptativity. Int J Numer Meth Engng 1992;33:207±24. [6] Wiberg NE, Abdulwahab F. A posteriori error estimation based on superconvergent derivatives and equilibrium. Int J Numer Meth Engng 1993;36:2703±24. [7] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery (SPR) and adaptative ®nite element re®nement. Comput Meth Appl Mech Engng 1992;101:207±24. [8] Zienkiewicz OC, Zhu JZ, Wu J. Superconvergent patch recovery techniques. Some further tests. Commun Numer Meth Engng 1993;9:251±8. [9] Rodriguez R, Some remarks on Zienkiewicz±Zhu estimators, Int J Numer Meth Engng (in press). [10] LadeveÁze P. Comparaison de modeÁles de milieux continus. TheÁse d'Etat, Universite P. et M. Curie, Paris 1975. [11] LadeveÁze P, Leguillon D. Error estimate procedure in the ®nite element method and applications. SIAM J Numer Anal 1983;20:485±509. [12] LadeveÁze P, Pelle JP, Rougeot P. Error estimation and mesh optimisation for classical ®nite elements. Engng Comput 1991;8:69±80. [13] Coorevits P, LadeveÁze P, Pelle JP. Mesh optimisation for problems with steep gradients. Engng Comput 1994;11:129±44. [14] Coorevits P, LadeveÁze P, Pelle JP. An automatic procedure with a control of accuracy for ®nite element analysis in 2D elasticity. Comput Meth Appl Mech Engng 1995;121:91±120. [15] Coorevits P, Dumeau JP, Pelle JP. Analyse eÂleÂments ®nis

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

627

adaptative pour les structures tridimensionnelles en eÂlasticiteÂ. Rev Eur eÂleÂments ®nis 1996;5:341±73. Coorevits P, Dumeau JP, Pelle JP. Adaptation et automatisation des maillages 3D. In: Actes du troisieÁme colloque national en calcul des structures, vol. 1. Presses acadeÂmiques de l'Ouest, 1997. p. 367±72. LadeveÁze P, Cognal G, Pelle JP. Accuracy of elastoplastic and dynamic analysis. In: Babuska I, Zienkiewicz OC, Cago JP, Oliovera de ER A, editors. Accuracy estimates and adaptativity for ®nite elements. New York: Wiley, 1986. p. 181±203. Gallimard L, LadeveÁze P, Pelle JP. Error estimation and adaptativity in elastoplasticity. Int J Numer Meth Engng 1996;39:189±217. Pelle JP, Ryckelynck D. Sur une strateÁgie adaptative pour les calculs non lineÂaires; application aux calculs en viscoplasticiteÂ. In: Actes du troisieÁme colloque national en calcul des structures, vol. 1. Presses acadeÂmiques de l'Ouest, 1997. p. 343±8. LadeveÁze P, Pelle JP. Accuracy in ®nite element computations for eigen frequencies. Int J Numer Meth Engng 1989;29:1929±49. Bussy P. Optimisation et ®abilite des calculs par eÂleÂments ®nis en non lineÂarite geÂomeÂtrique. TheÁse DocteurIngeÂnieur, Universite Paris VI 1984. Boisse P, Perrin S. Calcul de l'erreur en relation de comportement pour les modeÂlisations par des eÂleÂments ®nis de plaques C 0. C R Acad Sci Paris 1995;SeÂr. IIb 320:289±94. Boisse P, Daniel JL, Gelin JC. A simple isoparametric three node shell ®nite element. Comput Struct 1992;44:1263±74. Boisse P, Daniel JL, Gelin JC. A C 0 three node shell ®nite element for non-linear analyses. Int J Numer Meth Engng 1994;37:2239±364. Boisse P. EleÂments ®nis de Coques C 0 triangulaire aÁ interpolation lineÂaire pour les calculs de structures nonlineÂaires. In: Proceedings of 2eÁme Colloque National en Calcul des Structures, Giens, France, 1995. HermeÁs, Paris, 1995. 1. p. 299±305. Belytschko T, Stolarski H, Carpenter N. A C 0 triangular plate element with one-point quadrature. Int J Numer Meth Engng 1984;20:787±802. Hadjeb K. Analyse de la qualite des calculs de structures statiques et dynamiques par eÂleÂments ®nis de plaques et coques, Ph.D. Thesis, ENSAM Paris, France (in press). Watwood VB, Hartz BJ. An equilibrium stress ®eld model for ®nite element solutions of two dimensional elastoplastic problems. Int J Solid Struct 1968;4:857±73. Boisse P, Daniel JL, Gelin JC. Numerical simulations of 3D sheet metal forming employing new Ahmad shell elements. J Mater Process Technol 1992;34:117±24. Boisse P, Cognal G. Error through the constitutive relation for beam or C 0 plate ®nite elements: static and vibration analyses. In: LadeveÁze P, Oden JT, editors. Advances in adaptive computational methods in mechanics. Amsterdam: Elsevier, 1998. Cognal G. Fiabilite et optimisation des calculs eÂlastoplasiques par eÂleÂments ®nis. TheÁse d'eÂtat, Universite P. et M. Curie, Paris 1987.