Chapter 8

 

Theory and Implementation of the Finite Element Method

 

 

 

The derivation and implementation of the finite element method outlined in the previous chapter is simple and easy to follow, but it gives the misleading impression that the finite element method relies on the principle of minimum potential energy, and so is applicable only to linear elastic solids.  This is not the case, of course MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  the finite element method can solve problems involving very complex materials and large shape changes.

 

This chapter contains

  1. A more rigorous derivation of the finite element equations, based on the principle of virtual work, which was derived in Section 2.4.
  2. A more sophisticated implementation of finite element method for static linear elasticity. This includes more accurate element interpolation schemes, and also extends the finite element method to three dimensions.
  3. A discussion of time integration schemes that are used in finite element simulations of dynamic problems, and a discussion of modal techniques for dynamic linear elasticity problems;
  4. An extension of the finite element method to nonlinear materials, using the hypoelastic material model described in Section 3.2 as a representative nonlinear material;
  5. An extension of the finite element method to account for large shape changes, using finite strain elasticity as a representative example;
  6. A discussion of finite element procedures for history dependent solids, using small strain viscoplasticity as a representative example.
  7. A discussion of the phenomenon of `locking’ that can cause the standard finite element method to fail in certain circumstances.   Several techniques for avoiding locking are presented.

 

In addition, a set of sample finite element codes (implemented in MAPLE and MATLAB) are provided to illustrate the how the various finite element procedures are implemented in practice.

 

 

8.1 Generalized FEM for static linear elasticity

 

This section gives a more general derivation and implementation of the finite element method for static linear elastic solids than the energy-based derivation given in Chapter 7.

 

8.1.1 Review of the principle of virtual work

 

Governing equations: We begin by summarizing the usual governing equations of linear elasticity, which must be solved by the FEA code.  

 

Given:

1.      The shape of the solid in its unloaded condition R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOuaaaa@31B4@

2.      The initial stress field in the solid (we will take this to be zero in setting up our FEM code)

3.      The elastic constants for the solid C ijkl MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaam4qamaaBaaaleaacaWGPbGaamOAai aadUgacaWGSbaabeaaaaa@3581@

4.      The thermal expansion coefficients for the solid, and temperature distribution (we will take this to be zero for our FEM code, for simplicity)

5.      A body force distribution b MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCOyaaaa@31C8@  acting on the solid (Note that in this section we will use b to denote force per unit volume rather than force per unit mass, to avoid having to write out the mass density all the time)

6.      Boundary conditions, specifying displacements u * (x) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCyDamaaCaaaleqabaGaaiOkaaaaki aacIcacaWH4bGaaiykaaaa@350A@  on a portion 1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIy7aaSbaaSqaaiaaigdaaeqaaO GaamOuaaaa@33FD@  or tractions on a portion 2 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIy7aaSbaaSqaaiaaikdaaeqaaO GaamOuaaaa@33FE@  of the boundary of R

 

Calculate displacements, strains and stresses u i , ε ij , σ ij MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadwhadaWgaaWcbaGaamyAaaqabaGcca GGSaGaeqyTdu2aaSbaaSqaaiaadMgacaWGQbaabeaakiaacYcacqaH dpWCdaWgaaWcbaGaamyAaiaadQgaaeqaaaaa@3B8C@  satisfying the governing equations of static linear elasticity

1.      The strain-displacement equation ε ij = 1 2 ( u i / x j + u j / x i ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqyTdu2aaSbaaSqaaiaadMgacaWGQb aabeaakiabg2da9maalaaabaGaaGymaaqaaiaaikdaaaGaaiikaiab gkGi2kaadwhadaWgaaWcbaGaamyAaaqabaGccaGGVaGaeyOaIyRaam iEamaaBaaaleaacaWGQbaabeaakiabgUcaRiabgkGi2kaadwhadaWg aaWcbaGaamOAaaqabaGccaGGVaGaeyOaIyRaamiEamaaBaaaleaaca WGPbaabeaakiaacMcaaaa@48CF@

2.      The elastic stress-strain law σ ij = C ijkl ε kl MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4Wdm3aaSbaaSqaaiaadMgacaWGQb aabeaakiabg2da9iaadoeadaWgaaWcbaGaamyAaiaadQgacaWGRbGa amiBaaqabaGccqaH1oqzdaWgaaWcbaGaam4AaiaadYgaaeqaaaaa@3E1B@

3.      The equation of static equilibrium for stresses σ ij / x i + b j =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIyRaeq4Wdm3aaSbaaSqaaiaadM gacaWGQbaabeaakiaac+cacqGHciITcaWG4bWaaSbaaSqaaiaadMga aeqaaOGaey4kaSIaamOyamaaBaaaleaacaWGQbaabeaakiabg2da9i aaicdaaaa@3EF3@

4.      The boundary conditions on displacement and stress u i = u i * on 1 R σ ij n i = t j * on 2 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaBaaaleaacaWGPbaabeaaki abg2da9iaadwhadaqhaaWcbaGaamyAaaqaaiaacQcaaaGccaaMc8Ua aGPaVlaab+gacaqGUbGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7cq GHciITdaWgaaWcbaGaaGymaaqabaGccaWGsbGaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8Uaeq4Wdm3aaSbaaSqaaiaadMgacaWGQbaabeaakiaad6ga daWgaaWcbaGaamyAaaqabaGccqGH9aqpcaWG0bWaa0baaSqaaiaadQ gaaeaacaGGQaaaaOGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaab+gacaqGUbGaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7cqGHciITdaWgaaWcbaGaaGOmaaqabaGccaWGsbaaaa@7A78@

 

The principle of virtual work: As we discussed in section 2.4, the principle of virtual work can be used to replace the stress equilibrium equations.  

 

To express the principle, we define a kinematically admissible virtual displacement field δv(x) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaaCODaiaacIcacaWH4bGaai ykaaaa@35CB@ , satisfying  δv=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaaCODaiabg2da9iaaicdaaa a@3531@  on 1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIy7aaSbaaSqaaiaaigdaaeqaaO GaamOuaaaa@33FD@ .  You can visualize this field as a small change in the displacement of the solid, if you like, but it is really just an arbitrary differentiable vector field.  The term `kinematically admissible’ is just a complicated way of saying that δv=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaaCODaiabg2da9iaaicdaaa a@3531@  on 1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIy7aaSbaaSqaaiaaigdaaeqaaO GaamOuaaaa@33FD@  - that is to say, if you perturb the displacement slightly, the boundary conditions on displacement are still satisfied.

 

In addition, we define an associated virtual strain field

