VIEWS: 2 PAGES: 88 POSTED ON: 1/20/2014
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