Chapter 8

 

Theory and Implementation of the Finite Element Method

 

 

 

8.2 The Finite Element Method for dynamic linear elasticity

 

In this section, we show how to extend the finite element method to solve dynamic problems. Specifically:

1.      The principle of virtual work is used to derive a discrete system of equations for the (time varying) nodal displacements

2.      Three methods for integrating the equation of motion with respect to time are presented: (i) explicit time-stepping; (ii) implicit time stepping; and (iii) modal analysis.  The properties of each scheme are illustrated by solving simple problems.

3.      As always example codes are provided so you can see how the method is actually implemented, and explore its predictions for yourself.

 

8.2.1 Review of the governing equations of dynamic linear elasticity

 

As before, we begin by summarizing the governing equations for a standard dynamic linear elasticity problem.   Given:

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

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.      The initial displacement field in the solid u o (x) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaCbiaeaacaWH1baaleqabaGaam4Baa aakiaacIcacaWH4bGaaiykaaaa@356E@ , and the initial velocity field v o (x) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaCbiaeaacaWH2baaleqabaGaam4Baa aakiaacIcacaWH4bGaaiykaaaa@356F@

6.      A body force distribution b(x,t) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCOyaiaacIcacaWH4bGaaiilaiaads hacaGGPaaaaa@35BB@  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)

7.      Boundary conditions, specifying displacements u * (x,t) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCyDamaaCaaaleqabaGaaiOkaaaaki aacIcacaWH4bGaaiilaiaadshacaGGPaaaaa@36B5@  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 t * (x,t) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCiDamaaCaaaleqabaGaaiOkaaaaki aacIcacaWH4bGaaiilaiaadshacaGGPaaaaa@36B4@  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 (time varying) 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 dynamic 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 motion for stresses σ ij / x i + b j =ρ 2 u j / t 2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIyRaeq4Wdm3aaSbaaSqaaiaadM gacaWGQbaabeaakiaac+cacqGHciITcaWG4bWaaSbaaSqaaiaadMga aeqaaOGaey4kaSIaamOyamaaBaaaleaacaWGQbaabeaakiabg2da9i abeg8aYjabgkGi2oaaCaaaleqabaGaaGOmaaaakiaadwhadaWgaaWc baGaamOAaaqabaGccaGGVaGaeyOaIyRaamiDamaaCaaaleqabaGaaG Omaaaaaaa@486C@

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@

 

 

8.2.2 Expressing the governing equations using the principle of virtual work

 

Just as for static problems, 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+ R ρ 2 u i t 2 δ 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 biaabeqaaiqabaWaaaGceaqabeaadaWdrbqaaiaadoeadaWgaaWcba GaamyAaiaadQgacaWGRbGaamiBaaqabaGcdaWcaaqaaiabgkGi2kaa dwhadaWgaaWcbaGaam4AaaqabaaakeaacqGHciITcaWG4bWaaSbaaS qaaiaadYgaaeqaaaaakmaalaaabaGaeyOaIyRaeqiTdqMaamODamaa BaaaleaacaWGPbaabeaaaOqaaiabgkGi2kaadIhadaWgaaWcbaGaam OAaaqabaaaaOGaamizaiaadAfacqGHsisldaWdrbqaaiaadkgadaWg aaWcbaGaamyAaaqabaGccqaH0oazcaWG2bWaaSbaaSqaaiaadMgaae qaaOGaamizaiaadAfacqGHRaWkdaWdrbqaaiabeg8aYnaalaaabaGa eyOaIy7aaWbaaSqabeaacaaIYaaaaOGaamyDamaaBaaaleaacaWGPb aabeaaaOqaaiabgkGi2kaadshadaahaaWcbeqaaiaaikdaaaaaaOGa eqiTdqMaamODamaaBaaaleaacaWGPbaabeaakiaadsgacaWGwbaale aacaWGsbaabeqdcqGHRiI8aOGaeyOeI0Yaa8quaeaacaWG0bWaa0ba aSqaaiaadMgaaeaacaGGQaaaaOGaeqiTdqMaamODamaaBaaaleaaca WGPbaabeaakiaadsgacaWGbbGaeyypa0JaaGimaaWcbaGaeyOaIy7a aSbaaWqaaiaaikdaaeqaaSGaamOuaaqab0Gaey4kIipaaSqaaiaadk faaeqaniabgUIiYdaaleaacaWGsbaabeqdcqGHRiI8aaGcbaGaamyD amaaBaaaleaacaWGPbaabeaakiabg2da9iaadwhadaqhaaWcbaGaam yAaaqaaiaacQcaaaGccaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8Uaae4Baiaab6gacaaMc8UaaG PaVlabgkGi2oaaBaaaleaacaaIXaaabeaakiaadkfaaaaa@971F@

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 equation of motion and boundary conditions, so all the field equations and boundary conditions will be satisfied.

 

You can derive this result by following the method outlined in Section 8.1.2.

 

 

8.2.3 The finite element equations of motion for linear elastic solids


We now introduce a finite element approximate following exactly the same procedure that we used for static problems. We choose to calculate the displacement field (now a function of time) at a set of n discrete nodes.  We 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@ .  We then set up the FEA equations for dynamic problems as follows:

1.      As before, the displacement and virtual velocity fields are interpolated between nodal values 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@        δ 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@

2.      Substitute into the virtual work equation to see that

R ρ N b N a 2 u i b t 2 δ v i a dV + R C ijkl N b x l N a x j u k b δ v i a dV R b i N a δ v i a dV 2R t i * N a δ v i a dA =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaa8quaeaacqaHbpGCcaWGobWaaWbaaS qabeaacaWGIbaaaOGaamOtamaaCaaaleqabaGaamyyaaaakmaalaaa baGaeyOaIy7aaWbaaSqabeaacaaIYaaaaOGaamyDamaaDaaaleaaca WGPbaabaGaamOyaaaaaOqaaiabgkGi2kaadshadaahaaWcbeqaaiaa ikdaaaaaaOGaeqiTdqMaamODamaaDaaaleaacaWGPbaabaGaamyyaa aakiaadsgacaWGwbaaleaacaWGsbaabeqdcqGHRiI8aOGaey4kaSYa a8quaeaacaWGdbWaaSbaaSqaaiaadMgacaWGQbGaam4AaiaadYgaae qaaOWaaSaaaeaacqGHciITcaWGobWaa0baaSqaaaqaaiaadkgaaaaa keaacqGHciITcaWG4bWaaSbaaSqaaiaadYgaaeqaaaaakmaalaaaba GaeyOaIyRaamOtamaaDaaaleaaaeaacaWGHbaaaaGcbaGaeyOaIyRa amiEamaaBaaaleaacaWGQbaabeaaaaGccaWG1bWaa0baaSqaaiaadU gaaeaacaWGIbaaaOGaeqiTdqMaamODamaaDaaaleaacaWGPbaabaGa amyyaaaakiaadsgacaWGwbaaleaacaWGsbaabeqdcqGHRiI8aOGaey OeI0Yaa8quaeaacaWGIbWaaSbaaSqaaiaadMgaaeqaaOGaamOtamaa CaaaleqabaGaamyyaaaakiabes7aKjaadAhadaqhaaWcbaGaamyAaa qaaiaadggaaaGccaWGKbGaamOvaaWcbaGaamOuaaqab0Gaey4kIipa kiabgkHiTmaapefabaGaamiDamaaDaaaleaacaWGPbaabaGaaiOkaa aakiaad6eadaahaaWcbeqaaiaadggaaaGccqaH0oazcaWG2bWaa0ba aSqaaiaadMgaaeaacaWGHbaaaOGaamizaiaadgeaaSqaaiabgkGi2k aaikdacaWGsbaabeqdcqGHRiI8aOGaeyypa0JaaGimaaaa@8A7C@

3.      This is once again a linear system which may be expressed in the form