δ ε ij = 1 2 ( δ v i x j + δ v j x i ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rk0le9 v8qqaqFD0xXdHaVhbbf9v8qqaqFr0xc9pk0xbba9q8WqFfea0=yr0R Yxir=Jbba9q8aq0=yq=He9q8qqQ8frFve9Fve9Ff0dmeaabaqaciGa caGaaeqabaGadeaadaaakeaacqaH0oazcqaH1oqzdaWgaaWcbaGaam yAaiaadQgaaeqaaOGaeyypa0ZaaSaaaeaacaaIXaaabaGaaGOmaaaa daqadaqaamaalaaabaGaeyOaIyRaeqiTdqMaamODamaaBaaaleaaca WGPbaabeaaaOqaaiabgkGi2kaadIhadaWgaaWcbaGaamOAaaqabaaa aOGaey4kaSYaaSaaaeaacqGHciITcqaH0oazcaWG2bWaaSbaaSqaai aadQgaaeqaaaGcbaGaeyOaIyRaamiEamaaBaaaleaacaWGPbaabeaa aaaakiaawIcacaGLPaaaaaa@4EA1@

The principle of virtual work states that if the stress field σ ij MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeo8aZnaaBaaaleaacaWGPbGaamOAaa qabaaaaa@3434@  satisfies

R σ ij δ ε ij d V 0 R b i δ v i d V 0 2 R t i δ v i dA=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaKaaGgaacqaHdpWCkmaaBaaale aacaWGPbGaamOAaaqabaGccqaH0oazcqaH1oqzdaWgaaWcbaGaamyA aiaadQgaaeqaaKaaGkaaykW7caWGKbGaamOvaOWaaSbaaKqaGfaaca aIWaaabeaajaaOcqGHsislkmaapefajaaObaGaamOyaOWaaSbaaKqa GfaacaWGPbaabeaakiabes7aKjaadAhadaWgaaWcbaGaamyAaaqaba qcaaQaamizaiaadAfakmaaBaaajeaybaGaaGimaaqabaqcaaQaeyOe I0IcdaWdrbqcaaAaaiaadshakmaaBaaajeaybaGaamyAaaqabaGccq aH0oazcaWG2bWaaSbaaSqaaiaadMgaaeqaaaqcbaAaaiabgkGi2UWa aSbaaKGaGgaacaaIYaaabeaajeaycaWGsbaajeaObeqcdaQaey4kIi paaKqaGgaacaWGsbaabeqcdaQaey4kIipaaKqaGfaacaWGsbaabeqc daMaey4kIipakiaadsgacaWGbbqcaaMaeyypa0JaaGimaaaa@693C@

for all possible virtual displacement fields and corresponding virtual strains, it will automatically satisfy the equation of stress equilibrium σ ij / x i + b j =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIyRaeq4Wdm3aaSbaaSqaaiaadM gacaWGQbaabeaakiaac+cacqGHciITcaWG4bWaaSbaaSqaaiaadMga aeqaaOGaey4kaSIaamOyamaaBaaaleaacaWGQbaabeaakiabg2da9i aaicdaaaa@3EF3@ , and also the traction boundary condition σ ij n i = t j * on 2 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4Wdm3aaSbaaSqaaiaadMgacaWGQb aabeaakiaad6gadaWgaaWcbaGaamyAaaqabaGccqGH9aqpcaWG0bWa a0baaSqaaiaadQgaaeaacaGGQaaaaOGaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaab+gacaqGUbGaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7cqGHciITdaWgaaWcbaGaaGOmaaqabaGccaWGsb aaaa@5225@

 

 

8.1.2 Integral (weak) form of the governing equations of linear elasticity

 

The principle of virtual work can be used to write the governing equation for the displacement field in a linear elastic solid in an integral form (called the `weak form’).  Instead of solving the governing equations listed in the preceding section, the displacements, strains and stresses u i , ε ij , σ ij MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadwhadaWgaaWcbaGaamyAaaqabaGcca GGSaGaeqyTdu2aaSbaaSqaaiaadMgacaWGQbaabeaakiaacYcacqaH dpWCdaWgaaWcbaGaamyAaiaadQgaaeqaaaaa@3B8C@  are calculated as follows. 

1.      Find a displacement field u i ( x j ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadwhadaWgaaWcbaGaamyAaaqabaGcca GGOaGaamiEamaaBaaaleaacaWGQbaabeaakiaacMcaaaa@3601@  satisfying

R C ijkl u k x l δ v i x j dV R b i δ v i dV 2 R t i * δ v i dA=0, u i = u i * on 1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaeaacaWGdbWaaSbaaSqaaiaadM gacaWGQbGaam4AaiaadYgaaeqaaOWaaSaaaeaacqGHciITcaWG1bWa aSbaaSqaaiaadUgaaeqaaaGcbaGaeyOaIyRaamiEamaaBaaaleaaca WGSbaabeaaaaGcdaWcaaqaaiabgkGi2kabes7aKjaadAhadaWgaaWc baGaamyAaaqabaaakeaacqGHciITcaWG4bWaaSbaaSqaaiaadQgaae qaaaaakiaadsgacaWGwbGaeyOeI0Yaa8quaeaacaWGIbWaaSbaaSqa aiaadMgaaeqaaOGaeqiTdqMaamODamaaBaaaleaacaWGPbaabeaaki aadsgacaWGwbGaeyOeI0Yaa8quaeaacaWG0bWaa0baaSqaaiaadMga aeaacaGGQaaaaOGaeqiTdqMaamODamaaBaaaleaacaWGPbaabeaaki aadsgacaWGbbGaeyypa0JaaGimaiaacYcacaaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8oaleaacqGHciITdaWgaaadbaGaaG OmaaqabaWccaWGsbaabeqdcqGHRiI8aaWcbaGaamOuaaqab0Gaey4k IipaaSqaaiaadkfaaeqaniabgUIiYdGccaaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa dwhadaWgaaWcbaGaamyAaaqabaGccqGH9aqpcaWG1bWaa0baaSqaai aadMgaaeaacaGGQaaaaOGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caqGVbGaaeOBaiaaykW7caaMc8UaeyOaIy7aaSbaaSqaaiaaigdaae qaaOGaamOuaaaa@98A9@

for all virtual velocity fields δ v i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaamODamaaBaaaleaacaWGPb aabeaaaaa@3489@  satisfying δ v i =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaamODamaaBaaaleaacaWGPb aabeaakiabg2da9iaaicdaaaa@3653@  on 1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIy7aaSbaaSqaaiaaigdaaeqaaO GaamOuaaaa@33FD@ .

2.      Compute the strains from the definition ε ij = 1 2 ( u i / x j + u j / x i ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqyTdu2aaSbaaSqaaiaadMgacaWGQb aabeaakiabg2da9maalaaabaGaaGymaaqaaiaaikdaaaGaaiikaiab gkGi2kaadwhadaWgaaWcbaGaamyAaaqabaGccaGGVaGaeyOaIyRaam iEamaaBaaaleaacaWGQbaabeaakiabgUcaRiabgkGi2kaadwhadaWg aaWcbaGaamOAaaqabaGccaGGVaGaeyOaIyRaamiEamaaBaaaleaaca WGPbaabeaakiaacMcaaaa@48CF@

3.      Compute the stresses from the stress-strain law σ ij = C ijkl ε kl MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4Wdm3aaSbaaSqaaiaadMgacaWGQb aabeaakiabg2da9iaadoeadaWgaaWcbaGaamyAaiaadQgacaWGRbGa amiBaaqabaGccqaH1oqzdaWgaaWcbaGaam4AaiaadYgaaeqaaaaa@3E1B@

The stress will automatically satisfy the equilibrium equation and boundary conditions, so all the field equations and boundary conditions will be satisfied.

 

The significance of this result is that it replaces the derivatives in the partial differential equations of equilibrium with an equivalent integral, which is easier to handle numerically.  It is essentially equivalent to replacing the equilibrium equation with the principle of minimum potential energy, but the procedure based on the principle of virtual work is very easily extended to dynamic problems, other stress-strain laws, and even to problems involving large shape changes.

 

Derivation: start with the virtual work equation

R σ ij δ ε ij d V 0 R b i δ v i d V 0 2 R t i δ v i dA=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaKaaahaacqaHdpWCkmaaBaaale aacaWGPbGaamOAaaqabaGccqaH0oazcqaH1oqzdaWgaaWcbaGaamyA aiaadQgaaeqaaKaaalaaykW7caWGKbGaamOvaOWaaSbaaSqaaiaaic daaeqaaKaaalabgkHiTOWaa8quaKaaahaacaWGIbGcdaWgaaWcbaGa amyAaaqabaqcaaSaeqiTdqMaamODaOWaaSbaaSqaaiaadMgaaeqaaK aaalaadsgacaWGwbGcdaWgaaWcbaGaaGimaaqabaqcaaSaeyOeI0Ic daWdrbqcaaCaaiaadshakmaaBaaaleaacaWGPbaabeaajaaWcqaH0o azcaWG2bGcdaWgaaWcbaGaamyAaaqabaaajeaWbaGaeyOaIy7cdaWg aaqccaCaaiaaikdaaeqaaSGaamOuaaqcbaCabKWaalabgUIiYdaaje aWbaGaamOuaaqabKWaalabgUIiYdaaleaacaWGsbaabeqdcqGHRiI8 aOGaamizaiaadgeacqGH9aqpcaaIWaaaaa@696C@

Recall that σ ij = σ ji MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeo8aZnaaBaaaleaacaWGPbGaamOAaa qabaGccqGH9aqpcqaHdpWCdaWgaaWcbaGaamOAaiaadMgaaeqaaaaa @3910@ , and that δ ε ij =(δ v i / x j +δ v j / x i )/2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaeqyTdu2aaSbaaSqaaiaadM gacaWGQbaabeaakiabg2da9iaacIcacqGHciITcqaH0oazcaWG2bWa aSbaaSqaaiaadMgaaeqaaOGaai4laiabgkGi2kaadIhadaWgaaWcba GaamOAaaqabaGccqGHRaWkcqGHciITcqaH0oazcaWG2bWaaSbaaSqa aiaadQgaaeqaaOGaai4laiabgkGi2kaadIhadaWgaaWcbaGaamyAaa qabaGccaGGPaGaai4laiaaikdaaaa@4DA8@ , so that

σ ij δ ε ij =( σ ij δ v i / x j + σ ji δ v j / x i )/2= σ ij δ v i / x j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4Wdm3aaSbaaSqaaiaadMgacaWGQb aabeaakiabes7aKjabew7aLnaaBaaaleaacaWGPbGaamOAaaqabaGc cqGH9aqpcaGGOaGaeq4Wdm3aaSbaaSqaaiaadMgacaWGQbaabeaaki abgkGi2kabes7aKjaadAhadaWgaaWcbaGaamyAaaqabaGccaGGVaGa eyOaIyRaamiEamaaBaaaleaacaWGQbaabeaakiabgUcaRiabeo8aZn aaBaaaleaacaWGQbGaamyAaaqabaGccqGHciITcqaH0oazcaWG2bWa aSbaaSqaaiaadQgaaeqaaOGaai4laiabgkGi2kaadIhadaWgaaWcba GaamyAaaqabaGccaGGPaGaai4laiaaikdacqGH9aqpcqaHdpWCdaWg aaWcbaGaamyAaiaadQgaaeqaaOGaeyOaIyRaeqiTdqMaamODamaaBa aaleaacaWGPbaabeaakiaac+cacqGHciITcaWG4bWaaSbaaSqaaiaa dQgaaeqaaaaa@6761@

Finally, recall that σ ij = C ijkl ε kl = C ijkl ( u k / x l + u l / x k ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeo8aZnaaBaaaleaacaWGPbGaamOAaa qabaGccqGH9aqpcaWGdbWaaSbaaSqaaiaadMgacaWGQbGaam4Aaiaa dYgaaeqaaOGaeqyTdu2aaSbaaSqaaiaadUgacaWGSbaabeaakiabg2 da9iaadoeadaWgaaWcbaGaamyAaiaadQgacaWGRbGaamiBaaqabaGc caGGOaGaeyOaIyRaamyDamaaBaaaleaacaWGRbaabeaakiaac+cacq GHciITcaWG4bWaaSbaaSqaaiaadYgaaeqaaOGaey4kaSIaeyOaIyRa amyDamaaBaaaleaacaWGSbaabeaakiaac+cacqGHciITcaWG4bWaaS baaSqaaiaadUgaaeqaaOGaaiykaaaa@5541@  and that the elastic compliances must satisfy C ijkl = C ijlk MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadoeadaWgaaWcbaGaamyAaiaadQgaca WGRbGaamiBaaqabaGccqGH9aqpcaWGdbWaaSbaaSqaaiaadMgacaWG QbGaamiBaiaadUgaaeqaaaaa@3ADC@  so that σ ij = C ijkl ( u k / x l ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabeo8aZnaaBaaaleaacaWGPbGaamOAaa qabaGccqGH9aqpcaWGdbWaaSbaaSqaaiaadMgacaWGQbGaam4Aaiaa dYgaaeqaaOGaaiikaiabgkGi2kaadwhadaWgaaWcbaGaam4Aaaqaba GccaGGVaGaeyOaIyRaamiEamaaBaaaleaacaWGSbaabeaakiaacMca aaa@431C@ .  Finally, this shows that

σ ij δ ε ij = C ijkl u i x j δ v i x j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4Wdm3aaSbaaSqaaiaadMgacaWGQb aabeaakiabes7aKjabew7aLnaaBaaaleaacaWGPbGaamOAaaqabaGc cqGH9aqpcaWGdbWaaSbaaSqaaiaadMgacaWGQbGaam4AaiaadYgaae qaaOWaaSaaaeaacqGHciITcaWG1bWaaSbaaSqaaiaadMgaaeqaaaGc baGaeyOaIyRaamiEamaaBaaaleaacaWGQbaabeaaaaGcdaWcaaqaai abgkGi2kabes7aKjaadAhadaWgaaWcbaGaamyAaaqabaaakeaacqGH ciITcaWG4bWaaSbaaSqaaiaadQgaaeqaaaaaaaa@4F9A@

Substituting into the virtual work equation gives the result we need.

 

 

 

8.1.3 Interpolating the displacement field and the virtual velocity field

 

To solve the integral form of the elasticity equations given in 8.1.2, we discretize the displacement field.  That is to say, we choose to calculate the displacement field at a set of n discrete points in the solid (called `nodes’ in finite element terminology).  We will denote the coordinates of these special points by x i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiEamaaDaaaleaacaWGPbaabaGaam yyaaaaaaa@33CD@ , where the superscript a ranges from 1 to N.  The unknown displacement vector at each nodal point will be denoted by u i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaDaaaleaacaWGPbaabaGaam yyaaaaaaa@33CA@ .

 

The displacement field at an arbitrary point within the solid will be specified by interpolating between nodal values in some convenient way.  An efficient and robust implementation of the finite element method requires a careful choice of interpolation scheme, but for now we will denote the interpolation in a general way as

u i (x)= a=1 n N a (x) u i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaBaaaleaacaWGPbaabeaaki aacIcacaWH4bGaaiykaiabg2da9maaqahabaGaamOtamaaCaaaleqa baGaamyyaaaakiaacIcacaWH4bGaaiykaiaadwhadaqhaaWcbaGaam yAaaqaaiaadggaaaaabaGaamyyaiabg2da9iaaigdaaeaacaWGUbaa niabggHiLdaaaa@4363@

Here, x denotes the coordinates of an arbitrary point in the solid.  The interpolation functions N a (x) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOtamaaCaaaleqabaGaamyyaaaaki aacIcacaWH4bGaaiykaaaa@3519@  are functions of position only, which must have the property that

u i b = a=1 n N a ( x b ) u i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaDaaaleaacaWGPbaabaGaam Oyaaaakiabg2da9maaqahabaGaamOtamaaCaaaleqabaGaamyyaaaa kiaacIcacaWH4bWaaWbaaSqabeaacaWGIbaaaOGaaiykaiaadwhada qhaaWcbaGaamyAaaqaaiaadggaaaaabaGaamyyaiabg2da9iaaigda aeaacaWGUbaaniabggHiLdaaaa@430F@

for all b=1…N(This is to make sure that the displacement field has the correct value at each node).  Recently developed meshless finite element methods use very complex interpolation functions, but the more traditional approach is to choose them so that

N a ( x b )={ 1a=b 0ab MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOtamaaCaaaleqabaGaamyyaaaaki aacIcacaWH4bWaaWbaaSqabeaacaWGIbaaaOGaaiykaiabg2da9maa ceaabaqbaeqabiqaaaqaaiaaigdacaaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aadggacqGH9aqpcaWGIbaabaGaaGimaiaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaamyyaiabgcMi5kaadkgaaaaacaGL7baaaaa@7DF8@

The simple constant strain triangle elements introduced in 7.1 are one example of this type of interpolation scheme.  We will define more complicated interpolation functions shortly.

 

We can obviously interpolate the virtual velocity field in exactly the same way (since the principle of virtual work must be satisfied for all virtual velocities, it must certainly be satisfied for an interpolated velocity field…) so that

δ v i (x)= a=1 n N a (x)δ v i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaamODamaaBaaaleaacaWGPb aabeaakiaacIcacaWH4bGaaiykaiabg2da9maaqahabaGaamOtamaa CaaaleqabaGaamyyaaaakiaacIcacaWH4bGaaiykaiabes7aKjaadA hadaqhaaWcbaGaamyAaaqaaiaadggaaaaabaGaamyyaiabg2da9iaa igdaaeaacaWGUbaaniabggHiLdaaaa@46AF@

where δ v i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8skY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaamODamaaDaaaleaacaWGPb aabaGaamyyaaaaaaa@3580@  are arbitrary nodal values of virtual velocity.

 

 

8.1.4 Finite element equations

 

Substituting the interpolated fields into the virtual work equation, we find that

R C ijkl N b (x) x l u k b N a (x) x j δ v i a dV R b i N a (x)δ v i a dV 2R t i * N a (x)δ v i a dA=0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaWaa8quaeaacaWGdbWaaSbaaSqaaiaadM gacaWGQbGaam4AaiaadYgaaeqaaOWaaSaaaeaacqGHciITcaWGobWa aWbaaSqabeaacaWGIbaaaOGaaiikaiaahIhacaGGPaaabaGaeyOaIy RaamiEamaaBaaaleaacaWGSbaabeaaaaGccaWG1bWaa0baaSqaaiaa dUgaaeaacaWGIbaaaOWaaSaaaeaacqGHciITcaWGobWaaWbaaSqabe aacaWGHbaaaOGaaiikaiaahIhacaGGPaaabaGaeyOaIyRaamiEamaa BaaaleaacaWGQbaabeaaaaGccqaH0oazcaWG2bWaa0baaSqaaiaadM gaaeaacaWGHbaaaOGaaGPaVlaadsgacaWGwbaaleaacaWGsbaabeqd cqGHRiI8aOGaeyOeI0Yaa8quaeaacaWGIbWaaSbaaSqaaiaadMgaae qaaOGaamOtamaaCaaaleqabaGaamyyaaaakiaacIcacaWH4bGaaiyk aiabes7aKjaadAhadaqhaaWcbaGaamyAaaqaaiaadggaaaaabaGaam Ouaaqab0Gaey4kIipakiaadsgacaWGwbGaeyOeI0Yaa8quaeaacaWG 0bWaa0baaSqaaiaadMgaaeaacaGGQaaaaOGaamOtamaaCaaaleqaba GaamyyaaaakiaacIcacaWH4bGaaiykaiabes7aKjaadAhadaqhaaWc baGaamyAaaqaaiaadggaaaaabaGaeyOaIyRaaGOmaiaadkfaaeqani abgUIiYdGccaWGKbGaamyqaiaaykW7cqGH9aqpcaaIWaaaaa@7DF1@

where summation on a and b is implied, in addition to the usual summation on i,j,k,l

 

Note that the interpolation functions are known functions of position. We can therefore re-write the virtual work equation in matrix form as

( K aibk u k b F i a )δ v i a =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaeWaaeaacaWGlbWaaSbaaSqaaiaadg gacaWGPbGaamOyaiaadUgaaeqaaOGaamyDamaaDaaaleaacaWGRbaa baGaamOyaaaakiabgkHiTiaadAeadaqhaaWcbaGaamyAaaqaaiaadg gaaaaakiaawIcacaGLPaaacqaH0oazcaWG2bWaa0baaSqaaiaadMga aeaacaWGHbaaaOGaeyypa0JaaGimaaaa@443F@

where

K aibk = C ijkl N a (x) x j N b (x) x l dV F i a = b i N a (x) dV+ 2 t i * N a (x) dA MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamyAai aadkgacaWGRbaabeaakiabg2da9maapefabaGaam4qamaaBaaaleaa caWGPbGaamOAaiaadUgacaWGSbaabeaakmaalaaabaGaeyOaIyRaam OtamaaCaaaleqabaGaamyyaaaakiaacIcacaWH4bGaaiykaaqaaiab gkGi2kaadIhadaWgaaWcbaGaamOAaaqabaaaaOWaaSaaaeaacqGHci ITcaWGobWaaWbaaSqabeaacaWGIbaaaOGaaiikaiaahIhacaGGPaaa baGaeyOaIyRaamiEamaaBaaaleaacaWGSbaabeaaaaGccaWGKbGaam OvaaWcbaGaeyihHimabeqdcqGHRiI8aOGaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaadAeadaqhaaWcbaGaamyAaaqaaiaadggaaaGccqGH 9aqpdaWdrbqaaiaadkgadaWgaaWcbaGaamyAaaqabaGccaWGobWaaW baaSqabeaacaWGHbaaaOGaaiikaiaahIhacaGGPaaaleaacqGHCeIW aeqaniabgUIiYdGccaWGKbGaamOvaiabgUcaRmaapefabaGaamiDam aaDaaaleaacaWGPbaabaGaaiOkaaaakiaad6eadaahaaWcbeqaaiaa dggaaaGccaGGOaGaaCiEaiaacMcaaSqaaiabgkGi2kaaikdacqGHCe IWaeqaniabgUIiYdGccaWGKbGaamyqaiaaykW7aaa@888D@

Here K is known as the `stiffness matrix’ and f is known as the force vector.  K is a function only of the elastic properties of the solid, its geometry, and the interpolation functions and nodal positions. It is therefore a known matrix.  Similarly, f is a function only of the known boundary loading and body force field, and the interpolation scheme and nodal positions.  Observe that the symmetry of the elasticity tensor implies that K also has some symmetry MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  specifically K aibk = K bkai MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamyAai aadkgacaWGRbaabeaakiabg2da9iaadUeadaWgaaWcbaGaamOyaiaa dUgacaWGHbGaamyAaaqabaaaaa@3B2D@ .

 

The virtual work equation must be satisfied for all possible sets of δ v i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaamODamaaDaaaleaacaWGPb aabaGaamyyaaaaaaa@3570@  with δ v i a =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaamODamaaDaaaleaacaWGPb aabaGaamyyaaaakiabg2da9iaaicdaaaa@3738@  for nodes a that lie on 1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIy7aaSbaaSqaaiaaigdaaeqaaO GaamOuaaaa@340B@ .   At these nodes, the displacements must satisfy u i a = u i * MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaDaaaleaacaWGPbaabaGaam yyaaaakiabg2da9iaadwhadaqhaaWcbaGaamyAaaqaaiaacQcaaaaa aa@379D@ .   Evidently, this requires

K aibk u k b = F i a {a,i}: x k a not on  1 R u i a = u i * ( x i a ){a,i}: x k a on  1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGlbWaaSbaaSqaaiaadggaca WGPbGaamOyaiaadUgaaeqaaOGaamyDamaaDaaaleaacaWGRbaabaGa amOyaaaakiabg2da9iaadAeadaqhaaWcbaGaamyAaaqaaiaadggaaa GccaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaeyiaIiIaai4EaiaadggacaGGSaGaamyA aiaac2hacaaMc8UaaGPaVlaaykW7caGG6aGaaGPaVlaaykW7caaMc8 UaamiEamaaDaaaleaacaWGRbaabaGaamyyaaaakiaaykW7caqGUbGa ae4BaiaabshacaqGGaGaae4Baiaab6gacaqGGaGaeyOaIy7aaSbaaS qaaiaabgdaaeqaaOGaamOuaaqaaiaadwhadaqhaaWcbaGaamyAaaqa aiaadggaaaGccqGH9aqpcaWG1bWaa0baaSqaaiaadMgaaeaacaGGQa aaaOGaaiikaiaadIhadaqhaaWcbaGaamyAaaqaaiaadggaaaGccaGG PaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaeyiaIiIaai4Eaiaadg gacaGGSaGaamyAaiaac2hacaaMc8UaaGPaVlaaykW7caGG6aGaaGPa VlaaykW7caWG4bWaa0baaSqaaiaadUgaaeaacaWGHbaaaOGaaGPaVl aaykW7caqGVbGaaeOBaiaabccacqGHciITdaWgaaWcbaGaaeymaaqa baGccaWGsbaaaaa@BC05@

This is a system of n linear equations for the n nodal displacements.

 

 

8.1.5 Simple 1D Implementation of the finite element method

 

Before describing a fully general 3D implementation of the finite element method, we will illustrate all the key ideas using a simple 1-D example. Consider a long linear elastic bar, as shown in the picture.  Assume

1.      The bar has shear modulus μ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiVd0gaaa@3285@  and Poisson’s ratio ν MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqyVd4gaaa@3287@

2.      The bar has cross section h×h MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiAaiabgEna0kaadIgaaaa@34C0@  and length L

3.      It is constrained on all its sides so that u 2 = u 3 =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaBaaaleaacaaIYaaabeaaki abg2da9iaadwhadaWgaaWcbaGaaG4maaqabaGccqGH9aqpcaaIWaaa aa@376E@

4.      The bar is subjected to body force b=b( x 1 ) e 1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCOyaiabg2da9iaadkgacaGGOaGaam iEamaaBaaaleaacaaIXaaabeaakiaacMcacaWHLbWaaSbaaSqaaiaa igdaaeqaaaaa@38C3@ ,

5.      The bar is either loaded or constrained at its ends, so that the boundary conditions are either  t 1 (0)= t * (0), t 1 (L)= t * (L) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiDamaaBaaaleaacaaIXaaabeaaki aacIcacaaIWaGaaiykaiabg2da9iaadshadaahaaWcbeqaaiaacQca aaGccaGGOaGaaGimaiaacMcacaGGSaGaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caWG0bWaaSbaaSqaaiaaigdaaeqaaOGaaiikaiaadYea caGGPaGaeyypa0JaamiDamaaCaaaleqabaGaaiOkaaaakiaacIcaca WGmbGaaiykaaaa@4B4C@  or displacement u 1 (0)= u * (0) u 1 (L)= u * (L) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaBaaaleaacaaIXaaabeaaki aacIcacaaIWaGaaiykaiabg2da9iaadwhadaahaaWcbeqaaiaacQca aaGccaGGOaGaaGimaiaacMcacaaMc8UaaGPaVlaaykW7caaMc8Uaam yDamaaBaaaleaacaaIXaaabeaakiaacIcacaWGmbGaaiykaiabg2da 9iaaykW7caWG1bWaaWbaaSqabeaacaGGQaaaaOGaaiikaiaadYeaca GGPaaaaa@4AA0@  at x=0 and x=L

 

For this 1-D example, then, the finite element equations reduce to

K ab u 1 b = F a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamOyaa qabaGccaWG1bWaa0baaSqaaiaaigdaaeaacaWGIbaaaOGaeyypa0Ja amOramaaDaaaleaaaeaacaWGHbaaaaaa@3959@

where

K ab = h 2 0 L 2μ(1ν) 12ν N a ( x 1 ) x 1 N b ( x 1 ) x 1 d x 1 F a = h 2 0 L b N a ( x 1 )d x 1 + h 2 t 1 * (0) N a (0)+ h 2 t * (L) N a (L) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGceaqabeaacaWGlbWaaSbaaSqaaiaadggaca WGIbaabeaakiabg2da9iaadIgadaahaaWcbeqaaiaaikdaaaGcdaWd XbqaamaalaaabaGaaGOmaiabeY7aTjaacIcacaaIXaGaeyOeI0Iaeq yVd4MaaiykaaqaaiaaigdacqGHsislcaaIYaGaeqyVd4gaamaalaaa baGaeyOaIyRaamOtamaaCaaaleqabaGaamyyaaaakiaacIcacaWG4b WaaSbaaSqaaiaaigdaaeqaaOGaaiykaaqaaiabgkGi2kaadIhadaWg aaWcbaGaaGymaaqabaaaaOWaaSaaaeaacqGHciITcaWGobWaaWbaaS qabeaacaWGIbaaaOGaaiikaiaadIhadaWgaaWcbaGaaGymaaqabaGc caGGPaaabaGaeyOaIyRaamiEamaaBaaaleaacaaIXaaabeaaaaGcca WGKbGaamiEamaaBaaaleaacaaIXaaabeaaaeaacaaIWaaabaGaamit aaqdcqGHRiI8aaGcbaGaamOramaaDaaaleaaaeaacaWGHbaaaOGaey ypa0JaamiAamaaCaaaleqabaGaaGOmaaaakmaapehabaGaamOyamaa BaaaleaaaeqaaOGaamOtamaaCaaaleqabaGaamyyaaaakiaacIcaca WG4bWaaSbaaSqaaiaaigdaaeqaaOGaaiykaiaadsgacaWG4bWaaSba aSqaaiaaigdaaeqaaaqaaiaaicdaaeaacaWGmbaaniabgUIiYdGccq GHRaWkcaWGObWaaWbaaSqabeaacaaIYaaaaOGaamiDamaaDaaaleaa caaIXaaabaGaaiOkaaaakiaacIcacaaIWaGaaiykaiaad6eadaahaa WcbeqaaiaadggaaaGccaGGOaGaaGimaiaacMcacqGHRaWkcaWGObWa aWbaaSqabeaacaaIYaaaaOGaamiDamaaDaaaleaaaeaacaGGQaaaaO GaaiikaiaadYeacaGGPaGaamOtamaaCaaaleqabaGaamyyaaaakiaa cIcacaWGmbGaaiykaaaaaa@84DF@

We could obviously choose any interpolation scheme, evaluate the necessary integrals and solve the resulting system of equations to compute the solution.  It turns out to be particularly convenient, however, to use a piecewise-Lagrangian interpolation scheme, and to evaluate the integrals numerically using a Gaussian quadrature scheme.

 

To implement the Lagrangian interpolation scheme, we sub-divide the region 0 x 1 L MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaGimaiabgsMiJkaadIhadaWgaaWcba GaaGymaaqabaGccqGHKjYOcaWGmbaaaa@37B2@  into a series of elements, as illustrated in the figure.  Each element is bounded by two nodal points, and may also contain one or more interior nodes. The displacement field within the element is interpolated between the nodes attached to the element.  So, we would use a linear interpolation between the nodes on a 2-noded element, a quadratic interpolation between the nodes on a 3 noded element, and so on.

 

Linear 1-D element

N 1 (ξ)=0.5(1ξ) N 2 (ξ)=0.5(1+ξ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaaiikaiabe67a4jaacMcacqGH9aqpcaaIWaGaaiOlaiaaiwda caGGOaGaaGymaiabgkHiTiabe67a4jaacMcaaeaacaWGobWaaWbaaS qabeaacaaIYaaaaOGaaiikaiabe67a4jaacMcacqGH9aqpcaaIWaGa aiOlaiaaiwdacaGGOaGaaGymaiabgUcaRiabe67a4jaacMcaaaaa@4A78@

Quadratic 1-D element

N 1 (ξ)=0.5ξ(1ξ) N 2 (ξ)=0.5ξ(1+ξ) N 3 (ξ)=(1ξ)(1+ξ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaaiikaiabe67a4jaacMcacqGH9aqpcqGHsislcaaIWaGaaiOl aiaaiwdacqaH+oaEcaGGOaGaaGymaiabgkHiTiabe67a4jaacMcaae aacaWGobWaaWbaaSqabeaacaaIYaaaaOGaaiikaiabe67a4jaacMca cqGH9aqpcaaIWaGaaiOlaiaaiwdacqaH+oaEcaGGOaGaaGymaiabgU caRiabe67a4jaacMcaaeaacaWGobWaaWbaaSqabeaacaaIZaaaaOGa aiikaiabe67a4jaacMcacqGH9aqpcaGGOaGaaGymaiabgkHiTiabe6 7a4jaacMcacaGGOaGaaGymaiabgUcaRiabe67a4jaacMcaaaaa@5E52@

Generic linear and quadrilateral 1-D elements are illustrated in the table.  The local nodes on the element are numbered 1 and 2 for the linear element, and 1,2,3 for the quadratic element as shown.  We suppose that the element lies in the region 1 ξ 1 1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOeI0IaaGymaiabgsMiJkabe67a4n aaBaaaleaacaaIXaaabeaakiabgsMiJkaaigdaaaa@3950@ .  The displacements within the element are then interpolated as

u 1 ( ξ 1 )= a=1 N e N a ( ξ 1 ) u 1 a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaBaaaleaacaaIXaaabeaaki aacIcacqaH+oaEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaeyypa0Za aabCaeaacaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiikaiabe67a4n aaBaaaleaacaaIXaaabeaakiaacMcacaWG1bWaa0baaSqaaiaaigda aeaacaWGHbaaaaqaaiaadggacqGH9aqpcaaIXaaabaGaamOtamaaBa aameaacaWGLbaabeaaa0GaeyyeIuoaaaa@477A@

where N e MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOtamaaBaaaleaacaWGLbaabeaaaa a@32B8@  denotes the number of nodes on the element, u 1 a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaDaaaleaacaaIXaaabaGaam yyaaaaaaa@3397@  denotes the value of the displacement at each node, and the shape functions are given in the table.                

                             
Of course, the actual nodal coordinates do not lie at MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@ 1, +1 and 0 for all the elements.  For a general element, we map this special one to the region of interest.  A particularly convenient way to do this is to set

x 1 ( ξ 1 )= a=1 N e N a (ξ) x 1 a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiEamaaBaaaleaacaaIXaaabeaaki aacIcacqaH+oaEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaeyypa0Za aabCaeaacaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiikaiabe67a4j aacMcacaWG4bWaa0baaSqaaiaaigdaaeaacaWGHbaaaaqaaiaadgga cqGH9aqpcaaIXaaabaGaamOtamaaBaaameaacaWGLbaabeaaa0Gaey yeIuoaaaa@468F@

where x 1 a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiEamaaDaaaleaacaaIXaaabaGaam yyaaaaaaa@33A8@  denotes the coordinate of each node 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 or 3).  Elements that interpolate displacements and position using the same shape functions are called isoparametric elements.

 

Next, we need to devise a way to do the integrals in the expressions for the stiffness matrix and force vector.  We can evidently divide up the integral so as to integrate over each element in turn

K ab = l=1 N lmn h 2 x0 x1 2μ(1ν) 12ν N a ( x 1 ) x 1 N b ( x 1 ) x 1 d x 1 F a = l=1 N lmn h 2 x0 x1 b N a ( x 1 )d x 1 + h 2 t i * N a (0)+ h 2 t i * N a (L) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGceaqabeaacaWGlbWaaSbaaSqaaiaadggaca WGIbaabeaakiabg2da9maaqahabaGaamiAamaaCaaaleqabaGaaGOm aaaakmaapehabaWaaSaaaeaacaaIYaGaeqiVd0Maaiikaiaaigdacq GHsislcqaH9oGBcaGGPaaabaGaaGymaiabgkHiTiaaikdacqaH9oGB aaWaaSaaaeaacqGHciITcaWGobWaaWbaaSqabeaacaWGHbaaaOGaai ikaiaadIhadaWgaaWcbaGaaGymaaqabaGccaGGPaaabaGaeyOaIyRa amiEamaaBaaaleaacaaIXaaabeaaaaGcdaWcaaqaaiabgkGi2kaad6 eadaahaaWcbeqaaiaadkgaaaGccaGGOaGaamiEamaaBaaaleaacaaI XaaabeaakiaacMcaaeaacqGHciITcaWG4bWaaSbaaSqaaiaaigdaae qaaaaakiaadsgacaWG4bWaaSbaaSqaaiaaigdaaeqaaaqaaiaadIha caaIWaaabaGaamiEaiaaigdaa0Gaey4kIipaaSqaaiaadYgacqGH9a qpcaaIXaaabaGaamOtamaaBaaameaacaWGSbGaamyBaiaad6gaaeqa aaqdcqGHris5aaGcbaGaamOramaaDaaaleaaaeaacaWGHbaaaOGaey ypa0ZaaabCaeaacaWGObWaaWbaaSqabeaacaaIYaaaaOWaa8qCaeaa caWGIbWaaSbaaSqaaaqabaGccaWGobWaaWbaaSqabeaacaWGHbaaaO GaaiikaiaadIhadaWgaaWcbaGaaGymaaqabaGccaGGPaGaamizaiaa dIhadaWgaaWcbaGaaGymaaqabaaabaGaamiEaiaaicdaaeaacaWG4b GaaGymaaqdcqGHRiI8aaWcbaGaamiBaiabg2da9iaaigdaaeaacaWG obWaaSbaaWqaaiaadYgacaWGTbGaamOBaaqabaaaniabggHiLdGccq GHRaWkcaWGObWaaWbaaSqabeaacaaIYaaaaOGaamiDamaaDaaaleaa caWGPbaabaGaaiOkaaaakiaad6eadaahaaWcbeqaaiaadggaaaGcca GGOaGaaGimaiaacMcacqGHRaWkcaWGObWaaWbaaSqabeaacaaIYaaa aOGaamiDamaaDaaaleaacaWGPbaabaGaaiOkaaaakiaad6eadaahaa WcbeqaaiaadggaaaGccaGGOaGaamitaiaacMcaaaaa@973F@

where N lmn MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8skY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaad6eadaWgaaWcbaGaamiBaiaad2gaca WGUbaabeaaaaa@346D@  is the total number of elements, and x0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiEaiaaicdaaaa@3294@  and x1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiEaiaaigdaaaa@3295@  denote the coordinates of the ends of the lth element.  We now notice an attractive feature of our interpolation scheme.  The integral over the lth element depends only on the shape functions associated with the nodes on the lth element, since the displacement in this region is completely determined by its values at these nodes.  We can therefore define element stiffness matrix, and element force matrix

k ab = h 2 x0 x1 2μ(1ν) 12ν N a ( x 1 ) x 1 N b ( x 1 ) x 1 d x 1 f a = h 2 x0 x1 b N a ( x 1 )d x 1 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4AamaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpcaWGObWaaWbaaSqabeaacaaIYaaaaOWaa8qCaeaa daWcaaqaaiaaikdacqaH8oqBcaGGOaGaaGymaiabgkHiTiabe27aUj aacMcaaeaacaaIXaGaeyOeI0IaaGOmaiabe27aUbaadaWcaaqaaiab gkGi2kaad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGaamiEamaaBa aaleaacaaIXaaabeaakiaacMcaaeaacqGHciITcaWG4bWaaSbaaSqa aiaaigdaaeqaaaaakmaalaaabaGaeyOaIyRaamOtamaaCaaaleqaba GaamOyaaaakiaacIcacaWG4bWaaSbaaSqaaiaaigdaaeqaaOGaaiyk aaqaaiabgkGi2kaadIhadaWgaaWcbaGaaGymaaqabaaaaOGaamizai aadIhadaWgaaWcbaGaaGymaaqabaaabaGaamiEaiaaicdaaeaacaWG 4bGaaGymaaqdcqGHRiI8aOGaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaadAgadaqhaaWcbaaabaGaamyyaaaakiabg2da 9iaadIgadaahaaWcbeqaaiaaikdaaaGcdaWdXbqaaiaadkgadaWgaa Wcbaaabeaakiaad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGaamiE amaaBaaaleaacaaIXaaabeaakiaacMcacaWGKbGaamiEamaaBaaale aacaaIXaaabeaaaeaacaWG4bGaaGimaaqaaiaadIhacaaIXaaaniab gUIiYdaaaa@7D4E@

for each element, which depend on the geometry, interpolation functions and material properties of the element.  The first and last elements have additional contributions to the element force vector from the boundary terms h 2 t i * N a (0), h 2 t i * N a (L) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaamiAamaaCaaaleqabaGaaGOmaaaaki aadshadaqhaaWcbaGaamyAaaqaaiaacQcaaaGccaWGobWaaWbaaSqa beaacaWGHbaaaOGaaiikaiaaicdacaGGPaGaaiilaiaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaamiAamaaCaaaleqabaGaaGOmaaaakiaa dshadaqhaaWcbaGaamyAaaqaaiaacQcaaaGccaWGobWaaWbaaSqabe aacaWGHbaaaOGaaiikaiaadYeacaGGPaaaaa@4AAA@ .  The global stiffness matrix is computed by summing all the element stiffness matrices

K ab = l=1 N lmn k ab F a = l=1 N lmn f a + h 2 t i * N a (0)+ h 2 t i * N a (L) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpdaaeWbqaaiaadUgadaWgaaWcbaGaamyyaiaadkga aeqaaaqaaiaadYgacqGH9aqpcaaIXaaabaGaamOtamaaBaaameaaca WGSbGaamyBaiaad6gaaeqaaaqdcqGHris5aOGaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaadAeadaqhaaWcbaaabaGaamyyaaaaki abg2da9maaqahabaGaamOzamaaCaaaleqabaGaamyyaaaaaeaacaWG SbGaeyypa0JaaGymaaqaaiaad6eadaWgaaadbaGaamiBaiaad2gaca WGUbaabeaaa0GaeyyeIuoakiabgUcaRiaadIgadaahaaWcbeqaaiaa ikdaaaGccaWG0bWaa0baaSqaaiaadMgaaeaacaGGQaaaaOGaamOtam aaCaaaleqabaGaamyyaaaakiaacIcacaaIWaGaaiykaiabgUcaRiaa dIgadaahaaWcbeqaaiaaikdaaaGccaWG0bWaa0baaSqaaiaadMgaae aacaGGQaaaaOGaamOtamaaCaaaleqabaGaamyyaaaakiaacIcacaWG mbGaaiykaaaa@8359@

 

Finally we need to devise a way to compute the integrals for each element stiffness matrix.  It is convenient to map the domain of integration to [-1,+1] and integrate with respect to the normalized coordinate ξ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdGhaaa@32A0@  - thus

k ab = h 2 1 +1 2μ(1ν) 12ν N a (x) x N b (x) x Jdξ f a = h 2 1 +1 b N a ( x 1 )Jdξ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4AamaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpcaWGObWaaWbaaSqabeaacaaIYaaaaOWaa8qCaeaa daWcaaqaaiaaikdacqaH8oqBcaGGOaGaaGymaiabgkHiTiabe27aUj aacMcaaeaacaaIXaGaeyOeI0IaaGOmaiabe27aUbaadaWcaaqaaiab gkGi2kaad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGaamiEaiaacM caaeaacqGHciITcaWG4baaamaalaaabaGaeyOaIyRaamOtamaaCaaa leqabaGaamOyaaaakiaacIcacaWG4bGaaiykaaqaaiabgkGi2kaadI haaaGaamOsaiaadsgacqaH+oaEaSqaaiabgkHiTiaaigdaaeaacqGH RaWkcaaIXaaaniabgUIiYdGccaaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaamOzamaaDa aaleaaaeaacaWGHbaaaOGaeyypa0JaamiAamaaCaaaleqabaGaaGOm aaaakmaapehabaGaamOyamaaBaaaleaaaeqaaOGaamOtamaaCaaale qabaGaamyyaaaakiaacIcacaWG4bWaaSbaaSqaaiaaigdaaeqaaOGa aiykaiaadQeacaWGKbGaeqOVdGhaleaacqGHsislcaaIXaaabaGaey 4kaSIaaGymaaqdcqGHRiI8aaaa@8D2C@

where J=| x/dξ | MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOsaiabg2da9maaemaabaGaeyOaIy RaamiEaiaac+cacaWGKbGaeqOVdGhacaGLhWUaayjcSdaaaa@3B86@  is the Jacobian associated with the mapping, which may be computed as

x ξ = ξ a=1 N e N a ( ξ ) x a = a=1 N e N a ξ x a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaSaaaeaacqGHciITcaWG4baabaGaey OaIyRaeqOVdGhaaiabg2da9maalaaabaGaeyOaIylabaGaeyOaIyRa eqOVdGhaamaaqahabaGaamOtamaaCaaaleqabaGaamyyaaaakiaacI cacqaH+oaEdaWgaaWcbaaabeaakiaacMcacaWG4bWaa0baaSqaaaqa aiaadggaaaaabaGaamyyaiabg2da9iaaigdaaeaacaWGobWaaSbaaW qaaiaadwgaaeqaaaqdcqGHris5aOGaeyypa0ZaaabCaeaadaWcaaqa aiabgkGi2kaad6eadaahaaWcbeqaaiaadggaaaaakeaacqGHciITcq aH+oaEaaGaamiEamaaDaaaleaaaeaacaWGHbaaaaqaaiaadggacqGH 9aqpcaaIXaaabaGaamOtamaaBaaameaacaWGLbaabeaaa0GaeyyeIu oaaaa@5AC1@

1-D integration points and weights

M=1

ξ (1) =0 w 1 =2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aaWbaaSqabeaacaGGOaGaaG ymaiaacMcaaaGccqGH9aqpcaaIWaGaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 Uaam4DamaaBaaaleaacaaIXaaabeaakiabg2da9iaaikdaaaa@4CD0@

M=2

ξ (1) =0.5773502691 w 1 =1.0 ξ (2) =0.5773502691 w 2 =1.0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacqaH+oaEdaahaaWcbeqaaiaacI cacaaIXaGaaiykaaaakiabg2da9iabgkHiTiaaicdacaGGUaGaaGyn aiaaiEdacaaI3aGaaG4maiaaiwdacaaMc8UaaGimaiaaikdacaaI2a GaaGyoaiaaigdacaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8Uaam4DamaaBaaaleaacaaI Xaaabeaakiabg2da9iaaigdacaGGUaGaaGimaaqaaiabe67a4naaCa aaleqabaGaaiikaiaaikdacaGGPaaaaOGaeyypa0JaaGimaiaac6ca caaI1aGaaG4naiaaiEdacaaIZaGaaGynaiaaykW7caaIWaGaaGOmai aaiAdacaaI5aGaaGymaiaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caWG3bWaaSbaaSqaaiaaikdaaeqaaOGaeyypa0JaaGymaiaac6 cacaaIWaaaaaa@9A34@

M=3

ξ (1) =0.7745966692 w 1 =0.55555555555 ξ (2) =0. w 2 =0.88888888888 ξ (3) =0.7745966692 w 3 =0.55555555555 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacqaH+oaEdaahaaWcbeqaaiaacI cacaaIXaGaaiykaaaakiabg2da9iabgkHiTiaaicdacaGGUaGaaG4n aiaaiEdacaaI0aGaaGynaiaaiMdacaaI2aGaaGOnaiaaiAdacaaI5a GaaGOmaiaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caWG3bWaaSbaaSqaaiaaigdaaeqaaOGaey ypa0JaaGimaiaac6cacaaI1aGaaGynaiaaiwdacaaI1aGaaGynaiaa iwdacaaI1aGaaGynaiaaiwdacaaI1aGaaGynaaqaaiabe67a4naaCa aaleqabaGaaiikaiaaikdacaGGPaaaaOGaeyypa0JaaGimaiaac6ca caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caWG3bWaaSbaaSqaaiaaikdaaeqaaOGaeyypa0JaaGim aiaac6cacaaI4aGaaGioaiaaiIdacaaI4aGaaGioaiaaiIdacaaI4a GaaGioaiaaiIdacaaI4aGaaGioaaqaaiabe67a4naaCaaaleqabaGa aiikaiaaiodacaGGPaaaaOGaeyypa0JaaGimaiaac6cacaaI3aGaaG 4naiaaisdacaaI1aGaaGyoaiaaiAdacaaI2aGaaGOnaiaaiMdacaaI YaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaadEhadaWgaaWcbaGa aG4maaqabaGccqGH9aqpcaaIWaGaaiOlaiaaiwdacaaI1aGaaGynai aaiwdacaaI1aGaaGynaiaaiwdacaaI1aGaaGynaiaaiwdacaaI1aaa aaa@DE74@

 

Note that the mapping also enables us to calculate the shape function derivatives in the element stiffness matrix as

N a x = N a ξ ξ x = N a ξ [ x ξ ] 1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaSaaaeaacqGHciITcaWGobWaaWbaaS qabeaacaWGHbaaaaGcbaGaeyOaIyRaamiEaaaacqGH9aqpdaWcaaqa aiabgkGi2kaad6eadaahaaWcbeqaaiaadggaaaaakeaacqGHciITcq aH+oaEaaWaaSaaaeaacqGHciITcqaH+oaEaeaacqGHciITcaWG4baa aiabg2da9maalaaabaGaeyOaIyRaamOtamaaCaaaleqabaGaamyyaa aaaOqaaiabgkGi2kabe67a4baadaWadaqaamaalaaabaGaeyOaIyRa amiEaaqaaiabgkGi2kabe67a4baaaiaawUfacaGLDbaadaahaaWcbe qaaiabgkHiTiaaigdaaaaaaa@54C1@

 

Finally, note that integrals may be computed numerically using a quadrature formula, as follows

1 +1 g(ξ)dξ I=1 M w I g( ξ (I) ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8qCaeaacaWGNbGaaiikaiabe67a4j aacMcacaWGKbGaeqOVdGhaleaacqGHsislcaaIXaaabaGaey4kaSIa aGymaaqdcqGHRiI8aOGaeyisIS7aaabCaeaacaWG3bWaaSbaaSqaai aadMeaaeqaaOGaam4zaiaacIcacqaH+oaEdaahaaWcbeqaaiaacIca caWGjbGaaiykaaaakiaacMcaaSqaaiaadMeacqGH9aqpcaaIXaaaba GaamytaaqdcqGHris5aaaa@4CF5@

where  ξ (I) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aaWbaaSqabeaacaGGOaGaam ysaiaacMcaaaaaaa@34E6@   I=1…M denotes a set of integration points in the region [-1,+1], and w I MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaam4DamaaBaaaleaacaWGjbaabeaaaa a@32C3@  is a set of integration weights, which are chosen so as to make the approximation as accurate as possible. Values are given in the table to the right for M=1,2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaad2eacqGH9aqpcaaIXaGaaiilaiaaik daaaa@3467@  and 3. Higher order integration schemes exist but are required only for higher order elements.  For the linear 1-D element described earlier a single integration point is sufficient to evaluate the stiffness exactly. Similarly, for the quaratic element, two integration points will suffice.

 

 

8.1.6 Summary of the 1D finite element procedure

 

To summarize, then, the finite element solution requires the following steps:

1.      For each element, compute the element stiffness matrix as follows:

k ab = h 2 I=1 M w I 2μ(1ν) 12ν N a x N b x J( ξ I ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4AamaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpcaWGObWaaWbaaSqabeaacaaIYaaaaOWaaabCaeaa caWG3bWaaSbaaSqaaiaadMeaaeqaaOWaaSaaaeaacaaIYaGaeqiVd0 MaaiikaiaaigdacqGHsislcqaH9oGBcaGGPaaabaGaaGymaiabgkHi TiaaikdacqaH9oGBaaWaaSaaaeaacqGHciITcaWGobWaaWbaaSqabe aacaWGHbaaaaGcbaGaeyOaIyRaamiEaaaadaWcaaqaaiabgkGi2kaa d6eadaahaaWcbeqaaiaadkgaaaaakeaacqGHciITcaWG4baaaiaadQ eacaGGOaGaeqOVdG3aaSbaaSqaaiaadMeaaeqaaOGaaiykaaWcbaGa amysaiabg2da9iaaigdaaeaacaWGnbaaniabggHiLdaaaa@5A43@

where

J=| x ξ |=| a=1 N e N a ( ξ I ) ξ x a | MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOsaiabg2da9maaemaabaWaaSaaae aacqGHciITcaWG4baabaGaeyOaIyRaeqOVdGhaaaGaay5bSlaawIa7 aiabg2da9maaemaabaWaaabCaeaadaWcaaqaaiabgkGi2kaad6eada ahaaWcbeqaaiaadggaaaGccaGGOaGaeqOVdG3aaSbaaSqaaiaadMea aeqaaOGaaiykaaqaaiabgkGi2kabe67a4baacaWG4bWaa0baaSqaaa qaaiaadggaaaaabaGaamyyaiabg2da9iaaigdaaeaacaWGobWaaSba aWqaaiaadwgaaeqaaaqdcqGHris5aaGccaGLhWUaayjcSdaaaa@533B@             N a x = N a ( ξ I ) ξ ξ x = N a ( ξ I ) ξ [ x ξ ] 1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaSaaaeaacqGHciITcaWGobWaaWbaaS qabeaacaWGHbaaaaGcbaGaeyOaIyRaamiEaaaacqGH9aqpdaWcaaqa aiabgkGi2kaad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGaeqOVdG 3aaSbaaSqaaiaadMeaaeqaaOGaaiykaaqaaiabgkGi2kabe67a4baa daWcaaqaaiabgkGi2kabe67a4bqaaiabgkGi2kaadIhaaaGaeyypa0 ZaaSaaaeaacqGHciITcaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiik aiabe67a4naaBaaaleaacaWGjbaabeaakiaacMcaaeaacqGHciITcq aH+oaEaaWaamWaaeaadaWcaaqaaiabgkGi2kaadIhaaeaacqGHciIT cqaH+oaEaaaacaGLBbGaayzxaaWaaWbaaSqabeaacqGHsislcaaIXa aaaaaa@5D01@

and the integration points ξ I , w I I=1...M MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aaSbaaSqaaiaadMeaaeqaaO GaaiilaiaadEhadaWgaaWcbaGaamysaaqabaGccaaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caWGjbGaeyypa0JaaGymaiaac6caca GGUaGaaiOlaiaad2eaaaa@44FD@  are tabulated above, and shape functions N (a) (ξ),a=1... N e MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOtamaaCaaaleqabaGaaiikaiaadg gacaGGPaaaaOGaaiikaiabe67a4jaacMcacaGGSaGaaGPaVlaaykW7 caaMc8UaaGPaVlaadggacqGH9aqpcaaIXaGaaiOlaiaac6cacaGGUa GaamOtamaaBaaaleaacaWGLbaabeaaaaa@44B6@  were listed earlier.

2.      Assemble the contribution from each element to the global stiffness

K ab = l=1 N lmn k ab MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpdaaeWbqaaiaadUgadaWgaaWcbaGaamyyaiaadkga aeqaaaqaaiaadYgacqGH9aqpcaaIXaaabaGaamOtamaaBaaameaaca WGSbGaamyBaiaad6gaaeqaaaqdcqGHris5aaaa@406F@

3.      Similarly, if there is a non-zero body force, then compute for each element

f a = h 2 I=1 M w I b N a ( ξ I )J ( ξ I ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaamOzamaaDaaaleaaaeaacaWGHbaaaO Gaeyypa0JaamiAamaaCaaaleqabaGaaGOmaaaakmaaqahabaGaam4D amaaBaaaleaacaWGjbaabeaakiaadkgadaWgaaWcbaaabeaakiaad6 eadaahaaWcbeqaaiaadggaaaGccaGGOaGaeqOVdG3aaSbaaSqaaiaa dMeaaeqaaOGaaiykaiaadQeaaSqaaiaadMeacqGH9aqpcaaIXaaaba GaamytaaqdcqGHris5aOGaaiikaiabe67a4naaBaaaleaacaWGjbaa beaakiaacMcaaaa@49A5@

and assemble the global force vector

F a = l=1 N lmn f a MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaamOramaaDaaaleaaaeaacaWGHbaaaO Gaeyypa0ZaaabCaeaacaWGMbWaaWbaaSqabeaacaWGHbaaaaqaaiaa dYgacqGH9aqpcaaIXaaabaGaamOtamaaBaaameaacaWGSbGaamyBai aad6gaaeqaaaqdcqGHris5aaaa@3E99@

4.      Add contributions to the force vector from prescribed traction boundary conditions at x=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiEaiabg2da9iaaicdaaaa@339A@  and x=L MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiEaiabg2da9iaadYeaaaa@33B1@

F 1 = F 1 + h 2 t i * (0) F (L) = F (L) + h 2 t i * (L) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaamOramaaCaaaleqabaGaaGymaaaaki abg2da9iaadAeadaahaaWcbeqaaiaaigdaaaGccqGHRaWkcaWGObWa aWbaaSqabeaacaaIYaaaaOGaamiDamaaDaaaleaacaWGPbaabaGaai OkaaaakiaacIcacaaIWaGaaiykaiaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWGgbWa aWbaaSqabeaacaGGOaGaamitaiaacMcaaaGccqGH9aqpcaWGgbWaaW baaSqabeaacaGGOaGaamitaiaacMcaaaGccqGHRaWkcaWGObWaaWba aSqabeaacaaIYaaaaOGaamiDamaaDaaaleaacaWGPbaabaGaaiOkaa aakiaacIcacaWGmbGaaiykaaaa@6956@

where the (L) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaiikaiaadYeacaGGPaaaaa@3307@  superscript denotes the node that lies at x=L.

5.      Modify the stiffness matrix to enforce the constraints

u 1 = u * (0) u (L) = u * (L) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaCaaaleqabaGaaGymaaaaki abg2da9iaadwhadaahaaWcbeqaaiaacQcaaaGccaGGOaGaaGimaiaa cMcacaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caWG1bWaaWbaaSqabeaacaGGOaGaamitaiaacM caaaGccqGH9aqpcaWG1bWaaWbaaSqabeaacaGGQaaaaOGaaiikaiaa dYeacaGGPaaaaa@5BE3@

6.      Solve the system of linear equations

K ab u 1 b = F a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamOyaa qabaGccaWG1bWaa0baaSqaaiaaigdaaeaacaWGIbaaaOGaeyypa0Ja amOramaaDaaaleaaaeaacaWGHbaaaaaa@3959@

for the unknown displacements u 1 b MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadwhadaqhaaWcbaGaaGymaaqaaiaadk gaaaaaaa@3351@

 

 

 

 

8.1.7 Example FEM Code and solution

 

 

 

 

 

 

A simple example MAPLE code for this 1-D example can be found in the file FEM_1D_Static.mws

 

It is set up to solve for displacements for a bar with the following parameters:

1.      Length L=5, unit x-sect area,

2.      Shear modulus 50, Poisson’s ratio 0.3,

3.      Uniform body force magnitude 10,

4.      Displacement u=0 at x=0

5.      Traction t=2 at x=L.

 

The code computes the (1D) displacement distribution in the bar. The predicted displacement field is plotted on the right.

 

Of course in general we want to calculate more than just displacements MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  usually we want the stress field too. We can calculate the stress field anywhere within an element  by differentiating the displacements to calculate strains, and then substituting into the constitutive relation.  This gives

σ 11 = 2μ(1ν) (12ν) a=1 n N a x u a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4Wdm3aaSbaaSqaaiaaigdacaaIXa aabeaakiabg2da9maalaaabaGaaGOmaiabeY7aTjaacIcacaaIXaGa eyOeI0IaeqyVd4MaaiykaaqaaiaacIcacaaIXaGaeyOeI0IaaGOmai abe27aUjaacMcaaaWaaabCaeaadaWcaaqaaiabgkGi2kaad6eadaah aaWcbeqaaiaadggaaaaakeaacqGHciITcaWG4baaaaWcbaGaamyyai abg2da9iaaigdaaeaacaWGUbaaniabggHiLdGccaWG1bWaaWbaaSqa beaacaWGHbaaaaaa@4FB0@

This works well for a uniform body force with quadratic (3 noded elements) as the plot on the right shows.

 

However, if we switch to linear elements, the stress results are not so good (displacements are still calculated exactly).  In this case, the stress must be uniform in each element (because strains are constant for linear displacement field), so the stress plot looks like the figure to the right.


Notice that the stresses are most accurate near the center of each element (at the integration point).  For this reason, FEM codes generally output stress and strain data at integration points. 

 

It is interesting also to examine the stiffness matrix (shown below for 3 linear elements, before addition of the u=0 constraint for the first node )

[ 105 105 0 0 105 210 105 0 0 105 210 105 0 0 105 105 ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaamWaaeaafaqabeabeaaaaaqaaiaaig dacaaIWaGaaGynaaqaaiabgkHiTiaaigdacaaIWaGaaGynaaqaaiaa icdaaeaacaaIWaaabaGaeyOeI0IaaGymaiaaicdacaaI1aaabaGaaG OmaiaaigdacaaIWaaabaGaeyOeI0IaaGymaiaaicdacaaI1aaabaGa aGimaaqaaiaaicdaaeaacqGHsislcaaIXaGaaGimaiaaiwdaaeaaca aIYaGaaGymaiaaicdaaeaacqGHsislcaaIXaGaaGimaiaaiwdaaeaa caaIWaaabaGaaGimaaqaaiabgkHiTiaaigdacaaIWaGaaGynaaqaai aaigdacaaIWaGaaGynaaaaaiaawUfacaGLDbaaaaa@52CD@

Notice that stiffness is symmetric, as expected, and also banded.  A large FEM matrix is sparse MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  most of the elements are zero.  This allows the matrix to be stored in compact form MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  for very large matrices indexed storage (where only the nonzero elements together with their indices are stored) is the best approach; for smaller problems skyline storage or band storage (where only the central, mostly nonzero, band of the matrix is stored) may be preferable.  In this case equation numbers need to be assigned to each degree of freedom so as to minimize the bandwidth of the stiffness matrix.

 

 

 

8.1.8 Extending the 1D finite element method to 2 and 3 dimensions

 

It is straightforward to extend the 1-D case to more general problems.  All the basic ideas remain the same.   Specifically

1.      In both 2D and 3D we divide up our solid of interest into a number of elements, shown schematically for a 2D region in the picture on the right.

2.      We define interpolation functions N a ( ξ i ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaad6eadaahaaWcbeqaaiaadggaaaGcca GGOaGaeqOVdG3aaSbaaSqaaiaadMgaaeqaaOGaaiykaaaa@36B8@  for each element in terms of a local, dimensionless, coordinate system within the element.  The coordinates satisfy 1 ξ i +1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkHiTiaaigdacqGHKjYOcqaH+oaEda WgaaWcbaGaamyAaaqabaGccqGHKjYOcqGHRaWkcaaIXaaaaa@3A1E@ .  The displacement field and the position of a point inside an element are computed in terms of the interpolation functions as

u i = a=1 N e N a ( ξ j ) u 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 aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaamiEamaaBaaaleaacaWG Pbaabeaakiabg2da9maaqahabaGaamOtamaaCaaaleqabaGaamyyaa aakiaacIcacqaH+oaEdaWgaaWcbaGaamOAaaqabaGccaGGPaGaamiE amaaDaaaleaacaWGPbaabaGaamyyaaaaaeaacaWGHbGaeyypa0JaaG ymaaqaaiaad6eadaWgaaadbaGaamyzaaqabaaaniabggHiLdaaaa@61FC@

where N a ( ξ j ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOtamaaCaaaleqabaGaamyyaaaaki aacIcacqaH+oaEdaWgaaWcbaGaamOAaaqabaGccaGGPaaaaa@3700@  denote the shape functions, 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.

3.      We introduce an element stiffness matrix for each element by defining

k aibk (l) = V e (l) C ijkl N a (x) x j N b (x) x l dV f i a (l) = V e (l) b i N a (x) dV+ 2 V e (l) t i * N a (x) dA 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 PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaadA gadaqhaaWcbaGaamyAaaqaaiaadggaaaGcdaahaaWcbeqaaiaacIca caWGSbGaaiykaaaakiabg2da9maapefabaGaamOyamaaBaaaleaaca WGPbaabeaakiaad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGaaCiE aiaacMcaaSqaaiaadAfadaqhaaadbaGaamyzaaqaaiaacIcacaWGSb GaaiykaaaaaSqab0Gaey4kIipakiaadsgacaWGwbGaey4kaSYaa8qu aeaacaWG0bWaa0baaSqaaiaadMgaaeaacaGGQaaaaOGaamOtamaaCa aaleqabaGaamyyaaaakiaacIcacaWH4bGaaiykaaWcbaGaeyOaIy7a aSbaaWqaaiaaikdaaeqaaSGaamOvamaaDaaameaacaWGLbaabaGaai ikaiaadYgacaGGPaaaaaWcbeqdcqGHRiI8aOGaamizaiaadgeacaaM c8oaaa@9C54@

where k aibk (l) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadUgadaqhaaWcbaGaamyyaiaadMgaca WGIbGaam4AaaqaaiaacIcacaWGSbGaaiykaaaaaaa@379A@  denotes the element stiffness matrix for the (lth) element, and V e (l) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadAfadaqhaaWcbaGaamyzaaqaaiaacI cacaWGSbGaaiykaaaaaaa@34C4@  denotes the volume (in 3D) or area (in 2D) of the (lth) element, while 2 V e (l) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaWccqGHciITdaWgaaadbaGaaGOmaaqaba WccaWGwbWaa0baaWqaaiaadwgaaeaacaGGOaGaamiBaiaacMcaaaaa aa@3790@  denotes the surface of the (lth) element

4.      The volume integrals over each element are calculated by expressing the volume or surface integral in terms of the dimensionless coordinates 1 ξ i +1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkHiTiaaigdacqGHKjYOcqaH+oaEda WgaaWcbaGaamyAaaqabaGccqGHKjYOcqGHRaWkcaaIXaaaaa@3A1E@ , and then evaluating the integrals numerically, using a quadrature formula of the form

V e f( ξ i )d V ξ = I=1 N I w I f( ξ i I ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaeaacaWGMbGaaiikaiabe67a4n aaBaaaleaacaWGPbaabeaakiaacMcacaWGKbGaamOvamaaBaaaleaa cqaH+oaEaeqaaaqaaiaadAfadaWgaaadbaGaamyzaaqabaaaleqani abgUIiYdGccqGH9aqpdaaeWbqaaiaadEhadaWgaaWcbaGaamysaaqa baGccaWGMbGaaiikaiabe67a4naaDaaaleaacaWGPbaabaGaamysaa aakiaacMcaaSqaaiaadMeacqGH9aqpcaaIXaaabaGaamOtamaaBaaa meaacaWGjbaabeaaa0GaeyyeIuoaaaa@4DB2@

Here, w I MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadEhadaWgaaWcbaGaamysaaqabaaaaa@327E@  are a set of I=1... N I MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadMeacqGH9aqpcaaIXaGaaiOlaiaac6 cacaGGUaGaamOtamaaBaaaleaacaWGjbaabeaaaaa@36FA@  integration weights (just numbers), and ξ i (I) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabe67a4naaDaaaleaacaWGPbaabaGaai ikaiaadMeacaGGPaaaaaaa@358D@  are a set of coordinates that are selected to make the integration scheme as accurate as possible (also just numbers).

5.      The global stiffness matrix

K aibk = R C ijkl N a (x) x j N b (x) x l dV F i a = R b i N a (x) dV+ 2R t i * N a (x) dA MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamyAai aadkgacaWGRbaabeaakiabg2da9maapefabaGaam4qamaaBaaaleaa caWGPbGaamOAaiaadUgacaWGSbaabeaakmaalaaabaGaeyOaIyRaam OtamaaCaaaleqabaGaamyyaaaakiaacIcacaWH4bGaaiykaaqaaiab gkGi2kaadIhadaWgaaWcbaGaamOAaaqabaaaaOWaaSaaaeaacqGHci ITcaWGobWaaWbaaSqabeaacaWGIbaaaOGaaiikaiaahIhacaGGPaaa baGaeyOaIyRaamiEamaaBaaaleaacaWGSbaabeaaaaGccaWGKbGaam OvaaWcbaGaamOuaaqab0Gaey4kIipakiaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaadAeadaqh aaWcbaGaamyAaaqaaiaadggaaaGccqGH9aqpdaWdrbqaaiaadkgada WgaaWcbaGaamyAaaqabaGccaWGobWaaWbaaSqabeaacaWGHbaaaOGa aiikaiaahIhacaGGPaaaleaacaWGsbaabeqdcqGHRiI8aOGaamizai aadAfacqGHRaWkdaWdrbqaaiaadshadaqhaaWcbaGaamyAaaqaaiaa cQcaaaGccaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiikaiaahIhaca GGPaaaleaacqGHciITcaaIYaGaamOuaaqab0Gaey4kIipakiaadsga caWGbbGaaGPaVdaa@8E28@

is then computed by summing the contribution from each element as

K aibk = l=1 N lmn k aibk (l) F i a = l=1 N lmn f i a (l) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamyAai aadkgacaWGRbaabeaakiabg2da9maaqahabaGaam4AamaaDaaaleaa caWGHbGaamyAaiaadkgacaWGRbaabaGaaiikaiaadYgacaGGPaaaaa qaaiaadYgacqGH9aqpcaaIXaaabaGaamOtamaaBaaameaacaWGSbGa amyBaiaad6gaaeqaaaqdcqGHris5aOGaaGzaVlaaygW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaamOramaaDa aaleaacaWGPbaabaGaamyyaaaakiabg2da9maaqahabaGaamOzamaa BaaaleaacaWGPbaabeaakmaaCaaaleqabaGaamyyaaaakmaaCaaale qabaGaaiikaiaadYgacaGGPaaaaaqaaiaadYgacqGH9aqpcaaIXaaa baGaamOtamaaBaaameaacaWGSbGaamyBaiaad6gaaeqaaaqdcqGHri s5aaaa@7482@

6.      The stiffness matrix is modified to enforce any prescribed displacements

7.      The system of equations

K aibk u k b = F i a {a,i}: x k a not on  1 R u i a = u i * ( x i a ){a,i}: x k a on  1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGlbWaaSbaaSqaaiaadggaca WGPbGaamOyaiaadUgaaeqaaOGaamyDamaaDaaaleaacaWGRbaabaGa amOyaaaakiabg2da9iaadAeadaqhaaWcbaGaamyAaaqaaiaadggaaa GccaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaeyiaIiIaai4EaiaadggacaGGSaGaamyA aiaac2hacaaMc8UaaGPaVlaaykW7caGG6aGaaGPaVlaaykW7caaMc8 UaamiEamaaDaaaleaacaWGRbaabaGaamyyaaaakiaaykW7caqGUbGa ae4BaiaabshacaqGGaGaae4Baiaab6gacaqGGaGaeyOaIy7aaSbaaS qaaiaabgdaaeqaaOGaamOuaaqaaiaadwhadaqhaaWcbaGaamyAaaqa aiaadggaaaGccqGH9aqpcaWG1bWaa0baaSqaaiaadMgaaeaacaGGQa aaaOGaaiikaiaadIhadaqhaaWcbaGaamyAaaqaaiaadggaaaGccaGG PaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaeyiaIiIaai4Eaiaadg gacaGGSaGaamyAaiaac2hacaaMc8UaaGPaVlaaykW7caGG6aGaaGPa VlaaykW7caWG4bWaa0baaSqaaiaadUgaaeaacaWGHbaaaOGaaGPaVl aaykW7caqGVbGaaeOBaiaabccacqGHciITdaWgaaWcbaGaaeymaaqa baGccaWGsbaaaaa@BC05@

is solved for the unknown nodal displacements.

8.      The stresses and strains within each element are then deduced.

 

To implement this procedure, we must (a) Define the element interpolation functions; (b) Express the integrals for the element stiffness matrices and force vectors in terms of normalized coordinates; (c) Formulate a numerical integration scheme to evaluate the element stiffness matrices and force vectors.

These details are addressed in the sections to follow.

 

 

 

8.1.9 Interpolation functions for 2D elements

 

The 2D interpolation functions listed below are defined for the region

1 ξ 1 11 ξ 2 1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8skY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkHiTiaaigdacqGHKjYOcqaH+oaEda WgaaWcbaGaaGymaaqabaGccqGHKjYOcaaIXaGaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaeyOeI0IaaGymaiabgsMiJkabe67a4naaBaaaleaacaaI YaaabeaakiabgsMiJkaaigdaaaa@541F@

The numbers shown inside the element show the convention used to number the element faces.

 

2D interpolation functions

N 1 = ξ 1 N 2 = ξ 2 N 3 =1 ξ 1 ξ 2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaeyypa0JaeqOVdG3aaSbaaSqaaiaaigdaaeqaaOGaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaad6eadaahaaWcbe qaaiaaikdaaaGccqGH9aqpcqaH+oaEdaWgaaWcbaGaaGOmaaqabaaa keaacaWGobWaaWbaaSqabeaacaaIZaaaaOGaeyypa0JaaGymaiabgk HiTiabe67a4naaBaaaleaacaaIXaaabeaakiabgkHiTiabe67a4naa BaaaleaacaaIYaaabeaaaaaa@5172@

N 1 =( 2 ξ 1 1 ) ξ 1 N 2 =( 2 ξ 2 1 ) ξ 2 N 3 =( 2(1 ξ 1 ξ 2 )1 )(1 ξ 1 ξ 2 ) N 4 =4 ξ 1 ξ 2 N 5 =4 ξ 2 (1 ξ 1 ξ 2 ) N 6 =4 ξ 1 (1 ξ 1 ξ 2 ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaeyypa0ZaaeWaaeaacaaIYaGaeqOVdG3aaSbaaSqaaiaaigda aeqaaOGaeyOeI0IaaGymaaGaayjkaiaawMcaaiabe67a4naaBaaale aacaaIXaaabeaakiaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaamOt amaaCaaaleqabaGaaGOmaaaakiabg2da9maabmaabaGaaGOmaiabe6 7a4naaBaaaleaacaaIYaaabeaakiabgkHiTiaaigdaaiaawIcacaGL PaaacqaH+oaEdaWgaaWcbaGaaGOmaaqabaaakeaacaWGobWaaWbaaS qabeaacaaIZaaaaOGaeyypa0ZaaeWaaeaacaaIYaGaaiikaiaaigda cqGHsislcqaH+oaEdaWgaaWcbaGaaGymaaqabaGccqGHsislcqaH+o aEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaeyOeI0IaaGymaaGaayjk aiaawMcaaiaacIcacaaIXaGaeyOeI0IaeqOVdG3aaSbaaSqaaiaaig daaeqaaOGaeyOeI0IaeqOVdG3aaSbaaSqaaiaaikdaaeqaaOGaaiyk aaqaaiaad6eadaahaaWcbeqaaiaaisdaaaGccqGH9aqpcaaI0aGaeq OVdG3aaSbaaSqaaiaaigdaaeqaaOGaeqOVdG3aaSbaaSqaaiaaikda aeqaaOGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaamOtam aaCaaaleqabaGaaGynaaaakiabg2da9iaaisdacqaH+oaEdaWgaaWc baGaaGOmaaqabaGccaGGOaGaaGymaiabgkHiTiabe67a4naaBaaale aacaaIXaaabeaakiabgkHiTiabe67a4naaBaaaleaacaaIYaaabeaa kiaacMcaaeaacaWGobWaaWbaaSqabeaacaaI2aaaaOGaeyypa0JaaG inaiabe67a4naaBaaaleaacaaIXaaabeaakiaacIcacaaIXaGaeyOe I0IaeqOVdG3aaSbaaSqaaiaaigdaaeqaaOGaeyOeI0IaeqOVdG3aaS baaSqaaiaaikdaaeqaaOGaaiykaaaaaa@9BA3@

N 1 =0.25(1 ξ 1 )(1 ξ 2 ) N 2 =0.25(1+ ξ 1 )(1 ξ 2 ) N 3 =0.25(1+ ξ 1 )(1+ ξ 2 ) N 4 =0.25(1 ξ 1 )(1+ ξ 2 ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaeyypa0JaaGimaiaac6cacaaIYaGaaGynaiaacIcacaaIXaGa eyOeI0IaeqOVdG3aaSbaaSqaaiaaigdaaeqaaOGaaiykaiaacIcaca aIXaGaeyOeI0IaeqOVdG3aaSbaaSqaaiaaikdaaeqaaOGaaiykaaqa aiaad6eadaahaaWcbeqaaiaaikdaaaGccqGH9aqpcaaIWaGaaiOlai aaikdacaaI1aGaaiikaiaaigdacqGHRaWkcqaH+oaEdaWgaaWcbaGa aGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+oaEdaWgaa WcbaGaaGOmaaqabaGccaGGPaaabaGaamOtamaaCaaaleqabaGaaG4m aaaakiabg2da9iaaicdacaGGUaGaaGOmaiaaiwdacaGGOaGaaGymai abgUcaRiabe67a4naaBaaaleaacaaIXaaabeaakiaacMcacaGGOaGa aGymaiabgUcaRiabe67a4naaBaaaleaacaaIYaaabeaakiaacMcaae aacaWGobWaaWbaaSqabeaacaaI0aaaaOGaeyypa0JaaGimaiaac6ca caaIYaGaaGynaiaacIcacaaIXaGaeyOeI0IaeqOVdG3aaSbaaSqaai aaigdaaeqaaOGaaiykaiaacIcacaaIXaGaey4kaSIaeqOVdG3aaSba aSqaaiaaikdaaeqaaOGaaiykaaaaaa@7526@

N 1 =(1 ξ 1 )(1 ξ 2 )(1+ ξ 1 + ξ 2 )/4 N 2 =(1+ ξ 1 )(1 ξ 2 )( ξ 1 ξ 2 1)/4 N 3 =(1+ ξ 1 )(1+ ξ 2 )( ξ 1 + ξ 2 1)/4 N 4 =(1 ξ 1 )(1+ ξ 2 )( ξ 2 ξ 1 1)/4 N 5 =(1 ξ 1 2 )(1 ξ 2 )/2 N 6 =(1+ ξ 1 )(1 ξ 2 2 )/2 N 7 =(1 ξ 1 2 )(1+ ξ 2 )/2 N 8 =(1 ξ 1 )(1 ξ 2 2 )/2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaeyypa0JaeyOeI0IaaiikaiaaigdacqGHsislcqaH+oaEdaWg aaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+o aEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGHRaWk cqaH+oaEdaWgaaWcbaGaaGymaaqabaGccqGHRaWkcqaH+oaEdaWgaa WcbaGaaGOmaaqabaGccaGGPaGaai4laiaaisdaaeaacaWGobWaaWba aSqabeaacaaIYaaaaOGaeyypa0JaaiikaiaaigdacqGHRaWkcqaH+o aEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsisl cqaH+oaEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiabe67a4n aaBaaaleaacaaIXaaabeaakiabgkHiTiabe67a4naaBaaaleaacaaI YaaabeaakiabgkHiTiaaigdacaGGPaGaai4laiaaisdaaeaacaWGob WaaWbaaSqabeaacaaIZaaaaOGaeyypa0JaaiikaiaaigdacqGHRaWk cqaH+oaEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacq GHRaWkcqaH+oaEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiab e67a4naaBaaaleaacaaIXaaabeaakiabgUcaRiabe67a4naaBaaale aacaaIYaaabeaakiabgkHiTiaaigdacaGGPaGaai4laiaaisdaaeaa caWGobWaaWbaaSqabeaacaaI0aaaaOGaeyypa0Jaaiikaiaaigdacq GHsislcqaH+oaEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaa igdacqGHRaWkcqaH+oaEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaai ikaiabe67a4naaBaaaleaacaaIYaaabeaakiabgkHiTiabe67a4naa BaaaleaacaaIXaaabeaakiabgkHiTiaaigdacaGGPaGaai4laiaais daaeaacaWGobWaaWbaaSqabeaacaaI1aaaaOGaeyypa0Jaaiikaiaa igdacqGHsislcqaH+oaEdaqhaaWcbaGaaGymaaqaaiaaikdaaaGcca GGPaGaaiikaiaaigdacqGHsislcqaH+oaEdaWgaaWcbaGaaGOmaaqa baGccaGGPaGaai4laiaaikdacaaMc8UaaGPaVlaaykW7caaMc8Uaam OtamaaCaaaleqabaGaaGOnaaaakiabg2da9iaacIcacaaIXaGaey4k aSIaeqOVdG3aa0baaSqaaiaaigdaaeaaaaGccaGGPaGaaiikaiaaig dacqGHsislcqaH+oaEdaqhaaWcbaGaaGOmaaqaaiaaikdaaaGccaGG PaGaai4laiaaikdaaeaacaWGobWaaWbaaSqabeaacaaI3aaaaOGaey ypa0JaaiikaiaaigdacqGHsislcqaH+oaEdaqhaaWcbaGaaGymaaqa aiaaikdaaaGccaGGPaGaaiikaiaaigdacqGHRaWkcqaH+oaEdaWgaa WcbaGaaGOmaaqabaGccaGGPaGaai4laiaaikdacaaMc8UaaGPaVlaa ykW7caaMc8UaamOtamaaCaaaleqabaGaaGioaaaakiabg2da9iaacI cacaaIXaGaeyOeI0IaeqOVdG3aa0baaSqaaiaaigdaaeaaaaGccaGG PaGaaiikaiaaigdacqGHsislcqaH+oaEdaqhaaWcbaGaaGOmaaqaai aaikdaaaGccaGGPaGaai4laiaaikdaaaaa@E35A@

 

                                           

 

 

8.1.10 Interpolation Functions for 3D elements

 

The 3D interpolation functions listed below are defined for the region 1 ξ i 1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8skY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkHiTiaaigdacqGHKjYOcqaH+oaEda WgaaWcbaGaamyAaaqabaGccqGHKjYOcaaIXaaaaa@394C@

                                         

3D Interpolation Functions

N 1 = ξ 1 N 2 = ξ 2 N 3 = ξ 3 N 4 =1 ξ 1 ξ 2 ξ 3 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaeyypa0JaeqOVdG3aaSbaaSqaaiaaigdaaeqaaOGaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caWGobWaaWbaaSqabeaacaaIYaaaaOGaeyypa0JaeqOV dG3aaSbaaSqaaiaaikdaaeqaaaGcbaGaamOtamaaCaaaleqabaGaaG 4maaaakiabg2da9iabe67a4naaBaaaleaacaaIZaaabeaakiaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaamOtam aaCaaaleqabaGaaGinaaaakiabg2da9iaaigdacqGHsislcqaH+oaE daWgaaWcbaGaaGymaaqabaGccqGHsislcqaH+oaEdaWgaaWcbaGaaG OmaaqabaGccqGHsislcqaH+oaEdaWgaaWcbaGaaG4maaqabaaaaaa@6D0F@

N 1 =( 2 ξ 1 1 ) ξ 1 N 2 =( 2 ξ 2 1 ) ξ 2 N 3 =( 2 ξ 3 1 ) ξ 3 N 4 =( 2 ξ 4 1 ) ξ 4 N 5 =4 ξ 1 ξ 2 N 6 =4 ξ 2 ξ 3 N 7 =4 ξ 3 ξ 1 N 8 =4 ξ 4 ξ 1 ξ 4 =1 ξ 1 ξ 2 ξ 3 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaeyypa0ZaaeWaaeaacaaIYaGaeqOVdG3aaSbaaSqaaiaaigda aeqaaOGaeyOeI0IaaGymaaGaayjkaiaawMcaaiabe67a4naaBaaale aacaaIXaaabeaakiaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaad6eadaahaaWcbeqaaiaaikdaaaGccq GH9aqpdaqadaqaaiaaikdacqaH+oaEdaWgaaWcbaGaaGOmaaqabaGc cqGHsislcaaIXaaacaGLOaGaayzkaaGaeqOVdG3aaSbaaSqaaiaaik daaeqaaaGcbaGaamOtamaaCaaaleqabaGaaG4maaaakiabg2da9maa bmaabaGaaGOmaiabe67a4naaBaaaleaacaaIZaaabeaakiabgkHiTi aaigdaaiaawIcacaGLPaaacqaH+oaEdaWgaaWcbaGaaG4maaqabaGc caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWGobWaaWbaaS qabeaacaaI0aaaaOGaeyypa0ZaaeWaaeaacaaIYaGaeqOVdG3aaSba aSqaaiaaisdaaeqaaOGaeyOeI0IaaGymaaGaayjkaiaawMcaaiabe6 7a4naaBaaaleaacaaI0aaabeaaaOqaaiaad6eadaahaaWcbeqaaiaa iwdaaaGccqGH9aqpcaaI0aGaeqOVdG3aaSbaaSqaaiaaigdaaeqaaO GaeqOVdG3aaSbaaSqaaiaaikdaaeqaaOGaaGPaVlaaykW7caaMc8Ua aGPaVlaad6eadaahaaWcbeqaaiaaiAdaaaGccqGH9aqpcaaI0aGaeq OVdG3aaSbaaSqaaiaaikdaaeqaaOGaeqOVdG3aaSbaaSqaaiaaioda aeqaaaGcbaGaamOtamaaCaaaleqabaGaaG4naaaakiabg2da9iaais dacqaH+oaEdaWgaaWcbaGaaG4maaqabaGccqaH+oaEdaWgaaWcbaGa aGymaaqabaGccaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaad6eada ahaaWcbeqaaiaaiIdaaaGccqGH9aqpcaaI0aGaeqOVdG3aaSbaaSqa aiaaisdaaeqaaOGaeqOVdG3aaSbaaSqaaiaaigdaaeqaaaGcbaGaeq OVdG3aaSbaaSqaaiaaisdaaeqaaOGaeyypa0JaaGymaiabgkHiTiab e67a4naaBaaaleaacaaIXaaabeaakiabgkHiTiabe67a4naaBaaale aacaaIYaaabeaakiabgkHiTiabe67a4naaBaaaleaacaaIZaaabeaa aaaa@B9AA@

N 1 =(1 ξ 1 )(1 ξ 2 )(1 ξ 3 )/8 N 2 =(1+ ξ 1 )(1 ξ 2 )(1 ξ 3 )/8 N 3 =(1+ ξ 1 )(1+ ξ 2 )(1 ξ 3 )/8 N 4 =(1 ξ 1 )(1+ ξ 2 )(1 ξ 3 )/8 N 5 =(1 ξ 1 )(1 ξ 2 )(1+ ξ 3 )/8 N 6 =(1+ ξ 1 )(1 ξ 2 )(1+ ξ 3 )/8 N 7 =(1+ ξ 1 )(1+ ξ 2 )(1+ ξ 3 )/8 N 8 =(1 ξ 1 )(1+ ξ 2 )(1+ ξ 3 )/8 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaeyypa0JaaiikaiaaigdacqGHsislcqaH+oaEdaWgaaWcbaGa aGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+oaEdaWgaa WcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+oaE daWgaaWcbaGaaG4maaqabaGccaGGPaGaai4laiaaiIdacaaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaamOtamaaCaaaleqa baGaaGOmaaaakiabg2da9iaacIcacaaIXaGaey4kaSIaeqOVdG3aaS baaSqaaiaaigdaaeqaaOGaaiykaiaacIcacaaIXaGaeyOeI0IaeqOV dG3aaSbaaSqaaiaaikdaaeqaaOGaaiykaiaacIcacaaIXaGaeyOeI0 IaeqOVdG3aaSbaaSqaaiaaiodaaeqaaOGaaiykaiaac+cacaaI4aaa baGaamOtamaaCaaaleqabaGaaG4maaaakiabg2da9iaacIcacaaIXa Gaey4kaSIaeqOVdG3aaSbaaSqaaiaaigdaaeqaaOGaaiykaiaacIca caaIXaGaey4kaSIaeqOVdG3aaSbaaSqaaiaaikdaaeqaaOGaaiykai aacIcacaaIXaGaeyOeI0IaeqOVdG3aaSbaaSqaaiaaiodaaeqaaOGa aiykaiaac+cacaaI4aGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caWGobWaaWbaaSqabeaacaaI0aaaaOGaeyyp a0JaaiikaiaaigdacqGHsislcqaH+oaEdaWgaaWcbaGaaGymaaqaba GccaGGPaGaaiikaiaaigdacqGHRaWkcqaH+oaEdaWgaaWcbaGaaGOm aaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+oaEdaWgaaWcba GaaG4maaqabaGccaGGPaGaai4laiaaiIdaaeaacaWGobWaaWbaaSqa beaacaaI1aaaaOGaeyypa0JaaiikaiaaigdacqGHsislcqaH+oaEda WgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH +oaEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGHRa WkcqaH+oaEdaWgaaWcbaGaaG4maaqabaGccaGGPaGaaGPaVlaac+ca caaI4aGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWGobWaaWbaaS qabeaacaaI2aaaaOGaeyypa0JaaiikaiaaigdacqGHRaWkcqaH+oaE daWgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcq aH+oaEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGH RaWkcqaH+oaEdaWgaaWcbaGaaG4maaqabaGccaGGPaGaai4laiaaiI daaeaacaWGobWaaWbaaSqabeaacaaI3aaaaOGaeyypa0Jaaiikaiaa igdacqGHRaWkcqaH+oaEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaai ikaiaaigdacqGHRaWkcqaH+oaEdaWgaaWcbaGaaGOmaaqabaGccaGG PaGaaiikaiaaigdacqGHRaWkcqaH+oaEdaWgaaWcbaGaaG4maaqaba GccaGGPaGaaGPaVlaac+cacaaI4aGaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaamOtamaaCaaaleqabaGaaGioaaaakiabg2da9i aacIcacaaIXaGaeyOeI0IaeqOVdG3aaSbaaSqaaiaaigdaaeqaaOGa aiykaiaacIcacaaIXaGaey4kaSIaeqOVdG3aaSbaaSqaaiaaikdaae qaaOGaaiykaiaacIcacaaIXaGaey4kaSIaeqOVdG3aaSbaaSqaaiaa iodaaeqaaOGaaiykaiaac+cacaaI4aaaaaa@06B4@

N 1 =(1 ξ 1 )(1 ξ 2 )(1 ξ 3 )( ξ 1 ξ 2 ξ 3 2)/8 N 2 =(1+ ξ 1 )(1 ξ 2 )(1 ξ 3 )( ξ 1 ξ 2 ξ 3 2)/8 N 3 =(1+ ξ 1 )(1+ ξ 2 )(1 ξ 3 )( ξ 1 + ξ 2 ξ 3 2)/8 N 4 =(1 ξ 1 )(1+ ξ 2 )(1 ξ 3 )( ξ 1 + ξ 2 ξ 3 2)/8 N 5 =(1 ξ 1 )(1 ξ 2 )(1+ ξ 3 )( ξ 1 ξ 2 + ξ 3 2)/8 N 6 =(1+ ξ 1 )(1 ξ 2 )(1+ ξ 3 )(+ ξ 1 ξ 2 + ξ 3 2)/8 N 7 =(1+ ξ 1 )(1+ ξ 2 )(1+ ξ 3 )(+ ξ 1 + ξ 2 + ξ 3 2)/8 N 8 =(1 ξ 1 )(1+ ξ 2 )(1+ ξ 3 )( ξ 1 + ξ 2 + ξ 3 2)/8 N 9 =(1 ξ 1 2 )(1 ξ 2 )(1 ξ 3 )/4 N 10 =(1+ ξ 1 )(1 ξ 2 2 )(1 ξ 3 )/4 N 11 =(1 ξ 1 2 )(1+ ξ 2 )(1 ξ 3 )/4 N 12 =(1 ξ 1 )(1 ξ 2 2 )(1 ξ 3 )/4 N 13 =(1 ξ 1 2 )(1 ξ 2 )(1+ ξ 3 )/4 N 14 =(1+ ξ 1 )(1 ξ 2 2 )(1+ ξ 3 )/4 N 15 =(1 ξ 1 2 )(1+ ξ 2 )(1+ ξ 3 )/4 N 16 =(1 ξ 1 )(1 ξ 2 2 )(1+ ξ 3 )/4 N 17 =(1 ξ 1 )(1 ξ 2 )(1 ξ 3 2 )/4 N 18 =(1+ ξ 1 )(1 ξ 2 )(1 ξ 3 2 )/4 N 19 =(1+ ξ 1 )(1+ ξ 2 )(1 ξ 3 2 )/4 N 20 =(1 ξ 1 )(1+ ξ 2 )(1 ξ 3 2 )/4 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGobWaaWbaaSqabeaacaaIXa aaaOGaeyypa0JaaiikaiaaigdacqGHsislcqaH+oaEdaWgaaWcbaGa aGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+oaEdaWgaa WcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+oaE daWgaaWcbaGaaG4maaqabaGccaGGPaGaaiikaiabgkHiTiabe67a4n aaBaaaleaacaaIXaaabeaakiabgkHiTiabe67a4naaBaaaleaacaaI YaaabeaakiabgkHiTiabe67a4naaBaaaleaacaaIZaaabeaakiabgk HiTiaaikdacaGGPaGaai4laiaaiIdaaeaacaWGobWaaWbaaSqabeaa caaIYaaaaOGaeyypa0JaaiikaiaaigdacqGHRaWkcqaH+oaEdaWgaa WcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+oaE daWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcq aH+oaEdaWgaaWcbaGaaG4maaqabaGccaGGPaGaaiikaiabe67a4naa BaaaleaacaaIXaaabeaakiabgkHiTiabe67a4naaBaaaleaacaaIYa aabeaakiabgkHiTiabe67a4naaBaaaleaacaaIZaaabeaakiabgkHi TiaaikdacaGGPaGaai4laiaaiIdaaeaacaWGobWaaWbaaSqabeaaca aIZaaaaOGaeyypa0JaaiikaiaaigdacqGHRaWkcqaH+oaEdaWgaaWc baGaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHRaWkcqaH+oaEda WgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH +oaEdaWgaaWcbaGaaG4maaqabaGccaGGPaGaaiikaiabe67a4naaBa aaleaacaaIXaaabeaakiabgUcaRiabe67a4naaBaaaleaacaaIYaaa beaakiabgkHiTiabe67a4naaBaaaleaacaaIZaaabeaakiabgkHiTi aaikdacaGGPaGaai4laiaaiIdaaeaacaWGobWaaWbaaSqabeaacaaI 0aaaaOGaeyypa0JaaiikaiaaigdacqGHsislcqaH+oaEdaWgaaWcba GaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHRaWkcqaH+oaEdaWg aaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+o aEdaWgaaWcbaGaaG4maaqabaGccaGGPaGaaiikaiabgkHiTiabe67a 4naaBaaaleaacaaIXaaabeaakiabgUcaRiabe67a4naaBaaaleaaca aIYaaabeaakiabgkHiTiabe67a4naaBaaaleaacaaIZaaabeaakiab gkHiTiaaikdacaGGPaGaai4laiaaiIdaaeaacaWGobWaaWbaaSqabe aacaaI1aaaaOGaeyypa0JaaiikaiaaigdacqGHsislcqaH+oaEdaWg aaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcqaH+o aEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGHRaWk cqaH+oaEdaWgaaWcbaGaaG4maaqabaGccaGGPaGaaiikaiabgkHiTi abe67a4naaBaaaleaacaaIXaaabeaakiabgkHiTiabe67a4naaBaaa leaacaaIYaaabeaakiabgUcaRiabe67a4naaBaaaleaacaaIZaaabe aakiabgkHiTiaaikdacaGGPaGaai4laiaaiIdaaeaacaWGobWaaWba aSqabeaacaaI2aaaaOGaeyypa0JaaiikaiaaigdacqGHRaWkcqaH+o aEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsisl cqaH+oaEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacq GHRaWkcqaH+oaEdaWgaaWcbaGaaG4maaqabaGccaGGPaGaaiikaiab gUcaRiabe67a4naaBaaaleaacaaIXaaabeaakiabgkHiTiabe67a4n aaBaaaleaacaaIYaaabeaakiabgUcaRiabe67a4naaBaaaleaacaaI ZaaabeaakiabgkHiTiaaikdacaGGPaGaai4laiaaiIdaaeaacaWGob WaaWbaaSqabeaacaaI3aaaaOGaeyypa0JaaiikaiaaigdacqGHRaWk cqaH+oaEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacq GHRaWkcqaH+oaEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaa igdacqGHRaWkcqaH+oaEdaWgaaWcbaGaaG4maaqabaGccaGGPaGaai ikaiabgUcaRiabe67a4naaBaaaleaacaaIXaaabeaakiabgUcaRiab e67a4naaBaaaleaacaaIYaaabeaakiabgUcaRiabe67a4naaBaaale aacaaIZaaabeaakiabgkHiTiaaikdacaGGPaGaai4laiaaiIdaaeaa caWGobWaaWbaaSqabeaacaaI4aaaaOGaeyypa0Jaaiikaiaaigdacq GHsislcqaH+oaEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaa igdacqGHRaWkcqaH+oaEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaai ikaiaaigdacqGHRaWkcqaH+oaEdaWgaaWcbaGaaG4maaqabaGccaGG PaGaaiikaiabgkHiTiabe67a4naaBaaaleaacaaIXaaabeaakiabgU caRiabe67a4naaBaaaleaacaaIYaaabeaakiabgUcaRiabe67a4naa BaaaleaacaaIZaaabeaakiabgkHiTiaaikdacaGGPaGaai4laiaaiI daaeaacaWGobWaaWbaaSqabeaacaaI5aaaaOGaeyypa0Jaaiikaiaa igdacqGHsislcqaH+oaEdaqhaaWcbaGaaGymaaqaaiaaikdaaaGcca GGPaGaaiikaiaaigdacqGHsislcqaH+oaEdaWgaaWcbaGaaGOmaaqa baGccaGGPaGaaiikaiaaigdacqGHsislcqaH+oaEdaWgaaWcbaGaaG 4maaqabaGccaGGPaGaai4laiaaisdacaaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caWGobWaaWbaaSqabeaacaaIXaGaaGimaaaaki abg2da9iaacIcacaaIXaGaey4kaSIaeqOVdG3aaSbaaSqaaiaaigda aeqaaOGaaiykaiaacIcacaaIXaGaeyOeI0IaeqOVdG3aa0baaSqaai aaikdaaeaacaaIYaaaaOGaaiykaiaacIcacaaIXaGaeyOeI0IaeqOV dG3aaSbaaSqaaiaaiodaaeqaaOGaaiykaiaac+cacaaI0aaabaGaam OtamaaCaaaleqabaGaaGymaiaaigdaaaGccqGH9aqpcaGGOaGaaGym aiabgkHiTiabe67a4naaDaaaleaacaaIXaaabaGaaGOmaaaakiaacM cacaGGOaGaaGymaiabgUcaRiabe67a4naaBaaaleaacaaIYaaabeaa kiaacMcacaGGOaGaaGymaiabgkHiTiabe67a4naaBaaaleaacaaIZa aabeaakiaacMcacaGGVaGaaGinaiaaykW7caaMc8UaaGPaVlaad6ea daahaaWcbeqaaiaaigdacaaIYaaaaOGaeyypa0Jaaiikaiaaigdacq GHsislcqaH+oaEdaWgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaa igdacqGHsislcqaH+oaEdaqhaaWcbaGaaGOmaaqaaiaaikdaaaGcca GGPaGaaiikaiaaigdacqGHsislcqaH+oaEdaWgaaWcbaGaaG4maaqa baGccaGGPaGaai4laiaaisdaaeaacaWGobWaaWbaaSqabeaacaaIXa GaaG4maaaakiabg2da9iaacIcacaaIXaGaeyOeI0IaeqOVdG3aa0ba aSqaaiaaigdaaeaacaaIYaaaaOGaaiykaiaacIcacaaIXaGaeyOeI0 IaeqOVdG3aaSbaaSqaaiaaikdaaeqaaOGaaiykaiaacIcacaaIXaGa ey4kaSIaeqOVdG3aaSbaaSqaaiaaiodaaeqaaOGaaiykaiaac+caca aI0aGaaGPaVlaaykW7caaMc8UaamOtamaaCaaaleqabaGaaGymaiaa isdaaaGccqGH9aqpcaGGOaGaaGymaiabgUcaRiabe67a4naaBaaale aacaaIXaaabeaakiaacMcacaGGOaGaaGymaiabgkHiTiabe67a4naa DaaaleaacaaIYaaabaGaaGOmaaaakiaacMcacaGGOaGaaGymaiabgU caRiabe67a4naaBaaaleaacaaIZaaabeaakiaacMcacaGGVaGaaGin aaqaaiaad6eadaahaaWcbeqaaiaaigdacaaI1aaaaOGaeyypa0Jaai ikaiaaigdacqGHsislcqaH+oaEdaqhaaWcbaGaaGymaaqaaiaaikda aaGccaGGPaGaaiikaiaaigdacqGHRaWkcqaH+oaEdaWgaaWcbaGaaG OmaaqabaGccaGGPaGaaiikaiaaigdacqGHRaWkcqaH+oaEdaWgaaWc baGaaG4maaqabaGccaGGPaGaai4laiaaisdacaaMc8UaaGPaVlaayk W7caWGobWaaWbaaSqabeaacaaIXaGaaGOnaaaakiabg2da9iaacIca caaIXaGaeyOeI0IaeqOVdG3aaSbaaSqaaiaaigdaaeqaaOGaaiykai aacIcacaaIXaGaeyOeI0IaeqOVdG3aa0baaSqaaiaaikdaaeaacaaI YaaaaOGaaiykaiaacIcacaaIXaGaey4kaSIaeqOVdG3aaSbaaSqaai aaiodaaeqaaOGaaiykaiaac+cacaaI0aaabaGaamOtamaaCaaaleqa baGaaGymaiaaiEdaaaGccqGH9aqpcaGGOaGaaGymaiabgkHiTiabe6 7a4naaBaaaleaacaaIXaaabeaakiaacMcacaGGOaGaaGymaiabgkHi Tiabe67a4naaBaaaleaacaaIYaaabeaakiaacMcacaGGOaGaaGymai abgkHiTiabe67a4naaDaaaleaacaaIZaaabaGaaGOmaaaakiaacMca caGGVaGaaGinaiaaykW7caaMc8UaaGPaVlaad6eadaahaaWcbeqaai aaigdacaaI4aaaaOGaeyypa0JaaiikaiaaigdacqGHRaWkcqaH+oaE daWgaaWcbaGaaGymaaqabaGccaGGPaGaaiikaiaaigdacqGHsislcq aH+oaEdaWgaaWcbaGaaGOmaaqabaGccaGGPaGaaiikaiaaigdacqGH sislcqaH+oaEdaqhaaWcbaGaaG4maaqaaiaaikdaaaGccaGGPaGaai 4laiaaisdaaeaacaWGobWaaWbaaSqabeaacaaIXaGaaGyoaaaakiab g2da9iaacIcacaaIXaGaey4kaSIaeqOVdG3aaSbaaSqaaiaaigdaae qaaOGaaiykaiaacIcacaaIXaGaey4kaSIaeqOVdG3aaSbaaSqaaiaa ikdaaeqaaOGaaiykaiaacIcacaaIXaGaeyOeI0IaeqOVdG3aa0baaS qaaiaaiodaaeaacaaIYaaaaOGaaiykaiaaykW7caGGVaGaaGinaiaa ykW7caaMc8UaamOtamaaCaaaleqabaGaaGOmaiaaicdaaaGccqGH9a qpcaGGOaGaaGymaiabgkHiTiabe67a4naaBaaaleaacaaIXaaabeaa kiaacMcacaGGOaGaaGymaiabgUcaRiabe67a4naaBaaaleaacaaIYa aabeaakiaacMcacaGGOaGaaGymaiabgkHiTiabe67a4naaDaaaleaa caaIZaaabaGaaGOmaaaakiaacMcacaGGVaGaaGinaaaaaa@79DD@

 

    

The element faces are numbered as follows.

 

Linear and quadratic tetrahedral

 

Face 1 has nodes 1,2,3

Face 2 has nodes 1,4,2
Face 3 has nodes 2,4,3
Face 4 has nodes 3,4,1

 

Linear and quadratic brick elements

 

Face 1 has nodes 1,2,3,4

Face 2 has nodes 5,8,7,6

Face 3 has nodes 1,5,6,3

Face 4 has nodes 2,6,7,3

Face 5 has nodes 3,7,8,4

Face 6 has nodes 4,8,5,1

 



 

8.1.11 Volume integrals for stiffness and force in terms of normalized coordinates

 

In this section we outline the procedure that is used to re-write the integrals for the element stiffness and force in terms of the normalized coordinates ξ i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aaSbaaSqaaiaadMgaaeqaaa aa@33BA@ .  The integrals are

k aibk (l) = V e (l) C ijkl N a (x) x j N b (x) x l dV f i a (l) = V e (l) b i N a (x) dV+ 2 V e (l) t i * N a (x) dA 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 PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaadA gadaqhaaWcbaGaamyAaaqaaiaadggaaaGcdaahaaWcbeqaaiaacIca caWGSbGaaiykaaaakiabg2da9maapefabaGaamOyamaaBaaaleaaca WGPbaabeaakiaad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGaaCiE aiaacMcaaSqaaiaadAfadaqhaaadbaGaamyzaaqaaiaacIcacaWGSb GaaiykaaaaaSqab0Gaey4kIipakiaadsgacaWGwbGaey4kaSYaa8qu aeaacaWG0bWaa0baaSqaaiaadMgaaeaacaGGQaaaaOGaamOtamaaCa aaleqabaGaamyyaaaakiaacIcacaWH4bGaaiykaaWcbaGaeyOaIy7a aSbaaWqaaiaaikdaaeqaaSGaamOvamaaDaaameaacaWGLbaabaGaai ikaiaadYgacaGGPaaaaaWcbeqdcqGHRiI8aOGaamizaiaadgeacaaM c8oaaa@9C54@

To evaluate them, we need to

1.      Find a way to calculate the derivatives of the shape functions in terms of ξ i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aaSbaaSqaaiaadMgaaeqaaa aa@33AA@

2.      Map the volume (or area) integral to the region 1 ξ i +1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkHiTiaaigdacqGHKjYOcqaH+oaEda WgaaWcbaGaamyAaaqabaGccqGHKjYOcqGHRaWkcaaIXaaaaa@3A1E@

 

Calculating the shape function derivatives. The shape function derivatives can be evaluated by writing

N a x j = N a ξ i ξ i x j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaSaaaeaacqGHciITcaWGobWaaWbaaS qabeaacaWGHbaaaaGcbaGaeyOaIyRaamiEamaaBaaaleaacaWGQbaa beaaaaGccqGH9aqpdaWcaaqaaiabgkGi2kaad6eadaahaaWcbeqaai aadggaaaaakeaacqGHciITcqaH+oaEdaWgaaWcbaGaamyAaaqabaaa aOWaaSaaaeaacqGHciITcqaH+oaEdaWgaaWcbaGaamyAaaqabaaake aacqGHciITcaWG4bWaaSbaaSqaaiaadQgaaeqaaaaaaaa@4851@

where the derivatives N a / ξ i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkGi2kaad6eadaahaaWcbeqaaiaadg gaaaGccaGGVaGaeyOaIyRaeqOVdG3aaSbaaSqaaiaadMgaaeqaaaaa @38D4@  are easy to compute (just differentiate the expressions given earlier…).  To compute ξ i / x j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkGi2kabe67a4naaBaaaleaacaWGPb aabeaakiaac+cacqGHciITcaWG4bWaaSbaaSqaaiaadQgaaeqaaaaa @3906@  recall that the coordinates of a point at position ξ j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabe67a4naaBaaaleaacaWGQbaabeaaaa a@3366@  within an element can be determined as

x i = a=1 N e N a ( ξ j ) x i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadIhadaWgaaWcbaGaamyAaaqabaGccq GH9aqpdaaeWbqaaiaad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGa eqOVdG3aaSbaaSqaaiaadQgaaeqaaOGaaiykaiaadIhadaqhaaWcba GaamyAaaqaaiaadggaaaaabaGaamyyaiabg2da9iaaigdaaeaacaWG obWaaSbaaWqaaiaadwgaaeqaaaqdcqGHris5aaaa@43A6@

where N e MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaad6eadaWgaaWcbaGaamyzaaqabaaaaa@3271@  denotes the number of nodes on the element. Therefore

x i ξ j = ξ j ( a=1 N e N a (ξ) x i a )= 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 aeqaaaaakiabg2da9maalaaabaGaeyOaIylabaGaeyOaIyRaeqOVdG 3aaSbaaSqaaiaadQgaaeqaaaaakmaabmaabaWaaabCaeaacaWGobWa aWbaaSqabeaacaWGHbaaaOGaaiikaiabe67a4jaacMcacaWG4bWaa0 baaSqaaiaadMgaaeaacaWGHbaaaaqaaiaadggacqGH9aqpcaaIXaaa baGaamOtamaaBaaameaacaWGLbaabeaaa0GaeyyeIuoaaOGaayjkai aawMcaaiabg2da9maaqahabaWaaSaaaeaacqGHciITcaWGobWaaWba aSqabeaacaWGHbaaaOGaaiikaiabe67a4jaacMcaaeaacqGHciITcq aH+oaEdaWgaaWcbaGaamOAaaqabaaaaOGaamiEamaaDaaaleaacaWG PbaabaGaamyyaaaaaeaacaWGHbGaeyypa0JaaGymaaqaaiaad6eada WgaaadbaGaamyzaaqabaaaniabggHiLdaaaa@657F@

Note that x i / ξ j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkGi2kaadIhadaWgaaWcbaGaamyAaa qabaGccaGGVaGaeyOaIyRaeqOVdG3aaSbaaSqaaiaadQgaaeqaaaaa @3906@  is a 2x2 matrix (in 2D) or a 3x3 matrix (in 3D).  Finally, ξ i / x j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkGi2kabe67a4naaBaaaleaacaWGPb aabeaakiaac+cacqGHciITcaWG4bWaaSbaaSqaaiaadQgaaeqaaaaa @3906@  follows as the inverse of this matrix

ξ 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@

 

Mapping the volume integral: To map the region of integration we define

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@

where the matrix x i / ξ j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabgkGi2kaadIhadaWgaaWcbaGaamyAaa qabaGccaGGVaGaeyOaIyRaeqOVdG3aaSbaaSqaaiaadQgaaeqaaaaa @3906@  was defined earlier.  Then the integral with respect to x is mapped into an integral with respect to ξ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8skY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaah67aaaa@31E2@  by setting

k aibk = Ω C ijkl N a (x) x j N b (x) x l Jd V ξ f i a = Ω b i N a (x) Jd V ξ + Ω t i * N a (x) J d A ξ ^ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaam4AamaaBaaaleaacaWGHbGaamyAai aadkgacaWGRbaabeaakiabg2da9maapefabaGaam4qamaaBaaaleaa caWGPbGaamOAaiaadUgacaWGSbaabeaakmaalaaabaGaeyOaIyRaam OtamaaCaaaleqabaGaamyyaaaakiaacIcacaWH4bGaaiykaaqaaiab gkGi2kaadIhadaWgaaWcbaGaamOAaaqabaaaaOWaaSaaaeaacqGHci ITcaWGobWaaWbaaSqabeaacaWGIbaaaOGaaiikaiaahIhacaGGPaaa baGaeyOaIyRaamiEamaaBaaaleaacaWGSbaabeaaaaGccaWGkbGaam izaiaadAfadaWgaaWcbaGaaCOVdaqabaaabaGaeuyQdCfabeqdcqGH RiI8aOGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caWGMbWaa0baaSqaaiaadMgaaeaacaWGHbaaaOGaey ypa0Zaa8quaeaacaWGIbWaaSbaaSqaaiaadMgaaeqaaOGaamOtamaa CaaaleqabaGaamyyaaaakiaacIcacaWH4bGaaiykaaWcbaGaeuyQdC fabeqdcqGHRiI8aOGaamOsaiaadsgacaWGwbWaaSbaaSqaaiaah67a aeqaaOGaey4kaSYaa8quaeaacaWG0bWaa0baaSqaaiaadMgaaeaaca GGQaaaaOGaamOtamaaCaaaleqabaGaamyyaaaakiaacIcacaWH4bGa aiykaaWcbaGaeyOaIyRaeuyQdCfabeqdcqGHRiI8aOGabmOsayaata GaamizaiaadgeadaWgaaWcbaGabCOVdyaajaaabeaakiaaykW7aaa@9548@

We note in passing that the boundary integral in the element force vector can be regarded as a 1-D line integral for 2D elements and a 2D surface integral for 3D elements.  So the procedures we developed in 8.1.5 (1D elements) can be used to evaluate the surface integral for a 2D element.  Similarly, the procedures we develop to integrate stiffness matrices for 2D elements can be used to evaluate the surface integral for a 3D element.

 

 

 

8.1.12 Numerical integration schemes for 2D and 3D elements

 

Finally, to evaluate the integrals, we once again adopt a quadrature scheme, so that

Ω f( ξ i )d V ξ = I=1 N I w I f( ξ i I ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaeaacaWGMbGaaiikaiabe67a4n aaBaaaleaacaWGPbaabeaakiaacMcacaWGKbGaamOvamaaBaaaleaa cqaH+oaEaeqaaaqaaiabfM6axbqab0Gaey4kIipakiabg2da9maaqa habaGaam4DamaaBaaaleaacaWGjbaabeaakiaadAgacaGGOaGaeqOV dG3aa0baaSqaaiaadMgaaeaacaWGjbaaaOGaaiykaaWcbaGaamysai abg2da9iaaigdaaeaacaWGobWaaSbaaWqaaiaadMeaaeqaaaqdcqGH ris5aaaa@4D23@

The integration points ξ j I MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiabe67a4naaDaaaleaacaWGQbaabaGaam ysaaaaaaa@3435@  and weights w I MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaadEhadaWgaaWcbaGaamysaaqabaaaaa@327E@  depend on the element geometry, and are listed below for a few common element types

 

Integration points for triangular elements

 

1 point 

ξ 1 1 =1/3 ξ 2 1 =1/3 w 1 =1/2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aa0baaSqaaiaaigdaaeaaca aIXaaaaOGaeyypa0JaaGymaiaac+cacaaIZaGaaGzaVlaaygW7caaM b8UaaGzaVlaaygW7caaMb8UaaGzaVlaaygW7caaMb8UaaGzaVlaayg W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7cqaH+oaEdaqhaaWcbaGaaGOmaaqaaiaaigdaaaGccqGH9a qpcaaIXaGaai4laiaaiodacaaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWG3b WaaSbaaSqaaiaaigdaaeqaaOGaeyypa0JaaGymaiaac+cacaaIYaaa aa@7484@

3 point 

ξ 1 1 =1/2 ξ 2 1 =0 w 1 =1/6 ξ 1 2 =0 ξ 2 2 =1/2 w 2 =1/6 ξ 1 3 =1/2 ξ 2 3 =1/2 w 3 =1/6 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacqaH+oaEdaqhaaWcbaGaaGymaa qaaiaaigdaaaGccqGH9aqpcaaIXaGaai4laiaaikdacaaMb8UaaGza VlaaygW7caaMb8UaaGzaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7cqaH+oaEdaqhaaWcbaGaaGOmaaqaaiaaigdaaaGc cqGH9aqpcaaIWaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaadEhadaWgaaWcbaGaaGymaaqa baGccqGH9aqpcaaIXaGaai4laiaaiAdaaeaacqaH+oaEdaqhaaWcba GaaGymaaqaaiaaikdaaaGccqGH9aqpcaaIWaGaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaeqOVdG3aa0baaSqaaiaaikdaaeaacaaIYaaaaOGaeyyp a0JaaGymaiaac+cacaaIYaGaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caWG3bWaaSbaaSqaaiaaikdaaeqaaOGaeyypa0JaaGymaiaac+ca caaI2aaabaGaeqOVdG3aa0baaSqaaiaaigdaaeaacaaIZaaaaOGaey ypa0JaaGymaiaac+cacaaIYaGaaGzaVlaaygW7caaMb8UaaGzaVlaa ygW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Uaeq OVdG3aa0baaSqaaiaaikdaaeaacaaIZaaaaOGaeyypa0JaaGymaiaa c+cacaaIYaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWG3bWaaS baaSqaaiaaiodaaeqaaOGaeyypa0JaaGymaiaac+cacaaI2aaaaaa@BB84@  or ξ 1 1 =0.6 ξ 2 1 =0.2 w 1 =1/6 ξ 1 2 =0.2 ξ 2 2 =0.6 w 2 =1/6 ξ 1 3 =0.2 ξ 2 3 =0.2 w 3 =1/6 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacqaH+oaEdaqhaaWcbaGaaGymaa qaaiaaigdaaaGccqGH9aqpcaaIWaGaaiOlaiaaiAdacaaMb8UaaGza VlaaygW7caaMb8UaaGzaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlabe67a4naaDaaaleaacaaIYaaabaGaaGymaaaakiabg2da 9iaaicdacaGGUaGaaGOmaiaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8Uaam4DamaaBaaaleaacaaIXaaabeaakiab g2da9iaaigdacaGGVaGaaGOnaaqaaiabe67a4naaDaaaleaacaaIXa aabaGaaGOmaaaakiabg2da9iaaicdacaGGUaGaaGOmaiaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlabe67a4naaDaaaleaacaaIYa aabaGaaGOmaaaakiabg2da9iaaicdacaGGUaGaaGOnaiaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWG3bWaaSbaaSqaai aaikdaaeqaaOGaeyypa0JaaGymaiaac+cacaaI2aaabaGaeqOVdG3a a0baaSqaaiaaigdaaeaacaaIZaaaaOGaeyypa0JaaGimaiaac6caca aIYaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaeqOVdG3a a0baaSqaaiaaikdaaeaacaaIZaaaaOGaeyypa0JaaGimaiaac6caca aIYaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa dEhadaWgaaWcbaGaaG4maaqabaGccqGH9aqpcaaIXaGaai4laiaaiA daaaaa@AD6C@

(the first scheme here is optimal, but has some disadvantages for quadratic elements because the integration points coincide with the midside nodes.  The second scheme is less accurate but more robust).

4 point

ξ 1 1 =0.6 ξ 2 1 =0.2 w 1 =25/96 ξ 1 2 =0.2 ξ 2 2 =0.6 w 2 =25/96 ξ 1 3 =0.2 ξ 2 3 =0.2 w 3 =25/96 ξ 1 4 =1/3 ξ 2 4 =1/3 w 4 =27/96 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacqaH+oaEdaqhaaWcbaGaaGymaa qaaiaaigdaaaGccqGH9aqpcaaIWaGaaiOlaiaaiAdacaaMb8UaaGza VlaaygW7caaMb8UaaGzaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlabe67a4naaDaaaleaacaaIYaaa baGaaGymaaaakiabg2da9iaaicdacaGGUaGaaGOmaiaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaadEhadaWgaaWcbaGaaGymaaqabaGccqGH9aqpca aIYaGaaGynaiaac+cacaaI5aGaaGOnaaqaaiabe67a4naaDaaaleaa caaIXaaabaGaaGOmaaaakiabg2da9iaaicdacaGGUaGaaGOmaiaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa Vlabe67a4naaDaaaleaacaaIYaaabaGaaGOmaaaakiabg2da9iaaic dacaGGUaGaaGOnaiaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Uaam4DamaaBaaaleaaca aIYaaabeaakiabg2da9iaaikdacaaI1aGaai4laiaaiMdacaaI2aaa baGaeqOVdG3aa0baaSqaaiaaigdaaeaacaaIZaaaaOGaeyypa0JaaG imaiaac6cacaaIYaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaeqOVdG3aa0baaSqaaiaaikdaaeaaca aIZaaaaOGaeyypa0JaaGimaiaac6cacaaIYaGaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caWG3bWaaSbaaSqaaiaaiodaaeqaaOGaeyypa0JaaGOmaiaaiwda caGGVaGaaGyoaiaaiAdaaeaacqaH+oaEdaqhaaWcbaGaaGymaaqaai aaisdaaaGccqGH9aqpcaaIXaGaai4laiaaiodacaaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7cqaH+oaEda qhaaWcbaGaaGOmaaqaaiaaisdaaaGccqGH9aqpcaaIXaGaai4laiaa iodacaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8Uaam4DamaaBaaaleaacaaI0aaabeaakiabg2da 9iabgkHiTiaaikdacaaI3aGaai4laiaaiMdacaaI2aaaaaa@048D@

 

 Integration points for tetrahedral elements

 

1 point

 

ξ 1 1 =1/4 ξ 2 1 =1/4 ξ 3 1 =1/4 w 1 =1/6 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aa0baaSqaaiaaigdaaeaaca aIXaaaaOGaeyypa0JaaGymaiaac+cacaaI0aGaaGzaVlaaygW7caaM b8UaaGzaVlaaygW7caaMb8UaaGzaVlaaygW7caaMb8UaaGzaVlaayg W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7cqaH+oaEdaqhaaWcbaGaaGOmaaqaaiaaigdaaaGccqGH9a qpcaaIXaGaai4laiaaisdacaaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7cqaH+o aEdaqhaaWcbaGaaG4maaqaaiaaigdaaaGccqGH9aqpcaaIXaGaai4l aiaaisdacaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWG3bWaaSbaaSqaaiaa igdaaeqaaOGaeyypa0JaaGymaiaac+cacaaI2aaaaa@8DB2@

4 point

 

ξ 1 1 =α ξ 2 1 =β ξ 3 1 =β w 1 =1/24 ξ 1 2 =β ξ 2 1 =α ξ 3 2 =β w 2 =1/24 ξ 1 3 =β ξ 2 3 =β ξ 3 3 =α w 3 =1/24 ξ 1 4 =β ξ 2 4 =β ξ 3 4 =β w 4 =1/24 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacqaH+oaEdaqhaaWcbaGaaGymaa qaaiaaigdaaaGccqGH9aqpcqaHXoqycaaMb8UaaGzaVlaaygW7caaM b8UaaGzaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7cqaH+oaEdaqh aaWcbaGaaGOmaaqaaiaaigdaaaGccqGH9aqpcqaHYoGycaaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7cqaH+oaEdaqhaaWcbaGaaG4maaqaaiaaigdaaa GccqGH9aqpcqaHYoGycaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWG3bWaaS baaSqaaiaaigdaaeqaaOGaeyypa0JaaGymaiaac+cacaaIYaGaaGin aaqaaiabe67a4naaDaaaleaacaaIXaaabaGaaGOmaaaakiabg2da9i abek7aIjaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaeqOVdG 3aa0baaSqaaiaaikdaaeaacaaIXaaaaOGaeyypa0JaeqySdeMaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlabe67a4naaDaaaleaacaaIZaaabaGaaGOmaaaakiabg2da 9iabek7aIjaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaadEhadaWgaaWcbaGa aGOmaaqabaGccqGH9aqpcaaIXaGaai4laiaaikdacaaI0aaabaGaeq OVdG3aa0baaSqaaiaaigdaaeaacaaIZaaaaOGaeyypa0JaeqOSdiMa aGzaVlaaygW7caaMb8UaaGzaVlaaygW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlabe67a4naaDaaaleaacaaIYaaabaGaaG4maa aakiabg2da9iabek7aIjaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaeqOVdG3aa0baaS qaaiaaiodaaeaacaaIZaaaaOGaeyypa0JaeqySdeMaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8Uaam4DamaaBaaaleaacaaIZaaabeaakiabg2da9iaa igdacaGGVaGaaGOmaiaaisdaaeaacqaH+oaEdaqhaaWcbaGaaGymaa qaaiaaisdaaaGccqGH9aqpcqaHYoGycaaMc8UaaGPaVlaaykW7caaM c8UaaGzaVlaaygW7caaMb8UaaGzaVlaaygW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaeqOV dG3aa0baaSqaaiaaikdaaeaacaaI0aaaaOGaeyypa0JaeqOSdiMaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7cqaH+oaEdaqhaaWcbaGaaG4maaqaaiaaisdaaa GccqGH9aqpcqaHYoGycaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caWG3bWaaS baaSqaaiaaisdaaeqaaOGaeyypa0JaaGymaiaac+cacaaIYaGaaGin aaaaaa@8B5C@  

where α=0.58541020,β=0.13819660 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqySdeMaeyypa0JaaGimaiaac6caca aI1aGaaGioaiaaiwdacaaI0aGaaGymaiaaicdacaaIYaGaaGimaiaa cYcacaaMc8UaaGPaVlaaykW7caaMc8UaeqOSdiMaeyypa0JaaGimai aac6cacaaIXaGaaG4maiaaiIdacaaIXaGaaGyoaiaaiAdacaaI2aGa aGimaaaa@4BAA@

 

 

Quadrilateral and hexahedral elements

 

For quadrilateral elements we can simply regard the integral over 2 spatial dimensions as successive 1-D integrals

1 +1 1 +1 f( ξ 1 , ξ 2 )d ξ 1 d ξ 2 = I=1 N J=1 N w I w J f( ξ 1 I , ξ 2 J ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8qCaeaadaWdXbqaaiaadAgadaqada qaaiabe67a4naaBaaaleaacaaIXaaabeaakiaacYcacqaH+oaEdaWg aaWcbaGaaGOmaaqabaaakiaawIcacaGLPaaacaWGKbGaeqOVdG3aaS baaSqaaiaaigdaaeqaaOGaamizaiabe67a4naaBaaaleaacaaIYaaa beaaaeaacqGHsislcaaIXaaabaGaey4kaSIaaGymaaqdcqGHRiI8aO Gaeyypa0ZaaabCaeaadaaeWbqaaiaadEhadaWgaaWcbaGaamysaaqa baGccaWG3bWaaSbaaSqaaiaadQeaaeqaaOGaamOzamaabmaabaGaeq OVdG3aa0baaSqaaiaaigdaaeaacaWGjbaaaOGaaiilaiabe67a4naa DaaaleaacaaIYaaabaGaamOsaaaaaOGaayjkaiaawMcaaaWcbaGaam Osaiabg2da9iaaigdaaeaacaWGobaaniabggHiLdaaleaacaWGjbGa eyypa0JaaGymaaqaaiaad6eaa0GaeyyeIuoaaSqaaiabgkHiTiaaig daaeaacqGHRaWkcaaIXaaaniabgUIiYdaaaa@6664@

which gives rise to the following 2D quadrature scheme:  Let η I MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4TdG2aaSbaaSqaaiaadMeaaeqaaa aa@3375@  and v I MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamODamaaBaaaleaacaWGjbaabeaaaa a@32C4@  for I=1…M denote 1-D quadrature points and weights listed below. Then in 2D, an N=M×M MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOtaiabg2da9iaad2eacqGHxdaTca WGnbaaaa@3663@  quadrature scheme can be generated as follows:

for J=1…M and K=1…M    let   ξ 1 I = η J ξ 2 I = η K w I = v J v K ,I=M(K1)+J MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aa0baaSqaaiaaigdaaeaaca WGjbaaaOGaeyypa0Jaeq4TdG2aaSbaaSqaaiaadQeaaeqaaOGaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlabe67a4naaDaaaleaacaaIYaaabaGaamysaaaakiabg2da 9iabeE7aOnaaBaaaleaacaWGlbaabeaakiaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaadEhadaWgaaWcbaGaamysaaqabaGccqGH9aqpcaWG2bWaaS baaSqaaiaadQeaaeqaaOGaamODamaaBaaaleaacaWGlbaabeaakiaa cYcacaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caWGjbGaeyypa0JaamytaiaacIcacaWGlbGaeyOeI0Ia aGymaiaacMcacqGHRaWkcaWGkbaaaa@7EF9@

Similarly, in 3D, we generate an N=M×M×M MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOtaiabg2da9iaad2eacqGHxdaTca WGnbGaey41aqRaamytaaaa@394C@  scheme as:

for J=1…M , K=1…M  L=1…M   let

ξ 1 I = η J ξ 2 I = η K ξ 3 I = η L w I = v J v K v L ,I= M 2 (L1)+M(K1)+J MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aa0baaSqaaiaaigdaaeaaca WGjbaaaOGaeyypa0Jaeq4TdG2aaSbaaSqaaiaadQeaaeqaaOGaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaeqOVdG3aa0baaSqaai aaikdaaeaacaWGjbaaaOGaeyypa0Jaeq4TdG2aaSbaaSqaaiaadUea aeqaaOGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl abe67a4naaDaaaleaacaaIZaaabaGaamysaaaakiabg2da9iabeE7a OnaaBaaaleaacaWGmbaabeaakiaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caWG3bWaaSbaaSqaaiaadMeaaeqaaOGaeyyp a0JaamODamaaBaaaleaacaWGkbaabeaakiaadAhadaWgaaWcbaGaam 4saaqabaGccaWG2bWaaSbaaSqaaiaadYeaaeqaaOGaaiilaiaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaadMeacqGH9aqpcaWGnb WaaWbaaSqabeaacaaIYaaaaOGaaiikaiaadYeacqGHsislcaaIXaGa aiykaiabgUcaRiaad2eacaGGOaGaam4saiabgkHiTiaaigdacaGGPa Gaey4kaSIaamOsaaaa@86FB@

M=1

 

η 1 =0 v 1 =2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4TdG2aaSbaaSqaaiaaigdaaeqaaO Gaeyypa0JaaGimaiaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaadAhadaWgaa WcbaGaaGymaaqabaGccqGH9aqpcaaIYaaaaa@4B5E@

M=2

 

η 1 =0.5773502691 v 1 =1.0 η 2 =0.5773502691 v 2 =1.0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacqaH3oaAdaWgaaWcbaGaaGymaa qabaGccqGH9aqpcqGHsislcaaIWaGaaiOlaiaaiwdacaaI3aGaaG4n aiaaiodacaaI1aGaaGPaVlaaicdacaaIYaGaaGOnaiaaiMdacaaIXa GaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaadAhadaWgaaWcbaGaaGymaaqabaGccqGH 9aqpcaaIXaGaaiOlaiaaicdaaeaacqaH3oaAdaWgaaWcbaGaaGOmaa qabaGccqGH9aqpcaaIWaGaaiOlaiaaiwdacaaI3aGaaG4naiaaioda caaI1aGaaGPaVlaaicdacaaIYaGaaGOnaiaaiMdacaaIXaGaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaadAhadaWgaaWcbaGaaGOm aaqabaGccqGH9aqpcaaIXaGaaiOlaiaaicdaaaaa@9750@

M=3

 

η 1 =0.7745966692 v 1 =0.55555555555 η 2 =0. v 2 =0.88888888888 η 3 =0.7745966692 v 3 =0.55555555555 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacqaH3oaAdaWgaaWcbaGaaGymaa qabaGccqGH9aqpcqGHsislcaaIWaGaaiOlaiaaiEdacaaI3aGaaGin aiaaiwdacaaI5aGaaGOnaiaaiAdacaaI2aGaaGyoaiaaikdacaaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caWG2bWaaSbaaSqaaiaaigdaaeqaaOGaeyypa0JaaGimaiaac6caca aI1aGaaGynaiaaiwdacaaI1aGaaGynaiaaiwdacaaI1aGaaGynaiaa iwdacaaI1aGaaGynaaqaaiabeE7aOnaaBaaaleaacaaIYaaabeaaki abg2da9iaaicdacaGGUaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaadAhadaWgaaWcbaGaaGOmaaqabaGccqGH9a qpcaaIWaGaaiOlaiaaiIdacaaI4aGaaGioaiaaiIdacaaI4aGaaGio aiaaiIdacaaI4aGaaGioaiaaiIdacaaI4aaabaGaeq4TdG2aaSbaaS qaaiaaiodaaeqaaOGaeyypa0JaaGimaiaac6cacaaI3aGaaG4naiaa isdacaaI1aGaaGyoaiaaiAdacaaI2aGaaGOnaiaaiMdacaaIYaGaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaamODamaaBaaaleaacaaIZaaabeaaki abg2da9iaaicdacaGGUaGaaGynaiaaiwdacaaI1aGaaGynaiaaiwda caaI1aGaaGynaiaaiwdacaaI1aGaaGynaiaaiwdaaaaa@D3F2@

 

Choosing the number of integration points: There are two considerations.  If too many integration points are used, time is wasted without gaining any accuracy).  If too few integration points are used, the stiffness matrix may be singular, or else the rate of convergence to the exact solution with mesh refinement will be reduced.  The following schemes will avoid both

 

Number of integration points for fully integrated elements

Linear triangle (3 nodes):  1 point

Quadratic triangle (6 nodes): 4 points

Linear quadrilateral (4 nodes): 4 points

Quadratic quadrilateral (8 nodes): 9 points

Linear tetrahedron (4 nodes): 1 point

Quadratic tetrahedron (10 nodes): 4 points

Linear brick (8 nodes): 8 points

Quadratic brick (20 nodes): 27 points

 

There are situations where it is preferable to use fewer integration points and purposely make the stiffness singular.   These are discussed in more detail in Section 8.5.

 

 

 

8.1.13 Summary of formulas for element stiffness and force matrices

 

With these definitions, then, we write the element stiffness matrix as

k aibk = I=1 N I w I C ijkl N a ( ξ i I ) ξ p ξ p x j N b ( ξ i I ) x q ξ q x l J( ξ i I ) f i a = I=1 N w I b i N a ( ξ i I )J + I=1 N ^ w I t i * N a ( ξ i I ) J MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGceaqabeaacaWGRbWaaSbaaSqaaiaadggaca WGPbGaamOyaiaadUgaaeqaaOGaeyypa0ZaaabCaeaacaWG3bWaaSba aSqaaiaadMeaaeqaaOGaam4qamaaBaaaleaacaWGPbGaamOAaiaadU gacaWGSbaabeaakmaalaaabaGaeyOaIyRaamOtamaaCaaaleqabaGa amyyaaaakiaacIcacqaH+oaEdaqhaaWcbaGaamyAaaqaaiaadMeaaa GccaGGPaaabaGaeyOaIyRaeqOVdG3aaSbaaSqaaiaadchaaeqaaaaa kmaalaaabaGaeyOaIyRaeqOVdG3aaSbaaSqaaiaadchaaeqaaaGcba GaeyOaIyRaamiEamaaBaaaleaacaWGQbaabeaaaaGcdaWcaaqaaiab gkGi2kaad6eadaahaaWcbeqaaiaadkgaaaGccaGGOaGaeqOVdG3aa0 baaSqaaiaadMgaaeaacaWGjbaaaOGaaiykaaqaaiabgkGi2kaadIha daWgaaWcbaGaamyCaaqabaaaaOWaaSaaaeaacqGHciITcqaH+oaEda WgaaWcbaGaamyCaaqabaaakeaacqGHciITcaWG4bWaaSbaaSqaaiaa dYgaaeqaaaaakiaadQeacaGGOaGaeqOVdG3aa0baaSqaaiaadMgaae aacaWGjbaaaOGaaiykaaWcbaGaamysaiabg2da9iaaigdaaeaacaWG obWaaSbaaWqaaiaadMeaaeqaaaqdcqGHris5aaGcbaGaamOzamaaDa aaleaacaWGPbaabaGaamyyaaaakiabg2da9maaqahabaGaam4Damaa BaaaleaacaWGjbaabeaakiaadkgadaWgaaWcbaGaamyAaaqabaGcca WGobWaaWbaaSqabeaacaWGHbaaaOGaaiikaiabe67a4naaDaaaleaa caWGPbaabaGaamysaaaakiaacMcacaWGkbaaleaacaWGjbGaeyypa0 JaaGymaaqaaiaad6eaa0GaeyyeIuoakiabgUcaRmaaqahabaGaam4D amaaBaaaleaacaWGjbaabeaakiaadshadaqhaaWcbaGaamyAaaqaai aacQcaaaGccaWGobWaaWbaaSqabeaacaWGHbaaaOGaaiikaiabe67a 4naaDaaaleaacaWGPbaabaGaamysaaaakiaacMcaceWGkbGbambaaS qaaiaadMeacqGH9aqpcaaIXaaabaGabmOtayaajaaaniabggHiLdGc caaMc8oaaaa@9CFA@

where

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@

 

 

 

8.1.14 Sample 2D/3D linear elastostatic FEM code

 

You can find a MAPLE implementation of a simple 2D/3D static linear elasticity code in the file   FEM_2Dor3D_linelast_standard.mws

 

The code reads an input file.  Several examples are provided :

1.      Linear_elastic_triangles.txt: Simple 2D plane strain problem with two triangular elements

2.      Linear_elastic_quad4.txt: Simple 2D plane strain problem with eight 4 noded quadrilateral elements

3.      Linear_elastic_quad8.txt: Simple 2D plane strain problem with two 8 noded quadrilateral elements

4.      Linear_elastic_brick4.txt: Simple 3D problem with 8 noded brick elements

No._material_props:    3

    Shear_modulus:   10.

    Poissons_ratio:  0.3
    Plane strain/stress: 1

No._coords_per_node:   2

No._DOF_per_node:      2

No._nodes:             4

Nodal_coords:

    0.0   0.0

    1.0   0.0

    1.0   1.0

    0.0   1.0

No._elements:                       2

Max_no._nodes_on_any_one_element:   3

element_identifier; no._nodes_on_element; connectivity:

    1  3  1 2 4

    1  3  2 3 4

No._nodes_with_prescribed_DOFs:  3

Node_#, DOF#, Value:

   1 1 0.0

 

 

As an example, we show how to run the program with the first input file. The file sets up the problem illustrated in the figure above. The elements are linear elastic plane strain with μ=10,ν=0.3 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiVd0Maeyypa0JaaGymaiaaicdaca GGSaGaaGPaVlaaykW7cqaH9oGBcqGH9aqpcaaIWaGaaiOlaiaaioda aaa@3DAD@

 

The program input file is listed on the right. Here is a brief explanation of the data in the file

1.      The first part of the input file specifies material properties.  A number ‘1’ on the Plane strain/stress line indicates a plane strain analysis; a number ‘0’ indicates plane stress.

2.      The second part specifies properties and coordinates of the nodes.  For a 2D problem each node has 2 coordinates and 2 DOF; for a 3D problem each node has 3 coordinates and 3DOF.  Then enter nodal coordinates for each node.

3.      The third part lists the element properties.  Here, you must specify the number of elements, and the maximum number of nodes on any one element (you can mix element types if you like).  Then you must specify the nodes connected to each element (known as element connectivity).   For each element, you must specify the number of nodes attached to the element; an identifier that specifies the element type (you can enter any number in this version of the code MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  the identifier is provided to allow addition of more sophisticated element types such as reduced integration elements), then enter the nodes on each element following the convention shown earlier.

4.      The fourth part of the file specifies boundary constraints.  For any constrained displacements, enter the node number, the displacement component to be prescribed, and its value.

5.      The last part of the file specifies distributed loading acting on the element faces.  The loading is assumed to be uniform.  For each loaded boundary, you should specify the element number, the face of the element (the face numbering convention was described in section 7.2.9 and 7.2.10 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  note that you must be consistent in numbering nodes and faces on each element), and the components of traction acting on the element face, as a vector with 2 or 3 components.

 

Note that the program performs absolutely no error checking on the input file.  If you put in a typo, you will get some bizzarre error message from MAPLE MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  often during element stiffness assembly.

 

For the input file shown, the program produces an output file that looks like this

The code prints the displacements at each node in the mesh, and also the strains and stresses at each integration point (where these quantities are most accurate)  for each element.

 

To run the code, you must complete the following steps

1.      Open the maple executable file;

1.      Edit the code to insert the full path for the input file in the line near the top of the code that reads

> # Change the name of the file below to point to your input file

> infile :=fopen(`D:/fullpathoffile/Linear_elastic_triangles.txt`,READ):

2.      Scroll down near the bottom to the line that reads

> #      Print nodal displacements, element strains and stresses to a file

> #

> outfile := fopen(`path/Linear_elastic_triangles.out`,WRITE):

and enter a name for the output file.

3.      Return to the top of the file, and press <enter> to execute each MAPLE block. You will see the code plot the undeformed and deformed finite element mesh at the end.  The stresses and strains in the elements are printed to the output file.