Chapter 8

 

Theory and Implementation of the Finite Element Method

 

 

 

8.6 Advanced element formulations MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqabKqzaiaeaa aaaaaaa8qacaWFtacaaa@3788@  Incompatible modes; reduced integration; and hybrid elements

 

Techniques for interpolating the displacement field within 2D and 3D finite elements were discussed in Section 8.1.9 and 8.1.10.  In addition, methods for evaluating the volume or area integrals in the principle of virtual work were discussed in Section 8.1.11.   These procedures work well for most applications, but there are situations where the simple element formulations can give very inaccurate results.  In this section

  1. We illustrate some of the unexpected difficulties that can arise in apparently perfectly well designed finite element solutions to boundary value problems;
  2. We describe a few more sophisticated elements that have been developed to solve these problems.

 

We focus in particular on `Locking’ phenomena.  Finite elements are said to `lock’ if they exhibit an unphysically stiff response to deformation.  Locking can occur for many different reasons. The most common causes are (i) the governing equations you are trying to solve are poorly conditioned, which leads to an ill conditioned system of finite element equations; (ii) the element interpolation functions are unable to approximate accurately the strain distribution in the solid, so the solution converges very slowly as the mesh size is reduced; (iii) in certain element formulations (especially beam, plate and shell elements) displacements and their derivatives are interpolated separately. Locking can occur in these elements if the interpolation functions for displacements and their derivatives are not consistent.

 

 

8.6.1 Shear locking and incompatible mode elements

 

Shear locking can be illustrated by attempting to find a finite element solution to the simple boundary value problem illustrated in the picture. Consider a cantilever beam, with length L, height 2a and out-of-plane thickness b, as shown in the figure.  The top and bottom of the beam x 2 =±a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacaWG4bWaaSbaaSqaaiaaikdaaeqaaO Gaeyypa0JaeyySaeRaamyyaaaa@388F@  are traction free, the left hand end is subjected to a resultant force P, and the right hand end is clamped.  Assume that b<<a, so that a state of plane stress is developed in the beam.  The analytical solution to this problem is given in Section 5.2.4.

 

The figures below compare this result to a finite element solution, obtained with standard 4 noded linear quadrilateral plane stress elements.  Results are shown for two different ratios of a/L MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rrps0l bbf9q8WrFfeuY=Hhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacaWGHbGaai4laiaadYeaaaa@3419@ , with P L 3 /E a 3 b=2.22 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadcfacaWGmbWaaWbaaSqabeaacaaIZa aaaOGaai4laiaadweacaWGHbWaaWbaaSqabeaacaaIZaaaaOGaamOy aiabg2da9iaaikdacaGGUaGaaGOmaiaaikdaaaa@3B2C@  for both cases. 

   

 

