Dynamic Analysis of Solid Structures based on Space-Time Finite Element Analysis by Iyandri_TilukWahyono

VIEWS: 2 PAGES: 88

									           UNIVERSITY OF CINCINNATI

                                         December 18, 2008
                                    Date:___________________


   David Neil Alpert
I, _________________________________________________________,
hereby submit this work as part of the requirements for the degree of:
Master of Science
in:
Mechanical Engineering
It is entitled:
Dynamic Analysis of Solid Structures based on Space-Time
Finite Element Analysis




                    This work and its defense approved by:


                           Dr. Dong Qian
                    Chair: _______________________________
                           Dr. Randall Allemang
                           _______________________________
                           Dr. Thomas Eason, AFRL
                           _______________________________
                            Dr. Yijun Liu
                            _______________________________
                            Dr. Stephen Spottswood, AFRL
                            _______________________________
Dynamic Analysis of Solid Structures based on Space-Time Finite Element Analysis




                             A thesis submitted to the

                    Division of Research and Advanced Studies

                           of the University of Cincinnati

                             in partial fulfillment of the

                           requirements for the degree of




                             MASTER OF SCIENCE

                   in the Department of Mechanical Engineering

                           of the College of Engineering




                                December 18, 2008

                                          by

                                 David Neil Alpert

     Bachelor of Science in Mechanical Engineering with Minor in Mathematics

               University of Cincinnati, Cincinnati, Ohio, USA, 2004.

               Thesis Advisor and Committee Chair: Dr. Dong Qian

      Committee Members: Dr. Randall Allemang, Dr. Thomas Eason (AFRL),

                Dr. Yijun Liu and Dr. Stephen Spottswood (AFRL)
Abstract

       Accurate prediction of structural responses under combined, extreme environments is

traditionally analyzed using semi-discrete finite element methods. These methods pose

difficulties when simulating responses in the high frequency range and having long duration, and

capturing sharp gradients and discontinuities. These difficulties motivate the extension of the

finite element discretization scheme to the temporal domain through time-discontinuous Galerkin

(or space-time finite element) methods.

       By establishing approximations in both the spatial and temporal domains, this thesis

demonstrates the robustness of this approach in handling temporal variations in the loading

conditions. Both 1D and 2D space-time finite element codes are developed and applied to four

benchmark problems. Convergence studies with the use of different error estimators are

conducted. The developed space-time formulation is shown to be both robust and accurate.

Comparing with the semi-discrete scheme, it has unique advantages in tracing structural

responses under a variety of time-dependent loading conditions.




                                               iii
iv
Acknowledgements

       I would like to thank my advisor, Dr. Dong Qian, for his help and guidance in completing

this thesis. I also want to thank Dr. Randall Allemang for all of his support and guidance over the

past few years. I want to thank Dr. Yijun Liu, Dr. Thomas Eason and Dr. Stephen Spottswood

for serving on my committee and for their help and input towards making this thesis a success. I

greatly appreciate the assistance all of the above have provided in helping to open up

opportunities for me to do research at the Air Force Research Laboratory at Wright Patterson Air

Force Base over the past couple of years.

       I also want to thank my parents, Bob and Merle Alpert, along with the rest of my family,

Theta Tau brothers, and friends for their continuous support and encouragement throughout my

studies. Without their help, I never would have made it this far.




                                                 v
Table of Contents
Abstract ........................................................................................................................................................ iii
Acknowledgements ....................................................................................................................................... v
List of Figures ............................................................................................................................................ viii
Chapter 1.             Motivation, Introduction, and Literature Review ................................................................. 1
Chapter 2.             Space-Time Formulation for Linear Elastodynamics ........................................................... 6
   2.1.        Linear Elastodynamics .................................................................................................................. 6
   2.2.        Space-Time Finite Element Formulation ...................................................................................... 7
   2.3.        A Time-Discontinuous Galerkin Approach ................................................................................ 11
Chapter 3.             Discretization ...................................................................................................................... 13
   3.1.        Discretization of BDG ................................................................................................................... 16
   3.2.        Discretization of LDG ................................................................................................................... 20
   3.3.        Combination of Discretized BDG and LDG ................................................................................... 24
   3.4.        Temporal Mesh Limitations ........................................................................................................ 29
Chapter 4.             Mesh Refinement and Convergence Rate ........................................................................... 30
Chapter 5.             Analytical Space-Time Matrix Formulation ....................................................................... 33
   5.1.        Quadratic Lagrangian Temporal Shape Functions ...................................................................... 33
Chapter 6.             Implementation of Space-Time Formulation ...................................................................... 37
   6.1.        Axially Loaded Bar ..................................................................................................................... 37
      6.1.1. General Analytical Solution – Modal Superposition Method .................................................... 37
      6.1.2. General Analytical Solution – Acoustic Wave Propagation Method ......................................... 41
      6.1.3. Numerical Solution ..................................................................................................................... 44
      6.1.4. Convergence Results .................................................................................................................. 46
      6.1.5. Numerical Results....................................................................................................................... 49
   6.2.        Transversely Loaded Straight Beam ........................................................................................... 52
      6.2.1. Analytical Solution ..................................................................................................................... 53
      6.2.2. Numerical Solution ..................................................................................................................... 57
      6.2.3. Convergence Results .................................................................................................................. 59
      6.2.4. Numerical Results....................................................................................................................... 61
   6.3.        Transversely Loaded Thick Curved Beam .................................................................................. 63
   6.4.        Transversely Loaded Slender Curved Beam ............................................................................... 67
      6.4.1. Convergence Results .................................................................................................................. 69


                                                                                 vi
      6.4.2. Numerical Results....................................................................................................................... 70
Chapter 7.            Conclusions and Future Work............................................................................................. 74
Bibliography ............................................................................................................................................... 76




                                                                             vii
List of Figures
Figure 1. Space-Time Slabs Emphasizing Temporal Discontinuities ........................................................... 8
Figure 2. Normalized Quadratic Temporal Shape Functions...................................................................... 34
Figure 3. Bar Subject to Dynamic Axial Load, Analytical Example .......................................................... 38
Figure 4. Section of Bar Subject to Dynamic Axial Load, Analytical Example......................................... 39
Figure 5. Bar Subject to Dynamic Axial Load, Numerical Example .......................................................... 44
Figure 6. Cross Section A-A of Bar Subject to Dynamic Axial Load, Numerical Example ...................... 45
Figure 7. 1D Bar Spatial Mesh ................................................................................................................... 46
Figure 8. Strain Energy Norm at 0.00405 seconds ..................................................................................... 47
Figure 9. L2 Residual Displacement Norm ................................................................................................. 48
Figure 10. Total Energy Norm .................................................................................................................... 49
Figure 11. Total Energy Across Temporal Jump ........................................................................................ 49
Figure 12. 1D Bar Axial Displacement (Δt = 5.29×10-5 s), in m ................................................................ 50
Figure 13. 1D Bar Axial Displacement (Δt = 2.65×10-4 s), in m ................................................................ 50
Figure 14. Axial Displacement of Free End of Bar, in m ........................................................................... 51
Figure 15. Axial Stress at Midpoint of Bar, in Pa ....................................................................................... 51
Figure 16. Axial Stress at Midpoint of Bar, in Pa ....................................................................................... 52
Figure 17. Beam Subject to Dynamic Transverse Load, Analytical Example ............................................ 53
Figure 18. Section of Beam in Figure 17 .................................................................................................... 54
Figure 19. Beam Subject to Dynamic Transverse Load, Numerical Example ........................................... 57
Figure 20. Cross Section A-A of Beam in Figure 19 .................................................................................. 57
Figure 21. 2D Beam Spatial Mesh .............................................................................................................. 59
Figure 22. Zoomed In 2D Beam Spatial Mesh ........................................................................................... 59
Figure 23. L2 Residual Displacement Norm ............................................................................................... 60
Figure 24. Total Energy Norm .................................................................................................................... 61
Figure 25. Mesh of Transverse Displacement of Straight Beam, in m ....................................................... 62
Figure 26. Transverse Displacement of Midpoint of Straight Beam, in m ................................................. 63
Figure 27. Thick Curved Beam Subject to Dynamic Transverse Load, Numerical Example .................... 63
Figure 28. Cross Section A-A of Thick Curved Beam in Figure 27 ........................................................... 64
Figure 29. Spatial Mesh of Thick Curved Beam......................................................................................... 65
Figure 30. Zoomed In Spatial Mesh of Thick Curved Beam ...................................................................... 65
Figure 31. Mesh of X Displacement of Thick Curved Beam, in m ............................................................ 66
Figure 32. Mesh of Y Displacement of Thick Curved Beam, in m ............................................................ 66
Figure 33. Work Done to Thick Curved Beam, in J ................................................................................... 67
Figure 34. Slender Curved Beam Subject to Dynamic Transverse Load, Numerical Example.................. 68
Figure 35. Cross Section A-A of Slender Curved Beam in Figure 34 ........................................................ 68
Figure 36. Spatial Mesh of Slender Curved Beam ...................................................................................... 69
Figure 37. Zoomed in Spatial Mesh of Slender Curved Beam ................................................................... 69
Figure 38. L2 Residual Displacement Norm ............................................................................................... 70
Figure 39. Mesh of X Displacement of Slender Curved Beam, in m ......................................................... 71
Figure 40. Mesh of Y Displacement of Slender Curved Beam, in m ......................................................... 71
Figure 41. Stress in the XX-Direction, in Pa .............................................................................................. 72
Figure 42. Stress in the YY-Direction, in Pa .............................................................................................. 72

                                                                       viii
Figure 43. Stress in the XY-Direction, in Pa .............................................................................................. 73




                                                                     ix
Chapter 1.             Motivation, Introduction, and Literature Review

       Accurate prediction of structural responses under combined, extreme environments is of

critical interest to many industrial applications. In many of these applications, a wide range of

spatial and temporal scales are involved. In the traditional analysis of structural response

problems, such as those based on finite element methods, time dependent problems are generally

solved using the so-called semi-discrete scheme, i.e., the spatial domain is discretized using

finite elements and the temporal domain is resolved using a finite difference approach.

Consequently, although the spatial feature of the system can be solved at any time step, the time

history of the system has to be resolved in a sequential manner. Such an approach works well

when the response of the system is generally in the low frequency range and short term, but

poses stringent limits when simulating the responses in the high frequency range and for long

duration. The time step that can be taken, if an explicit scheme is chosen, is limited by the

critical time step that is related to both the numerical discretization and the physical wave speed.

Another disadvantage of semi-discrete analysis is that sharp gradients and discontinuities are

often difficult to accurately capture (Hulbert & Hughes, 1990, p. 327). Such sharp gradients and

discontinuities are frequently encountered in simulation of shocks and material failure. As a

result, simulating long term responses of engineering systems is still a challenging issue with the

use of semi-discrete schemes.

       Motivated by the difficulties associated with the semi-discrete schemes, it has been

suggested that the finite element discretization scheme can also be extended to the temporal

domain. The resulting formulation, which is termed “space-time formulation,” was first put forth

by Argyris and Scharpf (1969), Fried (1969), and Oden (1969). The initial research assumed that

the unknown fields of interest, such as displacement and velocity, were continuous in time giving


                                               1 / 79
a time-continuous Galerkin (TCG) formulation. Some of these methods have been shown to be

conditionally stable. More recently, the space-time approach has been extended by incorporating

discontinuities in the temporal domain. This formulation was originally developed for first order

hyperbolic equations and later extended to the fields of fluid dynamics and heat conduction. It

has been shown that this time-discontinuous Galerkin formulation is both stable and higher-order

accurate (Hulbert & Hughes, 1990, pp. 330-331).

       Time-discontinuous Galerkin (TDG) methods were extended to problems involving

second order hyperbolic equations like elastodynamics and structural dynamics by Hulbert and

Hughes (Hulbert & Hughes, 1990, p. 327; Hulbert, 1992, p. 307) and to acoustics by Thompson

and Pinsky (1996, p. 3297). In the work of Hulbert and Hughes, two types of general space-time

finite element formulations were derived. The first is termed the single field formulation. In this

formulation, the displacements alone are the basic unknowns. The second is called a two field

formulation and it uses both displacements and velocities as the basic unknowns (Li & Wiberg,

1998, p. 211). Johnson went on to prove a priori and a posteriori error estimates for solving

linear second order hyperbolic equations using TDG methods (1993). Additional studies have

focused on applying adaptive (Cavin, Gravouil, Lubrecht, & Combescure, 2005), stabilization

(Oñate & Manzán, 1999), and multiscale methods (Hughes & Stewart, 1996) to the TDG

approach.

       Costanzo and Huang (2005) expanded on Hulbert and Hughes’s (1990) formulations to a

more general case allowing unstructured finite element grids and provided proofs showing the

unconditional stability of the new methods. Huang and Costanzo’s work built on this research

demonstrating the use of space-time finite elements to solve dynamic solid/solid phase

transitions as well as dynamic fracture (Huang & Costanzo, 2002; Huang & Costanzo, 2004).



                                              2 / 79
Idesman developed an implicit approach for structured and unstructured meshes using both TCG

and TDG methods showing an increase in accuracy (Idesman A. V., 2007a; Idesman A. V.,

2007b; Idesman, Schmidt, & Sierakowski, 2008). Recently, Chen et al have derived a single field

TDG method using velocity as the basic unknowns (Chen, Steeb, & Diebels, 2008).

       Besides elastodynamics, the research has also gone into numerous other fields of study

such as viscoelasticity (Shaw & Whiteman, 2000), rotating systems (Gudla & Ganguli, 2006),

and contact/impact problems (Karaoğlan & Noor, 1997; Pattillo & Tortorelli, 2004). In the field

of fluid mechanics, space-time finite elements have been applied to the convection-diffusion

(Varoḡlu & Finn, 1980) (Varoḡlu & Finn, 1978), Burgers (Froncioni, Labbe, Garon, &

Camarero, 1997), and Navier-Stokes (Masud & Hughes, 1997) equations and to fluid flow

