Three-dimensional variational theory of laminated composite plates and its implementation with Bernstein basis functions

Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

www.elsevier.com/locate/cma

Three-dimensional variational theory of laminated composite plates and its implementation with Bernstein basis functions Alexander E. Bogdanovich * 3TEX, Inc., 109 MacKenan Drive, Cary, NC 27511, USA Accepted 26 March 1999

Abstract A new three-dimensional (3-D) variational theory aimed at the stress analysis of thick laminated rectangular plates with anisotropic layers is presented. The developed theory is versatile, it allows one to use various types of basis functions and apply them independently to each of the three coordinate directions. Also, the theory is material-adaptive, enabling us to impose the conditions of displacement continuity between all of the discretization elements (3-D bricks) within the body and, in addition, the conditions of stress continuity between adjacent bricks made from the same material. The theory can be applied, in principle, to any boundary value problem of linear elasticity, considering arbitrarily distributed static surface forces and external kinematic boundary conditions. One speci®c realization of the theory, elaborated in this work, applies Bernstein basis functions of an arbitrary degree for the displacement approximation in the three coordinate directions. This type of basis function, which are a rarity in the structural analysis, deserve more attention. They provide computational ef®ciency and unique analytical elegance to the mathematical algorithms. The developed theory can be also used for generating new types of higher-order hexahedral ®nite element, as illustrated here on the example of a novel 8-node ``Bernstein ®nite element''. Numerical examples given in this work show that, following the concept of material-adaptive inter-element continuity conditions, a higher smoothness and accuracy of the computed stresses can be achieved by incorporating additional inter-element constraints between the same material bricks. However, it is also easy to spoil the solution by adding the same type constraints between adjacent bricks having distinct material properties. Finally, numerical comparison for the benchmark 3-D problem of transverse bending of simply supported 3-layer cross-ply laminated plate validates that the present analysis is signi®cantly more accurate and computationally ef®cient than ANSYS SOLID 46 element. Ó 2000 Elsevier Science S.A. All rights reserved.

1. Introduction Development of three-dimensional (3-D) solid ®nite elements, since the mid-1960s, has been one of the major directions of the ®nite element methodology. A variety of 3-D elements were derived based on the earlier 1-D and 2-D ®nite element approaches. The most noticeable types of solid ®nite elements are: (i) Tetrahedral elements developed by Gallagher et al. [1], Melosh [2] and Argyris [3±5], and (ii) Hexahedral elements which are of a particular interest in the context of this work. Various elements of the latter type have been developed by Argyris et al. [6±8], Ergatoudis et al. [9±11] and other authors [12±16]. Among the most popular elements are: 8-node (linear displacement) element with 24 degrees of freedom (d.o.f.), 20node (quadratic displacement) element with 60 d.o.f., 32-node (incomplete cubic displacement) element with 96 d.o.f., 64-node (complete cubic displacement) element with 192 d.o.f. (the isoparametric version of this element has been termed LUMINA [6]), 8-node (incomplete quintic displacement ®elds) element with 96 d.o.f. and Hermite element with both displacements and displacement derivatives considered as

*

Tel.: +1-919-755-1818; fax: +1-919-755-1914. E-mail address: [email protected] (A.E. Bogdanovich).

0045-7825/00/$ see front matter Ó 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 2 6 3 - 7

280

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

independent d.o.f. A detailed description of the aforementioned, as well as other particular 3-D solid ®nite elements, can be found in [17±20] and other books on ®nite element analysis. In spite of a considerable amount of analytical work devoted to 3-D solid ®nite elements and their utilization in commercial ®nite element packages like MSC NASTRAN, ABAQUS, ANSYS, DYNA 3D, ANIDA, etc., their practical applications are still far less common than the use of 1-D and 2-D elements. This is especially obvious in the case of composite structures. The major arguments against 3-D elements are: (i) computer RAM and CPU time required for a full 3-D analysis are incomparably higher than for the respective 1-D and 2-D analysis and (ii) the necessity and bene®ts of a full 3-D structural analysis are often not obvious. Owing to the rapidly increasing a€ordability of powerful personal computers, the ®rst argument will be obsolete very soon while the second one needs convincing demonstrations that in many practical structural analysis problems it is imperative using 3-D stress and failure analysis. Yet, most of the applications of 3-D ®nite elements are still related to the analysis and design of bulk metal and concrete structures. A few publications are known on 3-D analysis of laminated composite plates and shells (see reviews in [21±25,27]). It has been well recognized since the early works (examples are [7,8,11,13,27]) that one of the primary challenges of the ®nite element methodology is to correctly formulate compatibility of the elements. As noted in [27] ``The problems of discrete element analysis which have been most vexing are: (1) the identi®cation of geometric admissibility conditions for the element and (2) generation of local displacement patterns which will satisfy these conditions. The elements must be based upon modes that cause them to be geometrically compatible with adjacent elements along their entire common boundary, and satisfy everywhere the geometric boundary conditions.'' Thus, the ®rst requirement to the shape functions is that they should assure a C0 class inter-element continuity of the displacements. Further, in the case of homogeneous materials, where there are no physical interfaces, enforcing inter-element continuity of the corresponding ®rst derivatives of the shape functions allows one to satisfy the strain (in this case, also the stress) continuity conditions. For example, the authors of [8] have stated that ``If physically admissible, HERMES 8 hexahedral element allows one to stipulate continuity of strains and stresses.'' At the same time it has been recognized that ``in case of discontinuity of the material properties of the elements meeting at a nodal point, the imposition of common strains at such points involves, in general, additional constraints. To relax these one may either assign di€erent derivatives to the element in question or introduce stresses and rotations as freedoms.'' When characterizing HERMES 8 element in [13] it has been emphasized that ``this involves assumption of derivative continuity, which is not true in nonhomogeneous situations''. The above methodological remarks are especially important in the context of the development of numerical analysis methods for structurally inhomogeneous materials and, particularly, laminated composites. Obviously, in the case of laminated structures, which are characterized by a step-wise material property variation, it is impossible to simultaneously satisfy the interfacial strain and stress continuity. From the mechanics standpoint, obtaining continuous interfacial stresses (which would assure equilibrium of the interface) should be a correct aim, which is, however, automatically deprived if the inter-laminar strain continuity is enforced. Nevertheless, some ®nite element developments for laminated structures (the examples are [28] and [29]) have missed the point and enforced inter-element strain continuity at the interfaces between distinct materials. A discussion of this issue can be found in [24,25]. Illustrative example showing that the strain continuity conditions enforced at the interfaces between distinct layers totally spoil the stress solution was presented in [30]. Based on the above considerations, it has been suggested in [24,25,31] that, if wishing to develop a more accurate and computationally ecient analysis tools for composite structures, a material-adaptive concept should be considered. Speci®cally, in the context of developing innovative ®nite element analysis approaches, this means that properties of the element shape functions have to be tied to the material properties. From a glance, such an idea may seem a little unrealistic, though our earlier analytical developments and numerical results presented in [31±34] show that this is analytically and computationally viable. A utilization of 3-D mosaic body model introduced in [31] allows one to enforce different types of connectivity between the discretization elements (3-D bricks), depending on their material properties. In the case of laminated composite structures this means that the element shape functions should be different for the in-plane coordinate directions on one side, and through thickness coordinate direction on the other side.

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

281

Ideally, the ``in-plane'' shape functions should possess continuous derivatives up to the third-order (an explanation for that can be found in classical text [35]; see also Section 1.5 of [36]) while the ``transverse'' shape functions must have discontinuous ®rst derivatives at all of the interfaces within the laminated structure. This issue has been discussed in [24] and [25]. In the cases of a more complex type of the composite material inhomogeneity, where there are property variations in the two or three coordinate directions, the 3-D mosaic body concept [31] seems to be an adequate one. Present work can be viewed as further development of 3-D variational analysis approach, di€erent versions of which have been reported in our earlier publications, [24,25,30,31,36]. A novel type of 3-D displacement approximation, using products of three independent sets of Bernstein basis functions, is applied here. Various analytical advantages of this novel approach are demonstrated. A ®nite element interpretation of the developed variational analysis is given in terms of a 3-D higher-order hexahedral ®nite element. Numerical examples of 3-D stress analysis of rectangular homogeneous and laminated composite bars and plates illustrate the applications. Solving the benchmark problem [37] validates the accuracy of stress predictions. It is also shown that present analysis is signi®cantly more accurate and computationally ecient than a popular ANSYS SOLID 46 ®nite element analysis, which has been applied to the same problem in [36].

2. Three-dimensional variational theory of laminated plates Consider a rectangular laminated plate consisting of some number of anisotropic layers each characterized by its individual 21 elastic constants. As shown in Fig. 1, the plate is discretized by three sets of orthogonal planes x ˆ xl , y ˆ ym , z ˆ zn …l ˆ 2; . . . ; L, m ˆ 2; . . . ; M, n ˆ 2; . . . ; N † into LMN parallelepipeds (further called ``bricks''). In the case under consideration, all bricks belonging to the same layer have identical properties. Thus, the entire plate is represented as a brick-type mosaic structure composed from linearly elastic anisotropic solids. The stress±strain relations are written in the following concise form: …s†

…s† …s†

ri …r† ˆ Cij ej …r†; …s†

i; j ˆ 1; . . . ; 6; …s†

…1† …s†

where ri are stresses and ej strains in the sth brick, Cij are components of its sti€ness matrix, and r ˆ fx; y; zg is a position vector. In this theory the strains are related to the displacements u1 …x; y; z†, u2 …x; y; z† and u3 …x; y; z†, corresponding to the x-, y- and z-directions, through the equations of linear 3-D elasticity

Fig. 1. Three-dimensional mosaic model of a rectangular laminated plate.

282

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304 …s†

…s†

e1 ˆ

ou1 ; ox

…s†

…s†

e2 ˆ

ou2 ; oy

…s†

e3 ˆ

…s†

ou3 ; oz

…s†

e4 ˆ

…s†

…s†

ou2 ou ‡ 3 ; oz oy

…s†

…s†

e5 ˆ

…s†

ou1 ou ‡ 3 ; oz ox

…s†

e6 ˆ

…s†

…s†

ou1 ou ‡ 2 : oy ox …2†

The principle of minimum total potential energy is implemented for the entire plate dP ˆ d…P ÿ W † ˆ 0:

…3†

Here, P is strain energy and W is work of external surface forces acting on the plate. The total potential energy P is further represented as the sum of potential energies P…s† of all bricks in the structure Pˆ

S X

P…s† ˆ

S X ÿ

sˆ1


P …s† ÿ W …s† ;

…4†

sˆ1

where S ˆ LMN . In order to evaluate expression (4), it is necessary to use the stress±strain equation (1), the strain-displacement equation (2), and also to assume some form of the displacement approximation in all three coordinate directions. Here we adopt the following form of the displacement approximation for the sth brick: u…s† a …x; y; z† ˆ

I X J X K X iˆ0

jˆ0

kˆ0

a…s†

Uijk Xi …x†Yj …y†Zk …z†;

a ˆ 1; 2; 3;

…5† a…s†

