Ghostscript wrapper for D:\Digitalizacja\MTS88_t26_z1_4_PDF_artykuly\04mts88_t26_zeszyt_4.pdf 576 A .  POCESKJ,  G .  KOKA.LAN0V equivalent  nodal  forces  and  nodal  deformations,  which  represent  the  element  matrix. This  concept  of development  of the  elements was applied  in the problems  of plane  stress and  plate  bending  elements  [2, 3]. The  three­dimensional  element  which we consider  in this  paper, is based  on the plane stress  element  (Fig. la).  This  plane  stress  element  gives  very  good  results.  For  instance * • N Fig.  1.  Mixed  plane  stress  element  (a), cantilever  represented  by  one element  (b) in  the analysis  of a  cantilever  (Fig. lb) even  one element,  with  totaly  5 equations, gives exactly the beam, solution! Such results were obtained  primarily due to particular  distribu­ tion  of the  unknowns  (degrees of freedom).  Note that  the stresses and the displacements, were taken  at the midside nodes. In such  cases the degrees  of  freedom,  the  displacements and the stresses, are not independent of each other. For instance by variation of the stresses we can develope  the equations  of the compatibility  of  displacements,  which  are already satisfied  due to the presence  of the  displacements  as the degrees  of freedom  at the same nodes. The  element  was developed  by  application  of the  direct  method  of  development  of finite  elements  [6]. In  this  way the parasitic  shear  stresses  Sxy  =f(NXi,Nyi),  functions of the axial stresses, were automatically excluded. The exclusion of the stresses was another reason  for the very  good  behaviour  of the  element.  The element  developed in this way is the same as the element developed by the assumption of independent stresses and  deforma­ tions. It means  that, in the case of energetic approach, the stresses and the displacements have to be independently  assumed!  The correctness  of such  assumption  for the first  time was  proved  in  Ref. [4]. The  element  which  we discuss  is a  simple  prismatic  three­dimensional  element. The element  is our first  step  in the development  of mixed  three­dimensional  elements.  The development  of the element has to  show the way of  development  of the  mixed  elements and  the accuracy  which  could  be expected.  The element  presented  here  gives  very  good accuracy  and is  very  promising  for  further  development  of  the mixed FEM. 2.  Prismatic  element  with  36 d.o.f. The  element  is presented  on Fig. 2. At the corner  nodes the unknowns  are the stress components  (8 x 3), and at the edge nodes are the displacements in the direction of the edges (12 x 1),  it  makes  36 d.o.f.  in  total. M I X E D  FINITE  ELEMENTS 577 Fig.  2.  Prismatic  mixed  element  with  36 d.o.f. 2.1.  The deformation  shape  function.  The  number  of  parameters  defining  the  D S h F  is equal to the number  of the degrees of freedom — 36. It  means that the deformation  shape function  should be defined  by 3 polynomials with 12 parameters each, for the 3 deformation components  u,  v  and  w  separately.  However,  an  assumption  of  12  terms  polynomial, as  the  standard  procedure  in  the  FEM,  can  be  misleading.  We  shall  discuss  this  below. Here  somewhat  different  procedure  will  be  considered. If an edge of the element, for  instance the edge  1­4,  Fig. 2, is considered  as  an  axially loaded  rod,  the  deformations  in  the  rod  will  be  defined  by  the  following  expression: 2 These  deformations  are  translated  in  the  xy  plane,  in  such  a  way  that  the  deformations at the edge 2 ­ 3  vanishe, and next translated along z axis so that the side 5­^8  the  deforma­ tions vanish. The contribution  of the nodal parameters  U9,  slx  and £4x to the  deformation U  in  the  element  is  defined  as  follows: (2) In  this  way  the  complete  DShF  can  be  defined  and  represented  as  follows: V W x  ®u  0  0  0  0 10  0  $v  &„ 0  0 0 0 0 0 0, sxi (3) where  we  have: =  Y  (i+£o/2) 1(1+no) a  + Co), W,  t ­  1,8 (4) ­o­ (1 +  f 0) (1 + Vo) (1 + C 578  A.  POCESKI,  G.  KOKALANOV { ( l  + W,  t=  9,12 (1 + Co),  i =  13,  16 ,  * = 1 7 , 20 It  is interesting  to note  that  the polynomial  by which for instance  deformations  U  are defined  is  of the following  type: U  —  a1 + a2X + a3y  + a4.z + asx z  + a6xy+a­!xz+a8y2+a9x 2y  + alL0x 2z+ + anxyz+a12x 2yz.  (6) The  polynomial  is not a  "complete"  one,  as usually  is required.  It is evident  that  such a polynomial is very difficult  to assume in advance. The polynomial corresponds to a Iline­ ar  variation  of  strains  (stresses),  similarly  to the stresses  at the element  in Fig.  2. 2.2. The element  matrix. The element  matrix  can be represented  as  follows: • "n 'uvw  I FLw K J - The  vector  of nodal  parameters  is the following: rf'=  [Nlx  ...N8x...Nly...NSy  ...Nlz  ...N8Z...  U9 ... Ui2V1% ... V16WX ... W20].  (8) The  first  row submatrices  in the element  matrix  represents  the flexibility  submatrices. The submatrix  Fn gives nodal displacements in direction of the axial stresses due to the same • stresses. Submatrix  FH1)W affords  the  displacements due to the midedge node displacements. The  second row represents the stiffness  submatrices.  Submatrix  F£„w gives the nodal  forces in the  direction  of the  edge nodes  due to the axial stresses, and  submatrix  K­nodal  forces in  the  same  direction  due to  the  edge  node  displacements. 2.2.1. Flexibility  submatrices  F„  and  Fuow.  From  the  DShF  (3 ­ 5) for |  = ­ 1  the  deforma­ tions  on the surface  1,2, 6,5 are  derived.  These  deformations, for instance  due to slx are as  follows: U =  exl®xi  ­  ­ ^ | 1  exl(l  ­tj)  (1 ­ 0 .  (9) The  volume  of  these  deformations is: / /  Udydz = —  abcsxl. The volume is equal to the nodal deformations  1,2,6, 5 due to the same strain. The distribu­ tion  due to  the particular  node  is as follows: FnlJ^fJu8U(et)dydz, where  U is as Exp. (9) and  dU(si) is the variation  of the  displacement  on the particular strains.  In this  way the derived  nodal  displacements  due to exl  take  the following  form: [U.UzUeUsV  = ^^­[4  2  1 2). M I XE D   FIN ITE ELEMENTS 579 F rom  £«!  we  obtain  deformations  also  on  the surface  3, 4, 8,  7.  The nodal  deformations derived  in  the  same  way  for  these  nodes  are  as  follows: 2  4  2  1]. These two  submatrices  define  the first  column of  submatrix  F B. In  this  way  the complete submatrix  Fn  can be derived. The submatrix can be derived  by  an energetic way  as  follows: F«i,= The  values  of the same nodal displacements  due to  e xl   derived  in this way  are as  follows: [U ±   U 2   U 6   U s   U 3   U 4   U 8   U n ]'  =   ^ ^ -   [8  4  2  4  2  4  2  1]. The  coefficients  derived in both ways are not the same, bu t their sum in one row  or column is  the  same,  equal  to  1.  Therefore  submatrix  F„  defined  in  both  ways  finał y  would  give the  same  results.  H ere it  is  convenient to- use  F„  derived  by  the energetic  way.  Submatrix FB  derived  in  this  way  is  th e  following: ' n F ° = F as 1  n  • nx ~~ —  abc 27 hx0 0 4  2 8  4 8 0 0 4 2 4 8 Sym. 4 2 1 2 8 0 0 F«. 2 4 2 1 4 8 I 2 4 2 2 4 8 2 1 2 4 4 2 4 (10) — F = F If  instead  of  strains  e  we  introduce  stresses  N ,  for  instance  for  strains  s x : subm atrix  F„   becom es: F„°  - v?°„  - vF° n ­v?°n  F°  ~vF°„ ­ <  ~v?°„  F„° (12) In  the  same  way  submatrix  Foow  is  derived.  This  submatrix  can  be  defined  as  follows: ' w m v Fu 0 0 " 0 F„ 0 0 0 Fw. (13) 580 A .  POCESKI,  G .  KOKALANOV The  values  of  the  submatrices  take  the  form: Ftn — ac 4  2 ­ 2 ­ 4 2 1 ­ 1 ­ 2 2  4 ­ 4 ­ 2 1 2 ­ 2 ­ 1 2  1  ­ 1 ­ 2  4  2  ­ 2  ­ 4 1 2 ­ 2 ­ 1 2 4 ­ 4  ­ 2 . 4  ­ 4 ­ 2  2  2  ­ 2  ­ 1  l" 2 ­ 2 ­ 4 4 1 ­ 1 ­ 2 2 2 ­ 2 ­ 1 1 4 ­ 4 ­ 2 2 1  ­ 1 ­ 2  2  2 ­ 2  ­ 4  4. 4  2  1 2 ­ 4 ­ 2  ­ 2  ­ 1 2  4  2  1 ­ 2 ­ 4 ­ 1 ­ 2 1 2  4  2 ­ 1 ­ 2 ­ 2 ­ 4 2  1 2  4 ­ 2 ­ 1  ­ 4 ­ 2 (14) (15) (16) 2.2.2. Stiffness  submatrix  K. The stiffness  submatrix  K  gives the nodal  shear forces  in the direction  of  the  edge  node  displacements,  due  to  the  displacements  of  the  nodes.  For instance  the  forces  due  to  U9 =  1  are  as  follows: ­G 4b ­G (1­0, The  shear  force  due to the shear  stresses  acting  on the surface  1, 4, 8, 5 takes  the  form: 2a2c_ -rac -GT. Since the distribution  of the shear stresses from  node 9 to 11 is triangular, 2/3 of this shear force  is applied  to  node  9 and  1/3  to node  11. In this way the defined  nodal  forces  give the  coefficients  of the stiffness  submatrix  K. The submatrix can be represented  as follows: =  ~ |  K2  K23 [Sym.  K3 The  values  of  submatrices  in  this  expression  are  the  following: ac­2/?6  ­ac­fib 2(ac+/S6) (17) K1== Sym. (18) a  =  a/b,  /3 =  a/c, MIXED  FINITE  ELEMENTS 581 K1 3 2  1 ­ 1 2  ­ 1  ­ 2 2  ­ 2 2 2  1  ­ 2  ­ I " 2  ­ 1  ­ 2 2  1 . Sym.  2 b. (19) (20) Submatrix  K,  is  obtained  if  instead  of  «  in  K±  we  substitute  a" 1. K_9 ­i =r: 2 ­ 2  1 ­ 1 1 ­ 1  2 ­ 2 ­ 2  2 ­ 1  1 ­ 1  1 ­ 2  2 a, (21) (22) 2.3. Numerical examples. 2.3.1. Simply supported square plate.  The  plate is  subdivided  in  4  elements (Fig. 3a). Since there  is the  symmetry,  only  a  quarter  of  the  plate  is  analyzed,  (Fig. 3b). There are only  3 unknowns:  N6,  U10  and  W18.  The  analysis  gives  the  following  results: a ~ 1 ­ 2 a  —a"1 —a  —2a""1 + a 2(a~1 + a)  —2a"14­a  —a"1 —a 2(a~1 +  a)  a " 1 ­ 2 a Fig.  3. Analysis  of  simply  supported  square  plate  (a), the  quarter  of  the  plate  represented by  one  element  (b) l6Ec2  v ­3PL2 192D " " ° ' 0 1 5 6 2 (Thin  Plate  Th.0,0116PL2/D). 582 A.  POCESKI,  G.  KAKALANOV The  results  for  the  case  of  constant  distributed  load  q can  be  derived  by the  substitution of  qaz  instead  of  P/4.  The  results  are  as  follows: wia = 768D (Thin  plate  theory:  0,00406  qL4/D). The  following  bending  moment  corresponds  to  stress  N6: M6  =  N6 (1 +v)  =  0,0609«L 2,  fy  ­  0,3), (Thin  plate  theory:  0,0479 The  results  for  the  displacement  WiB  contain  two  components:  the  thin  plate  theory component  and  the  contribution  of  the  shear  forces  component.  The  fist  component is the  same  as is  the result  obtained  by  our  very first  rectangular  mixed bending  element [7], with the assumption  that the moments and displacements are independent. The second component  will  be  analyzed  in  the  next  chapter. 2.3.2. Cantilever.  The  analysis  of  the  cantilever  in  Fig.  4  loaded  by  concentrated  edge force  gives  the  following  results: 2bc2 P :  2Ec N6 Fig.  4.  Analysis  of  a  cantilever  as  one  finite  element These  results  are  exactly  the  beam  solution,  with  cross  section  shape  coefficient  k  =  1. As could  be exepected,  the  three­dimensional  element gives exactly the same results as the plane  stress  element  from  Fig.  1  [2]. In  the case of  a  cantilever  loaded  by a  moment,  the element  also  gives  also  the  beam  solution. MIXED  HNTTE  ELEMENTS 583 3.  Reduced  three­dimensional  element  for  plate  bending  analysis 3.1. Element matrix. The  reduced  three­dimensional  element  is presented  in Fig.  5.  The primary  unknowns  for  the element  are the  following: d< =  [Mlx  ...  M4X  ... Mly  ... MAy  ...  09  ... <912,  Wt  ... It  means  that  the  element  is  reduced  to  16 d.o.f. (23) Fig.  5.  Reduced  plate  bending  element The strains for the three­dimensional  element can be substituted  by the bending  moments as  follows: ex = ­v2) c D(l­vz) (Mx­vMy), (­vMx+My). The  displacements  u and v  can be  substituted  by  the rotations,  for  instance  according to  to  the following  relation: The upper stresses and displacements, and the lower stresses and displacements, in the case of  plate  bending  are of  the  same  intensity,  but  with  different  signes,  for  instance: U9 ­  ­  Uu The  element  matrix  can be  represented  as  follows: Fk = ­vFm FQy 0 0 JSym. (24) Submatrix  Fm derived in this  way is the same as for the element  with  the assumption  that the  moments and displacements are independent. The derivation  of the  other  submatrices is  very  simple  and will  not be given. 584 A .  POCESKI,  G .  KOKALANOV 3.2. Analysis of  thick  plates.  As  an  illustration  of  the  accuracy  of  the  reduced  three­ dimensional element, and consequently of the original element, the simply supported square plate in Fig. 3a will be analyzed. The results of the analysis are presented in the table below and  in  Figs,  6,  7.  , Table 1. Results of the analysis of thick simply supported square plate Thickness Span WL) 0.0125 0.075 0.125 0.200 0.250 Constant  distributed  load wtwS* 1.0059 1.031 1.077 1.187 1.290 M/M?x 1.054 1.054 1.054 1.054 1.054 Concent.  Force W1WS, 1.1055 1.160 1.256 1.492 1.709 The  results  presented  in  the  table  are  the  ratio  of  the  computed  deflection  W  at  the center  of the plate versus the theoretical  thin plate deflection,  and the similar ratio of the computed  M  and theoretical M:%x  moments. In  the case of very thin plates  (h/L  =  0.0125) the  element gives the  same results  as the  plate bending element with the assumption that the moments and deflections are independent. The moments do not depend on the thickness of  the  plate.  The  results  are  derived  by  subdivision  of  the  quarter  of  the  plate  on  2 x 2 elements. 0,05 0,1 t h e o r . UH [10] 6*6 n-147 0,15 0,2 0,25 PR [ 8 ] 6«6 n-204 — Mw 2»2 n - 8 7 Fig. 6.  Central  displacement  of  simply  supported  square thick  plate subjected  to  constant  distributed  load With  the  increase  of  the  thickness  of  the  plate  the  relative  deflections  of  the  plate increase  (Figs.  6  and  7).  Besides  the  results  obtained  by  our  element  Mw,  we  present in the figures the results obtained by Prior et al.  [8] (Pr) with stiffness  element, Chang­Chun Wu  [10]  (UH) —with  a  hybrid  element,  and  Rao  et  al.  [9]  (Ra) —with  a  triangular stiffness  element.  The  number  n  besides  the  particular  results  denotes  the  number  of equations  by  means  of  which  the  results  were  obtained. From  Fig.  6  one  can  see  that  our  mixed  element  (Mw)  gives  excellent  results,  with values  somewhat  higher  than  the  theoretical  ones,  while  the  stiffnes  element  (Pr)  gives MIXED  FINITE  ELEMENTS 585 similar accuracy, with values below the exact. The results in the case of concentrated  force (Fig. 7) take the values  between  those  obtained  by the other.  In region  of thin  plates the results are not so good,  since the rough  mesh was used  (2 x 2). If the refined  mesh  (4 x 4) is  used  then  the results  fall  down  close  to  the exact  ones  (doted  line  in  the  figure). We Mw 2' 2 n­87 Fig.  7.  Central  displacement  of  simply  supported  square  thick  plate  subjected  to  concentrated  force A  disadvantage  of  the stiffness  elements for  analysis  of thick  plates  is that  they  give bad results for thin plates. In the case of thin plates there is the so called  "locking" pheno­ menon.  Contrary  to that,  our mixed  element  gives good  results for  thick  and thin  plates, regardless  of  their  thickness. 4.  Further  development  and  conclusions The three­dimensional  element  presented  here  represents  our first  step in the develop­ ment  of  mixed  three­dimensional  elements.  Due to  its  shape,  the  practical  application of the element is limited. The purpose of the  development  of the element was to show the way  of  development  and  expected  accuracy  of  the  mixed  elements. The accuracy of the element is very good. In the case of plane stress problem the element gives the same results as the corresponding plane stress element  [2]. The plane stress element is one of the best elements  available at present. In the analysis of plate bending the three­ dimensional  element  gives the same results  as the plate  bending  element  with  assumption that the moments and displacements  are independent  [7]. The plate bending element  gives also  very  good  results. Such  results  obtained  by means  of the presented  element  are encouraging  for  further research  of  the three­dimensional  mixed  elements.  As the next  step  in  the  development of mixed  elements  should be the development  of a general  hexahedron  element  (Fig.  8a) with  36 d.o.f.  It  would  be convenient  if such an  element to be taken  with  displacements along  the edges  as the degrees  of freedom.  The  same  element  could  be  developed  with curved  boundaries.  The development  of the  element  could  be based  on the  assumption that the stresses and displacements  are independent. It was shown that such an assumption is correct and has advantage in the exclusion  of the parasitic  stresses. The final  aim would be  to  develope  such  elements  explicitly. The best mixed element that can be developed is the in Fig. 8b, with curved  boundaries and  60 d.o.f.  Such  an  element  corresponds  to  ourkplane  stress  isoparametric  element 2  Mech.  Teoret.  i  Stos.  4/88 586 A .  POCESKI,  G .  KOKALAN OV 1 + 8 ( N x , N y ) r s !z ) 9 *2 0 ( u , v, w) F ig.  8.  Possible  three- dimensional  elements:  (a)  with  36  d.o.f.  (b)  with  60  d.o.f. with  16 d.o.f.  [2]. This plane stress  element seems to be the best  plane stress  element  avail- able  at  present.  Thus,  the  corresponding  three- dimensional  element  should  be  expected to  give  excellent  results  also.  The  development  of  such  an  element  is  in  progress. The  reduced  three- dimensional  element  with  16  d.o.f.  succesfully  can  be  applied  for analysis  of  plate bending. The element gives results  of very good  accuracy, unconditionally stable,  regardless  the stiffness  of  the plate.  N ow  we  will  try  to  eliminate  the rotations  as d.o.f., although they are internal d.o.f., and in such a way  develope a 12 d.o.f. plate bending element  for  analysis  of  thick  plates. References 1.  R .  W.  C LO U Q H ,  Companion  of  T hree- dimensional  Finite Elements, Symp. Application  F E M   in  Civil. Engineering,  N ashville,  Ten n .  1969. 2.  A.  M .  P OCESKI,  G .  KOKALAN OV,  Mixed plane stress finite  elements,  I n t . Conf.  Computer  aided  Anal. D esign  C on crete  Str.,  Split,  1984,  Pineridge  press,  part  I,  707- 719. 3.  A.  P OCESKI,  G .  KOKALAN OV,  Mixed  plate  bending  elements, I n t .  Conf.  Comp.  Aidede  an al.  D esign C on e.  Str.,  Split  1984,  Pineridge  I ,  721 -  34. 4.  A.  POCESKI,  G .  KOKALAN OV,  Ploci- ramninska  sostojba  na  napregnajata,  primena na  MKE,  G radeź en F akultet,  Skopje  1981,  Izv.  2.3. 5.  G .  KOKALAN OV,  Mesoviti  izoparametrijski elementi za  analiza  na ploć i,  D r.  D is.,  G radezen  F acultet, Skopje,  1983. 6.  A.  POCESKI,  T he  direct method of  development  of finite  elements, J.  Theor.  App.  M echanics,  Belgrade 1985,  N o .  11. 7.  A.  POCESKI,  Mesovit  metod  na  konecni  elementi  (III),  12  Ju.  Kongres  Teor.  Prim.  M ehanike,  Ohrid 1974. 8.  C H .  P R YOR ,  R .  BAKER,  D .  F RED ERIC, Finite  element bending analysis  of  Reisner's plates,  ASCE  96, N o .  E M 6,  1970,  967  -  Sd. 9.  G .  R AO  et al., A  high precision triangular plate  bending element for  analysis of  thick plates, N uclear  Eng. D esign  30,  1974,  408  -  12. 10.  Wu.  C H AN G - C H U N,  Some  problems  of  a plate  bending  hybrid model  with shear effects, I n t .  J.  N um . M eth .  E n g.  18,  1982,  755 -  64. M I XE D   F I N I TE  ELEMENTS  587 P  e  3  K>  M  e T P fiXM E P H BI E  C M EIU AH H BIE KOH E^LH BIE 3J I E M E H T B I n apaiueTpw  H BJI H I OTC JI  Hen3BecTHbiMH  B  ojieiueH Te:  KOMnoHenTH  H an pji> Kem n ł   B  yrcro- Bbix  y3Jiax  ( 8 x 3 )  H  n epeM em ein ra:  r p a m n i  ( 1 2 x 1 ) .  SjieineH T  pa3BHM>rii  n p a  n oM om ji  H en ocpe«crrBeH - H oro  iweTo.ua,  6e3  npHMeHenHH   BapiiauH OH H oro  n p i n m i i n a .  H 3BecrH bie  3jnapa3H TH Bie  H a u p jjjK e in r a " HcraiKMaioTCH   aBTOMaTH^recKH.  C Beflem ibiS  TpexM epH tift  ajieiweHT  flaeT  oyeH B  x o p o m n e  pe3yjibTaTBi B  aH ajiH 3e  TOJICTLIX  H  TOH KH X  njiH T.  B  p a S o ie  noKa3aH   nyTŁ  pa3BH TH «  cjien iaH H tix  ajiejueiiTOB  u  OH < H - flaeM an  To^raocrb.  flanee  paccMOTpeH bi  H 3onapaMeTpH tiecKH e  3JieMeHTW  J I I O SM X  cbopM   ( p u c . S3 a - 6 ). H aH jiyquiH M   o>KHflaeMMM  ojieivieHTOM   H BJI H C TC H  noKa3aH biH  Ha p u c .  86 aJieiweirr  c 60 CTeneHHMH   CBo6oflbi, KOTopbifi  T e n e p t  pa3pa6oTH BaeTCfl:. S t r e s z c z e n i e TR ÓJWYM I AR OWE  M IESZAN E  E LE M E N T Y SK O Ń C Z O NE N ieznanymi  param etram i  elementu  są   skł adowe  naprę ż enia  w  wę zł ach  naroż nych  ( 8x3)  i  przemiesz- czenia  krawę dzi  ( 12x1) .  Elementy  został y  rozwinię te  z  pomocą   metody  bezpoś redniej  bez  zastosowania zasady  wariacyjnej.  „ N aprę ż enia  param etryczne"  został y  wyeliminowane  automatycznie.  Z redukowane trójwymiarowe  elementy  dają   znakom ite  wyniki  w  analizie  pł yt  grubych  i  cienkich. Celem  pracy  jest  pokazanie  sposobu  rozwoju  mieszanych  elementów  i  spodziewanej  dokł adnoś ci ich  stosowania.  Oprócz  tego rozpatrujemy  izoparametryczne  elementy  o  dowolnych  kształ tach  (rys.  8a -  b). N ajlepszym  jest  element  pokazany  n a  rys.  8b  z  60  stopniami  swobody.  Jest  on  obecnie  rozwijany. Praca  wpł ynę ł a do  Redakcji  dnia  29  czerwca  1987  roku.