(Surana, Allu, Tenpas, & Reddy, 2007; Feng & Perić, 2000; N'dri, Garon, & Fortin, 2001). The

formulation of a space-time least-squares approach has been developed for fluid dynamics to

achieve better stability and reduce the effect of spurious modes (Hughes, Franca, & Hulbert,

1989; Shakib & Hughes, 1991; Nguyen & Reynen, 1984). TDG methods were also extended to

dynamics and wave propagation in nonlinear solids and porous media (Li, Yao, & Lewis, 2003;

Chen, Steeb, & Diebels, 2006). Chien, et al have applied the TDG approach to boundary element

methods (2003). Chessa and Belytschko explored space-time finite elements using enriched

finite element methods (2004; 2006). Kurdi and Beran (2008) developed space-time finite

element method by applying a spectral element method for the discretization of the temporal

domain.

       The research presented in this thesis is largely motivated by the idea of establishing

approximations for both the spatial and temporal domains in the space-time formulation. While a

number of studies have proven the convergence properties of the method in both the spatial and



                                            3 / 79
temporal domains, the loading conditions that are employed in many of these cases are mostly

linear or constant functions of time. The capability of the method of treating loading conditions

with strong temporal variations has not been fully demonstrated. Therefore, one of the main

objectives of this thesis is to explore the robustness of the space-time approach in handling

dynamic loading conditions. By following largely the theoretical background laid in the work of

Hulbert and Hughes (1990), both 1D and 2D space-time finite element codes are developed in

this thesis. Furthermore, a number of benchmark problems are developed and it is shown that the

space-time formulation has unique advantages in tracing dynamic structural responses.

       The rest of the thesis is organized as follows: In Chapter 2 the mathematical background

behind elastodynamics and the space-time finite methods using a single field formulation is

discussed. Chapter 3 introduces the discretization schemes used to approximate the spatial and

temporal domains. Limitations to the temporal mesh are also briefly discussed in this chapter.

Chapter 4 focuses on the refinement of the discretization meshes and the convergence of the

finite element results towards an analytical solution. In Chapter 5 the space-time formulation

using quadratic Lagrangian shape functions for the approximations in the temporal domain is

analytically derived. In Chapter 6, the codes developed for this thesis using the space-time

formulation are presented. The code is applied to a 1D axially loaded bar problem in Section 6.1

and convergence studies are performed using an analytical solution for the axial displacement. In

Section 6.2 the code is applied to a 2D transversely loaded beam problem with a time varying

forcing function. An analytical solution for the transverse displacement is derived for this

problem and used for comparison with the finite element solution. Convergence studies are also

performed on this problem. In Section 6.3 and Section 6.4 the code is applied to the problem of a

curved beam with a temporally periodic forcing function in the vertical direction. In Chapter 7



                                             4 / 79
the thesis is concluded and the work that possibly could be extended in the future in this field is

discussed.




                                              5 / 79
Chapter 2.             Space-Time Formulation for Linear Elastodynamics

2.1.   Linear Elastodynamics

       Consider a spatial region Ω∈        d
                                                bounded by Γ , where d is the number of spatial

dimensions. The boundary can be divided up into the essential and traction boundary conditions,

Γg and Γh respectively. These two subregions do not overlap, thus

                                               Γ = Γ g ∪ Γh                                     [2-1]

                                               Γ g ∩ Γh = ∅ .                                   [2-2]

Let the displacement vector be u ( x , t ) , where x ∈Ω and t ∈[ 0, T ] , with time interval of length

T > 0 . The stress σ is a function of the gradient of the displacement u according to the

generalized Hooke’s law:

                                          σ (∇u ) = C ⋅ ∇u                                      [2-3]


in which C is the elasticity tensor and ∇ is the gradient operator. The equation of motion of the


elastodynamics problem is

                                        ρ u = ∇ ⋅σ (∇ u ) + f     on Q ≡ Ω× ]0, T [ .           [2-4]

The essential boundary conditions are prescribed values of the unknown field(s) (Cook, Malkus,

Plesha, & Witt, 2001, p. 137). In this case, the essential (displacement) boundary condition is

                                                  u=g             on ϒ g ≡ Γ g × ]0, T [ .      [2-5]

The traction boundary conditions are prescribed values of higher derivatives of the unknown

field quantities than are usually used as nodal degrees of freedom (Cook, Malkus, Plesha, &

Witt, 2001, p. 137). The traction boundary condition of consists of prescribed forces, moments,

and pressures. The traction boundary condition in the present case is

                                                  6 / 79
                                              n ⋅σ (∇u ) = h         on ϒh ≡ Γh × ]0, T [ .        [2-6]

The higher derivative described above is shown in Eq. [2-6] by the dependency of the stress

function on the gradient of the displacement field. The initial conditions are

                                             u ( x ,0 ) = u0 ( x )   for x ∈Ω                      [2-7]

                                             u ( x ,0 ) = v0 ( x )   for x ∈Ω                      [2-8]

where ρ = ρ ( x ) > 0 is volumetric mass density, a superposed dot indicates partial

differentiation with respect to time, f is the body force (as a function of time and space) per unit

volume,   g     is the prescribed boundary condition, h is the prescribed traction (pressure) on the

boundary (as a function of time and space) with units of force per unit length, u0 is the initial

displacement, v0 is the initial velocity, n is the unit outward normal to Γ .

2.2.    Space-Time Finite Element Formulation

        Let the temporal domain, I = ]0, T [ , be split up into N segments, such that

0 = t0 < t1 <     < tN = T . Considering a particular segment, let I n = ]tn−1, tn [ and Δtn = tn − tn−1

where 1 ≤ n ≤ N as shown in Figure 1.




                                                   7 / 79
                                 t


                         −
                        tN

                          +
                        t N −1

                          −
                        t N −1
                         +
                        tn +1

                         −
                        tn +1

                         +
                        tn

                         −
                        tn

                         +
                        tn −1

                         −
                        tn −1

                         +
                        tn − 2

                         −
                        tn − 2
                        t1+

                        t1−

                         +
                        t0
                                                                 x

             Figure 1. Space-Time Slabs Emphasizing Temporal Discontinuities



                                        Qn = Ω× I n                                   [2-9]

                                         ϒn = Γ× I n                                 [2-10]

                                      (ϒ )
                                         g n   = Γg × In                             [2-11]


                                      ( ϒ h )n = Γ h × I n                           [2-12]

                                      th
From Eq. [2-9], Qn is considered the n space-time slab. Let space-time slab Qn be broken up

into ( nel )n space-time elements, with the domain (interior) of the e element defined as
                                                                      th




                                            8 / 79
Qn ⊂ Qn and its boundary is ϒe . The domain and boundary of the interior of the slab are defined
 e
                             n


as

                                                                     ( nel )n
                                                      Σ
                                                     Qn =                ∪Q
                                                                         e =1
                                                                                        e
                                                                                        n               (element interiors)    [2-13]


                                                              ( nel )n
                                                ϒ =  Σ
                                                     n            ∪ϒ
                                                                  e =1
                                                                                e
                                                                                n   − ϒn                (interior boundary).   [2-14]


To proceed with the derivation, a number of bilinear operators are introduced, given as

                                      (w ,u ) = ∫
                                            h        h
                                                            Ω             Ω
                                                                                wh ⋅ u h d Ω                                   [2-15]


                               a ( wh , u h ) = ∫ ( ∇wh ) ⋅ σ ( ∇u h ) d Ω
                                                Ω             Ω
                                                                                                                               [2-16]
                                                     = ∫ ( ∇wh ) ⋅ σ d Ω
                                                              Ω



                                      (w    h
                                                ,uh )
                                                           Qn
                                                                  = ∫ w h ⋅ u h dQ
                                                                          Qn
                                                                                                                               [2-17]


                              a ( wh , u h )         =∫  ( ∇w ) ⋅ σ ( ∇u ) dQ       h               h
                                                Qn            Qn
                                                                                                                               [2-18]
                                                     = ∫ ( ∇w ) ⋅ σ dQ              h
                                                              Qn



                                      (w ,u )
                                            h        h
                                                            Σ
                                                           Qn
                                                                  = ∫ Σ w h ⋅ u h dQ
                                                                          Qn
                                                                                                                               [2-19]


                                      (w ,u )
                                            h        h
                                                           ϒΣ
                                                            n
                                                                  = ∫ Σ wh ⋅ u hd ϒ
                                                                          ϒn
                                                                                                                               [2-20]


                                   ( w , u )(
                                       h        h
                                                         ϒ h )n
                                                                  =∫
                                                                          ( ϒ h )n
                                                                                        wh ⋅ u h d ϒ                           [2-21]


where

                                       ∫   Qn
                                                     dQ = ∫
                                                                         In   ∫ Ω
                                                                                            d Ωdt                              [2-22]

                                                                     ( nel )n

                                       ∫    Σ
                                           Qn
                                                     dQ =                ∑∫
                                                                         e =1
                                                                                     e
                                                                                    Qn
                                                                                             dQ                                [2-23]




                                                                  9 / 79
                                        ∫(   ϒ h )n
                                                      dϒ = ∫
                                                                 In   ∫
                                                                      Γh
                                                                           d Γdt             [2-24]


and a ( ⋅, ⋅)Ω is the strain-energy inner product.

        A typical space-time element encompasses domain Qe with boundary ϒ e . The unit

vector n is normal to the boundary ϒ e ∩ {t} and pointed outward in the spatial hyperplane

Q e ∩ {t} . If Q equals the product of the spatial domain Ω and a time interval, then n is the
                e                                          e




usual spatial unit outward normal vector to the spatial boundary Γ e (Hulbert & Hughes, 1990, p.

335).

        For two adjacent space-time elements, designated by + and − , let n + and n − be the

spatial unit outward normal vectors along the boundary between them. The jump operator           is

defined by

                                             w ( tn ) = w ( tn ) − w ( tn )
                                                             +          −
                                                                                             [2-25]


                                         w( x ) = w( x + ) − w( x − ) ,                      [2-26]

where

                                              w ( tn ) = lim w ( tn ± ε )
                                                   ±
                                                            ±
                                                                                             [2-27]
                                                          ε →0


                                         w ( x ± ) = lim± w ( x + ε n )                      [2-28]
                                                         ε →0



                                                      n = n+ = −n− ;                         [2-29]

w ( x , t ) is discontinuous at time tn and point x in Eq. [2-25] and [2-26], respectively; ε is an

infinitesimal perturbation; and n is the normal to the surface of discontinuity. Then

                           σ ( ∇w, x ) ⋅ n = σ ( ∇w, x + ) ⋅ n − σ ( ∇w, x − ) ⋅ n
                                                                                         .   [2-30]
                           σ ( ∇w, x ) ⋅ n = σ ( ∇w, x + ) ⋅ n + + σ ( ∇w, x − ) ⋅ n −

                                                         10 / 79
2.3.    A Time-Discontinuous Galerkin Approach

        Using an updated Lagrangian time-discontinuous Galerkin formulation, the derivation of

the weak form of the differential equations of motion is based on the references (Hulbert &

Hughes, 1990) and (Belytschko, Liu, & Moran, 2000).

        First of all, the trial functions u ( x , t ) and test functions δ u ( x , t ) are assumed to be C 0

continuous and satisfy the essential boundary conditions.

                         u ( x, t ) ∈U        U = {u ( x, t ) u ∈ C 0 ( x ) , u = g on Γ g }                  [2-31]


                         δ u ( x , t ) ∈U0    U0 = {δ u ( x, t ) δ u ∈ C 0 ( x ) , δ u = 0 on Γ g }           [2-32]

where C 0 gives the space of continuous functions. The stress σ ( x , t ) is assumed to be a C −1

function in space. The conservation of momentum is denoted by

                                                   ∇⋅σ + f = ρu .                                             [2-33]

        The weak form is developed by pre-multiplying the momentum equation in Eq. [2-33] by

the test function δ u ( x , t ) and integrating over the current domain. The integration is performed

over the time period ⎡tn −1 , tn ⎤ .
                     ⎣
                       −       −
                                 ⎦

                         δ u ⋅ ( ρ u − ∇ ⋅ σ − f ) dQ + ∫               δ u ⋅ ( ρ u − ∇ ⋅ σ − f ) d Ωdt = 0
                                                             +
                                                            tn −1
                   ∫Qn                                    −
                                                         tn −1      ∫
                                                                    Ω
                                                                                                              [2-34]

        The momentum equation here is broken up into a discontinuous portion that is integrated

over the time jump between the immediately previous and current time slabs where t ∈ ⎡tn −1 , tn −1 ⎤
                                                                                     ⎣
                                                                                       −       +
                                                                                                    ⎦

and a continuous portion integrated over the current time slab where t ∈ ⎡tn−1, tn ⎤ . Based on the
                                                                         ⎣
                                                                           +     −
                                                                                   ⎦

work of Hughes and Hulbert (1990), the final weak form is given as




                                                        11 / 79
                       (δ u , ρ u )Q                                  +
                                                                         (
                                             + a (δ u , u )Q + δ u ( tn −1 ) , ρ u ( tn −1 )
                                                                                      +
                                                                                                                   )   Ω
                                                                                                                                  (
                                                                                                                           + a δ u ( tn −1 ) , u ( tn −1 )
                                                                                                                                      +             +
                                                                                                                                                                          )   Ω
                                                                                                                                                                                  =
                                                                                                                                                                                      .               [2-35]
                                         n                       n



                       ( δ u , f )Q + ( δ u , h ) ( ϒ
                                     n
                                                            )
                                                           h n
                                                                     (
                                                                 + δ u ( tn −1 ) , ρ u ( tn −1 )
                                                                          +               −
                                                                                                              )    Ω
                                                                                                                              (
                                                                                                                       + a δ u ( tn −1 ) , u ( tn −1 )
                                                                                                                                  +             −
                                                                                                                                                                      )   Ω



Equation [2-35] can be equivalently expressed in bilinear form

                                                                     BDG (δ u , u ) n = LDG (δ u ) n                                                                                                  [2-36]

for n = 1, 2,… , N , where

              BDG (δ u, u )n = (δ u, ρ u )Q + a (δ u, u )Q + δ u ( tn−1 ) , ρ u ( tn −1 )
                                                       n
                                                                    +              +
                                                                                     n
                                                                                                (                                         )   Ω
                                                                                                                                                         (
                                                                                                                                                  + a δ u ( tn −1 ) , u ( tn−1 )
                                                                                                                                                             +             +
                                                                                                                                                                                              )   Ω
                                                                                                                                                                                                      [2-37]


               LDG (δ u ) n = (δ u , f ) + (δ u , h )( ϒ
                                                  Qn                           )
                                                                              h n
                                                                                           (
                                                                                    + δ u ( tn −1 ) , ρ u ( tn −1 )
                                                                                             +               −
                                                                                                                                  )   Ω
                                                                                                                                                     (
                                                                                                                                          + a δ u ( tn −1 ) , u ( tn −1 )
                                                                                                                                                     +             −
                                                                                                                                                                                          )   Ω
                                                                                                                                                                                                  . [2-38]

In integral form, this is


                ∫ (δ u ⋅ ρ u + ∇ (δ u ) ⋅σ ) dQ + ∫ (δ u ( t ) ⋅ ρ u ( t ) + ∇δ u ( t ) ⋅σ ( t ) ) d Ω =
                                                                                               +              +                                +               +
                 Qn                                                       Ω                    n−1            n−1                              n−1             n−1
                                                                                                              .                                                                                       [2-39]
                ∫( )
                  ϒh
                     (δ u ⋅ h ) d ϒ + ∫ (δ u ⋅ f ) dQ + ∫ (δ u ( t ) ⋅ ρ u ( t ) + ∇δ u ( t ) ⋅ σ ( t ) ) d Ω
                                                  Qn                                 Ω
                                                                                                     +
                                                                                                     n −1
                                                                                                                           −
                                                                                                                           n−1
                                                                                                                                                         +
                                                                                                                                                         n−1
                                                                                                                                                                          −
                                                                                                                                                                          n −1
                       n




Substituting for the stress using Hooke’s Law in Eq. [2-3]


∫Qn    (δ u ⋅ ρ u + ∇ (δ u ) ⋅ C ⋅ ∇ ( u )) dQ + ∫ (δ u (t ) ⋅ ρ u (t ) + ∇δ u (t ) ⋅ C ⋅ ∇ ( u (t ))) d Ω = .[2-40]
                                                                      Ω
                                                                                         +
                                                                                         n −1
                                                                                                            +
                                                                                                            n −1
                                                                                                                                              +
                                                                                                                                              n −1
                                                                                                                                                                          +
                                                                                                                                                                          n −1



∫(       (δ u ⋅ h ) d ϒ + ∫ (δ u ⋅ f ) dQ + ∫ (δ u ( t ) ⋅ ρ u ( t ) + ∇δ u ( t ) ⋅ C ⋅ ∇ ( u ( t ) ) ) d Ω
     ϒ h )n                     Qn                               Ω
                                                                                    +
                                                                                    n −1
                                                                                                       −
                                                                                                       n −1
                                                                                                                                      +
                                                                                                                                      n −1
                                                                                                                                                                     −
                                                                                                                                                                     n −1



Equation [2-40] represents the final weak form for elastodynamics in terms of the space-time

formulation.




                                                                                           12 / 79
Chapter 3.              Discretization

         For simplicity, a multiplicative form of the space-time approximation is developed. In

other words, the approximations in the spatial and temporal domains are established

independently. The spatial dimension is approximated by a discretization using shape functions

                  th
[ Ns ]   for the ns    space-time element. The temporal dimension is approximated by a


discretization using shape functions [ Nt ] = ⎡ Nt1
                                                                                                    th
                                              ⎣                       Nti            Ntk ⎤ for the nt space-time slab
                                                                                         ⎦

with k nodes. Combining both, the overall shape functions for the space-time slab are

                             N ( x, t ) = ⎡ Nt1 N s
                                          ⎣                   Nti N s             N tk N s ⎤ .
                                                                                           ⎦                    [3-1]

The approximated displacement at any point ( x , t ) is given by

                         u = N ( x, t ) d = ⎡ Nt1 N s
                                            ⎣                      Nti N s           N tk N s ⎤ d
                                                                                              ⎦                 [3-2]

where d is the nodal displacement vector of the space-time slab. The velocity and acceleration

approximations at any point are given by

                                v =u
                                  = N ( x, t ) d                                                                [3-3]
                                  = ⎡ N t1 N s
                                    ⎣                   N ti N s             N tk N s ⎤ d
                                                                                      ⎦

                               a=u
                                  = N ( x, t ) d                                            .                   [3-4]
                                  = ⎡ N t1 N s
                                    ⎣                   N ti N s             N tk N s ⎤ d
                                                                                      ⎦

The test function at any point is given by

                             δ u = δ ( N ( x, t ) d )
                                  = N ( x, t ) δ d                                              .               [3-5]
                                  = ⎡ N t1 N s
                                    ⎣                   N ti N s             N tk N s ⎤ δ d
                                                                                      ⎦


                                                      13 / 79
        For the discontinuous portion, the displacements and velocities from the previous space-

                                                                                     −
time slab are used as initial conditions for the current space-time slab. For time t n −1 , the nodal

                                                                                          +
displacement vector from the previous space-time slab is used, namely d n −1 . For time t n −1 , the

nodal displacement vector for the current space-time slab ( d n = d ) is used. The test function

                                         +
only for the current time step at time t n −1 is used.

                   u ( tn−1 ) = N ( x, tn−1 ) d n−1
                        −               −

                                                                                                                          [3-6]
                             = ⎡ Nt1 ( tn−1 ) N s
                               ⎣
                                        −
                                                                Nti ( tn−1 ) N s
                                                                       −
                                                                                             Ntk ( tn−1 ) N s ⎤ d n−1
                                                                                                    −
                                                                                                              ⎦

                    u ( tn−1 ) = N ( x, tn−1 ) d
                         +               +

                                                                                                                          [3-7]
                              = ⎡ Nt1 ( tn−1 ) N s
                                ⎣
                                         +
                                                                  Nti ( tn−1 ) N s
                                                                         +
                                                                                               Ntk ( tn−1 ) N s ⎤ d
                                                                                                      +
                                                                                                                ⎦

                         +
                                   (
                  δ u ( tn−1 ) = δ N ( x , tn−1 ) d
                                            +
                                                          )
                              = N ( x , tn−1 ) δ d
                                         +
                                                                                                                          [3-8]

                              = ⎡ N t1 ( tn−1 ) N s
                                ⎣
                                          +
                                                                  N ti ( tn−1 ) N s
                                                                          +
                                                                                               N tk ( tn−1 ) N s ⎤ δ d
                                                                                                       +
                                                                                                                 ⎦

The rate-of-deformation can be expressed in terms of the shape functions by

                                                  ∇u ( x , t ) = ∇ N ( x , t ) d
                                                                 = N,x ( x , t ) d                                        [3-9]
                                                                 = B ( x, t ) d

where

                        N,x ( x, t ) = ⎡ Nt1 N s , x
                                       ⎣   (          )               ( N N ),
                                                                          ti     s    x         ( N N ), ⎤
                                                                                                     tk  ⎦ s   x
                                                                                                                         [3-10]
                                       = ⎡ Nt1 N s , x
                                         ⎣                         Nti N s , x            N tk N s , x ⎤
                                                                                                       ⎦

and the strain displacement matrix is

                                                              Bs = N s , x .                                             [3-11]

Similarly,

                                                                14 / 79
                                                         ∇u ( x, t ) = ∇N ( x, t ) d
                                                                        = N,x ( x, t ) d                                   [3-12]
                                                                        = B ( x, t ) d

                                                      ∇δ u ( x , t ) = ∇ N ( x , t ) δ d
                                                                        = N,x ( x , t ) δ d .                              [3-13]
                                                                        = B ( x, t ) δ d

The weak form in Eq. [2-40] is approximated as


         ∫Qn( N ( x, t )δ d ⋅ ρ N ( x, t ) d + ∇ ( N ( x, t )δ d ) ⋅ C ⋅∇ ( N ( x, t ) d )) dQ
         + ∫ ( N ( x, t ) δ d ⋅ ρ N ( x, t ) d + ∇N ( x, t ) δ d ⋅ C ⋅∇ ( N ( x, t ) d ) ) d Ω
                                  +
                                  n −1
                                                       +
                                                       n −1
                                                                               +
                                                                               n −1
                                                                                                    +
                                                                                                    n −1
             Ω
                                                                                                                         . [3-14]
         =∫
               ( ϒ h )n
                          ( N ( x, t ) δ d ⋅ h ) d ϒ + ∫ Qn
                                                              (                       )
                                                                  N ( x, t ) δ d ⋅ f dQ


             Ω   (                                                                              (
         + ∫ N ( x, tn−1 ) δ d ⋅ ρ N ( x, tn−1 ) dn −1 + ∇N ( x, tn−1 ) δ d ⋅ C ⋅∇ N ( x, tn−1 ) d n−1 d Ω
                     +                     −                      +                        −
                                                                                                              ))
Using the substitutions from Eq. [3-9] through [3-13] and using the associative and commutative

properties


               ∫ Qn          (
                      δ d N ( x, t ) ⋅ ρ N ( x, t ) + B ( x, t ) ⋅ C ⋅ B ( x, t ) ddQ     )
                     Ω                   (
               + ∫ δ d ( tn−1 ) N ( x, tn−1 ) ⋅ ρ N ( x, tn−1 ) + B ( x, tn−1 ) ⋅ C ⋅ B ( x, tn−1 ) d ( tn−1 ) d Ω
                          +             +                 +               +                   +
                                                                                                     )   +

                                                                                                                     .     [3-15]
               =∫
                      ( ϒh )n
                                ( N ( x, t ) δ d ⋅ h ) d ϒ + ∫ ( N ( x, t ) δ d ⋅ f ) dQ
                                                                  Qn



                     Ω
                          +             +
                                         (                                                           )
               + ∫ δ d ( tn−1 ) N ( x, tn−1 ) ⋅ ρ N ( x, tn−1 ) + B ( x, tn−1 ) ⋅ C ⋅ B ( x, tn−1 ) d ( tn−1 ) d Ω
                                                          −               +                   −          −




Using the notation in Eq. [2-22]




                                                                       15 / 79
                ∫ ∫
                 In       Ω            (                                                     )
                              δ d N ( x, t ) ⋅ ρ N ( x, t ) + B ( x, t ) ⋅ C ⋅ B ( x, t ) dd Ωdt


                      Ω            (
                + ∫ δ d N ( x, tn−1 ) ⋅ ρ N ( x, tn−1 ) + B ( x, tn −1 ) ⋅ C ⋅ B ( x, tn−1 ) dd Ω
                                +                 +               +                    +
                                                                                                     )   .   [3-16]
                =∫
                      ( ϒ h )n
                                   ( N ( x, t ) δ d ⋅ h ) d ϒ + ∫
                                                                Qn
                                                                     (                   )
                                                                         N ( x, t ) δ d ⋅ f dQ

                                   (                                                                 )
                + ∫ δ d N ( x, tn−1 ) ⋅ ρ N ( x, tn−1 ) + B ( x, tn −1 ) ⋅ C ⋅ B ( x, tn−1 ) d n −1d Ω
                      Ω
                                +                 −               +                    −




3.1.   Discretization of BDG

       The expansion of the left hand or BDG side of Eq. [3-16] into matrix notation is


                          ∫ ∫              δ d T ( N T ρ N + BT [C ] B ) dd Ωdt
                              In       Ω
                                                                                                             [3-17]
                                              (                                                  )
                          + ∫ δ d T N T ( tn −1 ) ρ N ( tn −1 ) + BT ( tn −1 ) [C ] B ( tn −1 ) dd Ω
                                   Ω
                                           +             +              +                +




where the independent variables that B and N depend on ( x and t ) have been omitted except

for when the time is specified. More specifically, the expansion gives




                                                                    16 / 79
           ⎛ ⎡ N t N sT ⎤                                                      ⎞
           ⎜⎢ 1         ⎥                                                      ⎟
           ⎜⎢           ⎥                                                      ⎟
           ⎜ ⎢ N NT ⎥ ρ ⎡N N                     N ti N s         N tk N s ⎤ ⎟ ⎟
           ⎜ ⎢ ti s ⎥ ⎣ t1 s                                               ⎦
           ⎜⎢           ⎥                                                      ⎟
           ⎜⎢         T ⎥
                                                                               ⎟
           ⎜ ⎢ N tk N s ⎥
             ⎣          ⎦                                                      ⎟
∫In ∫Ω ⎜ ⎡ N BT ⎤
      δ dT ⎜                                                                   ⎟ dd Ωdt
                                                                               ⎟
                    t1 s
           ⎜ ⎢            ⎥                                                    ⎟
           ⎜ ⎢            ⎥                                                    ⎟
           ⎜ + ⎢ N t BsT ⎥ [C ] ⎡ N t Bs             N ti Bs         N tk Bs ⎤ ⎟
           ⎜ ⎢ i ⎥              ⎣ 1                                          ⎦⎟
           ⎜ ⎢            ⎥                                                    ⎟
           ⎜ ⎢ N BT ⎥                                                          ⎟
           ⎝ ⎢ tk s ⎥
               ⎣          ⎦                                                    ⎠
         ⎛ ⎡ N t ( tn −1 ) N sT ⎤
                     +                                                                                        ⎞
         ⎜⎢ 1                   ⎥                                                                             ⎟
         ⎜⎢                     ⎥                                                                             ⎟
         ⎜⎢                   T ⎥
                                                                                                              ⎟
         ⎜ ⎢ N ti ( tn −1 ) N s ⎥ ρ ⎡ N t1 ( tn −1 ) N s       N ti ( tn −1 ) N s      N tk ( tn −1 ) N s ⎤ ⎟
                     +                        +                        +                       +
                                    ⎣                                                                     ⎦
         ⎜⎢                     ⎥                                                                             ⎟
         ⎜⎢                     ⎥                                                                             ⎟
           ⎣ tk ( n −1 ) s ⎦
         ⎜ ⎢N t+ NT ⎥                                                                                         ⎟
+∫ δ d T ⎜
         ⎜ ⎡
                                                                                                              ⎟ dd Ω
                                                                                                              ⎟
         ⎜ ⎢ N t1 ( tn −1 ) Bs ⎤
  Ω                      +      T

                                  ⎥                                                                           ⎟
         ⎜ ⎢                      ⎥                                                                           ⎟      .[3-18]
         ⎜ ⎢                                                                                                  ⎟
         ⎜ + ⎢ N t ( tn −1 ) BsT ⎥ [C ] ⎡ N t ( tn −1 ) Bs
                         +
                                  ⎥
                                                    +
                                                                   N ti ( tn −1 ) Bs
                                                                           +
                                                                                          N tk ( tn −1 ) Bs ⎤ ⎟
                                                                                                   +

         ⎜ ⎢       i                    ⎣ 1                                                                 ⎦⎟
         ⎜ ⎢                      ⎥                                                                           ⎟
         ⎜ ⎢                      ⎥                                                                           ⎟
         ⎜ N tk ( tn −1 ) BsT ⎥
                         +
                                                                                                              ⎟
         ⎝ ⎣                      ⎦                                                                           ⎠




                                                       17 / 79
                   ⎛ ⎡ N t N sT ρ N t N s          N t1 N sT ρ N tk N s ⎤       ⎞
                   ⎜⎢ 1              1
                                                                        ⎥       ⎟
                   ⎜⎢                                                   ⎥       ⎟
                   ⎜ ⎢N NT ρ N N                   N tk N s ρ N t k N s ⎥
                                                           T                    ⎟
                   ⎜ ⎣ tk s         t1  s
                                                                        ⎦       ⎟
       ∫In ∫Ω δ dT ⎜                                                            ⎟ dd Ωdt
                   ⎜ ⎡ N t1 Bs [C ] N t1 Bs             N t1 BsT [C ] N tk Bs ⎤ ⎟
                                T

                   ⎜+⎢                                                        ⎥⎟
                   ⎜ ⎢                                                        ⎥⎟
                   ⎜ ⎢ N tk BsT [C ] N t1 Bs            N tk BsT [C ] N tk Bs ⎥ ⎟
                   ⎝ ⎣                                                        ⎦⎠
                ⎛ ⎡ N t ( tn −1 ) N sT ρ N t ( tn −1 ) N s
                            +                    +
                                                                  N t1 ( tn −1 ) N sT ρ N tk ( tn −1 ) N s ⎤
                                                                          +                     +                  ⎞
                ⎜⎢     1                    1
                                                                                                           ⎥       ⎟
                ⎜⎢                                                                                         ⎥       ⎟
                ⎜⎢                                                                                         ⎥       ⎟
                ⎜ ⎢ N tk ( tn −1 ) N s ρ N t1 ( tn −1 ) N s       N tk ( tn −1 ) N s ρ N tk ( tn −1 ) N s ⎥
                            +                    +                        +                     +
                                                                                                                   ⎟
                                     T                                              T
                  ⎣                                                                                        ⎦
       +∫ δ d T ⎜                                                                                                  ⎟ dd Ω
         Ω
                ⎜ ⎡ N t1 ( tn −1 ) BsT [C ] N t1 ( tn −1 ) Bs
                               +                        +
                                                                       N t1 ( tn −1 ) BsT [C ] N tk ( tn −1 ) Bs ⎤ ⎟
                                                                               +                       +

                ⎜   ⎢                                                                                            ⎥⎟
                ⎜+⎢                                                                                              ⎥⎟
                ⎜ ⎢                                                                                              ⎥⎟
                ⎜ ⎢ N tk ( tn −1 ) BsT [C ] N t1 ( tn −1 ) Bs
                                +                       +
                                                                       N tk ( tn −1 ) BsT [C ] N tk ( tn −1 ) Bs ⎦ ⎟
                                                                               +                       +
                                                                                                                 ⎥⎠
                ⎝ ⎣                                                                                                       [3-19]

The test function is arbitrary and can be moved outside of the integral, and the nodal

displacement vector d is constant and can also be moved outside of the integral. The spatial

integral can also be moved inside the matrices to act on the individual elements. With some

rearrangement the BDG of Eq. [3-16] now becomes




                                                          18 / 79
             ⎛ ⎡ N t N t N sT ρ N s d Ω                                                         ⎞
                 ∫Ω 1 1                                  ∫Ω Nt1 Ntk N s ρ N s d Ω ⎤ ⎟
                                                                        T
             ⎜⎢                                                                       ⎥
             ⎜⎢                                                                       ⎥         ⎟
             ⎜⎢                                                                       ⎥         ⎟
             ⎜ ⎢ ∫ N tk N t1 N sT ρ N s d Ω             ∫Ω  N tk N tk N sT ρ N s d Ω ⎥          ⎟
    δd ∫
       T     ⎜⎣ Ω                                                                     ⎦         ⎟ dtd
          In ⎜
                  ⎡ N t N t BsT [C ] Bs d Ω                                                   ⎤⎟
             ⎜ ⎢ ∫Ω 1 1                                      ∫Ω Nt1 Ntk Bs [C ] Bs d Ω ⎥ ⎟
                                                                            T


             ⎜+⎢                                                                              ⎥⎟
             ⎜ ⎢                                                                              ⎥⎟
             ⎜ ⎢ ∫ N tk N t1 BsT CBs d Ω                     ∫Ω Ntk Ntk Bs [C ] Bs d Ω ⎥ ⎠      ⎟
                                                                            T
             ⎝    ⎣ Ω                                                                         ⎦
           ⎛ ⎡ N t ( tn −1 ) N t ( tn −1 ) N sT ρ N s d Ω
               ∫Ω 1                                                      ∫Ω Nt1 ( tn−1 ) Ntk ( tn−1 ) N s ρ N s d Ω ⎤ ⎞
                         +              +                                          +              +         T
           ⎜⎢                      1
                                                                                                                        ⎥        ⎟
           ⎜⎢                                                                                                           ⎥        ⎟
           ⎜⎢                                                                                                           ⎥        ⎟
           ⎜ ⎢ ∫ N tk ( tn −1 ) N t1 ( tn −1 ) N sT ρ N s d Ω
                         +              +
                                                                         ∫Ω Ntk ( tn−1 ) Ntk ( tn−1 ) N s ρ N s d Ω ⎥ ⎟
                                                                                   +              +         T
             ⎣ Ω                                                                                                        ⎦
    +δ d T ⎜                                                                                                                     ⎟d
           ⎜ ⎡ N ( t + ) N ( t + ) BT [ C ] B d Ω
           ⎜ ⎢ ∫Ω t1 n −1 t1 n −1 s                                           ∫Ω Nt1 ( tn−1 ) Ntk ( tn−1 ) Bs [C ] Bs d Ω ⎤ ⎟
                                                                                         +              +       T

                                                                                                                               ⎥⎟
                                                             s

           ⎜+⎢                                                                                                                 ⎥⎟
           ⎜ ⎢                                                                                                                 ⎥⎟
           ⎜ ⎢ ∫ N tk ( tn −1 ) N t1 ( tn −1 ) BsT [C ] Bs d Ω
                             +              +
                                                                             ∫Ω  N tk ( tn −1 ) N tk ( tn −1 ) BsT [C ] Bs d Ω ⎥ ⎟ .[3-20]
                                                                                         +              +
           ⎝ ⎣ Ω                                                                                                               ⎦⎠

Taking the temporal shape functions out of the spatial integrals because they are constant in

space leads to

             ⎛ ⎡ N t N t N sT ρ N s d Ω                N t1 N tk ∫ N sT ρ N s d Ω ⎤       ⎞
             ⎜ ⎢ 1 1 ∫Ω                                           Ω               ⎥       ⎟
             ⎜⎢                                                                   ⎥       ⎟
             ⎜⎢                                                                   ⎥       ⎟
             ⎜ ⎢ N tk N t1 ∫ N sT ρ N s d Ω            N tk N tk ∫ N sT ρ N s d Ω ⎥       ⎟
               ⎣            Ω                                     Ω               ⎦
    δ dT ∫ ⎜                                                                              ⎟ dtd
          In ⎜
                 ⎡N N           BT [C ] Bs d Ω              N t1 N tk ∫ Bs [C ] Bs d Ω ⎤ ⎟
             ⎜ ⎢ t1 t1 ∫Ω s
                                                                         T
                                                                       Ω                ⎥⎟
             ⎜+⎢                                                                        ⎥⎟
             ⎜ ⎢                                                                        ⎥⎟
             ⎜ ⎢ N tk N t1 ∫ BsT [C ] Bs d Ω                N tk N tk ∫ BsT [C ] Bs d Ω ⎥ ⎟
             ⎝ ⎣              Ω                                        Ω                ⎦⎠
            ⎛ ⎡ N t ( tn −1 ) N t ( tn −1 ) N sT ρ N s d Ω            N t1 ( tn −1 ) N tk ( tn −1 ) ∫ N sT ρ N s d Ω ⎤       ⎞
                                              ∫Ω
                        +              +                                      +              +
            ⎜⎢ 1                  1                                                                  Ω               ⎥       ⎟
            ⎜⎢                                                                                                       ⎥       ⎟
            ⎜⎢                                                                                                       ⎥       ⎟
            ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) ∫ N sT ρ N s d Ω
                        +              +
                                                                      N tk ( tn −1 ) N tk ( tn −1 ) ∫ N sT ρ N s d Ω ⎥
                                                                              +              +
                                                                                                                             ⎟
              ⎣                                Ω                                                     Ω               ⎦
     +δ d T ⎜                                                                                                                ⎟d
            ⎜ ⎡ N ( t + ) N ( t + ) BT [ C ] B d Ω                         N t1 ( tn −1 ) N tk ( tn −1 ) ∫ BsT [C ] Bs d Ω ⎤ ⎟
            ⎜ ⎢ t1 n −1 t1 n −1 ∫Ω s
                                                                                   +              +

                                                                                                                           ⎥⎟
                                                              s
                                                                                                          Ω
            ⎜+⎢                                                                                                            ⎥⎟
            ⎜ ⎢                                                                                                            ⎥⎟
            ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) ∫ BsT [C ] Bs d Ω
                            +              +
                                                                           N tk ( tn −1 ) N tk ( tn −1 ) ∫ BsT [C ] Bs d Ω ⎥ ⎟ .[3-21]
                                                                                   +              +
            ⎝ ⎣                                  Ω                                                        Ω                ⎦⎠


                                                             19 / 79