where integer values I, J and K de®ne number of d.o.f. in the model; Uijk are undetermined coef®cients; Xi …x†, Yj …y† and Zk …z† are three sets of basis functions. At this point we do not specify them and perform further derivations in a general form. Note that total number of d.o.f. in the model is 3S…I ‡ 1†…J ‡ 1†…K ‡ 1†. Having the displacement approximation de®ned by Eq. (5) one can further proceed in evaluating P …s† and …s† W entering in Eq. (4). One important algorithmic step is to de®ne a kinematically admissible displacement ®eld. This should allow one to eliminate those undetermined coef®cients from (5), which become known due to the prescribed external kinematic (geometric) boundary conditions. This procedure is, generally similar to the ``condensation'' in the conventional ®nite element analysis. The ``reduced'' set of undetera…s† mined coef®cients, remaining in the displacement approximation, is denoted U~ijk . After the elimination has been accomplished, the running indices i, j and k can be rearranged for the convenience, in order to maintain their count from zero. The upper summation limits in (5) are now some different numbers, say Ia…s† 6 I, Ja…s† 6 J , Ka…s† 6 K …a ˆ 1; 2; 3†. Accordingly, the modi®ed displacement approximation is now represented in the form …s†

u…s† a …x; y; z†

ˆ

…s†

…s†

Ia Ja Ka X X X iˆ0

jˆ0

kˆ0

a…s† e ijk U Xi …x†Yj …y†Zk …z† ‡

I X

J X

…s†

…s†

K X

a…s† U^ijk Xi …x†Yj …y†Zk …z†;

…6†

…s†

iˆIa ‡1 jˆJa ‡1 kˆKa ‡1

where U^ijk are already known from the external kinematic boundary conditions. Accordingly, total number of d.o.f. in the model is reduced to a…s†

Dˆ

S X 3 X ÿ …s†
ÿ
ÿ
Ia ‡ 1 Ja…s† ‡ 1 Ka…s† ‡ 1 : sˆ1

…7†

aˆ1

Importantly, the modi®ed displacement approximation (6) may contain di€erent number of undetermined coecients for di€erent bricks while the primary approximation equation (5) had assumed equal number of them for all of the bricks. Using the primary displacement approximation (5), the total potential energy function (4) is written in the form Pˆ

S X sˆ1

S h   X    i a…s† a…s† a…s† P…s† Uijk ˆ P …s† Uijk ÿ W …s† Uijk ; sˆ1

…8†

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

283

where it is emphasized by the arguments in parentheses that both the strain energy and work of external surface forces are fully de®ned in terms of the displacement approximation coecients. The following expressions of P …s† and W …s† are then obtained: T 1 P …s† ˆ U…s† A…s† U…s† ; 2 T …s† W ˆ U…s† Q…s† ;

where

U…s†

…9† …10†

o9 8n 1…s† > > U > > ijk > > > = 2…s† ; ˆ Uijk > > > > n o > > > ; : U 3…s† > ijk

Q…s†

o9 8n 1…s† > > Q > > ijk > > > = 2…s† ; ˆ Qijk > > > > n o > > > ; : Q3…s† > ijk

ab…s†

2h

ih i3 13…s† aijk;pqr 6h 7 6 21…s† i h 22…s† i h 23…s† i 7 7 ˆ6 a a a ijk;pqr ijk;pqr 7: 6 ijk;pqr 4h ih ih i5 31…s† 32…s† 33…s† aijk;pqr aijk;pqr aijk;pqr

A…s†

11…s†

aijk;pqr

ih

12…s†

aijk;pqr

…11†

a…s†

Explicit expressions of aijk;pqr and Qijk are given in Appendices A and B. To accomplish the condensation procedure, vectors U…s† , Q…s† and matrix A…s† in (11) are substructured as following: ( ) ( ) " # …s† U…s† Q…s† A…s† Aab a a …s† …s† aa …s† ; Q ˆ ; A ˆ : …12† U ˆ …s† …s† …s† …s† Ub Qb Aba Abb …s† Here, U…s† a and Qa correspond to

i ˆ 0; 1; . . . ; Ia…s† ; …s†

j ˆ 0; 1; . . . ; Ja…s† ;

k ˆ 0; 1; . . . ; Ka…s† ;

…13†

…s†

while Ub and Qb correspond to i ˆ Ia…s† ‡ 1; . . . ; I; …s†

j ˆ Ja…s† ‡ 1; . . . ; J ;

…s†

k ˆ Ka…s† ‡ 1; . . . ; K:

…14†

…s†

Matrices A…s† aa , Aab , Aba and Abb are de®ned for the following index combinations: A…s† aa …s†

Aab

for i; p ˆ 0; 1; . . . ; Ia…s† ; for i ˆ 0; 1; . . . ; Ia…s† ;

q ˆ Ja…s† ‡ 1; . . . ; J ; …s†

Aba

j; q ˆ 0; 1; . . . ; Ja…s† ; p ˆ Ia…s† ‡ 1; . . . ; I;

k ˆ 0; 1; . . . ; Ka…s† ;

for i ˆ Ia…s† ‡ 1; . . . ; I;

…s†

Abb

r ˆ 0; 1; . . . ; Ka…s† ;

The following expressions are then obtained from (9) and (10) with the use of (12):  T 1  …s†T …s† …s† …s† …s† …s†T …s† …s† …s†T …s† …s† Ua Aaa Ua ‡ U…s† ; P …s† ˆ a Aab Ub ‡ Ub Aba Ua ‡ Ub Abb Ub 2 …s†

T

…s†

…s† W …s† ˆ U…s† a Qa ‡ Ub Qb :

…16†

j ˆ Ja…s† ‡ 1; . . . ; J ;

for i; p ˆ Ia…s† ‡ 1; . . . ; I; j; q ˆ Ja…s† ‡ 1; . . . ; J ; k; r ˆ Ka…s† ‡ 1; . . . ; K:

T

…15†

j ˆ 0; 1; . . . ; Ja…s† ;

r ˆ Ka…s† ‡ 1; . . . ; K;

p ˆ 0; 1; . . . ; Ia…s† ;

q ˆ 0; 1; . . . ; Ja…s† ; k ˆ Ka…s† ‡ 1; . . . ; K;

k; r ˆ 0; 1; . . . ; Ka…s† ;

…17† …18†

…19† …20†

It is worth mentioning that the ®rst term in (19) is a quadratic form of the undetermined coecients, the second and third terms are their linear forms, and the fourth term is a constant. The ®rst term in (20) is a linear form of the undetermined coecients, and the second term is a constant. It is known that one of the consequences of the variational principle (3) are Ritz-type equations which are, in our case, written in the following form: ! S X oP …s† oW …s† ÿ ˆ 0: …21† oU…s† oU…s† sˆ1 a a

284

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

Here, P …s† and W …s† are de®ned by (19) and (20), respectively. Before imposing connectivity of the bricks, P (21) provides S separate systems of linear simultaneous equations, each containing 3aˆ1 …Ia…s† ‡ 1†…Ja…s† ‡ 1† …Ka…s† ‡ 1† unknown displacement approximation coef®cients. Matrix form of the sth representative system, derived from Eqs. (19)±(21) with account for the symmetry of matrix A…s† (this property is proven in Appendix A), is …s†

…s†

…s† …s† A…s† aa Ua ˆ Qa ÿ Aab Ub :

…22†

Note that the right-hand side of the above equation contains both the external surface forces and external displacements applied to the brick. The above theory allows one to solve any 3-D elasticity problem for an anisotropic brick with imposed external displacements, surface forces or any appropriate combination of thereof. It has to be pointed out that no a priori assumptions have yet been made regarding speci®c type of basis functions in the displacement approximation (5) or regarding their number, which is determined by the indices I, J, and K. In particular, the following types of basis functions could be used: power series, Lagrange, Hermite, Tschebyscheff, Legendre, and Bernstein polynomials, B-splines, trigonometric functions, etc. Importantly, this theory is allowing choosing them independently for each of the coordinate directions. For example, a combination of trigonometric functions in the x- and y-directions along with the Bernstein polynomials in the z-direction have been utilized in [30]. Speci®cs of the solution related to the basis functions will be mainly exposed in the integrals given in Appendices A and B. A particular choice of basis functions can be ruled by their preferable order of continuity, analytical convenience, higher convergence rate, computational stability, and other considerations. Further development of the theory requires formulating a full set of internal kinematic and static boundary conditions between the bricks i.e., establishing their connectivity. 3. Connectivity of 3-D bricks When using nomenclature of the bricks introduced in Fig. 1, the full set of displacement continuity conditions in the mosaic body can be formulated as follows. For the inter-brick boundaries perpendicular to the x, y and z axis, respectively, it is obtained uas…l;m;n† …xl ; y; z† ˆ uas…l‡1;m;n† …xl ; y; z†;

…23†

uas…l;m;n† …x; ym ; z† ˆ uas…l;m‡1;n† …x; ym ; z†; uas…l;m;n† …x; y; zn † ˆ uas…l;m;n‡1† …x; y; zn †;

…24† …25†

where a ˆ 1; 2; 3 in all the above equations, and s…l; m; n† ˆ 1 ‡ …l ÿ 2† ‡ …m ÿ 2†L ‡ …n ÿ 2†LM:

…26†

The following values of l, m and n take place l ˆ 2; . . . ; L;

m ˆ 2; . . . ; M ‡ 1;

n ˆ 2; . . . ; N ‡ 1;

…27†

n ˆ 2; . . . ; N ‡ 1;

…28†

in Eq. (23), l ˆ 2; . . . ; L ‡ 1;

m ˆ 2; . . . ; M;

in Eq. (24), l ˆ 2; . . . ; L ‡ 1;

m ˆ 2; . . . ; M ‡ 1;

n ˆ 2; . . . ; N ;

…29†

in Eq. (25). Further, considering continuity of transverse stresses between adjacent bricks, one obtains the following relations (double-index notations r11 ˆ r1 , r22 ˆ r2 , r33 ˆ r3 , r23 ˆ r4 , r31 ˆ r5 , r12 ˆ r6 are used here for conciseness). For the inter-brick boundaries perpendicular to the x axis: s…l;m;n†

r1a

s…l‡1;m;n†

…xl ; y; z† ˆ r1a

…xl ; y; z†:

…30†

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

285

For the inter-brick boundaries perpendicular to the y axis: s…l;m;n†

r2a

s…l;m‡1;n†

…x; ym ; z† ˆ r2a

…x; ym ; z†:

…31†

And for the inter-brick boundaries perpendicular to the z axis: s…l;m;n†

r3a

s…l;m;n‡1†

…x; y; zn † ˆ r3a

…x; y; zn †:

…32†

As before, a ˆ 1; 2; 3. Values of s are de®ned by Eq. (26); values of l, m and n de®ned by Eqs. (27)±(29), apply to Eqs. (30)±(32), respectively. Eqs. (23)±(32) cover all possible cases of the displacement and stress continuity conditions between the bricks in a mosaic body shown in Fig. 1. Next step of the solution procedure is to express these conditions in terms of the displacement approximation coecients in Eq. (5). However, in order to accomplish this step one has to fully specify all basis functions in (5). Section 4 illustrates one particular implementation of the developed theory.