For the thick beam, the finite element and exact solutions agree nearly perfectly.   For the thin beam, however, the finite element solution is very poor MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  even though the mesh resolution is unchanged. The error in the finite element solution occurs because the standard 4 noded quadrilateral elements cannot accurately approximate the strain distribution associated with bending.  The phenomenon is known as `shear locking’ because the element interpolation functions give rise to large, unphysical shear strains in bent elements.  The solution would eventually converge if a large number of elements were added along the length of the beam. The elements would have to be roughly square, which would require about 133 elements along the length of the beam, giving a total mesh size of about 500 elements. 

 

Shear locking is therefore relatively benign, since it can be detected by refining the mesh, and can be avoided by using a sufficiently fine mesh.  However, finite element analysts sometimes cannot resist the temptation to reduce computational cost by using elongated elements, which can introduce errors.

 

Shear locking can also be avoided by using more sophisticated element interpolation functions that can accurately approximate bending. `Incompatible Mode’ elements do this by adding an additional strain distribution to the element.  The elements are called `incompatible’ because the strain is not required to be compatible with the displacement interpolation functions. The approach is conceptually straightforward:

1.      The displacement fields in the element are interpolated using the standard scheme, by setting

u i = a=1 N e N a ( ξ j ) u i a δ v i (x)= a=1 n N a (x)δ v i a x i = a=1 N e N a ( ξ j ) x i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaBaaaleaacaWGPbaabeaaki abg2da9maaqahabaGaamOtamaaCaaaleqabaGaamyyaaaakiaacIca cqaH+oaEdaWgaaWcbaGaamOAaaqabaGccaGGPaGaamyDamaaDaaale aacaWGPbaabaGaamyyaaaaaeaacaWGHbGaeyypa0JaaGymaaqaaiaa d6eadaWgaaadbaGaamyzaaqabaaaniabggHiLdGccaaMc8UaaGPaVl aaykW7caaMc8UaeqiTdqMaamODamaaBaaaleaacaWGPbaabeaakiaa cIcacaWH4bGaaiykaiabg2da9maaqahabaGaamOtamaaCaaaleqaba GaamyyaaaakiaacIcacaWH4bGaaiykaiabes7aKjaadAhadaqhaaWc baGaamyAaaqaaiaadggaaaaabaGaamyyaiabg2da9iaaigdaaeaaca WGUbaaniabggHiLdGccaaMc8UaaGPaVlaaykW7caWG4bWaaSbaaSqa aiaadMgaaeqaaOGaeyypa0ZaaabCaeaacaWGobWaaWbaaSqabeaaca WGHbaaaOGaaiikaiabe67a4naaBaaaleaacaWGQbaabeaakiaacMca caWG4bWaa0baaSqaaiaadMgaaeaacaWGHbaaaaqaaiaadggacqGH9a qpcaaIXaaabaGaamOtamaaBaaameaacaWGLbaabeaaa0GaeyyeIuoa aaa@77E6@  

where N a ( ξ j ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOtamaaCaaaleqabaGaamyyaaaaki aacIcacqaH+oaEdaWgaaWcbaGaamOAaaqabaGccaGGPaaaaa@3700@  are the shape functions listed in Sections 8.1.9 or 8.1.10, ξ i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaH+oaEdaWgaaWcbaGaamyAaaqaba aaaa@35A3@  are a set of local coordinates in the element, u i a , x i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadwhadaqhaaWcbaGaamyAaaqaaiaadg gaaaGccaGGSaGaamiEamaaDaaaleaacaWGPbaabaGaamyyaaaaaaa@373B@  denote the displacement values and coordinates of the nodes on the element, and N e MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaad6eadaWgaaWcbaGaamyzaaqabaaaaa@3271@  is the number of nodes on the element.

2.      The Jacobian matrix for the interpolation functions, its determinant, and its inverse are defined in the usual way

x i ξ j = a=1 N e N a (ξ) ξ j x i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaSaaaeaacqGHciITcaWG4bWaaSbaaS qaaiaadMgaaeqaaaGcbaGaeyOaIyRaeqOVdG3aaSbaaSqaaiaadQga aeqaaaaakiabg2da9maaqahabaWaaSaaaeaacqGHciITcaWGobWaaW baaSqabeaacaWGHbaaaOGaaiikaiabe67a4jaacMcaaeaacqGHciIT cqaH+oaEdaWgaaWcbaGaamOAaaqabaaaaOGaamiEamaaDaaaleaaca WGPbaabaGaamyyaaaaaeaacaWGHbGaeyypa0JaaGymaaqaaiaad6ea daWgaaadbaGaamyzaaqabaaaniabggHiLdaaaa@4E50@     J=det( x i ξ j ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOsaiabg2da9iGacsgacaGGLbGaai iDamaabmaabaWaaSaaaeaacqGHciITcaWG4bWaaSbaaSqaaiaadMga aeqaaaGcbaGaeyOaIyRaeqOVdG3aaSbaaSqaaiaadQgaaeqaaaaaaO GaayjkaiaawMcaaaaa@3EDD@          ξ j x i = ( x ξ ) ji 1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaSaaaeaacqGHciITcqaH+oaEdaWgaa WcbaGaamOAaaqabaaakeaacqGHciITcaWG4bWaaSbaaSqaaiaadMga aeqaaaaakiabg2da9maabmaabaWaaSaaaeaacqGHciITcaWG4baaba GaeyOaIyRaeqOVdGhaaaGaayjkaiaawMcaamaaDaaaleaacaWGQbGa amyAaaqaaiabgkHiTiaaigdaaaaaaa@4491@

3.      The usual expression for displacement gradient in the element is replaced by

u i x j =( a=1 N e N a (ξ) ξ m u i a + k=1 p J(0) J(ξ) α i (k) δ km ξ k ) ξ m x j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaadaWcaaqaaiabgkGi2kaadwhadaWgaa WcbaGaamyAaaqabaaakeaacqGHciITcaWG4bWaaSbaaSqaaiaadQga aeqaaaaakiabg2da9maabmaabaWaaabCaeaadaWcaaqaaiabgkGi2k aad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGaeqOVdGNaaiykaaqa aiabgkGi2kabe67a4naaBaaaleaacaWGTbaabeaaaaaabaGaamyyai abg2da9iaaigdaaeaacaWGobWaaSbaaWqaaiaadwgaaeqaaaqdcqGH ris5aOGaamyDamaaDaaaleaacaWGPbaabaGaamyyaaaakiabgUcaRm aaqahabaWaaSaaaeaacaWGkbGaaiikaiaaicdacaGGPaaabaGaamOs aiaacIcacqaH+oaEcaGGPaaaaaWcbaGaam4Aaiabg2da9iaaigdaae aacaWGWbaaniabggHiLdGccqaHXoqydaqhaaWcbaGaamyAaaqaaiaa cIcacaWGRbGaaiykaaaakiabes7aKnaaBaaaleaacaWGRbGaamyBaa qabaGccqaH+oaEdaWgaaWcbaGaam4AaaqabaaakiaawIcacaGLPaaa daWcaaqaaiabgkGi2kabe67a4naaBaaaleaacaWGTbaabeaaaOqaai abgkGi2kaadIhadaWgaaWcbaGaamOAaaqabaaaaaaa@7254@

where p=2 for a 2D problem and p=3 for a 3D problem, α i (k) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaHXoqydaqhaaWcbaGaamyAaaqaai aacIcacaWGRbGaaiykaaaaaaa@37C9@  are a set of unknown displacement gradients in the element, which must be determined as part of the solution.

4.      Similarly, the virtual displacement gradient is written as

δ v i x j =( a=1 N e N a (ξ) ξ m δ v i a + k=1 p J(0) J(ξ) δ α i (k) δ km ξ k ) ξ m x j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaadaWcaaqaaiabgkGi2kabes7aKjaadA hadaWgaaWcbaGaamyAaaqabaaakeaacqGHciITcaWG4bWaaSbaaSqa aiaadQgaaeqaaaaakiabg2da9maabmaabaWaaabCaeaadaWcaaqaai abgkGi2kaad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGaeqOVdGNa aiykaaqaaiabgkGi2kabe67a4naaBaaaleaacaWGTbaabeaaaaaaba Gaamyyaiabg2da9iaaigdaaeaacaWGobWaaSbaaWqaaiaadwgaaeqa aaqdcqGHris5aOGaeqiTdqMaamODamaaDaaaleaacaWGPbaabaGaam yyaaaakiabgUcaRmaaqahabaWaaSaaaeaacaWGkbGaaiikaiaaicda caGGPaaabaGaamOsaiaacIcacqaH+oaEcaGGPaaaaaWcbaGaam4Aai abg2da9iaaigdaaeaacaWGWbaaniabggHiLdGccqaH0oazcqaHXoqy daqhaaWcbaGaamyAaaqaaiaacIcacaWGRbGaaiykaaaakiabes7aKn aaBaaaleaacaWGRbGaamyBaaqabaGccqaH+oaEdaWgaaWcbaGaam4A aaqabaaakiaawIcacaGLPaaadaWcaaqaaiabgkGi2kabe67a4naaBa aaleaacaWGTbaabeaaaOqaaiabgkGi2kaadIhadaWgaaWcbaGaamOA aaqabaaaaaaa@7744@

where δ α i (k) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaH0oazcqaHXoqydaqhaaWcbaGaam yAaaqaaiaacIcacaWGRbGaaiykaaaaaaa@396E@  is a variation in the internal displacement gradient field for the element.  

5.      These expressions are then substituted into the virtual work equation, which must now be satisfied for all possible values of virtual nodal displacements δ v i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaH0oazcaWG2bWaa0baaSqaaiaadM gaaeaacaWGHbaaaaaa@3767@  and virtual displacement gradients δ α i (k) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaH0oazcqaHXoqydaqhaaWcbaGaam yAaaqaaiaacIcacaWGRbGaaiykaaaaaaa@396E@ .  At first sight this procedure appears to greatly increase the size of the global stiffness matrix, because a set of unknown displacement gradient components must be calculated for each element.  However, the unknown α i (k) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaHXoqydaqhaaWcbaGaamyAaaqaai aacIcacaWGRbGaaiykaaaaaaa@37C9@  are local to each element, and can be eliminated while computing the element stiffness matrix.  The procedure to do this can be shown most clearly in a sample code.

 

A sample small-strain, linear elastic code with incompatible mode elements is provided in Femlinelast_incompatible_modes.mws.  When run with the input file shear_locking_demo.txt it produces the results shown in the figure.  The incompatible modes clearly give a spectacular improvement in the performance of the element.

 

HEALTH WARNINGS: (i) The procedure outlined here only works for small-strain problems.  Finite strain versions exist but are somewhat more complicated.  (ii)  Adding strain variables to elements can dramatically improve their performance, but this procedure must be used with great care to ensure that the strain and displacement degrees of freedom are independent variables. For further details see Simo,  J. C., and M. S. Rifai, Int J. Num. Meth in Eng.  29, pp. 1595 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@ 1638, 1990. Simo,  J. C., and F. Armero, Int J. Num. Meth in Eng., 33, pp. 1413 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@ 1449, 1992.

 

 

8.6.2 Volumetric locking and reduced integration elements

 

Volumetric locking can be illustrated using a simple boundary value problem. Consider a long hollow cylinder with internal radius a and external radius b as shown in the figure.  The solid is made from a linear elastic material with Young’s modulus E MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacaWGfbaaaa@3390@  and Poisson’s ratio ν MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaH9oGBaaa@347E@ .  The cylinder is loaded by an internal pressure p a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacaWGWbWaaSbaaSqaaiaadggaaeqaaa aa@34CD@  and deforms in plane strain. 

The analytical solution to this problem is given in Section 4.1.9.

 

The figures below compare the analytical solution to a finite element solution with standard 4 noded plane strain quadrilateral elements.  Results are shown for two values of Poisson’s ratio ν MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaH9oGBaaa@347E@ .  The dashed lines show the analytical solution, while the solid line shows the FEA solution.

 

           

The two solutions agree well for ν=0.3 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaH9oGBcqGH9aqpcaaIWaGaaiOlai aaiodaaaa@37AD@ , but the finite element solution grossly underestimates the displacements as Poisson’s ratio is increased towards 0.5 (recall that the material is incompressible in the limit ν=0.5 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaH9oGBcqGH9aqpcaaIWaGaaiOlai aaiwdaaaa@37AF@  ).  In this limit, the finite element displacements tend to zero MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  this is known as `volumetric locking’

 

