* Your assessment is very important for improving the workof artificial intelligence, which forms the content of this project

# Download STABILISED FINITE ELEMENT SOLUTION OF

Survey

Document related concepts

Analytical mechanics wikipedia , lookup

N-body problem wikipedia , lookup

Theoretical and experimental justification for the Schrödinger equation wikipedia , lookup

Viscoelasticity wikipedia , lookup

Classical central-force problem wikipedia , lookup

Finite element method wikipedia , lookup

Relativistic quantum mechanics wikipedia , lookup

Hooke's law wikipedia , lookup

Derivations of the Lorentz transformations wikipedia , lookup

Deformation (mechanics) wikipedia , lookup

Routhian mechanics wikipedia , lookup

Wave packet wikipedia , lookup

Equations of motion wikipedia , lookup

Transcript

STABILISED FINITE ELEMENT SOLUTION OF LAGRANGIAN EXPLICIT DYNAMIC SOLID MECHANICS by Dhrubajyoti Mukherjee Thesis submitted to Swansea University in fulfilment of the requirements for the degree of Master of Science in Computational Mechanics May 2012 Abstract Since the advent of computational mechanics, the numerical modelling of transient phenomena has been a major field of interest in industry, including applications such as crash simulation, impact, forging and many others. Traditionally, a Lagrangian formulation is employed for the numerical simulation of these problems and low order spatial interpolation is preferred for computational workload convenience. For fast dynamics applications, the use of explicit time integrators is regarded as efficient in the majority of cases. The well-known second order solid dynamics formulation, where the primary variable is the displacement, is typically discretised in space by using the finite element method and discretised in the time domain by means of a Newmark time integrator. However, the resulting space-time discretised formulation presents a series of shortcomings. From the time discretisation point of view, the Newmark method has a tendency to introduce high frequency noise in the solution field, especially in the vicinity of sharp spatial gradients. From the space discretisation point of view, the use of isoparametric linear finite elements leads to second order convergence in displacements, but only first order convergence for stresses and strains. It is also known that these elements exhibit locking phenomena in incompressible or nearly incompressible scenarios. This work proposes a novel approach to obtain stabilised finite element solution of a new system of first order hyperbolic equations, which aims to alleviate locking by introducing the deformation gradient tensor as a conservation variable. A series of numerical experiments are performed in order to determine the feasibility of the one step Taylor-Galerkin and the Streamline Upwind Petrov-Galerkin finite element formulations. A comparative study is also conducted with finite volume method. Declarations and Statements DECLARATION This work has not previously been accepted in substance for any degree and is not being concurrently submitted in candidature for any degree. Signed ................................................ Date........................................... STATEMENT 1 This dissertation is the result of my own independent work/investigation, except where otherwise stated. Other sources are acknowledged by footnotes giving explicit references. A bibliography is appended. Signed ................................................ Date .......................................... STATEMENT 2 I hereby give my consent for my dissertation, if relevant and accepted, to be available for photocopying and for inter-library loan, and for the title and summary to be made available to outside organisations. Signed ................................................ Date ........................................... Acknowledgements This work would not have come into existence without the help and support of some of the best minds I have ever come across. I am grateful to Dr. Antonio Gil for his guidance and dedication to this work. Every conversation with him has been highly inspiring and intellectually stimulating. I would also like to extend my gratitude to my supervisor Prof. Javier Bonet. He has helped me to keep an open mind and encouraged me to explore the directions of research, which are otherwise considered as very unconventional. I am indebted to Dr. Chun Hean Lee and Mr. Miquel Aguirre for the lively discussion sessions. Those meetings have certainly helped me to gain an overall perspective of the research area and its key aspects. I would like to thank Prof. Djordje Peric, Dr. Rubén Sevilla, Dr. Aurelio Arranz Carreño and Mr. Rogelio Ortigosa Martínez for taking an avid interest in the work and giving invaluable suggestions for improvement. I am grateful to Prof. Carlos Agelet De Saracibar and Dr. Sergio Zlotnik at UPC, Barcelona for providing me with required expertise. Special thanks to Mr. Rodolfo Fleury and Mr. Renjith Krishnan for helping me gather relevant literatures. I wish to thank my friends, especially Caner, Hannah, Astor, Anais, Georgina, Fanny, Edurne, Natalie, Selina, Tiphaine, Anna, Elise, Alba, Mike, Nikolay, Sylvain, Manon and Regina for supporting me during this research period. I am also grateful to European Commission for granting Erasmus Mundus scholarship. Finally, I would like to express my gratitude to my parents and brother for their unconditional support and compassion. Contents LIST OF TABLES vi LIST OF FIGURES x ABBREVIATION xi NOTATION xii 1 2 INTRODUCTION 1 1.1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 EARLY DEVELOPMENTS IN FEM . . . . . . . . . . . . . . . . . . . . 3 1.3 MOTIVATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 OBJECTIVES AND SCOPE OF WORK . . . . . . . . . . . . . . . . . . 9 1.5 LAYOUT OF THESIS . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 LARGE STRAIN STRUCTURAL DYNAMICS 12 2.1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2 KINEMATICS PRELIMINARIES . . . . . . . . . . . . . . . . . . . . . . 13 2.2.1 The Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.2 Deformation Gradient Tensor . . . . . . . . . . . . . . . . . . . . 14 i ii Contents 2.2.3 Volume Change . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Velocity Measures . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.2 Rate of Deformation . . . . . . . . . . . . . . . . . . . . . . . . . 17 GOVERNING EQUATIONS . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4.1 Conservation of Linear Momentum . . . . . . . . . . . . . . . . . 18 2.4.2 Conservation of Deformation Gradient . . . . . . . . . . . . . . . 19 2.4.3 Conservation of Energy . . . . . . . . . . . . . . . . . . . . . . . 20 2.4.4 Conservative Law Formulation . . . . . . . . . . . . . . . . . . . . 21 2.4.5 Interface Flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 HYPERELASTICITY . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.5.1 Nearly Incompressible Neo-Hookean Material . . . . . . . . . . . 25 2.5.2 Linear Elasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.6 EIGENSTRUCTURE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.7 SHOCK AND CONTACT CONDITIONS . . . . . . . . . . . . . . . . . . 31 2.8 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 2.3 2.4 2.5 3 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 35 3.1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2 TAYLOR-GALERKIN SCHEME . . . . . . . . . . . . . . . . . . . . . . 36 3.3 DISCRETISATION IN SPACE: SUPG . . . . . . . . . . . . . . . . . . . . 38 3.4 DISCRETISATION IN TIME . . . . . . . . . . . . . . . . . . . . . . . . 41 3.4.1 Forward Euler Method . . . . . . . . . . . . . . . . . . . . . . . . 42 3.4.2 TVD Runge-Kutta Scheme . . . . . . . . . . . . . . . . . . . . . . 43 iii Contents 3.5 CONSISTENCY ANALYSIS . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.5.1 Consistency Analysis of The Taylor-Galerkin Scheme . . . . . . . 45 3.5.2 Consistency Analysis of The SUPG-FE Scheme . . . . . . . . . . . 46 VON NEUMANN STABILITY ANALYSIS . . . . . . . . . . . . . . . . . 48 3.6.1 Stability Analysis of The Taylor-Galerkin Scheme . . . . . . . . . 49 3.6.2 Stability Analysis of The SUPG-FE Scheme . . . . . . . . . . . . 51 SPECTRAL ANALYSIS . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.7.1 Spectral Analysis of The Taylor-Galerkin Scheme . . . . . . . . . . 54 3.7.2 Spectral Analysis of SUPG-FE Scheme . . . . . . . . . . . . . . . 55 NUMERICAL EXAMPLE: COSINE WAVE . . . . . . . . . . . . . . . . 57 3.8.1 Results & Discussions . . . . . . . . . . . . . . . . . . . . . . . . 59 3.8.2 Convergence Study . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.8.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 RIEMANN PROBLEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.10 NUMERICAL EXAMPLE: STEP WAVE . . . . . . . . . . . . . . . . . . 66 3.10.1 Results & Discussions . . . . . . . . . . . . . . . . . . . . . . . . 66 3.11 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.6 3.7 3.8 3.9 4 SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 70 4.1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.2 ANALYTICAL SOLUTION . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.2.1 Numerical Example: Propagation of Acoustic Wave . . . . . . . . 74 NUMERICAL SOLUTION . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.3 iv Contents 4.4 4.5 4.6 5 4.3.1 Taylor-Galerkin Scheme . . . . . . . . . . . . . . . . . . . . . . . 75 4.3.2 SUPG-RK Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.3.3 Numerical Example: Propagation of Acoustic Wave . . . . . . . . 81 CONSERVATIVE FRAMEWORK . . . . . . . . . . . . . . . . . . . . . . 82 4.4.1 Taylor-Galerkin Scheme in Conservative Form . . . . . . . . . . . 84 4.4.2 SUPG-RK Scheme in Conservative Form . . . . . . . . . . . . . . 86 SELECTION OF STABILISATION PARAMETER . . . . . . . . . . . . . 87 4.5.1 Numerical Experiment: Propagation of Acoustic Wave . . . . . . . 89 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 1D FAST STRUCTURAL DYNAMICS PROBLEM 92 5.1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.2 GOVERNING EQUATION . . . . . . . . . . . . . . . . . . . . . . . . . . 93 5.3 BAR SUBJECTED TO SINUSOIDAL TRACTION . . . . . . . . . . . . . 95 5.4 CONVERGENCE STUDY . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.5 BAR SUBJECTED TO CONSTANT TRACTION . . . . . . . . . . . . . . 103 5.6 5.5.1 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.5.2 Effect of Varying Stabilisation Parameter . . . . . . . . . . . . . . 106 5.5.3 Comparison of Different Schemes . . . . . . . . . . . . . . . . . . 109 SHOCK CAPTURING SCHEME . . . . . . . . . . . . . . . . . . . . . . 111 5.6.1 Y Zβ Shock Capturing . . . . . . . . . . . . . . . . . . . . . . . . 112 5.6.2 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . 113 5.7 NON-LINEAR ELASTICITY . . . . . . . . . . . . . . . . . . . . . . . . 115 5.8 LUMPED MASS FORMULATION . . . . . . . . . . . . . . . . . . . . . 117 v Contents 5.9 6 2D FAST STRUCTURAL DYNAMICS PROBLEM 123 6.1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.2 CURL-FREE INVOLUTION . . . . . . . . . . . . . . . . . . . . . . . . . 124 6.3 NUMERICAL EXAMPLES . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.4 7 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 6.3.1 Short Column . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.3.2 Tensile Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 CONCLUSIONS 129 7.1 SUMMARY AND CONCLUSION . . . . . . . . . . . . . . . . . . . . . 130 7.2 FUTURE RESEARCH DIRECTIONS . . . . . . . . . . . . . . . . . . . . 132 BIBLIOGRAPHY 144 List of Tables 5.1 Properties of axially-loaded bar . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Properties of short column . . . . . . . . . . . . . . . . . . . . . . . . . . 125 6.2 Properties of plate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 vi 97 List of Figures 1.1 Collapse of a building under earthquake [1] . . . . . . . . . . . . . . . . . 5 1.2 Collision of two cars [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 Welding of a metal [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Motion of a deformable body . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Surface discontinuity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.3 Contact motion generated shock waves . . . . . . . . . . . . . . . . . . . . 33 2.4 Different boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . 34 3.1 Polar representation of the amplification factor for the Taylor-Galerkin scheme 51 3.2 Polar representation of the amplification factor for the SUPG-FE scheme . . 3.3 Cartesian representation of the diffusion error as a function of phase angle, in degrees, for Taylor-Galerkin scheme . . . . . . . . . . . . . . . . . . . . 3.4 56 Cartesian representation of the diffusion error as a function of phase angle, in degrees, for SUPG-FE scheme . . . . . . . . . . . . . . . . . . . . . . . 3.6 55 Cartesian representation of the dispersion error as a function of phase angle, in degrees, for Taylor-Galerkin scheme . . . . . . . . . . . . . . . . . . . . 3.5 53 56 Cartesian representation of the dispersion error as a function of phase angle, in degrees, for SUPG-FE scheme . . . . . . . . . . . . . . . . . . . . . . . vii 57 viii List of Figures 3.7 Propagation of the wave using Bubnov-Galerkin (BG) scheme . . . . . . . 58 3.8 Algorithm for solving 1D advection equation . . . . . . . . . . . . . . . . 59 3.9 Propagation of a cosine wave using Taylor-Galerkin scheme . . . . . . . . 60 3.10 Propagation of a cosine wave using SUPG-FE scheme . . . . . . . . . . . . 62 3.11 Propagation of a cosine wave using SUPG-RK scheme . . . . . . . . . . . 63 3.12 Convergence of stabilised numerical schemes based on L1 norm . . . . . . 64 3.13 Convergence of stabilised numerical schemes based on L2 norm . . . . . . 65 3.14 Propagation of a step wave using Taylor-Galerkin scheme . . . . . . . . . . 67 3.15 Propagation of a step wave using SUPG-RK scheme . . . . . . . . . . . . 68 4.1 Evolution of pressure and velocity waves . . . . . . . . . . . . . . . . . . . 76 4.2 Evolution of pressure and velocity waves using Taylor-Galerkin Scheme at t = 0.5, C = 0.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Evolution of pressure and velocity waves using SUPG-RK Scheme at t = 0.5, C = 0.5 4.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Evolution of pressure and velocity waves using SUPG-RK scheme with varying τ at C = 0.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 82 89 Evolution of pressure and velocity waves using SUPG-RK scheme with varying τ at C = 0.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 5.1 Algorithm for solving 1D solid mechanics problem . . . . . . . . . . . . . 95 5.2 Algorithm for TVD Runge-Kutta time integration scheme . . . . . . . . . . 96 5.3 Algorithm for Taylor-Galerkin time integration scheme . . . . . . . . . . . 96 5.4 Axially-loaded cantilever bar . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.5 Time history of a bar under sinusoidal traction using TG scheme (form A) . 98 List of Figures 5.6 Time history of a bar under sinusoidal traction using TG scheme (form B) . 5.7 Time history of a bar under sinusoidal traction using SUPG-RK scheme ix 99 (form A) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.8 Time history of a bar under sinusoidal traction using SUPG-RK scheme (form B) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.9 Convergence of the Taylor-Galerkin scheme based on L2 norm for 1D fast dynamics problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.10 Convergence of SUPG-RK scheme based on L2 norm for 1D fast dynamics problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.11 Time history of a bar under constant traction using TG scheme (form A) . . 105 5.12 Time history of a bar under constant traction using SUPG-RK scheme (form A) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.13 Time history of a bar under constant traction using TG scheme (form B) . . 107 5.14 Time history of a bar under constant traction using SUPG-RK scheme (form B) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 5.15 Effect of varying τ at C = 0.4 . . . . . . . . . . . . . . . . . . . . . . . . 109 5.16 Comparison of numerical schemes . . . . . . . . . . . . . . . . . . . . . . 110 5.17 Time history of a bar under constant traction using TG scheme with shock capturing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.18 Time history of a bar under constant traction using SUPG-RK scheme with shock capturing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.19 Convergence of stabilised finite element schemes in the presence of shock capturing term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 5.20 Comparison of shock capturing schemes with Newmark scheme and Slope limiter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 List of Figures x 5.21 Time history of a bar under sinusoidal traction using TG scheme (NeoHookean material) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.22 Time history of a bar under sinusoidal traction using SUPG-RK scheme (Neo-Hookean material) . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.23 Time history of a bar under constant traction using TG scheme and lumped mass formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.24 Time history of a bar under constant traction using SUPG-RK scheme and lumped mass formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.1 Evolution of pressure distribution at different time instant for a short column 126 6.2 A tensile test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 6.3 Evolution of pressure distribution at different time instant for a tensile test . 128 Abbreviation WTW Water treatment works BG Bubnov-Galerkin w.r.t With respect to CFL Courant-Friedrichs-Lewy D.O.F Degrees of freedom FE Forward Euler FEM Finite element method L.H.S. Left hand side R.H.S. Right sand side RK Runge-Kutta SUPG Streamline Upwind Petrov-Galerkin TG Taylor-Galerkin TVD Total variation diminishing xi Notation English Symbols a Advection velocity h Length of element m Linear density p Pressure distribution in gas in case of acoustic wave propagation s Wave speed t Time u Unknown in 1d Hyperbolic equation Velocity distribution in gas in case of acoustic wave propagation A Cross section area C Courant-Friedrichs-Lewy number E Young’s modulus ET Total energy of a system G Amplification factor I Imaginary number J Determinant of deformation gradient tensor K0 Bulk modulus of compressibility N Unit normal direction Pe Peclet number R Spectral radius of Flux-jacobian matrix xii List of Figures T Traction force T0 Amplitude of traction force Up Longitudinal wave speed X Displacement Z non-dimensional parameter for selecting τ A Assembly operator b Left Cauchy-Green deformation tensor Body force d Rate of deformation tensor l Velocity gradient tensor p Linear momentum vector x Spatial coordinates of particles C Right Cauchy-Green deformation tensor E Green-Lagrange strain tensor F Deformation gradient tensor MeT G Consistent element mass matrix for the Taylor-Galerkin formulation MT G Consistent global mass matrix for the Taylor-Galerkin formulation MeSU P G Consistent element mass matrix for the SUPG formulation MSU P G Consistent global mass matrix for the SUPG formulation P First Piola-Kirchhoff stress tensor W Galerkin weighting function WSU P G Petrov-Galerkin weighting function X Material coordinates of particles Y Scaling matrix in Y Zβ Shock capturing scheme Z Residual in Y Zβ Shock capturing scheme A Flux-Jacobian matrix B Tensor containing derivative of element shape functions F Vector of flux xiii List of Figures L1 , L2 Left eigenvectors b L(•), L(•) Operator acting on • N Tensor containing element shape functions O(•) Order of truncation error of • R1 , R2 Right eigenvectors R (•) Residual of • S Source term U Vector of unknowns Uh Discretised version of U •t Partial time derivative of • J•K Jump in • •x Partial derivative of • w.r.t x Greek Symbols ∆t Time step size κ Shear modulus ν Poison’s ratio χ Stability parameter for selecting τ ∇0 Gradient with respect to the material configuration ∇ Gradient with respect to the spatial configuration ρ0 Constant material density λ Lame constant Λ Eigen values ζ Wavelength µ Lame constant φ Motion of a particle ξ, η Parameters denoting evolution of Amplification factor in complex plane xiv List of Figures ρ0 Constant material density τ Scalar stabilisation parameter τ Stabilisation parameter tensor ϑ, $ Parameters used in consistency analysis ω Frequency ϕ Phase angle D Diffusion error ϕ Dispersion error ωn Natural frequency β Parameter in Y Zβ Shock capturing scheme δ Shock capturing parameter γ1 , γ2 , γ3 Parameters to compute time step size in non-linear analysis xv Chapter 1 INTRODUCTION 1 INTRODUCTION 1.1 2 INTRODUCTION Study of behaviour of Solid continuum under external and body forces has been a challenging and active field of research. Early developments in the field of solid mechanics stemmed from civil engineering applications and currently has a wide range of applications from nano-scale to macro-scale. An important consideration in developing a numerical scheme is to able to represent the behaviour of a continuum. Fundamentally, the laws of continuum mechanics adheres to two types of description of motion, namely Lagrangian and Eulerian descriptions [12, 19, 35, 84]. The Lagrangian description is mainly used in the field of solid mechanics. In this description, the mesh motion coincides with the material motion. Hence, the history-dependent constitutive relation can be resolved naturally [12, 38] and the movement of free surface and material interfaces can be easily tracked. Lagrangian description naturally conserves mass. However, large distortions may emerge as a consequence of large strains in this description. The Eulerian description is widely popular in fluid mechanics and geodynamics applications. Here, the nodes and elements remain fixed in space, so that the continuum moves or deforms with respect to the computational mesh. In this description, the material interfaces and free surfaces can lose their accurate definitions [107]. This approach requires higher mesh density to capture the material response, making the method very computationally expensive. This method is usually unsuitable for solids as this description provides information only pertaining to current configuration. With both Lagrangian and Eulerian descriptions certain difficulties and advantages occur and on occasion it is possible to provide an alternative which attempts to secure the best features of both Lagrangian and Eulerian description by combining these. Such methods are known as Arbitrary-Lagrangian-Eulerian methods [106] and mostly used in fluid mechanics. There are two broad classes of external loads, namely static and dynamical loading [81]. Static forces are those that are applied slowly to a structure and influence of inertia forces INTRODUCTION 3 are not taken into account. In contrast, dynamic forces are time-varying forces which can trigger vibration of structures. Many engineering problems in which the dynamic effects are of particular importance are transportation, manufacturing and civil engineering structures under environmental loadings (i.e. wind and snow load). A detailed discussion on the various possible scenarios where dynamic loading is considered will be discussed in later part of this chapter. Traditionally, finite element method (FEM) is extensively used for problems in Computational Solid Mechanics [105]. Whereas, Finite Volume schemes are popular within the field of Computational Fluid Dynamics [47]. Both schemes can be considered as methods of weighted residuals where they differ in the choice of weighting functions [74]. The finite element Galerkin method treats the shape function as the weighting function and can be easily extended to higher order by increasing the order of polynomial interpolation. In contrast, the finite volume method results by selecting the weighting function as element piecewise unit constant. These two numerical schemes are equivalent in many applications [54]. 1.2 EARLY DEVELOPMENTS IN FEM The development of finite element method can be contributed to pioneering works of many researchers around the world. Perhaps, the first idea of finite elements was proposed by Henrikoff [45]. Henrikoff suggested an approach for solving continuum problems with proper boundary conditions which is fairly similar to the strategy that was used to resolve truss problems. Until 1952, finite element applications were mostly limited to elements that were connected by two points in space (e.g. rods, beams). In 1952, Ray Clough tackled the problem of modelling membranes or plates that were part of bending regions of a structure [31]. In 1956, Ioannis Argyris presented stress analysis results of an aircraft fuselages with many cut-outs, openings and severe irregularities [9]. This work demonstrated systematic implementation of finite element method, which later on paved way for immense popularity of this scheme in computer applications. In the same year, a complete solution framework INTRODUCTION 4 using the finite element method was reported [100] where triangular elements were considered to find solution to a plane stress problem. The well known direct stiffness method is proposed in this work. Further developments of this method can be credited to Zienkiewicz and Cheung [104] where they demonstrated that FEM is applicable to all problems that can be recast into a variational form. Earlier development of finite element method were mostly limited to infinitesimal or small deformation problems [41, 51, 68, 70]. Since, many problems arising in solid mechanics involve very large displacements, rotations and strains, the computational development has to incorporate both large strain continuum mechanics and non-linear material behaviour. Large strain theory, or finite strain theory also known as large deformation theory, deals with deformations in which both rotations and strains are arbitrarily large. In this case, the initial and deformed configurations of the continuum are significantly different and the stress distribution due to varying geometry can not be neglected. For problems dealing with dynamics, where the external loading increases or changes rapidly, the solution process are even more complex and inertia effects needs to be taken into account. Many applications in engineering require the solution of the dynamic response of structures or deformable solids. In the field of civil engineering, the structural response of a building under earthquakes loading (Figure 1.1), where, prior to the collapse, the structure undergoes large displacements and displays elasto-plastic response. Developing a finite element scheme capable of modelling this collapse accurately still poses a real challenge to engineers and researchers. However, a considerable amount of research work have been carried out in this direction. Dynamic response of systems of flexible beams and plates undergoing finite strains has been considered in [36, 55, 89, 90]. A stress resultant finite element formulation is proposed for the dynamic plastic analysis of plates and shells undergoing moderate to large deformations [28, 66]. Crash worthiness is one of the main research interest of automotive industry. The impact of a vehicle (Figure 1.2) involves fast dynamics and triggers geometrical non-linearity. Large INTRODUCTION 5 Figure 1.1: Collapse of a building under earthquake [1] strain finite element analysis of an impact can lead to accurate future design modifications that could ensure the safe design of the vehicle. Similarly, the metal forming processes (Figure 1.3), including rolling, forging, welding and extrusion of the material is relatively complex. These processes involve important changes in the original geometry, large strains, isochoric plastic flow, interaction (contact) with the forming tools and in many cases, selfcontact and a strong thermo-mechanical coupling [25, 29]. In field of biomedical engineering, the effect of impact loading human body and biological soft tissues using hyperelasticity and plasticity constitutive formulations have been studied [23, 69, 102] 6 INTRODUCTION Figure 1.2: Collision of two cars [1] 1.3 MOTIVATION A wide range of commercial finite element packages such as LS-DYNA3D [42], ANSYS [3], ABAQUS [5], SAP 2000 [6] and PAM-CRASH [4] have been developed for dynamic transient analysis. Majority of the packages use the Lagrangian description of motion. Lagrangian hydrocodes are also extensively used in automotive industry for resolving high velocity impact problems [14]. Various versions of the LS-DYNA code use an Updated Lagrangian description. Traditionally, Updated Lagrangian solid or structural mechanics codes such as ABAQUS and NASTRAN [2] use an implicit time integration scheme for evolving in time and are found not to be suitable for dynamic analyses in which shock waves are dominant [67]. Since, the variables in the Updated Lagrangian are evaluated at current configuration, which is constantly changing, finite element mesh coordinates need to be recomputed at every time step. This leads to a significant rise in computer storage and computational time [72]. In current work, the Total Lagrangian description is adopted, 7 INTRODUCTION Figure 1.3: Welding of a metal [1] where formulations are always with respect to the reference configuration. The traditional solid dynamics formulation, where displacement field is the primary variable, is solved by standard finite element spatial discretisation coupled with a family of Newmark time integration schemes. However, the resulting space-time discretised formulation presents a series of shortcomings. Firstly, Newmark’s method has a tendency for high frequency noise to persist in the solution and most importantly, its accuracy is remarkably degraded once artificial damping is implemented. Some minor modifications were introduced to improve the accuracy of numerical dissipation without the inclusion of a discontinuity sensor, which consequently made the Newmark scheme unsuitable for problems with shock discontinuity [7, 30, 46, 103]. Moreover, it is well known that using linear elements INTRODUCTION 8 in displacement-based FEM leads to second order convergence for displacements but one order less for strains and stresses. Constant stress elements exhibit volumetric locking in incompressible or nearly incompressible applications; for instance, plastic flows involving large isochoric strains [20]. In order to alleviate the locking phenomena, different approaches have been proposed. One of the approaches recommends adoption of p-refinement along with high order interpolating functions [37]. Another general strategy is to introduce a multi-field Veubeke-Hu-Washizu type variational principle, which enables the use of independent kinematic descriptions for the volumetric and deviatoric deformations [19]. The mean dilatation technique, in which a constant interpolation for volumetric variables over an element is involved, is widely accepted. This specific technique can be identified as a particular case of Selective Reduced Integration, where the volumetric stress components are integrated strategically using less number of gauss points. Unfortunately, this scheme is not applicable for low order elements (i.e. linear triangles and linear tetrahedrons) as these elements have already used the simplest Gaussian quadrature rule. Bonet and Burton [16] suggested that the volumetric strain energy is approximated by evaluating averaged nodal pressures in terms of nodal volumes while the deviatoric component is treated in a standard manner. However, the resulting solution performed poorly in bending dominated cases. To circumvent this issue, Dohrmann et al. [37] proposed a new linear tetrahedron by employing nodal averaging process to the whole small strain tensor. Subsequently, Bonet et al. [18] extended this application to large strain regime with the idea of implementing an averaged nodal deformation gradient tensor as the main kinematic variable. It has been recently been observed that the explicit solid dynamics formulation can be expressed as a system of first order conservation laws similar to the system in wave propagation problems [17, 62, 73, 75]. An effort has been made to solve the first order conservation laws using Discontinuous Galerkin method [75]. Additionally, the equation for the rate of change of the deformation gradient tensor has also been added. The same system of first order conservation laws has also been solved using two-step Taylor-Galerkin method [58, 79] INTRODUCTION 9 and cell-centerd finite volume method [61] for Lagrangian solid dynamics. 1.4 OBJECTIVES AND SCOPE OF WORK The objectives of the current study can be summarised as follows: • Propose stabilised finite element schemes to solve first order hyperbolic system of conservation laws. • Investigate the consistency and stability of the schemes. • Perform Spectral analysis of the schemes. • Verify the order of convergence of the stabilised schemes. • Investigate the suitability of different constitutive models, for instance, linear elastic, non-linear Neo-Hookean models. • Investigate the suitability of lumped mass formulation. • Implement shock capturing scheme in order to minimise spurious oscillations in the vicinity of shocks. • Compare the accuracy of the stabilised finite element schemes with first and second order finite volume schemes. • Address issues regarding multi-dimensional implementation of the solution of the conservation laws. The effect of plasticity has not been studied in current work. However, the proposed numerical framework permits the inclusion of any suitable plasticity models. This work briefly introduces challenges arising from multi-dimensional implementation. However, further study is needed in order to eliminate all the problems satisfactorily. INTRODUCTION 1.5 10 LAYOUT OF THESIS This work is divided into seven chapters including this chapter which reviews the earlier research in FEM and introduces the motivation and objectives of current research. It also provides the scope of the ongoing work. Chapter 2 reviews the concept and terminologies used in Kinematics, which are used in later part of the thesis. It also provides the formulation of the governing first order conservation laws. The constitutive equations for the material models used in this study are also presented. Finally, the contact and shock conditions corresponding to the Riemann solution on the boundaries as well as interfaces between elements are described. Chapter 3 proposes two stabilised finite element schemes e.g. one-step Taylor-Galerkin scheme and the Streamline Upwind Petrov Galerkin scheme (SUPG) for solving the pure advection equation. The stability and consistency of the resulting space-time schemes are evaluated. Moreover, Spectral analysis is performed in order to determine the nature of the dispersion and diffusion error. Finally, the Riemann problem is introduced and few numerical simulations are conducted. In chapter 4, the stabilised finite element formulations are extended in order to solve one dimensional (1D) fast dynamics problems. Two different finite element framework, namely, non-conservative and conservative framework, are introduced. The formulations for two possible ways to obtain numerical solution, namely, “form A” and “form B”, for each framework, are presented. The analytical and numerical solutions for propagation of an acoustic wave in an 1D tube is discussed. Finally, few numerical experiments are performed with regards to the SUPG scheme. Chapter 5 focuses on the numerical implementation to solve 1D fast dynamics problems. It briefly introduces the governing equations for the problem and presents results obtained for an 1D example for both “form A” and “form B”. The convergence of the proposed schemes are also verified in this chapter. Moreover, it proposes an shock capturing scheme for ensuring better solution in the presence of shock discontinuity. The accuracy of the INTRODUCTION 11 stabilised finite element schemes is compared with first and second order finite volume schemes. The suitability of the formulation for incorporating different constitutive models is investigated. Finally, the feasibility of adopting a lumped mass formulation is reported. Chapter 6, briefly introduces the challenges arising from multi-dimensional implementation of fast dynamics problems. The initial steps to circumvent these issues are presented. A numerical example is presented and compared with finite volume schemes. Finally, chapter 7 summarises the conclusions of the current work and addresses areas of future research concerning the total Lagrangian solid dynamics finite element formulation. Chapter 2 LARGE STRAIN STRUCTURAL DYNAMICS 12 LARGE STRAIN STRUCTURAL DYNAMICS 2.1 13 INTRODUCTION This chapter presents an overview of large strain solid dynamics and introduces some of the notations and terminologies that will be used in the later part of this thesis. The chapter begins with discussing few fundamental concepts of kinematics. The discussion is followed by the new mixed formulation for fast dynamics. Later, the constitutive models utilised in this work have been presented. Afterwards, the eigenstructure of the governing equation is introduced. Finally, the contact and shock conditions are discussed. 2.2 KINEMATICS PRELIMINARIES Kinematics can be defined as the science behind motion and deformation without any reference to the forces causing it [19]. Traditionally, behaviour of finite deformation is associated with two alternative description settings, namely, Material or Lagrangian description and Spatial or Eulerian description. This section introduces some of the basic concepts of kinematics. An extensive discussion on kinematics is available in references [12, 19, 35, 48, 84]. 2.2.1 The Motion Let B be a body with volume V , which consists of particles X at time t = 0. The motion can be defined as the continuous time sequence of displacements which carries the body B into current time t with volume v(t), containing the particles at x. In mathematical terms x(t) = φ (X, t) . (2.1) Alternatively, it can also be stated that φ : V → v (t) where ∀X ∈ V, ∃x ∈ v(t)/x(t) = φ (X, t) . (2.2) 14 LARGE STRAIN STRUCTURAL DYNAMICS The motion of B is shown in Figure 2.1. If the physical and kinematic quantities are exX 3 ,x3 v t time t V time t 0 x X X 2 ,x2 X1 ,x1 Figure 2.1: Motion of a deformable body pressed in terms of where the body was before deformation, it is known as material or Lagrangian description. Alternatively, if the relevant quantities are described where the body is after the deformation, it is called as spatial or Eulerian description. 2.2.2 Deformation Gradient Tensor The deformation gradient tensor F is key to the description of deformation and hence strain. It is expressed as the parital derivative of the mapping φ (X, t) with respect to the initial configuration as, F = ∂x ∂φ (X, t) = = ∇0 x, ∂X ∂X (2.3) where ∇0 indicates the gradient with respect to the material configuration. The elemental vector in spatial configuration dx can be obtained in terms of elemental vectors in initial configuration dX as dx = F dX. (2.4) LARGE STRAIN STRUCTURAL DYNAMICS 15 Since, F transforms vectors in initial configuration to vectors in spatial configuration, it is known as a two-point tensor. 2.2.3 Volume Change The elemental spatial volume can be expressed as dv = JdV, (2.5) where dV is the volume at material configuration and Jacobian J is defined as the determinant of the deformation gradient tensor as, J = det(F ). (2.6) The element mass dm can be written as dm = ρ0 V = ρv, (2.7) where ρ0 and ρ indicate the densities in material and spatial configurations, respectively. From equations (2.5) and (2.7), it can be deduced that ρ0 = Jρ. Equation (2.8) is known as conservation of mass or continuity equation. (2.8) LARGE STRAIN STRUCTURAL DYNAMICS 2.3 16 Strain In order to measure strains in large deformation, a suitable strain tensor E referred to as the Green-Lagrange strain tensor, is introduced E= 1 (C − I) , 2 (2.9) where C is the right Cauchy-Green deformation tensor, which is expressed in terms of F as C = FTF. (2.10) It is worth noticing that C operates on material elemental vectors and hence, it is a material tensor. The right Cauchy-Green deformation tensor can also be written as C= 3 X Λ2α Nα ⊗ Nα , (2.11) α=1 where Nα represents the principal material directions of C and their corresponding eigenvalues are denoted by Λα . The spatial measure of strain is known as left Cauchy-Green deformation tensor and is related to F as b = FFT. (2.12) Similar to C, b can also be expressed in terms of principal spatial directions nα and associated eigenvalues Λα as b= 3 X Λ2α nα ⊗ nα . (2.13) α=1 Therefore, the relationship between Nα and nα can be written as F Nα = Λα nα . (2.14) 17 LARGE STRAIN STRUCTURAL DYNAMICS 2.3.1 Velocity Measures The Lagrangian description of velocity of a particle is defined as v (X, t) = ∂φ (X, t) . ∂t (2.15) The velocity can be written in a spatial position x, by inverting equation (2.15) as v (x, t) = v φ−1 (x, t) , t . (2.16) The partial derivative of equation (2.16) with respect to the spatial configuration leads to the definition of velocity gradient tensor l as l= ∂v (φ−1 (x, t) , t) = ∇v, ∂x (2.17) where ∇ indites the gradient with respect to the spatial configuration. Therefore, the time derivative of deformation gradient tensor in terms of the velocity gradient can be written as Ḟ = 2.3.2 ∂v ∂φ (X, t) = ∇vF = lF . ∂x ∂X (2.18) Rate of Deformation The time derivative of Lagrangian strain tensor is known as the material strain rate tensor and can be easily obtained in terms of Ḟ as 1 1 T Ė = Ċ = Ḟ F + F T Ḟ . 2 2 (2.19) The rate of deformation tensor d can be expressed as d= 1 l + lT , 2 (2.20) 18 LARGE STRAIN STRUCTURAL DYNAMICS which is the symmetric part of velocity gradient tensor. It can also be observed that d = F −T ĖF −1 . 2.4 (2.21) GOVERNING EQUATIONS In this section, the governing conservation equations or balance laws, in the context of Lagrangian dynamics, are presented. The obvious advantage of this formulation is that all derivatives with respect to spatial coordinates are calculated based upon an original undeformed configuration. The new mixed formulation combines conservation of linear momentum, deformation gradient and energy [17, 58, 61, 63, 62, 73, 79]. Similar formulations exist in an Eulerian configuration [71, 99]. 2.4.1 Conservation of Linear Momentum For a continuum, conservation of linear momentum dictates that the rate of change of linear momentum of particles in material configuration is equal to resultant forces applied to these particles. Mathematically, this principle can be written as d dt Z Z t dA, ρ0 b(X, t) dV + p(X, t) dV = V Z V (2.22) ∂V where p(X, t) = ρ0 v(X, t) is the linear momentum per unit of material volume, ρ0 represents the material density, v is the velocity field, b stands for the body force per unit mass and t denotes the nominal traction vector. The traction force t can be related to the first Piola-Kirchhoff stress tensor as t = P N, (2.23) 19 LARGE STRAIN STRUCTURAL DYNAMICS where N indicates the unit normal in reference configuration. Applying the divergence theorem on the last term at right hand side yields Z V dp(X, t) dV = dt Z Z ∇0 · P dV. ρ0 b(X, t) dV + V (2.24) V It is worth mentioning that the derivative of the integral term at the L.H.S. of equation (2.22) can rewritten as the integral of a derivative term in equation (2.24) as the integration takes place over material volume which does not change over time. Since equation (2.24) is valid for any arbitrary volume, the local form of conservation of linear momentum is dp(X, t) = ρ0 b(X, t) + ∇0 · P . dt 2.4.2 (2.25) Conservation of Deformation Gradient In order to alleviate shear locking as well as volumetric locking, it is useful to treat F as an independent variable with the aim of increasing the degrees of freedom (or flexibility) of the problem [18]. Conservation of deformation gradient tensor can be easily derived by noting that the time derivative of F (X, t) is related to the linear momentum p(X, t) as ∂F ∂ = ∂t ∂X ∂φ (X, t) ∂t = ∇0 p ρ0 . (2.26) With the help of the identity tensor I, equation (2.26) can be alternatively written as ∂F = ∇0 · ∂t p ⊗I . ρ0 (2.27) Considering the integral form of equation (2.27) over volume V and applying divergence theorem to the right hand side Z V ∂F dV = ∂t Z ∂V p ⊗N ρ0 dA. (2.28) 20 LARGE STRAIN STRUCTURAL DYNAMICS Equation (2.28) can be considered as the generalization of continuity equation usually employed in fluid mechanics. 2.4.3 Conservation of Energy The conservation of the total energy dictates that, in a continuum, the change in total energy of a material volume should be equal to sum of the work done and heat exchanged in the system. Mathematically, it can be expressed as d dt Z Z Z t · v dA − ET dV = V ∂V Q · N dA, (2.29) ∂V where ET is the total energy per unit of undeformed volume, t describes the traction vector, v stands for the velocity vector, Q denotes the heat flow vector and N represents the outward-pointing unit normal vector in reference configuration. For simplicity, the heat source term is ignored. Applying divergence theorem to convert surface integral in to volume integral leads to d dt Z Z ET dV = V V T p ∇0 · P − Q dV. ρ0 (2.30) The local differential form of equation (2.30) is given as ∂ET = ∇0 · ∂t 1 T P p−Q . ρ0 (2.31) The total energy ET can be expressed as ET = e + Ψext + 1 p·p 2ρ0 (2.32) where e indicates the internal energy in unit material volume consisting of heat and elastic energy components, Ψext denotes the potential energy resulting from external forces ρ0 b and the last term expresses kinetic energy. Assuming that the external forces remain constant, 21 LARGE STRAIN STRUCTURAL DYNAMICS the potential energy can be written as Ψext = −ρ0 b · x, (2.33) Ψ̇ext = −ρ0 b · v. (2.34) and the time derivative leads to Therefore, equation (2.31) can modified by combining with equation (2.27) as ∂e ∂F =P : − ∇0 · Q. ∂t ∂t 2.4.4 (2.35) Conservative Law Formulation The physical laws for the linear momentum and the deformation gradient tensor are summarised here for convenience: ∂p − ∇0 · P = ρ0 b, ∂t p ∂F − ∇0 · ( ⊗ I) = 0, ∂t ρ 0 ∂ET 1 T − ∇0 · P p − Q = 0. ∂t ρ0 (2.36) (2.37) (2.38) These conservation laws can then be grouped into a single system of first order hyperbolic equations as ∂U ∂F I + = S, ∂t ∂XI ∀ I = 1, 2, 3; (2.39) 22 LARGE STRAIN STRUCTURAL DYNAMICS where their components are illustrated as follows p1 p 2 p 3 F 11 F 12 F 13 U = F21 , F 22 F 23 F 31 F 32 F 33 E T FI = −P1I (F ) −P2I (F ) −P3I (F ) −δI1 p1 /ρ0 −δI2 p1 /ρ0 −δI3 p1 /ρ0 −δI1 p2 /ρ0 −δI2 p2 /ρ0 −δI3 p2 /ρ0 −δI1 p3 /ρ0 −δI2 p3 /ρ0 −δI3 p3 /ρ0 QI − (PiI pi ) /ρ0 , S= ρ 0 b 1 ρ b 0 2 ρ b 0 3 0 0 0 0 For a reversible process, the energy term can be solved independently. 0 0 0 0 0 0 . (2.40) 23 LARGE STRAIN STRUCTURAL DYNAMICS 2.4.5 Interface Flux At a given interface defined by the outward unit normal vector material configuration N = [N1 N2 N3 ]T , the interface flux without the energy term is expressed as FN −t1 (F ) −t (F ) 2 −t (F ) 3 −N p /ρ 1 1 0 −N p /ρ 2 1 0 −N3 p1 /ρ0 = F I NI = , −N p /ρ 1 2 0 −N p /ρ 2 2 0 −N p /ρ 3 2 0 −N1 p3 /ρ0 −N p /ρ 2 3 0 −N3 p3 /ρ0 ∀ I = 1, 2, 3, (2.41) with the help of t = P N . 2.5 HYPERELASTICITY Traditionally, materials for which constitutive equations depend solely on the current state of deformation are known as elastic. Hence, the stress at a particle X in body B can be measured as a function of the current deformation gradient F associated with that particle. Consequently, the first Piola-Kirchhoff stress P can be measured in terms of its conjugate F as P = P (F (X), X) , (2.42) 24 LARGE STRAIN STRUCTURAL DYNAMICS where the direct dependence on X takes into account the possible non-homogeneous nature of B. Hyperelasticity occurs when work done by the stresses during a deformation process is dependent only on the initial state at time t0 and the final configuration at time t. This path-independent behaviour allows the stored strain energy function or elastic potential per undeformed volume Ψ to be captured in terms of P and its work conjugate Ft as Z t P (F (X), X) : Ft dt; Ψt = P : Ft , Ψ (F (X), X) = (2.43) t0 which leads to P (F (X), X) = ∂Ψ (F (X), X) . ∂F (2.44) Behaviour of all hyperelastic materials are governed by equation (2.43) and equation (2.44). This general constitutive equation can be further extended by observing that as a consequence of the objectivity requirement, Ψ must be independent of the rotation component of the deformation gradient. Therefore, Ψ can be expressed as a function of the right CauchyGreen tensor C as Ψ (F (X), X) = Ψ (C(X), X) . (2.45) Since, 12 Ct = Et is work conjugate to the second Piola-Kirchoff stress S, Lagrangian constitutive equation for hyperelastic materials can be formulated as Ψt = ∂Ψ 1 ∂Ψ ∂Ψ : Ct = S : Ct ; S(C(X), X) = 2 = . ∂C 2 ∂C ∂E (2.46) For isotropic materials, Ψ can be expressed in terms of the principal invariants of F as Ψ (F (X), X) = Ψ (IF , IIF , IIIF , X) , (2.47) where IF = tr(F ), IIF = F : F , IIIF = det(F ). Therefore, P can be written in terms 25 LARGE STRAIN STRUCTURAL DYNAMICS of the invariants as P = ∂Ψ ∂Ψ ∂IF ∂Ψ ∂IIF ∂Ψ ∂IIIF = + + , ∂F ∂IF ∂F ∂IIF ∂F ∂IIIF ∂F (2.48) where the derivatives of the invariants are computed as ∂IF = I; ∂F ∂IIF = 2F ; ∂F ∂IIIF = det(F )F −T . ∂F (2.49) The stored energy function Ψ can be decomposed into the summation of deviatoric Ψdev (J −1/3 F ) and volumetric Ψvol (J) components as Ψ(F ) = Ψdev (J −1/3 F ) + Ψvol (J), (2.50) which in turn, leads to P = Pdev + Pvol ; 2.5.1 Pdev = ∂Ψdev ; ∂F Pvol = ∂Ψvol . ∂F (2.51) Nearly Incompressible Neo-Hookean Material The simplest model satisfying the conditions described in previous section is the nearly incompressible Neo-Hookean (NH) material. Its deviatoric and volumetric parts are described as Ψdev = 1 µ [J −2/3 (F : F ) − 3]; 2 1 Ψvol = κ(J − 1)2 . 2 (2.52) Here, κ is the bulk modulus which only appears in the volumetric term whereas the shear modulus µ, on the other hand, appears in the deviatoric counterpart. The corresponding stress components can be evaluated as 1 Pdev = µJ −2/3 [F − (F : F )F −T ]; 3 Pvol = dΨvol ∂J = κ(J − 1)JF −T . dJ ∂F (2.53) 26 LARGE STRAIN STRUCTURAL DYNAMICS The fourth-order constitutive tensor is defined as C= 2.5.2 ∂P . ∂F (2.54) Linear Elasticity A linearised elastic constitutive relationship is considered as an excellent model to describe small deformation behaviour for engineering materials such as, concrete, steel and metal. In this type of material, Ψ is defined as 1 ψ (ε) = λ (tr(ε))2 + µ (ε : ε) , 2 (2.55) where µ and λ are the so-called Lamé constants. It is worth mentioning that Saint-Venant Kirchhoff material is recovered if ε is replaced by the Green-Lagrange strain tensor E. Deformation gradient tensor can be split into a displacement gradient H = ∂u/∂X and identity tensor I, i.e. F = I + H. In the context of infinitesimal strain, an assumption is made such that only linear contributions of H are considered. In what follows, the engineering strain ε and its trace can be further developed as ε= 1 1 H + HT = F + F T − 2I ; 2 2 tr(ε) = tr(H) = tr(F ) − 3. (2.56) In the absence of deformation (F = I), the stored energy functional vanishes as expected (Ψ(ε) = 0). Based on equation (2.44), after few algebraic manipulations, the stress tensor is obtained as 2 P (F ) = µ F + F − tr(F )I + κ (tr(F ) − 3) I. 3 T (2.57) 27 LARGE STRAIN STRUCTURAL DYNAMICS 2.6 EIGENSTRUCTURE The mixed conservation law expressed in equation (2.39) can be rewritten in non-conservation form as ∂U ∂U I + AI = S, ∂t ∂XI ∀I = 1, 2, 3, (2.58) where Flux Jacobian matrix AI can be expressed as AI = ∂F I . ∂U (2.59) The Flux Jacobian matrix at a given interface can be written as AN = AI NI = ∂F I ∂F N NI = . ∂U ∂U (2.60) In order to comprehend the eigenstructure of the Flux Jacobian matrix, it is useful to separate the momentum and the deformation gradient tensor of U and F N as U= p and F N = −t , (2.61) −v ⊗ N F where the term corresponding to energy is neglected for simplicity (i.e., no heat flow takes place). In equation (2.61), F should be interpreted as column vectors of 9 entries corresponding to each of its components. Consequently, AI can be expressed as AN = − − ∂(P∂pN ) ∂ ρ1 p ⊗N 0 ∂p − N) − ∂(P ∂F ∂ ρ1 p ⊗N 0 ∂F 03×3 −C N = , − ρ10 IN 09×9 (2.62) 28 LARGE STRAIN STRUCTURAL DYNAMICS where the normal component of fourth order constitutive tensor C = [C N ]ijJ = ∂PiI NI ∂FjJ and ∂P ∂F is indicated as [IN ]iIk = δik NI . (2.63) The right and left eigenvectors of AN , namely Rα and Lα , and their corresponding eigenvalues Uα can be related as AN Rα = Uα Rα (2.64a) LTα AN = Uα LTα . (2.64b) The orthogonality condition between the left and right eigenvectors dictates that 12 X Rα LTα . Uα T AN = Rα Lα α=1 (2.65) In order to derive expressions for these eigenvectors, it is crucial to separate their components into Rα = pR α FαR , Lα = pL α . (2.66) FαL Substituting the explicit expression for AN in equation (2.62) into equation (2.64a) leads to − −C N : FαR = Uα pR α (2.67a) 1 R p ⊗ N = Uα FαR . ρ0 α (2.67b) Eliminating FαR by inserting equation (2.67b) into equation (2.67a) yields a symmetric eigenvalue problem for pR as 2 R C N N pR α = ρ0 Uα pα , (2.68) 29 LARGE STRAIN STRUCTURAL DYNAMICS where the symmetric 3 × 3 tensor C N N is given as [C N N ]ij = 3 X CiIjJ NI NJ . (2.69) I,J In the nonlinear elastic context, the eigenproblem discussed above yields 3 pairs of wave speeds, which correspond to the volumetric (or P-wave) Up and shear (or S-wave) Us , where s Up = U1,2 = ±Up , (2.70a) U3,4 = U5,6 = ±Us , (2.70b) γ2 + γ1 Λ2 + 2γ3 ρ0 r ; Us = γ2 , ρ0 (2.71) where 5 γ1 = κJ 2 + µJ −2/3 (F : F ) , 9 (2.72a) γ2 = µJ −2/3 , (2.72b) 2 γ3 = − µJ −2/3 , 3 1 Λ= . −T kF N k (2.72c) (2.72d) Expression (2.70) concludes that the remaining six eigenvalues of matrix AN are zero. In linear elasticity context, since F ≈ I and J ≈ 1, the longitudinal and shear waves can be simplified as s Up = λ + 2µ ; ρ0 r Us = µ . ρ0 (2.73) The matrix AN can therefore be reconstructed in terms of non-zero wave speeds as 6 X Rα LTα . AN = Uα T Rα Lα α=1 (2.74) Moreover, the eigenvalue structure, also leads to 3 pairs of orthogonal eigenvectors, where 30 LARGE STRAIN STRUCTURAL DYNAMICS the first one n corresponds to the outward unit normal vector in spatial configuration associated to material vector N and the remaining two are arbitrary tangential vectors t1,2 orthogonal to n. These orthogonal eigenvectors are written as R1,2 = n ± ; R3,4 = 1 n ρ0 Up R5,6 = ⊗ N t2 1 t ρ 0 Us 2 ⊗ N ± ± t1 1 t ρ 0 Us 1 ⊗ N ; . (2.75) Subsequently, the set of left eigenvectors is obtained as L1,2 = n ± ; L3,4 = t1 ; 1 C Up L5,6 ± 1 C : (t1 ⊗ N ) : (n ⊗ N ) Us t2 . = ± 1 C : (t2 ⊗ N ) Us (2.76) Noting that RTα Lα = 2 for α = 1, 2, . . . , 6, the Flux Jacobian matrix AN can now be rewritten as 6 1X AN = Uα Rα LTα 2 α=1 (2.77a) 0 0 0 0 0 Up 0 −U 0 0 0 0 p T L1 0 0 Us 0 0 0 1 .. = {R1 , . . . , R6 } . 2 0 0 0 −U 0 0 s T L6 0 0 0 0 Us 0 0 0 0 0 0 −Us = 1 Up (R1 LT1 − R2 LT2 ) + Us (R3 LT3 − R4 LT4 + R5 LT5 − R6 LT6 ) . 2 (2.77b) (2.77c) 31 LARGE STRAIN STRUCTURAL DYNAMICS This expression suggests existence of 3 sets of waves being transmitted in opposite direction: one set of longitudinal waves travelling with speeds Up and −Up in the direction n as well as 2 sets of shear waves moving with speeds Us and −Us in the directions t1 and t2 . 2.7 SHOCK AND CONTACT CONDITIONS The conservation law in equation (2.39) can be written in integral form (neglecting source terms) as d dt Z Z U dV + V V ∂F I dV = 0. ∂XI (2.78) F N d∂V = 0. (2.79) Application of divergence theorem leads to d dt Z Z U dV + V ∂V These conservation laws allow solutions with discontinuities or shocks propagating through the medium at certain speed. Jump conditions are derived by considering an arbitrary volume of deformable body shown in Figure 2.2. The domain is divided by a jump discontinu- U+ V+ N ∂V − U− V− ∂V + U Γ Figure 2.2: Surface discontinuity ity across surface that separates it into two regions V + and V − with boundaries ∂V + and ∂V − , respectively. The speed of the moving surface Γ with unit normal N is denoted as U . 32 LARGE STRAIN STRUCTURAL DYNAMICS Application of the Reynold’s transport theorem yields Z Z Z d ∂U U dV = dV + −U U + d∂V dt V + ∂t ZΓ Z ZV + d ∂U dV + U U − d∂V. U dV = dt V − ∂t Γ V− (2.80) Summation of both equations leads to the Reynold’s transport theorem with jump discontinuity as d dt Z Z U dV = V V ∂U dV + ∂t Z U J U K d∂V, Γ (2.81) where J U K indicates jump in U as J U K = U + − U −. (2.82) Subsequently, the expression for Flux term in equation (2.78) can be written as Z V ∂F I dV = ∂XI Z Z F N d∂V + ∂V Γ J F N K d∂V, (2.83) Adding equation (2.81) and (2.81) yields U J U K = J F N K. (2.84) This is so called Rankine-Hugoniot condition [80, 98]. The generalised jump condition derived in equation (2.84) can be utilised to formulate the jump conditions for momentum, deformation gradient and energy as U J p K = −J P KN , 1 J p K ⊗ N, ρ0 1 U J ET K = − J P T p K · N . ρ0 U JF K = − (2.85) 33 LARGE STRAIN STRUCTURAL DYNAMICS Contact conditions play a major role in analysing dynamic problems. Figure 2.3 displays motion of an arbitrary body which comes into contact with another body. Physically, this φ− , p − , F − v− (t ) time t U s− U p− U p+ U s+ n v+ (t ) V− N− X 2 ,x2 φ+ , p+ ,F + N+ V+ time t = 0 N = N − = −N + n = n − = −n+ X 1 ,x1 Figure 2.3: Contact motion generated shock waves can be the result of impact between two bodies or two parts of the same body. Numerically, contacts may arise from the use of discontinuous interpolations of the problem variables, such as in Godunov’s type of methods or discontinuous Galerkin. In current formulation, A general interface flux at a contact point can be expressed as FC N = −t C (2.86) − 1 pC ⊗ N ρ0 where the term corresponding to energy has been neglected. The traction and linear mo- 34 LARGE STRAIN STRUCTURAL DYNAMICS mentum term dependant on various contact conditions (Figure 2.4) can be expressed as pC = 0, − pC t = pt , pC n = 0, pC = p− , tC = t− (sticking surface) B tC t = tt , − tC n = tn (sliding surface) (2.88) (free surface). (2.89) tC = tB tB ttB ି ݒሺݐሻ (a) Sticking surface (2.87) ି ݒሺݐሻ (b) Sliding surface ି ݒሺݐሻ (c) Free surface Figure 2.4: Different boundary conditions 2.8 CONCLUSION Various concepts pertaining to Kinematics are introduced in this chapter. Among these, the terms pertaining to deformation gradient tensor plays a major role in later part of this thesis. Moreover, the conservation laws of deformation gradient tensor, linear momentum and energy are formulated. The constitutive models of nearly incompressible Neo-Hookean and linear elastic materials are also presented in this chapter. Furthermore, the expressions for longitudinal and shear wave velocities for different material models are obtained. The Rankine-Hugoniot condition is derived and extended to the conservation variables. Finally, the expression for contact for sticking surface, sliding surface and free surface are presented. Chapter 3 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 35 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 3.1 36 INTRODUCTION From the discussion of the previous chapter, it is known that the set of equations arising from the new mixed Lagrangian solid dynamics problems are hyperbolic in nature. This chapter aims to envisage different finite element formulations to solve the one-dimensional (1D) hyperbolic equation. Later, these formulations can be extended to multi-dimensions. This chapter first introduces the one step Taylor-Galerkin formulation for solving the pure advection equation. This discussion is followed by Streamline Upwind Petrov-Galerkin (SUPG) formulation for spatial discretisation. Later, various time integration schemes are employed for obtaining the solution of the 1D linear advection equation, which is hyperbolic in nature. Various numerical analyses such as, consistency analysis, stability analysis and spectral analysis are performed on the schemes. Few numerical examples are presented in order to demonstrate the capabilities and efficiency of the schemes. Finally, some conclusions are reached on the proper choice of numerical schemes as well as their consistency and stability. 3.2 TAYLOR-GALERKIN SCHEME In this particular scheme, the advection equation is discretised in space using Galerkin approximation and in time with a second order Taylor series discretisation in time. The advection equation (3.12) at time step n can be written as unt = −aunx . (3.1) The second derivative with respect to time yields utt = −aunxt = −auntx = a2 unxx . (3.2) 37 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION It is known from a Taylor series expansion that ∆t n ∆u un+1 − un := = unt + utt + O ∆t2 ∆t ∆t 2 (3.3) Substituting expressions for first and second order time derivatives (equations (3.1) and (3.2)) into equation (3.3) ∆u a2 ∆t n = −aunx + uxx + O ∆t2 . ∆t 2 (3.4) Application of Galerkin discretisation to equation (3.4) results into ∆uh W dΩ = − ∆t Ω Z Z W auhx Ω Z dΩ + W Ω ∆t 2 h a uxx dΩ. 2 (3.5) Application of the divergence theorem to the R.H.S. leads to ∆uh W dΩ = ∆t Ω Z Z h Z Wx au dΩ − Ω Ω ∆t 2 a Wx uhx dΩ + 2 Z W Z∂Ω − ∆t 2 h a ux d∂Ω 2 W auh d∂Ω. (3.6) ∂Ω The discretised version of above equation is nel X Z nel X ∆ũ dΩ = W̃ · NN W̃ · aBN T ũ dΩ ∆t e e Ω Ω e=1 e=1 Z n el X a2 ∆t ∆t W̃ · − BBT ũ dΩ + W̃N a2 uxN − W̃N auN . 2 2 Ωe e=1 Z T (3.7) The discretised form of equation (3.7) can be written as MT G ∆ũ + KT G ũ = 0, ∆t (3.8) 38 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION where nel MT G = A nel e=1 (MeT G ) , KT G = A e=1 nnode (KeT G ) and ũ = A i=1 (ũi ) . (3.9) In equation (3.9), A represents the assembly operator, ũi = [ui , ui+1 ]T and consistent element mass matrix MeT G is expressed as MeT G and KeT G h 2 1 = 6 1 2 a −1 −1 ∆ta2 =− + 2 1 2h 1 (3.10) 1 −1 . −1 1 (3.11) The one step Taylor-Galerkin scheme is second order accurate in time and space. 3.3 DISCRETISATION IN SPACE: SUPG The one-dimensional linear advection equation can be expressed as ut + aux = 0 in Ω, (3.12) with boundary and initial conditions u = ů in (Ω, t0 ), u = uD on (∂ΩD , t), ux = uN on (∂ΩN , t), (3.13) where a is the advection velocity, u denotes an unknown scalar variable, ut indicates the partial derivative of u with respect to time (t), t0 indicates initial time, ux represents the partial derivative of u with respect to space, Ω is a bounded open domain in < with Lipschitz SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 39 continuous boundary ∂Ω = ∂ΩD ∪ ∂ΩN with Neumann and Dirichlet boundary conditions denoted by uN , uD , respectively. As discussed earlier, the standard Bubnov-Galerkin solution of hyperbolic equations leads problems such as spurious oscillations. In order to obtain a stabilised finite element solution, time-space elements are a natural choice [56, 85, 86]. Here, upwinding effect can be incorporated by combining the time derivative and the advective term into a ”single material” derivative, which leads to Petrov-Galerkin weighting function (WSU P G ) [15, 56] in one-dimension as Duh = u̇h + auhx , WSU P G uh = Dt (3.14) Nevertheless, the most effective and popular approach is to first discretise the hyperbolic equation in space, which leads to simultaneous first order ordinary differential equations dependent only on time. Next, the equations are solved using appropriate time integration scheme. As a result, these algorithms allow usage of existing spatial finite element frameworks and apply a suitable time integration scheme without significant development of new software [59, 24, 34, 44, 49]. This technique is also known as Method of Lines [65]. Moreover, in practical applications the increased computational cost due to higher number unknowns for coupled time-space formulations is a significant drawback. A popular and efficient way to circumvent the problem arising from spurious oscillations and locking, is to augment the Bubnov-Galerkin form of equation (3.12) with terms which add numerical dissipation but diminish for all sufficiently smooth solutions. Resulting schemes (e.g., SUPG method, Taylor-Galerkin method) are called consistently stabilised methods [15] as the order of accuracy of Galerkin approximation is not disturbed. For 1D advection equation, the SUPG scheme is expressed as W , uht + auhx Ω = 0, (3.15) where h , i indicates L2 inner product in the Ω. The weighting function W can be written SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 40 as W = W + WSU P G , (3.16) where W uh is the Bubnov-Galerkin weighting function and the Petrov-Galerkin weighting function (WSU P G ) is formulated as WSU P G = τ aWx , (3.17) where τ is the stabilisation parameter. It can be noticed that the Bubnov-Galerkin formulation is recovered for τ = 0. The optimal value of τ for a convection-diffusion equation in 1D, is given as (Donea and Huerta, 2003) h = 2|a| τopt 1 coth|P e| − , Pe (3.18) where P e denotes the Peclet number and h indicates the length of the element. For a pure advection problem, P e → ∞ and equation (3.18) results into τopt = h . 2|a| (3.19) For an implicit time-integration scheme, an alternative expression for τopt is obtained as [15, 21, 22, 76] h τopt = √ . 15|a| (3.20) This maximises the phase accuracy of the semi discrete equation. Therefore, the SUPG formulation of equation (3.12) in terms of τ can be written as Z Ω W uht dΩ + nel Z X e=1 Ωe τ aWx uht + auhx Z dΩ = − Ω aW uhx dΩ. (3.21) 41 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION Applying the divergence theorem on the first term at right hand side (R.H.S.), Z W uht dΩ + Ω nel Z X e=1 Ωe τ aWx uht + auhx Z Z h aWx u dΩ − dΩ = Ω aW uh d∂Ω. (3.22) ∂Ω The finite element approximation for uh and W within a subdomain can be expressed as function of a shape function vector N as uh ≈ ũ = Nj ũj = N T ũ; W ≈ W̃ = Nj W̃j = N T W̃ , where j increments up to 2 for linear elements N = xi+1 −x h x−xi T h (3.23) . Similarly, ux is obtained as uhx ≈ ũx = Bj ũj = BT ũ, (3.24) where B represents the vector field containing the first derivative of the shape functions and T B = − h1 h1 for linear elements. Therefore, the semi-discretised form of equation (3.22) can be expressed as nel X e=1 Z W̃ · T N N ũt dΩ + Ωe − nel X e=1 nel X e=1 Z T W̃ · τ aBN ũt dΩ = Ωe Z W̃ · nel X Z W̃ · e=1 aBN T ũ dΩ Ωe τ a2 BBT ũ dΩ − W̃N auN , (3.25) Ωe where the constant term represents the Neumann boundary condition effect. 3.4 DISCRETISATION IN TIME Once equation (3.12) is discretised in space, a suitable time integration scheme is used to obtain the final solution. Due to the nature of the problems to be solved subsequently, explicit schemes are preferred for implementation. In current section, Forward Euler and Total Variation Diminishing (TVD) Runge-Kutta schemes are employed. 42 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 3.4.1 Forward Euler Method According to this scheme, ut is expressed as n un+1 ∆u ∂u − uni i = + O(∆t) ≈ , ∂t i ∆t ∆t (3.26) where un+1 indicates magnitude of variable at ith node and n + 1 time step. The discretised i form of equation (3.25) can be written as MSU P G ∆ũ + KSU P G ũ = 0, ∆t (3.27) where nel MSU P G = A e=1 nel (MeSU P G ) , In equation (3.28), KSU P G = A e=1 nnode (KeSU P G ) and ũ = A represents the assembly operator, ũi A i=1 (ũi ) . (3.28) = [ui , ui+1 ]T and consistent element mass matrix MeSU P G , constructed using linear shape functions, is expressed as MeSU P G and KeSU P G h 2 1 1 τ a −1 −1 = + , 6 1 2 2 ∆t 1 1 a −1 −1 τ a2 =− + 2 1 h 1 1 −1 . −1 1 (3.29) (3.30) The numerical scheme using SUPG discretisation in space and Forward Euler time integration (SUPG-FE) leads to second order accuracy in space and first order accuracy in time. SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 3.4.2 43 TVD Runge-Kutta Scheme The total variation of a numerical scheme is given as T V (u) = X |ui+1 − ui |. (3.31) i where i indicates the node number. TVD Runge-Kutta time integration scheme [8, 33, 40, 87, 88] ensures T V (un+1 ) ≤ T V (un ). (3.32) where n indicates the time step number. However, it is worth mentioning that a numerical scheme is not entirely TVD unless the spatial derivatives are also approximated by a TVD finite difference or finite element scheme [43, 47]. The semi-discrete advection equation can be recast as M du = L(u) dt (3.33) where M represents the global consistent mass matrix, u represents the global vector of unknowns, L indicates the right side global vector. Rearranging terms du b = M−1 L(u) = L(u) dt (3.34) Two stage second order TVD Runge-Kutta scheme may be written as Step 1: u(1) = un + ∆tLbn , (3.35) b n ). Lbn = L(u (3.36) where 44 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION Step 2: u(2) = u(1) + ∆tLb(1) , (3.37) Lb(1) = Lb u(1) . (3.38) where The solution at (n + 1) time step is given as un+1 = 1 n u + u(2) . 2 (3.39) If observed carefully, the solution un+1 in equation (3.39) can be viewed as the average of solution at time step n and n + 2. The solution at n + 2 time step is obtained using Forward Euler two times consecutively. The TVD Runge-Kutta schme coupled with the SUPG formulation (SUPG-RK) leads to second order accuracy both in space and time. 3.5 CONSISTENCY ANALYSIS The consistency condition dictates that when time and space steps tend to zero, the numerical scheme must tend to the original differential equation. This is one of the requirements imposed on a numerical scheme. In this section, the consistency condition for the second order Taylor-Galerkin scheme and the SUPG-FE scheme are derived. Subsequently, the expression for the stabilisation parameter τ in SUPG scheme is derived by means of a dimensional analysis. 45 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 3.5.1 Consistency Analysis of The Taylor-Galerkin Scheme The resulting discretised equation for an internal node i can be written as a2 ∆t n a h [∆ui−1 + 4∆ui + ∆ui+1 ] = − −uni−1 + uni+1 − −ui−1 + 2uni − uni+1 . 6∆t 2 2h (3.40) Applying Taylor discretisation in time to the left hand side (L.H.S.) results into a2 ∆t n a h n ϑi−1 + 4ϑni + ϑni+1 = − −uni−1 + uni+1 − −ui−1 + 2uni − uni+1 , (3.41) 6∆t 2 2h where ϑ = ∆t (ut ) + ∆t2 ∆t3 (utt ) + (uttt ) . 2 6 (3.42) Grouping terms with same order of derivatives at L.H.S yields a2 ∆t n h ∆t2 a ∆t $tt + $tt = − −uni−1 + uni+1 − −ui−1 + 2uni − uni+1 , $t + 6 2 6 2 2h (3.43) where $ = (u)ni−1 + 4 (u)ni + (u)ni+1 . (3.44) Applying Taylor discretisation in space h (ut )ni h2 ∆t h2 ∆t n n n 2 + (utxx )i + (utt )i + (uttxx )i + O ∆t = 6 2 12 h3 a2 ∆t h2 n n n −a h (ux )i + (uxxx )i − −2 (uxx )i . 6 2h 2 (3.45) 46 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION Converting second order time derivative at L.H.S into second order space derivative using relationship at equation (3.2) ah2 a2 ∆t n a2 ∆t n n n 2 2 (uxxx )i + + (uxx )i . (uxx )i + O h , ∆t = −a (ux )i − 2 6 2 (3.46) (ut + aux )ni + O h2 , ∆t2 = 0. (3.47) (ut )ni Finally Hence, it is proved that the Taylor-Galerkin scheme is second order accurate in time and space. 3.5.2 Consistency Analysis of The SUPG-FE Scheme SUPG space discretisation along with Forward Euler time integration scheme, should result into a solution first order accurate in time and second order accurate in space (SUPG-FE). The discretised equation for an internal node i can be expressed as h τa a [∆ui−1 + 4∆ui + ∆ui+1 ] + [∆ui−1 − ∆ui+1 ] = − −uni−1 + uni+1 6∆t 2∆t 2 τ a2 n − −ui−1 + 2uni − uni+1 . h (3.48) Some of the terms of equation (3.48) is similar to those of equation (3.40) and hence, can be simplified as h (ut )ni a2 ∆t n 2 2 (uxx )i + O h , ∆t + + 2 h3 n n −a h (ux )i + (uxxx )i − 6 τa [∆ui−1 − ∆ui+1 ] = 2∆t τ a2 h2 n −2 (uxx )i . h 2 (3.49) SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 47 The second term on the L.H.S of equation (3.49) can be further manipulated using a Taylor series expansion as a2 ∆t n 2 2 (uxx )i + O h , ∆t h + 2 τa ∆t ∆t n n n n 2 (ut )i−1 + (utt )i−1 − (ut )i+1 − (utt )i+1 + O ∆t + 2 2 2 τ a2 h2 h3 n n n (uxxx )i − −2 (uxx )i . = −a h (ux )i + 6 h 2 (ut )ni (3.50) Applying Taylor discretisation in space for the same term h (ut )ni a2 ∆t n 2 2 + − τ ah (utx )ni + (uttx )ni + O ∆t2 (uxx )i + O h , ∆t 2 h3 τ a2 2 n n = −a h (ux )i + (uxxx )i + h (uxx )ni . (3.51) 6 h Converting the time derivatives into space derivatives leads to h (ut )ni a2 ∆t n n n 2 2 2 + (uxx )i + O h , ∆t − τ a2 h − (u xx )i + (uttx )i + O ∆t 2 h3 n n n 2 ) (uxxx )i + τ a h(u (3.52) = −a h (ux )i + xx i . 6 Finally (ut + aux )ni + O h2 , ∆t = 0. (3.53) Hence, the convergence order holds (i.e. first order in time and second order in space). However, it can be seen from equation (3.52) that the diffusion introduced by BubnovGalerkin scheme still prevails and the stabilisation parameter τ has failed to suppress it. Therefore, this scheme seems to be unstable. An in-depth stability analysis can confirm this speculation. It is also worth noticing that the stabilisation term is similar to Taylor-Galerkin scheme if τ = ∆t . 2 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 3.6 48 VON NEUMANN STABILITY ANALYSIS Once the consistency of the numerical scheme has been established, the following step is to verify its stability. The Von Neumman method [10, 47, 57, 64, 101] can analyse the stability of a scheme in relation to a partial differential equation when considering periodic boundary conditions. The novelty of this stability analysis is the expansion of the error, or the numerical solution, in a finite Fourier series in the spatial frequency. The stability of a scheme is guaranteed on the condition that the amplitude of any harmonic must not increase indefinitely in time. The procedure to perform this stability analysis can be summarised in following steps: n+k I(i+m)ϕ 1. In the numerical scheme, replace all the terms of the form un+k e , i+m with V √ where ϕ represents the phase angle and I indicates −1. It can be noticed that this decomposition separates the time and space dependence of the solution. The time behaviour is denoted by amplitudes V n+k , while the Fourier modes contain the full space dependence. 2. Simplify all the terms by eliminating the factor eIiϕ . 3. Obtain the expression for the amplification factor given as G= V n+1 . Vn (3.54) 4. Ensure the stability is maintained by imposing |G| ≤ 1 for ϕ ∈ [−π, π]. Given the wavelength ζ of a wave pattern, phase angle ϕ can be computed as ϕ= 2πh . ζ (3.55) SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 3.6.1 49 Stability Analysis of The Taylor-Galerkin Scheme The equation for an internal node i can be rewritten from equation (3.40) as a n h n+1 h n ui−1 + 4un+1 ui−1 + 4uni + uni+1 − −ui−1 + uni+1 + un+1 i i+1 = 6∆t 6∆t 2 2 a ∆t (3.56) −uni−1 + 2uni − uni+1 . − 2h Introducing Fourier decomposition h n+1 I(i−1)ϕ V e + 4V n+1 eI(i)ϕ + V n+1 eI(i+1)ϕ = 6∆t a h n I(i−1)ϕ V e + 4V n eI(i)ϕ + V n eI(i+1)ϕ − −V n eI(i−1)ϕ + V n eI(i+1)ϕ 6∆t 2 a2 ∆t −V n eI(i−1)ϕ + 2V n eI(i)ϕ − V n eI(i+1)ϕ . − 2h (3.57) Further simplification yields h I(i−1)ϕ h I(i−1)ϕ e + 4eI(i)ϕ + eI(i+1)ϕ V n+1 = e + 4eI(i)ϕ + eI(i+1)ϕ V n 6∆t 6∆t a I(i−1)ϕ − −e + eI(i+1)ϕ V n 2 a2 ∆t I(i−1)ϕ −e + 2eI(i)ϕ − eI(i+1)ϕ V n . − 2h (3.58) Further manipulation leads to h −Iϕ a −Iϕ h −Iϕ e + 4 + eIϕ V n+1 = e + 4 + eIϕ V n − −e + eIϕ V n 6∆t 6∆t 2 a2 ∆t −Iϕ − −e + 2 − eIϕ V n . (3.59) 2h 50 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION Introducing the expansion eIϕ = cos ϕ + I sin ϕ, e−Iϕ = cos ϕ − I sin ϕ, equation (3.59) can be expressed as h h a2 ∆t [cos ϕ + 2] V n+1 = [cos ϕ + 2] V n − a [I sin ϕ] V n + [cos ϕ − 1] V n . 3∆t 3∆t h (3.60) Hence, the amplification factor (G) can be written as (cos ϕ + 2) − V n+1 = G= n V 3a∆t h (I sin ϕ) + cos ϕ + 2 3a2 ∆t2 h2 (cos ϕ − 1) . (3.61) Introducing the Courant-Friedrichs-Lewy (CFL) number as C= a∆t , h (3.62) equation (3.61) can be expressed as (cos ϕ + 2) − 3C (I sin ϕ) + 3C 2 (cos ϕ − 1) , cos ϕ + 2 (cos ϕ + 2) + 3C 2 (cos ϕ − 1) 3C sin ϕ = −I . cos ϕ + 2 cos ϕ + 2 G= (3.63) In order to satisfy, stability criteria i.e. |G| < 1, in the complex plane, G needs to be contained in the unit circle. Figure 3.1 shows the evolution of G in the complex plane with an increasing CFL number. It is found that the one step Taylor-Galerkin scheme is stable up to a CFL number 0.5. SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 51 Figure 3.1: Polar representation of the amplification factor for the Taylor-Galerkin scheme 3.6.2 Stability Analysis of The SUPG-FE Scheme SUPG-FE scheme is first order in time and second order in space. Considering the particular case where τ= ∆t , 2 (3.64) equation (3.48) can be expressed as a a h [∆ui−1 + 4∆ui + ∆ui+1 ] + [∆ui−1 − ∆ui+1 ] = − −uni−1 + uni+1 6∆t 4 2 a2 ∆t n −ui−1 + 2uni − uni+1 . − 2h (3.65) Separating time steps, a n+1 h n+1 h n n+1 n+1 n n ui−1 + 4un+1 + u + u − u = u + 4u + u i i+1 i i+1 i+1 6∆t 4 i−1 6∆t i−1 a n a + uni−1 − uni+1 − −ui−1 + uni+1 4 2 a2 ∆t n − −ui−1 + 2uni − uni+1 . (3.66) 2h 52 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION Introducing Fourier decomposition, equation (3.66) is simplified as h a h a [2 cos ϕ + 4] V n+1 − [2I sin ϕ] V n+1 = [2 cos ϕ + 4] V n − [2I sin ϕ] V n 6∆t 4 6∆t 4 2 a a ∆t − [2I sin ϕ] V n − [−2 cos ϕ + 2] V n , 2 2h (3.67) Further simplifying, amplification factor (G) is obtained as 2 2 (I sin ϕ) + 3ah∆t (cos ϕ − 1) (cos ϕ + 2) − 9a∆t V n+1 2 2h G= , = 3a∆t n V (cos ϕ + 2) − 2h (I sin ϕ) (cos ϕ + 2) − 9C (I sin ϕ) + 3C 2 (cos ϕ − 1) 2 , (cos ϕ + 2) − 3C (I sin ϕ) 2 3C 2 (cos ϕ + 2) − 9C (I sin ϕ) + 3C (cos ϕ − 1) (cos ϕ + 2) + (I sin ϕ) 2 2 , = 3C (I sin ϕ) (cos ϕ + 2) + (I sin ϕ) (cos ϕ + 2) − 3C 2 2 = = ξ + Iη, (3.68) where ξ= η= (cos ϕ + 2)2 + 3C 2 (cos ϕ + 2) (cos ϕ − 1) + 27 2 C 4 sin2 ϕ (cos ϕ + 2)2 + 49 C 2 sin2 ϕ 3 C 2 , sin ϕ (cos ϕ + 2) + 92 C 2 sin ϕ (cos ϕ − 1) − 29 C sin ϕ (cos ϕ + 2) (cos ϕ + 2)2 + 94 C 2 sin2 ϕ (3.69) . (3.70) Figure 3.2 shows the evolution of G in complex plane with increasing CFL number. It can be clearly seen that SUPG-FE scheme is unstable. This conclusion agrees perfectly with the observations made during consistency analysis in subsection 3.5.2. Efforts have been made in the past to make SUPG-FE scheme stable by neglecting the effect of τ in mass matrix [34]. The stability of the scheme is further improved by lumping the mass matrix. However, the resulting scheme is fairly similar to an one-step Taylor-Galerkin scheme and is not in the interest of current work. 53 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION Figure 3.2: Polar representation of the amplification factor for the SUPG-FE scheme 3.7 SPECTRAL ANALYSIS While performing a Von Neumann stability analysis, the expression of the amplification factor G is obtained as a function of the phase angle ϕ and other parameters of the scheme. Therefore, it can be concluded that G should contain the complete information concerning the behaviour of a numerical scheme and, in particular, the information regarding the generated numerical errors. The error in amplitude of numerical solution is called diffusion or dissipation error. It is defined as the ratio of the computed amplitude to the exact amplitude and can be expressed when particularised for a linear hyperbolic equation as [47] D = |G| (3.71) Accuracy of a numerical scheme requires the modulus of G to be as close to one as possible, but stability requires |G| to be lower than one.The error on the phase of the solution is defined as dispersion or phase error and can be written when particularized for a linear 54 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION hyperbolic equation as [47] ϕ = ϕnum Cϕ (3.72) where −1 ϕnum = tan Im G − Re G (3.73) Dispersion error can be interpreted as the ratio between the numerical and physical convection speeds. When the dispersion error is larger than one, ϕ > 1, the phase error is said to be a leading error and the numerical convection velocity is larger than the exact one. This means that the computed solution moves faster than the physical one. On the contrary, when ϕ < 1, the phase error is a lagging error and the computed solution propagates at a lower velocity than the physical one. 3.7.1 Spectral Analysis of The Taylor-Galerkin Scheme The diffusion error can be calculated using equations (3.63) and (3.71) " D = (cos ϕ + 2) + 3C 2 (cos ϕ − 1) cos ϕ + 2 2 + 3C sin ϕ cos ϕ + 2 2 # 21 . (3.74) The diffusion error is represented in Figure 3.3 in a Cartesian diagram as function of the phase angle ϕ for constant values of the CFL number. It can be observed from Figure 3.3 that the diffusion error is more than one for C = 0.75 at higher phase angle, which renders the numerical scheme unstable. This is quite expected as the stability limit of the Taylor-Galerkin scheme, obtained from Von Neumann stability analysis, is C = 0.5. Using the relationships given in equation (3.72) and equation (3.73), the dispersion error can be 55 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION Figure 3.3: Cartesian representation of the diffusion error as a function of phase angle, in degrees, for Taylor-Galerkin scheme computed as 1 ϕ = tan−1 Cϕ 3C sin ϕ (cos ϕ + 2) + 3C 2 (cos ϕ − 1) . (3.75) Figure 3.4 shows the variation of the dispersion error with increasing phase angle, in a cartesian diagram ϕ for constant values of the CFL number. It is noticed from Figure 3.4 that for lower values values of CFL the numerical solution lags behind the analytical solution. However, for higher magnitude of CFL, the numerical solution leads the analytical solution at lower frequency followed by a sudden phase shift at higher frequency. 3.7.2 Spectral Analysis of SUPG-FE Scheme The diffusion error of the SUPG-FE scheme can be evaluated using the relationship provided in equation (3.68) as D = p ξ 2 + η2. (3.76) SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 56 Figure 3.4: Cartesian representation of the dispersion error as a function of phase angle, in degrees, for Taylor-Galerkin scheme The change in the diffusion error with increasing phase angle is displayed in Figure 3.5. It Figure 3.5: Cartesian representation of the diffusion error as a function of phase angle, in degrees, for SUPG-FE scheme can be observed from Figure 3.5 that the amplitude of diffusive error is less than or equal to one only at extreme ends of the phase angle axis. Hence, it reaffirms that the SUPG-FE scheme is unstable. Similarly, dispersion error for this numerical scheme can be calculated SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 57 as, η 1 −1 − ϕ = tan . Cϕ ξ (3.77) Figure 3.6 represents the dispersion error as a function of phase angle. The numerical Figure 3.6: Cartesian representation of the dispersion error as a function of phase angle, in degrees, for SUPG-FE scheme solution leads the analytical solution at lower phase angle and lags at higher phase angle for low CFL number. However, for higher CFL number, the numerical solution initially lags the analytical solution followed by sudden phase shift. 3.8 NUMERICAL EXAMPLE: COSINE WAVE Hereafter, few numerical examples are discussed in order to demonstrate the accuracy and applicability of each scheme. To begin with, a sufficiently smooth wave pattern, consisting SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 58 of a cosine wave is considered, u(x) = 1 1 + cos π(x−x0 ) , if |x − x0 | ≤ σ. 2 σ 0, (3.78) otherwise. The magnitude of σ and x0 are taken as 0.12 and 0.2, respectively. The domain of the problem is [0, 1]. The speed of the wave is considered as 0.1. Therefore, the time needed for the wave to leave the domain is 9. The propagation of the wave using Bubnov-Galerkin (BG) solution is depicted in Figure 3.7. Number of elements to discretise the domain is 20. As expected, it can be seen that the magnitude of spurious oscillations rises uncontrollably as CFL number increases. The algorithms of various stabilised numerical schemes are Figure 3.7: Propagation of the wave using Bubnov-Galerkin (BG) scheme presented in Figure 3.8. The accuracy of these schemes are compared with the analytical solution. Finally, a detailed convergence study is conducted to verify the order of convergence of each scheme. SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 59 • INPUT number of elements and solution parameters • COMPUTE number of nodes • INITIALIZE vector containing nodal coordinates in domain • COMPUTE number of time steps using equation (3.62) • IF one-step Taylor-Galerkin scheme – COMPUTE global matrices based on equation (3.7) • ELSEIF SUPG scheme – COMPUTE global matrices based on equation (3.25) • INITIALIZE vector containing initial condition (u0 ) based on the wave profile • INCORPORATE Dirichlet boundary conditions • IF time integration scheme TVD Runge-Kutta – ALLOCATE MEMORY for u(1) , u(2) as described in equation (3.35) and equation (3.37) – LOOP over time steps ∗ COMPUTE unknown vector at next time step using equation (3.35)(3.39) – END LOOP • ELSE – LOOP over time steps ∗ COMPUTE unknown vector at next time step using un+1 = un + ∆tLbn – END LOOP • COMPUTE post-processing parameters. Figure 3.8: Algorithm for solving 1D advection equation 3.8.1 Results & Discussions Figure 3.9 shows the propagation of the wave using Taylor-Galerkin (TG) scheme. The value of CFL is taken as 0.5. The phase angles can be computed for different spatial discretisation using the relationship given in equation (3.55). The corresponding phase angles SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION (a) h=0.05, time=3 (b) h=0.05, time=5 (c) h=0.025, time=3 (d) h=0.025, time=5 (e) h=0.0125, time=3 (f) h=0.0125, time=5 60 Figure 3.9: Propagation of a cosine wave using Taylor-Galerkin scheme for the given wave are 45◦ , 22.5◦ and 11.25◦ for element size of 0.05, 0.025 and 0.0125, respectively. Therefore, the reduced amplitude of numerical approximation at Figure 3.9(a) SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 61 and Figure 3.9(b) can be explained through discussion in subsection 3.7.1. Similarly, the lagging in the numerical solution can be explained from an earlier discussion. The overall performance of the TG scheme improves with higher mesh density. The oscillations observed in the numerical solution of BG scheme diminishes effectively when TG scheme is implemented. As wave propagates towards end domain, the numerical solution becomes less diffusive. Figure 3.10 displays movement of the wave using SUPG-FE scheme. The value of CFL is considered as low as 0.1. The increment in diffusion of numerical approximation at Figure 3.10(a) and Figure 3.10(b) can be explained through discussion in subsection 3.7.2. Likewise, the lagging in the numerical solution can be linked to earlier discussion. Although, the numerical approximation improves with lower element size, the spurious oscillations can not be eliminated for the smooth problem. This supports the conclusion reached from the numerical analysis that the SUPG-FE scheme is unstable. Figure 3.11 shows propagation of the wave using a SUPG scheme in space and TVD RungeKutta time integration in time (SUPG-RK scheme). Due to the nature of the time integration scheme and the presence of a consistent mass matrix, it is increasingly difficult to perform analyses such as, Consistency analysis, Von-Neuman stability analysis and Spectral analysis, in an analytical manner. However, the stability condition and consistency of this scheme can be evaluated very accurately through multiple numerical simulations. SUPG-RK scheme displays the same diffusive pattern as the TG scheme. Although, careful observations reveal that the dispersion error in SUPG-RK scheme is slightly smaller than TG scheme. After multiple function evaluations, it is found that SUPG-RK scheme is stable up to CFL number 0.6. This gives a computational advantage over the TG scheme, which is stable only up to C = 0.5. SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION (a) h=0.05, time=3 (b) h=0.05, time=5 (c) h=0.025, time=3 (d) h=0.025, time=5 (e) h=0.0125, time=3 (f) h=0.0125, time=5 Figure 3.10: Propagation of a cosine wave using SUPG-FE scheme 62 SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION (a) h=0.05, time=3 (b) h=0.05, time=5 (c) h=0.025, time=3 (d) h=0.025, time=5 (e) h=0.0125, time=3 (f) h=0.0125, time=5 63 Figure 3.11: Propagation of a cosine wave using SUPG-RK scheme 3.8.2 Convergence Study A convergence study is performed keeping CFL number constant and varying mesh size as well as time step size, accordingly. Number of elements in the domain is taken as 10, SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 64 20, 40 and 80. CFL is considered as 0.5. Error is calculated and normalised with respect to analytical solutions. Convergence of stabilised numerical schemes based on L1 norm is shown in Figure 3.12. The slope is calculated based on least squares approximation. Figure (a) t=1 (b) t=2 (c) t=3 (d) t=4 Figure 3.12: Convergence of stabilised numerical schemes based on L1 norm 3.13 shows the convergence pattern obtained by calculating error based on L2 norm. In order to establish order of convergence, the slope needs to be 1.0 for SUPG-FE and 2.0 for TG and SUPG-RK schemes. The convergence pattern for SUPG-RK and TG is same as expected, in spite of minute differences in the magnitude of error. SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION (a) t=1 (b) t=2 (c) t=3 (d) t=4 65 Figure 3.13: Convergence of stabilised numerical schemes based on L2 norm 3.8.3 Conclusion The propagation of a cosine wave is simulated using different numerical schemes. Clearly, TG and SUPG-RK schemes have better predictability than BG scheme. Their stability and consistency criteria is established. It is observed that SUPG-FE scheme is unstable for solving hyperbolic problems and henceforth, will not be used. 3.9 RIEMANN PROBLEM A Riemann Problem is defined by means of hyperbolic equation together with special initial conditions. The initial data is piecewise constant with a single jump discontinuity at some SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 66 point, say x = 0. u(x) = uL , if x < 0. (3.79) uR , if x > 0. The Riemann Problem plays an important role for developing numerical fluxes for utilising finite volume methods [13, 65, 98]. 3.10 NUMERICAL EXAMPLE: STEP WAVE A step wave propagation problem can be categorised as a Riemann problem. Hence, in this section, the capability of the TG scheme and the SUPG-RK scheme to capture the propagation of a step wave is illustrated. The following wave is considered, u(x) = 1, if |x − x0 | ≤ σ. (3.80) 0, otherwise. The magnitude of σ and x0 are taken as 0.1 and 0.2, respectively. The domain of the problem is [0, 1]. The speed of wave is taken as 0.1. Therefore, the time consumed for the wave to leave the domain is 9. Both TG scheme and SUPG-RK schemes are used to simulate the propagation of this wave. The focus of this example is to investigate the prediction capability of numerical schemes with varying CFL. Number of elements is taken as 160. The accuracy of the approximation is compared against analytical solution at time t = 5. 3.10.1 Results & Discussions Figure 3.14 shows the propagation of the step wave using the TG scheme with CFL numbers varying from 0.2 to 0.5. Whereas, Figure 3.15 displays the movement of the wave using SUPG-RK scheme with different CFL numbers varying from 0.2 to 0.5. Both in Figure 3.14 and Figure 3.15, small wiggles can be noticed near the vicinity of SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION (a) C=0.2 (b) C=0.3 (c) C=0.4 (d) C=0.5 67 Figure 3.14: Propagation of a step wave using Taylor-Galerkin scheme discontinuity. These wiggles are linked to a property known as monotonicity [47]. The monotonicity condition for a numerical scheme can be stated as the requirement that no new extrema be created by the scheme, other than those eventually present in the initial solution. In short, it can be concluded that the numerical oscillations are the consequence of the non-monotone behaviour of the second order schemes. The amplitude of the wiggles in the SUPG-RK scheme is smaller than in the TG scheme. One of the reasons for this phenomenon is the total variation diminishing Runge-Kutta time integration formulation implemented in SUPG-RK scheme. Further observations from Figure 3.14 and Figure 3.15 reveal that the wiggles appears in the upstream end of the discontinuities at lower values of CFL and they move to the downstream SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION (a) C=0.2 (b) C=0.3 (c) C=0.4 (d) C=0.5 68 Figure 3.15: Propagation of a step wave using SUPG-RK scheme end as CFL increases. This can be explained through variation of dispersion error with CFL number. At lower CFL, the magnitude of the dispersion error is predominantly less than one. This means that the numerical solution moves slower than the analytical solution. Therefore, oscillations appear at the upstream side of discontinuity. On the other hand, at higher CFL, the dispersion error is larger than one, which triggers wiggles to appear on the downstream side. SOLUTION OF HYPERBOLIC EQUATION IN ONE-DIMENSION 3.11 69 CONCLUSION This chapter has presented three finite element formulations for solving the linear advection equation. One of them, SUPG-FE is found to be unstable.The remaining two numerical schemes, SUPG-RK and TG perform well and provide second order convergence both in space and time. Moreover, SUPG-RK scheme shows higher stability and smaller wiggles when compared with the TG solution. A constant value of stabilisation parameter τ for SUPG formulation is adopted throughout this chapter. However, in later part of this work the effect of different values of τ will be discussed. Chapter 4 SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 70 SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 4.1 71 INTRODUCTION The general form of a hyperbolic system of linear equations is given as U t + AI ∂U I = 0, ∂XI ∀I = 1, 2, 3 in Ω, (4.1) with boundary and initial conditions U = Ů in (Ω, t0 ), U = UD on (∂ΩD , t), Un = UN on (∂ΩN , t), (4.2) where U denotes the vector containing the unknown conservation variables, U t indicates the partial derivative of U with respect to time t, t0 indicates the initial time and Ω is a bounded open domain in <n , n = 1, 2, 3 with Lipschitz continuous boundary ∂Ω = ∂ΩD ∪ ∂ΩN . AI is so called flux-Jacobian matrix in the “I” spatial direction and AI NI possesses real eigenvalues ∀N = [N1 , N2 , N3 ] /||N || = 1. For 1D problems with 2 degrees of freedom (D.O.F), equation (4.1) can be written as U t + AU x = 0, (4.3) where u f a U= , A= . v b g (4.4) In equation (4.4), u, v denote the unknown scalar variables. A simplified expression for A with physical relevance is obtained by setting and by setting f, g = 0. In this chapter, the analytical approach for solving equation (4.3) is discussed. Next, the basic formulation for obtaining a numerical solution will be presented. Both analytical and SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 72 numerical approaches are illustrated with a numerical example of the propagation of an acoustic wave in a shock tube. This discussion is followed by the introduction of a finite element formulation in a conservative framework. Finally, concluding remarks are made regarding the suitability of different frameworks. 4.2 ANALYTICAL SOLUTION Looking at the nature of the hyperbolic equations presented in equations (4.3) and (4.4), the solution is expected to be composed of two waves propagating to the left and right of the domain. This suggests for a solution of the form U (x, t) = U (x − st) (4.5) for wave speed s, where U (η) is some function of a variable η. With this assumption, it can be seen that U t (x, t) = −sU 0 (x − st) , U x (x, t) = U 0 (x − st) . (4.6) Therefore, equation (4.3) reduces to AU 0 (x − st) = sU 0 (x − st) . (4.7) As s is a scalar and A is a matrix, this is only plausible if s is an eigenvalue of A and U 0 (η) is its corresponding right eigenvector of A for each value of η. For the matrix A given in equation (4.4) the eigenvalues are computed as Λ1 = −c0 and Λ2 = c0 , (4.8) 73 SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS where √ c0 = ab, (4.9) which represents the speed of the wave. As expected, waves propagate in both directions. The right eigenvectors of the flux-jacobian matrix A are, R1 = − c 0 b , R2 = 1 c0 b . (4.10) 1 Hence the general solution can be expressed as, U = w1 (x + c0 t) R1 + w2 (x − c0 t) R2 , (4.11) where w1 (η) , w2 (η) are scalar functions of some variable η also known as characteristic variables. The values of these functions will depend on the initial data given U (x, t0 ) = Ů (x) = ů (x) . (4.12) v̊ (x) For simplicity, t0 is considered as zero. Therefore, after implementing the initial condition equation (4.11) reduces to Ů (x) = w1 (x) R1 + w2 (x) R2 . (4.13) After solving equation (4.17), the values of the characteristic variables are obtained as b w = 2c0 b w2 = 2c0 1 h i c0 −ů (x) + v̊ (x) , b h i c0 ů (x) + v̊ (x) . b (4.14) SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 74 Therefore, the solution of equation (4.3) becomes 1 c0 [ů (x + c0 t) + ů (x − c0 t)] − [v̊ (x + c0 t) − v̊ (x − c0 t)] , 2 2b 1 b [ů (x + c0 t) − ů (x − c0 t)] + [v̊ (x + c0 t) + v̊ (x − c0 t)] . v(x, t) = − 2c0 2 u(x, t) = 4.2.1 (4.15) Numerical Example: Propagation of Acoustic Wave An acoustic wave is a small pressure disturbance which travels through a compressible gas resulting infinitesmall changes in density and pressure. The linearised equations modelling this phenomenon can be written as p u0 K0 p + = 0, 1 u u u 0 ρ0 t (4.16) x where p and u represent pressure and velocity distribution in the gas, respectively. K0 indicates the bulk modulus of compressibility, ρ0 denotes the initial density of the gas and u0 represents the velocity of the gas. A simplified set of equations are obtained for motionless state by setting u0 = 0, yielding pt + K0 ux = 0, ρ0 ut + px = 0. (4.17) Following the procedure discussed earlier, the final solution is obtained as [65] 1 Z0 [p̊ (x + c0 t) + p̊ (x − c0 t)] − [ů (x + c0 t) − ů (x − c0 t)] , 2 2 1 1 u(x, t) = − [p̊ (x + c0 t) − p̊ (x − c0 t)] + [ů (x + c0 t) + ů (x − c0 t)] , 2Z0 2 p(x, t) = (4.18) SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 75 where s Z0 = ρ0 c0 and c0 = K0 . ρ0 (4.19) Figure 4.1 shows the evolution of an initial pressure perturbation, concentrated near the origin, into distinct simple waves propagating with velocities −c0 and c0 and corresponding change in velocities for the following set of initial conditions [65] 1 2 p̊ = e−80x + S (x) , 2 ů = 0, where S (x) = 1, if −0.3 < x < −0.1. (4.20) (4.21) 0, otherwise. The magnitude of ρ0 and K0 is taken as 1 and 0.25, respectively. 4.3 NUMERICAL SOLUTION In this section, different numerical approaches to obtain the solution of equation (4.3) is discussed. 4.3.1 Taylor-Galerkin Scheme Equation (4.3) can be expressed as U t = −AU x , (4.22) U tt = (−AU x )t = −AU tx = A (AU x )x . (4.23) which leads to 76 SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS Velocity at time t=0.5 1.5 1 1 0.5 0.5 0 0 u p Pressure at time t=0.5 1.5 −0.5 −0.5 −1 −1 −1.5 −1.5 −1 −0.5 0 x 0.5 1 −1.5 −1.5 1.5 −1 −0.5 (a) 1 1 0.5 0.5 0 0 −0.5 −0.5 −1 −1 −0.5 0 x 1 1.5 1 1.5 Velocity at time t=1 1.5 u p Pressure at time t=1 −1 0.5 (b) 1.5 −1.5 −1.5 0 x 0.5 1 1.5 −1.5 −1.5 −1 (c) −0.5 0 x 0.5 (d) Figure 4.1: Evolution of pressure and velocity waves Applying Taylor discretisation in time, an expression similar to equation (3.4) is obtained as ∆U U n+1 − U n ∆t := = −AU nx + A (AU x )nx + O ∆t2 . ∆t ∆t 2 (4.24) The weak form of the above strong form can be written as ∆U h dΩ + W· ∆t Ω Z Z W· Ω AU hx Z dΩ = Ω ∆t W · A AU hx x dΩ. 2 (4.25) SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 77 After applying divergence theorem to the R.H.S., the weak formulation leads to ∆U h dΩ + W· ∆t Ω Z Z W· AU hx Ω Z ∆t Wx · AAU hx dΩ 2 ZΩ ∆t W · AAU hx d∂Ω, + ∂Ω 2 dΩ = − (4.26) hereafter, referred as “form A”. Further application of the divergence theorem to the second term of the L.H.S. results in ∆U h W· dΩ = ∆t Ω Z Z Z ∆t Wx · AU dΩ − Wx · AAU hx dΩ Ω Ω 2 Z Z ∆t h W · AU d∂Ω + W · AAU hx d∂Ω, − ∂Ω ∂Ω 2 h (4.27) henceforth, referred as “form B”. Following a standard finite element expansion, form A given in equation (4.26) leads to nel X Z nel X ∆Ũ dΩ = − W̃ · NN W̃ · N ABT Ũ dΩ ∆t e e Ω Ω e=1 e=1 Z Z n el X ∆t ∆t T BAAB Ũ dΩ + WN · AAU x d∂Ω. − W̃ · ∂ΩN 2 Ωe 2 e=1 Z T (4.28) Analogously, form B as expressed in equation (4.27) can be expressed as nel X Z nel X ∆Ũ W̃ · NN dΩ = W̃ · BAN T Ũ dΩ ∆t e e Ω Ω e=1 e=1 Z n el X ∆t − W̃ · BAABT Ũ dΩ Ωe 2 Ze=1 Z ∆t − WN · AU d∂Ω + WN · AAU x d∂Ω. ∂ΩN ∂ΩN 2 Z T (4.29) SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 78 For linear shape functions, x(i+1) −x h N = x−xi h 0 0 1 0 − h 0 1 0 0 h , B = . x(i+1) −x 0 −1 h h x−xi 1 0 h h (4.30) The discretised form of equation (4.29) can be written as MT G ∆Ũ + KT G Ũ = 0, ∆t (4.31) where nel MT G = A nel e=1 (MeT G ) , KT G = A nnode e=1 (KeT G ) and Ũ = A i=1 Ũi . (4.32) In equation (4.44), A represents the assembly operator, Ũi = [ui , ui+1 , vi , vi+1 ]T and consistent element mass matrix MeT G is expressed as MeT G and KeT G 0 1 0 =− 2 −b −b 2 h 1 = 6 0 0 1 0 2 0 0 2 0 1 0 0 , 1 2 (4.33) 0 −a a 0 1 −1 0 0 −a a −1 1 0 0 ab∆t + . 2h 0 b 0 0 0 1 −1 b 0 0 0 0 −1 1 (4.34) 79 SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 4.3.2 SUPG-RK Scheme The general form of SUPG scheme is expressed as h(W + WSU P G ), R U h iΩ = 0, (4.35) ∂U h h = U t + AI −S , ∂XI (4.36) where R U h ∀I = 1, 2, 3, and ∂W . ∂XI (4.37) WSU P G = τ AT Wx . (4.38) WSU P G = τ ATI For one dimensional problems The SUPG formulation of equation (4.3) can then be written as Z W· U ht dΩ + Ω nel Z X T τ A Wx · Ωe e=1 U ht + AU hx Z dΩ = − W · AU hx dΩ, (4.39) Ω to be referred to as form A. The application of divergence theorem at the R.H.S. yields Z W· Ω U ht dΩ + nel Z X e=1 Ωe T τ A Wx · U ht + AU hx Z Wx · AU h dΩ dΩ = Ω Z − ∂Ω W · AU h d∂Ω, (4.40) 80 SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS henceforth known as form B. After FE expansion, form A given in equation (4.39) can be expressed as Z nel X ∆Ũ ∆Ũ W̃ · τ BAN T dΩ + dΩ = W̃ · NN ∆t ∆t Ωe Ωe e=1 e=1 Z Z nel nel X X T τ BAABT Ũ dΩ. W̃ · W̃ · N AB Ũ dΩ − − nel X Z T Ωe e=1 (4.41) Ωe e=1 Whereas, form B of given in equation (4.40) leads to nel X Z nel X ∆Ũ ∆Ũ NN W̃ · τ BAN T W̃ · dΩ + dΩ = ∆t ∆t e e Ω Ω e=1 e=1 Z Z nel nel X X T BAN Ũ dΩ − W̃ · W̃ · Z T Ωe e=1 τ BAABT Ũ dΩ Ωe e=1 Z − WN · AU d∂Ω. (4.42) ∂ΩN Equation (4.42) can be rewritten as M dU = L(U) dt (4.43) where nel M= nnode A (M ) , e e=1 and Ũ = A i=1 Ũi . (4.44) where A represents the assembly operator, Ũi = [ui , ui+1 , vi , vi+1 ]T and consistent element mass matrix MeT G represents element consistent mass matrix. Performing simple operation, dU b = M−1 L(U) = L(U) dt (4.45) A two step second order TVD Runge-Kutta scheme can be sketched as Step 1: U(1) = Un + ∆tLbn , (4.46) SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 81 where b n ). Lbn = L(U (4.47) U(2) = U(1) + ∆tLb(1) , (4.48) Lb(1) = Lb U(1) . (4.49) Step 2: where The solution at (n + 1)th time step is given as Un+1 = 4.3.3 1 Un + U(2) . 2 (4.50) Numerical Example: Propagation of Acoustic Wave The acoustic wave propagation problem discussed in section 4.2.1 is solved using form A of both the one-step Taylor-Galerkin scheme and the SUPG-RK shecme. Homogeneous dirichlet boundary conditions are considered on the left boundary. Homogeneous Neumann boundary conditions are considered on the right boundary. Analyses are carried out with 300, 600 and 1200 elements. CFL is chosen as 0.5. Figure 4.2 displays the propagation of the pressure and velocity waves using the Taylor-Galerkin scheme at t = 0.5. The propagation of the acoustic wave is captured using the SUPG-RK scheme in Figure 4.3. The stabilisation parameter τ is chosen as τ = τ I, τ= ∆t 2 (4.51) It can be observed from Figure 4.2 and Figure 4.3, that for both the schemes numerical solutions converge to corresponding analytical solutions quickly. The wiggles in the solution appear as a result of the non-monotonicity. SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 82 Pressure at time t=0.5 using TG scheme 1.5 Analytical solution 300 elements 600 elements 1200 elements p 1 0.5 0 −0.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 x (a) Pressure Velocity at time t=0.5 using TG scheme 2 Analytical solution 300 elements 600 elements 1200 elements 1.5 1 u 0.5 0 −0.5 −1 −1.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 x (b) Velocity Figure 4.2: Evolution of pressure and velocity waves using Taylor-Galerkin Scheme at t = 0.5, C = 0.5 4.4 CONSERVATIVE FRAMEWORK In fast dynamics problems, shock discontinuities along with geometric and material nonlinearities are often encountered. Therefore, it is convenient to express the first order hyperbolic equation in equation (4.1) in a conservative form as Ut + ∂F I = 0 ∀I = 1, 2, 3, ∂XI (4.52) SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 83 Pressure at time t=0.5 using SUPG−RK scheme 1.5 Analytical solution 300 elements 600 elements 1200 elements p 1 0.5 0 −0.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 x (a) Velocity at time t=0.5 using SUPG−RK scheme 2 Analytical solution 300 elements 600 elements 1200 elements 1.5 1 u 0.5 0 −0.5 −1 −1.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 x (b) Figure 4.3: Evolution of pressure and velocity waves using SUPG-RK Scheme at t = 0.5, C = 0.5 where F I represents the Flux component in I direction and Flux Jacobian matrix AI is defined as AI = ∂F I . ∂U (4.53) It can be noticed that if AI 6= f (U ) i.e. the Flux Jacobian matrix is linear ∂F I ∂U = AI , ∂XI ∂XI (4.54) SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 84 both conservative and non-conservative framework of governing equations render the same solution. 4.4.1 Taylor-Galerkin Scheme in Conservative Form Equation (4.52) for 1D problems can be rewritten as U t = −F x , (4.55) Further simplifications lead to U tt = −F xt = −F tx = − ∂F Ut ∂U = (AF x )x (4.56) x Taylor discretisation in time leads to U n+1 − U n ∆t ∆U := = −F nx + (AF x )nx + O ∆t2 . ∆t ∆t 2 (4.57) The weak form for the above equation can be written as ∆U h W· dΩ + ∆t Ω Z Z W· F hx Ω Z dΩ = Ω ∆t W · AF hx x dΩ. 2 (4.58) Application of divergence theorem on the R.H.S results in form A as ∆U h W· dΩ + ∆t Ω Z Z W· Ω F hx Z dΩ = − Ω ∆t Wx · AF hx dΩ + 2 Z ∂Ω ∆t W · AF hx d∂Ω. 2 (4.59) 85 SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS Further application of the divergence theorem on the second term of the L.H.S. leads to form B as ∆U h dΩ = W· ∆t Ω Z Z Z h Wx · F dΩ − Ω Ω ∆t Wx · AF hx dΩ − 2 Z W · F h d∂Ω Z∂Ω + ∂Ω ∆t W · AF hx d∂Ω. 2 (4.60) After FE expansion, the equation (4.59) can be expressed as nel X Z nel X ∆Ũ W̃ · NN W̃ · N BT F̃ dΩ dΩ = − ∆t e e Ω Ω e=1 e=1 Z Z n el X ∆t ∆t − W̃ · BABT F̃ dΩ + WN · AF x d∂Ω, Ωe 2 ∂ΩN 2 e=1 Z T (4.61) and form B is presented as nel X n el X ∆Ũ W̃ · NN dΩ = W̃ · ∆t e Ω e=1 e=1 Z T − nel X Z Z W̃ · Ωe Ze=1 + ∂ΩN BN T F̃ dΩ Ωe ∆t BABT F̃ dΩ − 2 Z WN · F d∂Ω ∂ΩN ∆t WN · AF x d∂Ω. 2 (4.62) For linear shape functions, the expression in equation (4.61) can be simplified when particularized for an intermediate element as be written as MT G ∆Ũ + KT G Ũ = 0, ∆t (4.63) where nel nel MT G = A e=1 (MeT G ) , KT G = A e=1 nnode (KeT G ) and Ũ = A i=1 Ũi . (4.64) SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 86 In equation (4.64), A indicates the assembly operator, Ũi = [ui , ui+1 , vi , vi+1 ]T and consistent element mass matrix MeT G is expressed as MeT G and KeT G 2 h 1 = 6 0 0 0 0 , 1 2 1 0 2 0 0 2 0 1 (4.65) 0 1 −1 a −a 0 0 0 √ 0 −1 1 ab∆t 1 a −a 0 0 0 = + 2 0 0 b −b 2h 1 −1 0 0 0 0 b −b −1 1 0 0 (4.66) Similar expression can be computed for equation (4.62). 4.4.2 SUPG-RK Scheme in Conservative Form The SUPG formulation in conservative framework for form A equation can be written as Z W· U ht dΩ + Ω nel Z X e=1 T τ A Wx · Ωe U ht + F hx Z W · F hx dΩ. dΩ = − (4.67) Ω Using the divergence theorem on the R.H.S. form B can be expressed as Z W· Ω U ht dΩ+ nel Z X e=1 Ωe T τ A Wx · U ht + F hx Z h Z Wx · F dΩ− dΩ = Ω W · F h d∂Ω. ∂Ω (4.68) 87 SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS The spatial semi-discretised form for A leads to Z nel X ∆Ũ τ BAN T Ũ t dΩ = W̃ · W̃ · NN dΩ + ∆t e e Ω Ω e=1 e=1 Z Z nel n el X X T W̃ · W̃ · N B F̃ dΩ − − nel X Z T Ωe e=1 e=1 τ BABT F̃ dΩ, Ωe (4.69) and spatial semi-discretised form for B can be expressed as nel X Z nel X ∆Ũ ∆Ũ W̃ · NN τ BAN T W̃ · dΩ + dΩ = ∆t ∆t Ωe Ωe e=1 e=1 Z Z nel nel X X T BN F̃ dΩ − W̃ · W̃ · Z T Ωe e=1 e=1 τ BABT F̃ dΩ Ωe Z − WN · F d∂Ω. (4.70) ∂ΩN Final solution can be obtained by implementing TVD-RK scheme as sketched in subsection 4.3.2 (equations (4.69) and (4.70)). 4.5 SELECTION OF STABILISATION PARAMETER The stabilisation parameter τ proposed in original SUPG formulation is expressed as [22, 50] τ = τI (4.71) Here, two separate expressions for τ is employed based on spatial discretisation and temporal discretisation. They are given as follows: • Spatial criteria: τ = Zχh , R where Z indicates a non-dimensional parameter, χ is a parameter related to stability and R is the spectral radius of Flux-jacobian matrix AI . • Temporal criteria: τ = Zχ∆t. SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 88 It is worth pointing out that, in the original formulation [22] τ is proposed based only on spatial criteria. This is highly logical for an implicit scheme as τ needs to be independent of the size of time step ∆t. However, for an explicit scheme, h and ∆t are related through CFL number and therefore, τ can be defined based on either criteria. Considerable work has been carried out [22, 50] to observe the effect on stability and accuracy of the solution by considering Zχ = 12 , 14 , √115 . This definition of τ yields lower values for higher order elements [49, 94]. Subsequent minor modifications of τ are carried out considering the interaction between a shock-capturing term and the earlier version of τ [91]. Lately, new formulations of τ are proposed for higher order elements [39]. One of these method suggests calculation of τ based on element level matrices and vectors [95]. These definitions are generally written in terms of the ratios of the norms of the matrices or vectors e.g., for an advection dominated problem τ is expressed as τ= ∆t kCe k , 2 C̃e (4.72) where ∂FhI W· C : dΩ ∂XI Ωe e Z e Z and C̃ : Ωe ∂W ∂Uh · dΩ. ∂XI ∂t (4.73) These definitions automatically take into account the local length scales and advection field. Moreover, based on these definitions τ can be calculated for each element and each degree of freedom [26]. It is also well known that the stabilization parameters to be used in advancing the solution from time step n to n + 1 should be evaluated at time step n [93, 95]. In this way, further nonlinearity in solution scheme is avoided. In current work, τ is defined from a temporal perspective. Generally, τ = analyses using SUPG-RK scheme, unless otherwise stated. ∆t 2 for all the SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 4.5.1 89 Numerical Experiment: Propagation of Acoustic Wave An experiment is conducted in order to observe the effect of varying τ . Figure 4.5 depicts the evolution of pressure and velocity waves for the acoustic wave propagation problem varying τ at t = 0.5. Number of elements taken for the study is 600 and the magnitude of CFL is Pressure at time t=0.5 using SUPG−RK scheme Pressure at time t=0.5 using SUPG−RK scheme 1 Analytical solution τ=dt/2 τ=dt/3 τ=dt/4 p p 0.5 0 −0.5 −0.8 0.5 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0 x x (a) Pressure (b) Magnified view Velocity at time t=0.5 using SUPG−RK scheme Velocity at time t=0.5 using SUPG−RK scheme 2 Analytical solution τ=dt/2 τ=dt/3 τ=dt/4 1.5 1 u u 0.5 0 1 −0.5 −1 −1.5 −0.8 −0.6 −0.4 −0.2 0 x (c) Velocity 0.2 0.4 0.6 0 x (d) Magnified view Figure 4.4: Evolution of pressure and velocity waves using SUPG-RK scheme with varying τ at C = 0.5 considered as 0.5. Further investigations reveal that the solution can not be obtained beyond the value of τ = 0.5∆t when CFL is 0.5. However, by reducing CFL at 0.3, solution can be reached when τ = ∆t. Figure 4.4 displays the evolution of pressure and velocity waves using the SUPG-RK scheme with varying τ at t = 0.5 and C = 0.3. It can be observed both from Figure 4.4 and Figure 4.5, increasing magnitude of τ adds more diffusion to the solution and therefore, minimises the spurious oscillations. However, after a certain value SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS Pressure at time t=0.5 using SUPG−RK scheme 90 Pressure at time t=0.5 using SUPG−RK scheme 1 Analytical solution τ=dt/1 τ=dt/2 τ=dt/4 p p 0.5 0.5 0 −0.5 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0 x x (a) Pressure (b) Magnified view Velocity at time t=0.5 using SUPG−RK scheme Velocity at time t=0.5 using SUPG−RK scheme 2 Analytical solution τ=dt/1 τ=dt/2 τ=dt/4 1.5 1 u u 0.5 0 1 −0.5 −1 −1.5 −0.8 −0.6 −0.4 −0.2 0 x (c) Velocity 0.2 0.4 0.6 0 x (d) Magnified view Figure 4.5: Evolution of pressure and velocity waves using SUPG-RK scheme with varying τ at C = 0.3 of τ , the solution becomes over diffusive and renders unstable. Therefore, for a given CFL there exists an optimum value of τ . 4.6 CONCLUSION This chapter introduces the stabilised finite element formulations for solving 1D hyperbolic problems with 2 D.O.F. The numerical formulations for the one step Taylor-Galerkin scheme and the SUPG-RK scheme are presented both in a non-conservative and conservative framework. It is concluded that for a conservative framework is more suitable for solving fast dynamics problems. Moreover, the accuracy of the stabilised schemes is demon- SOLUTION OF A 1D TWO-EQUATION SYSTEM OF HYPERBOLIC EQUATIONS 91 strated through a numerical example of propagation of an acoustic wave in a shock-tube. Furthermore, two possible ways to solve the discretised equations are presented and they are indicated as “form A” and “form B”. An in depth investigation on the influence of these forms will be presented in next chapter. Finally, a discussion is included on selection of the stabilisation parameter τ for the SUPG-RK scheme. The effect on the solution by varying τ is demonstrated. It is observed that τ plays an important role in determining the accuracy of the solution and the optimum magnitude of τ is dependant on the CFL number. Chapter 5 1D FAST STRUCTURAL DYNAMICS PROBLEM 92 93 1D FAST STRUCTURAL DYNAMICS PROBLEM 5.1 INTRODUCTION This chapter addresses 1D fast structural dynamics problems. The nature of this type of problems is similar to the problems discussed in previous chapter. Therefore, the conservative finite element framework can be implemented to solve these problems. In the following sections, the governing equations for 1D fast dynamics problems is derived, their eigenstructure is discussed. Extensive numerical analyses are carried out in order to address the accuracy, convergence as well as suitability of form A and form B. Furthermore, the accuracy of the stabilised finite element schemes is compared with the finite volume schemes. Afterwards, a shock capturing scheme is introduced in order to improve the accuracy of the stabilised finite element schemes. The performance of the shock capturing scheme is compared with slope limiters implemented in finite volume schemes. Finally, some conclusions are reached regarding the suitability of numerical schemes. 5.2 GOVERNING EQUATION For a reversible process, the mixed formulation for one dimensional fast dynamics reduces to U t + F 1x = 0, (5.1) where U= p1 F11 , F1 = −P11 . (5.2) − p1 ρ0 In equation (5.2), p1 denotes linear momentum, F11 represents the deformation gradient and ρ0 indicates constant material density. In the small strain linear regime, the First Piola- 94 1D FAST STRUCTURAL DYNAMICS PROBLEM Kirchhoff stress is expressed as, P11 = (λ + 2µ)(F11 − 1), (5.3) where λ and µ are Lamé constants. Therefore, in a non-conservative form, the governing equations can be written as U t + A1 U x = 0, (5.4) where 0 −(λ + 2µ) A1 = . − ρ10 0 (5.5) Furthermore, imposition of a null Poisson’s ratio, ν = 0, leads to λ= νE = 0, and E = 2µ(1 + ν) = 2µ, (1 + ν)(1 − 2ν) (5.6) which results in a Flux Jacobian matrix A1 as 0 −E A1 = . 1 − ρ0 0 (5.7) The eigenvalues of A1 are computed as s Λ1 = − E and Λ2 = ρ0 s E , ρ0 (5.8) For a linear problem, time is computed based on the maximum eigenvalue as ∆t = Ch . Λ2 (5.9) 95 1D FAST STRUCTURAL DYNAMICS PROBLEM The algorithm to solve 1D solid mechanics problems is presented in Figure 5.1. • INPUT number of elements and solution parameters • COMPUTE number of nodes and nodal coordinates in domain [0, 10] • DEFINE material properties • COMPUTE number of time steps using equation (5.9) • ALLOCATE MEMORY for p1 , F11 and X • IF one-step Taylor-Galerkin scheme – COMPUTE global matrices based on equation (4.61) (form A) or (4.62) (form B) • ELSEIF SUPG-RK scheme – COMPUTE global matrices based on equation (4.69) or (4.70) (form B) • COMPUTE global matrices • IF time integration scheme TVD Runge-Kutta – IMPLEMENT time integration algorithm as presented in Figure 5.2 • ELSE – IMPLEMENT time stepping scheme as presented in Figure 5.3 • COMPUTE post-processing parameters. Figure 5.1: Algorithm for solving 1D solid mechanics problem 5.3 BAR SUBJECTED TO SINUSOIDAL TRACTION In order to demonstrate the capabilities of the one step Taylor-Galerkin and SUPG-RK schemes, a sinusoidal traction force is applied axially to a bar displayed in Figure 5.4 with properties as shown in Table 5.1. The bar is fixed at one end and the tensile traction loading applied to the free end is expressed as T = T0 sin(ωt) (5.10) 96 1D FAST STRUCTURAL DYNAMICS PROBLEM • ALLOCATE MEMORY for U(1) , U(2) as referred in equation (4.46) and equation (4.48) • LOOP over time steps – EVALUATE flux F n1 using equation (5.2) and (5.3) – CORRECT flux using Neumann boundary conditions (only for form A) R – EVALUATE boundary term ∂Ω W F (only for form B) – COMPUTE U(1) using equation (4.46) – IMPLEMENT Dirichlet boundary conditions and CORRECT F11 at nodes of known traction – COMPUTE F (1) and U(2) – CORRECT flux using Neumann boundary conditions (only for form A) R – EVALUATE boundary term ∂Ω W F (only for form B) – IMPLEMENT Dirichlet boundary conditions and CORRECT F11 at nodes of known traction – COMPUTE U(2) using equation (4.48) – COMPUTE X using trapezoidal rule • END LOOP Figure 5.2: Algorithm for TVD Runge-Kutta time integration scheme • LOOP over time steps – EVALUATE flux F n1 using equation (5.2) and (5.3) – CORRECT flux using Neumann boundary conditions (only for form A) R – EVALUATE boundary term ∂Ω W F (only for form B) R – COMPUTE ∂Ω ∆t W AF x d∂Ω 2 – COMPUTE Un+1 = Un + ∆tLn – IMPLEMENT Dirichlet boundary conditions and CORRECT F11 at nodes of known traction – COMPUTE X using trapezoidal rule • END LOOP Figure 5.3: Algorithm for Taylor-Galerkin time integration scheme where T0 = 1 × 10−3 N and ω = 0.1. (5.11) 97 1D FAST STRUCTURAL DYNAMICS PROBLEM y T=T0 f(t) x L Figure 5.4: Axially-loaded cantilever bar Table 5.1: Properties of axially-loaded bar Property Magnitude Linear density (m) 1 kg/m Young’s modulus (E) 1 N/m2 Poison’s ratio (ν) 0 Length (L) 10 m Cross section area (A) 1 m2 The objective of this study is to verify the accuracy the numerical schemes provide by comparing the numerical solutions with corresponding analytical solutions. The analytical solution is obtained through Fourier analysis. The displacement at mid-point of the bar is given as ∞ X 2T0 L X(x, t) = (−1) mL (ωn2 − ω 2 ) n=1 n (2n − 1)πx ω sin(ωn t) sin sin(ωt) − ωn 2L (5.12) where π ωn = (2n − 1) 2 r EA . mL2 (5.13) The domain is discretized into 50, 100 and 200 elements. CFL number is considered as 0.5. Figure 5.5 shows the time history of displacement, velocity and axial stress at the mid-span of the bar using form A of Taylor-Galerkin scheme as presented in equation (4.61). The re- 98 1D FAST STRUCTURAL DYNAMICS PROBLEM sults are compared against analytical solution provided in equation (5.11). The time history Solution at mid−span using TG scheme Solution at mid−span using TG scheme 0.03 Analytical solution 50 elements 100 elements 200 elements 0.025 Displacement (m) Displacement (m) 0.02 0.015 0.01 0.005 0 −0.005 −0.01 −0.015 −0.02 0 20 40 60 80 Time (s) 100 Time (s) (a) Displacement Solution at mid−span using TG scheme −3 3 x 10 Solution at mid−span using TG scheme Analytical solution 50 elements 100 elements 200 elements Velocity (m/s) 2 Velocity (m/s) (b) Magnified view 1 0 −1 −2 0 20 40 60 80 Time (s) 100 Time (s) (c) Velocity Solution at mid−span using TG scheme −3 5 (d) Magnified view x 10 4 2 Axial stress (N/m2) 3 Axial stress (N/m ) Solution at mid−span using TG scheme Analytical solution 50 elements 100 elements 200 elements 2 1 0 −1 −2 −3 0 20 40 60 80 100 Time (s) Time (s) (e) Axial stress (f) Magnified view Figure 5.5: Time history of a bar under sinusoidal traction using TG scheme (form A) of the same variables, calculated based on form B given in equation (4.62), is depicted in Figure 5.6. 99 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at mid−span using TG scheme Solution at mid−span using TG scheme 0.03 Analytical solution 50 elements 100 elements 200 elements 0.025 Displacement (m) Displacement (m) 0.02 0.015 0.01 0.005 0 −0.005 −0.01 −0.015 −0.02 0 20 40 60 80 Time (s) 100 Time (s) (a) Displacement x 10 Analytical solution 50 elements 100 elements 200 elements Velocity (m/s) 2 Velocity (m/s) Solution at mid−span using TG scheme Solution at mid−span using TG scheme −3 3 (b) Magnified view 1 0 −1 −2 0 20 40 60 80 Time (s) 100 Time (s) (c) Velocity x 10 Analytical solution 50 elements 100 elements 200 elements 4 2 2 Axial stress (N/m ) 3 Axial stress (N/m ) Solution at mid−span using TG scheme Solution at mid−span using TG scheme −3 5 (d) Magnified view 2 1 0 −1 −2 −3 0 20 40 60 80 100 Time (s) Time (s) (e) Axial stress (f) Magnified view Figure 5.6: Time history of a bar under sinusoidal traction using TG scheme (form B) Similarly, time history of displacement, velocity and axial stress at mid-span of the bar obtained using form A and form B of SUPG-RK scheme is plotted in Figure 5.7 and Figure 5.8, respectively. 100 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at mid−span using SUPG−RK scheme Solution at mid−span using SUPG−RK scheme 0.03 Analytical solution 50 elements 100 elements 200 elements 0.025 Displacement (m) Displacement (m) 0.02 0.015 0.01 0.005 0 −0.005 −0.01 −0.015 −0.02 0 20 40 60 80 Time (s) 100 Time (s) (a) Displacement x 10 Analytical solution 50 elements 100 elements 200 elements Velocity (m/s) 2 Velocity (m/s) Solution at mid−span using SUPG−RK scheme Solution at mid−span using SUPG−RK scheme −3 3 (b) Magnified view 1 0 −1 −2 0 20 40 60 80 Time (s) 100 Time (s) (c) Velocity Solution at mid−span using SUPG−RK scheme −3 5 (d) Magnified view x 10 Analytical solution 50 elements 100 elements 200 elements 4 2 Axial stress (N/m2) 3 Axial stress (N/m ) Solution at mid−span using SUPG−RK scheme 2 1 0 −1 −2 −3 0 20 40 60 80 100 Time (s) Time (s) (e) Axial stress (f) Magnified view Figure 5.7: Time history of a bar under sinusoidal traction using SUPG-RK scheme (form A) From the results, it seems that form A of both stabilised schemes provide slightly better solution than form B for stress and velocity. However, the solution for displacement is almost identical as it is evaluated through time integration of velocity. A trapezoidal time 101 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at mid−span using SUPG−RK scheme Solution at mid−span using SUPG−RK scheme 0.03 Analytical solution 50 elements 100 elements 200 elements 0.025 Displacement (m) Displacement (m) 0.02 0.015 0.01 0.005 0 −0.005 −0.01 −0.015 −0.02 0 20 40 60 80 Time (s) 100 Time (s) (a) Displacement x 10 Analytical solution 50 elements 100 elements 200 elements Velocity (m/s) 2 Velocity (m/s) Solution at mid−span using SUPG−RK scheme Solution at mid−span using SUPG−RK scheme −3 3 (b) Magnified view 1 0 −1 −2 0 20 40 60 80 Time (s) 100 Time (s) (c) Velocity Solution at mid−span using SUPG−RK scheme −3 5 (d) Magnified view x 10 Analytical solution 50 elements 100 elements 200 elements 4 Axial stress (N/m2) 3 Axial stress (N/m2) Solution at mid−span using SUPG−RK scheme 2 1 0 −1 −2 −3 0 20 40 60 80 100 Time (s) Time (s) (e) Axial stress (f) Magnified view Figure 5.8: Time history of a bar under sinusoidal traction using SUPG-RK scheme (form B) integration scheme is implemented to obtain displacement. 102 1D FAST STRUCTURAL DYNAMICS PROBLEM 5.4 CONVERGENCE STUDY The purpose of this numerical experiment is to verify the order convergence of the proposed stabilised finite element schemes. In order to obtain the convergence plot, the bar at section 5.3 is subjected to an exponential loading as it is C ∞ continuous. The properties and boundary conditions are considered same as before. The tensile traction loading applied at free end is expressed as 2 T = T0 e−ω(t−13) (5.14) T0 = 1 × 10−3 N and ω = 0.1. (5.15) where The convergence study for is conducted using 20, 40, 80, 160 and 320 elements at C = 0.5. Normalised error is calculated based on L2 norm at t = 10 s. Figure 5.9 depicts the convergence of displacement, velocity and stress for the Taylor-Galerkin formulation. Figure 5.10 shows the convergence of the same parameters for the SUPG-RK formulation. TG scheme 0 10 −1 −1 10 −2 L2 Norm Error L2 Norm Error 10 10 −3 10 Displacement (form A) Velocity (form A) Stress (form A) First−order Second−order −4 10 −5 10 TG scheme 0 10 −2 10 −1 10 Element Size (a) form A −2 10 −3 10 Displacement (form B) Velocity (form B) Stress (form B) First−order Second−order −4 10 −5 0 10 10 −2 10 −1 10 Element Size 0 10 (b) form B Figure 5.9: Convergence of the Taylor-Galerkin scheme based on L2 norm for 1D fast dynamics problem 103 1D FAST STRUCTURAL DYNAMICS PROBLEM The error is calculated based on interpolation of the nodal variables in a coarser mesh to SUPG−RK scheme −1 10 −2 −2 10 L2 Norm Error L2 Norm Error 10 −3 10 −4 −3 10 −4 10 10 Displacement (form A) Velocity (form A) Stress (form A) First−order Second−order −5 10 SUPG−RK scheme −1 10 −2 10 −1 10 Element Size Displacement (form B) Velocity (form B) Stress (form B) First−order Second−order −5 0 10 10 (a) form A −2 10 −1 10 Element Size 0 10 (b) form B Figure 5.10: Convergence of SUPG-RK scheme based on L2 norm for 1D fast dynamics problem nodal coordinates of very fine mesh. In current work, a domain of 640 elements is taken for computation of normalised error. It can be seen from Figure 5.9 and Figure 5.10, both schemes provide second order accuracy for displacement, velocity and stress. However, in a relatively coarser mesh the convergence of stress is faster in SUPG-RK scheme than Taylor-Galerkin scheme. The convergence pattern for form A is identical to that of form B. 5.5 BAR SUBJECTED TO CONSTANT TRACTION A step-function loading is applied to the same bar. The purpose for this experiment is to observe the ability of both form A and form B to capture the shock waves. The properties and boundary conditions are considered same as before. The tensile traction loading applied at free end is expressed as T = T0 (5.16) 104 1D FAST STRUCTURAL DYNAMICS PROBLEM where T0 = 1 × 10−3 N. (5.17) The displacement of the bar obtained using Fourier analysis is expressed as X(x, t) = ∞ X n+1 (−1) n=1 2T0 1 − cos(ωn t) sin mL ωn2 (2n − 1)πx 2L (5.18) where ωn is given in equation (5.13). The domain is discretized into 50, 100 and 200 elements. CFL number is considered as 0.5. Figure 5.11 and Figure 5.13 display the time history of displacement, velocity at free end and axial stress at fixed end of the bar using form A and form B of Taylor-Galerkin scheme, respectively. Variation of the same parameters is simulated in Figure 5.12 for form A and in Figure 5.14 for form B of the SUPG-RK scheme. 5.5.1 Discussion From Figure 5.5, Figure 5.6, Figure 5.7 and Figure 5.8, it can be noticed that the numerical solutions converge quickly to corresponding analytical results. The spurious oscillations arise for axial stress and velocity. Moreover, if observed carefully, it can be seen the accuracy of the form A is slightly higher than form B for both Taylor-Galerkin and SUPG-RK schemes. The reason is that the jump in flux in intermediate elements is neglected in form B. Nevertheless, it is still advantageous to implement form B for two reasons. Firstly, the underlying principle behind form A is based on the fact that flux can be linearly interpolated in terms of shape functions i.e. F = Ni Fi . However, this is true only in small-strain framework. Secondly, the Neumann boundary conditions appear naturally in the form B as a result of applying divergence theorem, hence, easier to implement. Therefore, form B is adopted in later part of this work. The predication of displacement is fairly similar of all the schemes and not of primary concern. 105 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at free end using TG scheme 0.03 Analytical solution 50 elements 100 elements 200 elements Displacement (m) Displacement (m) 0.025 0.02 0.015 0.01 0.02 0.005 0 0 20 40 60 80 60 Time (s) 100 Time (s) (a) Displacement −3 Solution at free end using TG scheme −3 x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 3 2.5 2 1.5 1 Velocity (m/s) Velocity (m/s) (b) Magnified view 1.5 1 0.5 0 0.5 0 −0.5 −0.5 −1 −1.5 −1 −2 0 20 40 60 80 60 100 Time (s) Time (s) (c) Velocity x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 4 3.5 3 2.5 2 2 Axial stress (N/m ) Axial stress (N/m2) −3 Solution at fixed end using TG scheme −3 4.5 (d) Magnified view 2.5 2 1.5 1 0.5 1.5 1 0.5 0 −0.5 −1 0 0 20 40 60 Time (s) (e) Axial stress 80 100 60 Time (s) (f) Magnified view Figure 5.11: Time history of a bar under constant traction using TG scheme (form A) 106 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at free end using SUPG−RK scheme 0.03 Analytical solution 50 elements 100 elements 200 elements Displacement (m) Displacement (m) 0.025 0.02 0.015 0.01 0.02 0.005 0 0 20 40 60 80 60 Time (s) 100 Time (s) (a) Displacement −3 Solution at free end using SUPG−RK scheme −3 x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 3 2.5 2 1.5 1 Velocity (m/s) Velocity (m/s) (b) Magnified view 1.5 1 0.5 0 0.5 0 −0.5 −0.5 −1 −1.5 −1 −2 0 20 40 60 80 60 100 Time (s) Time (s) (c) Velocity x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 4 3.5 3 2.5 Axial stress (N/m2) Axial stress (N/m2) −3 Solution at fixed end using SUPG−RK scheme −3 4.5 (d) Magnified view 2.5 2 1.5 1 0.5 2 1.5 1 0.5 0 −0.5 −1 0 0 20 40 60 80 100 Time (s) (e) Axial stress 60 Time (s) (f) Magnified view Figure 5.12: Time history of a bar under constant traction using SUPG-RK scheme (form A) 5.5.2 Effect of Varying Stabilisation Parameter A numerical experiment is conducted in order to observe the change in solution with varying stabilisation parameter τ . For this purpose, CFL number is reduced to 0.4 and magnitude 107 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at free end using TG scheme 0.03 Analytical solution 50 elements 100 elements 200 elements Displacement (m) Displacement (m) 0.025 0.02 0.015 0.01 0.02 0.005 0 0 20 40 60 80 60 Time (s) 100 Time (s) (a) Displacement (b) Magnified view x 10 Analytical solution 50 elements 100 elements 200 elements 3 2.5 2 1.5 1 Velocity (m/s) Velocity (m/s) −3 Solution at free end using TG scheme −3 x 10 1.5 1 0.5 0 0.5 0 −0.5 −0.5 −1 −1.5 −1 −2 0 20 40 60 80 60 100 Time (s) Time (s) (c) Velocity x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 4 3.5 3 2.5 Axial stress (N/m2) Axial stress (N/m2) −3 Solution at fixed end using TG scheme −3 4.5 (d) Magnified view 2.5 2 1.5 1 2 1.5 1 0.5 0.5 0 −0.5 −1 0 0 20 40 60 80 100 60 Time (s) Time (s) (e) Axial stress (f) Magnified view Figure 5.13: Time history of a bar under constant traction using TG scheme (form B) of τ is considered as ∆t, ∆t and 2 ∆t . 4 Number of elements is taken as 100. The result is shown in Figure 5.15. As expected, as τ increases more diffusion is added to the solution until it reaches an optimum value. When CFL number is equal to 0.4, ∆t is found to be the 108 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at free end using SUPG−RK scheme 0.03 Analytical solution 50 elements 100 elements 200 elements Displacement (m) Displacement (m) 0.025 0.02 0.015 0.01 0.02 0.005 0 0 20 40 60 80 60 Time (s) 100 Time (s) (a) Displacement (b) Magnified view −3 Analytical solution 50 elements 100 elements 200 elements 3 2.5 1.5 1 Velocity (m/s) 2 Velocity (m/s) x 10 Solution at free end using SUPG−RK scheme −3 x 10 1.5 1 0.5 0 0.5 0 −0.5 −0.5 −1 −1.5 −1 −2 0 20 40 60 80 60 100 Time (s) Time (s) (c) Velocity (d) Magnified view −3 x 10 Analytical solution 50 elements 100 elements 200 elements 4 3.5 3 2.5 Axial stress (N/m2) Axial stress (N/m2) x 10 Solution at fixed end using SUPG−RK scheme −3 4.5 2.5 2 1.5 1 0.5 2 1.5 1 0.5 0 −0.5 −1 0 0 20 40 60 Time (s) (e) Axial stress 80 100 60 Time (s) (f) Magnified view Figure 5.14: Time history of a bar under constant traction using SUPG-RK scheme (form B) optimum value of τ . However, at stability limit (C = 0.5), the optimum magnitude of τ is known to be ∆t . 2 109 1D FAST STRUCTURAL DYNAMICS PROBLEM −3 x 10 x 10 Analytical solution τ=dt/1 τ=dt/2 τ=dt/4 3 2.5 1.5 1 Velocity (m/s) 2 Velocity (m/s) −3 Solution at free end using SUPG−RK scheme 1.5 1 0.5 0 0.5 0 −0.5 −0.5 −1 −1.5 −1 −2 0 20 Time (s) 20 40 Time (s) (a) Velocity x 10 x 10 Analytical solution τ=dt/1 τ=dt/2 τ=dt/4 4 3.5 Axial stress (N/m2) −3 Solution at fixed end using SUPG−RK scheme 3 2.5 2 Axial stress (N/m2) −3 4.5 (b) Magnified view 2.5 2 1.5 1 0.5 1.5 1 0.5 0 0 −0.5 −1 0 20 Time (s) (c) Axial stress Time (s) 40 (d) Magnified view Figure 5.15: Effect of varying τ at C = 0.4 5.5.3 Comparison of Different Schemes A study is conducted to compare the accuracy of stabilised finite element schemes implemented in current work with cell-centered finite volume schemes [61, 63, 62]. The domain is discretised into 100 elements and analyses are carried out at C = 0.5. Figure 5.16 shows the solutions obtained through one step Taylor-Galerkin, SUPG, first and second finite volume schemes. The solutions are plotted against analytical solution. If observed carefully, it can be seen that the magnitude of wiggles in Taylor-Gaerkin scheme is higher than that of SUPG-RK scheme. This occurs as larger diffusion is added into the solution due to TVD Runge-Kutta time integration scheme associated with SUPG-RK scheme. Further observations reveal that monotonicity of first order finite volume scheme ensures a solution without 110 1D FAST STRUCTURAL DYNAMICS PROBLEM x 10 x 10 SUPG−RK Taylor Galerkin FV first order FV second order Analytical solution 2.5 2 Solution at free end 1.5 1 Velocity (m/s) 1.5 Velocity (m/s) −3 Solution at free end −3 3 1 0.5 0 0.5 0 −0.5 −0.5 −1 −1.5 −1 −2 0 20 Time (s) 20 Time (s) 40 (a) Velocity x 10 Axial stress (N/m2) 3 2.5 x 10 SUPG−RK Taylor Galerkin FV first order FV second order Analytical solution Solution at fixed end 2.5 2 Axial stress (N/m2) 3.5 −3 Solution at fixed end −3 4 (b) Magnified view 2 1.5 1 0.5 1.5 1 0.5 0 −0.5 0 −1 0 20 Time (s) (c) Axial stress Time (s) 40 (d) Magnified view Figure 5.16: Comparison of numerical schemes wiggles. However, the non-monotonic second order finite volume scheme solution includes spurious oscillations. If noticed carefully, the magnitude of the wiggles in second order finite volume method is less than finite element schemes. This is mainly because the flux is evaluated at the boundary of each element in finite volume scheme. Whereas, in finite element scheme flux is computed only at the boundary of the domain making it less diffusive. Additionally, higher number of D.O.F are present in the finite element scheme as compared to finite volume method. 111 1D FAST STRUCTURAL DYNAMICS PROBLEM 5.6 SHOCK CAPTURING SCHEME The purpose of introducing the shock capturing scheme is to remove spurious oscillations from the solution. This can be achieved usually done by adding diffusion near the vicinity of shock. The generalised semidiscretised form of SUPG-RK scheme in multi-dimension can be expressed as Z W ·R U h dΩ + Ω n nel Z X Ωe e=1 τ ATJ el X ∂W · R U h dΩ + ∂XJ e=1 Z δ Ωe ∂W ∂U h · dΩ = 0, ∂XI ∂XI (5.19) where R U h = ∂U h ∂F hI + ∂t ∂XI I, J = 1, 2, 3. (5.20) The semi discrete form when particularised for 1D fast structural dynamics problem leads to Z W· U ht dΩ + Ω + nel Z X e e=1 Ω nel Z X e=1 τ AT Wx · U ht + F hx dΩ δWx · Ωe U hx Z h Z W · F h d∂Ω. (5.21) Wx · F dΩ − dΩ = Ω ∂Ω Similarly, the semi discrete form for TG scheme with shock capturing can be written as n el X ∆U h W· dΩ + ∆t Ω e=1 Z Z δWx · Ωe U hx Z h Z ∆t Wx · AF hx dΩ 2 Wx · F dΩ − dΩ = Ω Z − Ω h Z W · F d∂Ω + ∂Ω ∂Ω ∆t W · AF hx d∂Ω. 2 (5.22) An earlier version of the shock capturing parameter δ is proposed by Le Beau and Tezduyar [59, 60]. The expression for this shock capturing term, used in the context of conservation variables, is derived from the shock capturing term expressed in terms of entropy variables [52, 53]. Moreover, computation of this term involves evaluation of the Jacobian matrix of 1D FAST STRUCTURAL DYNAMICS PROBLEM 112 the transformation from the entropy variables to the conservation variables [60]. Recently, a novel shock capturing scheme, known as Y Zβ shock capturing, has been proposed [96, 97]. This scheme has been widely used for capturing shocks in flow problems [11, 26, 27, 77, 78]. In addition to being much simpler than the earlier shock capturing scheme, Y Zβ Shock capturing scheme yields substantially better shock quality [96, 97]. 5.6.1 Y Zβ Shock Capturing Y Zβ shock capturing are formulated based on scaled residuals and are defined with options for smoother or sharper shocks. “YZ” version of Y Zβ shock capturing scheme can be expressed as [83, 96, 97] !β/2−1 β nsd −1 X −1 ∂U h 2 hshoc , δ= Y Z Y ∂Xi 2 (5.23) i=1 where Z= ∂U h ∂F hI + ∂t ∂XI hshoc = 2 nen X , Y = diag (U ref ) , !−1 |j · ∇Na | = h (for 1D problems). (5.24) a=1 In equation (5.24), Y is a diagonal scaling matrix containing reference values of components of U and j is defined in direction of density for flow problems. “YZU ” version of Y Zβ shock capturing scheme can be written as !β/2−1 nsd −1 X −1 ∂U h 2 −1 h 1−β hshoc β δ= Y Z Y U . Y ∂Xi 2 (5.25) i=1 In flow problems, different expressions for scaling matrix Y has been tested. One of these expressions recommends inflow condition for variables in Eulerian configuration [92], which corresponds to initial values of variables at Lagrangian configuration. Since, the sec- 1D FAST STRUCTURAL DYNAMICS PROBLEM 113 ond Piola-Kirchhoff stress is zero initially, this expression can not adopted. Therefore, the expression for Y is considered as Y = αI, (5.26) where α is a scalar parameter which needs to be tuned for a specific problem. Lowering the magnitude of α enables better shock capturing. However, there exists an optimum value beyond which α can not be reduced. In equations (5.23) and (5.25) β = 1 for smoother shocks and β = 2 for sharper shocks. In current work, the magnitude of β is taken as 2. 5.6.2 Numerical Experiments A series of numerical experiments are carried out for determining the efficiency of Y Zβ shock capturing scheme. To begin with, the elastic bar subjected to constant traction force, discussed in section , is analysed. The value of CFL is taken as 0.3 and α = 0.0006. Number elements to discretise the domain is taken as 50, 100 and 200. Figure 5.17 shows the variation in axial stress and velocity using Taylor-Galerkin scheme. The change in same parameters, using SUPG-RK scheme, is displayed in Figure 5.18. If observed carefully, Y Zβ shock capturing scheme performs slightly better along with SUPG-RK scheme. This can again be related to two stage TVD RK time integration. An investigation is conducted to determine the convergence trend for a sufficiently smooth problem in presence of shock capturing scheme. For this purpose, the bar subjected to exponential traction force (as discussed in section 5.4 is analysed. Figure 5.19 displays the convergence pattern for Taylor-Galerkin and SUPG-RK schemes. It can be observed that the velocity, displacement and axial stress of the bar are second order convergent in both schemes in spite of inclusion of shock capturing term in the stabilised finite element formulation. A further comparison is made between Y Zβ shock capturing scheme implemented in sta- 114 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at free end using TG scheme −3 x 10 Analytical solution 50 elements 100 elements 200 elements 3 2.5 2 1 0.5 1.5 Velocity (m/s) Velocity (m/s) −3 x 10 1 0.5 0 0 −0.5 −0.5 −1 −1.5 −1 −2 0 20 40 60 80 60 100 Time (s) Time (s) (a) Velocity Solution at fixed end using TG scheme −3 4.5 (b) Magnified view x 10 Analytical solution 50 elements 100 elements 200 elements 4 3.5 3 2 Axial stress (N/m2) Axial stress (N/m2) −3 x 10 2.5 2 1.5 1 0.5 1.5 1 0.5 0 −0.5 −1 0 0 20 40 60 80 100 Time (s) Time (s) (c) Axial stress (d) Magnified view Figure 5.17: Time history of a bar under constant traction using TG scheme with shock capturing bilised finite element formulations and slope limiters [61] applied in second order cell centred finite volume scheme (Figure 5.20). The solution of the SUPG-RK scheme is slightly better than the Taylor Galerkin scheme due to the inclusion of TVD RK scheme. The slope limiter seems to work well when implemented in more diffusive cell centred finite volume scheme. The solution of the finite volume and stabilised finite element schemes is also compared with Trapezoidal Newmark scheme [20]. The spurious oscillations in the Newmark scheme near discontinuity appear as expected. 115 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at free end using SUPG−RK scheme −3 x 10 Analytical solution 50 elements 100 elements 200 elements 3 2.5 2 1 0.5 1.5 Velocity (m/s) Velocity (m/s) −3 x 10 1 0.5 0 0 −0.5 −0.5 −1 −1.5 −1 −2 0 20 40 60 80 100 60 Time (s) Time (s) (a) Velocity x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 4 3.5 3 2 Axial stress (N/m2) Axial stress (N/m2) −3 Solution at fixed end using SUPG−RK scheme −3 4.5 (b) Magnified view 2.5 2 1.5 1 0.5 1.5 1 0.5 0 −0.5 −1 0 0 20 40 60 80 Time (s) 100 Time (s) (c) Axial stress (d) Magnified view Figure 5.18: Time history of a bar under constant traction using SUPG-RK scheme with shock capturing 5.7 NON-LINEAR ELASTICITY An investigation is carried out for observing the change in the solution when the constitutive model is changed from linear elastic to nearly incompressible Neo-Hookean material (discussed in subsection ). The imposition of null Poisson’s ratio leads to calculation of Lamé constants and shear modulus κ as λ = 0, µ= E 2 and κ = 2µ . 3 (5.27) The magnitude of E is considered as unity. The amplitude of the loading is taken as 0.1 i.e. 100 times more than elastic analysis in order to activate the non linearity of the constitutive 116 1D FAST STRUCTURAL DYNAMICS PROBLEM TG scheme 0 SUPG−RK scheme −1 10 10 −1 10 −2 −2 L2 Norm Error L2 Norm Error 10 10 −3 10 −3 10 −4 −5 10 10 Displacement Velocity Stress First−order Second−order −4 10 −2 10 −1 −5 0 10 Element Size Displacement Velocity Stress First−order Second−order 10 10 −2 −1 10 10 Element Size (a) 0 10 (b) Figure 5.19: Convergence of stabilised finite element schemes in the presence of shock capturing term model. The other parameters are kept same as in elastic analysis. Figure 5.21 shows the variation in displacement, stress and velocity at mid point of a bar while being subjected to a sinusoidal traction force using Taylor-Galerkin scheme. The variation in the same parameters using the SUPG-RK formulation is displayed in Figure 5.22. The elastic solution is obtained through Fourier analysis. The non linear behaviour can very well be captured using both schemes. It is worth noticing that the size of the time step in a non-liear analysis is not constant, but changes as the solution evolves in time. The steps for calculation of time step size are summarised hereafter. • For every node, compute parameters as 2 5 γ1 = κJ 2 + µJ − 3 (F : F ) , 9 2 γ2 = µJ − 3 2 2 and γ3 = − µJ − 3 3 (5.28) • Calculate longitudinal wave speed as s Up = γ2 + γ1 Λ2 ρ0 + 2γ3 , where Λ is eigenvalue related to F , ρ0 is taken as unity. (5.29) 117 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at free end −3 3 x 10 SUPG−RK Taylor Galerkin FV second order Newmark solution Analytical solution 2.5 2 Solution at free end 1.5 1 Velocity (m/s) 1.5 Velocity (m/s) −3 x 10 1 0.5 0 0.5 0 −0.5 −0.5 −1 −1.5 −1 −2 0 20 Time (s) 20 Time (s) 40 (a) Velocity Solution at fixed end −3 x 10 3.5 Axial stress (N/m2) 3 2.5 −3 x 10 Solution at fixed end 2.5 SUPG−RK Taylor Galerkin FV second order Newmark Solution Analytical solution 2 Axial stress (N/m2) 4 (b) Magnified view 2 1.5 1 0.5 1.5 1 0.5 0 0 −0.5 −1 0 20 Time (s) 40 (c) Axial stress Time (s) (d) Magnified view Figure 5.20: Comparison of shock capturing schemes with Newmark scheme and Slope limiter • Finally time step size can be computed as ∆t = 5.8 Ch max(Up ) (5.30) LUMPED MASS FORMULATION One of the concerns of current study is the range of stability of the stabilised finite element schemes. Since, adoption of consistent mass matrix in the finite element formulation leads to a stability limit of C = 0.5 for both the schemes. One of the ways to increase this limit is by lumping the mass matrix [82, 105]. Figure 5.23 shows the time history of a bar under 118 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at mid−span using TG scheme Solution at mid−span using TG scheme 2.5 Linear elastic solution 50 elements 100 elements 200 elements 2 Displacement (m) Displacement (m) 1.5 1 0.5 0 1 −0.5 −1 −1.5 0 20 40 60 80 Time (s) 100 Time (s) (a) Displacement (b) Magnified view Solution at mid−span using TG scheme Solution at mid−span using TG scheme 0.3 Linear elastic solution 50 elements 100 elements 200 elements Velocity (m/s) Velocity (m/s) 0.2 0.1 0 −0.1 −0.1 −0.2 0 20 40 60 80 Time (s) 100 Time (s) (c) Velocity (d) Magnified view Solution at mid−span using TG scheme Solution at mid−span using TG scheme 0.5 Linear elastic solution 50 elements 100 elements 200 elements 0.4 Axial stress (N/m2) Axial stress (N/m2) 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 0 20 40 60 80 100 Time (s) Time (s) (e) Axial stress (f) Magnified view Figure 5.21: Time history of a bar under sinusoidal traction using TG scheme (Neo-Hookean material) constant traction using the Taylor-Galerkin scheme with lumped mass matrix. Similarly, the time history obtained through SUPG-RK scheme is displayed in Figure 5.24. The CFL 119 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at mid−span using SUPG−RK scheme Solution at mid−span using SUPG−RK scheme 2.5 Linear elastic solution 50 elements 100 elements 200 elements 2 Displacement (m) Displacement (m) 1.5 1 0.5 0 1 −0.5 −1 −1.5 0 20 40 60 80 Time (s) 100 Time (s) (a) Displacement (b) Magnified view Solution at mid−span using SUPG−RK scheme Solution at mid−span using SUPG−RK scheme 0.3 Linear elastic solution 50 elements 100 elements 200 elements Velocity (m/s) Velocity (m/s) 0.2 0.1 0 −0.1 −0.1 −0.2 0 20 40 60 80 Time (s) 100 Time (s) (c) Velocity (d) Magnified view Solution at mid−span using SUPG−RK scheme Solution at mid−span using SUPG−RK scheme 0.5 Linear elastic solution 50 elements 100 elements 200 elements 0.4 Axial stress (N/m2) Axial stress (N/m2) 0.3 0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 0 20 40 60 80 100 Time (s) Time (s) (e) Axial stress (f) Magnified view Figure 5.22: Time history of a bar under sinusoidal traction using SUPG-RK scheme (NeoHookean material) number for these analyses is extended to 1.0 which is twice the stability limit obtained from consistent mass matrix. Although, adopting lumped mass matrix offers further advantages 120 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at free end using TG scheme 0.03 Analytical solution 50 elements 100 elements 200 elements Displacement (m) Displacement (m) 0.025 0.02 0.015 0.01 0.02 0.005 0 0 20 40 60 80 60 Time (s) 100 Time (s) (a) Displacement −3 Solution at free end using TG scheme −3 x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 3 2.5 1 0.5 Velocity (m/s) 2 Velocity (m/s) (b) Magnified view 1.5 1 0.5 0 0 −0.5 −0.5 −1 −1 −1.5 −1.5 −2 0 20 40 60 80 60 100 Time (s) Time (s) (c) Velocity x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 4 3.5 3 20 Axial stress (N/m2) Axial stress (N/m2) −4 Solution at fixed end using TG scheme −3 4.5 (d) Magnified view 2.5 2 1.5 1 15 10 5 0.5 0 0 −0.5 −1 0 −5 20 40 60 80 100 Time (s) Time (s) (e) Axial stress (f) Magnified view Figure 5.23: Time history of a bar under constant traction using TG scheme and lumped mass formulation such as faster convergence and smaller spurious oscillations, it also raises dispersion (or phase) error [105] in the solution. However, this phase error can be reduced by increasing 121 1D FAST STRUCTURAL DYNAMICS PROBLEM Solution at free end using SUPG−RK scheme 0.03 Analytical solution 50 elements 100 elements 200 elements Displacement (m) Displacement (m) 0.025 0.02 0.015 0.01 0.02 0.005 0 0 20 40 60 80 60 Time (s) 100 Time (s) (a) Displacement −3 Solution at free end using SUPG−RK scheme −3 x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 3 2.5 1 0.5 Velocity (m/s) 2 Velocity (m/s) (b) Magnified view 1.5 1 0.5 0 0 −0.5 −0.5 −1 −1 −1.5 −1.5 −2 0 20 40 60 80 60 100 Time (s) Time (s) (c) Velocity x 10 x 10 Analytical solution 50 elements 100 elements 200 elements 4 3.5 3 20 Axial stress (N/m2) Axial stress (N/m2) −4 Solution at fixed end using SUPG−RK scheme −3 4.5 (d) Magnified view 2.5 2 1.5 1 15 10 5 0.5 0 0 −0.5 −1 0 −5 20 40 60 80 100 Time (s) Time (s) (e) Axial stress (f) Magnified view Figure 5.24: Time history of a bar under constant traction using SUPG-RK scheme and lumped mass formulation the mesh density. Since, one of the objectives of the proposed scheme is to obtain solution with moderately 1D FAST STRUCTURAL DYNAMICS PROBLEM 122 coarse mesh, consistent mass matrix has been adopted throughout this work. The stability limit can also be increased by using a multi-stage time stepping scheme at the expense of higher computational cost. 5.9 CONCLUSION This chapter presents the implementation of the stabilised finite element formulation for solving 1D fast dynamics problems. The capabilities of these schemes have been demonstrated through several numerical experiments. It is found that the stresses in the solution converge at same rate as displacements and velocities. Although, form A provides slightly better result than form B, it is found the implementation of form B is easier and therefore, adopted in the rest of the work. The finite element framework also provides opportunities to implement different constitutive models. The SUPG-RK scheme performs better than the one-step Taylor-Galerkin scheme. The solution obtained through finite element schemes are compared to that of finite volume schemes. It is found that spurious oscillations in finite element scheme and second order finite volume schemes arise in presence of discontinuity due to their non-monotonic nature. A shock capturing scheme developed in Eulerian finite element formulation is implemented in the Lagrangian formulation in order to alleviate those wiggles and their efficiency is compared to slope limiters applied in second order finite volume formulation. The effect of lumped mass matrix on the solution is also investigated. Although, lumped mass matrix drastically increases the stability limit, it is not adopted to reduce dispersion error in the solution. Chapter 6 2D FAST STRUCTURAL DYNAMICS PROBLEM 123 2D FAST STRUCTURAL DYNAMICS PROBLEM 6.1 124 INTRODUCTION This chapter briefly introduces the two-dimensional (2D) implementation of fast dynamics problems. It addresses the curl-fee involution arising from compatibility condition. Moreover, the numerical results obtained from preliminary analysis is presented. 6.2 CURL-FREE INVOLUTION In this SUPG scheme, the deformation gradient F is a primary variable which contrasts with the usual finite element formulation where the spatial coordinate x is the primary variable. For convenience, the expression for F is recalled here as F = ∂x ∂φ (X, t) = = ∇0 x, ∂X ∂X (6.1) ∂xi . ∂XI (6.2) or following index notation FiI = The gradient of deformation gradient tensor satisfies ∂ 2 xi ∂ 2 xi = . ∂XI ∂XJ ∂XJ ∂XI (6.3) In the case, where F is primary variable, the compatibility condition must be satisfied as curlF = 0, (6.4) in order to guarantee the existence of a single-valued continuous displacement field. These conditions are met under exact integration provided they are satisfied by the initial condition, which in turn implies that the curl preservation is an inherent analytical property of the evolution operator. A robust numerical scheme must have the ability to control curlerrors under long-term response analysis. These errors usually accumulate and lead to a 2D FAST STRUCTURAL DYNAMICS PROBLEM 125 breakdown of the scheme. An initial effort to ensure compatibility is made by adopting the discretised “form B” for the conservation of linear momentum equation and “form A” for conservation of deformation gradient tensor. This seems logical ans the curl-error arises from the second equation and “form A” ensures that the order of continuity of v is not altered. This is useful to recall that “form A” does not invoke the Divergence theorem unlike “form B”. Additionally, adopting “form B” for the conservation of linear momentum allows the traction force in the boundary to appear naturally as a result of employing the Divergence theorem. 6.3 6.3.1 NUMERICAL EXAMPLES Short Column In order to evaluate the performance of the proposed scheme in large strain bending dominated problems, the vibration of a short column is considered. The column is being subjected to a linear distributed initial velocity V0 = 0.1 Y /L m/s. The properties of the column is listed in Table 6.1. The magnitude of CFL is taken as 0.2 and value of stabilisation paTable 6.1: Properties of short column Property Magnitude Density (ρ0 ) 1.1e03 kg/m3 Young’s modulus (E) 1.7e07 N/m2 Poison’s ratio (ν) 0.45 Length (L) 6m Width (b) 1m rameter is considered as τ = 0.2∆t. The domain is discretised in to 4 × 24 elements.Figure 6.1 displays the evolution of pressure in the column at different time instant. 2D FAST STRUCTURAL DYNAMICS PROBLEM 126 Figure 6.1: Evolution of pressure distribution at different time instant for a short column 6.3.2 Tensile Test For further investigating, a large deformation of punching test is considered. A square plate of size 1 m × 1 m is constrained at the bottom, left and right boundary as shown in Figure 6.2. The plate is pulled uniformly with an initial velocity V0 = 500 m/s. Nearly incompressible constitutive material is implemented for the numerical computation. The 127 2D FAST STRUCTURAL DYNAMICS PROBLEM V0=500 m/s Y L X L Figure 6.2: A tensile test properties of the material is show in Table 6.2. The domain is discretised in to 10 × 10 Table 6.2: Properties of plate Property Density (ρ0 ) Young’s modulus (E) Poison’s ratio (ν) Length (L) Width (b) Magnitude 7e03 kg/m3 2.1e07 N/m2 0.3 1m 1m elements. The magnitude of CFL is kept at 0.2 and value of stabilisation parameter is considered as τ = 0.2∆t. Figure 6.3 displays pressure distribution in the plate. If observed carefully, the solution leads to singularity in the solution at t = 0.0016131 s. 6.4 CONCLUSION This chapter presents some preliminary results obtained from 2D fast dynamics problems. The analysis results need to be compared with corresponding finite element solution and other numerical schemes. The solution obtained from the vibration of the short column seems to be agreeable. However, from the solution of punching test, it seems that the curlerror has not been diminished completely. Therefore, additional measures need to be taken 2D FAST STRUCTURAL DYNAMICS PROBLEM 128 Figure 6.3: Evolution of pressure distribution at different time instant for a tensile test to ensure compatibility condition. Chapter 7 CONCLUSIONS 129 CONCLUSIONS 7.1 130 SUMMARY AND CONCLUSION This work begins with a discussion on various terminologies and concepts in Kinematics. The terms associated to deformation gradient tensor is emphasised particularly as they play crucial role in the later part of the thesis. The first order conservation laws formulated from conservation of linear momentum, deformation gradient and total energy is derived. Moreover, different constitutive models, namely, nearly incompressible Neo-Hookean material model and linear elastic materials are described. The eigen structure of the governing equations is presented, as well. The eigenvalues represent the shear and longitudinal wave velocities and are important for obtaining time step size for the evolution of the solution in time. Contact and shock conditions are an integral part of the fast dynamics problems. The boundary of the problem domain is treated as contact. Hence, formulations related to free surface, sliding surface and sticking surface are derived. In this work, two particular stabilised schemes, e.g. one step Taylor-Galerkin scheme and SUPG scheme, are employed. Their similarity in structure makes them a favourable choice. In order to determine the suitability of these schemes, first they are employed to solve a pure advection equation. Extensive numerical experiments are performed on the schemes to evaluate their stability and consistency. For the one step Taylor-Galerkin scheme the stability limit is found to be at CFL number of 0.5 and it is second order convergent both in time and space. Two time integration schemes, namely, Forward Euler and TVD Runge-Kutta, are chosen for updating the SUPG solution in time. The Forward Euler time integration is found to be highly unstable and hence, discarded. The resulting space-time discretised SUPG formulation using TVD Runge-Kutta scheme displays second order convergence and has a stability limit of C = 0.6. A spectral analysis is also carried out in order to determine the dispersion and diffusion error of the schemes. It is found that for lower values values of CFL the numerical solution lags behind the analytical solution. However, for higher magnitude of CFL, the numerical solution leads the analytical solution at lower frequency followed by a sudden phase shift at higher frequency. The stabilised schemes are imple- CONCLUSIONS 131 mented to propagate a cosine wave and a step wave. Wiggles appear in the solution while propagating the step wave because of non-monotonic nature of both schemes. In later part of this work, the stabilised finite element formulations are extended to solve 1D hyperbolic equation with two degrees of freedom. Two different computational framework, namely, non-conservative and conservative framework is presented. For a linear problem, both conservative and non-conservative frameworks yield same solution. However, for a non-linear problem conservative framework is preferred. In each framework, two different formulations are presented which are indicated as “form A” and “form B”. The difference in those formulations arise from application of the Divergence theorem. The proposed formulations are employed to solve propagation of acoustic wave in a 1D gas tube. The numerical solution is compared against analytical solution. As expected, spurious oscillations arise in the solution due to presence of shock waves. One of the flexibilities of SUPG formulation is the choice of stabilisation parameter. A detailed study is performed on effect of stabilisation parameter in the solution. It is observed that for a given CFL there exists an optimum magnitude of stabilisation parameter. The analyses of one dimensional solid mechanics problems show that the SUPG-RK scheme performs better than the one step Taylor-Galerkin scheme. The improved accuracy in SUPG-TK scheme arises from the implementation of TVD-RK time integration. An indepth analysis reveals that “form A” yields slightly better results than “form B”. However, in 1D case, it is easier to work with “form B”, since the Neumann boundary conditions appear naturally. The formulation performs adequately when nearly incompressible NeoHookean material model is employed. The implementation of lumped mass matrix increases the stability of the scheme up to two times, but introduces dispersion error in the solution. Dispersion error can be minimised by refining the mesh. However, usage of consistent mass matrix is preferred since one of the aims of this study is to use relatively coarser mesh. The performance of the stabilised finite element schemes is compared to that of finite volume scheme. In the presence of shock discontinuity, spurious oscillations arise only in second order finite volume scheme as the first order finite volume scheme is monotonic. A shock CONCLUSIONS 132 capturing scheme is also employed for diminishing the wiggles appearing in the finite element solution. The performance of the shock capturing scheme is found to be satisfactory when compared with finite volume solution with slope limiters. The current work also includes preliminary analysis results of two dimensional fast dynamics problems. A novel idea is presented to ensure the satisfaction of compatibility condition. This strategy seems to work well for most of the cases, with few exceptions. Therefore, further improvements need to be made for ensuring curl-free involution. The solutions also need to be compared with finite volume results. 7.2 FUTURE RESEARCH DIRECTIONS This work has a tremendous potential for research in future. Some of the ideas are as follows: • Establish and reliable and robust framework for two-dimensional and three-dimensional applications. • Include plasticity in the formulations. This can be easily achieved by means of a sub-routine. • Incorporate different type of elements. Especially, implementation of linear triangular elements will offer significant advantage over contemporary finite element solvers • Introduce formulation for unstructured mesh. • Improve shock capturing scheme. This work can also be extended to other stabilised finite element schemes such as discontinuousGalerkin method [32]. Bibliography [1] http://www.wikipedia.org. [2] MSC/NASTRAN. User’s Manual - Version 66. The MacNeal-Schwendler Corporation, 1990. [3] ANSYS user’s guide. ANSYS Incorporation, 1997. [4] PAM-CRASH/PAM-SAFE Solver reference with manual. PAM SYSTEM International, 2002. [5] ABAQUS/Explicit. User’s Manual, Version 6.8. Dassault Systémes Simulia Corporation, 2008. [6] CSI Analysis Reference Manual For SAP2000. Computers and Structures Incorporation, 2009. [7] A DAMS , D. D., AND W OOD , W. L. Comparison of hilber-hughes-taylor and bossak α methods for the numerical integration of vibration equations. International Journal for Numerical Methods in Engineering 19, 5 (1983), 765–771. [8] A NGUELOV, R., L UBUMA , J. M.-S., AND M INANI , F. Total variation diminish- ing nonstandard finite difference schemes for conservation laws. Mathematical and Computer Modelling 51 (2010), 160 – 166. [9] A RGYRIS , J. Energy theorems and structural analysis: A generalized discourse with applications on energy principles of structural analysis including the effects of tem133 134 Bibliography perature and non-linear stress-strain relations. Aircraft Engineering and Aerospace Technology 26 (1954), 347 – 356. [10] BALSARA , D. S. von neumann stability analysis of smoothed particle hydrodynamics suggestions for optimal algorithms. Journal of Computational Physics 121, 2 (1995), 357 – 372. [11] BAZILEVS , Y., C ALO , V. M., T EZDUYAR , T. E., AND H UGHES , T. J. R. Y Zβ discontinuity capturing for advection-dominated processes with application to arterial drug delivery. International Journal for Numerical Methods in Fluids 54, 6-8 (2007), 593–608. [12] B ELYTSCHKO , T., L IU , W. K., AND M ORAN , B. Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons, Ltd, 2000. [13] B EN -A RTZI , M., AND FALCOVITZ , J. Generalized Riemann Problems in Compu- tational Fluid Dynamics. Cambridge University Press, 2003. [14] B ENSON , D. J. Computational methods in lagrangian and eulerian hydrocodes. Comput. Methods Appl. Mech. Eng. 99, 2-3 (Sept. 1992), 235–394. [15] B OCHEV, P. B., G UNZBURGER , M. D., AND S HADID , J. N. Stability of the SUPG finite element method for transient advection-diffusion problems. Computer Methods in Applied Mechanics and Engineering 193 (2004), 2301 – 2323. [16] B ONET, J., AND B URTON , A. J. A simple average nodal pressure tetrahedral el- ement for incompressible and nearly incompressible dynamic explicit applications. Communications in Numerical Methods in Engineering 14, 5 (1998), 437–449. [17] B ONET, J., AND G IL , A. J. Two step taylor-galerkin solution of lagrangian explicit dynamic solid mechanics. In Proceedings of the 8th World Congress on Computational Mechanics (2008). 135 Bibliography [18] B ONET, J., M ARRIOTT, H., AND H ASSAN , O. An averaged nodal deformation gradient linear tetrahedral element for large strain explicit dynamic applications. Communications in Numerical Methods in Engineering 17, 8 (2001), 551–561. [19] B ONET, J., AND W OOD , R. D. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press, 1997. [20] B OWER , A. F. Applied mechanics of solids. CRC Press, 2010. [21] B RADFORD , S. F., AND K ATOPODES , N. D. The anti-dissipative, non-monotone behavior of petrov-galerkin upwinding. International Journal for Numerical Methods in Fluids 33, 4 (2000), 583–608. [22] B ROOKS , A. N., AND H UGHES , T. J. Streamline upwind/petrov-galerkin formula- tions for convection dominated flows with particular emphasis on the incompressible navier-stokes equations. Computer Methods in Applied Mechanics and Engineering 32 (1982), 199 – 259. [23] B URTON , A. Explicit, Large Strain, Dynamic Finite Element Analysis with Applications to Human Body Impact Problems. PhD thesis, Swansea University, 1996. [24] C ALLE , J. L., D EVLOO , P. R., AND G OMES , S. M. Stabilized discontinuous galerkin method for hyperbolic equations. Computer Methods in Applied Mechanics and Engineering 194, 17 (2005), 1861 – 1874. [25] C ASTELLÓ , W. B., AND F LORES , F. G. A triangular finite element with local remeshing for the large strain analysis of axisymmetric solids. Computer Methods in Applied Mechanics and Engineering 198, 2 (2008), 332 – 343. [26] C ATABRIGA , L., C OUTINHO , A., AND T EZDUYAR , T. Compressible flow SUPG stabilization parameters computed from degree-of-freedom submatrices. Computational Mechanics 38 (2006), 334–343. 136 Bibliography [27] C ATABRIGA , L., DE S OUZA , D. A. F., C OUTINHO , A. L. G. A., AND T EZDUYAR , T. E. Three-dimensional edge-based SUPG computation of inviscid compressible flows with Y Zβ shock-capturing. Journal of Applied Mechanics 76 (2009), 021208. [28] C ESCOTTO , S., AND F ONDER , G. A finite element approach for large strains of nearly incompressible rubber-like materials. International Journal of Solids and Structures 15, 8 (1979), 589 – 605. [29] C ESCOTTO , S., AND Z HU , Y. Large strain dynamic analysis using solid and contact finite elements based on a mixed formulation; application to metal forming. Journal of Materials Processing Technology 45 (1994), 657 – 663. [30] C HUNG , J., AND H ULBERT, G. M. A time integration algorithm for structural dy- namics with improved numerical dissipation: The generalized-alpha method. Journal of Applied Mechanics 60, 2 (1993), 371–375. [31] C LOUGH , R. W. Original formulation of the finite element method. Finite Elem. Anal. Des. 7, 2 (Nov. 1990), 89–101. [32] C OCKBURN , B. Discontinuous galerkin methods for convection dominated problems. Tech. rep., University of Minnesota, 1999. [33] C OCKBURN , B., AND S HU , C.-W. The runge-kutta discontinuous galerkin method for conservation laws v: Multidimensional systems. Journal of Computational Physics 141, 2 (1998), 199 – 224. [34] C ODINA , R. Stability analysis of the forward euler scheme for the convectiondiffusion equation using the SUPG formulation in space. International Journal for Numerical Methods in Engineering 36, 9 (1993), 1445–1464. [35] C RISFIELD , M. A. Non-Linear Finite Element Analysis of Solids and Structures, Volume 1, Essentials. John Wiley & Sons, Ltd, 1996. 137 Bibliography [36] D IMAGGIO , S., AND B IENIEK , M. Nonlinear dynamics of flexible structures: A finite element approach. International Journal of Solids and Structures 32 (1995), 1179 – 1193. [37] D OHRMANN , C., H EINSTEIN , M. W., J UNG , J., K EY, S. W., AND W ITKOWSKI , W. R. Node-based uniform strain elements for three-node triangular and four-node tetrahedral meshes. International Journal for Numerical Methods in Engineering 47, 9 (2000), 1549–1568. [38] D ONEA , J., AND H UERTA ., A. Finite element methods for flow problems. John Wiley & Sons, Ltd, 2003. [39] F RANCA , L. P., F REY, S. L., AND H UGHES , T. J. Stabilized finite element meth- ods: I. application to the advective-diffusive model. Computer Methods in Applied Mechanics and Engineering 95, 2 (1992), 253 – 276. [40] G OTTLIEB , S., AND S HU , C.-W. Total variation diminishing runge-kutta schemes. Mathematics of Computation 67 (1998), 73–85. [41] G URTIN , M. E. An introduction to continuum mechanics. Academic Press, 1981. [42] H ALLQUIST, J. O., S TILLMAN , D. W., AND L IN ., T. L. LS-DYNA3D theoritical manual. Livermore Software Technology Corporation, 1991. [43] H ARTEN , A. High resolution schemes for hyperbolic conservation laws. Journal of Computational Physics 49, 3 (1983), 357 – 393. [44] H AUKE , G., AND D OWEIDAR , M. Fourier analysis of semi-discrete and space-time stabilized methods for the advective-diffusive-reactive equation: I. SUPG. Computer Methods in Applied Mechanics and Engineering 194, 1 (2005), 45 – 81. [45] H ENRIKOFF , A. Solution of problems in elasticity by the framework method. Journal of Applied Mechanics 8 (1941), 169–175. 138 Bibliography [46] H ILBER , H. M., H UGHES , T. J. R., AND TAYLOR , R. L. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engineering & Structural Dynamics 5, 3 (1977), 283–292. [47] H IRSCH , C. Numerical Computation of Internal and External Flows. Elsevier, 2007. [48] H OLZAPFEL , G. A. Nonlinear Solid Mechanics: A Continuum Approach for Engineering. John Wiley & Sons, Ltd, 2000. [49] H UGHES , T., S COVAZZI , G., AND T EZDUYAR , T. Stabilized methods for compressible flows. Journal of Scientific Computing 43 (2010), 343–368. [50] H UGHES , T., AND T EZDUYAR , T. Finite element methods for first-order hyperbolic systems with particular emphasis on the compressible euler equations. Computer Methods in Applied Mechanics and Engineering 45, 1-3 (1984), 217 – 284. [51] H UGHES , T. J. The finite element method: linear static and dynamic finite element analysis. Prentice-Hall, 1987. [52] H UGHES , T. J., F RANCA , L. P., AND M ALLET, M. A new finite element formula- tion for computational fluid dynamics: VI. convergence analysis of the generalized SUPG formulation for linear time-dependent multidimensional advective-diffusive systems. Computer Methods in Applied Mechanics and Engineering 63, 1 (1987), 97 – 112. [53] H UGHES , T. J., AND M ALLET, M. A new finite element formulation for compu- tational fluid dynamics: Iv. a discontinuity-capturing operator for multidimensional advective-diffusive systems. Computer Methods in Applied Mechanics and Engineering 58, 3 (1986), 329 – 336. [54] I DELSOHN , S. R., AND O ÑATE , E. Finite volumes and finite elements: Two ‘good friends’. International Journal for Numerical Methods in Engineering 37, 19 (1994), 3323–3341. 139 Bibliography [55] I URA , M., AND ATLURI , S. Dynamic analysis of finitely stretched and rotated threedimensional space-curved beams. Computers & Structures 29, 5 (1988), 875 – 889. [56] J OHNSON , C., NÄVERT, U., AND P ITKÄRANTA , J. Finite element methods for lin- ear hyperbolic problems. Computer Methods in Applied Mechanics and Engineering 45 (1984), 285 – 312. [57] J ONES , E. Stability analysis for maccormack’s time-splitting technique applied to a model conduction problem. Journal of Computational Physics 30, 3 (1979), 389 – 406. [58] K ARIM , I. A. A Two-step Taylor-Galerkin Formulation For Explicit Solid Dynamics Large Strain Problems. PhD thesis, Swansea University, 2011. [59] L E B EAU , G. J., R AY, S., A LIABADI , S., AND T EZDUYAR , T. SUPG finite element computation of compressible flows with the entropy and conservation variables formulations. Computer Methods in Applied Mechanics and Engineering 104, 3 (1993), 397 – 422. [60] L E B EAU , G. J., AND T EZDUYAR , T. E. Finite element computation of compress- ible flows with the SUPG formulation. Advances in finite element analysis in fluid dynamics FED 123 (1991), 21–27. [61] L EE , C. H. Development of a Cell Centred Upwind Finite Volume Algorithm for a New Conservation Law Formulation in Structural Dynamics. PhD thesis, Swansea University, 2012. [62] L EE , C. H., G IL , A. J., AND B ONET, J. Development of a finite volume algorithm for a new conservation law formulation in structural dynamics. In Proceedings of the 19th UK Conference of the Association for Computational Mechanics in Engineering (2011). 140 Bibliography [63] L EE , C. H., G IL , A. J., AND B ONET, J. Development of a cell centred upwind finite volume algorithm for a new conservation law formulation in structural dynamics. Computers & Structures (Under Review). [64] L EONARD , B. Note on the von neumann stability of explicit one-dimensional advection schemes. Computer Methods in Applied Mechanics and Engineering 118 (1994), 29 – 46. [65] L E V EQUE , R. J. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, 2002. [66] L UKKUNAPRASIT, P., AND K ELLY, J. M. Dynamic plastic analysis using stress resultant finite element formulation. International Journal of Solids and Structures 15, 3 (1979), 221 – 240. [67] M AIR , H. U. Hydrocode methodologies for underwater explosion structure medium interaction. In 66th Shock and Vibration Symposium SAVIAC (1995). [68] M ALVERN , L. E. Introduction to the mechanics of a continuous medium. PrenticeHall, 1969. [69] M ARRIOTT, H. Unstructured tetrahedral elements for explicit large strain, dynamic impact on human body Components. PhD thesis, Swansea University, 2000. [70] M ASE , G. T., AND M ASE , G. E. Continuum mechanics for engineers. CRC Press, 1999. [71] M ILLER , G., AND C OLELLA , P. A high-order eulerian godunov method for elasticplastic flow in solids. Journal of Computational Physics 167, 1 (2001), 131 – 176. [72] M ILLER , K., J OLDES , G., L ANCE , D., AND W ITTEK , A. Total lagrangian explicit dynamics finite element algorithm for computing soft tissue deformation. Communications in Numerical Methods in Engineering 23, 2 (2007), 121–134. 141 Bibliography [73] M UKHERJEE , D., G IL , A. J., AND B ONET, J. Streamline upwind petrov-galerkin solution of lagrangian explicit dynamic solid mechanics. In Proceedings of the 20th UK Conference of the Association for Computational Mechanics in Engineering (2012). [74] O ÑATE , E., C ERVERA , M., AND Z IENKIEWICZ , O. C. A finite volume format for structural mechanics. International Journal for Numerical Methods in Engineering 37, 2 (1994), 181–201. [75] P ERAIRE , J., P ERSSON , P., V IDAL , Y., B ONET, J., AND H UERTA ., A. A discontinuos galerkin formulation for lagrangian dynamic analysis of hyperelastic materials. In VII World Congress of computational Mechanics (2006). [76] R AYMOND , W. H., AND G ARDER , A. Selective damping in a galerkin method for solving wave problems with variable grids. Monthly Weather Review 104 (1976), 1583–1590. [77] R ISPOLI , F., S AAVEDRA , R., C ORSINI , A., AND T EZDUYAR , T. E. Computation of inviscid compressible flows with the V-SGS stabilization and Y Zβ shock-capturing. International Journal for Numerical Methods in Fluids 54, 6-8 (2007), 695–706. [78] R ISPOLI , F., S AAVEDRA , R., M ENICHINI , F., AND T EZDUYAR , T. E. Computation of inviscid supersonic flows around cylinders and spheres with the V-SGS stabilization and Y Zβ shock-capturing. Journal of Applied Mechanics 76, 2 (2009), 021209. [79] RODRÍGUEZ , M. C. A two-step taylor-galerkin algorithm applied to lagrangian dynamics. Master’s thesis, Swansea University, 2006. [80] ROE , P. Approximate riemann solvers, parameter vectors, and difference schemes. Journal of Computational Physics 135, 2 (1997), 250 – 258. [81] S CHODEK , D. L. Structures. Prentice-Hall, 2000. [82] S EGERLIND , L. J. Applied Finite Element Analysis. John Wiley & Sons, Ltd, 1884. 142 Bibliography [83] S ENGA , M. Improved SUPG formulations for compressible flows. PhD thesis, Rice University, 2005. [84] S HABANA , A. A. Computational Continuum Mechanics. Cambridge University Press, 2008. [85] S HAKIB , F., AND H UGHES , T. J. A new finite element formulation for computa- tional fluid dynamics: IX. fourier analysis of space-time galerkin/least-squares algorithms. Computer Methods in Applied Mechanics and Engineering 87, 1 (1991), 35 – 58. [86] S HAKIB , F., H UGHES , T. J., AND J OHAN , Z. A new finite element formulation for computational fluid dynamics: X. the compressible euler and navier-stokes equations. Computer Methods in Applied Mechanics and Engineering 89, 1-3 (1991), 141 – 219. [87] S HU , C.-W., AND O SHER , S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics 77, 2 (1988), 439 – 471. [88] S HU , C.-W., AND O SHER , S. Efficient implementation of essentially non-oscillatory shock-capturing schemes, ii. Journal of Computational Physics 83, 1 (1989), 32 – 78. [89] S IMO , J. A finite strain beam formulation. the three-dimensional dynamic problem. part i. Computer Methods in Applied Mechanics and Engineering 49, 1 (1985), 55 – 70. [90] S IMO , J., AND V U -Q UOC , L. On the dynamics in space of rods undergoing large motions - a geometrically exact approach. Computer Methods in Applied Mechanics and Engineering 66, 2 (1988), 125 – 161. [91] T EZDUYAR , T., AND PARK , Y. Discontinuity-capturing finite element formulations for nonlinear convection-diffusion-reaction equations. Computer Methods in Applied Mechanics and Engineering 59, 3 (1986), 307 – 325. 143 Bibliography [92] T EZDUYAR , T., S ENGA , M., AND V ICKER , D. Computation of inviscid supersonic flows around cylinders and spheres with the SUPG formulation and Y Zβ shockcapturing. Computational Mechanics 38 (2006), 469–481. [93] T EZDUYAR , T. E. Computation of moving boundaries and interfaces and stabilization parameters. International Journal for Numerical Methods in Fluids 43, 5 (2003), 555–575. [94] T EZDUYAR , T. E. Comments on adiabatic shock capturing in perfect gas hypersonic flows. International Journal for Numerical Methods in Fluids 66, 7 (2011), 935–938. [95] T EZDUYAR , T. E., AND O SAWA , Y. Finite element stabilization parameters com- puted from element matrices and vectors. Computer Methods in Applied Mechanics and Engineering 190, 3-4 (2000), 411 – 430. [96] T EZDUYAR , T. E., AND S ENGA , M. Stabilization and shock-capturing parameters in SUPG formulation of compressible flows. Computer Methods in Applied Mechanics and Engineering 195 (2006), 1621 – 1632. [97] T EZDUYAR , T. E., AND S ENGA , M. SUPG finite element computation of inviscid supersonic flows with Y Zβ shock-capturing. Computers & Fluids 36 (2007), 147 – 159. [98] T ORO , E. F. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer Berlin / Heidelberg, 2009. [99] T RANGENSTEIN , J. A., AND C OLELLA , P. A higher-order godunov method for modeling finite deformation in elastic-plastic solids. Communications on Pure and Applied Mathematics 44, 1 (1991), 41–100. [100] T URNER , M. J., C LOUGH , R. W., M ARTIN , H. C., AND T OPP, L. P. Stiffness and deflection analysis of complex structures. J. Aeronautical Society 23 (1956). 144 Bibliography [101] VON N EUMANN , J., AND R ICHTMYER , R. A method for the numerical calculation of hydrodynamic shocks. Journal of Applied Physics 21, 3 (1950), 232–237. [102] W EINBERG , E. J., AND K AAZEMPUR -M OFRAD , M. R. A large-strain finite ele- ment formulation for biological tissues with application to mitral valve leaflet tissue mechanics. Journal of Biomechanics 39, 8 (2006), 1557 – 1561. [103] W OOD , W. L., B OSSAK , M., AND Z IENKIEWICZ , O. C. An alpha modification of newmark’s method. International Journal for Numerical Methods in Engineering 15, 10 (1980), 1562–1566. [104] Z IENKIEWICZ , O. C., AND C HEUNG , Y. K. Finite element solution of field prob- lems. The Engineer (1965), 507–510. [105] Z IENKIEWICZ , O. C., AND M ORGAN , K. Finite elements and approximation. John Wiley & Sons, Ltd, 1983. [106] Z IENKIEWICZ , O. C., TAYLOR , R. L., AND N ITHIARASU ., P. The Finite Element Method for Fluid Dynamics. Elsevier, 2005. [107] Z IENKIEWICZ , O. C., TAYLOR , R. L., method. McGraw Hill, 2006. AND N ITHIARASU , P. The finite element