4. Implementation of the theory with Bernstein basis functions Theory of Bernstein polynomials is well developed and widely used in several areas of applied mathematics. The applications are common in computational geometry, computer graphics, computer-aided geometric design (see, for example [38±40]), also in probability theory (examples can be found in [41]), but are very rare in computational mechanics. The author's opinion is that this type of approximation polynomials could be of a great value in advanced structural analysis. In particular, it has been shown in [30] that the use of Bernstein polynomials for the through thickness approximations of displacements and stresses in the displacement-assumed and mixed variational analyses of laminated composite plates provides uniquely accurate results. Bernstein basis functions of degree N are de®ned as BNn …t† ˆ

N! N ÿn tn …1 ÿ t† ; n!…N ÿ n†!

t 2 ‰0; 1Š;

n ˆ 0; 1; . . . ; N ;

…33†

where N may be any integer number. For any value of n, Eq. (33) de®nes the N th degree polynomial. It is well-known that all algebraic polynomials of degree 6 N constitute a linear space of dimension N ‡ 1. Since the N ‡ 1 functions BNn …t†, where n ˆ 0; 1; . . . ; N are linearly independent, they form a basis. Bernstein basis functions posses certain valuable properties, which make them suitable for developing rather simple and elegant mathematical solutions and computational algorithms for various boundary value problems of elasticity. One of their properties, which could be of a particular interest in the structural analysis applications, is that N X nˆ0

BNn …t† ˆ 1 for any t 2 ‰0; 1Š:

…34†

This property is perfectly suited for one important class of boundary conditions, which requires satisfying a uniformly distributed displacement along the exterior surface or its part. Another valuable property of Bernstein polynomials is that in contrast to the other types of approximation polynomials (Tschebysche€ polynomials, for example), they yield smooth approximants [42]. This means that they provide simultaneous approximation of the function and its derivatives. Now, considering sth brick located between the planes xl and xl‡1 , ym and ym‡1 , zn and zn‡1 (see, Fig. 1), basis functions entering in the displacement approximation (5) are de®ned as Xi …x† ˆ

8 < BI …x† ˆ i

:

0

   x ÿ xl i xl‡1 ÿ x Iÿi I! xl‡1 ÿ xl i!…I ÿ i†! xl‡1 ÿ xl

for x 2 ‰xl ; xl‡1 Š; else;

…35†

286

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

Yj …y† ˆ

Zk …z† ˆ

8 < BJ …y† ˆ

   y ÿ ym j ym‡1 ÿ y J ÿj J! ym‡1 ÿ ym j!…J ÿ j†! ym‡1 ÿ ym

j

:

0

8 < BK …z† ˆ k

:

   z ÿ zn k zn‡1 ÿ z Kÿk K! zn‡1 ÿ zn k!…K ÿ k†! zn‡1 ÿ zn

0

for y 2 ‰ym ; ym‡1 Š;

…36†

else; for z 2 ‰zn ; zn‡1 Š;

…37†

else:

It is important to note that, with the use of basis functions (35)±(37), all integrals introduced in Appendix A, see Eq. (A.2), can be evaluated analytically. Their analytical expressions are presented in Appendix C. Thus, a numerical integration (usually, Gaussian quadrature) can be avoided. Actually, the Gauss points, which are of some signi®cance in conventional ®nite element approaches, will never appear in the present analysis approach. Also, in some speci®c cases of the surface load distributions, is would be possible performing an analytical evaluation of the integrals introduced in Appendix B, Eqs. (B.4). It is easy to show that with the use of property (34), it is possible to exactly satisfy the uniformly distributed displacement boundary condition. Speci®cally, let us assume that uniformly distributed displacement u01 is imposed along the entire exterior surface x ˆ x1 (see, Fig. 1), so the boundary condition u1 …x1 ; y; z† ˆ u01 should be identically satis®ed for all y and z values. Now, using Eqs. (5) and (35) one obtains I X J X K X iˆ0

jˆ0

kˆ0

1…s†

Uijk BIi …x1 †BJj …y†BKk …z† ˆ u01 :

…38†

This equation has to be satis®ed for all those bricks which are exposed to the surface x ˆ x1 , and for all values of the coordinates y and z. Now, using Eq. (35), one obtains that the only nonzero basis function at x ˆ x1 is BI0 …x1 † ˆ 1. Hence, Eq. (38) yields J X K X jˆ0

kˆ0

1…s†

U0jk BJj …y†BKk …z† ˆ u01 :

…39†

Taking into account (34), the above equation is further written in the following form: J X K   X 1…s† U0jk ÿ u01 BJj …y†BKk …z† ˆ 0: jˆ0

…40†

kˆ0

Now, due to the linear independence of basis functions, the above equation is reduced to the following simple relation: 1…s†

U0jk ˆ u01 :

…41†

This has to be satis®ed for j ˆ 0; 1; . . . ; J and k ˆ 0; 1; . . . ; K and all the respective values of s. In a similar way, the uniformly distributed displacement boundary condition can be formulated for any other surface element of the mosaic body. Next step is to express the displacement continuity conditions (23)±(25) in terms of the undetermined coecients in (5). By substituting (5) together with (35)±(37) in (23)±(25), it is obtained J X K h X jˆ0

kˆ0

I X K h X iˆ0

kˆ0

I X J h X iˆ0

jˆ0

a;s…l;m;n†

ÿ U0jk

a;s…l;m;n†

ÿ Ui0k

a;s…l;m;n†

ÿ Uij0

UIjk

UiJk

UijK

a;s…l‡1;m;n†

a;s…l;m‡1;n†

a;s…l;m;n‡1†

i

i i

BJj …y†BKk …z† ˆ 0; BIi …x†BKk …z† ˆ 0; BIi …x†BJj …y† ˆ 0:

…42†

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

287

Again, using linear independence of basis functions, these equations are reduced to the following: a;s…l;m;n†

UIjk

a;s…l‡1;m;n†

ˆ U0jk

;

a;s…l;m;n†

UiJk

a;s…l;m‡1;n†

ˆ Ui0k

;

a;s…l;m;n†

UijK

a;s…l;m;n‡1†

ˆ Uij0

:

…43†

Each of the above three groups of equations should be satis®ed for their respective values of i, j and k. Values of s, l, m and n in (43) are de®ned by (26)±(29), respectively. If all of the Eq. (43) are satis®ed, then the solution of any boundary value problem will provide continuous displacement ®eld within the entire mosaic body, regardless of the material properties of individual bricks, external kinematic and static boundary conditions. Next, consider conditions of the transverse stress continuity between adjacent bricks in the mosaic body. Those are presented by Eqs. (30)±(32). The corresponding equations, formulated through the undetermined coecients of the displacement approximation (5), are derived in Appendix D. It seems virtually impossible to satisfying those conditions exactly between the bricks having distinct material properties. Yet, it is shown in Appendix D that the strain and stress continuity conditions can be exactly satis®ed between the same material bricks. These conditions are of the following form …a ˆ 1; 2; 3 in all of these equations): xl‡1 ÿ xlÿ1 a;s…l;m;n† xl‡1 ÿ xl a;s…l;m;n† a;s…l‡1;m;n† ˆ U ÿ U ; j ˆ 0; 1; . . . J ; k ˆ 0; 1; . . . ; K; …44† U1jk xl ÿ xlÿ1 Ijk xl ÿ xlÿ1 Iÿ1jk ym‡1 ÿ ymÿ1 a;s…l;m;n† ym‡1 ÿ ym a;s…l;m;n† a;s…l;m‡1;n† ˆ U ÿ U ; i ˆ 0; 1; . . . I; k ˆ 0; 1; . . . ; K; …45† Ui1k ym ÿ ymÿ1 iJk ym ÿ ymÿ1 iJ ÿ1k zn‡1 ÿ znÿ1 a;s…l;m;n† zn‡1 ÿ zn a;s…l;m;n† a;s…l;m; n‡1† ˆ U ÿ U ; i ˆ 0; 1; . . . I; j ˆ 0; 1; . . . ; J : …46† Uij1 zn ÿ znÿ1 ijK zn ÿ znÿ1 ijKÿ1 In the derivation of Eqs. (44)±(46) it has been assumed that each layer in the laminated plate is a homogeneous 3-D elastic body, so the stress and strain continuity conditions are identical within the layer. Obviously, this does not apply to the stress and strain continuity conditions between distinct layers. Accordingly, Eq. (46) does not apply for those n values, which correspond to the interfaces between distinct layers. Symbolically, this is indicated in Eq. (46) by using index n instead of n. Thus, in the present analysis the transverse stress continuity conditions are not enforced between layers, but are enforced at all of the other inter-brick boundaries. This feature of the developed analysis approach supports the statement that this analysis is a material-adaptive. Note that if, despite the above warning, one would apply Eq. (46) to the layer interfaces (this means enforcing the transverse strain continuity between the layers), the results for stresses would be absolutely wrong. This kind of numerical experiment has been exercised in [30], and the issue will be also illustrated later in this work. It should be added that, though the transverse stress continuity conditions between the layers are not enforced in the present analysis, those are expected to be accurately satis®ed in a ``soft'' sense, as the natural internal static boundary conditions. After a full set of the connectivity equations for the mosaic body has been expressed in terms of the displacement approximation coecients, it remains to incorporate Eqs. (43)±(46) in the previously obtained separate systems of linear simultaneous Eqs. (22), as to couple the systems. Details of this procedure, which is very tedious, are omitted here. It can only be mentioned that, after this step is accomplished, the system matrix still retains some distinct blocks, though their inter-connection pattern is very complex. Also, it is worth mentioning that total number of equations in the coupled system is, usually, considerably lower than total number of equations in S uncoupled systems. Numerical solutions of the aforementioned, coupled system of linear simultaneous equations may be performed using routine procedures. The result of numerical solution provides some part of the disa…s† placement approximation coecients U~ijk . The other part of them is then obtained from Eqs. (43)±(46). a…s† ^ Recall that Uijk have been already de®ned from the imposed external kinematic boundary conditions. Thus, the displacement ®eld, which is represented in the form of triple series (6), is now fully de®ned. In order to perform summation of the series, one needs to prescribe speci®c values of the x, y and z coordinates. The summation provides displacement values at the corresponding point within the plate. Further, the strains are directly obtained by substituting (6) in Eq. (2). The respective derivatives of the basis

288

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

functions are analytically evaluated from Eqs. (35)±(37), with the use of formula (C.3) from Appendix C. Finally, the stresses are obtained directly from the constitutive equation (1). Note that the whole solution algorithm utilizes only one approximate numerical procedure, namely solution of the coupled system of linear simultaneous equations. There are no extrapolation, smoothing, or additional postprocessing procedures (like utilization of the equilibrium equations for the transverse stress evaluation) involved in this analysis. 5. Bernstein family of 8-noded elements The developed variational analysis approach with Bernstein basis functions used in the displacement approximation (5), can be reformulated in terms of the ®nite element analysis. Indeed, the primary displacement approximation (5) can be represented in the following equivalent form: Ua…s† …x; y; z† ˆ

F X f ˆ1

a…s†

Uf Nf …x; y; z†;

a ˆ 1; 2; 3;

…47†

where Nf …x; y; z† ˆ Nijk …x; y; z† ˆ BIi …x†BJj …y†BKk …z†;

…48†