The error in the finite element solution occurs because the finite element interpolation functions are unable to properly approximate a volume preserving strain field. In the incompressible limit, a nonzero volumetric strain at any of the integration points gives rise to a very large contribution to the virtual power.  The interpolation functions can make the volumetric strain vanish at some, but not all, the integration points in the element. 

 

Volumetric locking is a much more serious problem than shear locking, because it cannot be avoided by refining the mesh.  In addition, all the standard fully integrated finite elements will lock in the incompressible limit; and some elements show very poor performance even for Poisson’s ratios as small as 0.45 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacaaIWaGaaiOlaiaaisdacaaI1aaaaa@35AF@ .  Fortunately, most materials have Poisson’s ratios around 0.3 or less, so the standard elements can be used for most linear elasticity and small-strain plasticity problems.   To model rubbers, or to solve problems involving large plastic strains, the elements must be redesigned to avoid locking.

 

Reduced Integration is the simplest way to avoid locking.   The basic idea is simple: since the fully integrated elements cannot make the strain field volume preserving at all the integration points, it is tempting to reduce the number of integration points so that the constraint can be met.  `Reduced integration’ usually means that the element stiffness is integrated using an integration scheme that is one order less accurate than the standard scheme.   The number of reduced integration points for various element types is listed in the table below.  The coordinates of the integration points are listed in the tables in Section 8.1.12.

 

Number of integration points for reduced integration schemes

Linear triangle (3 nodes) 1 point

Quadratic triangle (6 nodes): 3 points

Linear quadrilateral (4 nodes): 1 point

Quadratic quadrilateral (8 nodes): 4 points

Linear tetrahedron (4 nodes): 1 point

Quadratic tetrahedron (10 nodes): 4 points

Linear brick (8 nodes): 1 point

Quadratic brick (20 nodes): 8 points

 

HEALTH WARNING: Notice that the integration order cannot be reduced for the linear triangular and tetrahedral elements MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  these elements should not be used to model near incompressible materials, although in desperation you can a few such elements in regions where the solid cannot be meshed using quadrilaterals.

 

Remarkably, reduced integration completely resolves locking in some elements (the quadratic quadrilateral and brick), and even improves the accuracy of the element.  As an example, the figure below shows the solution to the pressurized cylinder problem, using both full and reduced integration for 8 noded quadrilaterals.  With reduced integration, the analytical and finite element results are indistinguishable.

                

 

 

Reduced integration does not work in 4 noded quadrilateral elements or 8 noded brick elements.  For example, the figure on the right shows the solution to the pressure vessel problem with linear (4 noded) quadrilateral elements with reduced integration.  The solution is clearly a disaster.  The error occurs because the stiffness matrix is nearly singular MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@   the system of equations includes a weakly constrained deformation mode. This   phenomenon is known as `hourglassing’ because of the characteristic shape of the spurious deformation mode. 

 

Selectively Reduced Integration can be used to cure hourglassing.   The procedure is illustrated most clearly by modifying the formulation for static linear elasticity. To implement the method: 

1.      The volume integral in the virtual work principle is separated into a deviatoric and volumetric part by writing

R σ ij δ ε ij d V 0 = R ( σ ij δ ε ij σ kk 3 δ ε qq )d V 0 + R σ kk 3 δ ε qq d V 0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaKaaGgaacqaHdpWCkmaaBaaale aacaWGPbGaamOAaaqabaGccqaH0oazcqaH1oqzdaWgaaWcbaGaamyA aiaadQgaaeqaaKaaGkaaykW7caWGKbGaamOvaOWaaSbaaKqaGfaaca aIWaaabeaaaeaacaWGsbaabeqcdaMaey4kIipajaaycqGH9aqpkmaa pefajaaObaGcdaqadaqaaiabeo8aZnaaBaaaleaacaWGPbGaamOAaa qabaGccqaH0oazcqaH1oqzdaWgaaWcbaGaamyAaiaadQgaaeqaaOGa eyOeI0YaaSaaaeaacqaHdpWCdaWgaaWcbaGaam4AaiaadUgaaeqaaa GcbaGaaG4maaaacqaH0oazcqaH1oqzdaWgaaWcbaGaamyCaiaadgha aeqaaaGccaGLOaGaayzkaaqcaaQaaGPaVlaadsgacaWGwbGcdaWgaa qcbawaaiaaicdaaeqaaaqaaiaadkfaaeqajmaycqGHRiI8aOGaey4k aSYaa8quaKaaGgaakmaalaaajaaObaGaeq4WdmNcdaWgaaqcbaAaai aadUgacaWGRbaabeaaaKaaGgaacaaIZaaaaiabes7aKjabew7aLPWa aSbaaKqaGgaacaWGXbGaamyCaaqabaqcaaQaaGPaVlaadsgacaWGwb GcdaWgaaqcbawaaiaaicdaaeqaaaqaaiaadkfaaeqajmaycqGHRiI8 aaaa@7B0F@

 Here, the first integral on the right hand side vanishes for a hydrostatic stress.

2.      Substituting the linear elastic constitutive equation and the finite element interpolation functions into the virtual work principle, we find that the element stiffness matrix can be reduced to

k aibk (l) = V e (l) ( C ijkl N a (x) x j N b (x) x l 1 3 C ppkl N a (x) x i N b (x) x l )dV + V e (l) 1 3 C ppkl N a (x) x i N b (x) x l dV MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4AamaaDaaaleaacaWGHbGaamyAai aadkgacaWGRbaabaGaaiikaiaadYgacaGGPaaaaOGaeyypa0Zaa8qu aeaadaqadaqaaiaadoeadaWgaaWcbaGaamyAaiaadQgacaWGRbGaam iBaaqabaGcdaWcaaqaaiabgkGi2kaad6eadaahaaWcbeqaaiaadgga aaGccaGGOaGaaCiEaiaacMcaaeaacqGHciITcaWG4bWaaSbaaSqaai aadQgaaeqaaaaakmaalaaabaGaeyOaIyRaamOtamaaCaaaleqabaGa amOyaaaakiaacIcacaWH4bGaaiykaaqaaiabgkGi2kaadIhadaWgaa WcbaGaamiBaaqabaaaaOGaeyOeI0YaaSaaaeaacaaIXaaabaGaaG4m aaaacaWGdbWaaSbaaSqaaiaadchacaWGWbGaam4AaiaadYgaaeqaaO WaaSaaaeaacqGHciITcaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiik aiaahIhacaGGPaaabaGaeyOaIyRaamiEamaaBaaaleaacaWGPbaabe aaaaGcdaWcaaqaaiabgkGi2kaad6eadaahaaWcbeqaaiaadkgaaaGc caGGOaGaaCiEaiaacMcaaeaacqGHciITcaWG4bWaaSbaaSqaaiaadY gaaeqaaaaaaOGaayjkaiaawMcaaiaadsgacaWGwbaaleaacaWGwbWa a0baaWqaaiaadwgaaeaacaGGOaGaamiBaiaacMcaaaaaleqaniabgU IiYdGccaaMc8UaaGPaVlabgUcaRmaapefabaWaaSaaaeaacaaIXaaa baGaaG4maaaacaWGdbWaaSbaaSqaaiaadchacaWGWbGaam4AaiaadY gaaeqaaOWaaSaaaeaacqGHciITcaWGobWaaWbaaSqabeaacaWGHbaa aOGaaiikaiaahIhacaGGPaaabaGaeyOaIyRaamiEamaaBaaaleaaca WGPbaabeaaaaGcdaWcaaqaaiabgkGi2kaad6eadaahaaWcbeqaaiaa dkgaaaGccaGGOaGaaCiEaiaacMcaaeaacqGHciITcaWG4bWaaSbaaS qaaiaadYgaaeqaaaaakiaadsgacaWGwbaaleaacaWGwbWaa0baaWqa aiaadwgaaeaacaGGOaGaamiBaiaacMcaaaaaleqaniabgUIiYdaaaa@9949@