The spatial consistent mass and stiffness matrices are defined as

                                                                  M = ∫ N sT ρ N s d Ω                                                        [3-22]
                                                                             Ω



                                                                 K = ∫ BsT [C ] Bs d Ω .                                                      [3-23]
                                                                         Ω


Both of these are assumed to be constant in time because the volumetric mass and elastic moduli

tensor are both assumed to be constant in time. Inserting these into Eq. [3-21] gives

                          ⎛ ⎡ Nt Nt M                                  Nt1 Ntk M ⎤ ⎡ N t1 N t1 K                   Nt1 N tk K ⎤ ⎞
                          ⎜⎢ 1 1                                                 ⎥ ⎢                                           ⎥⎟
            BDG   =δd ∫ ⎜⎢
                     T
                                                                                 ⎥+⎢                                           ⎥ ⎟ dtd
                       In
                          ⎜ ⎢N N M
                          ⎜ tk t1                                      Ntk Ntk M ⎥ ⎢ N tk N t1 K                   N tk N tk K ⎥ ⎟
                                                                                                                                 ⎟
                          ⎝⎣                                                     ⎦ ⎣                                           ⎦⎠
                        ⎛ ⎡ N t ( tn −1 ) Nt ( tn −1 ) M
                                    +              +
                                                                                      Nt1 ( tn −1 ) Ntk ( tn −1 ) M ⎤ ⎞
                                                                                              +              +

                        ⎜⎢ 1                  1
                                                                                                                      ⎥⎟
                        ⎜⎢                                                                                            ⎥⎟                  .   [3-24]
                        ⎜⎢                                                                                            ⎥⎟
                        ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M                           N tk ( tn −1 ) N tk ( tn −1 ) M ⎥ ⎟
                                    +              +                                          +              +
                       T ⎣                                                                                            ⎦
                   +δ d ⎜                                                                                               ⎟d
                        ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K
                                       +              +
                                                                                       N t1 ( tn −1 ) Ntk ( tn −1 ) K ⎤ ⎟
                                                                                                 +              +

                        ⎜   ⎢                                                                                          ⎥⎟
                        ⎜+⎢                                                                                            ⎥⎟
                        ⎜ ⎢                                                                                            ⎥⎟
                        ⎜ ⎢ N tk ( tn −1 ) Nt1 ( tn −1 ) K
                                        +              +
                                                                                       Ntk ( tn −1 ) Ntk ( tn −1 ) K ⎥ ⎟
                                                                                                 +              +
                        ⎝ ⎣                                                                                            ⎦⎠

3.2.    Discretization of LDG

        The expansion of the right hand or LDG side of Eq. [3-16] is


                  ∫(            δ d ( N ( x, t ) ⋅ h ) d ϒ + ∫ δ d ( N ( x, t ) ⋅ f ) dQ
                       ϒ h )n                                          Qn



                                     (                                                                                )
                                                                                                                                      .       [3-25]
                  + ∫ δ d N ( x, tn −1 ) ⋅ ρ N ( x, tn −1 ) + B ( x, tn −1 ) ⋅ C ⋅ B ( x, tn −1 ) d n−1 d Ω
                                  +                  −                +                    −
                        Ω



Using the notations in Eq. [2-24] and [2-22] while also taking the test function and the nodal

displacement vectors outside of the integrations yields

                  δd∫
                            In   ∫ ( N ( x , t ) ⋅ h ) d Γdt + δ d ∫ ∫ ( N ( x , t ) ⋅ f ) d Ωdt
                                     Γh                                          In   Ω
                                                                                                                                      .       [3-26]
                  +δ d ∫
                                 Ω   (                                                                                )
                                          N ( x , t n −1 ) ⋅ ρ N ( x , t n −1 ) + B ( x , t n −1 ) ⋅ C ⋅ B ( x , t n −1 ) d Ωd n −1
                                                    +                    −                  +                      −




In matrix notation, this is

                                                                             20 / 79
                                  δ dT ∫                  ∫ ( Nh ) d Γdt + δ d ∫ ∫ ( Nf ) d Ωdt
                                                                                           T
                                                     In       Γh                                   In   Ω

                                  +δ d               ∫ ( N ( t ) ρ N ( t ) + B ( t ) [C ] B ( t ) ) d Ωd                                                                 [3-27]
                                                 T                 T    +               −                   T    +               −
                                                          Ω             n −1            n −1                     n −1            n −1          n −1



where the independent variables that B and N depend on ( x and t ) have been omitted expect

for when the time is specified. Expanding Eq. [3-27] into full matrix notation yields

      δ dT ∫
               In
                  (⎡N N
                   ⎣∫    Γh
                                  t1         s                     N ti N s             N tk N s ⎤ [ h ] d Γdt
                                                                                                 ⎦               )
      +δ d T   ∫ ∫ (⎡N N
                     ⎣
                    In        Ω        t1        s                  N ti N s             N tk N s ⎤ [ f ] d Ωdt
                                                                                                  ⎦                  )
                ⎡ N t ( tn −1 ) N sT ⎤
                         +

                ⎢    1
                                     ⎥
                ⎢                    ⎥
                ⎢                    ⎥
      +δ d T ∫ ⎢ N ti ( tn −1 ) N sT ⎥ ρ ⎡ N t1 ( tn −1 ) N s
                         +
                                         ⎣
                                                   −
                                                                                                          N ti ( tn −1 ) N s
                                                                                                                  −
                                                                                                                                        N tk ( tn −1 ) N s ⎤ d Ωd n −1
                                                                                                                                                −
                                                                                                                                                           ⎦
              Ω
                ⎢                    ⎥
                ⎢                    ⎥
                ⎢ N (t + ) N T ⎥
                ⎣ tk n −1 s ⎦
                ⎡ N t ( tn −1 ) BsT ⎤
                         +

                ⎢ 1                 ⎥
                ⎢                   ⎥
                ⎢                   ⎥
      +δ d T ∫ ⎢ N ti ( tn −1 ) BsT ⎥ [C ] ⎣ N t1 ( tn −1 ) Bs
                         +
                                           ⎡         −
                                                                                                            N ti ( tn −1 ) Bs
                                                                                                                    −
                                                                                                                                        N tk ( tn −1 ) Bs ⎦ d Ωd n −1
                                                                                                                                                −
                                                                                                                                                          ⎤
              Ω
                ⎢                   ⎥
                ⎢                   ⎥
                ⎢ N ( t + ) BT ⎥
                ⎣ tk n −1 s ⎦                                                                                                                                            .[3-28]

      δ dT ∫
               In∫ ( ⎡ N N [h]
                     ⎣   Γh        t1        s                                      ) d Γdt
                                                                        N t k N s [ h ]⎤
                                                                                       ⎦
                                                                                           T




                ∫ ∫ (⎡N N [ f ]                                           N N [ f ]⎤ ) d Ωdt
                                                                                               T
      +δ d T           ⎣                t1       s                             tk  ⎦s
                    In        Ω

              ⎛ ⎡ N t ( tn −1 ) N sT ρ N t ( tn −1 ) N s
                          +                    −
                                                                                                        N t1 ( tn −1 ) N sT ρ N tk ( tn −1 ) N s ⎤
                                                                                                                +                     −                     ⎞
              ⎜⎢ 1                        1
                                                                                                                                                 ⎥          ⎟
              ⎜⎢                                                                                                                                 ⎥          ⎟
              ⎜⎢                                                                                                                                 ⎥          ⎟
              ⎜ ⎢ N tk ( tn −1 ) N s ρ N t1 ( tn −1 ) N s                                               N tk ( tn −1 ) N sT ρ N tk ( tn −1 ) N s ⎥
                          +                    −                                                                +                     −
                                                                                                                                                            ⎟
                                   T
                ⎣                                                                                                                                ⎦
      +δ d ∫ ⎜
          T
                                                                                                                                                            ⎟ d Ωd n −1
            Ω
              ⎜ ⎡ N t1 ( tn −1 ) BsT [C ] N t1 ( tn −1 ) Bs
                             +                        −
                                                                                                                N t1 ( tn −1 ) BsT [C ] N tk ( tn −1 ) Bs ⎤ ⎟
                                                                                                                        +                       −

              ⎜ ⎢                                                                                                                                         ⎥⎟
              ⎜+⎢                                                                                                                                         ⎥⎟
              ⎜ ⎢                                                                                                                                         ⎥⎟            [3-29]
              ⎜ ⎢ N tk ( tn −1 ) BsT [C ] N t1 ( tn −1 ) Bs
                              +                       −
                                                                                                                N tk ( tn −1 ) BsT [C ] N tk ( tn −1 ) Bs ⎥ ⎟
                                                                                                                        +                       −
              ⎝ ⎣                                                                                                                                         ⎦⎠




                                                                                           21 / 79
The spatial integral can be moved inside the matrices to act on the individual elements. The

temporal shape functions can be taken out of the spatial integrals because they are constant in