f ˆ f …i; j; k† ˆ …k ‡ 1† ‡ j…K ‡ 1† ‡ i…J ‡ 1†…K ‡ 1†;

…49†

F ˆ …I ‡ 1†…J ‡ 1†…K ‡ 1†:

…50†

Now it is recognized that Bernstein basis functions (35)±(37) provide the following properties to the functions Nf …x; y; z†: N1 …xl ; ym ; zn † ˆ 1; NK‡1 …xl ; ym ; zn‡1 † ˆ 1; N1‡J …K‡1† …xl ; ym‡1 ; zn † ˆ 1; N…J ‡1†…K‡1† …xl ; ym‡1 ; zn‡1 † ˆ 1; N1‡I…J ‡1†…K‡1† …xl‡1 ; ym ; zn † ˆ 1;

…51†

N…K‡1†‰1‡I…J ‡1†Š …xl‡1 ; ym ; zn‡1 † ˆ 1; N1‡…K‡1†‰J ‡I…J ‡1†Š …xl‡1 ; ym‡1 ; zn † ˆ 1; N…I‡1†…J ‡1†…K‡1† …xl‡1 ; ym‡1 ; zn‡1 † ˆ 1: For each of the above eight combinations of l, m and n (geometrically they represent eight corner points of the brick: A, B, C, D, E, F, G and H, as shown in Fig. 2) all functions Nf , except those listed in Eq. (51), take zero values. Thus, Bernstein basis functions (35)±(37) can be viewed as speci®c type of ``shape functions'', which de®ne the eight-node hexahedral ®nite element. Though, the element is quite unusual: it may have only eight nodes (which are always its corners), independently of the degrees I,J and K. It is easy to verify that there are no other points in the brick where any of the functions (48) take zero value. Therefore, no edge nodes, side nodes or internal nodes (which are typical, for example, with higher-order Lagrange ®nite element), can exist in this ``Bernstein element''. An explanation of this fundamental difference is simple: Lagrange polynomials are interpolation polynomials while Bernstein polynomials are approximation polynomials. Nevertheless, it is possible to directly compare the number of d.o.f. associated with the ``interpolation'' and ``approximation'' hexahedral ®nite elements A comparison of their computational effectiveness and accuracy is an interesting issue for the future study. However, one advantage of this novel family of hexahedral Bernstein ®nite elements has been already shown in this work: their analytical elegance. Numerical examples illustrating applications of the developed variational approach are presented in the next two sections.

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

289

Fig. 2. An eight-node ®nite element.

6. 3-D analysis of laminated bar under uniaxial extension Consider a laminated rectangular bar composed of an arbitrary number of orthotropic layers perpendicular to the x-axis, as shown in Fig. 3. The layers may have identical or distinct properties. Dimensions of the bar in the x-, y- and z-directions are a, b and c, respectively. In all numerical studies performed here, a single-brick was used in the y±z cross-section. When relating the problem under consideration to the mosaic body model shown in Fig. 1, it is seen that x1 ˆ 0, y1 ˆ 0, z1 ˆ 0, xL ˆ a, y2 ˆ b, z2 ˆ c. In this section we will use subscripts x, y and z in relation with the displacement, strains and stress components, instead of the subscripts 1, 2 and 3, which were used in the previous sections. Thus, the very left brick is fully ®xed at its left end, x ˆ 0, and the very right brick is exposed to the uniformly distributed displacement u0x at its right end x ˆ a. The corresponding boundary conditions are formulated as follows ux…1† …y; z† ˆ 0;

uy…1† …y; z† ˆ 0;

0 u…L† x …y; z† ˆ ux

at x ˆ xL‡1 :

u…1† at x ˆ x1 ; z …y; z† ˆ 0

…52† …53†

The above conditions should be satis®ed for all y 2 ‰0; bŠ and z 2 ‰0; cŠ. First, consider the variants of two and four-bricks made from the same material and having equal length. Accordingly, in the case of two-bricks it is taken x2 ˆ a=2, and in the case of four-bricks it is taken x2 ˆ a=4 x3 ˆ a=2, and x4 ˆ 3a=4. The following material properties are used in the analysis: E1 ˆ 164 GPa; E2 ˆ E3 ˆ 9:82 GPa; v12 ˆ v13 ˆ 0:24; v23 ˆ 0:34;

G12 ˆ G13 ˆ 6:78 GPa;

G23 ˆ 3:66 GPa;

Fig. 3. Illustration of a laminated rectangular bar under uniaxial extension.

…54†

290

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

where the principal material axes 1, 2 and 3 correspond to the coordinate directions x, y and z. Third degree Bernstein basis functions have been used here for the three coordinate directions. The displacements are normalized by the value u0x while the strains are normalized by e0x ˆ …u0x =a† and stresses by r0x ˆ E1 …u0x =a†. The computed distribution of the axial displacement at the loaded end of the bar is shown in Fig. 4. This con®rms that the imposed kinematic boundary condition (53) is exactly satis®ed in the solution for all y and z values. The same is true for the imposed zero value of ux at x ˆ x1 . When considering the computed normal stresses shown in Fig. 5, it is hard to see any di€erence between the results obtained with two- and four-bricks. Also, the single-brick solutions (not shown), which were obtained using third, ®fth and seventh degrees of the basis functions, gave practically the same results. Thus, the discretization error is negligible for the normal stress. Variations of the shear stresses are shown in Fig. 6. In this case the discretization error is considerable. The analysis conducted without enforcing stress continuity between the bricks (dashed lines) shows clear jumps in the stress components sxy and sxz at the border between Bricks 1 and 2. This is seen for both the two-brick and four-brick models. However, the jumps become signi®cantly smaller when using the fourbrick model. Results of the analysis performed with enforced stress continuity between the bricks, e.g., with the use of Eq. (44) (those are shown by solid lines), indicate that the jumps are no more present, yet some ``humps'' remain between Bricks 1 and 2. The humps are, certainly, an indication that there is still some discretization error in the solution. This suggestion is supported by the fact that the humps are signi®cantly smaller in the case of a four-brick model. It is also recognized that the humps are most pronounced between those bricks, which are the closest to the ®xed end, x ˆ x1 . So, the solution becomes less accurate when moving closer to the stress concentrator.

Fig. 4. Distribution of the normalized displacement ux =u0x in the y±z plane at the loaded end of the bar x ˆ a.

Fig. 5. Variations of the normalized stresses rx =r0x (O), ry =r0x (D), and rz =r0x (r) along the x-coordinate at y ˆ b, z ˆ c for the twobrick (a) and four-brick (b) models.

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

291

Fig. 6. Variations of the normalized stresses syz =r0x (O), sxz =r0x (D), and sxy =r0x (r) along the x-coordinate at y ˆ b, z ˆ c for the twobrick (a) and four-brick (b) models.

Fig. 7. Variations of the normalized displacements ux =u0x (O), uy =u0x …D† and uz =u0x (r) (a), normalized strains ex =e0x (O), ey =e0x (D) and ez =e0x (r) (b), and normalized stresses rx =r0x (O), ry =r0x (D) and rz =r0x (r) (c) along the x-coordinate at y ˆ b, z ˆ c, computed without ( ± ±) and with (- - - - -) strain continuity enforced at the interface x ˆ a=2.

The above numerical results lead to the following conclusion. When solving homogeneous structures using any discretization method, the enforcement of stress continuity between the discretization elements (e.g., continuity of the respective ®rst derivatives of the displacements) should signi®cantly improve the solution. However, unless continuity of the second and third derivatives of the displacements is not enforced, the stresses will inevitably have some humps at the boundaries between the discretization elements. Those may be not signi®cant in the zones of smooth stress variation, but become more pronounced in the zones of high stress concentration. Though applying ®ner brick meshes helps to reduce the humps, it is not possible to remove them completely, unless continuity of the second and third derivatives of the displacements is enforced. A theoretical

292

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

Fig. 8. Variations of the normalized stress rx =r0x along the interface x ˆ a=2, computed in the left (a) and right (b) layers.

explanation to this effect can be found in classical book of Novozhilov [35]. The problem is that the strain compatibility equations of 3-D elasticity may not be satis®ed if the second and third derivatives of the displacements are not continuous within the homogeneous solid. It should be pointed out that the above statement does not apply to the case of a laminated or other type inhomogeneous structures possessing discontinuities of the material properties. Further discussion of this issue can be found in [36]. Numerical example illustrating the issue is presented next. Consider a two-layer rectangular bar with the layers having equal size in the x-direction (i.e., x2 ˆ a=2) but distinct material properties. Elastic properties of the right orthotropic layer are de®ned by (54). The left orthotropic layer has all its elastic and shear moduli by a factor of 10 higher than the right layer. Poisson's ratios of the right and left layers are the same. Like in the previous example, one brick is taken in the y±z cross-section. The ®fth degree Bernstein basis functions are used in this example for the three coordinate directions. In this numerical experiment, two radically di€erent versions of the solution have been obtained. In the ®rst version, only the displacement continuity conditions (43) were enforced while in the second one the strain continuity conditions (44) at the interface x ˆ a=2 were added. The displacement variations shown in Fig. 7(a) indicate, however, that there is very little difference between the results. Their closeness may cause a suggestion that it does not matter, if the strain continuity is enforced or not. However, such a suggestion is readily denied after looking at the results for the strains and stresses shown in Figs. 7(b) and 7(c). Indeed, the effect of enforcing strain continuity is dramatic: the strain ex is continuous in this case while rx suffers a huge jump at the interface. Note that in the analysis version without enforcing continuity of ex , the stress rx is practically continuous at the interface (recall that the stress continuity has not been enforced in this analysis). Thus, enforcing interfacial strain continuity results in a qualitatively wrong solution. Another illustration shown in Fig. 8 veri®es that continuity of rx is very accurately satis®ed in the present analysis as the natural internal boundary condition for all values of y and z.

7. Transverse bending of [0°=90°=0°] simply supported plate Next example considers transverse bending of a simply supported 3-layer [0°=90°=0°] rectangular plate having length a, width b and thickness h. External surface load is sinusoidally distributed over the top surface, as shown in Fig. 9. The closed-form 3-D solution of this problem and benchmark results were presented in [37]. The following kinematic boundary conditions are imposed: …uy †xˆ0 ˆ …uy †xˆa ˆ …uz †xˆ0 ˆ …uz †xˆa ˆ 0;

…55†

…ux †yˆ0 ˆ …ux †yˆb ˆ …uz †yˆ0 ˆ …uz †yˆb ˆ 0:

…56†

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

293

Fig. 9. Illustration of the transverse bending problem.