3.      When selectively reduced integration is used, the first volume integral is evaluated using the full integration scheme; the second integral is evaluated using reduced integration points.

 

Selective reduced integration has been implemented in the sample program fem_selective_reduced_integration.mws.  When this code is run with the input file volumetric_locking_demo.txt it produces the results shown in the figure.  The analytical and finite element solutions agree, and there are no signs of hourglassing.

 

In many commercial codes, the `fully integrated’ elements actually use selective reduced integration.

 

 

The ‘B-bar’ method:  Like selective reduced integration, the B-bar method works by treating the volumetric and deviatoric parts of the stiffness matrix separately.  Instead of separating the volume integral into two parts, however, the B-bar method modifies the definition of the strain in the element.  Its main advantage is that the concept can easily be generalized to finite strain problems.   Here, we will illustrate the method by applying it to small-strain linear elasticity.  The procedure starts with the usual virtual work principle 

R σ ij [ ε ij ( u k )]δ ε ij dV R b i δ v i dV 2 R t i * δ v i dA=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaeaacqaHdpWCdaWgaaWcbaGaam yAaiaadQgaaeqaaOGaai4waiabew7aLnaaBaaaleaacaWGPbGaamOA aaqabaGccaGGOaGaamyDamaaBaaaleaacaWGRbaabeaakiaacMcaca GGDbGaeqiTdqMaeqyTdu2aaSbaaSqaaiaadMgacaWGQbaabeaakiaa dsgacaWGwbGaeyOeI0Yaa8quaeaacaWGIbWaaSbaaSqaaiaadMgaae qaaOGaeqiTdqMaamODamaaBaaaleaacaWGPbaabeaakiaadsgacaWG wbGaeyOeI0Yaa8quaeaacaWG0bWaa0baaSqaaiaadMgaaeaacaGGQa aaaOGaeqiTdqMaamODamaaBaaaleaacaWGPbaabeaakiaadsgacaWG bbGaeyypa0JaaGimaaWcbaGaeyOaIy7aaSbaaWqaaiaaikdaaeqaaS GaamOuaaqab0Gaey4kIipaaSqaaiaadkfaaeqaniabgUIiYdaaleaa caWGsbaabeqdcqGHRiI8aaaa@63E4@

In the B-bar method

1.      We introduce a new variable to characterize the volumetric strain in the elements.  Define

ω= 1 3 V e V e ε kk dV = B bk (vol) u k b B bk (vol) = 1 3 V e V e N b x k dV MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=fYJH8sipiYdHaVhbbf9v8qqaqFr0xc9pk0xbb a9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9 Ff0dmeaabaqaciGacaGaaeqabaWaaeaaeaaakeaacqaHjpWDcqGH9a qpdaWcaaqaaiaaigdaaeaacaaIZaGaamOvamaaBaaaleaacaWGLbaa beaaaaGcdaWdrbqaaiabew7aLnaaBaaaleaacaWGRbGaam4Aaaqaba GccaWGKbGaamOvaaWcbaGaamOvamaaBaaameaacaWGLbaabeaaaSqa b0Gaey4kIipakiabg2da9iaadkeadaqhaaWcbaGaamOyaiaadUgaae aacaGGOaGaamODaiaad+gacaWGSbGaaiykaaaakiaadwhadaqhaaWc baGaam4AaaqaaiaadkgaaaGccaaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWGcbWaa0baaSqaai aadkgacaWGRbaabaGaaiikaiaadAhacaWGVbGaamiBaiaacMcaaaGc cqGH9aqpdaWcaaqaaiaaigdaaeaacaaIZaGaamOvamaaBaaaleaaca WGLbaabeaaaaGcdaWdrbqaamaalaaabaGaeyOaIyRaamOtamaaCaaa leqabaGaamOyaaaaaOqaaiabgkGi2kaadIhadaWgaaWcbaGaam4Aaa qabaaaaOGaamizaiaadAfaaSqaaiaadAfadaWgaaadbaGaamyzaaqa baaaleqaniabgUIiYdaaaa@845B@

where the integral is taken over the volume of the element.

2.      The strain variation in each element is replaced by the approximation ε ¯ ij = ε ij +(ω ε kk /3) δ ij MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=fYJH8sipiYdHaVhbbf9v8qqaqFr0xc9pk0xbb a9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9 Ff0dmeaabaqaciGacaGaaeqabaWaaeaaeaaakeaacuaH1oqzgaqeam aaBaaaleaacaWGPbGaamOAaaqabaGccqGH9aqpcqaH1oqzdaWgaaWc baGaamyAaiaadQgaaeqaaOGaey4kaSIaaiikaiabeM8a3jabgkHiTi abew7aLnaaBaaaleaacaWGRbGaam4AaaqabaGccaGGVaGaaG4maiaa cMcacqaH0oazdaWgaaWcbaGaamyAaiaadQgaaeqaaaaa@4B22@

3.      Similarly, the virtual strain in each element is replaced by δ ε ¯ ij =δ ε ij +(δωδ ε kk /3) δ ij MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=fYJH8sipiYdHaVhbbf9v8qqaqFr0xc9pk0xbb a9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9 Ff0dmeaabaqaciGacaGaaeqabaWaaeaaeaaakeaacqaH0oazcuaH1o qzgaqeamaaBaaaleaacaWGPbGaamOAaaqabaGccqGH9aqpcqaH0oaz cqaH1oqzdaWgaaWcbaGaamyAaiaadQgaaeqaaOGaey4kaSIaaiikai abes7aKjabeM8a3jabgkHiTiabes7aKjabew7aLnaaBaaaleaacaWG RbGaam4AaaqabaGccaGGVaGaaG4maiaacMcacqaH0oazdaWgaaWcba GaamyAaiaadQgaaeqaaaaa@51B6@

This means that the volumetric strain in the element is everywhere equal to its mean value.  The virtual work principle is then written in terms of ε ¯ ij MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=fYJH8sipiYdHaVhbbf9v8qqaqFr0xc9pk0xbb a9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9 Ff0dmeaabaqaciGacaGaaeqabaWaaeaaeaaakeaacuaH1oqzgaqeam aaBaaaleaacaWGPbGaamOAaaqabaaaaa@3888@  and δ ε ¯ ij MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=fYJH8sipiYdHaVhbbf9v8qqaqFr0xc9pk0xbb a9q8WqFfea0=yr0RYxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9 Ff0dmeaabaqaciGacaGaaeqabaWaaeaaeaaakeaacqaH0oazcuaH1o qzgaqeamaaBaaaleaacaWGPbGaamOAaaqabaaaaa@3A2D@  as

R σ ij [ ε ¯ ij ( u k )]δ ε ¯ ij dV R b i δ v i dV 2 R t i * δ v i dA=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaeaacqaHdpWCdaWgaaWcbaGaam yAaiaadQgaaeqaaOGaai4waiqbew7aLzaaraWaaSbaaSqaaiaadMga caWGQbaabeaakiaacIcacaWG1bWaaSbaaSqaaiaadUgaaeqaaOGaai ykaiaac2facqaH0oazcuaH1oqzgaqeamaaBaaaleaacaWGPbGaamOA aaqabaGccaWGKbGaamOvaiabgkHiTmaapefabaGaamOyamaaBaaale aacaWGPbaabeaakiabes7aKjaadAhadaWgaaWcbaGaamyAaaqabaGc caWGKbGaamOvaiabgkHiTmaapefabaGaamiDamaaDaaaleaacaWGPb aabaGaaiOkaaaakiabes7aKjaadAhadaWgaaWcbaGaamyAaaqabaGc caWGKbGaamyqaiabg2da9iaaicdaaSqaaiabgkGi2oaaBaaameaaca aIYaaabeaaliaadkfaaeqaniabgUIiYdaaleaacaWGsbaabeqdcqGH RiI8aaWcbaGaamOuaaqab0Gaey4kIipaaaa@6414@

Finally, introducing the finite element interpolation and using the constitutive equation yields the usual system of linear finite element equations, but with a modified element stiffness matrix given by

k aibk (l) = V e (l) C pjql B ¯ pji a B ¯ qlk b dV B ¯ pji a = N a x j δ ip +( B ia (vol) 1 3 N a x i ) δ pj MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4AamaaDaaaleaacaWGHbGaamyAai aadkgacaWGRbaabaGaaiikaiaadYgacaGGPaaaaOGaeyypa0Zaa8qu aeaacaWGdbWaaSbaaSqaaiaadchacaWGQbGaamyCaiaadYgaaeqaaO GabmOqayaaraWaa0baaSqaaiaadchacaWGQbGaamyAaaqaaiaadgga aaGcceWGcbGbaebadaqhaaWcbaGaamyCaiaadYgacaWGRbaabaGaam OyaaaakiaadsgacaWGwbaaleaacaWGwbWaa0baaWqaaiaadwgaaeaa caGGOaGaamiBaiaacMcaaaaaleqaniabgUIiYdGccaaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UabmOqayaaraWaa0baaSqaaiaadcha caWGQbGaamyAaaqaaiaadggaaaGccqGH9aqpdaWcaaqaaiabgkGi2k aad6eadaahaaWcbeqaaiaadggaaaaakeaacqGHciITcaWG4bWaaSba aSqaaiaadQgaaeqaaaaakiabes7aKnaaBaaaleaacaWGPbGaamiCaa qabaGccqGHRaWkdaqadaqaaiaadkeadaqhaaWcbaGaamyAaiaadgga aeaacaGGOaGaamODaiaad+gacaWGSbGaaiykaaaakiabgkHiTmaala aabaGaaGymaaqaaiaaiodaaaWaaSaaaeaacqGHciITcaWGobWaaWba aSqabeaacaWGHbaaaaGcbaGaeyOaIyRaamiEamaaBaaaleaacaWGPb aabeaaaaaakiaawIcacaGLPaaacqaH0oazdaWgaaWcbaGaamiCaiaa dQgaaeqaaaaa@98DA@

This expression can be integrated using a standard, full integration scheme.

 

The B-bar method has been implemented in the sample code FEM_Bbar.mws.  The code can be run with the input file volumetric_locking_demo.txt.  Run the code for yourself to verify that the analytical and finite element solutions agree, and there are no signs of hourglassing.

Hourglass base vectors

Linear quadrilateral

Γ a(1) =(+1,1,+1,1) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabfo5ahnaaCaaaleqabaGaamyyaiaacI cacaaIXaGaaiykaaaakiabg2da9iaacIcacqGHRaWkcaaIXaGaaiil aiabgkHiTiaaigdacaGGSaGaey4kaSIaaGymaiaacYcacqGHsislca aIXaGaaiykaaaa@3FFA@

Linear brick

Γ a(1) =(+1,+1,1,1,1,1,+1,+1) Γ a(2) =(+1,1,1,+1,1,+1,+1,1) Γ a(3) =(+1,1,+1,1,+1,1,+1,1) Γ a(4) =(1,+1,1,+1,+1,1,+1,1) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOabaeqabaGaeu4KdC0aaWbaaSqabeaacaWGHb GaaiikaiaaigdacaGGPaaaaOGaeyypa0JaaiikaiabgUcaRiaaigda caGGSaGaey4kaSIaaGymaiaacYcacqGHsislcaaIXaGaaiilaiabgk HiTiaaigdacaGGSaGaeyOeI0IaaGymaiaacYcacqGHsislcaaIXaGa aiilaiabgUcaRiaaigdacaGGSaGaey4kaSIaaGymaiaacMcaaeaacq qHtoWrdaahaaWcbeqaaiaadggacaGGOaGaaGOmaiaacMcaaaGccqGH 9aqpcaGGOaGaey4kaSIaaGymaiaacYcacqGHsislcaaIXaGaaiilai abgkHiTiaaigdacaGGSaGaey4kaSIaaGymaiaacYcacqGHsislcaaI XaGaaiilaiabgUcaRiaaigdacaGGSaGaey4kaSIaaGymaiaacYcacq GHsislcaaIXaGaaiykaaqaaiabfo5ahnaaCaaaleqabaGaamyyaiaa cIcacaaIZaGaaiykaaaakiabg2da9iaacIcacqGHRaWkcaaIXaGaai ilaiabgkHiTiaaigdacaGGSaGaey4kaSIaaGymaiaacYcacqGHsisl caaIXaGaaiilaiabgUcaRiaaigdacaGGSaGaeyOeI0IaaGymaiaacY cacqGHRaWkcaaIXaGaaiilaiabgkHiTiaaigdacaGGPaaabaGaeu4K dC0aaWbaaSqabeaacaWGHbGaaiikaiaaisdacaGGPaaaaOGaeyypa0 JaaiikaiabgkHiTiaaigdacaGGSaGaey4kaSIaaGymaiaacYcacqGH sislcaaIXaGaaiilaiabgUcaRiaaigdacaGGSaGaey4kaSIaaGymai aacYcacqGHsislcaaIXaGaaiilaiabgUcaRiaaigdacaGGSaGaeyOe I0IaaGymaiaacMcaaaaa@93E7@

 

Reduced integration with hourglass control: Hourglassing in 4 noded quadrilateral and 8 noded brick elements can also be cured by adding an artificial stiffness to the element that acts to constrain the hourglass mode.   The stiffness must be carefully chosen so as to influence only the hourglass mode.  Only the final result will be given here MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  for details see D.P. Flanagan and T. Belytschko, International J. Num Meth in Engineering, 17, pp. 679 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@ 706, (1981).  To compute the corrective term:

1.      Define a series of `hourglass base vectors’  Γ a(i) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabfo5ahnaaCaaaleqabaGaamyyaiaacI cacaWGPbGaaiykaaaaaaa@352A@  which specify the displacements of the ath node in the ith hourglass mode.   The 4 noded quadrilateral element has only one hourglass mode; the 8 noded brick has 4 modes, listed in the table.