space.

    δ dT ∫
          In   ∫Γh      (   ⎡ Nt1 N s [ h]
                            ⎣                                              ) d Γdt
                                                            N t k N s [ h ]⎤
                                                                           ⎦
                                                                                     T




    +δ d T ∫        ∫       ( ⎡⎣ N N [ f ]                       N N [ f ]⎤ ) d Ωdt
                                                                                         T

               In    Ω           t1   s                            tk     ⎦  s


           ⎛ ⎡ N t ( tn −1 ) Nt ( tn −1 ) N sT ρ N s d Ω                                           Nt1 ( tn −1 ) Ntk ( tn −1 ) ∫ N sT ρ N s d Ω ⎤       ⎞
                                           ∫Ω
                      +             −                                                                     +             −
           ⎜⎢     1            1                                                                                                Ω               ⎥       ⎟
           ⎜⎢                                                                                                                                   ⎥       ⎟
           ⎜⎢                                                                                                                                   ⎥       ⎟
           ⎜ ⎢ Ntk ( tn −1 ) Nt1 ( tn −1 ) ∫ N sT ρ N s d Ω
                      +             −
                                                                                                   Ntk ( tn−1 ) Ntk ( tn −1 ) ∫ N sT ρ N s d Ω ⎥
                                                                                                          +             −
                                                                                                                                                        ⎟
         T ⎜⎣                               Ω                                                                                   Ω               ⎦       ⎟ d n −1
    +δ d
           ⎜ ⎡ N ( t + ) N ( t − ) BT [ C ] B d Ω                                                        Nt1 ( tn −1 ) Ntk ( tn −1 ) ∫ Bs [C ] Bs d Ω ⎤⎟
           ⎜ ⎢ t1 n −1 t1 n −1 ∫Ω s
                                                                                                                +             −         T

                                                                                                                                                      ⎥⎟
                                                           s                                                                          Ω
           ⎜+⎢                                                                                                                                        ⎥⎟
           ⎜ ⎢                                                                                                                                        ⎥⎟         [3-30]
           ⎜ ⎢ Ntk ( tn−1 ) Nt1 ( tn −1 ) ∫ BsT [C ] Bs d Ω
                          +             −
                                                                                                         Ntk ( tn −1 ) Ntk ( tn −1 ) ∫ Bs [C ] Bs d Ω ⎥ ⎟
                                                                                                                +             −         T
           ⎝   ⎣                              Ω                                                                                       Ω               ⎦⎠

Using the substitutions in equations Eq. [3-22] and [3-23] for the spatial mass and stiffness

matrices gives

                              LDG = δ d T ∫
                                                   In   ∫
                                                        Γh       ( ⎡⎣ N N [h]
                                                                        t1       s                 N t k N s [ h ]⎤
                                                                                                                  ⎦
                                                                                                                      T
                                                                                                                          ) d Γdt
                              +δ d T ∫
                                          In   ∫ (⎡N N [ f ]
                                               Ω  ⎣         t1    s                          N tk N s [ f ]⎤
                                                                                                           ⎦
                                                                                                               T
                                                                                                                   ) d Ωdt
                                   ⎛ ⎡ N t ( tn −1 ) N t ( tn −1 ) M
                                               +              −
                                                                                                       N t1 ( tn −1 ) N tk ( tn −1 ) M ⎤
                                                                                                               +              −           ⎞
                                   ⎜⎢     1              1
                                                                                                                                       ⎥  ⎟
                                   ⎜⎢                                                                                                  ⎥  ⎟
                                   ⎜⎢                                                                                                  ⎥  ⎟
                                   ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M                                 N tk ( tn −1 ) N tk ( tn −1 ) M ⎥
                                               +              −                                                +              −
                                  T ⎣                                                                                                  ⎦  ⎟
                              +δ d ⎜                                                                                                      ⎟ d n −1
                                   ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K
                                                  +              −
                                                                                                        N t1 ( tn −1 ) N tk ( tn −1 ) K ⎤ ⎟
                                                                                                                +              −

                                   ⎜ ⎢                                                                                                  ⎥⎟
                                   ⎜+⎢                                                                                                  ⎥⎟
                                   ⎜ ⎢                                                                                                  ⎥⎟         .       [3-31]
                                   ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) K
                                                   +              −
                                                                                                        N tk ( tn −1 ) N tk ( tn −1 ) K ⎥ ⎟
                                                                                                                +              −
                                   ⎝ ⎣                                                                                                  ⎦⎠

The surface traction and body forces can often be assumed to be functions that can be separated

into the multiplication of distinct spatial and temporal functions.

                                                                             h ( x , t ) = hs ( x ) ht ( t )                                               [3-32]


                                                                                             22 / 79
                                                          f ( x, t ) = f s ( x ) ft ( t )                                       [3-33]

Inserting these into Eq. [3-31] gives

                 LDG = δ d T ∫
                                       In   ∫
                                            Γh   ( ⎡⎣ N N [h h ]
                                                     t1    s    s t             N tk N s [ hs ht ]⎤
                                                                                                  ⎦
                                                                                                        T
                                                                                                            ) d Γdt
                 +δ d T ∫
                          In   ∫
                               Ω   (   ⎡ N t1 N s [ f s f t ]
                                       ⎣                                N t k N s [ f s f t ]⎤
                                                                                             ⎦
                                                                                                 T
                                                                                                     ) d Ωdt
                      ⎛ ⎡ N t ( tn −1 ) N t ( tn −1 ) M
                                  +              −
                                                                              N t1 ( tn −1 ) N tk ( tn −1 ) M ⎤
                                                                                      +              −             ⎞
                      ⎜⎢     1              1
                                                                                                              ⎥    ⎟
                      ⎜⎢                                                                                      ⎥    ⎟
                      ⎜⎢                                                                                      ⎥    ⎟
                      ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M                     N tk ( tn −1 ) N tk ( tn −1 ) M ⎥
                                  +              −                                    +              −
                     T ⎣                                                                                      ⎦    ⎟
                 +δ d ⎜                                                                                            ⎟ d n −1
                      ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K
                                     +              −
                                                                                 N t1 ( tn −1 ) N tk ( tn −1 ) K ⎤ ⎟
                                                                                         +              −

                      ⎜ ⎢                                                                                        ⎥⎟
                      ⎜+⎢                                                                                        ⎥⎟
                      ⎜ ⎢                                                                                        ⎥⎟         .   [3-34]
                      ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) K
                                      +              −
                                                                                 N tk ( tn −1 ) N tk ( tn −1 ) K ⎥ ⎟
                                                                                         +              −
                      ⎝ ⎣                                                                                        ⎦⎠

The functions ft and ht can be taken out of the body force and traction vectors. The temporal

shape functions and temporal components of the traction and body forces can be taken out of the

spatial integrals because they are constant in space.

                                 ⎛                                                          T
                                                                                              ⎞
                LDG = δ d T ∫ ⎜ ⎡ Nt1 ht1 ∫ N s [ hs ] d Γ       Ntk htk ∫ N s [ hs ] d Γ ⎤ ⎟ dt
                               In ⎢                                                       ⎥ ⎠
                                 ⎝⎣         Γh                            Γh              ⎦

                +δ d T ∫ ⎛ ⎡ N t1 ft1 ∫ N s [ f s ] d Ω    Ntk ftk ∫ N s [ f s ] d Ω ⎤ ⎞ dt
                                                                                      T
                           ⎜
                        In ⎝ ⎣                                                       ⎦ ⎟
                                       Ω                            Ω                   ⎠
                     ⎛ ⎡ N t ( tn −1 ) N t ( tn −1 ) M
                                 +              −
                                                                             Nt1 ( tn −1 ) Ntk ( tn −1 ) M ⎤
                                                                                     +             −            ⎞
                     ⎜⎢     1              1
                                                                                                            ⎥   ⎟
                     ⎜⎢                                                                                     ⎥   ⎟
                     ⎜⎢                                                                                     ⎥   ⎟
                     ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M                     N tk ( tn −1 ) Ntk ( tn −1 ) M ⎥
                                 +              −                                    +             −
                    T ⎣                                                                                     ⎦   ⎟
                +δ d ⎜                                                                                          ⎟ d n −1
                     ⎜ ⎡ N t1 ( tn −1 ) Nt1 ( tn −1 ) K
                                    +              −
                                                                               Nt1 ( tn −1 ) Ntk ( tn −1 ) K ⎤ ⎟
                                                                                       +             −

                     ⎜ ⎢                                                                                      ⎥⎟
                     ⎜+⎢                                                                                      ⎥⎟
                     ⎜ ⎢                                                                                      ⎥⎟                [3-35]
                     ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) K
                                     +              −
                                                                               N tk ( tn −1 ) Ntk ( tn −1 ) K ⎥ ⎟
                                                                                       +             −
                     ⎝ ⎣                                                                                      ⎦⎠

Defining the spatial body force and surface traction amplitude vectors (with units of force) as

                                                           Φ = ∫ N s [ hs ] d Γ                                                 [3-36]
                                                                   Γh




                                                                    23 / 79
                                                      F = ∫ Ns [ fs ] d Ω                                                    [3-37]
                                                             Ω


and substituting into Eq. [3-35] yields

                 LDG = δ d T ∫
                                  In    ( ⎡⎣ N h Φ
                                              t1 t1           N tk htk Φ ⎤
                                                                         ⎦
                                                                                T
                                                                                    ) dt
                 +δ d T ∫
                          In   ( ⎡⎣ N f F
                                   t1    t1            N t k f tk F ⎤
                                                                    ⎦
                                                                        T
                                                                            ) dt
                      ⎛ ⎡ N t ( tn −1 ) N t ( tn −1 ) M
                                  +              −
                                                                             N t1 ( tn −1 ) N tk ( tn −1 ) M ⎤
                                                                                     +              −           ⎞
                      ⎜⎢     1              1
                                                                                                             ⎥  ⎟
                      ⎜⎢                                                                                     ⎥  ⎟
                      ⎜⎢                                                                                     ⎥  ⎟
                      ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M                    N tk ( tn −1 ) N tk ( tn −1 ) M ⎥
                                  +              −                                   +              −
                     T ⎣                                                                                     ⎦  ⎟
                 +δ d ⎜                                                                                         ⎟ d n −1
                      ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K
                                     +              −
                                                                              N t1 ( tn −1 ) N tk ( tn −1 ) K ⎤ ⎟
                                                                                      +              −

                      ⎜ ⎢                                                                                     ⎥⎟
                      ⎜+⎢                                                                                     ⎥⎟
                      ⎜ ⎢                                                                                     ⎥⎟         .   [3-38]
                      ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) K
                                      +              −
                                                                              N tk ( tn −1 ) N tk ( tn −1 ) K ⎥ ⎟
                                                                                      +              −
                      ⎝ ⎣                                                                                     ⎦⎠


                                                 (
                 LDG = δ d T ∫ ⎛ ⎡ Nt1 ht1 Φ + f t1 F            )                     (                  )
                                                                                    N tk htk Φ + f tk F ⎤ ⎞ dt
                                                                                                         T
                                 ⎜
                              In ⎝ ⎣                                                                    ⎦ ⎟⎠
                      ⎛ ⎡ N t ( tn −1 ) N t ( tn −1 ) M
                                  +              −
                                                                             N t1 ( tn −1 ) N tk ( tn −1 ) M ⎤
                                                                                     +              −            ⎞
                      ⎜⎢     1              1
                                                                                                             ⎥   ⎟
                      ⎜⎢                                                                                     ⎥   ⎟
                      ⎜⎢                                                                                     ⎥   ⎟
                      ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M                    N tk ( tn −1 ) N tk ( tn −1 ) M ⎥               [3-39]
                                  +              −                                   +              −
                     T ⎣                                                                                     ⎦   ⎟
                 +δ d ⎜                                                                                          ⎟ d n −1
                      ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K
                                     +              −
                                                                               N t1 ( tn −1 ) N tk ( tn −1 ) K ⎤ ⎟
                                                                                       +              −

                      ⎜ ⎢                                                                                      ⎥⎟
                      ⎜+⎢                                                                                      ⎥⎟
                      ⎜ ⎢                                                                                      ⎥⎟
                      ⎜ ⎢ Ntk ( tn −1 ) N t1 ( tn −1 ) K
                                      +              −
                                                                               N tk ( tn −1 ) N tk ( tn −1 ) K ⎥ ⎟
                                                                                       +              −
                      ⎝ ⎣                                                                                      ⎦⎠

3.3.   Combination of Discretized BDG and LDG

Combining the updated BDG and LDG equations results in




                                                             24 / 79
                  ⎛ ⎛ ⎡ Nt Nt M                        N t1 N tk M ⎤ ⎡ N t1 N t1 K                  N t1 N tk K ⎤ ⎞ ⎞
                  ⎜ ⎜⎢ 1 1                                         ⎥ ⎢                                            ⎥⎟ ⎟
                  ⎜∫ ⎜⎢                                            ⎥+⎢                                            ⎥ ⎟ dt ⎟
                  ⎜ In ⎜ ⎢                                                                                               ⎟
                  ⎜ ⎜ ⎣ N tk N t1 M                   N tk N tk M ⎥ ⎢ N tk N t1 K                   N tk N tk K ⎥ ⎟ ⎟
                                                                                                                    ⎟
                        ⎝                                          ⎦ ⎣                                            ⎦⎠
                  ⎜                                                                                                      ⎟
                  ⎜ ⎛ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) M
                                   +              +
                                                                      N t1 ( tn −1 ) N tk ( tn −1 ) M ⎤ ⎞
                                                                               +              +
                                                                                                                         ⎟
                  ⎜ ⎜⎢                                                                                  ⎥ ⎟              ⎟
            δ dT ⎜ ⎜ ⎢                                                                                  ⎥ ⎟              ⎟d
                  ⎜    ⎜⎢                                                                               ⎥ ⎟              ⎟
                  ⎜ ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M               N tk ( tn −1 ) N tk ( tn −1 ) M ⎥ ⎟
                                   +              +                             +              +
                         ⎣                                                                              ⎦                ⎟
                  ⎜+⎜                                                                                      ⎟             ⎟
                  ⎜ ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K                 N t1 ( tn −1 ) N tk ( tn −1 ) K ⎤ ⎟
                                      +              +                            +              +
                                                                                                                         ⎟
                  ⎜ ⎜ ⎢                                                                                  ⎥⎟              ⎟
                  ⎜ ⎜+⎢                                                                                  ⎥⎟              ⎟
                  ⎜ ⎜ ⎢N t+ N t+ K
                  ⎜ ⎜ ⎢ tk ( n −1 ) t1 ( n −1 )
                                                                                                         ⎥⎟
                                                                        N tk ( tn −1 ) N tk ( tn −1 ) K ⎥ ⎟
                                                                                  +              +                       ⎟
                                                                                                                         ⎟
                  ⎝ ⎝ ⎣                                                                                  ⎦⎠              ⎠

                                  (
            = δ d T ∫ ⎛ ⎡ N t1 ht1 Φ + f t1 F       )                 (                  )
                                                                 N tk htk Φ + ftk F ⎤ ⎞ dt
                                                                                              T
                        ⎜
                     In ⎝ ⎣                                                                 ⎦ ⎟  ⎠
                    ⎛ ⎡ N t ( tn −1 ) N t ( tn −1 ) M
                                +              −
                                                                    N t1 ( tn −1 ) N tk ( tn −1 ) M ⎤ ⎞
                                                                            +              −

                    ⎜⎢ 1                  1
                                                                                                      ⎥ ⎟
                    ⎜⎢                                                                                ⎥ ⎟
                    ⎜⎢                                                                                ⎥ ⎟
                    ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M             N tk ( tn −1 ) N tk ( tn −1 ) M ⎥ ⎟
                                +              −                            +              −
                      ⎣                                                                               ⎦
            +δ d T ⎜                                                                                     ⎟ d n −1
                    ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K
                                   +              −
                                                                     N t1 ( tn −1 ) N tk ( tn −1 ) K ⎤ ⎟
                                                                               +              −

                    ⎜ ⎢                                                                                ⎥⎟
                    ⎜+⎢                                                                                ⎥⎟
                    ⎜    ⎢                                                                             ⎥⎟                   .   [3-40]
                    ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) K
                                    +              −
                                                                     N tk ( tn −1 ) N tk ( tn −1 ) K ⎥ ⎟
                                                                               +              −
                    ⎝ ⎣                                                                                ⎦⎠

Dropping the arbitrary test functions from both sides provides




                                                               25 / 79
               ⎛ ⎛ ⎡ Nt Nt M                        N t1 N tk M ⎤ ⎡ N t1 N t1 K                  N t1 N tk K ⎤ ⎞ ⎞
               ⎜ ⎜⎢ 1 1                                         ⎥ ⎢                                           ⎥⎟ ⎟
               ⎜∫ ⎜⎢                                            ⎥+⎢                                           ⎥ ⎟ dt ⎟
               ⎜ In ⎜ ⎢                                                                                              ⎟
               ⎜ ⎜ ⎣ N tk N t1 M                   N tk N tk M ⎥ ⎢ Ntk N t1 K                    N tk N t k K ⎥ ⎟ ⎟
                                                                                                                ⎟
                    ⎝                                           ⎦ ⎣                                           ⎦⎠
               ⎜                                                                                                     ⎟
               ⎜ ⎛ ⎡ N t1 ( tn −1 ) Nt1 ( tn −1 ) M
                                +              +
                                                                   Nt1 ( tn −1 ) N tk ( tn −1 ) M ⎤ ⎞
                                                                            +              +
                                                                                                                     ⎟
               ⎜ ⎜⎢                                                                                  ⎥ ⎟             ⎟
               ⎜ ⎜⎢                                                                                  ⎥ ⎟             ⎟d
               ⎜   ⎜⎢                                                                                ⎥ ⎟             ⎟
               ⎜ ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M               N tk ( tn −1 ) N tk ( tn −1 ) M ⎥ ⎟
                                +              +                            +              +
                     ⎣                                                                               ⎦               ⎟
               ⎜+⎜                                                                                      ⎟            ⎟
               ⎜ ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K                 N t1 ( tn −1 ) N tk ( tn −1 ) K ⎤ ⎟
                                   +              +                           +              +
                                                                                                                     ⎟
               ⎜ ⎜ ⎢                                                                                  ⎥⎟             ⎟
               ⎜ ⎜+⎢                                                                                  ⎥⎟             ⎟
               ⎜ ⎜ ⎢N t+ N t+ K
               ⎜ ⎜ ⎢ tk ( n −1 ) t1 ( n −1 )                         Ntk ( tn −1 ) Ntk ( tn −1 ) K ⎥ ⎟
                                                                              +              +        ⎥⎟             ⎟
                                                                                                                     ⎟
               ⎝ ⎝ ⎣                                                                                  ⎦⎠             ⎠

                               (
               = ∫ ⎛ ⎡ N t1 ht1 Φ + ft1 F        )                 (                  )
                                                               N tk htk Φ + ftk F ⎤ ⎞ dt
                                                                                            T
                      ⎜
                   In ⎝ ⎣                                                                 ⎦ ⎟  ⎠
                 ⎛ ⎡ N t ( tn −1 ) N t ( tn −1 ) M
                             +              −
                                                                 N t1 ( tn −1 ) N tk ( tn −1 ) M ⎤ ⎞
                                                                         +              −

                 ⎜⎢ 1                  1
                                                                                                   ⎥ ⎟
                 ⎜⎢                                                                                ⎥ ⎟
                 ⎜⎢                                                                                ⎥ ⎟
                 ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M             N tk ( tn −1 ) N tk ( tn −1 ) M ⎥ ⎟
                             +              −                            +              −
                   ⎣                                                                               ⎦
               +⎜                                                                                     ⎟ d n −1
                 ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K
                                 +              −
                                                                   N t1 ( tn −1 ) N tk ( tn −1 ) K ⎤ ⎟
                                                                            +              −

                 ⎜ ⎢                                                                                ⎥⎟
                 ⎜+⎢                                                                                ⎥⎟
                 ⎜    ⎢                                                                             ⎥⎟                  .   [3-41]
                 ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) K
                                 +              −
                                                                   N tk ( tn −1 ) N tk ( tn −1 ) K ⎥ ⎟
                                                                            +              −
                 ⎝ ⎣                                                                                ⎦⎠