No static boundary conditions are imposed in this analysis. It is expected that the loading condition  px   py  sin …57† …rz †zˆh ˆ q0 sin 2 a b should be accurately satis®ed as a ``natural'' boundary condition, through the respective term entering in the work of external surface forces (10). Material properties and the stress normalization factors are taken as in [37]. The example considered here corresponds to a ˆ b, a=h ˆ 4. In this analysis a single-brick is taken within each layer of the plate. Hence, total number of bricks is 3, and L ˆ 1, M ˆ 1, N ˆ 3, l ˆ 1; 2. According to the notations introduced in Fig. 1, in this case m ˆ 1; 2 and n ˆ 1; . . . ; 4. The sixth degree Bernstein basis functions were used in all the coordinate directions, I ˆ J ˆ K ˆ 6. Distributions of the displacements in the x±y, x±z and y±z planes are shown in Figs. 10 and 11. Results in Fig. 10 verify that the geometric boundary conditions (55) and (56) are exactly satis®ed, and that the solution provides smooth sinusoidal-type distribution of the displacements along the in-plane coordinates. Fig. 11 shows that the displacements are continuous at the interfaces. Besides, the distribution of ux indicates a typical higher-order ``zig±zag'' pattern through the thickness, though this pattern is not so obvious in the distribution of uy . This distinction emphasizes that the through-thickness variation of elastic

Fig. 10. Distributions of the normalized displacements ux (a), uy (b), and uz (c) in the x±y plane at z ˆ h=2.

294

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

Fig. 11. Distributions of the normalized displacements: ux in the x± z plane at y ˆ b=2 (a), uy , (b) and uz (c) in the y±z plane at x ˆ a=2.

properties signi®cantly affects the displacement variations. Also, distribution of uz indicates its minor variation in the thickness direction. Distributions of the in-plane stresses shown in Fig. 12 illustrate that the static boundary conditions at the edges of the plate, which were not enforced in the solution, are accurately satis®ed. Also, it is clearly seen that rx and ry have jumps at the interfaces. At the same time, stress distributions shown in Fig. 13 indicate that there are no visible jumps of the transverse stresses sxz , syz and rz . So, continuity of these stress components, which has not been enforced in the present solution, is satis®ed with a high accuracy. Overall, when comparing through-the-thickness distributions of displacements and stresses shown in Figs. 10±13 with the respective curves in Figs. 3±6 from [37], the closeness of results is obvious. Additional comparisons are presented in Table 1. Along with the present results and the benchmark data from [37], also the results obtained with ANSYS SOLID 46 element and another two groups of results are given. One of the groups has been published in [25], where cubic splines were used in all the coordinate directions. The other group has been published in [30], where trigonometric functions were used for the displacement approximation in the x- and y-directions, together with Bernstein polynomials in the z-direction. It is concluded from the comparison, that present results, as well as results of [25] and [30] are practically identical to the benchmark data from [37]. At the same time, results obtained with ANSYS SOLID 46 element are considerably less accurate, even with over 10,000 d.o.f. used. Besides, the values of ry seemingly converge to wrong values (possible explanation for that has been suggested in [30]). In addition, it is worth mentioning that present analysis gave the values     a b h a b h ; ; ; ; ˆ 0:99994; sxz ˆ ÿ6
10ÿ9 ; rz 2 2 2 2 2 2 and

 syz

a b h ; ; 2 2 2

 ˆ ÿ6
10ÿ9 :

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

Fig. 12. Distributions of the normalized stresses in the x±z plane: rx (a) and ry (b) at y ˆ b=2, and sxy (c) at y ˆ b=4.

Fig. 13. Distributions of the normalized stresses in the x±z plane: sxz (a) and syz (b) at y ˆ b=4, and rz (c) at y ˆ b=2.

295

296

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

Table 1 Comparison among the present results, benchmark data [37], results obtained with ANSYS SOLID 46 element [36], and results from previous versions of 3-D variational analysis [25,30], for the case of ‰0°=90°=0°Š square plate with a=h ˆ 4 Analysis approach Present [37] [25] [30] ANSYSa ANSYSb ANSYSc ANSYSd ANSYSe ANSYSf

ÿ
rx a2 ; b2 ; h2

ÿ
rx a2 ; b2 ; ÿ h2

ry

ÿa

0.8005 0:801 0:8008 0:8008 0:6849 0:6912 0:7751 0:7826 0:7867 0:7943

ÿ0.7545 ÿ0:755 ÿ0:7548 ÿ0:7548 ÿ0:6469 ÿ0:6530 ÿ0:7301 ÿ0:7371 ÿ0:7414 ÿ0:7585

0:5339 0:534 0:5341 0:5341 0:2609 0:3632 0:2742 0:2760 0:2781 0:2800

;b;h 2 2 4




ry

ÿa




; b ; ÿ h4 2 2

ÿ0:5560 ÿ0:556 ÿ0:55625 ÿ0:5563 ÿ0:2667 ÿ0:2689 ÿ0:2786 ÿ0:2804 ÿ0:2826 ÿ0:2844

ÿ
sxy 0; 0; h2

ÿ
sxz 0; b2 ; 0

syz

ÿ0:0511 ÿ0:0511 ÿ0:0511 ÿ0:0511 ÿ0:0461 ÿ0:0467 ÿ0:0490 ÿ0:0493 ÿ0:0502 ÿ0:0505

0:2559 0:256 0:2559 0:2559 0:2304 0:2324 0:2502 0:2506 0:2538 0:2544

0:2172 0:2172 0:2172 0:2171 0:1098 0:1105 0:2043 0:2116 0:2066 0:2140

ÿa 2

; 0; 0




Note: In ANSYS analysis the following numbers of d.o.f. were used: a 432, b 768, c 1404, d 2700, e 5616, and f 10,800. Present results were obtained with 3087 d.o.f. in the primary displacement approximation (5).

This shows, how accurately the external static boundary conditions at the top surface of the plate are satis®ed without enforcing them in the solution. The above example illustrates that all six stress components in thick laminated plates can be very accurately predicted by the developed analysis, even with a single-brick is used within each layer. Though, a higher degree Bernstein basis functions have to be used as to achieve such accuracy. Numerical examples considered here are rather simple and academic. Their purpose was illustrating some characteristic features, algorithmic advantages, and accuracy of the developed 3-D variational analysis approach. A more complex and practical problems have been analyzed using this approach in [43±45]. Some comparisons with experimental data are also available in the above works. 8. Conclusions The developed 3-D variational theory and numerical analysis approach provide novel analytical and computational tools for solving diverse boundary value problems of laminated anisotropic plates. The developed theory is versatile, allowing us to use any type of basis functions, selectively for each of the coordinate directions, without substantial modi®cation of the major mathematical algorithm. Also, the developed theory is material-adaptive, e.g., this provides a capability of imposing di€erent combinations of the inter-brick connectivity conditions, depending on their individual elastic properties. The approach is capable of solving practically any speci®c type of external kinematic and static boundary conditions. As shown in this work, the application of Bernstein basis functions brings an enchanting elegance to the mathematical algorithm. This is especially remarkable if the uniformly distributed displacement boundary conditions are imposed. Even more importantly, the displacement continuity conditions between the bricks and the stress continuity conditions between the same material bricks, can be expressed in a simple form of the linear relations. Besides, the developed analysis approach can be easily reformulated as an 8-noded hexahedral ``Bernstein ®nite element'' of an arbitrary order. The presented numerical examples show that if imposing transverse strain continuity between the same material bricks, the stress solution quality (its smoothness and accuracy) is considerably improved. However, if imposing transverse strain continuity between di€erent material bricks, the solution becomes totally wrong. A comparison with benchmark data for the problem of 3-D transverse bending of rectangular simply supported laminated plate validates high accuracy of the developed analysis approach and, at the same time, its superior accuracy to ANSYS SOLID 46 element. Some generalizations of the developed theory and analysis approach, aimed at solving complex geometrical con®gurations, are underway. Among those are: application of 2-D and 3-D uncoupled Bernstein

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

297

polynomials, which will provide new types of 3-D elements, and development of the curvilinear Bernstein elements, following fundamental methodologies of [6,7,11,13]. Acknowledgements This work has been performed under ®nancial support of the US Army Research Oce, project DAAH04-96-1-0057, Dr. Gary L. Anderson project monitor, and the US Air Force Materials Directorate, Wright-Patterson Air Force Base, contract F33615-95-C-5041, Dr. Steven L. Donaldson project manager. The author is especially grateful to Dr. Sergei P. Yushanov who developed computer code and provided numerical results for this paper. Appendix A First three groups of the matrix A…s† elements (11) are de®ned as following:   11…s† …s† …4† …1† …1† …s† …1† …4† …1† …s† …1† …1† …4† …s† …1† …2† …3† …1† …3† …2† aijk;pqr ˆ C11 Xip Yjq Zkr ‡ C66 Xip Yjq Zkr ‡ C55 Xip Yjq Zkr ‡ C56 Xip Yjq Zkr ‡ Xip Yjq Zkr     …s† …3† …2† …1† …2† …3† …1† …s† …3† …1† …2† …2† …1† …3† ‡ C16 Xip Yjq Zkr ‡ Xip Yjq Zkr ‡ C15 Xip Yjq Zkr ‡ Xip Yjq Zkr ; 12…s†

…s†

…2†

…3†

…1†

…s†

…3†

…2†

…1†

…s†

…1†

…4†

…1†

…s†

…4†

…1†

…1†

…s†

…1†

…1†

…4†

aijk;pqr ˆ C12 Xip Yjq Zkr ‡ C66 Xip Yjq Zkr ‡ C26 Xip Yjq Zkr ‡ C16 Xip Yjq Zkr ‡ C45 Xip Yjq Zkr …s†

…1†

…3†

…2†

…s†

…2†

…1†

…3†

…s†

…1†

…2†

…3†

…s†

…3†

…1†

…2†

‡ C25 Xip Yjq Zkr ‡ C14 Xip Yjq Zkr ‡ C46 Xip Yjq Zkr ‡ C56 Xip Yjq Zkr ; 13…s†

…s†

…2†

…1†

…3†

…s†

…3†

…1†

…2†

…s†

…1†

…2†

…3†

…s†

…1†

…3†

…2†

…A:1† …s†

…1†

…1†

…4†

aijk;pqr ˆ C13 Xip Yjq Zkr ‡ C55 Xip Yjq Zkr ‡ C36 Xip Yjq Zkr ‡ C45 Xip Yjq Zkr ‡ C35 Xip Yjq Zkr …s†

…2†

…3†

…1†

…s†

…1†

…4†

…1†

…s†

…4†

…1†

…1†

…s†

…3†

…2†

…1†

‡ C14 Xip Yjq Zkr ‡ C46 Xip Yjq Zkr ‡ C15 Xip Yjq Zkr ‡ C56 Xip Yjq Zkr ; where …1†

Xip ˆ …4†

Xip ˆ

Z

xl‡1 x

Xi …x†Xp …x† dx;

Z lxl‡1 xl

…2†

Z

Xip ˆ

xl‡1

xl

Xi …x†

dXp …x† dx; dx

…3†

Xip ˆ

Z

xl‡1

xl

dXi …x† Xp …x† dx; dx

dXi …x† dXp …x† dx: dx dx

…A:2†