( M ab u ¨ i b + 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 biaabeqaaiqabaWaaaGcbaWaaeWaaeaacaWGnbWaaSbaaSqaaiaadg gacaWGIbaabeaakiqadwhagaWaamaaDaaaleaacaWGPbaabaGaamOy aaaakiabgUcaRiaadUeadaWgaaWcbaGaamyyaiaadMgacaWGIbGaam 4AaaqabaGccaWG1bWaa0baaSqaaiaadUgaaeaacaWGIbaaaOGaeyOe I0IaamOramaaDaaaleaacaWGPbaabaGaamyyaaaaaOGaayjkaiaawM caaiabes7aKjaadAhadaqhaaWcbaGaamyAaaqaaiaadggaaaGccqGH 9aqpcaaIWaaaaa@4B06@

where

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 PaVlaaykW7caWGgbWaa0baaSqaaiaadMgaaeaacaWGHbaaaOGaeyyp a0Zaa8quaeaacaWGIbWaaSbaaSqaaiaadMgaaeqaaOGaamOtamaaCa aaleqabaGaamyyaaaakiaacIcacaWH4bGaaiykaaWcbaGaamOuaaqa b0Gaey4kIipakiaadsgacaWGwbGaey4kaSYaa8quaeaacaWG0bWaa0 baaSqaaiaadMgaaeaacaGGQaaaaOGaamOtamaaCaaaleqabaGaamyy aaaakiaacIcacaWH4bGaaiykaaWcbaGaeyOaIyRaaGOmaiaadkfaae qaniabgUIiYdGccaWGKbGaamyqaiaaykW7aaa@8671@

are the finite element stiffness matrix and force vector introduced in Section 7.2, and

M ab = R ρ N a N b dV MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpdaWdrbqaaiabeg8aYjaad6eadaahaaWcbeqaaiaa dggaaaGccaWGobWaaWbaaSqabeaacaWGIbaaaOGaamizaiaadAfaaS qaaiaadkfaaeqaniabgUIiYdaaaa@3F30@

is known as the finite element mass matrix.  Procedures for computing the mass matrix are given below.

4.      Again, since the principal of virtual power holds for all virtual velocity fields such that δ v i a =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqiTdqMaamODamaaDaaaleaacaWGPb aabaGaamyyaaaakiabg2da9iaaicdaaaa@373A@  for x i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8skY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiEamaaDaaaleaacaWGPbaabaGaam yyaaaaaaa@33DD@  on 1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIy7aaSbaaSqaaiaaigdaaeqaaO GaamOuaaaa@340B@  we conclude that this requires

( M ab 2 u i b t 2 + K aibk u k b F i a )=0(i,a): x k a not on 1 R u i a = u i * (i,a): x k a  on 1 R MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaadaqadaqaaiaad2eadaWgaaWcba GaamyyaiaadkgaaeqaaOWaaSaaaeaacqGHciITdaahaaWcbeqaaiaa ikdaaaGccaWG1bWaa0baaSqaaiaadMgaaeaacaWGIbaaaaGcbaGaey OaIyRaamiDamaaCaaaleqabaGaaGOmaaaaaaGccqGHRaWkcaWGlbWa aSbaaSqaaiaadggacaWGPbGaamOyaiaadUgaaeqaaOGaamyDamaaDa aaleaacaWGRbaabaGaamOyaaaakiabgkHiTiaadAeadaqhaaWcbaGa amyAaaqaaiaadggaaaaakiaawIcacaGLPaaacqGH9aqpcaaIWaGaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7cqGHaiIicaGGOaGaamyAaiaacYcacaWG HbGaaiykaiaacQdacaaMc8UaaGPaVlaadIhadaqhaaWcbaGaam4Aaa qaaiaadggaaaGccaaMc8UaaGPaVlaaykW7caqGUbGaae4Baiaabsha caqGGaGaae4Baiaab6gacaaMc8UaaGPaVlabgkGi2oaaBaaaleaaca qGXaaabeaakiaadkfaaeaacaWG1bWaa0baaSqaaiaadMgaaeaacaWG HbaaaOGaeyypa0JaamyDamaaDaaaleaacaWGPbaabaGaaiOkaaaaki aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7cqGHaiIicaGGOa GaamyAaiaacYcacaWGHbGaaiykaiaacQdacaaMc8UaaGPaVlaadIha daqhaaWcbaGaam4AaaqaaiaadggaaaGccaaMc8UaaGPaVlaaykW7ca qGGaGaae4Baiaab6gacaaMc8UaaGPaVlabgkGi2oaaBaaaleaacaqG Xaaabeaakiaadkfaaaaa@2064@

Nodal velocities must also satisfy initial conditions. This is a set of coupled second-order linear differential equations that may be integrated with respect to time to determine u i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaDaaaleaacaWGPbaabaGaam yyaaaaaaa@33CA@  as a function of time.  Observe that M ab MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaaaaa@339A@  and K aibk MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamyAai aadkgacaWGRbaabeaaaaa@3576@  are constant matrices for linear elastic problems, but F i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOramaaDaaaleaacaWGPbaabaGaam yyaaaaaaa@339B@  will in general be a function of time.

 

To find the solution, we now need to integrate these equations with respect to time.  For linear problems, there are two choices

1.      Brute-force time stepping methods (e.g. Newmark time integration).  This technique can be used for both linear and nonlinear problems.

2.      Modal methods, which integrate the equations of motion exactly.  This method only works for linear problems.

 

Both methods are described in detail below.

 

 

8.2.4 Newmark time integration for elastodynamics

 

The general idea is simple.  Given values of u i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaDaaaleaacaWGPbaabaGaam yyaaaaaaa@33CA@  and u i a /t MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeyOaIyRaamyDamaaDaaaleaacaWGPb aabaGaamyyaaaakiaac+cacqGHciITcaWG0baaaa@384C@  at some time t, we wish to determine values at a slightly later time t+Δt MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiDaiabgUcaRiabfs5aejaadshaaa a@3509@ , using the governing equations of motion.  There are of course many ways to do this.  Here, we will present the popular ‘Newmark’ time integration scheme.

 

Newmark method for a 1 DOF system: It is simplest to illustrate the procedure using a 1-DOF system (e.g. a forced,  undamped spring-mass system) first.   For this case the equation of motion is

m u ¨ +kuf=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyBaiqadwhagaWaaiabgUcaRiaadU gacaWG1bGaeyOeI0IaamOzaiabg2da9iaaicdaaaa@3927@

where f(t) and the initial conditions u(0), u ˙ (0) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDaiaacIcacaaIWaGaaiykaiaacY cacaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlqadwhagaGaaiaacIca caaIWaGaaiykaaaa@3F59@  are given.

 

Suppose that we are able to get estimates for the acceleration u ¨ (t) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabmyDayaadaGaaiikaiaadshacaGGPa aaaa@3425@  and u ¨ (t+Δt) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabmyDayaadaGaaiikaiaadshacqGHRa WkcqqHuoarcaWG0bGaaiykaaaa@3766@  both at the start and end of a general time step Δt MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeuiLdqKaamiDaaaa@332E@ .  We could then use a Taylor series expansion to obtain estimates of displacement and velocity at time t+Δt MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiDaiabgUcaRiabfs5aejaadshaaa a@3509@

u(t+Δt)u(t)+Δt u ˙ (t)+ Δ t 2 2 [ (1 β 2 ) u ¨ (t)+ β 2 u ¨ (t+Δt) ] u ˙ (t+Δt) u ˙ (t)+Δt[ (1 β 1 ) u ¨ (t)+ β 1 u ¨ (t+Δt) ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWG1bGaaiikaiaadshacqGHRa WkcqqHuoarcaWG0bGaaiykaiabgIKi7kaadwhacaGGOaGaamiDaiaa cMcacqGHRaWkcqqHuoarcaWG0bGabmyDayaacaGaaiikaiaadshaca GGPaGaey4kaSYaaSaaaeaacqqHuoarcaWG0bWaaWbaaSqabeaacaaI YaaaaaGcbaGaaGOmaaaadaWadaqaaiaacIcacaaIXaGaeyOeI0Iaeq OSdi2aaSbaaSqaaiaaikdaaeqaaOGaaiykaiqadwhagaWaaiaacIca caWG0bGaaiykaiabgUcaRiabek7aInaaBaaaleaacaaIYaaabeaaki qadwhagaWaaiaacIcacaWG0bGaey4kaSIaeuiLdqKaamiDaiaacMca aiaawUfacaGLDbaaaeaaceWG1bGbaiaacaGGOaGaamiDaiabgUcaRi abfs5aejaadshacaGGPaGaeyisISRabmyDayaacaGaaiikaiaadsha caGGPaGaey4kaSIaeuiLdqKaamiDamaadmaabaGaaiikaiaaigdacq GHsislcqaHYoGydaWgaaWcbaGaaGymaaqabaGccaGGPaGabmyDayaa daGaaiikaiaadshacaGGPaGaey4kaSIaeqOSdi2aaSbaaSqaaiaaig daaeqaaOGabmyDayaadaGaaiikaiaadshacqGHRaWkcqqHuoarcaWG 0bGaaiykaaGaay5waiaaw2faaaaaaa@80A1@

Here, β 1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdi2aaSbaaSqaaiaaigdaaeqaaa aa@3357@  and β 2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdi2aaSbaaSqaaiaaikdaaeqaaa aa@3358@  are two adjustable parameters that determine the nature of the time integration scheme.  If we set β 1 = β 2 =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdi2aaSbaaSqaaiaaigdaaeqaaO Gaeyypa0JaeqOSdi2aaSbaaSqaaiaaikdaaeqaaOGaeyypa0JaaGim aaaa@38BA@  the acceleration is estimated based on its value at time t.  This is known as an explicit time integration scheme.  Alternatively, if we put β 1 = β 2 =1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdi2aaSbaaSqaaiaaigdaaeqaaO Gaeyypa0JaeqOSdi2aaSbaaSqaaiaaikdaaeqaaOGaeyypa0JaaGym aaaa@38BB@ , the acceleration is estimated from its value at time t+Δt MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiDaiabgUcaRiabfs5aejaadshaaa a@3507@ .  This is known as implicit time integration.

 

The acceleration at the start and end of the time step is computed using the equation of motion.  At time t we have

u ¨ (t)= 1 m [ ku(t)+f(t) ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabmyDayaadaGaaiikaiaadshacaGGPa Gaeyypa0ZaaSaaaeaacaaIXaaabaGaamyBaaaadaWadaqaaiabgkHi TiaadUgacaWG1bGaaiikaiaadshacaGGPaGaey4kaSIaamOzaiaacI cacaWG0bGaaiykaaGaay5waiaaw2faaaaa@4220@

while at time t+Δt MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiDaiabgUcaRiabfs5aejaadshaaa a@3507@

m u ¨ (t+Δt)+ku(t+Δt)f(t+Δt)=0 m u ¨ (t+Δt)+k{ u(t)+Δt u ˙ (t)+ Δ t 2 2 [ (1 β 2 ) u ¨ (t)+ β 2 u ¨ (t+Δt) ] }f(t+Δt)=0 u ¨ (t+Δt)= 1 m+k β 2 Δ t 2 { k[ u(t)+Δt u ˙ (t)+ Δ t 2 2 (1 β 2 ) u ¨ (t) ]+f(t+Δt) } MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWGTbGabmyDayaadaGaaiikai aadshacqGHRaWkcqqHuoarcaWG0bGaaiykaiabgUcaRiaadUgacaWG 1bGaaiikaiaadshacqGHRaWkcqqHuoarcaWG0bGaaiykaiabgkHiTi aadAgacaGGOaGaamiDaiabgUcaRiabfs5aejaadshacaGGPaGaeyyp a0JaaGimaaqaaiabgkDiElaad2gaceWG1bGbamaacaGGOaGaamiDai abgUcaRiabfs5aejaadshacaGGPaGaey4kaSIaam4AamaacmaabaGa amyDaiaacIcacaWG0bGaaiykaiabgUcaRiabfs5aejaadshaceWG1b GbaiaacaGGOaGaamiDaiaacMcacqGHRaWkdaWcaaqaaiabfs5aejaa dshadaahaaWcbeqaaiaaikdaaaaakeaacaaIYaaaamaadmaabaGaai ikaiaaigdacqGHsislcqaHYoGydaWgaaWcbaGaaGOmaaqabaGccaGG PaGabmyDayaadaGaaiikaiaadshacaGGPaGaey4kaSIaeqOSdi2aaS baaSqaaiaaikdaaeqaaOGabmyDayaadaGaaiikaiaadshacqGHRaWk cqqHuoarcaWG0bGaaiykaaGaay5waiaaw2faaaGaay5Eaiaaw2haai abgkHiTiaadAgacaGGOaGaamiDaiabgUcaRiabfs5aejaadshacaGG PaGaeyypa0JaaGimaaqaaiabgkDiElqadwhagaWaaiaacIcacaWG0b Gaey4kaSIaeuiLdqKaamiDaiaacMcacqGH9aqpdaWcaaqaaiaaigda aeaacaWGTbGaey4kaSIaam4Aaiabek7aInaaBaaaleaacaaIYaaabe aakiabfs5aejaadshadaahaaWcbeqaaiaaikdaaaaaaOWaaiWaaeaa cqGHsislcaWGRbWaamWaaeaacaWG1bGaaiikaiaadshacaGGPaGaey 4kaSIaeuiLdqKaamiDaiqadwhagaGaaiaacIcacaWG0bGaaiykaiab gUcaRmaalaaabaGaeuiLdqKaamiDamaaCaaaleqabaGaaGOmaaaaaO qaaiaaikdaaaGaaiikaiaaigdacqGHsislcqaHYoGydaWgaaWcbaGa aGOmaaqabaGccaGGPaGabmyDayaadaGaaiikaiaadshacaGGPaaaca GLBbGaayzxaaGaey4kaSIaamOzaiaacIcacaWG0bGaey4kaSIaeuiL dqKaamiDaiaacMcaaiaawUhacaGL9baaaaaa@BD62@

 

This then gives us the following time stepping scheme: Given: k, m, f(t), u(0), u ˙ (0) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDaiaacIcacaaIWaGaaiykaiaacY cacaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlqadwhagaGaaiaacIca caaIWaGaaiykaaaa@3F59@

1.      At t=0 compute u ¨ (0)= 1 m [ ku(0)+f(0) ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabmyDayaadaGaaiikaiaaicdacaGGPa Gaeyypa0ZaaSaaaeaacaaIXaaabaGaamyBaaaadaWadaqaaiabgkHi TiaadUgacaWG1bGaaiikaiaaicdacaGGPaGaey4kaSIaamOzaiaacI cacaaIWaGaaiykaaGaay5waiaaw2faaaaa@4165@ .

2.      Then for successive time steps, compute

u ¨ (t+Δt)= 1 m+k β 2 Δ t 2 { k[ u(t)+Δt u ˙ (t)+ Δ t 2 2 (1 β 2 ) u ¨ (t) ]+f(t+Δt) } MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabmyDayaadaGaaiikaiaadshacqGHRa WkcqqHuoarcaWG0bGaaiykaiabg2da9maalaaabaGaaGymaaqaaiaa d2gacqGHRaWkcaWGRbGaeqOSdi2aaSbaaSqaaiaaikdaaeqaaOGaeu iLdqKaamiDamaaCaaaleqabaGaaGOmaaaaaaGcdaGadaqaaiabgkHi TiaadUgadaWadaqaaiaadwhacaGGOaGaamiDaiaacMcacqGHRaWkcq qHuoarcaWG0bGabmyDayaacaGaaiikaiaadshacaGGPaGaey4kaSYa aSaaaeaacqqHuoarcaWG0bWaaWbaaSqabeaacaaIYaaaaaGcbaGaaG OmaaaacaGGOaGaaGymaiabgkHiTiabek7aInaaBaaaleaacaaIYaaa beaakiaacMcaceWG1bGbamaacaGGOaGaamiDaiaacMcaaiaawUfaca GLDbaacqGHRaWkcaWGMbGaaiikaiaadshacqGHRaWkcqqHuoarcaWG 0bGaaiykaaGaay5Eaiaaw2haaaaa@670C@

u(t+Δt)u(t)+Δt u ˙ (t)+ Δ t 2 2 [ (1 β 2 ) u ¨ (t)+ β 2 u ¨ (t+Δt) ] u ˙ (t+Δt) u ˙ (t)+Δt[ (1 β 1 ) u ¨ (t)+ β 1 u ¨ (t+Δt) ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWG1bGaaiikaiaadshacqGHRa WkcqqHuoarcaWG0bGaaiykaiabgIKi7kaadwhacaGGOaGaamiDaiaa cMcacqGHRaWkcqqHuoarcaWG0bGabmyDayaacaGaaiikaiaadshaca GGPaGaey4kaSYaaSaaaeaacqqHuoarcaWG0bWaaWbaaSqabeaacaaI YaaaaaGcbaGaaGOmaaaadaWadaqaaiaacIcacaaIXaGaeyOeI0Iaeq OSdi2aaSbaaSqaaiaaikdaaeqaaOGaaiykaiqadwhagaWaaiaacIca caWG0bGaaiykaiabgUcaRiabek7aInaaBaaaleaacaaIYaaabeaaki qadwhagaWaaiaacIcacaWG0bGaey4kaSIaeuiLdqKaamiDaiaacMca aiaawUfacaGLDbaaaeaaceWG1bGbaiaacaGGOaGaamiDaiabgUcaRi abfs5aejaadshacaGGPaGaeyisISRabmyDayaacaGaaiikaiaadsha caGGPaGaey4kaSIaeuiLdqKaamiDamaadmaabaGaaiikaiaaigdacq GHsislcqaHYoGydaWgaaWcbaGaaGymaaqabaGccaGGPaGabmyDayaa daGaaiikaiaadshacaGGPaGaey4kaSIaeqOSdi2aaSbaaSqaaiaaig daaeqaaOGabmyDayaadaGaaiikaiaadshacqGHRaWkcqqHuoarcaWG 0bGaaiykaaGaay5waiaaw2faaaaaaa@80A1@

 

 

Newmark applied to FEM equations: This can immediately be extended to the general n DOF case as follows. Given K aibk MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaam4samaaBaaaleaacaWGHbGaamyAai aadkgacaWGRbaabeaaaaa@3576@ , M ab MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaaaaa@339A@ , F i a (t) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOramaaDaaaleaacaWGPbaabaGaam yyaaaakiaacIcacaWG0bGaaiykaaaa@35F7@ , u i a (0), u ˙ i a (0) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaDaaaleaacaWGPbaabaGaam yyaaaakiaacIcacaaIWaGaaiykaiaacYcacaaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7ceWG1bGbaiaadaqhaaWcbaGaamyAaaqaai aadggaaaGccaGGOaGaaGimaiaacMcaaaa@44FA@

1.      At  t=0, compute u ¨ i a (0)= M ab 1 [ K bick u k c (0)+ F i b (0) ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabmyDayaadaWaa0baaSqaaiaadMgaae aacaWGHbaaaOGaaiikaiaaicdacaGGPaGaeyypa0JaamytamaaDaaa leaacaWGHbGaamOyaaqaaiabgkHiTiaaigdaaaGcdaWadaqaaiabgk HiTiaadUeadaWgaaWcbaGaamOyaiaadMgacaWGJbGaam4AaaqabaGc caWG1bWaa0baaSqaaiaadUgaaeaacaWGJbaaaOGaaiikaiaaicdaca GGPaGaey4kaSIaamOramaaDaaaleaacaWGPbaabaGaamOyaaaakiaa cIcacaaIWaGaaiykaaGaay5waiaaw2faaaaa@4DEF@

2.      Then for successive time steps, solve

M ab u ¨ i b (t+Δt)+ β 2 Δ t 2 K aibk u ¨ k b = K aibk [ u k b (t)+Δt u ˙ k b (t)+ Δ t 2 2 (1 β 2 ) u ¨ k b (t) ]+ F i a (t+Δt) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaGcceWG1bGbamaadaqhaaWcbaGaamyAaaqaaiaadkgaaaGccaGG OaGaamiDaiabgUcaRiabfs5aejaadshacaGGPaGaey4kaSIaeqOSdi 2aaSbaaSqaaiaaikdaaeqaaOGaeuiLdqKaamiDamaaCaaaleqabaGa aGOmaaaakiaadUeadaWgaaWcbaGaamyyaiaadMgacaWGIbGaam4Aaa qabaGcceWG1bGbamaadaqhaaWcbaGaam4AaaqaaiaadkgaaaGccqGH 9aqpcqGHsislcaWGlbWaaSbaaSqaaiaadggacaWGPbGaamOyaiaadU gaaeqaaOWaamWaaeaacaWG1bWaa0baaSqaaiaadUgaaeaacaWGIbaa aOGaaiikaiaadshacaGGPaGaey4kaSIaeuiLdqKaamiDaiqadwhaga GaamaaDaaaleaacaWGRbaabaGaamOyaaaakiaacIcacaWG0bGaaiyk aiabgUcaRmaalaaabaGaeuiLdqKaamiDamaaCaaaleqabaGaaGOmaa aaaOqaaiaaikdaaaGaaiikaiaaigdacqGHsislcqaHYoGydaWgaaWc baGaaGOmaaqabaGccaGGPaGabmyDayaadaWaa0baaSqaaiaadUgaae aacaWGIbaaaOGaaiikaiaadshacaGGPaaacaGLBbGaayzxaaGaey4k aSIaamOramaaDaaaleaacaWGPbaabaGaamyyaaaakiaacIcacaWG0b Gaey4kaSIaeuiLdqKaamiDaiaacMcaaaa@7AA8@

for u ¨ i a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabmyDayaadaWaa0baaSqaaiaadMgaae aacaWGHbaaaaaa@33D4@

3.      Then compute

u i a (t+Δt) u i a (t)+Δt u ˙ i a (t)+ Δ t 2 2 [ (1 β 2 ) u ¨ i a (t)+ β 2 u ¨ i a (t+Δt) ] u ˙ i a (t+Δt) u ˙ i a (t)+Δt[ (1 β 1 ) u ¨ i a (t)+ β 1 u ¨ i a (t+Δt) ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGceaqabeaacaWG1bWaa0baaSqaaiaadMgaae aacaWGHbaaaOGaaiikaiaadshacqGHRaWkcqqHuoarcaWG0bGaaiyk aiabgIKi7kaadwhadaqhaaWcbaGaamyAaaqaaiaadggaaaGccaGGOa GaamiDaiaacMcacqGHRaWkcqqHuoarcaWG0bGabmyDayaacaWaa0ba aSqaaiaadMgaaeaacaWGHbaaaOGaaiikaiaadshacaGGPaGaey4kaS YaaSaaaeaacqqHuoarcaWG0bWaaWbaaSqabeaacaaIYaaaaaGcbaGa aGOmaaaadaWadaqaaiaacIcacaaIXaGaeyOeI0IaeqOSdi2aaSbaaS qaaiaaikdaaeqaaOGaaiykaiqadwhagaWaamaaDaaaleaacaWGPbaa baGaamyyaaaakiaacIcacaWG0bGaaiykaiabgUcaRiabek7aInaaBa aaleaacaaIYaaabeaakiqadwhagaWaamaaDaaaleaacaWGPbaabaGa amyyaaaakiaacIcacaWG0bGaey4kaSIaeuiLdqKaamiDaiaacMcaai aawUfacaGLDbaaaeaaceWG1bGbaiaadaqhaaWcbaGaamyAaaqaaiaa dggaaaGccaGGOaGaamiDaiabgUcaRiabfs5aejaadshacaGGPaGaey isISRabmyDayaacaWaa0baaSqaaiaadMgaaeaacaWGHbaaaOGaaiik aiaadshacaGGPaGaey4kaSIaeuiLdqKaamiDamaadmaabaGaaiikai aaigdacqGHsislcqaHYoGydaWgaaWcbaGaaGymaaqabaGccaGGPaGa bmyDayaadaWaa0baaSqaaiaadMgaaeaacaWGHbaaaOGaaiikaiaads hacaGGPaGaey4kaSIaeqOSdi2aaSbaaSqaaiaaigdaaeqaaOGabmyD ayaadaWaa0baaSqaaiaadMgaaeaacaWGHbaaaOGaaiikaiaadshacq GHRaWkcqqHuoarcaWG0bGaaiykaaGaay5waiaaw2faaaaaaa@9304@

 

 

 

8.2.5 1-D implementation of a Newmark scheme.

 

It is straightforward to extend a static FEM code to dynamics.  As before, we illustrate the method in 1D before developing a 3D code.

 

Consider a long linear elastic bar.  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@   and mass density ρ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqyWdihaaa@328F@ ;

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.      The bar is at rest and strain free at t=0 with u 1 =d u 1 /dt=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamyDamaaBaaaleaacaaIXaaabeaaki abg2da9iaadsgacaWG1bWaaSbaaSqaaiaaigdaaeqaaOGaai4laiaa dsgacaWG0bGaeyypa0JaaGimaaaa@3AE9@  

4.      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@

5.      The bar is subjected to body force b=b( x 1 ,t) e 1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCOyaiabg2da9iaadkgacaGGOaGaam iEamaaBaaaleaacaaIXaaabeaakiaacYcacaWG0bGaaiykaiaahwga daWgaaWcbaGaaGymaaqabaaaaa@3A6C@ ,

6.      The bar is either loaded or constrained at its ends, so that the boundary conditions are either  t 1 = t * (0,t), t 1 = t * (L,t) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamiDamaaBaaaleaacaaIXaaabeaaki abg2da9iaadshadaahaaWcbeqaaiaacQcaaaGccaGGOaGaaGimaiaa cYcacaWG0bGaaiykaiaacYcacaaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaadshadaWgaaWcbaGaaGymaaqabaGccqGH9aqpcaWG0bWaaWba aSqabeaacaGGQaaaaOGaaiikaiaadYeacaGGSaGaamiDaiaacMcaaa a@4A61@  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

M ab u 1 b + K ab u 1 b = F a MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaGccaWG1bWaa0baaSqaaiaaigdaaeaacaWGIbaaaOGaey4kaSIa am4samaaBaaaleaacaWGHbGaamOyaaqabaGccaWG1bWaa0baaSqaai aaigdaaeaacaWGIbaaaOGaeyypa0JaamOramaaDaaaleaaaeaacaWG Hbaaaaaa@3FE3@

where

M ab = h 2 0 L ρ N a ( x 1 ) N b ( x 1 )d x 1 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 aiaabeqaaiqacaWaaaGceaqabeaacaWGnbWaaSbaaSqaaiaadggaca WGIbaabeaakiabg2da9iaadIgadaahaaWcbeqaaiaaikdaaaGcdaWd Xbqaaiabeg8aYjaad6eadaahaaWcbeqaaiaadggaaaGccaGGOaGaam iEamaaBaaaleaacaaIXaaabeaakiaacMcacaWGobWaaWbaaSqabeaa caWGIbaaaOGaaiikaiaadIhadaWgaaWcbaGaaGymaaqabaGccaGGPa GaamizaiaadIhadaWgaaWcbaGaaGymaaqabaaabaGaaGimaaqaaiaa dYeaa0Gaey4kIipakiaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Uaam4samaaBaaaleaa caWGHbGaamOyaaqabaGccqGH9aqpcaWGObWaaWbaaSqabeaacaaIYa aaaOWaa8qCaeaadaWcaaqaaiaaikdacqaH8oqBcaGGOaGaaGymaiab gkHiTiabe27aUjaacMcaaeaacaaIXaGaeyOeI0IaaGOmaiabe27aUb aadaWcaaqaaiabgkGi2kaad6eadaahaaWcbeqaaiaadggaaaGccaGG OaGaamiEamaaBaaaleaacaaIXaaabeaakiaacMcaaeaacqGHciITca WG4bWaaSbaaSqaaiaaigdaaeqaaaaakmaalaaabaGaeyOaIyRaamOt amaaCaaaleqabaGaamOyaaaakiaacIcacaWG4bWaaSbaaSqaaiaaig daaeqaaOGaaiykaaqaaiabgkGi2kaadIhadaWgaaWcbaGaaGymaaqa baaaaOGaamizaiaadIhadaWgaaWcbaGaaGymaaqabaaabaGaaGimaa qaaiaadYeaa0Gaey4kIipaaOqaaiaadAeadaqhaaWcbaaabaGaamyy aaaakiabg2da9iaadIgadaahaaWcbeqaaiaaikdaaaGcdaWdXbqaai aadkgadaWgaaWcbaaabeaakiaad6eadaahaaWcbeqaaiaadggaaaGc caGGOaGaamiEamaaBaaaleaacaaIXaaabeaakiaacMcacaWGKbGaam iEamaaBaaaleaacaaIXaaabeaaaeaacaaIWaaabaGaamitaaqdcqGH RiI8aOGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7cqGHRaWkcaWGOb WaaWbaaSqabeaacaaIYaaaaOGaamiDamaaDaaaleaacaaIXaaabaGa aiOkaaaakiaacIcacaaIWaGaaiykaiaad6eadaahaaWcbeqaaiaadg gaaaGccaGGOaGaaGimaiaacMcacqGHRaWkcaWGObWaaWbaaSqabeaa caaIYaaaaOGaamiDamaaDaaaleaaaeaacaGGQaaaaOGaaiikaiaadY eacaGGPaGaamOtamaaCaaaleqabaGaamyyaaaakiaacIcacaWGmbGa aiykaaaaaa@B639@

 

We adopt exactly the same interpolation scheme as for the static solution, and calculate stiffness matrix and force vector accordingly.  In addition, we need to determine the mass matrix.  This can also be done by evaluating integrals for each element by numerical quadrature, and then assembling element mass matrices into a global matrix.  Specifically,

 

   For each element, compute the element mass matrix. The mass matrix can either be evaluated by numerical qudrature, by evaluating the formula

m ab = h 2 I=1 n I w I ρ N a ( ξ I ) N b ( ξ I )J( ξ I ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=xi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaamyBamaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpcaWGObWaaWbaaSqabeaacaaIYaaaaOWaaabCaeaa caWG3bWaaSbaaSqaaiaadMeaaeqaaOGaeqyWdiNaamOtamaaCaaale qabaGaamyyaaaakmaabmaabaGaeqOVdG3aaSbaaSqaaiaadMeaaeqa aaGccaGLOaGaayzkaaGaamOtamaaCaaaleqabaGaamOyaaaakmaabm aabaGaeqOVdG3aaSbaaSqaaiaadMeaaeqaaaGccaGLOaGaayzkaaGa amOsaiaacIcacqaH+oaEdaWgaaWcbaGaamysaaqabaGccaGGPaaale aacaWGjbGaeyypa0JaaGymaaqaaiaad6gadaWgaaadbaGaamysaaqa baaaniabggHiLdaaaa@52B8@

 

where

J=| x ξ |=| a=1 N e N a ( ξ I ) ξ x a | MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamOsaiabg2da9maaemaabaWaaSaaae aacqGHciITcaWG4baabaGaeyOaIyRaeqOVdGhaaaGaay5bSlaawIa7 aiabg2da9maaemaabaWaaabCaeaadaWcaaqaaiabgkGi2kaad6eada ahaaWcbeqaaiaadggaaaGccaGGOaGaeqOVdG3aaSbaaSqaaiaadMea aeqaaOGaaiykaaqaaiabgkGi2kabe67a4baacaWG4bWaa0baaSqaaa qaaiaadggaaaaabaGaamyyaiabg2da9iaaigdaaeaacaWGobWaaSba aWqaaiaadwgaaeqaaaqdcqGHris5aaGccaGLhWUaayjcSdaaaa@531B@

and the integration points ξ I , w I I=1... n I MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOVdG3aaSbaaSqaaiaadMeaaeqaaO GaaiilaiaadEhadaWgaaWcbaGaamysaaqabaGccaaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caWGjbGaeyypa0JaaGymaiaac6caca GGUaGaaiOlaiaad6gadaWgaaWcbaGaamysaaqabaaaaa@461A@  , 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 in Section 7.2.5. Note that that a higher order integration scheme is required to integrate the mass matrix than we needed for the stiffness matrix.  Specifically, for quadratic elements we need a cubic

Mass matrices for 1D elements

Linear element

[m]= ρL 6 [ 2 1 1 2 ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaai4waiaad2gacaGGDbGaeyypa0ZaaS aaaeaacqaHbpGCcaWGmbaabaGaaGOnaaaadaWadaqaauaabeqaciaa aeaacaaIYaaabaGaaGymaaqaaiaaigdaaeaacaaIYaaaaaGaay5wai aaw2faaaaa@3CD8@

 

Quadratic element

(midside node must be at center of element)

[ m ]= ρL 15 [ 4 1 2 1 4 2 2 2 16 ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaamWaaeaacaWGTbaacaGLBbGaayzxaa Gaeyypa0ZaaSaaaeaacqaHbpGCcaWGmbaabaGaaGymaiaaiwdaaaWa amWaaeaafaqabeWadaaabaGaaGinaaqaaiabgkHiTiaaigdaaeaaca aIYaaabaGaeyOeI0IaaGymaaqaaiaaisdaaeaacaaIYaaabaGaaGOm aaqaaiaaikdaaeaacaaIXaGaaGOnaaaaaiaawUfacaGLDbaaaaa@4414@  

 

 integration scheme (with 3 integration points) while for linear elements we need a quadratic scheme (2 integration points).  If the mass matrix is under-integrated it will be singular (you can check this out using the example code provided below).  Alternatively, the mass matrix can be integrated analytically (a table is given on the right)

   Assemble the contribution from each element to the global mass matrix M ab = l=1 L m ab MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa aiaabeqaaiqacaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpdaaeWbqaaiaad2gadaWgaaWcbaGaamyyaiaadkga aeqaaaqaaiaadYgacqGH9aqpcaaIXaaabaGaamitaaqdcqGHris5aa aa@3D4E@

   Modify the stiffness and mass matrices  to enforce the constraints u ¨ 1 = u ¨ * (0,t) u ¨ (L) = u ¨ * (L,t) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabmyDayaadaWaaWbaaSqabeaacaaIXa aaaOGaeyypa0JabmyDayaadaWaaWbaaSqabeaacaGGQaaaaOGaaiik aiaaicdacaGGSaGaamiDaiaacMcacaaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlqadwha gaWaamaaCaaaleqabaGaaiikaiaadYeacaGGPaaaaOGaeyypa0Jabm yDayaadaWaaWbaaSqabeaacaGGQaaaaOGaaiikaiaadYeacaGGSaGa amiDaiaacMcaaaa@5490@

   Run the Newmark time-stepping scheme outlined above.

 

 

 

 

8.2.6 Example 1D dynamic FEM code and solution

 

A simple example MAPLE code for this 1D example is given in the file FEM1D_newmark.mws.  It is set up to solve for displacements for a bar with the following properties:

1.      Length 5 and unit unit x-sect area,

2.      Mass density 100; shear modulus 50, Poisson’s ratio 0., (so the wave speed is 1);

3.       No body force,  displacement u=0 at x=0, and traction t=10 applied at x=L suddenly at t=0 and then held constant for t>0

 

As an example, we plot the variation of displacement at the end of the bar (x=L) with time.  You should be able to verify for yourself that the exact solution is a saw-tooth, amplitude 0.5 and period 20 (the time for a plane wave to make two trips from one end of the bar to the other)

 

Notice how the peaks of the saw-tooth are blunted MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  this is because the FEM interpolation functions cannot resolve the velocity discontinuity associated with a propagating plane wave.

 

The MAPLE file produces animations of the displacement and velocity field in the bar.  In the exact solution, a plane wave propagates down the bar, and repeatedly reflects off the free and fixed ends of the bar.  You can see this quite clearly in the simulation results, but the wave-front is not sharp, as it is supposed to be.  This gives rise to spurious high frequency oscillations in the velocity fields.                 

 

It is interesting to investigate the effects of the two adjustable parameters in the time integration scheme. 

 

First, we look at the solution with β 1 = β 2 =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdi2aaSbaaSqaaiaaigdaaeqaaO Gaeyypa0JaeqOSdi2aaSbaaSqaaiaaikdaaeqaaOGaeyypa0JaaGim aaaa@38BA@  (both velocity and displacement update are fully explicit).  The plots below show the predicted displacement at the end of the bar for two values of time step Δt MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeuiLdqKaamiDaaaa@333C@

       

The result with smaller time step is good, but with a larger time step the solution is oscillatory.  This is an example of a numerical instability, which is a common problem in explicit time stepping schemes.  Eventually the solution blows up completely, because the oscillations grow exponentially.  With a small enough time step an explicit scheme will usually work, but the critical time step for stability is proportional to the time required for a wave to propagate through the smallest element in the mesh.

 

With larger values of β MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8YjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdigaaa@327E@  the problem disappears.  In fact, one can show that for the equation of motion considered here, setting

β 2 β 1 β 1 1/2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdi2aaSbaaSqaaiaaikdaaeqaaO GaeyyzImRaeqOSdi2aaSbaaSqaaiaaigdaaeqaaOGaaGPaVlaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlabek7aInaaBaaaleaacaaIXaaabeaakiab gwMiZkaaigdacaGGVaGaaGOmaaaa@524B@

makes the time stepping scheme is unconditionally stable MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  no oscillations will occur even for very large time steps.

 

 

Stability does not necessarily mean accuracy, however.  Results with a fully implicit integration scheme ( β 1 = β 2 =1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdi2aaSbaaSqaaiaaigdaaeqaaO Gaeyypa0JaeqOSdi2aaSbaaSqaaiaaikdaaeqaaOGaeyypa0JaaGym aaaa@38BB@  ) and a large time step are shown on the right (set dt=0.5 in the example code).  This shows a different problem MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  energy is dissipated due to the numerical time integration scheme.   So, larger values of β MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdigaaa@326E@  buy  stability by introducing artificial damping, at the expense of a loss of accuracy.  A good compromise is to set β 1 = β 2 =1/2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdi2aaSbaaSqaaiaaigdaaeqaaO Gaeyypa0JaeqOSdi2aaSbaaSqaaiaaikdaaeqaaOGaeyypa0JaaGym aiaac+cacaaIYaaaaa@3A2A@ .

 

 

8.2.7 Lumped mass matrices

 

In the Newmark time integration scheme, accelerations are computed by solving a set of linear equations

M ab u ¨ i b (t+Δt)+ β 2 Δ t 2 K aibk u ¨ k b = K aibk [ u k b (t)+Δt u ˙ k b (t)+ Δ t 2 2 (1 β 2 ) u ¨ k b (t) ]+ F i a (t+Δt) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaGcceWG1bGbamaadaqhaaWcbaGaamyAaaqaaiaadkgaaaGccaGG OaGaamiDaiabgUcaRiabfs5aejaadshacaGGPaGaey4kaSIaeqOSdi 2aaSbaaSqaaiaaikdaaeqaaOGaeuiLdqKaamiDamaaCaaaleqabaGa aGOmaaaakiaadUeadaWgaaWcbaGaamyyaiaadMgacaWGIbGaam4Aaa qabaGcceWG1bGbamaadaqhaaWcbaGaam4AaaqaaiaadkgaaaGccqGH 9aqpcqGHsislcaWGlbWaaSbaaSqaaiaadggacaWGPbGaamOyaiaadU gaaeqaaOWaamWaaeaacaWG1bWaa0baaSqaaiaadUgaaeaacaWGIbaa aOGaaiikaiaadshacaGGPaGaey4kaSIaeuiLdqKaamiDaiqadwhaga GaamaaDaaaleaacaWGRbaabaGaamOyaaaakiaacIcacaWG0bGaaiyk aiabgUcaRmaalaaabaGaeuiLdqKaamiDamaaCaaaleqabaGaaGOmaa aaaOqaaiaaikdaaaGaaiikaiaaigdacqGHsislcqaHYoGydaWgaaWc baGaaGOmaaqabaGccaGGPaGabmyDayaadaWaa0baaSqaaiaadUgaae aacaWGIbaaaOGaaiikaiaadshacaGGPaaacaGLBbGaayzxaaGaey4k aSIaamOramaaDaaaleaacaWGPbaabaGaamyyaaaakiaacIcacaWG0b Gaey4kaSIaeuiLdqKaamiDaiaacMcaaaa@7AA8@

This is the most time-consuming part of the procedure.  But notice that if we set β 2 =0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqOSdi2aaSbaaSqaaiaaikdaaeqaaO Gaeyypa0JaaGimaaaa@3522@  and somehow find a way to make the mass matrix diagonal, then equation solution is trivial.

 

The so-called `lumped’ mass matrix is a way to achieve this.  Instead of using the correct element mass matrix

M ab = R ρ N a N b dV MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpdaWdrbqaaiabeg8aYjaad6eadaahaaWcbeqaaiaa dggaaaGccaWGobWaaWbaaSqabeaacaWGIbaaaOGaamizaiaadAfaaS qaaiaadkfaaeqaniabgUIiYdaaaa@3F30@

we replace it with a diagonal approximation.  Various procedures can be used to do this. 

 

1.      Integrate the mass matrix using integration points located at the element nodes

2.      Diagonal scaling procedure: Set M aa =c M aa MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamyyaa qabaGccqGH9aqpcaWGJbGaamytamaaBaaaleaacaWGHbGaamyyaaqa baaaaa@385B@  and M ab =0ab MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpcaaIWaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caWGHbGaeyiyIKRaamOyaaaa@40AF@ , with the constant c selected to satisfy a M aa = V e ρdV MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaabCaeaacaWGnbWaaSbaaSqaaiaadg gacaWGHbaabeaaaeaacaWGHbaabaaaniabggHiLdGccqGH9aqpdaWd rbqaaiabeg8aYjaadsgacaWGwbaaleaacaWGwbWaaSbaaWqaaiaadw gaaeqaaaWcbeqdcqGHRiI8aaaa@3F91@

3.      Row-sum method: Set M aa = b M ab MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamyyaa qabaGccqGH9aqpdaaeWbqaaiaad2eadaWgaaWcbaGaamyyaiaadkga aeqaaaqaaiaadkgaaeaaa0GaeyyeIuoaaaa@3A92@  and M ab =0ab MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaamytamaaBaaaleaacaWGHbGaamOyaa qabaGccqGH9aqpcaaIWaGaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caWGHbGaeyiyIKRaamOyaaaa@40AF@ .

 

In general, the three approaches give different answers.  For some element types, it turns out that some of these procedures give zero or negative masses, which is clearly undesirable.   Use whichever method gives the best answers.

 

Explicit time integration with a lumped mass matrix has been implemented in the 1D FEM code FEM1D_newmark.mws. To use a lumped mass matrix, find the line in the file that reads

>        lumpedmass := false:

and change it to

>        lumpedmass := true:

As an example, the figure on the right shows the predicted displacement at the end of the bar with relatively large timestep (set dt=0.1, beta1=1 and lumpedmass=true in the code to reproduce the result)

 

In general lumped mass matrices introduce some additional damping which can help stabilize the time stepping scheme, but lead to some accuracy loss.  You can see this in the result shown in the figure.  Using a smaller timestep reduces the energy loss.

 

Explicit time integration is cheap, easy to implement and is therefore a very popular technique.  Its disadvantages are that it is conditionally stable and can require very small timesteps (the critical step size for stability is proportional to the size of the smallest element in the mesh and inversely proportional to the wave speed). 

 

 

8.2.8 Example 2D and 3D dynamic linear elastic code and solution

 

It is straightforward to extend a 3D static linear elasticity code to time-domain dynamics.  As an example, the Maple code outlined in Section 8.1.14 has been extended for this purpose.

1.      The code is in a file called FEM_2Dor3D_linelast_dynamic.mws.  You will need to modify the line of the code that reads the input file to point to wherever you store the input file, as described in 8.1.14.

2.      An input file that sets up a simulation of the vibration of a cantilever beam that is subjected to a transverse force at its end is provided in the file Linear_elastic_dynamic_beam.txt.

 

The format for the input file is described in detail in Section 8.1.14.  The mesh is assumed to be at rest at time t=0, and any loads are applied as a step and held constant for t>0.  You can always extend the code to do more complex loadings if you need to.  The motion of the beam is animated at the end of the file.

 

 

 

8.2.9 Modal method of time integration

 

The finite element equations

( M ab u ¨ i b + K aibk u k b F i a )=0 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaWaaeWaaeaacaWGnbWaaSbaaSqaaiaadg gacaWGIbaabeaakiqadwhagaWaamaaDaaaleaacaWGPbaabaGaamOy aaaakiabgUcaRiaadUeadaWgaaWcbaGaamyyaiaadMgacaWGIbGaam 4AaaqabaGccaWG1bWaa0baaSqaaiaadUgaaeaacaWGIbaaaOGaeyOe I0IaamOramaaDaaaleaacaWGPbaabaGaamyyaaaaaOGaayjkaiaawM caaiabg2da9iaaicdaaaa@465B@

are a standard set of coupled, second order linear differential equations such as one would meet in analyzing forced vibrations in discrete systems.  They can therefore be solved using the usual modal techniques.  The procedure is to rearrange the system of equations to yield n decoupled 2nd order differential equations, by means of the following substitutions.  We will drop index notation in favor of matrix notation.

1.      Write the governing equation as M u ¨ +Ku=F(t) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaah2eaceWH1bGbamaacqGHRaWkcaWHlb GaaCyDaiabg2da9iaahAeacaGGOaGaamiDaiaacMcaaaa@3921@

2.      Let

q= M 1/2 u MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCyCaiabg2da9iaah2eadaahaaWcbe qaaiaaigdacaGGVaGaaGOmaaaakiaahwhaaaa@3702@               H= M 1/2 K M 1/2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCisaiabg2da9iaah2eadaahaaWcbe qaaiabgkHiTiaaigdacaGGVaGaaGOmaaaakiaahUeacaWHnbWaaWba aSqabeaacqGHsislcaaIXaGaai4laiaaikdaaaaaaa@3BB6@

and substitute for u

q ¨ +Hq= M 1/2 F MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabCyCayaadaGaey4kaSIaaCisaiaahg hacqGH9aqpcaWHnbWaaWbaaSqabeaacqGHsislcaaIXaGaai4laiaa ikdaaaGccaWHgbaaaa@3A79@

3.      Note that H is a symmetric, positive definite matrix.  Consequently, we may perform a spectral decomposition and write

H=QΛ Q T MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCisaiabg2da9iaahgfacqqHBoatca WHrbWaaWbaaSqabeaacaWGubaaaaaa@36D5@

where Q is an orthogonal matrix ( Q T Q=I MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCyuamaaCaaaleqabaGaamivaaaaki aahgfacqGH9aqpcaWHjbaaaa@356B@  ) and Λ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeu4MdWeaaa@3242@  is diagonal.  This spectral decomposition is accomplished as follows: compute the eigenvalues λ i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4UdW2aaSbaaSqaaiaadMgaaeqaaa aa@339D@  and corresponding normalized eigenvectors r (i) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCOCamaaCaaaleqabaGaaiikaiaadM gacaGGPaaaaaaa@343E@  of H (normalized to satisfy r (i) r (i) =1 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqaq=JfrVkFH e9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr0=vqpWqaaeaabiGaciaa caqabeaacmqaamaaaOqaaiaahkhadaahaaWcbeqaaiaacIcacaWGPb GaaiykaaaakiabgwSixlaahkhadaahaaWcbeqaaiaacIcacaWGPbGa aiykaaaakiabg2da9iaaigdaaaa@3B65@  ), then set

Λ=[ λ 1 λ 2 λ n ]Q=[ r 1 (1) r 1 (2) r 2 (1) r 2 (2) r n (1) r n (n) ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeu4MdWKaeyypa0ZaamWaaeaafaqabe abeaaaaaqaaiabeU7aSnaaBaaaleaacaaIXaaabeaaaOqaaaqaaaqa aaqaaaqaaiabeU7aSnaaBaaaleaacaaIYaaabeaaaOqaaaqaaaqaaa qaaaqaaiablgVipbqaaaqaaaqaaaqaaaqaaiabeU7aSnaaBaaaleaa caWGUbaabeaaaaaakiaawUfacaGLDbaacaaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaahgfacqGH9aqpdaWadaqaauaabeqaeqaaaaaabaGaamOCamaa DaaaleaacaaIXaaabaGaaiikaiaaigdacaGGPaaaaaGcbaGaamOCam aaDaaaleaacaaIXaaabaGaaiikaiaaikdacaGGPaaaaaGcbaaabaaa baGaamOCamaaDaaaleaacaaIYaaabaGaaiikaiaaigdacaGGPaaaaa GcbaGaamOCamaaDaaaleaacaaIYaaabaGaaiikaiaaikdacaGGPaaa aaGcbaaabaaabaGaeSO7I0eabaaabaGaeSy8I8eabaaabaGaamOCam aaDaaaleaacaWGUbaabaGaaiikaiaaigdacaGGPaaaaaGcbaaabaaa baGaamOCamaaDaaaleaacaWGUbaabaGaaiikaiaad6gacaGGPaaaaa aaaOGaay5waiaaw2faaaaa@7EC0@

4.      Next let w= Q T q MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaC4Daiabg2da9iaahgfadaahaaWcbe qaaiaadsfaaaGccaWHXbaaaa@35B9@  and rearrange the equation of motion a second time to get

w ¨ +Λw= Q T M 1/2 F MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGabC4DayaadaGaey4kaSIaeu4MdWKaaC 4Daiabg2da9iaahgfadaahaaWcbeqaaiaadsfaaaGccaWHnbWaaWba aSqabeaacqGHsislcaaIXaGaai4laiaaikdaaaGccaWHgbaaaa@3D13@

Since Λ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8XjY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeu4MdWeaaa@3242@  is diagonal this is now a set of n uncoupled ODEs for w.

5.      The initial conditions for w follow as

w= Q T M 1/2 u(t=0) w ˙ = Q T M 1/2 u ˙ (t=0) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaC4Daiabg2da9iaahgfadaahaaWcbe qaaiaadsfaaaGccaWHnbWaaWbaaSqabeaacaaIXaGaai4laiaaikda aaGccaWH1bGaaiikaiaadshacqGH9aqpcaaIWaGaaiykaiaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7ceWH3bGbaiaacqGH9aqpcaWHrbWaaW baaSqabeaacaWGubaaaOGaaCytamaaCaaaleqabaGaaGymaiaac+ca caaIYaaaaOGabCyDayaacaGaaiikaiaadshacqGH9aqpcaaIWaGaai ykaaaa@5D5E@

6.      Solve the decoupled ODEs for w, then  recover the displacements through

u= M 1/2 Qw MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCyDaiabg2da9iaah2eadaahaaWcbe qaaiabgkHiTiaaigdacaGGVaGaaGOmaaaakiaahgfacaWH3baaaa@38D1@

You can use any method you like to solve the ODEs for w MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  for harmonic forcing you could use the exact solution.  For general F(t) you could use Fourier transforms (the FFT works for periodic F).  You could even use the Newmark algorithm if you want, although this would be rather perverse!

 

Effectively, this procedure reduces the general forced FEM problem to an equivalent one that involves solving the equations of motion for n forced, spring-mass systems.

 

 

8.2.10 Natural frequencies and mode shapes

 

The procedure outlined in 8.2.9 can also be used to extract the natural frequencies and mode shapes for a vibrating, linear elastic solid.  To see how to do this, note that

1.      By definition, the natural frequencies and mode shapes of a continuous solid or structure are special deflections and frequencies for which the solid will vibrate harmonically.

2.      The governing equations for w can be regarded as describing the vibration of n uncoupled spring-mass systems;

3.      If we solve the problem with F=0, and excite only the ith spring-mass system (e.g. by choosing appropriate initial conditions), the solution will be a harmonic vibration at the natural frequency of the ith spring-mass system. 

4.      The natural frequencies of the n decoupled systems are therefore the natural frequencies of vibration of the structure.  The corresponding eigenvectors r (i) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCOCamaaCaaaleqabaGaaiikaiaadM gacaGGPaaaaaaa@343E@  are related to the deflections u (i) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCyDamaaCaaaleqabaGaaiikaiaadM gacaGGPaaaaaaa@3441@  associated with the ith vibration mode through u (i) = M 1/2 r (i) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCyDamaaCaaaleqabaGaaiikaiaadM gacaGGPaaaaOGaeyypa0JaaCytamaaCaaaleqabaGaeyOeI0IaaGym aiaac+cacaaIYaaaaOGaaCOCamaaCaaaleqabaGaaiikaiaadMgaca GGPaaaaaaa@3CE4@ .

 

The following procedure can be used to determine the natural frequencies and mode shapes of a linear elastic solid:

1.      Compute the finite element mass and stiffness matrices M and K

2.      Find H= M 1/2 K M 1/2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCisaiabg2da9iaah2eadaahaaWcbe qaaiabgkHiTiaaigdacaGGVaGaaGOmaaaakiaahUeacaWHnbWaaWba aSqabeaacqGHsislcaaIXaGaai4laiaaikdaaaaaaa@3BB8@

3.      Compute the eigenvalues λ i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4UdW2aaSbaaSqaaiaadMgaaeqaaa aa@339D@  and corresponding eigenvectors r (i) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCOCamaaCaaaleqabaGaaiikaiaadM gacaGGPaaaaaaa@343E@  of H

4.      The natural frequencies ω i MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeqyYdC3aaSbaaSqaaiaadMgaaeqaaa aa@33B6@ , are related to the eigenvalues of H by λ i = ω i 2 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaeq4UdW2aaSbaaSqaaiaadMgaaeqaaO Gaeyypa0JaeqyYdC3aa0baaSqaaiaadMgaaeaacaaIYaaaaaaa@3851@  .

5.      The deflections u (i) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCyDamaaCaaaleqabaGaaiikaiaadM gacaGGPaaaaaaa@3441@  associated with the ith vibration mode follow as u (i) = M 1/2 r (i) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaebbnrfifHhDYfgasaacH8rkY=vi pgYlH8Gipec8Eeeu0xXdbba9frFj0=OqFfea0dXdd9vqai=hGuQ8ku c9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=xfr=xb9adbaqaaeGaciGa biaabeqaaiqabaWaaaGcbaGaaCyDamaaCaaaleqabaGaaiikaiaadM gacaGGPaaaaOGaeyypa0JaaCytamaaCaaaleqabaGaeyOeI0IaaGym aiaac+cacaaIYaaaaOGaaCOCamaaCaaaleqabaGaaiikaiaadMgaca GGPaaaaaaa@3CE4@ . The deflections calculated in this way will have some random magnitude MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  if you like, you can normalize them appropriately.

 

The procedure is simple on paper (as always) but there are two problems in practice.  First, we need to find the square root of the mass matrix.  For a general matrix, we would have to compute a spectral decomposition of M to do this.  However, if we use a lumped mass matrix, its square root is easily found MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  just take the square root of all the diagonal elements!  Secondly, computing all the eigenvalues and eigenvectors of a general dynamical matrix H is an exceedingly time consuming process.  To make this a viable scheme, we normally extract only the lowest few eigenvalues and the corresponding eigenvectors (say 10 to 100), and discard the rest. 

 

 

8.2.11 Example 1D code with modal dynamics

 

As an example, a version of the 1D code with modal dynamics is provided in the file FEM_1D_modal.mws

 

The code is set up to plot the displacements associated with the fundamental (lowest frequency) vibration mode, as well as to calculate the displacement of the end of the bar as a function of time.  The code will handle computations with a full mass matrix as well as a lumped matrix, so you can compare the results with the two cases. 

 

The plots below show the mode shape (displacements) for the fundamental vibration mode, and the predicted motion at the end of the bar.  The code will also animate the displacement in the bar (not shown here).

 

 

        

Here is a brief comparison of the Newmark and modal time integration schemes:

1.      For linear problems the modal decomposition method usually beats direct time integration in both accuracy and speed.

2.      If you are interested in looking in detail at transient wave propagation through the solid there is not much difference between the two methods MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  in this case you need to retain a huge number of terms in the modal expansion for accurate results.

3.      For most vibration problems, however, it is generally only necessary to retain a small number of modes in the expansion, in which case the modal decomposition method will be vastly preferable to direct time integration.  Moreover, knowledge of the vibration modes can itself be valuable information. You can use also use the modal approach to handle very complex forcing, including random vibrations, directly.

4.      The main limitation of the modal decomposition approach is that it can only handle undamped, linear systems (one can introduce modal damping for vibration applications, but this does not model any real energy dissipation mechanism).  This is of course far too restrictive for any real application except the simplest possible problems in vibration analysis.

 

8.2.12 Example 2D and 3D FEM code to compute mode shapes and natural frequencies

 

It is straightforward to extend a 3D static linear elasticity code to modal dynamics. 

1.      An example program is provided in the file FEM_2Dor3D_modeshapes.mws;

2.      The code can be run with the file Linear_elastic_dynamic_beam.txt MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiFKI8=feu0dXdh9vqqj=hEeeu0xXdbba9frFj0=OqFf ea0dXdd9vqaq=JfrVkFHe9pgea0dXdar=Jb9hs0dXdbPYxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaGqaaKqzGfaeaa aaaaaaa8qacaWFtacaaa@37E6@  it computes and plots mode shapes for the beam as shown below (you need to edit the code to select which mode is displayed)

 

Mode shapes for a 1D bar clamped at one end