Using substitutions from Eq. [3-2], [3-3], and [3-6]




                                                              26 / 79
          ⎛ ⎛ ⎡ Nt Nt M                      N t1 N tk M ⎤ ⎡ N t1 N t1 K                N t1 N tk K ⎤ ⎞ ⎞
          ⎜ ⎜⎢ 1 1                                       ⎥ ⎢                                         ⎥⎟ ⎟
          ⎜∫ ⎜⎢                                          ⎥+⎢                                         ⎥ ⎟ dt ⎟
          ⎜ In ⎜ ⎢                                                                                          ⎟
          ⎜ ⎜ ⎣ N tk N t1 M                 N tk N tk M ⎥ ⎢ N tk N t1 K                 N tk N t k K ⎥ ⎟ ⎟
                                                                                                       ⎟
               ⎝                                         ⎦ ⎣                                         ⎦⎠
          ⎜                                                                                                 ⎟
          ⎜ ⎛ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) M
                        +              +
                                                          N t1 ( tn −1 ) N tk ( tn −1 ) M ⎤ ⎞
                                                                  +              +
                                                                                                            ⎟
          ⎜ ⎜⎢                                                                              ⎥ ⎟             ⎟
          ⎜ ⎜⎢                                                                              ⎥ ⎟             ⎟d
          ⎜  ⎜⎢                                                                             ⎥ ⎟             ⎟
          ⎜ ⎜ ⎢ N tk ( tn −1 ) N t1 ( tn −1 ) M           N tk ( tn −1 ) N tk ( tn −1 ) M ⎥ ⎟
                        +              +                          +              +
               ⎣                                                                            ⎦               ⎟
          ⎜+⎜                                                                                 ⎟             ⎟
          ⎜ ⎜ ⎡ N t1 ( tn −1 ) N t1 ( tn −1 ) K            N t1 ( tn −1 ) N tk ( tn −1 ) K ⎤ ⎟
                            +              +                         +              +
                                                                                                            ⎟
          ⎜ ⎜ ⎢                                                                              ⎥⎟             ⎟
          ⎜ ⎜+⎢                                                                              ⎥⎟             ⎟
          ⎜ ⎜ ⎢N t+ N t+ K
          ⎜ ⎜ ⎢ tk ( n −1 ) t1 ( n −1 )
                                                                                             ⎥⎟
                                                           N tk ( tn −1 ) N tk ( tn −1 ) K ⎥ ⎟
                                                                     +              +                       ⎟
                                                                                                            ⎟
          ⎝ ⎝ ⎣                                                                              ⎦⎠             ⎠

                ⎜⎢ 1 1
                          (
                ⎛ ⎡ N t ht Φ + ft F ⎤ ⎞
                                  1
                                       ⎥⎟
                                           )    ⎛ ⎡ N t ( tn −1 ) M ⎤
                                                ⎜⎢ 1
                                                            +

                                                                     ⎥
                                                                                     ⎡ N t ( tn −1 ) K ⎤
                                                                                     ⎢    1
                                                                                               +

                                                                                                        ⎥
                                                                                                                      ⎞
                                                                                                                      ⎟
          = ∫ ⎜⎢                       ⎥ ⎟ dt + ⎜ ⎢                  ⎥ u ( tn −1 ) + ⎢                  ⎥ u ( tn −1 ) ⎟
                                                                            −                                  −
             In
                ⎜⎢                     ⎥⎟       ⎜⎢                   ⎥               ⎢                  ⎥             ⎟
                ⎝⎣       (
                ⎜ ⎢ N tk htk Φ + ftk F ⎥ ⎟
                                       ⎦⎠   )   ⎜ ⎢ N tk ( tn −1 ) M ⎥
                                                ⎝⎣
                                                            +
                                                                     ⎦               ⎢ N tk ( tn −1 ) K ⎥
                                                                                     ⎣
                                                                                               +
                                                                                                        ⎦
                                                                                                                      ⎟
                                                                                                                      ⎠.   [3-42]

This equation is of the form

                                                              Kd = F                                                       [3-43]

where the space-time stiffness matrix and force vector are

                   ⎛ ⎛ ⎡ Nt Nt M                    N t1 N tk M ⎤ ⎡ N t1 N t1 K              N t1 N tk K ⎤ ⎞ ⎞
                   ⎜ ⎜⎢ 1 1                                     ⎥ ⎢                                      ⎥⎟ ⎟
                   ⎜∫ ⎜⎢                                        ⎥+⎢                                      ⎥ ⎟ dt ⎟
                   ⎜ In ⎜ ⎢                                                                                     ⎟
                   ⎜ ⎜ ⎣ N tk N t1 M                N tk N tk M ⎥ ⎢ N tk N t1 K              N tk N tk K ⎥ ⎟ ⎟
                                                                                                           ⎟
                        ⎝                                       ⎦ ⎣                                      ⎦⎠
                   ⎜                                                                                            ⎟
                   ⎜ ⎛ ⎡ N t1 ( tn−1 ) N t1 ( tn−1 ) M
                                 +             +
                                                                 N t1 ( tn−1 ) N tk ( tn−1 ) M ⎤ ⎞
                                                                         +             +
                                                                                                                ⎟
                   ⎜ ⎜⎢                                                                         ⎥⎟              ⎟
                K =⎜ ⎜⎢                                                                         ⎥⎟              ⎟          [3-44]
                   ⎜ ⎜⎢                                                                         ⎥⎟              ⎟
                   ⎜ ⎜ ⎢ N tk ( tn−1 ) N t1 ( tn−1 ) M           N tk ( tn−1 ) N tk ( tn−1 ) M ⎥ ⎟
                                 +             +                         +             +
                         ⎣                                                                      ⎦               ⎟
                   ⎜ +⎜                                                                            ⎟            ⎟
                   ⎜ ⎜ ⎡ N t1 ( tn−1 ) N t1 ( tn−1 ) K            N t1 ( tn−1 ) N tk ( tn−1 ) K ⎤ ⎟
                                    +             +                        +             +
                                                                                                                ⎟
                   ⎜ ⎜ ⎢                                                                         ⎥⎟             ⎟
                   ⎜ ⎜+⎢                                                                         ⎥⎟             ⎟
                   ⎜ ⎜ ⎢N t+ N t+ K
                   ⎜ ⎜ ⎢ tk ( n−1 ) t1 ( n−1 )
                                                                                                 ⎥⎟
                                                                  N tk ( tn−1 ) N tk ( tn−1 ) K ⎥ ⎟
                                                                           +             +                      ⎟
                                                                                                                ⎟
                   ⎝ ⎝ ⎣                                                                         ⎦⎠             ⎠

and d is the nodal displacement vector for the current space-time slab. The space-time force

vector is of the form


                                                               27 / 79
                                                        F = F ext + F kin + F int                  [3-45]

with

                                                         ⎡ Nt ( tn−1 ) Mu ( tn−1 ) ⎤
                                                                  +            −

                                                         ⎢   1
                                                                                      ⎥
                                                F kin   =⎢                            ⎥            [3-46]
                                                         ⎢                            ⎥
                                                         ⎢ Ntk ( tn −1 ) Mu ( tn −1 ) ⎥
                                                                  +            −
                                                         ⎣                            ⎦

and

                                                         ⎡ Nt ( tn−1 ) Ku ( tn−1 ) ⎤
                                                                  +           −

                                                         ⎢   1
                                                                                    ⎥
                                                  Fint = ⎢                          ⎥              [3-47]
                                                         ⎢                          ⎥
                                                         ⎢ Ntk ( tn−1 ) Ku ( tn−1 ) ⎥
                                                                  +           −
                                                         ⎣                          ⎦

where F ext is the external space-time force vector, F kin is the inertial or kinetic space-time force

vector, and F int is the internal space-time force vector. In its general form, the external space-

time force vector is

                         F ext = ∫
                                 In ( ⎡⎣ N N [h]
                                          ∫Γh            t1   s                        ) d Γdt
                                                                            N t k N s [ h ]⎤
                                                                                           ⎦
                                                                                               T


                                                                                                   [3-48]
                               +∫ ∫ ( ⎡ N N [ f ]                           N N [ f ]⎤ ) d Ωdt
                                                                                               T

                                      ⎣
                                     In       Ω         t1    s               tk    s⎦

but using the assumption of separation of variables described in Eq. [3-32] and [3-33] it becomes


                                                        ⎜⎢ 1 1
                                                                   (
                                                        ⎛ ⎡ N t ht Φ + f t F ⎤ ⎞
                                                                           1
                                                                                ⎥⎟
                                                                                    )
                                     F    ext
                                                  = ∫ ⎜⎢                        ⎥ ⎟ dt .           [3-49]
                                                     In ⎜
                                                          ⎢                     ⎥⎟
                                                        ⎝⎣         (
                                                        ⎜ ⎢ N tk htk Φ + f tk F ⎥ ⎟
                                                                                ⎦⎠  )
       Note that the K matrix depends on the spatial mass and stiffness matrices and the

temporal shape functions. The F matrix also depends on the spatial mass and stiffness matrices

and the temporal shape functions plus the initial conditions for the current space-time slab (i.e.

the displacements and velocities at the end of the previous time step). Also, the order of the



                                                                  28 / 79
temporal shape functions needs to be at least quadratic. This requirement ensures the existence of

the second time derivative in Eq. [3-44].

3.4.   Temporal Mesh Limitations

       Semi-discrete analysis using explicit integration methods presents a limitation to the size

of the time step used. The solutions to these problems are stable as long as the time step is

limited to a critical upper bound:

                                                                 l
                                            Δt ≤ Δt critical =                             [3-50]
                                                                 c

where c is the wave propagation speed (speed of sound) in the material and l is the

characteristic element length. The speed of sound is

                                                        E
                                                 c=                                        [3-51]
                                                        ρ

where E is Young’s modulus of elasticity and ρ is the volumetric mass density of the material.

The limitation implied by Eq. [3-50] is called the CFL condition after Courant, Friedrichs, and

Levy (Cook, Malkus, Plesha, & Witt, 2001, pp. 411-414). Basically, the explicit method requires

that the wave front should not skip any spatial element during the numerical integration. In

contrast to the semi-discrete formulation, the space-time finite element formulations are

unconditionally stable for linear elastodynamics (Hulbert & Hughes, 1990, p. 331). Thus, any

time step can be used with this formulation without becoming unstable.




                                                  29 / 79
Chapter 4.             Mesh Refinement and Convergence Rate

        Discretization error is an important measure of the robustness of finite element analysis.

It is defined as the difference between an analytical (mathematical) model and its discretized

(finite element) model. The refinement of the mesh (also referred to as h refinement), which

reduces (refines) the characteristic length of an element, should lead to smaller discretization

error. Semi-discrete methods often employ displacement residual L2 error estimates and global

energy error norms. Error norms are typically of the form

                                          error = d app − d exact                                                     [4-1]
                                                                                2


where

                                                                        1

                                                         ⎛ n    k ⎞
                                                                    k

                                               a       = ⎜ ∑ ai ⎟ .                                                   [4-2]
                                                         ⎝ i =1   ⎠
                                                   k




The L2 norm of the displacement on a domain is usually calculated as


                                                       (∫                               )
                                                                                            1

                                    u ( x)                  u i ( x ) ui ( x ) d Ω
                                                                                            2
                                                   =                                                                  [4-3]
                                              L2        Ω



Similar estimates can be employed in space-time formulations.

        Wiberg and Li (1999, pp. 1788-1789) use a Zienkiewicz-Zhu error in the total energy

norm to measure the spatial discretization error at specific times:

                                                                                                              1
                            ⎡ ⎛ ⎡ *⎤
                                      (                      )     ( ⎡ε ⎤ − [ε ] )                  d Ω ⎞⎥
                                                                                                          ⎤       2
                                         ε ⎦ i − [ε ]i [ E ]
                                                       T

                            ⎢ ⎜ ∫Ω ⎣
                                                                            *
                               ne
                                                                     ⎣ ⎦                                ⎟
                       es = ⎢∑ ⎜
                                                                                    i           i
                                                                                                        ⎟⎥            [4-4]
                                          (
                            ⎢ i =1 ⎜ + ∫ ⎡u * ⎤ − [u ]i ρ        ) ( ⎡u ⎤ − [ u ] )
                                                         T
                                                                            *
                                                                                                    d Ω ⎟⎥
                            ⎣ ⎝ Ω ⎣ ⎦i                               ⎣ ⎦        i               i       ⎠⎦

or equivalently




                                                        30 / 79
                                                                                                                        1
                               ⎡ ⎛ ⎡ *⎤
                                                  (                 ) (           ⎤ 2
                                       σ ⎦ i − [σ ]i [ E ] ⎡σ * ⎤ i − [σ ]i d Ω ⎞ ⎥                         )
                                                    T     −1
                               ⎢ ⎜ ∫Ω ⎣
                                  ne
                                                              ⎣ ⎦               ⎟
                          es = ⎢ ∑ ⎜                                            ⎟⎥ .                                                    [4-5]
                                                      (              ) (
                               ⎢ ⎜ + ∫ ⎡u ⎤ − [u ]i ρ ⎡u ⎤ − [u ]i d Ω                              )
                                                     T
                                 i =1    *                   *
                                                                                ⎟⎥
                               ⎣ ⎝ Ω ⎣ ⎦i                 ⎣ ⎦i                  ⎠⎦

In these norms, ⎡ε * ⎤ , ⎡σ * ⎤ and ⎡u * ⎤ are a smoothed (or analytical) strain, stress and velocity
                ⎣ ⎦i ⎣ ⎦            ⎣ ⎦

field and the strains, stresses and velocities calculated from the finite element analysis are

[ε ] = [ B ][ d ] , [σ ] = [ E ][ε ] and [u ] = ⎡ N ⎤ [ d ] , respectively. The velocities usually converge at a
                                                ⎣ ⎦

faster rate than the stresses, so the kinetic energy terms are often not included in actual

implementation. The kinetic energy terms are needed though when the strain energy is very

small or the kinetic energy dominates the error measurement (Wiberg & Li, 1999, p. 1789).

        The L2 residual error is defined as

                                                                                                                                    1
             u h ( x, t ) − u a ( x, t )         = ⎡ ∫ ( u h ( x, t ) − u a ( x, t ) ) ( u h ( x, t ) − u a ( x, t ) ) d Ω ⎤ .
                                                                                                                             2
                                                                                                                                        [4-6]
                                            L2     ⎣ Ω                                                                     ⎦

The L2 norm defined by Eq. [4-6] provides a measure on the mean error over all degrees of

freedom (Belytschko, Liu, & Moran, 2000, p. 332)

         Li and Wiberg (1998, p. 216) also use a total energy norm along the temporal jumps in

the displacements and velocities to estimate the error in the time discretization.

                                                                                            1

                             ⎣
                                  h
                                    (
                        et = ⎡ a un , un
                                       h
                                                          ) +( v
                                                          Ω
                                                                    h
                                                                    n   ,ρ v   h
                                                                               n   )    ⎤
                                                                                       Ω⎦
                                                                                                2


                                                                                                                1
                             = ⎡ ∫ ∇un ⋅ σ ∇un d Ω + ∫ vn ⋅ ρ vn d Ω ⎤
                                                                                                                    2
                                     h       h          h      h
                                                                                                                                        [4-7]
                               ⎣ Ω                    Ω              ⎦

                                            (                        ) (                                            )
                                                                                                                            1
                               ⎡ B u h ( tn ) − u h ( tn ) ⋅ E ⋅ B u h ( tn ) − u h ( tn ) d Ω ⎤
                               ⎢ ∫Ω
                                           +            −                 +            −                                        2

                             =                                                                 ⎥
                               ⎣ ∫Ω      n (
                               ⎢+ vh (t + ) − vh (t − ) ⋅ ρ vh (t + ) − vh (t − ) d Ω
                                                      n            n) (         n
                                                                                               ⎥
                                                                                               ⎦        )
Their results show that the discontinuous Galerkin method is third order accurate in time.




                                                                  31 / 79
        Hughes and Hulbert (1988, pp. 350-351) have proposed a different norm for space-time

formulations including least squares terms that is stronger than the traditional norms used in

semi-discrete analysis. When simplified for single field space-time formulations, this norm is

given as:

                                                                                   N −1
                         ||| wh |||2 =E   (          )
                                              wh (T − ) +E     (               )
                                                                   wh ( 0+ ) + ∑E
                                                                                   n =1
                                                                                          (   wh ( t n )   )   [4-8]


where

                                  E   (w ) = 1 (w , ρw )
                                             2
                                              h          h         h
                                                                       Ω
                                                                            1
                                                                           + a ( wh , wh ) .
                                                                            2             Ω
                                                                                                               [4-9]


        In addition to the h refinement, another method of refinement, known as p refinement,

changes the element type by increasing the highest complete polynomial in the element field

quantity. A simple example of p refinement is changing a typical one-dimensional bar finite

element from a linear Lagrangian element with two degrees of freedom (two nodes) to a

quadratic Lagrangian element with three degrees of freedom (three nodes).


        The error estimates provided above will be implemented in the sample problems in

Chapter 6 to test the numerical performance of the method. The resulting log-log plots show the

relationship between the mesh size and the error for convergence using h refinement. Refinement

methods that converge well should appear roughly linear in these log-log plots with positive

slopes (i.e. an increase in the characteristic element length increase the error). The slopes of these

straight lines are proportional to the rate of convergence (Belytschko, Liu, & Moran, 2000, p.

459). A greater rate of convergence indicates that the finite element solution will converge to an

exact solution faster.




                                                             32 / 79
Chapter 5.                   Analytical Space-Time Matrix Formulation

       In this chapter, the analytical form of the space-time formulation is derived for certain

choices of space and time shape functions. In finite element analysis, shape functions are used to

interpolate field quantities. In the case of Lagrangian shape functions, it is possible to derive the

analytical form of the space-time formulation. The Lagrangian shape function has the following

general pattern:


                                  Ni =
                                           ( x1 − x )( x2 − x ) [ xi − x ] ( xk − x )                                                                 [5-1]
                                         ( x1 − xi )( x2 − xi ) [ xi − xi ] ( xk − xi )

for a finite element having k nodes and the bracketed terms are omitted from the products in the

numerator and denominator to obtain the i th shape function (Cook, Malkus, Plesha, & Witt,

2001, p. 86). Lagrangian shape functions are useful in finite element analysis because the

simplicity of its formula makes them easy to program for any size and dimension of element.

5.1.   Quadratic Lagrangian Temporal Shape Functions

       For the nt
                       th
                            space-time slab with temporal nodes at tnt , tnt − 12 = 1 tnt −1 + tnt , and
                                                                                    2                                           (                 )
tnt −1 = tnt − δ t , where δ t is the time step used, let the quadratic temporal shape functions be


          Nt = ⎢
                       (
               ⎡ tn − 1 − t tn − t
                       t 2
                                 )(
                                  t
                                           )            (t   nt −1      )(
                                                                     − t tnt − t   )        (t    nt −1       )(
                                                                                                          − t tnt − 12 − t          )        ⎤
                                                                                                                                             ⎥.       [5-2]
                   (             )(
               ⎢ tn − 1 − tn −1 tn − tn −1
               ⎣ t 2 t            t    t
                                               ) (t   nt −1 − tnt − 1 2 )(
                                                                        tnt − tnt − 12   ) (t   nt −1 − tnt   ) (t   nt   −1   2
                                                                                                                                 − tnt   )   ⎥
                                                                                                                                             ⎦

Figure 2 shows the quadratic temporal shape functions normalized between 0 and 1 time unit.




                                                           33 / 79
                     1.2
                                                                                                 Nt1
                         1                                                                       Nt2
                                                                                                 Nt3
                     0.8

                     0.6

                     0.4

                     0.2

                         0

                    -0.2
                        0           0.2           0.4        0.6                   0.8                 1
                                                 Normalized Time

                  Figure 2. Normalized Quadratic Temporal Shape Functions



The overall shape functions of the space-time slab are

                                  N ( x, t ) = ⎡ Nt1 N x
                                               ⎣                  N t2 N x   Nt3 N x ⎤
                                                                                     ⎦                       [5-3]

and the space-time nodal displacement vector is

                                                     ⎡ xn −1 ⎤
                                                     ⎢ t ⎥
                                                 d = ⎢ xnt − 1 2 ⎥ .                                         [5-4]
                                                     ⎢           ⎥
                                                     ⎢ xnt ⎥
                                                     ⎣           ⎦

The derivatives of the temporal shape functions are

                          1 ⎡
                  Nt =
                         δt2 ⎣
                              (           )
                               δ t − 4tnt + 4t       (8t   nt          )
                                                                − 4δ t − 8t     ( −4t   nt         )
                                                                                             + 3δ t + 4t ⎤
                                                                                                         ⎦
                                                                                                             [5-5]


                                                      1
                                              Nt =        [ 4 −8 4 ] .                                       [5-6]
                                                     δ t2