The y-direction integrals, which are analogous to the above ones, can be obtained from the respective expressions (A.2) by substituting y for x, Y for X, ym for xl , ym‡1 for xl‡1 , j for i, and q for p. The z-direction integrals can be obtained from the respective expressions (A.2) by substituting z for x, Z for X, zn for xl , zn‡1 for xl‡1 , k for i, and r for p. 22…s† 23…s† Another six groups of the matrix A…s† elements (11) are de®ned as follows. Coecients aijk;pqr ; aijk;pqr and 21…s† 11…s† 12…s† 13…s† aijk;pqr are obtained from aijk;pqr ; aijk;pqr and aijk;pqr , respectively, by the following index permutations: 2 for 1, …c† …c† …c† …c† …c† …c† 3 for 2, 1 for 3, 5 for 4, 6 for 5, 4 for 6, and substitutions of Yjq for Xip , Zkr for Yjq and Xip for Zkr , where 33…s† 31…s† 32…s† 22…s† 23…s† 21…s† c ˆ 1; 2; 3; 4. After that, coecients aijk;pqr ; aijk;pqr and aijk;pqr can be obtained from aijk;pqr ; aijk;pqr and aijk;pqr , respectively, by applying the same index permutations and substitutions. ab…s†

It is easy to check that coecients aijk;pqr possess the following symmetry properties: aa…s†

aa…s†

aijk;pqr ˆ apqr;ijk ;

ab…s†

ba…s†

aijk;ijk ˆ aijk;ijk ;

ab…s†

ba…s†

aijk;pqr ˆ apqr;ijk ;

for any a; b ˆ 1; 2; 3:

…A:3†

Properties (A.3) assure that matrix A…s† is symmetric. Appendix B a…s†

Expressions of the components Qijk of vector Q…s† (11) are derived as follows. Consider sth brick having exterior surfaces x ˆ xl , x ˆ xl‡1 , y ˆ ym , y ˆ ym‡1 , z ˆ zn and z ˆ zn‡1 . Each of the above surfaces may be exposed to three arbitrarily distributed surface tractions

298

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304 …s†l

qa1 …y; z† …s†l‡1

qa1

applied at

x ˆ xl ;

…y; z† applied at x ˆ xl‡1 ;

…s†m

qa2 …x; z†

applied at

…s†m‡1 qa2 …x; z† …s†n

qa3 …x; y†

y ˆ ym ;

…B:1†

applied at y ˆ ym‡1 ; applied at z ˆ zn ;

…s†n‡1 qa3 …x; y†

applied at z ˆ zn‡1 :

Now, total work of the external surface forces is written in the form "Z Z Z Z 3 X …s† …s†l‡1 …s† qa1 …y; z†u…s† qa1 …y; z†ua…s† …xl‡1 ; y; z† dy dz …x ; y; z† dy dz ‡ W ˆ l a …s†

Z

Z

‡

…s†

Z ‡

…s†

Fl

aˆ1

Fm

Z

…s†

Fn

…s†m qa2 …x; z†u…s† a …x; ym ; z†

…s†n qa3 …x; y†u…s† a …x; y; zn †

Fl‡1

Z

Z dx dz ‡ Z dx dy ‡

…s†m‡1

qa2

…s†

Fm‡1

…x; z†u…s† a …x; ym‡1 ; z† dx dz #

Z

…s†n‡1 qa3 …x; y†u…s† a …x; y; zn‡1 †

…s†

Fn‡1

dx dy :

…B:2†

After applying Eq. (5), expression (B.2) is transformed to the form (10) with vector Q…s† de®ned by (11). Its components are expressed as following a…s†

…s†l

…s†l‡1

…s†m

…s†m‡1

…s†n

Qijk ˆ Xi …xl †Qa1jk ‡ Xi …xl‡1 †Qa1jk ‡ Yj …ym †Qa2ik ‡ Yj …ym‡1 †Qa2ik

…s†n‡1

‡ Zk …zn †Qa3ij ‡ Zk …zn‡1 †Qa3ij ; …B:3†

where …s†l Qa1jk

Z ˆ

Z

…s† Fl

Z

…s†m

Qa2ik ˆ Qa3ij ˆ

…s†l‡1 Qa1jk

dy dz;

…s†m

qa2 …x; z†Xi …x†Zk …z† dx dz;

…s†

Fm

Z

Z

…s†n

Z

…s†l qa1 …y; z†Yj …y†Zk …z†

…s†n

qa3 …x; y†Xi …x†Yj …y† dx dy;

…s† Fn

Z ˆ

…s†n‡1

…s† Fl‡1

Z

…s†m‡1

Qa2ik Qa3ij

Z

ˆ

Z

Fm‡1

…y; z†Yj …y†Zk …z† dy dz;

…s†m‡1

qa2

…s†

Z ˆ

…s†l‡1

qa1

Z

…s† Fn‡1

…s†n‡1

qa3

…x; z†Xi …x†Zk …z† dx dz;

…B:4†

…x; y†Xi …x†Yj …y† dx dy:

By the use of integrals (B.4), one can evaluate work of external surface forces acting on the exterior of the brick for any given distributions of the surface tractions. Appendix C Four integrals de®ned by Eq. (A.2) are written in the following form using parametric representation of the Bernstein basis functions: …1†

I…1†

Xip ˆ …xl‡1 ÿ xl †Bip ; where I…1† Bip I…4† Bip

Z ˆ Z ˆ

1 0 1 0

BIi …t†BIp …t† dBIi …t† dt

dt;

dBIp …t† dt

…2†

I…2†

…3†

Xip ˆ Bip ;

I…2† Bip

dt:

Z ˆ

0

I…3†

Xip ˆ Bip ;

1

BIi …t†

dBIp …t† dt; dt

…4†

Xip ˆ

I…3† Bip

1 I…4† B ; xl‡1 ÿ xl ip Z ˆ

0

1

dBIi …t† I Bp …t† dt; dt

…C:1†

…C:2†

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

299

To perform further derivations, we ®rst make use of the di€erentiation formula ÿ N ÿ1
dBNn ˆ N Bnÿ1 ÿ BnN ÿ1 ; dt

…C:3†

which is valid for n ˆ 1; . . . ; N ÿ 1. The cases n ˆ 0 and n ˆ N require a special consideration. The following expressions are derived: Z 1   I…2† Iÿ1 BIi Bpÿ1 ÿ BIÿ1 dt; Bip ˆ I p 0

Z

I…3†

Bip ˆ I I…4†

Bip ˆ I 2

1

ÿ

0

Z

1 0


Iÿ1 ÿ BiIÿ1 BIp dt; Biÿ1

ÿ

Iÿ1 Biÿ1 ÿ BiIÿ1

…C:4†


 Iÿ1 Bpÿ1 ÿ BIÿ1 dt; p

which are valid for i; p ˆ 1; 2; . . . ; I ÿ 1. The following speci®c case of the binomial formula is useful when evaluating integrals (C.4): N

…1 ÿ t† ˆ

N X N !…ÿ1†n tn : n!…N ÿ n†! nˆ0

…C:5†

Analytical expressions of the integrals (C.2) are then obtained as following I…1†

Bip ˆ I…2†

Bip ˆ

2Iÿiÿp m I!I!…2I ÿ i ÿ p†! X …ÿ1† ; i!…I ÿ i†!p!…I ÿ p†! mˆ0 m!…2I ÿ i ÿ p ÿ m†!…i ‡ p ‡ m ‡ 1†

I!I!…2I ÿ i ÿ p ÿ 1†! i!…I ÿ i†!p!…I ÿ p ÿ 1†! i ˆ 0; 1; . . . ; I;

I…3†

mˆ0

m

…ÿ1† ; m!…2I ÿ i ÿ p ÿ m ÿ 1†!…i ‡ p ‡ m ‡ 1†

p ˆ 1; . . . ; I ÿ 1;

I!I!…2I ÿ i ÿ p ÿ 1†! i!…I ÿ i ÿ 1†!p!…I ÿ p†! i ˆ 1; . . . ; I ÿ 1;

I…4†

2Iÿiÿpÿ1 X

…C:7†

2Iÿiÿp m X I!I!…2I ÿ i ÿ p†! …ÿ1† …i ÿ 1†!…I ÿ i†!p!…I ÿ p†! mˆ0 m!…2I ÿ i ÿ p ÿ m†!…i ‡ p ‡ m†

ÿ

Bip ˆ

…C:6†

2Iÿiÿp m X I!I!…2I ÿ i ÿ p†! …ÿ1† i!…I ÿ i†!…p ÿ 1†!…I ÿ p†! mˆ0 m!…2I ÿ i ÿ p ÿ m†!…i ‡ p ‡ m†

ÿ

Bip ˆ

i; p ˆ 0; 1; . . . ; I;

2Iÿiÿpÿ1 X mˆ0

…ÿ1†m ; m!…2I ÿ i ÿ p ÿ m ÿ 1†!…i ‡ p ‡ m ‡ 1†

p ˆ 0; 1; . . . ; I;

…C:8†

2Iÿiÿp X I!I!…2I ÿ i ÿ p†! …ÿ1†m …i ÿ 1†!…I ÿ i†!…p ÿ 1†!…I ÿ p†! mˆ0 m!…2I ÿ i ÿ p ÿ m†!…i ‡ p ‡ m ÿ 1†

ÿ ÿ ‡

I!I!…2I ÿ i ÿ p ÿ 1†! …i ÿ 1†!…I ÿ i†!p!…I ÿ p ÿ 1†!

2Iÿiÿpÿ1 X

I!I!…2I ÿ i ÿ p ÿ 1†! i!…I ÿ i ÿ 1†!…p ÿ 1†!…I ÿ p†!

2Iÿiÿpÿ1 X

I!I!…2I ÿ i ÿ p ÿ 2†! i!…I ÿ i ÿ 1†!p!…I ÿ p ÿ 1†! i; p ˆ 1; . . . ; I ÿ 1;

mˆ0

mˆ0

2Iÿiÿpÿ2 X mˆ0

m

…ÿ1† m!…2I ÿ i ÿ p ÿ m ÿ 1†!…i ‡ p ‡ m† m

…ÿ1† m!…2I ÿ i ÿ p ÿ m ÿ 1†!…i ‡ p ‡ m† m

…ÿ1† ; m!…2I ÿ i ÿ p ÿ m ÿ 2†!…i ‡ p ‡ m ‡ 1† …C:9†

300

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

I…2†

Bi0 ˆ I…2†

BiI

ˆ

m X I
I!…2I ÿ i ÿ 1†! 2Iÿiÿ1 …ÿ1† ; m!…2I ÿ i ÿ m ÿ 1†!…i ‡ m ‡ 1† i!…I ÿ i†! mˆ0 m Iÿi I
I! X …ÿ1† ; i! mˆ0 m!…I ÿ i ÿ m†!…I ‡ i ‡ m†

I…3†

B0p ˆ ÿ I…3†

BIp ˆ I…4†

‡ I…4†

Bi0 ˆ ÿ ‡ I…4†

BiI

ˆ

I…4†

I…4†

I…4†

BII ˆ

p ˆ 0; 1; . . . ; I;

p ˆ 0; 1; . . . ; I;

…C:12†

…C:13†

2Iÿpÿ1 m I
I!…2I ÿ p ÿ 1†! X …ÿ1† …p ÿ 1†!…I ÿ p†! mˆ0 m!…2I ÿ p ÿ m ÿ 1†!…p ‡ m† 2Iÿpÿ2 I
I!…2I ÿ p ÿ 2†! X …ÿ1†m ; p!…I ÿ p ÿ 1†! m!…2I ÿ p ÿ m ÿ 2†!…p ‡ m ‡ 1† mˆ0

p ˆ 1; . . . ; I ÿ 1;

…C:14†