2.      Calculate the `hourglass shape vectors’ for each mode as follows

γ a(i) = Γ a(i) N a (ξ=0) x j b=1 n e Γ b(i) x j b MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeo7aNnaaCaaaleqabaGaamyyaiaacI cacaWGPbGaaiykaaaakiabg2da9iabfo5ahnaaCaaaleqabaGaamyy aiaacIcacaWGPbGaaiykaaaakiabgkHiTmaalaaabaGaeyOaIyRaam OtamaaCaaaleqabaGaamyyaaaakiaacIcacaWH+oGaeyypa0JaaGim aiaacMcaaeaacqGHciITcaWG4bWaaSbaaSqaaiaadQgaaeqaaaaakm aaqahabaGaeu4KdC0aaWbaaSqabeaacaWGIbGaaiikaiaadMgacaGG PaaaaOGaamiEamaaDaaaleaacaWGQbaabaGaamOyaaaaaeaacaWGIb Gaeyypa0JaaGymaaqaaiaad6gadaWgaaadbaGaamyzaaqabaaaniab ggHiLdaaaa@5639@

where n e MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaad6gadaWgaaWcbaGaamyzaaqabaaaaa@3271@  denotes the number of nodes on the element.

3.      The modified stiffness matrix for the element is written as

k aibk (l) = V e (l) C ijkl N a (x) x j N b (x) x l dV +κ V e (l) m γ a(m) γ b(m) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4AamaaDaaaleaacaWGHbGaamyAai aadkgacaWGRbaabaGaaiikaiaadYgacaGGPaaaaOGaeyypa0Zaa8qu aeaacaWGdbWaaSbaaSqaaiaadMgacaWGQbGaam4AaiaadYgaaeqaaO WaaSaaaeaacqGHciITcaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiik aiaahIhacaGGPaaabaGaeyOaIyRaamiEamaaBaaaleaacaWGQbaabe aaaaGcdaWcaaqaaiabgkGi2kaad6eadaahaaWcbeqaaiaadkgaaaGc caGGOaGaaCiEaiaacMcaaeaacqGHciITcaWG4bWaaSbaaSqaaiaadY gaaeqaaaaakiaadsgacaWGwbaaleaacaWGwbWaa0baaWqaaiaadwga aeaacaGGOaGaamiBaiaacMcaaaaaleqaniabgUIiYdGccaaMc8UaaG PaVlabgUcaRiabeQ7aRjaadAfadaqhaaWcbaGaamyzaaqaaiaacIca caWGSbGaaiykaaaakmaaqafabaGaeq4SdC2aaWbaaSqabeaacaWGHb Gaaiikaiaad2gacaGGPaaaaOGaeq4SdC2aaWbaaSqabeaacaWGIbGa aiikaiaad2gacaGGPaaaaaqaaiaad2gaaeqaniabggHiLdaaaa@6FCE@

where V e MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadAfadaWgaaWcbaGaamyzaaqabaaaaa@3259@  denotes the volume of the element, and κ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeQ7aRbaa@321A@  is a numerical parameter that controls the stiffness of the hourglass resistance.  Taking κ=0.01μ( N a / x j )( N a / x j ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeQ7aRjabg2da9iaaicdacaGGUaGaaG imaiaaigdacqaH8oqBdaqadaqaaiabgkGi2kaad6eadaahaaWcbeqa aiaadggaaaGccaGGVaGaeyOaIyRaamiEamaaBaaaleaacaWGQbaabe aaaOGaayjkaiaawMcaamaabmaabaGaeyOaIyRaamOtamaaCaaaleqa baGaamyyaaaakiaac+cacqGHciITcaWG4bWaaSbaaSqaaiaadQgaae qaaaGccaGLOaGaayzkaaaaaa@49EB@  where μ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeY7aTbaa@321E@  is the elastic shear modulus works well for most applications.  If κ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeQ7aRbaa@321A@  is too large, it will seriously over-stiffen the solid.

 

Sample code: Reduced integration with hourglass control has been implemented in the sample code Fem_hourglasscontrol.mws.   When run with the input file volumetric_locking_demo.txt it produces the results shown in the picture.  Hourglassing has clearly been satisfactorily eliminated.

 

HEALTH WARNING: Hourglass control is not completely effective: it can fail for finite strain problems and can also cause problems in a dynamic analysis, where the low stiffness of the hourglass modes can introduce spurious low frequency vibration modes and low wave speeds.

 

 

8.6.3 Hybrid elements for modeling near-incompressible materials

 

The bulk elastic modulus is infinite for a fully incompressible material, which leads to an infinite stiffness matrix in the standard finite element formulation (even if reduced integration is used to avoid locking).   This behavior can cause the stiffness matrix for a nearly incompressible material to become ill-conditioned, so that small rounding errors during the computation result in large errors in the solution.

 

Hybrid elements are designed to avoid this problem.  They work by including the hydrostatic stress distribution as an additional unknown variable, which must be computed at the same time as the displacement field.  This allows the stiff terms to be removed from the system of finite element equations.  The procedure is illustrated most easily using isotropic linear elasticity, but in practice hybrid elements are usually used to simulate rubbers or metals subjected to large plastic strains.

 

Hybrid elements are based on a modified version of the principle of virtual work, as follows. The virtual work equation (for small strains) is re-written as

R ( S ij δ ε ij )d V 0 + R pδ ε qq d V 0 + R δp( σ kk 3K p K )d V 0 R b i δ v i d V 0 2 R t i δ v i dA =0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaKaaGgaakmaabmaabaGaam4uam aaBaaaleaacaWGPbGaamOAaaqabaGccqaH0oazcqaH1oqzdaWgaaWc baGaamyAaiaadQgaaeqaaaGccaGLOaGaayzkaaqcaaQaaGPaVlaads gacaWGwbGcdaWgaaqcbawaaiaaicdaaeqaaaqaaiaadkfaaeqajmay cqGHRiI8aOGaey4kaSYaa8quaKaaGgaacaWGWbGaeqiTdqMaeqyTdu McdaWgaaqcbaAaaiaadghacaWGXbaabeaajaaOcaaMc8Uaamizaiaa dAfakmaaBaaajeaybaGaaGimaaqabaaabaGaamOuaaqabKWaGjabgU IiYdGccqGHRaWkdaWdrbqcaaAaaiabes7aKjaadchakmaabmaabaWa aSaaaeaacqaHdpWCdaWgaaWcbaGaam4AaiaadUgaaeqaaaGcbaGaaG 4maiaadUeaaaGaeyOeI0YaaSaaaeaacaWGWbaabaGaam4saaaaaiaa wIcacaGLPaaajaaOcaaMc8UaamizaiaadAfakmaaBaaajeaybaGaaG imaaqabaaabaGaamOuaaqabKWaGjabgUIiYdGccqGHsisldaWdrbqa aiaadkgadaWgaaWcbaGaamyAaaqabaGccqaH0oazcaWG2bWaaSbaaS qaaiaadMgaaeqaaOGaamizaiaadAfadaWgaaWcbaGaaGimaaqabaGc cqGHsisldaWdrbqaaiaadshadaWgaaWcbaGaamyAaaqabaGccqaH0o azcaWG2bWaaSbaaSqaaiaadMgaaeqaaOGaamizaiaadgeaaSqaaiab gkGi2oaaBaaameaacaaIYaaabeaaliaadkfaaeqaniabgUIiYdaale aacaWGsbaabeqdcqGHRiI8aOGaeyypa0JaaGimaaaa@8876@

Here

1.      S ij = σ ij σ kk δ ij /3 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadofadaWgaaWcbaGaamyAaiaadQgaae qaaOGaeyypa0Jaeq4Wdm3aaSbaaSqaaiaadMgacaWGQbaabeaakiab gkHiTiabeo8aZnaaBaaaleaacaWGRbGaam4AaaqabaGccqaH0oazda WgaaWcbaGaamyAaiaadQgaaeqaaOGaai4laiaaiodaaaa@421D@  is the deviatoric stress, determined from the displacement field;

2.      σ kk /3 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeo8aZnaaBaaaleaacaWGRbGaam4Aaa qabaGccaGGVaGaaG4maaaa@35B1@  is the hydrostatic stress, again determined from the displacement field;

3.      p is an additional degree of freedom that represents the (unknown) hydrostatic stress in the solid;

4.      δp MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabes7aKjaadchaaaa@3302@  is an arbitrary variation in the hydrostatic stress;

5.      K=(1/3) σ kk / ε pp MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadUeacqGH9aqpcaGGOaGaaGymaiaac+ cacaaIZaGaaiykaiabgkGi2kabeo8aZnaaBaaaleaacaWGRbGaam4A aaqabaGccaGGVaGaeyOaIyRaeqyTdu2aaSbaaSqaaiaadchacaWGWb aabeaaaaa@40D7@  is the bulk modulus of the solid.

The modified virtual work principle states that, if the virtual work equation is satisfied for all kinematically admissible variations in displacement δ v i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabes7aKjaadAhadaWgaaWcbaGaamyAaa qabaaaaa@3422@  and strain δ ε ij =( δ v i / x j +δ v j / x i ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabes7aKjabew7aLnaaBaaaleaacaWGPb GaamOAaaqabaGccqGH9aqpdaqadaqaaiabgkGi2kabes7aKjaadAha daWgaaWcbaGaamyAaaqabaGccaGGVaGaeyOaIyRaamiEamaaBaaale aacaWGQbaabeaakiabgUcaRiabgkGi2kabes7aKjaadAhadaWgaaWc baGaamOAaaqabaGccaGGVaGaeyOaIyRaamiEamaaBaaaleaacaWGPb aabeaaaOGaayjkaiaawMcaaaaa@4C02@  and all possible variations in pressure δp MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabes7aKjaadchaaaa@3302@ , the stress field will satisfy the equilibrium equations and traction boundary conditions, and the pressure variable p will be equal to the hydrostatic stress in the solid.

  

The finite element equations are derived in the usual way

1.      The displacement field, virtual displacement field and position in the each element are interpolated using the standard interpolation functions defined in Sections 8.1.9 and 8.1.10  as

u i = a=1 N e N a ( ξ j ) u i a δ v i = a=1 N e N a ( ξ j )δ v i a x i = a=1 N e N a ( ξ j ) x i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaBaaaleaacaWGPbaabeaaki abg2da9maaqahabaGaamOtamaaCaaaleqabaGaamyyaaaakiaacIca cqaH+oaEdaWgaaWcbaGaamOAaaqabaGccaGGPaGaamyDamaaDaaale aacaWGPbaabaGaamyyaaaaaeaacaWGHbGaeyypa0JaaGymaaqaaiaa d6eadaWgaaadbaGaamyzaaqabaaaniabggHiLdGccaaMc8UaaGPaVl aaykW7caaMc8UaeqiTdqMaamODamaaBaaaleaacaWGPbaabeaakiab g2da9maaqahabaGaamOtamaaCaaaleqabaGaamyyaaaakiaacIcacq aH+oaEdaWgaaWcbaGaamOAaaqabaGccaGGPaGaeqiTdqMaamODamaa DaaaleaacaWGPbaabaGaamyyaaaaaeaacaWGHbGaeyypa0JaaGymaa qaaiaad6eadaWgaaadbaGaamyzaaqabaaaniabggHiLdGccaaMc8Ua aGPaVlaaykW7caWG4bWaaSbaaSqaaiaadMgaaeqaaOGaeyypa0Zaaa bCaeaacaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiikaiabe67a4naa BaaaleaacaWGQbaabeaakiaacMcacaWG4bWaa0baaSqaaiaadMgaae aacaWGHbaaaaqaaiaadggacqGH9aqpcaaIXaaabaGaamOtamaaBaaa meaacaWGLbaabeaaa0GaeyyeIuoaaaa@786A@

2.      Since the pressure is now an independent variable, it must also be interpolated.  We write

p(x)= a=1 N p M a ( ξ j ) p a δp(x)= a=1 N p M a ( ξ j )δ p a x i = a=1 N e N a ( ξ j ) x i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiCaiaacIcacaWH4bGaaiykaiabg2 da9maaqahabaGaamytamaaCaaaleqabaGaamyyaaaakiaacIcacqaH +oaEdaWgaaWcbaGaamOAaaqabaGccaGGPaGaamiCamaaCaaaleqaba GaamyyaaaaaeaacaWGHbGaeyypa0JaaGymaaqaaiaad6eadaWgaaad baGaamiCaaqabaaaniabggHiLdGccaaMc8UaaGPaVlaaykW7caaMc8 UaeqiTdqMaamiCaiaacIcacaWH4bGaaiykaiabg2da9maaqahabaGa amytamaaCaaaleqabaGaamyyaaaakiaacIcacqaH+oaEdaWgaaWcba GaamOAaaqabaGccaGGPaGaeqiTdqMaamiCamaaCaaaleqabaGaamyy aaaaaeaacaWGHbGaeyypa0JaaGymaaqaaiaad6eadaWgaaadbaGaam iCaaqabaaaniabggHiLdGccaaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caWG4bWaaSbaaSqaaiaadMgaaeqaaOGaeyypa0ZaaabCae aacaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiikaiabe67a4naaBaaa leaacaWGQbaabeaakiaacMcacaWG4bWaa0baaSqaaiaadMgaaeaaca WGHbaaaaqaaiaadggacqGH9aqpcaaIXaaabaGaamOtamaaBaaameaa caWGLbaabeaaa0GaeyyeIuoaaaa@7D99@

where p a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadchadaahaaWcbeqaaiaadggaaaaaaa@3270@  are a discrete set of pressure variables, δ p a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabes7aKjaadchadaahaaWcbeqaaiaadg gaaaaaaa@3415@  is an arbitrary change in these pressure variables, M a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaad2eadaahaaWcbeqaaiaadggaaaaaaa@324D@  are a set of interpolation functions for the pressure, and N p MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaad6eadaWgaaWcbaGaamiCaaqabaaaaa@325C@  is the number of pressure variables associated with the element.  The pressure need not be continuous across neighboring elements, so that independent pressure variables can be added to each element.  The following schemes are usually used

a.       In linear elements (the 3 noded triangle, 4 noded quadrilateral, 5 noded tetrahedron or 8 noded brick) the pressure is constant.  The pressure is defined by its value at the centroid of each element, and the interpolation functions are constant.

b.      In quadratic elements (6 noded triangle, 8 noded quadrilateral, 10 noded tetrahedron or 20 noded brick) the pressure varies linearly in the element.  Its value can be defined by the pressure at the corners of each element, and interpolated using the standard linear interpolation functions.

3.      For an isotropic, linear elastic solid with shear modulus μ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeY7aTbaa@321E@  and Poisson ratio ν MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabe27aUbaa@3220@  the deviatoric stress is related to the displacement field by  S ij =μ( δ il δ jk + δ ik δ jl 2 δ ij δ kl /3)( u k / x l ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadofadaWgaaWcbaGaamyAaiaadQgaae qaaOGaeyypa0JaeqiVd0Maaiikaiabes7aKnaaBaaaleaacaWGPbGa amiBaaqabaGccqaH0oazdaWgaaWcbaGaamOAaiaadUgaaeqaaOGaey 4kaSIaeqiTdq2aaSbaaSqaaiaadMgacaWGRbaabeaakiabes7aKnaa BaaaleaacaWGQbGaamiBaaqabaGccqGHsislcaaIYaGaeqiTdq2aaS baaSqaaiaadMgacaWGQbaabeaakiabes7aKnaaBaaaleaacaWGRbGa amiBaaqabaGccaGGVaGaaG4maiaacMcadaqadaqaaiabgkGi2kaadw hadaWgaaWcbaGaam4AaaqabaGccaGGVaGaeyOaIyRaamiEamaaBaaa leaacaWGSbaabeaaaOGaayjkaiaawMcaaaaa@5B0B@ , while the hydrostatic stress is σ kk /3=K( u k / x k ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeo8aZnaaBaaaleaacaWGRbGaam4Aaa qabaGccaGGVaGaaG4maiabg2da9iaadUeacaGGOaGaeyOaIyRaamyD amaaBaaaleaacaWGRbaabeaakiaac+cacqGHciITcaWG4bWaaSbaaS qaaiaadUgaaeqaaOGaaiykaaaa@40A2@ , where K=E/3(12ν) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadUeacqGH9aqpcaWGfbGaai4laiaaio dacaGGOaGaaGymaiabgkHiTiaaikdacqaH9oGBcaGGPaaaaa@39ED@  is the bulk modulus.

4.      Substituting the linear elastic equations and the finite element interpolation functions into the virtual work principle leads to a system of equations for the unknown displacements u k b MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadwhadaqhaaWcbaGaam4Aaaqaaiaadk gaaaaaaa@3386@  and pressures p b MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadchadaahaaWcbeqaaiaadkgaaaaaaa@3291@  of the form

K aibk u k b + Q aib p b = F i a {a,i}: x k a not on  1 R Q akb u k b Π ab p b =0 u i a = u i * ( x i a ){a,i}: x k a on  1 R MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGlbWaaSbaaSqaaiaadggaca WGPbGaamOyaiaadUgaaeqaaOGaamyDamaaDaaaleaacaWGRbaabaGa amOyaaaakiabgUcaRiaadgfadaWgaaWcbaGaamyyaiaadMgacaWGIb aabeaakiaadchadaahaaWcbeqaaiaadkgaaaGccqGH9aqpcaWGgbWa a0baaSqaaiaadMgaaeaacaWGHbaaaOGaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlabgc GiIiaacUhacaWGHbGaaiilaiaadMgacaGG9bGaaGPaVlaaykW7caaM c8UaaiOoaiaaykW7caaMc8UaaGPaVlaadIhadaqhaaWcbaGaam4Aaa qaaiaadggaaaGccaaMc8UaaeOBaiaab+gacaqG0bGaaeiiaiaab+ga caqGUbGaaeiiaiabgkGi2oaaBaaaleaacaqGXaaabeaakiaadkfaae aacaWGrbWaaSbaaSqaaiaadggacaWGRbGaamOyaaqabaGccaWG1bWa a0baaSqaaiaadUgaaeaacaWGIbaaaOGaeyOeI0IaeuiOda1aaSbaaS qaaiaadggacaWGIbaabeaakiaadchadaahaaWcbeqaaiaadkgaaaGc cqGH9aqpcaaIWaaabaGaamyDamaaDaaaleaacaWGPbaabaGaamyyaa aakiabg2da9iaadwhadaqhaaWcbaGaamyAaaqaaiaacQcaaaGccaGG OaGaamiEamaaDaaaleaacaWGPbaabaGaamyyaaaakiaacMcacaaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGzaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7cqGHaiIicaGG7bGaamyyaiaacYcaca WGPbGaaiyFaiaaykW7caaMc8UaaGPaVlaacQdacaaMc8UaaGPaVlaa dIhadaqhaaWcbaGaam4AaaqaaiaadggaaaGccaaMc8UaaGPaVlaab+ gacaqGUbGaaeiiaiabgkGi2oaaBaaaleaacaqGXaaabeaakiaadkfa aaaa@F6DA@

where the global stiffness matrices K,Q,Π MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaahUeacaGGSaGaaCyuaiaacYcacaWHGo aaaa@34A2@  are obtained by summing the following element stiffness matrices

k aibk = V e C ijkl N a (x) x j N b (x) x l K N a (x) x i N b (x) x k dV q aib = V N a (x) x i M b (x)dV π ab = V e 1 K M a (x) M b (x)dV MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGceaqabeaacaWGRbWaaSbaaSqaaiaadggaca WGPbGaamOyaiaadUgaaeqaaOGaeyypa0Zaa8quaeaacaWGdbWaaSba aSqaaiaadMgacaWGQbGaam4AaiaadYgaaeqaaOWaaSaaaeaacqGHci ITcaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiikaiaahIhacaGGPaaa baGaeyOaIyRaamiEamaaBaaaleaacaWGQbaabeaaaaGcdaWcaaqaai abgkGi2kaad6eadaahaaWcbeqaaiaadkgaaaGccaGGOaGaaCiEaiaa cMcaaeaacqGHciITcaWG4bWaaSbaaSqaaiaadYgaaeqaaaaakiabgk HiTiaadUeadaWcaaqaaiabgkGi2kaad6eadaahaaWcbeqaaiaadgga aaGccaGGOaGaaCiEaiaacMcaaeaacqGHciITcaWG4bWaaSbaaSqaai aadMgaaeqaaaaakmaalaaabaGaeyOaIyRaamOtamaaCaaaleqabaGa amOyaaaakiaacIcacaWH4bGaaiykaaqaaiabgkGi2kaadIhadaWgaa WcbaGaam4AaaqabaaaaOGaamizaiaadAfaaSqaaiaadAfadaWgaaad baGaamyzaaqabaaaleqaniabgUIiYdGccaaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caWGXbWaaSbaaSqaaiaadggacaWGPbGaamOyaaqabaGccqGH9a qpdaWdrbqaamaalaaabaGaeyOaIyRaamOtamaaCaaaleqabaGaamyy aaaakiaacIcacaWH4bGaaiykaaqaaiabgkGi2kaadIhadaWgaaWcba GaamyAaaqabaaaaOGaamytamaaCaaaleqabaGaamOyaaaakiaacIca caWH4bGaaiykaiaadsgacaWGwbaaleaacaWGwbaabeqdcqGHRiI8aa GcbaGaeqiWda3aaSbaaSqaaiaadggacaWGIbaabeaakiabg2da9maa pefabaWaaSaaaeaacaaIXaaabaGaam4saaaacaWGnbWaaWbaaSqabe aacaWGHbaaaOGaaiikaiaahIhacaGGPaGaamytamaaCaaaleqabaGa amOyaaaakiaacIcacaWH4bGaaiykaiaadsgacaWGwbaaleaacaWGwb WaaSbaaWqaaiaadwgaaeqaaaWcbeqdcqGHRiI8aaaaaa@A772@

and the force F is defined in the usual way.  The integrals defining k iakb MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadUgadaWgaaWcbaGaamyAaiaadggaca WGRbGaamOyaaqabaaaaa@354F@  may be evaluated using the full integration scheme or reduced integration (hourglass control may be required in this case).  The remaining integrals must be evaluated using reduced integration to avoid element locking.

5.      Note that, although the pressure variables are local to the elements, they cannot be eliminated from the element stiffness matrix, since doing so would reduce the element stiffness matrix to the usual, non-hybrid form.  Consequently, hybrid elements increase the cost of storing and solving the system of equations.