The global space-time stiffness matrix becomes




                                                       34 / 79
                 ⎛ ⎡ Nt Nt M           N t1 N t2 M     N t1 N t3 M ⎤ ⎡ N t1 N t1 K       N t1 N t2 K     N t1 N t3 K ⎤ ⎞
                 ⎜⎢ 1 1                                            ⎥ ⎢                                               ⎥⎟
         K = ∫ ⎜ ⎢ N t2 N t1 M         N t2 N t2 M     N t2 N t3 M ⎥ + ⎢ Nt2 N t1 K      N t2 N t2 K     N t2 N t3 K ⎥ ⎟ dt
              In
                 ⎜⎢                                                ⎥ ⎢                                               ⎥⎟
                 ⎜ ⎢ Nt Nt M           N t3 N t2 M     N t3 N t3 M ⎥ ⎢ N t3 N t1 K       N t3 N t2 K                   ⎟
                                                                                                         N t3 N t3 K ⎥ ⎠
                 ⎝⎣ 3 1                                            ⎦ ⎣                                               ⎦
               ⎡ N t ( tn−1 ) Nt ( tn−1 ) M
                         +             +
                                                  N t1 ( tn−1 ) Nt2 ( tn−1 ) M
                                                          +            +
                                                                                  N t1 ( tn−1 ) N t3 ( tn−1 ) M ⎤
                                                                                          +             +

               ⎢    1             1
                                                                                                                ⎥
             + ⎢ N t2 ( tn−1 ) N t1 ( tn−1 ) M
                         +             +
                                                 N t2 ( tn−1 ) N t2 ( tn−1 ) M
                                                         +             +
                                                                                  N t2 ( tn−1 ) N t3 ( tn−1 ) M ⎥
                                                                                          +             +
                                                                                                                              .    [5-7]
               ⎢                                                                                                ⎥
               ⎢ N t ( tn−1 ) N t ( tn−1 ) M
                         +             +
                                                 Nt3 ( tn−1 ) N t2 ( tn−1 ) M
                                                        +             +
                                                                                  N t3 ( tn−1 ) N t3 ( tn−1 ) M ⎥
                                                                                          +             +
               ⎣ 3                1                                                                             ⎦
               ⎡ N t ( tn−1 ) N t ( tn−1 ) K
                         +             +
                                                 N t1 ( tn−1 ) N t2 ( tn−1 ) K
                                                         +             +
                                                                                 N t1 ( tn−1 ) N t3 ( tn−1 ) K ⎤
                                                                                         +             +

               ⎢    1             1
                                                                                                               ⎥
             + ⎢ N t2 ( tn−1 ) N t1 ( tn−1 ) K
                         +             +
                                                 N t2 ( tn−1 ) N t2 ( tn−1 ) K
                                                         +             +
                                                                                 Nt2 ( tn−1 ) N t3 ( tn−1 ) K ⎥
                                                                                         +             +

               ⎢                                                                                               ⎥
               ⎢ N t ( tn−1 ) N t ( tn−1 ) K
                         +             +
                                                 N t3 ( tn−1 ) N t2 ( tn−1 ) K
                                                         +             +
                                                                                 N t3 ( tn−1 ) Nt3 ( tn−1 ) K ⎥
                                                                                         +             +
               ⎣ 3                1                                                                            ⎦

By substituting the shape functions, the matrix reduces to

                        ⎡ 5                     4              1     ⎤ ⎡ 1                 2            1 ⎤
                        ⎢ δt2 M              −     M        −       M⎥ ⎢ K                − K             K
                                               δt2            δt  2
                                                                         2                 3            6 ⎥
                        ⎢                                            ⎥ ⎢                                    ⎥
                           12                 16               4         2                               2 ⎥
                    K = ⎢− 2 M                    M         − 2 M⎥+ ⎢ K                     0K         − K                         [5-8]
                        ⎢ δt                  δt2             δt     ⎥ ⎢ 3                               3 ⎥
                        ⎢                                            ⎥ ⎢                                    ⎥
                        ⎢ 72 M                 12
                                             − 2M
                                                              5           1
                                                                   M ⎥ ⎢− K
                                                                                            2
                                                                                              K
                                                                                                        1 ⎥
                                                                                                          K
                        ⎢ δt
                        ⎣                      δt            δt 2
                                                                     ⎥ ⎢ 6
                                                                     ⎦ ⎣                    3           2 ⎥ ⎦

in which the spatial global mass and stiffness matrices are

                                                      M = ∫ N sT ρ N s d Ω                                                         [5-9]
                                                               Ω



                                                     K = ∫ BsT [C ] Bs d Ω .                                                      [5-10]
                                                           Ω


The global space-time forcing vector becomes

                   ⎛ ⎡ N t ( ht Φ + f t F ) ⎤ ⎞      ⎡ N t ( tn −1 ) M u ( tn −1 ) ⎤ ⎡ N t ( tn −1 ) K u ( tn −1 ) ⎤
                                                              +             −                   +             −

                   ⎜⎢     1
                                            ⎥⎟       ⎢ 1                           ⎥ ⎢ 1                             ⎥
           F = ∫ ⎜⎢                         ⎥ ⎟ dt + ⎢ N t ( tn −1 ) M u ( tn −1 ) ⎥ + ⎢ N t ( tn −1 ) K u ( tn −1 ) ⎥ .
                                                              +             −                   +             −
                                                                                                                                  [5-11]
                In
                   ⎜ ⎢N h Φ + f F ⎥ ⎟                ⎢ 2                           ⎥ ⎢ 2                             ⎥
                   ⎜ tk ( t               ) ⎦ ⎟ ⎢ N t ( tn+−1 ) M u ( tn−−1 ) ⎥ ⎢ N t ( tn+−1 ) K u ( tn−−1 ) ⎥
                   ⎝⎣                 t
                                              ⎠      ⎣ 3                           ⎦ ⎣ 3                             ⎦

Assuming constant body forces and tractions in time, the space-time force vector becomes




                                                               35 / 79
                    ⎡ 3                  ⎤
                    ⎢ − δ t Mu ( tn −1 ) ⎥ ⎡
                                     −


                                         ⎥ ⎢ Ku ( tn −1 ) ⎥
   ⎡− ( Φ + F )⎤ ⎢
                                                     −
                                                          ⎤
   ⎢            ⎥ ⎢ 4
                ⎥ + ⎢ δ t Mu ( tn −1 ) ⎥ + ⎢ 0 ⎥
F=⎢       0                        −     ⎥
                                            ⎢             ⎥
   ⎣ (Φ + F ) ⎦ ⎢ 1
   ⎢            ⎥                        ⎥ ⎢ 0 ⎥
                                            ⎣             ⎦
                    ⎢ − Mu ( tn −1 ) ⎥
                                     −

                    ⎢ δt
                    ⎣                    ⎥
                                         ⎦
                                                            .   [5-12]
   ⎡                3                              ⎤
   ⎢ − ( Φ + F ) − δ t Mu ( tn −1 ) + Ku ( tn −1 ) ⎥
                              −             −


   ⎢                                               ⎥
                   4
 = ⎢                  Mu ( tn −1 )
                             −                     ⎥
   ⎢              δt                               ⎥
   ⎢                                               ⎥
                         1
   ⎢        ( Φ + F ) − Mu ( tn−−1 )               ⎥
   ⎢
   ⎣                    δt                         ⎥
                                                   ⎦




                          36 / 79
Chapter 6.             Implementation of Space-Time Formulation

       Based on the analytical formulation described in Chapter 5, a set of MATLAB codes

have been developed. The codes developed to implement the space-time formulation were

constructed in MATLAB using the Symbolic toolbox. All integrals are calculated using Gauss-

Legendre quadrature. The use of quadrature helps account for the time-varying traction applied

to the structures. Thus, a simulation using a relatively coarse mesh can account for the dynamic

behavior of the loading function by using a large number of quadrature points. The inputs to the

program include the structure geometry, physical properties, loading conditions, mesh

parameters, the order of the quadrature to be used, and an option to perform a convergence study.

The output of the codes includes the displacement, strain, stress, strain energy, kinetic energy,

and work histories as well as the error estimates used in the convergence studies. The use of the

code is demonstrated by evaluating four example structures: an axially loaded bar, a transversely

loaded straight beam, a transversely loaded thick curved beam (or arch), and a transversely

loaded thin curved beam. The specific implementations are outlined below.


6.1.   Axially Loaded Bar

       An axially loaded bar problem can been successfully solved analytically using the

techniques of classical structural dynamics and acoustics. This problem therefore makes a good

test case to compare the results of a one-dimensional space-time finite element simulation with

an analytical solution as well as a semi-discrete finite element solution.

6.1.1. General Analytical Solution – Modal Superposition Method

       Consider a uniform straight bar, without damping, initially undisplaced and at rest as

shown in Figure 3. The bar is considered to be continuous with distributed parameters. The left

hand side of the bar (at x = 0 ) is fixed rigidly. The axial stiffness is EA and the mass per unit
                                               37 / 79
length is m where E is Young’s modulus of elasticity and A is the cross section area of the bar.

The general loading applied to the bar is q ( x , t ) representing an arbitrary external distributed

axial load. The internal time-varying axial force distribution is P ( x, t ) . The end boundary

conditions are those used in the numerical example, but the analytical derivation could use

arbitrary boundary conditions. The forces acting on a differential section are shown in Figure 4.




                 AE, m(x)


P(0,t)                                                                       P(L,t)          x, u

                                                        q(x,t)
                       x               dx
                                         L
               Figure 3. Bar Subject to Dynamic Axial Load, Analytical Example




                                              38 / 79
                                                                                 ∂ 2u ( x, t )
                                                                    m ( x, t )                 dx
                                                                                     ∂t 2


                                                                                    ∂P ( x, t )
                        P ( x, t )                                   P ( x, t ) +               dx
                                                                                      ∂x



                                                                   q ( x, t ) dx




                                     x           dx


         Figure 4. Section of Bar Subject to Dynamic Axial Load, Analytical Example

       The axial displacement can be derived from the equation of motion as shown in general

by Clough and Penzien (1993, pp. 373-375, 391-393, 407-411). The resulting circular natural

frequencies from their analysis are

                                                            π c
                                           ωn = ( 2n − 1)                                            [6-1]
                                                            2 L

where the speed of sound in the material is

                                                       EA
                                                c=        .                                          [6-2]
                                                       m
The corresponding normal mode shapes are

                                                      ⎛           π x⎞
                                      Φ n ( x ) = sin ⎜ ( 2n − 1)    ⎟.                              [6-3]
                                                      ⎝           2 L⎠




                                                  39 / 79
These mode shapes are orthogonal to each other (Clough & Penzien, 1993, pp. 393-394). The

complementary solution is

                                   ∞
                                        ⎛          π x⎞
                     uc ( x, t ) = ∑sin ⎜ ( 2n −1)    ⎟ ( Ac,n cos (ωnt ) + Bc,n sin (ωnt ) )                  [6-4]
                                   n=1  ⎝          2 L⎠

where Ac,n and Bc,n are determined by the initial conditions. Let the loading on the bar be a

suddenly applied axial compressive force at the tip of magnitude P . The particular solution of
                                                                  0


the bar is

                                                             2
                              2P ∞      n⎛          π⎞     ⎛           π x⎞
                u p ( x, t ) = 0 ∑ ( −1) ⎜ ( 2n − 1) ⎟ sin ⎜ ( 2n − 1)    ⎟ (1 − cos (ωnt ) ) .                [6-5]
                              EA n=1     ⎝          2⎠     ⎝           2 L⎠

The general solution is the superposition of the complementary and particular solutions:

                                                  ⎛⎛                                   ⎞ ⎞
                                                                                         2
                                                                                                         ⎞
                                                                    n 2P ⎛      2
                                                  ⎜ ⎜ Ac ,n − ( −1)     0
                                                                          ⎜            ⎟   ⎟ cos (ωn t ) ⎟
                           ∞
                                ⎛           π x ⎞⎜⎝ ⎜                 EA ⎜ ( 2n − 1) π ⎟ ⎟
                                                                          ⎝            ⎠ ⎠               ⎟
             u ( x, t ) = ∑ sin ⎜ ( 2n − 1)     ⎟ ⎜                                                      ⎟.    [6-6]
                                            2 L ⎠⎜
                                ⎝                                                                  ⎞ ⎟
                                                                                                     2
                                                                              n 2P ⎛
                          n =1
                                                  ⎜ + B sin (ω t ) + ( −1)         0        2            ⎟
                                                                                     ⎜
                                                                                     ⎜             ⎟ ⎟
                                                                                                   ⎟
                                                  ⎜ c ,n                        EA ⎝ ( 2n − 1) π ⎠ ⎠
                                                                   n
                                                  ⎝

The initial conditions are

                                                  u ( x, 0 ) = ρ ( x )                                         [6-7]

                                                  u ( x, 0 ) = ψ ( x ) .                                       [6-8]

For 0 ≤ x ≤ L

                                              ∞
                                                        ⎛          π x⎞
                                  u ( x,0) = ∑ Ac,n sin ⎜ ( 2n −1)    ⎟ = ρ ( x)                               [6-9]
                                             n=1        ⎝          2 L⎠

                                                      ⎡                              ⎤
                                                      ⎢                              ⎥
                                                                         ( −1)
                                                                               n
                              ∞
                                    ⎛           π x ⎞⎢           2P
                 u ( x, 0 ) = ∑ sin ⎜ ( 2n − 1)     ⎟  ωn Bc ,n + 0                  ⎥ =ψ ( x) .              [6-10]
                                    ⎝           2 L⎠⎢            EA ⎛           π⎞ ⎥
                                                                                   2

                                                                    ⎜ ( 2n − 1) ⎟ ⎥
                              n =1
                                                      ⎢
                                                      ⎣             ⎝            2⎠ ⎦


                                                        40 / 79
Assuming that the bar is initially undeformed ( ρ ( x ) = 0 ) and at rest (ψ ( x ) = 0 ), then

                                                    Ac,n = 0                                     [6-11]


                                                         ( −1)
                                                                   n
                                               2P
                                     Bc ,n   =− 0                2
                                                                   .                             [6-12]
                                               EA ⎛          π⎞
                                                  ⎜ ( 2n − 1) ⎟ ωn
                                                  ⎝           2⎠

Inserting Eq. [6-11] and [6-12] into Eq. [6-6] gives

                                     ( −1) sin ⎛ 2n − 1 π x ⎞ ⎛1 − cos ω t − sin (ωnt ) ⎞
                                        n
                         2P ∞
             u ( x, t ) = 0 ∑                    (     ) ⎟⎜           ( n)              ⎟.       [6-13]
                         EA n=1 ⎛          π⎞ ⎜⎝        2 L ⎠⎝                  ωn ⎠
                                ⎜ ( 2n − 1) ⎟
                                ⎝          2⎠

6.1.2. General Analytical Solution – Acoustic Wave Propagation Method

        Another way of looking at the analytical solution is in the wave propagating behavior of

the stress and displacement. Clough and Penzien (1993, pp. 411-423) have also derived the wave

solution to a general axially loaded bar problem and the presentation here follows their example.

        The bar is initially unstressed and undisplaced until the force is applied at the free end

producing a compressive or tensile stress. The displacement and stress move as a wave through

the bar at the wave speed in Eq. [6-2] which is equivalent to

                                                            E
                                                    c=                                           [6-14]
                                                            ρ

where ρ is the volumetric mass density. The equation of motion becomes (Kinsler, Frey,

Coppens, & Sanders, 2000, p. 70)

                                         ∂ 2u ( x, t ) 1 ∂ 2u ( x, t )
                                                      = 2              .                         [6-15]
                                            ∂x 2       c    ∂t 2

The general solution has the form

                                    u ( x, t ) = f1 ( x − ct ) + f 2 ( x + ct )                  [6-16]


                                                     41 / 79
where f1 and f2 are arbitrary functional relationships of the parameters x − ct and x + ct . This

form of the solution denotes a pair of displacement waves propagating in the positive and

negative directions, along the axis of the bar. At the initial time t = 0 , the two waves are

functions of only position.

       Looking at the instant of time at t = Δt , the new position variable is x′ = x − ct . Thus

                                            f1 ( x − cΔt ) ≡ f1 ( x′ ) .                       [6-17]

The shape of the wave in relation to new coordinate x′ remains unchanged compared to the

shape relative to the initial coordinate x . Basically, the wave moves a distance of cΔt at a speed

of c during the time increment Δt without changing its shape. Similar arguments can show that

f2 in Eq. [6-16] is a waveform moving at speed c in the negative x -direction.

       The wave behavior of the bar can be extended to the stress distribution of the bar. The

basic strain-displacement relationship for one-dimensional problems is

                                                         ∂u
                                                    ε=      .                                  [6-18]
                                                         ∂x

The stress-strain relationship of Hooke’s Law reduces to

                                                   σ = Eε .                                    [6-19]

The stress wave of the bar becomes

                                           ∂u    ∂f               ∂f
                          σ ( x, t ) = E      = E 1 ( x − ct ) + E 2 ( x + ct ) .              [6-20]
                                           ∂x    ∂x               ∂x

                              ∂f1      ∂f
The stress wave functions E       and E 2 can be designated by g1 and g2 . Thus
                              ∂x       ∂x

                                  σ ( x, t ) = g1 ( x − ct ) + g 2 ( x + ct ) .                [6-21]

The stress wave also propagates with the velocity c and with unchanging shape.



                                                    42 / 79
       The shape of the wave propagating through a uniform bar is changed by the boundary

conditions at the ends of the bar. These boundary conditions impose equilibrium and

compatibility requirements on the waves. For the example shown in Figure 3, the right end of the

bar, at x = L , is free to move. This boundary condition results in the stress being zero at all times

at that material coordinate. This condition is accomplished by another stress wave being

superposed on top of the incident wave but propagating in the opposite (left) direction and in the

opposite sense. Thus the two cancel each other out at the end of the bar. This can be expressed

by the following equation

       If the right end ( x = L ) of the bar is free, as indicated in Figure 3, the condition of zero

stress must be maintained at all times at that end. This condition may be satisfied by a second

stress wave propagating toward the left, which when superposed on the incident wave, cancels

the end-section stresses. Expressing this concept mathematically by means of Eq. [6-20] leads to

                                               ∂f1               ∂f
                              σ x= L = 0 = E       ( L − ct ) + E 2 ( L + ct )                 [6-22]
                                               ∂x                ∂x

from which

                                     ∂f1               ∂f
                                         ( L − ct ) = − 2 ( L + ct ) .                         [6-23]
                                     ∂x                ∂x

                                                 ∂u
       This leads to the fact that the slope        of the left propagating wave is the negative slope
                                                 ∂x

of the right propagating wave as each part of the waves passes the end of the bar. Thus the two

cancel each other out at the end. When the wave front hits the free end, the stress is reflected but

in the opposite sense, so the stress switches between compressive and tensile, or vice versa. The

displacement wave is also reflected but the displacements of the reflected wave are the same

sense as the incident wave. Thus, the displacement is doubled by the superposition of the

reflected and incident waves at the free end while the stress is canceled.

                                                   43 / 79
        The other (left) side of the bar shown in Figure 3 is fixed rather than free, so the wave

reflections react differently to this boundary condition. The displacement wave has the form

                                 u x = 0 = 0 = f1 ( 0 − ct ) + f 2 ( 0 + ct ) .                [6-24]

The reflected wave (now moving in the positive direction) can be expressed by the incident wave

terms

                                       f1 ( 0 − ct ) = − f 2 ( 0 + ct ) .                      [6-25]

Similar to the case of stress at the fixed boundary condition, the displacement waves have

opposite signs. The incident and reflected stress waves at the fixed end, analogous to the

displacement waves at the free end, have the same sense. Thus, the fixed end has a canceling

effect on the displacement and a doubling effect on the stress.

6.1.3. Numerical Solution

        Two one-dimensional space-time finite element simulations were run for this example

using different time steps. The bar used in these examples is shown in Figure 5 and Figure 6.



                      A
                                                                              q=(7.5kN )H(t)
                                                                                     m            x

                      A

                        L=10m

               Figure 5. Bar Subject to Dynamic Axial Load, Numerical Example




                                                   44 / 79
                                         H =1m


                                                           T = 1m
   Figure 6. Cross Section A-A of Bar Subject to Dynamic Axial Load, Numerical Example



H ( t ) is the Heaviside step function in Figure 5. The material properties of the bar are:

volumetric mass density is ρ = 7860 kg m3 ; Young’s modulus of elasticity is E = 200GPa ; and

Poisson’s ratio is ν = 0.3 . The 2D finite element models are treated as a plane stress problem.

        All simulations use 75 elements to discretize the bar along its length. The temporal

domain is discretized using a time step of Δt = 5.29 × 10−5 s (with 150 slabs along the time axis)

in one simulation and a time step of Δt = 2.65 × 10−4 s (with 30 slabs along the time axis) in the

other space-time finite element simulation. The critical time step for the central difference

method in this problem is Δt critical = 2.65 ×10−5 s , thus the small time step used is twice the critical

time step and the large time step is ten times the critical time step. Quartic temporal shape

functions and quadratic spatial shape functions (3-node bar elements) were used. The

isoparametric form of the spatial shape functions is

                    N s = ⎡ 1 ( 0 − ξ )(1 − ξ ) − ( −1 − ξ )(1 − ξ )
                          ⎣2
                                                                       1
                                                                       2   ( −1 − ξ )( 0 − ξ )⎤
                                                                                              ⎦    [6-26]

with nodes at


                                             [ξ ] = [ −1   0 1] .                                  [6-27]