m X I
I!…2I ÿ i ÿ 1†! 2Iÿiÿ1 …ÿ1† …i ÿ 1†!…I ÿ i†! mˆ0 m!…2I ÿ i ÿ m ÿ 1†!…i ‡ m† m X I
I!…2I ÿ i ÿ 2†! 2Iÿiÿ2 …ÿ1† ; i!…I ÿ i ÿ 1†! m!…2I ÿ i ÿ m ÿ 2†!…i ‡ m ‡ 1† mˆ0

m X I
I! Iÿiÿ1 …ÿ1† ; i! mˆ0 m!…I ÿ i ÿ m ÿ 1†!…I ‡ i ‡ m†

Iÿpÿ1 m I
I! X …ÿ1† ; p! mˆ0 m!…I ÿ p ÿ m ÿ 1†!…I ‡ p ‡ m†

I…4†

B0I ˆ BI0 ˆ ÿI
I! B00 ˆ

…C:11†

p ˆ 1; . . . ; I ÿ 1;

…C:15†

i ˆ 1; . . . ; I ÿ 1;

…C:16†

Iÿp m I
I! X …ÿ1† …p ÿ 1†! mˆ0 m!…I ÿ p ÿ m†!…I ‡ p ‡ m ÿ 1†

ÿ I…4†

…C:10†

Iÿi I
I! X …ÿ1†m …i ÿ 1†! mˆ0 m!…I ÿ i ÿ m†!…I ‡ i ‡ m ÿ 1†

ÿ

BIp ˆ

i ˆ 0; 1; . . . ; I;

2Iÿ1ÿp m I
I!…2I ÿ p ÿ 1†! X …ÿ1† ; m!…2I ÿ p ÿ m ÿ 1†!…p ‡ m ‡ 1† p!…I ÿ p†! mˆ0

Iÿp m I
I! X …ÿ1† ; p! mˆ0 m!…I ÿ p ÿ m†!…I ‡ p ‡ m†

B0p ˆ ÿ

i ˆ 0; 1; . . . ; I;

Iÿ1 X mˆ0

p ˆ 1; . . . ; I ÿ 1;

…C:17†

m

…ÿ1† ; m!…I ÿ m ÿ 1†!…I ‡ m†

…C:18†

I2 ; 2I ÿ 1

…C:19†

I2 : 2I ÿ 1

…C:20†

The respective y- and z-direction integrals can be evaluated analogously. Appendix D Here we assume for simplicity that all bricks in the mosaic body are orthotropic. The derivation for a more general case of anisotropy can be performed in a similar way. Consider ®rst equation of (30), which

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

301

corresponds to a ˆ 1, By the use of stress±strain and strain-displacement relations (1) and (2), the equation is written in the following form: ! ! ! s…l;m;n† s…l;m;n† s…l;m;n† ou1 ou2 ou3 s…l;m;n† s…l;m;n† s…l;m;n† ‡ C12 ‡ C13 C11 ox oy oz xˆxl xˆxl xˆxl ! ! ! s…l‡1;m;n† s…l‡1;m;n† s…l‡1;m;n† ou1 ou2 ou3 s…l‡1;m;n† s…l‡1;m;n† s…l‡1;m;n† ˆ C11 ‡ C12 ‡ C13 …D:1† ox oy oz xˆxl

xˆxl

xˆxl

Further, by substituting (5) with (35)±(37) in (D.1), it is obtained J X K      X I I s…l;m;n† 1;s…l;m;n† 1;s…l;m;n† s…l‡1;m;n† 1;s…l‡1;m;n† 1;s…l‡1;m;n† C11 UIjk U1jk ÿ UIÿ1jk ÿ C11 ÿ U0jk xl ÿ xlÿ1 xl‡1 ÿ xl jˆ0 kˆ0 h i J s…l;m;n† 2;s…l;m;n† s…l‡1;m;n† 2;s…l‡1;m;n† dPj …y† K UIjk U0jk Pk …z† ÿ C12  PjJ …y†PkK …z† ‡ C12 dy  h i dP K …z† s…l;m;n† 3;s…l;m;n† s…l‡1;m;n† 3;s…l‡1;m;n† UIjk U0jk ˆ 0: ÿ C13 PjJ …y† k ‡ C13 dz

…D:2†

Assuming that elastic properties of the layers do not vary along the x-coordinate, e.g., that s…l;m;n† s…l‡1;m;n† s…l;m;n† s…l‡1;m;n† s…l;m;n† s…l‡1;m;n† ˆ C11 , C12 ˆ C12 and C13 ˆ C13 , Eq. (D.2) yields C11  J X K     X I I s…l;m;n† 1;s…l;m;n† 1;s…l;m;n† 1;s…l‡1;m;n† 1;s…l;m;n† C11 UIjk U1jk ÿ UIÿ1jk ÿ ÿ UIjk PjJ …y†PkK …z† ˆ 0: xl ÿ xlÿ1 xl‡1 ÿ xl jˆ0 kˆ0 …D:3† Further, taking into account that all basis functions are linearly independent, then accounting that I 6ˆ 0 s…l;m;n† 6ˆ 0, and incorporating ®rst equation of (43), into (D.3), the following equation is obtained: and C11 xl‡1 ÿ xlÿ1 1;s…l;m;n† xl‡1 ÿ xl 1;s…l;m;n† 1;s…l‡1;m;n† ˆ U ÿ U : …D:4† U1jk xl ÿ xlÿ1 Ijk xl ÿ xlÿ1 Iÿ1jk The above equation has to be satis®ed for j ˆ 0; 1; . . . ; J and k ˆ 0; 1; . . . ; K. If (D.4) is satis®ed, then continuity of the strain components e11 ; e12 ; e13 and stress components r11 ; r12 ; r13 is assured at all those borders between the same material bricks, which are perpendicular to the x-axis. After performing analogous derivations for the second (a ˆ 2) and third (a ˆ 3) equations from (30) and s…l;m;n† s…l‡1;m;n† s…l;m;n† s…l‡1;m;n† ˆ C66 and C55 ˆ C55 , it is obtained assuming that C66 x ÿ x x ÿ x l‡1 lÿ1 l‡1 l 2;s…l‡1;m;n† 2;s…l;m;n† 2;s…l;m;n† ˆ U ÿ U ; …D:5† U1jk xl ÿ xlÿ1 Ijk xl ÿ xlÿ1 Iÿ1jk 3;s…l‡1;m;n†

U1jk

ˆ

xl‡1 ÿ xlÿ1 3;s…l;m;n† xl‡1 ÿ xl 3;s…l;m;n† U ÿ U : xl ÿ xlÿ1 Ijk xl ÿ xlÿ1 Iÿ1jk

Now it is recognized that all three sets of Eqs. (D.4)±(D.6), can be written in a uni®ed form xl‡1 ÿ xlÿ1 a;s…l;m;n† xl‡1 ÿ xl a;s…l;m;n† a;s…l‡1;m;n† ˆ U ÿ U ; a ˆ 1; 2; 3: U1jk xl ÿ xlÿ1 Ijk xl ÿ xlÿ1 Iÿ1jk

…D:6†

…D:7†

When applying analogous derivations to the inter-brick boundaries perpendicular to y-axis, assuming that elastic properties of the layers do not vary in the y-direction, and accounting for the second equation of (43), the result obtained from (31) is ym‡1 ÿ ymÿ1 a;s…l;m;n† ym‡1 ÿ ym a;s…l;m;n† a;s…l;m‡1;n† ˆ U ÿ U ; a ˆ 1; 2; 3: …D:8† Ui1k ym ÿ ymÿ1 iJk ym ÿ ymÿ1 iJ ÿ1k This equation has to be satis®ed for i ˆ 0; 1; . . . ; I and k ˆ 0; 1; . . . ; K. If (D.8) is satis®ed, then continuity of the strain components e22 ; e12 ; e23 and stress components r22 ; r12 ; r23 is assured at all those borders between the same material bricks, which are perpendicular to the y-axis.

302

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