The isoparametric form of the temporal shape functions is




                                                   45 / 79
                                                       ⎡ 16 ( − 1 − ξ )( 1 − ξ ) (1 − ξ ) ⎤
                                                                                                  T
                                                            9
                                                                 3         3
                                                       ⎢ 27                                   ⎥
                                                       ⎢ − 16 ( −1 − ξ ) ( 3 − ξ ) (1 − ξ ) ⎥
                                                                             1
                                                  N t = 27                                                 [6-28]
                                                       ⎢ 16 ( −1 − ξ ) ( − 1 − ξ ) (1 − ξ ) ⎥
                                                                             3
                                                       ⎢ 9                                    ⎥
                                                       ⎢ − 16 ( −1 − ξ ) ( − 3 − ξ )( 3 − ξ ) ⎥
                                                                              1       1
                                                       ⎣                                      ⎦

with nodes at

                                                          [ξ ] = [ −1    −1
                                                                          3
                                                                                1
                                                                                3   1] .                   [6-29]

       The spatial finite element mesh is shown in Figure 7. Circles indicate nodes where

tractions are applied and squares indicate nodes where boundary conditions are applied. The

simulations were run for a total time period of 8 ×10−3 s so that the wave front could make one

complete cycle.




                                             3

                                             2
                  Initial Y-position, in m




                                             1

                                             0

                                             -1

                                             -2

                                             -3


                                              0      2              4            6                8   10
                                                              Initial X-position, in m

                                                    Figure 7. 1D Bar Spatial Mesh

6.1.4. Convergence Results

       The convergence of the spatial and temporal domains was simulated separately. For the

spatial convergence, 529 space-time slabs (with a time step of Δt = 15. × 10−6 s ) were used in


                                                                    46 / 79
each step of the convergence. The spatial domain was refined using 22+h elements with

characteristic element lengths of 10 × 1
                                       2  (           2+ h
                                                             ) m . Four simulations were run with differing types of
elements. Linear (SL) and quadratic (SQ) elements were used in the spatial discretization while

quadratic (TQ) and cubic (TC) were used in the temporal discretization.


       The strain energy norm es is shown in Figure 8 at time t = 4.05 ×10−3 s . The slopes of

the curves in these logarithmic plots are approximately 0.435.


                              -4
                            10




                              -5
                  ||e s||




                            10




                                                                                              SL-TQ
                                                                                              SL-TC
                                                                                              SQ-TQ
                              -6
                                                                                              SQ-TC
                            10 -0.8           -0.6               -0.4              -0.2               0
                              10           10                   10               10               10
                                                Characteristic Spatial Element Length, in m



                                 Figure 8. Strain Energy Norm at 0.00405 seconds




                                                               47 / 79
                             -8
                            10




                             -9
               ||e|| L 2    10




                             -10
                            10
                                                                                       SL-TQ
                                                                                       SL-TC
                                                                                       SQ-TQ
                                                                                       SQ-TC

                             -11
                            10        -0.8             -0.6     -0.4            -0.2            0
                                 10                10          10          10              10
                                              Characteristic Spatial Element Length, in m


                                             Figure 9. L 2 Residual Displacement Norm

       The energy norm u h                        used by Hulbert and Hughes below shows convergence with

slopes about 0.680.


                                 -3
                            10




                                               SL-TQ
                                               SL-TC
                                               SQ-TQ
                                               SQ-TC
                |||uh|||2




                                 -4
                            10




                                 -5
                            10 -0.8                    -0.6     -0.4        -0.2               0
                              10                   10          10          10              10
                                              Characteristic Spatial Element Length, in m




                                                              48 / 79
                                      Figure 10. Total Energy Norm




                         -4
                        10




                         -5
                        10
              ||e t||




                         -6
                        10



                                                                                   SL-TQ
                         -7                                                        SL-TC
                        10
                                                                                   SQ-TQ
                                                                                   SQ-TC
                                                 -2
                                               10
                                    Characteristic Temporal Element Length (cΔt), in m



                              Figure 11. Total Energy Across Temporal Jump

       Although a convergence study specifically on p refinement is not directly performed in

this thesis, the plots in Figure 9 and Figure 11 do provide an indication of convergence using p

refinement. The slopes, and thus the rates of convergence, increase from the lower order element

types to the higher order element types.


6.1.5. Numerical Results

       The axial displacement in the X-direction is shown in the meshes below in Figure 12 and

Figure 13.




                                                      49 / 79
                Figure 12. 1D Bar Axial Displacement (Δt = 5.29×10 -5 s), in m




                Figure 13. 1D Bar Axial Displacement (Δt = 2.65×10 -4 s), in m

       A comparison between the analytical model, the 1D space-time finite element models,

and a semi-discrete explicit finite element model using the central difference method is shown in

Figure 14 and Figure 15. The axial displacement of the free end of the bar is shown in Figure 14.

                                             50 / 79
Each method appears to accurately capture the displacement of the bar. Figure 15 shows the axial

stress at the midpoint of the bar.


                                                                      -7
                                                                   x 10
                                                              0


                                                              -1                   Analytical Model
                                                                                                                -5
                                                                                   Space-Time Model (Δt=5.29×10 s)
                                                                                                                -4
                                                              -2                   Space-Time Model (Δt=2.65×10 s)
                                                                                   Semi-Discrete Explicit Model
                                       X-displacement, in m




                                                              -3


                                                              -4


                                                              -5


                                                              -6


                                                              -7


                                                              -8
                                                                0          1   2          3       4         5             6   7          8
                                                                                               Time, in s                            -3
                                                                                                                                  x 10



                                                               Figure 14. Axial Displacement of Free End of Bar, in m



                                            2000
                                                                                     Analytical Model
                                                              0                                                      -5
                                                                                     Space-Time Model (Δt=5.29×10 s)
                                                                                                                     -4
                                       -2000                                         Space-Time Model (Δt=2.65×10 s)
                                                                                     Semi-Discrete Explicit Model
                                       -4000
                Axial stress, in Pa




                                       -6000

                                       -8000

                                      -10000

                                      -12000

                                      -14000

                                      -16000

                                      -18000
                                            0                              1   2          3           4     5             6   7          8
                                                                                              Time, in s                             -3
                                                                                                                                  x 10



                                                                   Figure 15. Axial Stress at Midpoint of Bar, in Pa

                                                                                              51 / 79
       The axially loaded bar problem exemplifies many of the difficulties inherent with semi-

discrete methods. The analytical stress history is discontinuous. The use of the central difference

method (an explicit semi-discrete method) is here limited by a relatively small time step (as

compared to the space-time model) and shows a great deal of oscillation when it encounters a

discontinuity. This is emphasized in the zoomed in plot of the axial stress in Figure 16. The

space-time model lifts this strict limitation on the time step and also exhibits greatly reduced

oscillation when it encounters a discontinuity.



                                       -5000
                                                                Analytical Model
                                       -5500                                                -5
                                                                Space-Time Model (Δt=5.29×10 s)
                                                                                            -4
                                                                Space-Time Model (Δt=2.65×10 s)
                                       -6000
                                                                Semi-Discrete Explicit Model
                 Axial stress, in Pa




                                       -6500

                                       -7000

                                       -7500

                                       -8000

                                       -8500

                                       -9000

                                       -9500

                                                5   5.2   5.4   5.6    5.8     6     6.2   6.4    6.6   6.8          7
                                                                        Time, in s                               -3
                                                                                                              x 10



                                               Figure 16. Axial Stress at Midpoint of Bar, in Pa

6.2.   Transversely Loaded Straight Beam

       The axially loaded bar problem in the previous section is a rather simple example. This

section uses a beam under dynamic transverse loading as a second test case. The numerical finite

element simulation is performed only with a two-dimensional spatial model.




                                                                        52 / 79
6.2.1. Analytical Solution

       Consider a uniform straight slender beam, without damping, initially undisplaced and at

rest as shown in Figure 17. The beam is assumed to be continuous with distributed parameters.

Both sides of the beam (at x = 0 and at x = L ) are fixed rigidly. The flexural stiffness is E I and

the mass per unit length is m where E is Young’s modulus of elasticity and I is the second

moment of inertia of the cross sectional area with respect to the neutral axis through the centroid.

The general loading applied to the beam is p ( x, t ) representing an arbitrary external distributed

transverse load. The transverse displacement response is v ( x , t ) . The end boundary conditions

are those used in the numerical example, but the analytical derivation could use arbitrary

boundary conditions. The forces acting on a differential section are shown in Figure 18.


  y ( x, t )
                                         p ( x, t )




                                                                                                 x
                                                EI ( x ) , m ( x )

                 x                         dx

                                            L

          Figure 17. Beam Subject to Dynamic Transverse Load, Analytical Example




                                              53 / 79
                                                     p( x,t ) dx


                                  V ( x,t )


                    M ( x,t )                                                                               ∂M ( x,t )
                                                                                              M ( x,t ) +              dx
                                                                                                              ∂x
                                                                                        ∂V ( x,t )
                                                                          V ( x,t ) +              dx
                                                                                          ∂x
                                                              f I ( x, t ) dx



                            x                            dx


                                Figure 18. Section of Beam in Figure 17

        The transverse displacement can be derived from the equation of motion as shown in

general by Clough and Penzien (1993, pp. 365-424) and Paz and Leigh (2004, pp. 527-552). Paz

uses Bernoulli-Euler and simple bending theories to setup the differential equation of motion for

a uniform beam. Shear force deformations and rotary inertia are neglected, as well as axial

effects. The frequency equation is

                                              cos an L cosh an L −1 = 0 .                                                   [6-30]

This relates to the natural frequency with the following

                                                                       EI
                                                ωn = ( an L )
                                                                   2
                                                                           .                                                [6-31]
                                                                       mL4

The corresponding normal mode shape is

                         Φ n ( x ) = cosh an x − cos an x − σ n ( sinh an x − sin an x )                                    [6-32]

where

                                                     cos an L − cosh an L
                                              σn =                        .                                                 [6-33]
                                                     sin an L − sinh an L



                                                          54 / 79
These mode shapes are assumed to be orthogonal to each other (Clough & Penzien, 1993, pp.

389-391). The complementary solution is

                                                 ∞
                                   yc ( x, t ) = ∑Φn ( x ) [ An cosωnt + Bn sin ωnt ]                             [6-34]
                                                n=1


where the An and Bn are determined by the initial conditions. Let the forcing function be

sinusoidal:

                                                     p ( x, t ) = p0 sin ω t .                                    [6-35]

The particular solution of the beam is
                                                  ∞
                                                                          In
                                    y p ( x, t ) = ∑ Φ n ( x )                     p0 sin ω t                     [6-36]
                                                  n =1                m (ω − ω 2 )
                                                                            2
                                                                            n



where

                                                                     Φ n ( x )dx
                                                                 L

                                                      In =
                                                             ∫ 0
                                                                                    .                             [6-37]
                                                                     Φ n ( x ) dx
                                                               L
                                                             ∫
                                                                       2
                                                               0



The general solution is the superposition of the complementary and particular solutions:
                          ∞                                                     ∞
                                                                                            In
              y ( x , t ) = ∑ Φ n ( x ) [ An cos ω n t + Bn sin ω n t ] + ∑ Φ n ( x )                p0 sin ω t   [6-38]
                          n =1                                                  n =1    m (ω − ω 2 )
                                                                                                2
                                                                                                n



                                 ∞          ⎡                                In                   ⎤
                   y ( x, t ) = ∑ Φ n ( x ) ⎢ An cos ωnt + Bn sin ωnt +                p0 sin ω t ⎥ .             [6-39]
                                n =1        ⎢
                                            ⎣                           m ( ωn − ω 2 )
                                                                             2
                                                                                                  ⎥
                                                                                                  ⎦

The initial conditions are

                                                         y ( x,0 ) = ρ ( x )                                      [6-40]

                                                         y ( x,0 ) = ψ ( x ) .                                    [6-41]

For 0 ≤ x ≤ L




                                                              55 / 79
                                                          ∞
                                               y ( x,0) = ∑ AnΦn ( x ) = ρ ( x )                           [6-42]
                                                          n=1


                                    ∞          ⎡             In             ⎤
                       y ( x,0 ) = ∑ Φ n ( x ) ⎢ Bnωn +                ω p0 ⎥ = ψ ( x ) .                  [6-43]
                                   n =1        ⎢
                                               ⎣        m ( ωn − ω 2 )
                                                             2
                                                                            ⎥
                                                                            ⎦

Assuming that the beam is initially undeformed ( ρ ( x ) = 0 ) and at rest (ψ ( x ) = 0 ), then

                                                              An = 0                                       [6-44]

                                                               In      ω p0
                                                Bn = −                      .                              [6-45]
                                                         ω n (ω n − ω ) m
                                                                2    2




Inserting Eq. [6-44] and [6-45] into Eq. [6-39] gives

                                   ∞
                                                           p0     ⎡ ω                 ⎤
                      y ( x, t ) = ∑ Φ n ( x ) I n             2 ⎢
                                                                   − sin ωnt + sin ωt ⎥ .                  [6-46]
                                  n=1                  m (ω − ω ) ⎣ ωn
                                                                2
                                                                n                     ⎦

The first nine modal frequencies and shapes were used to estimate the analytical model. The

numerical values for the modal information are taken and calculated from Blevins (1979, p. 108)

and are shown in Table 1 below.

                                Table 1. Modal Data for Beam Problem

 n         an L         ωn (in   rad
                                       s   )             σn                     In          Φ n ( x = 0.5 L )
  1    4.73004074           65.2                  0.982502215          0.83086152256645         1.5881
  2    7.85320462           179.5                 1.000777312                  0                0.0017
  3    10.9956079           291.3                 0.999966450          0.36376941224454        -1.2382
  4    14.1371655           581.6                 1.000001450                  0               -0.0041
  5    17.2787597           869.8                 0.999999937          0.23149808369867         1.4146
  6       6.5π             1213.4                      1.0                     0                0.0059
  7       7.5π             1616.9                      1.0             0.16976527261064        -1.4142
  8       8.5π             1972.8                      1.0                     0               -0.4693
  9       9.5π             2595.5                      1.0             0.13402521523525         1.4142




                                                           56 / 79
6.2.2. Numerical Solution

       A two-dimensional space-time finite element model was simulated for the transverse

beam problem as shown in Figure 19.


y

                        A                   p(x,t)=(75sin(0.011π t)kN )H(t)
                                                            2
                                                                    m

                                                                                                 x
                        A
                                        L=10m

          Figure 19. Beam Subject to Dynamic Transverse Load, Numerical Example




                                     H=0.2m



                                 T=0.005m
                      Figure 20. Cross Section A-A of Beam in Figure 19

 H ( t ) is the Heaviside step function in Figure 19. The material properties of the bar are:

volumetric mass density is ρ = 7860 kg m3 ; Young’s modulus of elasticity is E = 200GPa ; and

Poisson’s ratio is ν = 0.3 . The 2D finite element model is treated as a plane stress problem. The

mesh used to produce the results consisted of 65 elements to discretize the beam along its length,

3 elements along the height of the beam, and 460 space-time slabs along the time axis with a

time step of Δt = 0.5 ×10−3 s . Quartic temporal shape functions and quadratic spatial shape

functions (Q9 elements) were used. The isoparametric form of the 2D spatial shape functions is


                                             57 / 79
               ⎡ 1 ( 0 − ξ )(1 − ξ )( 0 − η )(1 − η )
                                                                                                          T
                     4                                                           0                    ⎤
               ⎢                                                                                      ⎥
                                     0                          4 ( 0 − ξ )(1 − ξ )( 0 − η )(1 − η )
                                                                1
               ⎢                                                                                      ⎥
               ⎢ 1 ( −1 − ξ )( 0 − ξ )( 0 − η )(1 − η )
                   4                                                             0                    ⎥
               ⎢                                                                                      ⎥
                                     0                        4 ( −1 − ξ )( 0 − ξ )( 0 − η )(1 − η ) ⎥
                                                              1
               ⎢
               ⎢ 1 ( −1 − ξ )( 0 − ξ )( −1 − η )( 0 − η )                        0                    ⎥
               ⎢4                                                                                     ⎥
                                     0                      4 ( −1 − ξ )( 0 − ξ )( −1 − η )( 0 − η ) ⎥
                                                            1
               ⎢
               ⎢ 1 ( 0 − ξ )(1 − ξ )( −1 − η )( 0 − η )                          0                    ⎥
               ⎢ 4                                                                                    ⎥
               ⎢                     0                        1
                                                              4  ( 0 − ξ )(1 − ξ )( −1 − η )( 0 − η ) ⎥
               ⎢ − 1 −1 − ξ 1 − ξ 0 − η 1 − η                                                         ⎥
          Ns = ⎢ 2(            )( )(            )( )                             0                    ⎥
               ⎢                     0                      − 1 ( −1 − ξ )(1 − ξ )( 0 − η )(1 − η ) ⎥
                                                                2
               ⎢ 1                                                                                    ⎥
               ⎢ − 2 ( −1 − ξ )( 0 − ξ )( −1 − η )(1 − η )                       0                    ⎥
               ⎢                     0                     − 1 ( −1 − ξ )( 0 − ξ )( −1 − η ) (1 − η ) ⎥
               ⎢ 1                                            2
                                                                                                      ⎥
               ⎢ − 2 ( −1 − ξ )(1 − ξ )( −1 − η )( 0 − η )                       0                    ⎥
               ⎢                     0                     − 2 ( −1 − ξ )(1 − ξ )( −1 − η )( 0 − η ) ⎥
                                                              1
               ⎢                                                                                      ⎥
               ⎢ − 1 ( 0 − ξ )(1 − ξ )( −1 − η )(1 − η )
                     2                                                           0                    ⎥
               ⎢                     0                      − 2 ( 0 − ξ )(1 − ξ )( −1 − η )(1 − η ) ⎥
                                                                1
               ⎢                                                                                      ⎥       [6-47]
               ⎢ ( −1 − ξ )(1 − ξ )( −1 − η )(1 − η )                            0                    ⎥
               ⎢
               ⎣                     0                        ( −1 − ξ )(1 − ξ )( −1 − η )(1 − η ) ⎥  ⎦

with nodes at


                              ⎡ξ ⎤ ⎡ −1 1 1 −1 0 1 0 −1 0⎤
                              ⎢η ⎥ = ⎢ −1 −1 1 1 −1 0 1 0 0⎥ .                                                [6-48]
                              ⎣ ⎦ ⎣                        ⎦

The spatial finite element mesh is shown in Figure 21 and Figure 22. Circles indicate nodes

where tractions are applied and squares indicate nodes where boundary conditions are applied.

The forcing function has a period of 0.011s (making its frequency about 90 Hz ). The behavior of

the displacement is strongly dominated by the first mode of vibration, which has a period near

0.1s (and a frequency near 10Hz ). The simulation was run for a total time period of 0.230s in

order to capture at least two periods of the first mode.




                                                     58 / 79
                                                             4


                                                             3


                                                             2



                                  Initial Y-position, in m
                                                             1


                                                             0


                                                             -1


                                                             -2


                                                             -3


                                                                 0       1         2         3      4        5        6      7         8         9         10
                                                                                                 Initial X-position, in m



                                                                             Figure 21. 2D Beam Spatial Mesh



                                                      0.4


                                                      0.3
                 Initial Y-position, in m




                                                      0.2


                                                      0.1


                                                             0


                                                -0.1


                                                -0.2

                                                                     0       0.1       0.2        0.3       0.4        0.5       0.6       0.7       0.8
                                                                                                 Initial X-position, in m


                                                                     Figure 22. Zoomed In 2D Beam Spatial Mesh

6.2.3. Convergence Results

      The analytical solution provides displacements for the transverse beam. The midpoint of

the beam exhibits the greatest displacements, so the Y-displacement of the midpoint of the


                                                                                                     59 / 79
neutral axis (the very center of the beam) is used to study the convergence rate of the temporal

mesh. The L2 norm of the residual displacement at this point is used taking a summation across

all the temporal nodes. The plot of this convergence is shown below in Figure 23. The slope of

this line is 4.010 .

                            -3
                           10




                            -4
                           10
                 ||e||L2




                            -5
                           10




                            -6
                           10




                            -7
                           10         1.2       1.3       1.4       1.5       1.6        1.7
                                    10        10        10        10        10          10
                                   Characteristic Temporal Element Length (cΔt), in m


                                 Figure 23. L 2 Residual Displacement Norm

The energy norm used by Hulbert and Hughes ||| u h ||| was also calculated and shown in Figure

24. The slope of this curve is 3.053.




                                                      60 / 79
                               3
                              10




                               2
                              10
                  |||uh|||2




                               1
                              10




                               0
                              10 0                            1                            2
                                10                          10                            10
                                     Characteristic Temporal Element Length (cΔt), in m



                                           Figure 24. Total Energy Norm

6.2.4. Numerical Results

       The Y-directional transverse displacement of the finite element model is shown in the

mesh of Figure 25. This displacement history shows the dominance of the first mode of

vibrations at about 10Hz (with a period of about 0.1s ). The period of the externally applied

forcing function is shown in the smaller perturbations.




                                                         61 / 79
             Figure 25. Mesh of Transverse Displacement of Straight Beam, in m

The analytical solution provides displacements for the transverse beam. The midpoint of the

beam exhibits the greatest displacements, so the Y-displacement at this midpoint is shown in

Figure 26.




                                          62 / 79
                                                          0.03




                 Y-displacement at middle of beam, in m
                                                          0.02


                                                          0.01


                                                               0                                           Analytical Model
                                                                                                           Space-Time Model
                                                          -0.01


                                                          -0.02


                                                          -0.03


                                                               0          0.05              0.1              0.15               0.2
                                                                                              Time, in s




             Figure 26. Transverse Displacement of Midpoint of Straight Beam, in m

6.3.    Transversely Loaded Thick Curved Beam

        One of the major goals of finite element analysis is to model the movement of a structure

whose analytical solution is unknown or highly complex. An arch (curved beam) provides more

complexity than the straight beam from the last section because of the coupling between the axial

and transverse or radial loadings. The arch that was simulated using the space-time finite element

code is shown in Figure 27.

                                                                                   L = 18 in [ 457 mm ]
                                                                                                                        ⎛          ⎛ 2π t ⎞ ⎞ lbf    ⎛      ⎛ 2π t ⎞ ⎞ kN
                                                                                                           p ( x, t ) = ⎜ 0.225sin ⎜         ⎟⎟     = 75sin ⎜         ⎟⎟
                                                                                                                        ⎝          ⎝ 0.0001s ⎠ ⎠  in ⎜
                                                                                                                                                     ⎝      ⎝ 0.0001s ⎠ ⎠ m

 y
                                                          A

         x                                                         H = 0.5 in [ 12.7 mm ]
                                                           A                                  R = 81.25 in [ 2.064 m ]


     Figure 27. Thick Curved Beam Subject to Dynamic Transverse Load, Numerical Example




                                                                                            63 / 79
                             T = 0.02 in [ 0.5 mm ]



                                H = 0.5 in [ 12.7 mm ]




               Figure 28. Cross Section A-A of Thick Curved Beam in Figure 27

       The arch is assumed to be made of isotropic steel with the following properties:

volumetric mass density ρ = 7860 kg m3 , Young’s modulus of elasticity E = 200GPa , and

Poisson’s ratio ν = 0.3 . The finite element model is discretized using quadratic spatial shape

function (Q9 elements) and quartic temporal shape functions. The relative thickness of this arch

makes the implementation quick as compared to a thin arched beam. The regularity of the spatial

elements in the mesh can be maintained without becoming extremely small. The initial mesh is

shown in Figure 29 below. It consists of 75 elements lengthwise and 3 elements height wise. The

temporal domain is discretized with a time step of Δt = 10−5 s and the simulation is run for a total

of 0.00455 s (using 455 space-time slabs).




                                              64 / 79
                                                      2.45


                                                       2.4


                                                      2.35
                           Initial Y-position, in m

                                                       2.3


                                                      2.25


                                                       2.2


                                                      2.15


                                                             0      0.05   0.1     0.15         0.2     0.25       0.3    0.35      0.4     0.45
                                                                                          Initial X-position, in m



                                                                 Figure 29. Spatial Mesh of Thick Curved Beam




                                                      2.29
                Initial Y-position, in m




                                              2.285



                                                      2.28



                                              2.275



                                                      2.27


                                                                     0     0.005          0.01        0.015        0.02     0.025         0.03
                                                                                          Initial X-position, in m


                                       Figure 30. Zoomed In Spatial Mesh of Thick Curved Beam

The X- and Y-directional displacements are shown in the meshes of the following figures.




                                                                                              65 / 79
Figure 33. The frequency of the applied traction is seen as the individual period of the work and

the frequency of the first mode of vibration dominates the amplitude.


                                                -3
                                             x 10


                                      1.5


                                        1


                                      0.5
                 Work Applied, in J




                                        0


                                      -0.5

                                       -1


                                      -1.5


                                       -2
                                                     0.5   1   1.5   2      2.5   3   3.5   4      4.5
                                                                     Time, in s                    -3
                                                                                                x 10



                                         Figure 33. Work Done to Thick Curved Beam, in J

       The curved beam problem example shows the ability of the space-time formulations to

hand excitations in the high frequency range. The excitation in this example has a period of

0.0001s (or a frequency of 10000 Hz )


6.4.   Transversely Loaded Slender Curved Beam

       Thin beams and arches are of practical use in application and can simplify the analytical

form of their displacements. A slender arch that was simulated using the space-time finite

element code is shown in Figure 34.




                                                                     67 / 79
                                                     ⎛          ⎛ 2π t ⎞ ⎞ lbf    ⎛      ⎛ 2π t ⎞ ⎞ kN
                                        p ( x, t ) = ⎜ 0.225sin ⎜         ⎟⎟     = 75sin ⎜         ⎟⎟
                                                     ⎝          ⎝ 0.0001s ⎠ ⎠  in ⎜
                                                                                  ⎝      ⎝ 0.0001s ⎠ ⎠ m
                                         L = 12 in [ 0.305 m ]
y
                    A

       x
                    A
                    H = 0.024 in [ 0.00061 m ]                   R = 71.25 in [ 1.810 m ]

    Figure 34. Slender Curved Beam Subject to Dynamic Transverse Load, Numerical Example


                                                                    T = 0.5 in [ 12.7 mm ]



                 H = 0.024 in [ 0.610 mm ]


                 Figure 35. Cross Section A-A of Slender Curved Beam in Figure 34

           The arch is similar to the arch used in the previous section in that it is also assumed to be

made of isotropic steel with the following properties: volumetric mass density ρ = 7860 kg m3 ,

Young’s modulus of elasticity E = 200GPa , and Poisson’s ratio ν = 0.3 . The finite element

model is discretized using quadratic spatial shape function (Q9 elements) and quartic temporal

shape functions. The relative thinness of this arch makes the implementation take an extended

period of time. In order to maintain some regularity of the spatial elements, 500 elements were

used to mesh the length of the arch and 2 along its height. The temporal domain is discretized

using a time step of Δt = 10−5 s and the simulation was run for a time length of 3.03 ×10−3 s

(using 303 space-time slabs). The initial spatial mesh is shown in Figure 36 with a zoomed in

view in Figure 37 below.




                                                  68 / 79
                                                                      2.1




                                                                     2.05



                                          Initial Y-position, in m
                                                                       2




                                                                     1.95




                                                                      1.9



                                                                        0         0.05     0.1             0.15          0.2         0.25        0.3
                                                                                                 Initial X-position, in m



                                                                            Figure 36. Spatial Mesh of Slender Curved Beam


                                                1.9795



                                                          1.979
               Initial Y-position, in m




                                                1.9785



                                                          1.978



                                                1.9775



                                                          1.977


                                                                              0     0.5    1             1.5          2        2.5          3
                                                                                                 Initial X-position, in m                          -3
                                                                                                                                                x 10


                                  Figure 37. Zoomed in Spatial Mesh of Slender Curved Beam

6.4.1. Convergence Results

       For convergence study, a fine mesh of 500 elements lengthwise and 2 elements

heightwise is used as an approximation of the analytical solution. The L2 norm shown in Figure


                                                                                                     69 / 79
38 has a slope of 2.815. The convergence study shows the robustness of the space-time method

in dealing with the pressure load which is sinusoidal function of time.

                              0
                             10
                  ||e||L 2




                              -1
                             10




                              -2
                             10     -2.9            -2.8                -2.7                -2.6
                                   10            10                    10                  10
                                                Characteristic Spatial Element Length, in m


                                        Figure 38. L 2 Residual Displacement Norm

6.4.2.   Numerical Results

         The displacement histories of the thin curved beam are shown in Figure 39 and Figure 40.

These histories are taken from the nodes at the midline of the curved beam. These two

directional displacements are coupled because of the geometry of the structure and the modes

present in the normal and tangential directions.




                                                              70 / 79
               Figure 31. Mesh of X Displacement of Thick Curved Beam, in m




               Figure 32. Mesh of Y Displacement of Thick Curved Beam, in m

       The behavior of the arched beam is similar to the straight beam in that it seems that the

displacement is dominated by the first mode of vibration at a frequency of about 500 Hz (and a

period near 0.002 s ). This behavior is clearly seen in plot of the work done to the system in

                                            66 / 79
              Figure 39. Mesh of X Displacement of Slender Curved Beam, in m




              Figure 40. Mesh of Y Displacement of Slender Curved Beam, in m

The stress histories are shown in Figure 41 through Figure 43. These meshes show the averages

between each row of elements in order to represent the stresses along the neutral axis.


                                              71 / 79
Figure 41. Stress in the XX-Direction, in Pa




Figure 42. Stress in the YY-Direction, in Pa




                  72 / 79
Figure 43. Stress in the XY-Direction, in Pa




                  73 / 79
Chapter 7.           Conclusions and Future Work

       In summary, this thesis presents the use of space-time finite element methods for

structural response prediction under dynamic loads. The method developed is then applied to

problems in both one and two spatial dimensions and convergence studies are performed to

ensure the consistency of the method. Based on the work performed, the following conclusions

can be drawn:


       1) The space-time finite element method allows for establishment of approximations in

          both space and time so that it is not restricted to a particular time interpolation as

          implied in a typical semi-discrete scheme. As such, the method is capable of refining

          the numerical resolutions in both the spatial and temporal domains by using different

          types of approximations in different domains. While the semi-discrete finite element

          scheme can also achieve this through the use of refined spatial mesh or subcycling

          algorithm, the time discretization incorporated in the space-time finite element

          scheme provides a direct approach to converge to the solutions that have strong time-

          dependent characteristics through systematic adaptivity. In contrast, it is difficult to

          obtain error estimates on time for a semi-discrete scheme since no mesh is used in the

          time domain and the switch between different time integration schemes can be

          tedious.

       2) The semi-discrete explicit analysis is limited in the way it resolves responses in the

          high frequency range and responses having long time duration. The ability to capture

          sharp gradients and discontinuities is also limited. It has been shown that most of

          these issues can be dealt with using the space-time finite element method by

          establishing a high-order accurate approximation in time. In addition, the time step

                                            74 / 79
           that limits explicit semi-discrete methods is lifted by this formulation since the

           algorithm is stable.

       3) The multiplicative form of the formulation using the Lagrangian type of shape

           function makes the implementation straightforward. It should be noted, however, the

           use of multiplicative form is not necessary and other types of finite element shape

           functions can also be incorporated in a straightforward manner.

       4) Based on the systematic convergence study for a number of problems involving

           dynamic structural response prediction, it is concluded that the space-time finite

           element code developed in this thesis is robust and accurate.


       For future research efforts, it is proposed that the research presented in this thesis can be

extended in various directions. Specifically,


       1) The space-time formulation and code can be extended to introduce multiscale

           approximation through the enrichment methodology based on the extended finite

           element method. This methodology will enhance the numerical resolution through

           incorporating analytical basis.

       2) The space-time framework can be extended to perform more reliable nonlinear

           analysis of complex thermo-mechanical problems. This is important to assess the

           effect of the coupling effect.

       3) The work can be extended to 3D structures, in particular the application of shell and

           plate structures will be of great interest to the research and industry community.

       4) The incorporation of damping into the space-time finite element formulation can lead

           to more accurate response prediction for problems where damping cannot be

           neglected.

                                                75 / 79
Bibliography

Argyris, J. H., & Scharpf, D. W. (1969). Finite elements in space and time. Nuclear Engineering
and Design , 10 (4), 456-464.

Belytschko, T., Liu, W. K., & Moran, B. (2000). Nonlinear Finite Elements for Continua and
Structures. Chichester, West Sussex, England: John Wiley & Sons.

Blevins, R. D. (1979). Formulas for Natural Frequency and Mode Shape. New York, New York,
United States of America: Van Nostrand Reinhold.

Cavin, P., Gravouil, A., Lubrecht, A. A., & Combescure, A. (2005). Automatic energy
conserving space-time refinement for linear dynamic structural problems. International Journal
for Numerical Methods in Engineering , 64 (3), 304-321.

Chen, Z., Steeb, H., & Diebels, S. (2008). A new hybrid velocity integration method applied to
elastic wave propagation. International Journal for Numerical Methods in Engineering , 74 (1),
56-79.

Chen, Z., Steeb, H., & Diebels, S. (2006). A time-discontinuous Galerkin method for dynamical
analysis of porous media. International Journal for Numerical and Analytial Methods in
Geomechanics , 30 (11), 1113-1134.

Chessa, J., & Belytschko, T. (2006). A local space-time discontinuous finite element method.
Computer Methods in Applied Mechanics and Engineering , 195 (13-16), 1325-1343.

Chessa, J., & Belytschko, T. (2004). Arbitrary discontinuities in space-time finite element by
level sets and X-FEM. International Journal for Numerical Methods in Engineering , 61 (15),
2595-2614.

Chien, C.-C., Chen, Y.-H., & Chuang, C.-C. (2003). Dual reciprocity BEM analysis of 2D
transient elastodynamic problems by time-discontinuous Galerkin FEM. Engineering Analysis
with Boundary Elements , 27 (6), 611-624.

Clough, R. W., & Penzien, J. (1993). Dynamics of Structures (2nd Edition ed.). (B. J. Clark, Ed.)
New Yok, New York, United States of America: McGraw-Hll.

Cook, R. D., Malkus, D. S., Plesha, M. E., & Witt, R. J. (2001). Concepts and Applications of
Finite Element Analysis (4th Edition ed.). John Wiley & Sons.

Feng, Y. T., & Peri , D. (2000). A time-adaptive space-time finite element method for
incompressible Lagrangian flows with free surfaces: computational issues. Computer Methods in
Applied Mechanics and Engineering , 190 (5-7), 499-518.



                                             76 / 79
Fried, I. (1969). Finite element analysis of time-dependent phenomena. AIAA Journal , 7 (6),
1170-1173.

Froncioni, A. M., Labbe, P., Garon, A., & Camarero, R. (1997). Interpolation-free space-time
remeshing for the Burgers equation. Communications in Numerical Methods in Engineering , 13
(11), 875-884.

Gudla, P. K., & Ganguli, R. (2006). Discontinuous Galerkin finite element in time for solving
periodic differential equations. Computer Methods in Applied Mechanics and Engineering , 196
(1-3), 682-696.

Huang, H., & Costanzo, F. (2004). On the use of space-time finite elements in the solution of
elasto-dynamic fracture problems. International Journal of Fracture , 127 (2), 119-146.

Huang, H., & Costanzo, F. (2002). On the use of space-time finite elements in the solution of
elasto-dynamic problems with strain discontinuities. Computer Methods in Applied Mechanics
and Engineering , 191 (46), 5315-5343.

Huang, H., & Costanzo, F. (2005). Proof of unconditional stability for a single-field
discontinuous Galerkin finite element formulation for linear elasto-dynamics. Computer Methods
in Applied Mechanics and Engineering , 194 (18-20), 2059-2076.

Hughes, T. J., & Hulbert, G. M. (1988). Space-time finite element methods for elastodynamics:
formulations and error estimates. Computer Methods in Applied Mechanics and Engineering , 66
(3), 339-363.

Hughes, T. J., & Stewart, J. R. (1996). A space-time formulation for multiscale phenomena.
Journal of Computational and Applied Mathematics , 74 (1-2), 217-229.

Hughes, T. J., Franca, L. P., & Hulbert, G. M. (1989). A new finite element formulation for
computational fluid dynamics: VIII. the Galerkin/least-squares method for advective-diffusive
equations. Computer Methods in Applied Mechanics and Engineering , 73 (2), 173-189.

Hulbert, G. M. (1992). Time finite element methods for structural dynamics. International
Journal for Numerical Methods in Engineering , 33 (2), 307-331.

Hulbert, G. M., & Hughes, T. J. (1990). Space-time finite element methods for second-order
hyperbolic equations. Computer Methods in Applied Mechanics and Engineering , 84 (3), 327-
348.

Idesman, A. V. (2007a). A new high-order accurate continuous Galerkin method for linear
elastodynamics problems. Computational Mechanics , 40 (2), 261-279.




                                           77 / 79
Idesman, A. V. (2007b). Solution of linear elastodynamics problems with space-time finite
elements on structured and unstructured meshes. Computer Methods in Applied Mechanics and
Engineering , 196 (9-12), 1787-1815.

Idesman, A. V., Schmidt, M., & Sierakowski, R. L. (2008). A new explicit predictor-
multicorrector high-order accurate method for linear elastodynamics. Journal of Sound and
Vibration , 310 (1-2), 217-229.

Johnson, C. (1993). Discontinuous Galerkin finite element methods for second order hyperbolic
problems. Computer Methods in Applied Mechanics and Engineering , 107 (1-2), 117-129.

Karao lan, L., & Noor, A. K. (1997). Space-time finite element methods for the sensitivity
analysis of contact/impact response of axisymmetric composite structures. Computer Methods in
Applied Mechanics and Engineering , 144 (3-4), 371-389.

Kinsler, L. E., Frey, A. R., Coppens, A. B., & Sanders, J. V. (2000). Fundamentals of Acoustics
(4th Edition ed.). New York, New York: John Wiley & Sons.

Kurdi, M. H., & Beran, P. S. (2008). Spectral element method in time for rapidly actuated
systems. Journal of Computational Physics , 227 (3), 1809-1835.

Li, X. D., & Wiberg, N. E. (1998). Implementation and adaptivity of a space-time finite element
method for structural dynamics. Computer Methods in Applied Mechanics and Engineering , 156
(1), 211-229.

Li, X., Yao, D., & Lewis, R. W. (2003). A discontinuous Galerkin finite element method for
dynamic and wave propagation problems in non-linear solids and saturated porous media.
International Journal for Numerical Methods in Engineering , 57 (12), 1775-1800.

Masud, A., & Hughes, T. J. (1997). A space-time Galerkin/least-squares finite element
formulation of the Navier-Stokes equations for moving domain problems. Computer Methods in
Applied Mechanics and Engineering , 146 (1-2), 91-126.

N'dri, D., Garon, A., & Fortin, A. (2001). A new stable space-time formulation for two-
dimensional and three-dimensional incompressible viscous flow. International Journal for
Numerical Methods in Fluids , 37 (8), 865-884.

Nguyen, H., & Reynen, J. (1984). A space-time least square finite element scheme for advection-
diffusion equations. Computer Methods in Applied Mechanics and Engineering , 42 (3), 331-
342.

Oden, J. T. (1969). A general theory of finite elements II. applications. International Journal for
Numerical Methods in Engineering , 1, 247-259.



                                              78 / 79
Oñate, E., & Manzán, M. (1999). A general procedure for deriving stabilized space-time finite
element methods for advective-diffusive problems. International Journal for Numerical Methods
in Fluids , 31 (1), 203-221.

Pattillo, P. D., & Tortorelli, D. A. (2004). A contact algorithm for the Signorini problem using
space-time finite elements. International Journal for Numerical Methods in Engineering , 60 (7),
1197-1213.

Paz, M., & Leigh, W. (2004). Structural Dynamics: Theory and Computation (5th Edition ed.).
Boston, Massachusetts, United States of America: Kluwer Academic Publishers.

Shakib, F., & Hughes, T. J. (1991). A new finite element formulation for computational fluid
dynamics: IX. Fourier analysis of space-time Galerkin/least-squares algorithms. Computer
Methods in Applied Mechanics and Engineering , 87 (1), 35-58.

Shaw, S., & Whiteman, J. R. (2000). Adaptive space-time finite element solution for Volterra
equations arising in viscoelasticity problems. Journal of Computational and Applied
Mathematics , 125 (1-2), 337-345.

Surana, K. S., Allu, S., Tenpas, P. W., & Reddy, J. N. (2007). k-version of finite element method
in gas dynamics: higher-order global differentiability numerical solutions. International Journal
for Numerical Methods in Engineering , 69 (6), 1109-1157.

Thompson, L. L., & Pinsky, P. M. (1996). A space-time finite element method for the exterior
acoustics problem. The Journal of the Acoustical Society of America , 99 (6), 3297-3311.

Varoḡlu, E., & Finn, W. D. (1978). A finite element method for the diffusion-convection
equation with constant coefficients. Advances in Water Resources , 1 (6), 337-343.

Varoḡlu, E., & Finn, W. D. (1980). Finite elements incorporating characteristics for one-
dimensional diffusion-convection equation. Journal of Computational Physics , 34 (3), 371-389.

Wiberg, N.-E., & Li, X. (1999). Adaptive finite element procedures for linear and non-linear
dynamics. International Journal for Numerical Methods in Engineering , 46, 1781-1802.




                                             79 / 79

								
To top