Further, the following result is obtained from (32) for the stress component r33 ( " # s…l;m;n†  I X J  C s…l;m;n‡1†   X C33 3;s…l;m;n† 3;s…l;m;n† 3;s…l;m;n‡1† 3;s…l;m;n‡1† 33 K UijK Uij1 ÿ UijKÿ1 ÿ ÿ Uij0 PiI …x†PjJ …y† z ÿ z z ÿ z n nÿ1 n‡1 n iˆ0 jˆ0 h i I s…l;m;n† 1;s…l;m;n† s…l;m;n‡1† 1;s…l;m;n‡1† dPi …x† J Pj …y† UijK Uij0 ÿ C13 ‡ C13 dx ) J h i dP …y† j s…l;m;n† 2;s…l;m;n† s…l;m;n‡1† 2;s…l;m;n‡1† ˆ0 UijK Uij0 ‡ C23 ÿ C23 PiI …x† dy

…D:9†

for the stress component r13 ( " s…l;m;n†  I X J h i I  X C55 s…l;m;n† 3;s…l;m;n† s…l;m;n‡1† 3;s…l;m;n‡1† Pi …x† J 1;s…l;m;n† 1;s…l;m;n† Pj …y† ‡ K C55 UijK Uij0 UijK ÿ C55 ÿ UijKÿ1 zn ÿ znÿ1 dx iˆ0 jˆ0 ) # s…l;m;n‡1†   C55 1;s…l;m;n‡1† 1;s…l;m;n‡1† I J Uij1 ÿ Uij0 Pi …x†Pj …y† ˆ 0 ÿ zn‡1 ÿ zn

…D:10†

and for the stress component r23 ( " s…l;m;n†  I X J h i  X dPjJ …y† C44 s…l;m;n† 3;s…l;m;n† s…l;m;n‡1† 3;s…l;m;n‡1† 2;s…l;m;n† 2;s…l;m;n† I ‡K C44 UijK Uij0 UijK ÿ C44 Pi …x† ÿ UijKÿ1 dy zn ÿ znÿ1 iˆ0 jˆ0 ) # s…l;m;n‡1†   C44 2;s…l;m;n‡1† 2;s…l;m;n‡1† I J Uij1 …D:11† ÿ Uij0 Pi …x†Pj …y† ˆ 0: ÿ zn‡1 ÿ zn s…l;m;n†

s…l;m;n‡1†

s…l;m;n†

s…l;m;n‡1†

s…l;m;n†

s…l;m;n‡1†

s…l;m;n†

s…l;m;n‡1†

6ˆ C13 , C23 6ˆ C13 , C33 6ˆ C33 , C44 6ˆ C44 , and It is anticipated that if C13 s…l;m;n† s…l;m;n‡1† 6ˆ C55 , then it will be hardly possible to analytically satisfy Eqs. (D.9)±(D.11). Thus, we do not C55 have a better choice now as to ignore these continuity equations and expect that there will be no signi®cant jumps in the computed transverse stresses at the interfaces. s…l;m; n† s…l;m; n‡1† ˆ C13 , Now assume that for some subset f ng of the full set fng the following takes place: C13 s…l;m; n† s…l;m; n‡1† s…l;m; n† s…l;m; n‡1† s…l;m; n† s…l;m; n‡1† s…l;m; n† s…l;m; n‡1† ˆ C13 , C33 ˆ C33 , C44 ˆ C44 , and C55 ˆ C55 . In the case of lamiC23 nated plates this applies to those inter-brick boundaries, perpendicular to the z-axis, which are located inside the layers. Consequently, the following equation is obtained from (D.9)±(D.11) with account for the third equation of (43) zn‡1 ÿ znÿ1 a;s…l;m;n† zn‡1 ÿ zn a;s…l;m;n† a;s…l;m; n‡1† ˆ U ÿ U : …D:12† Uij1 zn ÿ znÿ1 ijK zn ÿ znÿ1 ijKÿ1 The above equation should be satis®ed for i ˆ 0; 1; . . . ; I and j ˆ 0; 1; . . . ; J . If (D.12) is satis®ed, then continuity of the strain components e33 ; e13 ; e23 and stress components r33 ; r13 ; r23 is assured at all those borders between the same material bricks, which are perpendicular to the z-axis. Thus, if applied altogether, the derived Eqs. (D.7), (D.8), (D.12) assure continuity of all the strain components e11 ; e12 ; e13 ; e22 ; e23 ; e33 and stress components r11 ; r12 ; r13 ; r22 ; r23 , and r33 between all of the bricks belonging to the same homogeneous layer. Again, it should be emphasized that Eq. (D.12) do not apply to the interfaces between distinct layers. References [1] [2] [3] [4]

R.H. Gallagher, J. Padlog, P.P. Bijlaard, Stress analysis of heated complex shapes, J. Aero-Space Science 32 (1962) 700±707. R.J. Melosh, Structural analysis of solids, J. Structural Division, Proceedings of ASCE 89, No. ST4, 1963, pp. 205±223. J.H. Argyris, Matrix analysis of three-dimensional elastic media ± small and large displacements, AIAA J. 3 (1965) 45±51. J.H. Argyris, Three-dimensional anisotropic and inhomogeneous media ± matrix analysis for small and large displacements, Ingenieur Archiv 34 (1965) 33±35.

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

303

[5] J.H. Argyris, I. Fried, D.W. Scharpf, The TET 20 and TEA 8 elements for the matrix displacement method, The Aeronaut. J. Royal Aeronaut. Soc. 72 (691) (1968) 618±623. [6] J.H. Argyris, The LUMINA element for the matrix displacement method, The Aeronaut. J. Royal Aeronaut. Soc. 72 (690) (1968) 514±517. [7] J.H. Argyris, K.E. Buck, D.W. Scharpf, H.M. Hilber, G. Mareczek, Some new elements for the matrix displacement methods, in: Proceedings of 2nd Conference on Matrix Methods in Struct. Mech., Air Force Inst. of Techn., Wright Patterson AF Base, Ohio, 1968, pp. 333±397. [8] J.H. Argyris, I. Fried, D.W. Scharpf, The Hermes 8 element for the matrix displacement method, The Aeronaut. J. Royal Aeronaut. Soc. 72 (691) (1968) 613±617. [9] J.G. Ergatoudis, Isoparametric ®nite elements in two- and three-dimensional stress analysis, Ph.D. Thesis, University of Wales, Swansea, 1968. [10] J.G. Ergatoudis, B.M. Irons, O.C. Zienkiewicz, Three-dimensional analysis of arch dams and their foundations, in: Proceedings of Symposium on Arch Dams, Institution of Civil Engineers, London, 1968, pp. 20±21. [11] O.C. Zienkiewicz, B.M. Irons, J. Ergatoudis, S. Ahmad, F.C. Scott, Iso-parametric and associated element families for two- and three-dimensional analysis, in: I. Holand, K. Bell (Eds.), Finite Element Methods in Stress Analysis, chapter 13 TAPIR, The Techn. University of Norway, Trondheim, Norway, 1969, pp. 383±432. [12] R.W. Clough, Comparison of three-dimensional ®nite elements, in: Proceedings of Symposium on Application of Finite Element Methods in Civil Engineering, Vanderbilt University, Nashville, Tennessee, ASCE, New York, 1969, pp. 1±26. [13] O.C. Zienkiewicz, B.M. Irons, F.C. Scott, J.S. Campbell, Three-dimensional stress analysis, in: B. Fraeijs de Veubeke (Ed.), Proceedings of Symposium on High Speed Computing of Elastic Structures, August 23±28, 1970, vol. 1, University of Liege, Belgium, 1970, pp. 413±432. [14] Y. Rashid, Three-dimensional analysis of elastic solids, Int. J. Solids Structures, Pt. I 5 (1969) 1311±1332. [15] Y. Rashid, Three-dimensional analysis of elastic solids, Int. J. Solids Structures, Pt. II 6 (1970) 195±207. [16] G.L. Rigby, G.M. McNeice, A strain energy basis for studies of element sti€ness matrices, AIAA J. 10 (1972) 1490±1493. [17] O.C. Zienkiewicz, The Finite Element Method, 3rd ed., McGraw-Hill, London, 1977. [18] R.H. Gallagher, Finite Element Analysis: Fundamentals, Prentice-Hall, Englewood Cli€s, NJ, 1975. [19] T.Y. Yang, Finite Element Structural Analysis, Prentice-Hall, Englewood Cli€s, NJ, 1986. [20] D.J. Dawe, Matrix and Finite Element Displacement Analysis, Clarendon Press, Oxford, 1984. [21] A.K. Noor, W.S. Burton, J.M. Peters, Assessment of computational models for multilayered composite cylinders, in: A.K. Noor, T. Belytschko, J. Simo (Eds.), Analytical and Computational Models of Shells, CED-vol. 3, ASME, New York, 1989, pp. 419± 441. [22] A.K. Noor, W.S. Burton, Assessment of computational models for multilayered composite shells, Appl. Mech. Rev. 43 (4) (1990) 67±97. [23] A.K. Noor, W.S. Burton, Assessment of computational models for multilayered anisotropic plates, Composite Structures 14 (1990) 233±265. [24] A.E. Bogdanovich, Spline function aided analysis of inhomogeneous materials and structures, in: J.N. Reddy, K.L. Reifsnider (Eds.), Proceedings of IUTAM Symposium Local Mech. Concepts for Composite Material System, Blacksburg, Virginia, 1991, Springer, Berlin, Heidelberg, New York, 1992, pp. 355±382. [25] A.E. Bogdanovich, Three-dimensional analysis of laminated composite plates, in: T. Chandra, A.K. Dhingra (Eds.), Procceedings of International Conference on Advanced Composite Materials, University of Wollongong, Australia, 1993, TMS Minerals/ Metals/Materials, Warrendale, Pennsylvania, 1993, pp. 79±91. [26] J.N. Reddy, D.H. Robbins, Theories and computational models for composite laminates, Appl. Mech. Rev. 47(6), Pt 1 (1994) 147±169. [27] F.K. Bogner, R.L. Fox, L.A. Schmit, Jr., The generation of inter-element-compatible sti€ness and mass matrices by the use of interpolation formulas, in: Matrix Meth. in Struct. Mech., The Conference held at Wright Patterson AF Base, 1965, Dayton, Ohio, 1965, pp. 397±443. [28] A.N. Palazotto, W.P. Witt, Formulation of a non-linear compatible ®nite element for the analysis of laminated composites, Computers Structures 21 (1985) 1213±1234. [29] R.L. Hinrichsen, A.N. Palazotto, The non-linear ®nite element analysis of thick composite plates using cubic spline functions, AIAA J. 24 (1986) 1836±1842. [30] A.E. Bogdanovich, C.M. Pastore, B.P. Deepak, A comparison of various 3-D approaches for the analysis of laminated composite structures, Composites Engrg. 5 (1995) 1105±1126. [31] A.E. Bogdanovich, Three-dimensional analysis of anisotropic spatially reinforced structures, in: John D. Buckley (Ed.), Fiber-Tex 1992, The Sixth Conference on Adv. Engrg. Fibers and Textile Struct. for Composites, 1992, Philadelphia, Pennsylvania, NASA Conf. Publ. 3211, 1992, pp. 271±304; also published in Composites Manufacturing 4 (1993) 173±186. [32] A.E. Bogdanovich, A.B. Birger, Application of spline functions to the analysis of local boundary e€ects in laminated composite plates, in: D.H. Allen, D.C. Lagoudas (Eds.), Proceedings of ASME Winter Annual Meeting, Symp. on Damage Mech. in Compos., Anaheim, California, 1992, AMD-vol. 150/AD-vol. 32, ASME, New York, 1992, pp. 275±293. [33] A.E. Bogdanovich, A.B. Birger, Stress ®elds and initial failure in a cantilever laminated composite plate under transverse and inplane loads, in: L. Librescu (Ed.), Proceedings of the 1st Joint Mech. Meeting of ASME/ASCE/SES, Symposium on Non-Classical Probl. of the Theory and Behav. of Struct. Exposed to Complex Environ. Cond., Charlottesville, Virginia, 1993, AMD-vol. 164, ASME, New York, 1993, pp. 155±167.

304

A.E. Bogdanovich / Comput. Methods Appl. Mech. Engrg. 185 (2000) 279±304

[34] A.E. Bogdanovich, A.B. Birger, Three-dimensional stress ®eld analysis in uniformly loaded, simply supported composite plates, Computers Structures 52 (1994) 237±257. [35] V.V. Novozhilov, Theory of Elasticity, Sudpromgiz, Leningrad, 1958, in Russian; English translation by J.K. Lusher, Pergamon Press, New York, 1961. [36] A.E. Bogdanovich, C.M. Pastore, Mechanics of Textile and Laminated Composites, Chapman & Hall, London, 1996. [37] N.J. Pagano, Exact solutions for rectangular bidirectional composites and sandwich plates, J. Compos. Mater. 4 (1970) 20±34. [38] I.D. Faux, M.J. Pratt, Computational Geometry for Design and Manufacture, Ellis Harwood, England, 1981. [39] S. Bu-Qing, L. Ding-Yuan, Computational Geometry Curve and Surface Modeling, Academic Press, San Diego, NY, 1989. [40] G. Farin, Curves and Surfaces for Computer Aided Geometric Design, Academic Press, New York, 1988. [41] W. Feller, An Introduction to Probability Theory and Its Applications, 2nd ed., vol. II, Wiley, New York, 1974. [42] P.J. Davis, Interpolation and Approximation, Dover Publications, New York, 1975. [43] A.E. Bogdanovich, 3-D strength prediction of composite materials, Interim Techn. Rep. WLTR-98, Materials Directorate, Wright Patterson Air Force Base, Ohio, 1997, Part B, pp. 38±89. [44] A.E. Bogdanovich, S.P. Yushanov, 3-D progressive failure analysis of bonded composite joints, in: A Coll. of Techn. Pap. of the 39th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dyn., and Mater. Conf. and Exhibit, Long Beach, California, 1998, Pt. 3, AIAA, Reston, Virginia, 1998, pp. 1984±1993. [45] A.E. Bogdanovich, S.P. Yushanov, Progressive failure analysis of adhesive bonded joints with laminated composite adherends, in: Proceedings of the American Society for Composites, 13th Annual Techn. Conf., Baltimore, MD, September 21±23, 1998, Technomic Publishing Co. pp. 452±464.