Dynamic Finite Element Investigation of Wave
Attack on Sea Dikes
Dynamic Finite Element Investigation of Wave
Attack on Sea Dikes
A Coupled Approach using Plate and Volume Elements
PROEFSCHRIFT
ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,
op gezag van de Rector Magnificus prof. ir. K.C.A.M. Luyben, voorzitter van het College voor Promoties,
in het openbaar te verdedigen op woensdag 30 march 2016 om 15:00 uur
door
Shuhong TAN
Master of Engineering
China University of Geosciences (Beijing), China geboren te Hebei, China
promotor: Prof. dr. M. A. Hicks Samenstelling promotiecommissie: Rector Magnificus voorzitter
Prof. dr. M. A. Hicks Technische Universiteit Delft Onafhankelijke leden:
Prof. dr. C. Jommi Technische Universiteit Delft Prof. dr. A. Scarpas Technische Universiteit Delft Prof. dr. H. Di Benedetto ENTPE, University of Lyon Dr. P. J. Vardon Technische Universiteit Delft Overige leden:
Prof. dr. ir. P. A. Vermeer University of Stuttgart and Deltares
Dr. L. Beuth Deltares
Prof. dr. ir. P. A. Vermeer heeft in belangrijke mate aan de totstandkoming van het proefschrift bijgedragen.
Keywords: dynamic finite element method, plate element, elastic viscoplas-tcity, sea dikes, excess pore pressure
Printed by: Ipskamp Drukker
Author: Shuhong TAN (谭淑红)
Copyright© 2016 by S. TAN
Author email: shuhong.tan@gmail.com
All rights reserved. No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording or by any information storage and retrieval system, without the prior permission of the author.
ISBN 978-94-6186-618-9
An electronic version of this dissertation is available at http://repository.tudelft.nl/.
Contents
Summary xi
Samenvatting xiii
1 Introduction 1
1.1 Aims and objectives . . . 4
1.2 Thesis outline . . . 4
2 Single phase formulation of dynamic finite element method 7 2.1 Literature review . . . 7
2.2 Virtual work formulation for dynamic problems . . . 9
2.2.1 Mathematical model . . . 9
2.2.2 Conservation laws . . . 11
2.2.3 Virtual work formulation . . . 11
2.3 Space discretisation for elaboration of the virtual work formulation 12 2.3.1 Shape functions . . . 13
2.3.2 Numerical Gaussian integration . . . 15
2.3.3 Lumped mass matrix . . . 17
2.4 Euler-Cromer time discretisation . . . 18
2.5 Strain smoothing for low order elements . . . 19
2.6 Local damping . . . 21
2.7 Absorbing boundary . . . 22
2.8 Study of the boundary effect . . . 24
2.8.1 Rigid bottom boundary . . . 25
2.8.2 Absorbing bottom boundary . . . 26
2.8.3 Conclusions . . . 30
3.1 Basic theory of two-phase soil behaviour . . . 32
3.1.1 Conservation of mass . . . 32
3.1.2 Conservation of momentum . . . 33
3.2 v − w formulation in finite element method . . . 34
3.3 Euler-Cromer time discretisation . . . 36
3.4 Local damping for two-phase problems . . . 37
3.5 Absorbing boundary for two-phase problems . . . 38
3.6 Numerical validations of one-dimensional wave propagation . . . . 39
3.7 Further validation of the one-dimensional consolidation problem with rigid bottom . . . 42
3.8 Numerical investigation of the absorbing boundary condition . . . 45
3.8.1 Absorbing boundary with varying δ/L (αw= 1) . . . 46
3.8.2 Absorbing boundary with varying αw (δ/L = 3) . . . 49
3.8.3 Conclusions . . . 52
4 Elastic beam and plate elements 53 4.1 Literature review on beams and plates . . . 54
4.2 Dynamic Euler-Bernoulli beam theory and FE modelling . . . 56
4.2.1 Basic theory . . . 56
4.2.1.1 Geometry and deformation behaviour . . . 57
4.2.1.2 Constitutive equations . . . 58
4.2.1.3 Equilibrium equations . . . 58
4.2.2 Finite element space discretisation . . . 60
4.2.2.1 Shape functions . . . 60
4.2.2.2 Discretisation of the virtual work formulation . . . 62
4.2.2.3 Numerical Gaussian integration . . . 63
4.2.2.4 Lumped mass matrix . . . 64
4.2.3 Euler-Cromer time discretisation . . . 65
CONTENTS
4.3 Extended dynamic Kirchhoff plate theory and FE modelling . . . . 68
4.3.1 Basic theory . . . 69
4.3.1.1 Geometry and deformation behaviour . . . 69
4.3.1.2 Kinematics . . . 70
4.3.1.3 Stretching of plates . . . 71
4.3.1.4 Bending of plates . . . 72
4.3.1.5 Boundary conditions . . . 74
4.3.2 Finite element space discretisation . . . 75
4.3.2.1 Discretisation of the virtual work formulations . . 75
4.3.2.2 Numerical space integration . . . 77
4.3.3 Enhanced lumped mass matrix of plate elements . . . 78
4.3.4 Local damping of plate elements . . . 80
4.3.5 Euler-Cromer time discretisation . . . 80
4.4 Validation of simple plate analysis . . . 82
4.4.1 Cantilever plate . . . 82
4.4.2 Circular plate with simply supported edge . . . 84
4.5 Coupling of plate and volume elements . . . 86
4.6 Validation of plates on soil . . . 89
4.6.1 Plate on elastic soil (plane strain) . . . 89
4.6.2 Validation with two-layer theory . . . 92
4.7 Concluding remarks . . . 95
5 Non-linear material models 97 5.1 Non-linear soil models . . . 97
5.1.1 Triaxial compression test . . . 98
5.1.2 Elastic perfectly plasticity (Mohr-Coulomb model) . . . 99
5.1.2.1 Elastic perfectly plastic behaviour . . . 99
5.1.2.2 Formulation of the Mohr-Coulomb model . . . 100
5.2 Elastic viscoplastic model for plates . . . 107
5.2.1 Literature review . . . 107
5.2.1.1 Overview of elastic viscoplastic model . . . 107
5.2.1.2 Modelling of bituminous concrete . . . 108
5.2.2 Typical creep behaviour of bituminous concrete . . . 109
5.2.3 Theory of elastic viscoplasticity . . . 112
5.2.4 Viscoplastic model for plates . . . 113
5.2.5 Creep threshold parameters . . . 115
5.2.6 Viscosity coefficient η . . . 118
5.2.6.1 Viscosity observed for road bituminous concrete . 118 5.2.6.2 Viscosity observed for hydraulic bituminous con-crete . . . 119
5.2.7 Work by Bonnier on elasticity of bituminous concrete . . . 120
5.2.8 Influence of temperature on parameters of the elastic vi-coplastic model . . . 124
5.3 Concluding remarks . . . 125
6 Numerical validation and verification of elastic viscoplastic plate model 127 6.1 Uniform stretching . . . 127 6.1.1 Non-frictional material . . . 128 6.1.2 Frictional material . . . 130 6.2 Plate bending . . . 133 6.3 Model test . . . 136 6.3.1 Experimental description . . . 136 6.3.2 Experimental results . . . 138
6.3.3 Description of the numerical analysis . . . 139
6.3.4 Numerical analysis results . . . 141
CONTENTS
7 Dynamic analysis of wave attack on sea dikes 147
7.1 Background . . . 147
7.2 Description of numerical analysis . . . 149
7.3 Numerical results of a preliminary calculation . . . 152
7.3.1 Quasi-static analysis of the gravity loading phase . . . 152
7.3.2 Dynamic analysis of the primary wave loading phase . . . . 154
7.3.3 Discussion . . . 160
7.4 Numerical results of successive wave loading . . . 162
7.5 Concluding remarks . . . 169
8 Concluding remarks and further research 171 8.1 Main conclusions . . . 171
8.2 Recommendations for further research . . . 173
References 175
Appendices
A Space discretisation of beam 187
B Plate space discretisation 189
C Numerical integration of elastic viscoplastic model for plates 195
Notation 199
Acknowledgements 209
Summary
Approximately 400 kilometres of Dutch sea dikes are protected by bituminous concrete revetments to prevent damage from erosion and repeated wave attacks during storms. The numerical analysis of sea dikes subjected to cyclic wave load-ing needs to consider the behaviour of the bituminous concrete revetment, and the behaviour of the subsoil including the generation and dissipation of the excess pore pressures, as well as the interaction between the revetment and subsoil. This thesis develops and evaluates numerical methods for investigating this problem, using the dynamic Finite Element Method with the coupling of plate and volume elements.
In dynamic finite element analysis, the arbitrary selected boundary generates wave reflections which normally causes oscillation problems. Absorbing boundary conditions are therefore adopted to minimise the wave reflection at the artificial boundaries. The effects of the absorbing boundary conditions are investigated in detail for both solid and water phases, and appropriate sets of parameters are recommended.
In order to model the bituminous concrete revetment, which is a thin layer of material with a high stiffness, a very fine mesh is often needed in order to ob-tain accurate results. However, an explicit time integration algorithm is normally needed for dynamic analysis, which reduces the critical time step size required for numerical stability if the same element type is used for both the revetment and the soil, severely affecting simulation performance. In order to avoid this problem, in this thesis a structural plate element is developed based on the classical Kirch-hoff thin plate bending theory for modelling thin layer materials instead of the volume element. For consideration of both stretching and bending problems, the Kirchhoff thin plate bending theory has been extended. Specifically, a more con-ventional 3-noded triangular plate element, with 1 translational and 2 rotational degrees of freedom per node, has been extended to include 2 extra translational degrees of freedom per node, i.e. giving a total of 3 translational and 2 rotational degrees of freedom per node. Moreover, an enhanced lumped mass matrix for the plate element has been constructed by considering both translational and rota-tional degrees of freedom. Numerical validations indicate that the use of plate elements for simulating the thin layer material provides a great improvement on the accuracy of the results, as well as a significant decrease in the computational cost.
For the limit state design of some engineering problems, non-linear material mod-els provide more accurate results than elastic material modmod-els. Here an elastic viscoplastic model is investigated and implemented for the new plate element, using the Drucker-Prager yield condition and a non-associated flow rule. Typical creep behaviour of the bituminous concrete is presented and a parametric study of the viscoplastic model has been carried out for modelling the behaviour of
Delft bituminous concrete. It is concluded that the viscosity value η for a partic-ular material changes from different conditions, that it depends not only on the material model and stress level, but also on the temperature.
The developed numerical methods have been used to investigate the behaviour of sea dikes subjected to wave attacks. The results of the bending moments and stresses in the bituminous concrete, as well as the excess pore pressures and dis-placements in the soil, are presented. The tendency of a so-called “shake down” behaviour is observed after 10 cycles of wave attack during a severe storm, indi-cating that the displacement of the dike converges asymptotically to a particular final deformation.
Samenvatting
Ongeveer 400 kilometer van de Nederlandse zeedijken zijn beschermd met as-faltbetonbekledingen om schade door erosie en herhaaldelijke golfaanval tijdens stormen te voorkomen. Voor numerieke analyse van zeedijken, blootgesteld aan cyclische golfbelasting, dient rekening te worden gehouden met het gedrag van asfaltbetonbekleding en het gedrag van de ondergrond, inclusief opwekking en afstroming van wateroverspanningen en de interactie tussen bekleding en on-dergrond. In dit proefschrift worden numerieke methodes ontwikkeld en beo-ordeeld om dit probleem te onderzoeken door gebruik te maken van de dynamische Eindige Elementen Methode met de koppeling van plaat- en volume elementen. In een dynamische eindige elementen analyse worden door de willekeurig vast-gelegde randen golfreflecties veroorzaakt die normaliter tot oscillatie problemen leiden. Daarom zijn absorberende randvoorwaarden toegepast om golfreflectie aan de kunstmatige randen tot een minimum te beperken. De gevolgen van ab-sorberende randvoorwaarden zijn in detail onderzocht voor zowel de grond als de water fase, en geschikte parametersets worden aanbevolen.
Voor modellering van asfaltbetonbekleding, een dunne materiaallaag met een hoge stijfheid, is vaak een zeer fijn elementengrid nodig om nauwkeurige resultaten te verkrijgen. Als echter hetzelfde soort elementen gebruikt wordt voor zowel bek-leding als ondergrond, wordt de grootte van de kritische tijdstap gereduceerd die nodig is om numerieke stabiliteit te waarborgen in expliciete tijdintegratie schema’s, wat de rekenprestatie ernstig be¨ınvloedt. Om dit probleem te ver-mijden is in dit proefschrift voor de modellering van dunne materiaallagen een constructief plaatelement ontwikkeld, gebaseerd op de klassieke theorie van Kirch-hoff voor doorbuiging van dunnen platen, in plaats van volume elementen. Om met zowel extensie als doorbuiging rekening te kunnen houden is de plaattheorie van Kirchhoff uitgebreid. Specifiek is een conventioneler drieknoops driehoekig plaatelement, met 1 vrijheidsgraad per knoop voor translatie en 2 voor rotatie, uitgebreid door 2 extra vrijheidsgraden voor translatie toe te voegen. Hiermee komt het element in totaal uit op drie vrijheidsgraden per knoop voor translatie en twee voor rotatie. Bovendien is een verbeterde gediagonaliseerde massamatrix opgesteld door zowel de vrijheidsgraden van translatie als rotatie te beschouwen. Numerieke validaties laten zien dat het gebruik van plaatelementen voor simu-latie van de dunne materiaallaag zowel voor een aanzienlijke verbetering van de resultaten als ook een significante reductie van de rekentijd zorgen.
Bij het ontwerpen van de uiterste grenstoestand van sommige ingenieursproble-men zorgen niet-lineaire materiaalmodellen voor nauwkeurigere resultaten dan elastische materiaalmodellen. Voor het nieuwe plaatelement is hier een elastisch visco-plastisch model onderzocht en ge¨ımplementeerd, waarbij het Drucker-Prager bezwijkcriterium en een niet-associatieve vloeiregel gebruikt wordt. Typisch kruipgedrag van het asfaltbeton wordt gepresenteerd en een parameterstudie van
het visco-plastische model is uitgevoerd voor de modellering van het gedrag van Delfts asfaltbeton. Het wordt geconcludeerd dat de viscositeitswaarde η voor een bepaald materiaal niet alleen afhankelijk is van de temperatuur maar ook van het materiaalmodel en spanningsniveau.
De ontwikkelde numerieke methodes zijn toegepast om het gedrag van zeedijken, blootgesteld aan golfaanval, te onderzoeken. De resultaten van zowel de buigende momenten en spanningen in het asfaltbeton als ook de wateroverspanningen en verplaatsingen in de ondergrond worden gepresenteerd. De neiging tot het zo-genaamde ‘shake down’ gedrag wordt waargenomen na 10 cycli van golfaanval. Dit geeft aan dat de verplaatsing van de dijk asymptotisch convergeert naar een bepaalde eindvervorming.
1
Introduction
A large part of the Netherlands lies below sea level, and is protected from flooding by a system of dikes, levees and dams. Currently, more than 400 kilometres of Dutch sea dikes have a bituminous concrete revetment, to prevent damage from both repeated wave attack during storms and erosion of the soil. Figure 1.1 shows the failure of a bituminous concrete revetment of a sea dike, as the result of repeated wave loading under storm conditions. The behaviour of the revetment and the subsoil under wave attack is not well understood, nor the processes of generation and dissipation of water pressures in the subsoil. Therefore, it is very important to study the dynamic behaviour of sea dikes during a super storm. However, the analysis of such a problem, which includes the coupling of the
solid and water phases, is challenging, especially when considering its dynamic nature. For existing geotechnical software, it is not possible to provide reliable solutions. Either it is possible to solve dynamic problems with larger deformations without being able to take account of the pore fluid (e.g. ABAQUS, LS-Dyna), or the pore fluid is taken into account but solvers are limited to the consolidation case excluding wave propagation and inertial effects (e.g. PLAXIS, ABAQUS). In some cases, dynamics and fluid coupling are considered, but limited to the undrained situation where the pore fluid does not escape the pores during the simulation (e.g. PLAXIS dynamics).
The dynamic material point method (MPM) is an extension of the classical finite element method, and has been proven to be able to model large deformation problems. Furthermore, it is possible to solve dynamic problems involving the coupling of the solid and fluid phases. Therefore, MPM is a promising tool for the dynamic analysis of wave attacks on sea dikes, not only for cases involving larger deformations, but also for considering the dynamic effect of the pore water pressure, or a more severe situation, e.g. liquefaction in the sand.
In the current MPM code developed by Deltares, volume elements, i.e. 4-noded tetrahedral elements, are used for the numerical modelling. However, for typical Dutch sea dikes, most of the bituminous concrete revetments have a thickness of only 20 to 30 cm in the wave attack zone. Compared to the subsoil underneath, the bituminous concrete revetment is therefore a very thin layer of material with a much higher stiffness. Unfortunately, the numerical modelling of such a thin layer material using volume elements is computationally expensive and sometimes leads to inaccurate results, especially when using 4-noded tetrahedral elements which use linear shape functions.
As an example, a cantilever plate is considered as shown in Figure 1.2, with a unit width, a length of 5 m, and a thickness of 0.5 m. The plate has a Young’s modulus of 100 MPa and a Poisson’s ratio of 0. A uniform load of 1 kPa is applied on the top surface of the plate. It is modelled by multiple layers of volume elements across the thickness of the plate. The deflections of the free end are listed in Table 1.1, in which the errors indicate the difference between the numerical results and the analytical solution of 0.075 m. It is clear that, when modelling the cantilever
Figure 1.2: Elevation view of a cantilever plate modelled with 2 lays of volume elements
Table 1.1: Deflection results for cantilever using different numbers of layers of volume elements
number of layers deflection w error time step size CPU time
(-) (m) (%) (s) (s)
1 0.0327 56.5 3.56×10−2 22.6
2 0.0515 31.4 1.52×10−2 96.0
4 0.0683 8.9 8.6×10−3 1599.2
plate with relatively few layers of volume elements, the numerical results are far too inaccurate. Even though the result can be significantly improved by using more layers of elements across the thickness, it reduces the time step size (required for numerical stability) which significantly increases the computational (CPU) time as shown in the table. In order to prevent inaccurate results, as well as reducing significantly the required computational effort, an an alternative option, the structural plate element, will here be considered and implemented to model the thin layer bituminous concrete revetment for the sea dike.
For sea dike revetments, hardly any large deformations occur to the bituminous concrete before cracking. Therefore, as a start to this research topic, the sea dike with a bituminous concrete revetment can be considered as a small deformation problem before cracking occurs, as shown in Figure 1.1, whereupon the soil un-derneath may flow out. In order to retain the capability of the dynamic analysis with coupled water and solid phases, the MPM framework is still adopted, but with the focus only on the small deformation aspect. Note that, for the material point method, when a single material point for each element is used and is lo-cated at the Gaussian integration point, it reduces to a traditional dynamic finite element method (DFEM). Therefore, the dynamic finite element method in the MPM framework developed by Deltares will be used to model the behaviour of the sea dike under dynamic wave loading during severe storms.
The aim is to incorporate plate elements coupled with tetrahedral elements in a dynamic FEM computation code. This initial implementation within a small deformation framework has four advantages:
• the small deformation of structural plate elements in geotechnical applica-tions covers the majority of situaapplica-tions in which these elements are applied in engineering practice;
• it allows comparison to available theoretical results for structural elements; • adding structural elements to a coupled solid-fluid solver poses a variety of challenges and few validations of the code are possible for large deforma-tion problems; whereas, for small deformadeforma-tions, the coupling can be tested against known limiting cases that allow a better validation of the code;
• the plate elements, at a later stage, can be taken into account in the material point method, for problems that allow the plate element to move through the soil body with large deformations, e.g. anchors, sheet piles.
1.1 Aims and objectives
The main aim of the thesis is to study the behaviour of sea dikes with bitumi-nous concrete revetments under repeated wave loading during severe storms. It involves the viscoplastic behaviour of the bituminous concrete, the soil behaviour during the repeated wave loading, and the generation and dissipation of the excess pore pressures in the soil. Above all, it includes the efficient modelling of soil-structure interaction by developing and implementing enhanced plate elements to model the revetment.
The basic objectives of the thesis are to:
• study the basic mathematics of the dynamic finite element method and the two phase equations, and investigate the effects of the absorbing boundary conditions for dynamic problems;
• extend Kirchhoff thin plate theory so as to be suitable for both stretching and bending problems, and construct an enhanced lumped mass matrix with consideration of both translational and rotational degrees of freedom for implementation into the dynamic finite element method;
• develop an advanced constitutive model, i.e. elastic viscoplasticity, for plate elements and investigate the creep behaviour of hydraulic bituminous con-crete;
• apply plate elements to model the bituminous concrete revetment on top of a sea dike slope subjected to dynamic wave attack during a storm, and study the behaviour of the bituminous concrete and the soil underneath, i.e. the dynamic generation and dissipation of the excess pore pressures in the soil.
1.2 Thesis outline
This thesis is arranged in eight chapters.
The mathematics of the dynamic finite element method are presented in Chapter 2. Firstly, a literature review of the finite element method is presented. The basic
1.2. THESIS OUTLINE
equations of DFEM are then reviewed, including the virtual work formulation, space and time discretizations. The lumped mass matrix is presented, with the advantage of computational and storage cost over the consistent mass matrix. Strain smoothing is discussed to deal with locking in low order elements. Fur-thermore, local damping and absorbing boundary conditions are discussed. The effects of the absorbing boundary conditions are investigated in detail and an appropriate set of parameters is recommended.
Chapter 3 discusses coupled dynamic two-phase FE analyses using a v − w formu-lation. The governing equations are solved for both the water and solid phases. The dynamic finite element procedure for the v − w formulation is also presented in detail, as is the use of Euler-Cromer time discretisation. This v − w formu-lation can well simulate the generation and dissipation of pore pressures under dynamic loading conditions. Absorbing boundary conditions are applied to the water phase and the relevant parameters are investigated.
Chapter 4 discusses the structural elements. It starts with a short review relating to previous research about plates, followed by the Euler-Bernoulli beam theory and numerical validations. Special attention is paid to the elastic plate element. The classical Kirchhoff plate theory, for the case of a 3-noded triangular element with 3 degrees of freedom at each node, is extended to include five degrees of freedom per node. This allows the consideration of both bending and stretching problems. An enhanced scheme is presented to construct the lumped mass matrix of the presented plate element, considering not only the translational degrees of freedom, but also the rotational degrees of freedom. Furthermore, the plate elements are coupled with 4-noded tetrahedral elements. It is possible to solve simple plate problems, and also problems with a plate resting on top of soil. Finally, numerical validations are presented to validate the implemented plate element against numerical and analytical solutions.
Chapter 5 presents the non-linear material models. It gives an overview of the Mohr-Coulomb and hypoplastic models for modelling the soil. The parameter de-termination for Baskarp sand is considered. Then, particular attention is focused on the elastic viscoplastic model. The basic theory of the elastic viscoplastic model is discussed, especially for plate elements. The creep behaviour of bitumi-nous concrete is studied in detail. A parametric study of the viscoplastic model is presented for a particular hydraulic bituminous concrete used for revetments of sea dikes.
In Chapter 6, numerical simulations of plates modelled as an elastic viscoplastic material are presented. Uniform stretching and plate bending are considered and compared with analytical solutions. A model test of an experiment is presented and compared with experimental results, in which elastic viscoplastic material parameters are determined for the bituminous concrete.
The application of wave attack on the bituminous concrete revetment of a sea dike is described in Chapter 7. It includes analyses of the two-phase behaviour of the
soil, and the elastic viscoplastic behaviour of the bituminous concrete revetment, as well as absorbing boundary conditions. Attention is paid to the dynamic effect on the displacements of the revetment, effective stresses and excess pore pressures in the soil, and the “shake sown” behaviour due to the cyclic wave loading. The concluding remarks and recommended future work are presented in Chapter 8.
2
Single phase formulation of dynamic finite
element method
In this chapter, the basic theory of the dynamic finite element method (DFEM) using an explicit algorithm will be discussed. Section 2.1 presents a short liter-ature review about the finite element method. In Section 2.2, the weak form of the momentum balance is considered, and it leads to the well known virtual work equation which is used to develop the DFEM equations. Section 2.3 presents a detailed discretisation scheme for solving the virtual work formulation; special at-tention is paid to the lumping procedure of the mass matrix. Section 2.4 provides an overview of the Euler-Cromer time discretisation of a single time step for the finite element system. Section 2.5 discusses the approach of strain smoothing to deal with locking in low order elements. Artificial local damping is discussed in Section 2.6. A lot of attention is paid to the application of the absorbing bound-ary conditions in Section 2.7, where the well-known Kelvin-Voigt element is used to limit the continuous creep of the boundary. A numerical example for studying the boundary effect is given in Section 2.8, which contains a rigid boundary study and a parameter sensitivity analysis of the absorbing boundary.
2.1 Literature review
The finite element method (FEM) is a numerical technique for obtaining approx-imate solutions to partial differential equations for boundary-value problems. It has been successfully used to analyse a wide range of problems in engineering and applied mathematics. The principal applications arise both in structural mechanics and in continuum mechanics, e.g. solid mechanics and fluid flow. In geomechanics the pores of the solid may contain a fluid and, as a consequence, both solid and fluid mechanics may be involved. In this section, a general review of FEM development will be discussed with a view to geomechanics.
The modern development of the finite element method began in the 1940s. Hren-nikof (1943) and McHenry (1943) studied stresses in one-dimensional bars and beams. Courant (1943) introduced shape functions over triangular subregions, to interpolate the whole region, in order to get numerical solutions. This is the first application where the elements were pieces of a continuum rather than structural members. A collection of papers by Argyris (1954a; 1954b; 1954c; 1955a; 1955b; 1955c) and Argyris and Kelsey (1960) illustrated the important role of energy principles by developing the method for structures.
Turner et al. (1956) used line, triangular and rectangular elements to decompose a continuum, and solve 1- and 2-dimensional problems. Clough (1960) introduced the name “finite element” when using triangular and rectangular elements in plane-stress analyses.
In 1961 a flat, rectangular plate bending element was developed by Melosh. It was then followed by the development of the curved shell bending element for axisymmetrical shells and pressure vessels by Grafton and Strome (1963). Later, FEM was extended to three-dimensional problems (Martin 1961; Gallagher et al. 1962; Melosh 1963).
Up to the early 1960s, FEM was mostly used for problems with small strains, elastic behaviour and static loadings, until Turner et al. (1960) considered large deformation and thermal-mechanical coupling problems. Gallagher et al. (1962) studied material nonlinearities. With the development of the consistent mass matrix, dynamic analysis was considered by Archer (1965), which was applied to the analysis of distributed mass systems for bars and beams in structural analysis. Costantino (1967) used the finite element method for the dynamic analysis of stress wave problems. FEM was then extended to be applied in more field problems, such as in the determination of the torsion of a shaft (Zienkiewicz and Cheung 1965), fluid flow (Martin 1969) and heat conduction (Wilson and Nickell 1966).
In 1967, Zienkiewicz and Cheung wrote the first finite element textbook. Be-lytschko and Hsieh (1973) presented the use of the method for the transient anal-ysis of large displacement problems with material non-linearities. They presented the formulation of the constant strain triangular element and the Euler-Bernoulli beam element. They compared the results with available experimental data for several examples. After that, Belytschko (1976) considered a problem with large deformation and nonlinear dynamic behaviour, where he improved the numerical techniques for solving the system equations.
Readers interested in a more detailed development of the finite element method are referred to textbooks, e.g. by Cook et al. (1989), Bathe (1996), Zienkiewicz et al. (2005) and Logan (2007).
2.2. VIRTUAL WORK FORMULATION FOR DYNAMIC PROBLEMS
2.2 Virtual work formulation for dynamic problems
In the finite element method, the weak form of the balance equation is required. The weak form is an integral expression that implicitly contains the differential equation, i.e. the strong form (Cook et al. 1989). The strong form states the condition at any point, whereas the weak form states the condition in an average sense.
2.2.1
Mathematical model
Considering a continuum domain V , the coordinates of an arbitrary point are denoted by a vector x = [x y z]T in the Cartesian coordinate system. The
displacement, velocity and acceleration vectors at any point are respectively, ˆ u = [ˆux uˆy uˆz]T ˆ v = [ˆvx ˆvy vˆz]T= ˙uˆ ˆ a = [ˆax aˆy ˆaz]T= ˙vˆ (2.1)
Stress and stain
The stress state at a point is mostly expressed in the matrix notation:
σ = σx σxy σxz σyx σy σyz σzx σzy σz (2.2)
As the shear stress components are such that σxy = σyx, σyz = σzy and σxz =
σzx, there are only six independent stress components. Therefore, it can be
represented in the so-called Voigt notation (Belytschko et al. 2000), which has storage and computational advantages in numerical calculations:
σ = [σx σy σz σxy σyz σzx]T (2.3)
As for the stress tensor, the strain tensor is introduced in both matrix and vector representations. The most general matrix notation is
ε = εx εxy εxz εyx εy εyz εzx εzy εz (2.4)
where εxy = εyx, εyz = εzy and εxz = εzx. In FEM, the strain is often cast in
the vector notation:
where γxy= 2εxy.
Linear elastic model
For a linear elastic material, the constitutive relation relates the stress σ to the strain ε by Hooke’s law as
σ = Dε (2.6)
where D is an elastic coefficient matrix, which is related to the Young’s modulus E and Poisson’s ratio ν:
D = E (1 − 2ν)(1 + ν) 1 − ν ν ν 0 0 0 ν 1 − ν ν 0 0 0 ν ν 1 − ν 0 0 0 0 0 0 1−2ν2 0 0 0 0 0 0 1−2ν 2 0 0 0 0 0 0 1−2ν2 (2.7)
Note that throughout this thesis, tensile stresses are considered to be positive while compressive stresses are negative.
For a purely linear elastic material, all the energy is stored when there is no damping effect. On introducing the shear strains γxy, etc., the rate of work ˙W
and the strain energy density W can be simply written as ˙ W = σT˙ε W = Z t 0 ˙ W dt = Z t 0 σT˙ε dt = Z t 0 εTD ˙ε dt =1 2ε TDε = 1 2σ Tε (2.8)
where ˙ε denotes the rate of strain. It can be derived from Equation (2.7) that
p = −Kεv (2.9)
where p is the mean stress with p = −13(σx+ σy+ σz), which is positive for
pressure, and εv is the volumetric strain with εv = εx+ εy+ εz. K is the bulk
modulus, defined as
K = E
3(1 − 2ν) (2.10)
The shear modulus G is defined as
G = E
2(1 + ν) (2.11)
In later chapters of this thesis, a water-saturated porous material will be consid-ered. A distinction will be made between pore water pressure pw and the mean
effective stress p0. Here it should be realised that the mean effective stress p0 is taken to be positive for compression.
2.2. VIRTUAL WORK FORMULATION FOR DYNAMIC PROBLEMS
2.2.2
Conservation laws
The law of conservation of mass states that the mass of the system must remain constant over time, as the system mass cannot change when there is no mass added or removed from the system. The mathematical form of the conservation of mass can be written as
∂ρ ∂t + ρ
∂ ˆv
∂x= 0 (2.12)
where ρ is the material density and ˆv is the velocity.
The conservation of momentum represents the equation of motion of a continuum, i.e. Newton’s second law of motion. It relates the motion or the kinetic of a continuum to the internal and external forces acting upon it. The mathematical form of the conservation of momentum is
LTσ + ρg = ρ¨uˆ (2.13)
where ρg is the force due to gravity, which is referred to as the unit body force in this thesis, and g is the gravity acceleration, where g = [gx gy gz]T. L is a
matrix operator given by
LT= ∂ ∂x 0 0 ∂ ∂y 0 ∂ ∂z 0 ∂ ∂y 0 ∂ ∂x ∂ ∂z 0 0 0 ∂z∂ 0 ∂y∂ ∂x∂ (2.14)
2.2.3
Virtual work formulation
To apply FEM, the conservation of momentum has to be transformed into the weak form, i.e. the well-known virtual work formulation, which is obtained by multiplying Equation (2.13) by a virtual displacement field δ ˆu, and integrating over the domain V currently occupied by the body (Pilkey and Wunderlich 1994; Cook et al. 1989): Z V δ ˆuT(LTσ + ρg) dV = Z V δ ˆuTρ¨u dVˆ (2.15) The virtual displacement δ ˆu must be kinematically admissible and continuous over the domain V .
The first part of Equation (2.15) can be expanded by using integration by parts and the divergence theorem (Zienkiewicz and Cheung 1967). This yields
Z V δ ˆuTLTσ dV = − Z V (Lδ ˆu)Tσ dV + Z S δ ˆuTτ dS (2.16)
where τ is the surface traction vector. It acts on the external surface, S, with unit normal vector n. Thus the traction is given by τ = σn.
By substituting Equation (2.16) into Equation (2.15), the virtual work equation can be written as:
Z S δ ˆuTτ dS + Z V δ ˆuTρg dV − Z V (Lδ ˆu)Tσ dV = Z V δ ˆuTρ¨u dVˆ (2.17)
which will be used in the formulation of the discrete finite element equations. The boundary conditions are ˆu = ¯u or τ = ¯τ , where ¯u and ¯τ are the prescribed displacements and tractions on complementary parts of the body surface S. The initial conditions at t = 0 are ˆu = ˆu0 and ˙u = ˙ˆ uˆ0.
2.3 Space discretisation for elaboration of the virtual
work formulation
It is often difficult or even impossible to solve boundary value problems analyt-ically. In order to derive approximate equilibrium equations, the finite element method is used by first decomposing the domain into small interconnected com-ponents called “finite elements”, as shown in Figure 2.1 with triangular elements for a 2 dimensional case and tetrahedral elements for a 3 dimensional case. Every finite element is linked to adjacent elements through shared nodes. The method has been described by many researchers as stated in Section 2.1. The basics will be briefly described in this chapter.
(a) 2D domain with triangular elements (b) 3D domain with tetrahedral elements
2.3. SPACE DISCRETISATION OF VIRTUAL WORK FORMULATION
2.3.1
Shape functions
In this study, special attention will be paid to 4-noded tetrahedral elements as shown in Figure 2.1(b). In dynamic analyses, such low order elements have the advantage that the resulting mass matrix can be lumped as described in Section 2.3.3. Considering 3D problems, large matrices will have to be introduced, which tend to be computationally time consuming. However, by using low-order ele-ments and lumped linear matrices the computational cost will remain restricted. For 4-noded tetrahedral elements, a continuous displacement field ˆu for element e is written as ˆ u ≈ N ue= [N1 N2 N3 N4] u1 u2 u3 u4 (2.18)
where the vector uirepresents the displacement of node i, i.e. ui= [uxi uyi uzi]T.
N is the shape functions matrix of an element, and Ni is the shape function of
node i: Ni= Ni 0 0 0 Ni 0 0 0 Ni (2.19)
where Niis the linear shape function of the 4-noded tetrahedral element, defined
as Ni= ai+ bix + ciy + diz 6V and 4 X i=1 Ni= 1 (2.20)
where V is the volume of the tetrahedron,
V = 1 6det 1 x1 y1 z1 1 x2 y2 z2 1 x3 y3 z3 1 x4 y4 z4 (2.21)
and ai, bi, ci, di are defined as
ai= det xj yj zj xk yk zk xl yl zl bi= − det 1 yj zj 1 yk zk 1 yl zl ci= − det xj 1 zj xk 1 zk xl 1 zl di= − det xj yj 1 xk yk 1 xl yl 1 (2.22)
Written in the local coordinate system, i.e. ξ (ξ, η, ζ), the shape functions are specified as N1= 1 − ξ − η − ζ N2= ζ N3= ξ N4= η (2.23)
With all nodal displacements known, the strain can be computed from
ε = Bue (2.24)
where B is the strain-displacement matrix defined as,
B = [B1 B2 B3 B4] (2.25)
and for node i
Bi= LNi= ∂Ni ∂x 0 0 0 ∂Ni ∂y 0 0 0 ∂Ni ∂z ∂Ni ∂y ∂Ni ∂x 0 0 ∂Ni ∂z ∂Ni ∂y ∂Ni ∂z 0 ∂Ni ∂x (2.26)
The integrals of the virtual work Equation (2.17) are performed on each individ-ual element, when sweeping over the elements in the mesh. Element matrices resulting from the integration are assembled together for all nodes in the mesh. Substituting the shape functions element by element, with ˆu ≈ N ueand δ ˆuT≈
δuTNT, the virtual work Equation (2.17) can be written as Nel X e=1 Z Ve ρNTN ¨uedVe= Nel X e=1 Z Ve NTρg dVe+ Nel X e=1 Z Se NTτedSe− Nel X e=1 Z Ve BTσedVe (2.27) where Nel is the total number of elements in the mesh; Ve represents the volume
of element e; τeis the traction applied on the surface of the triangular element
Se with τe= σen.
By defining the integrals in vector or matrix form, the above equation can be written as
Mcu = F¨ body+ Ftrac− Fint= Fext− Fint (2.28) Note that the vectors in Equation (2.28) are for the whole system. u and u¨ are the nodal accelerations and displacements of the whole system. Fext is the vector of external forces, which consists of both the surface traction force Ftrac
2.3. SPACE DISCRETISATION OF VIRTUAL WORK FORMULATION
and the body force Fbody; Mc denotes the consistent mass matrix; Fint denotes the internal forces due to stresses. They are defined as
Mc = Nel X e=1 Z Ve ρNTN dVe (2.29a) Fbody = Nel X e=1 Z Ve NTρg dVe (2.29b) Ftrac= Nel X e=1 Z Se NTτedSe (2.29c) Fint = Nel X e=1 Z Ve BTσedVe (2.29d)
2.3.2
Numerical Gaussian integration
The integrations in Equation (2.29) can be performed either analytically or by numerical integration. In finite element analyses, Gaussian integration is gen-erally used for the numerical integration. It is usually an approximation of the definite integral of a function by using the weights of specified points within the domain of integration.
On denoting the number of Gaussian integration points in one element as Neq,
the integration of a function f (x) over the current volume Ve can be computed
from Z Ve f (x) dVe= Z ξ Z η Z ζ f (ξ)|J | dξdηdζ ≈ Neq X q=1 wqf (ξq)Ve (2.30)
where |J | denotes the determinant of the so-called Jacobian matrix which equals the element volume Ve; wq is the weight of the integration point q at the local
position ξq.
The fewer the number of integration points, the less the computational cost. Considerable savings can thus be achieved by choosing an appropriate number of integration points. For the considered 4-noded tetrahedral element with linear shape functions, which has constant strain and constant stress in one element, a single Gauss integration point (Neq = 1) is enough to fully capture the elastic
energy, i.e. to obtain the exact solution. The single Gauss point is located in the centre of the tetrahedral element. The local coordinates of the central point are ξ = η = ζ = 14 and its weight is wq = 1 (Hughes 1987). Because of the linear
shape functions, the integral computed with Gaussian integration yields the same result as the analytical solution.
Using “element-by-element” techniques, the numerical integrations in Equation (2.29) are evaluated at the Gaussian point,
Mc = Z V ρNTN dV ≈ Nel X e=1 wqρNT(ξq)N (ξq)Ve (2.31a) Fbody = Z V ρNTg dV = Nel X e=1 wqρNT(ξq)gVe (2.31b) Fint = Z V BTσ dV = Nel X e=1 wqBT(ξq)σe(ξq)Ve (2.31c)
where wq is the integration weight associated with the Gaussian point q; σeis the
stress at the local position of Gaussian point q in element e; B(ξq) and N (ξq)
are the matrices at the local position of Gaussian point q.
Note that B is a constant matrix at any point of one element because of the linear shape functions, so the internal force Fint is also constant in one element. Because the traction force is applied on the surface but not throughout the vol-ume, the traction force Ftrac is only integrated over the surface area on which the traction is applied. Even though all the nodes of the volume elements are included, the traction force of the interior nodes will be zero.
For a 4-noded tetrahedral element, the surface traction is prescribed on a 3-noded triangular element (Ntri = 3), as shown in Figure 2.2. The so-called Gaussian
integration method is explained in detail for the triangular element by Hughes (1987). Considering a linear integration of the integrand, full numerical integra-tion is obtained by using three integraintegra-tion points (Nseq= 3) in one element. The
positions of the integration points are shown in Figure 2.2(c). On denoting the
(a) tetrahedral element (b) triangular element in global domain
(c) triangular element in lo-cal domain with three Gauss points
2.3. SPACE DISCRETISATION OF VIRTUAL WORK FORMULATION
number of triangular elements at the traction surface as NSel, the traction force
can be calculated as Ftrac= Z S NTτ dS ≈ NSel X e=1 Nseq X q=1 wqNT(ξq)τe(ξq)Se = NSel X e=1 Nseq X q=1 wqNT(ξq) 3 X i=1 Ntri(ξq)τe(x) ! Se (2.32)
τe(ξq) denotes the traction vector at the local position of the integration point q
in element e and it can be obtained by the integration of the nodal tractions τe(x)
at the surface element; Ntri(ξq) is the shape function of the 3-noded triangular
element at local position ξq; Se is the surface area of the triangular element e
where the traction forces are applied; wq is the local integration weight.
2.3.3
Lumped mass matrix
In traditional finite element schemes, the consistent mass matrix is usually used, although it adds considerable complexity to the dynamic algorithm and requires considerable computational and storage costs, especially for a dynamic approach. Alternatively, a lumped mass matrix is sometimes used in practice, which has advantages with regards to computational and storage cost over the consistent mass matrix.
Several methods have been developed to construct the lumped mass matrix, but some have only limited mathematical support (Belytschko et al. 2000). A more detailed explanation of different methods for constructing the lumped mass ma-trix can be found in the textbook by Hughes (1987). The most common method is to construct the lumped mass matrix from the consistent one: each (diagonal) entry of the lumped mass matrix then corresponds to the row sum of the con-sistent mass matrix. For this procedure it is necessary to obtain the concon-sistent mass matrix first before being able to build up the lumped mass matrix, which is also time demanding.
Another method to construct the lumped mass matrix is called “direct lumping”. For this approach, the total mass of the element is distributed to the nodes without explicitly knowing the consistent mass matrix.
Let MLdenote the global lumped-mass matrix and Methe lumped mass matrix of the 4-noded tetrahedral element e,
Me= m1 0 0 0 0 m2 0 0 0 0 m3 0 0 0 0 m4 (2.33)
with sub-matrix mi corresponding to the mass of node i and 0 being a null
matrix. For a three-dimensional finite element, the nodal mass matrix has the following form: mi= mi 0 0 0 mi 0 0 0 mi (2.34)
where the nodal mass mi takes one quarter of the element mass and is defined as
mi= ρVe/4.
Direct lumping gives an identical lumped mass matrix to the common method, but requires much less computational time and is easy to implement. The draw-back of using a lumped-mass matrix is a slight numerical dissipation of kinetic energy (Burgess and Sulsky 1992). A key motivation for direct lumping is that a diagonal mass matrix offers computational and storage advantages in certain simulations, notably in explicit time integration.
2.4 Euler-Cromer time discretisation
The system in Equation (2.28) is discretised in space, but is still continuous in time. To get a fully discrete system of algebraic equations, a time integration algorithm is applied, which is based on the finite difference method.
Al-Kafaji (2013) discussed in detail the integration procedures for time discreti-sation and their energy conservation properties. The Euler-Cromer scheme is known in literature as a modified Euler method (MEM) (Hahn 1991) and widely used in modelling dynamic problems (Erhart et al. 2006). The procedure is as follows.
On solving Equation (2.28) and applying forward-Euler time integration, the velocity yields
˙
ut+∆t= ˙ut+ ¨ut∆t (2.35)
where ˙ut and ˙ut+∆t are the nodal velocities at time t and t + ∆t, respectively. The incremental nodal displacement is obtained by integrating the nodal velocity using the backward-Euler rule
∆ut+∆t= ˙ut+∆t∆t (2.36)
The strain increment follows from
∆εt+∆t= B∆ut+∆t (2.37)
Now, the stress σt+∆t can be updated according to the constitutive model of the
material. Consequently the internal force Fint will be obtained using Equation
(2.31c). Together with the external force Fextusing Equations (2.31b) and (2.32),
2.5. STRAIN SMOOTHING FOR LOW ORDER ELEMENTS
The Euler-Cromer scheme is conditionally stable. This means that if the time step, ∆t, exceeds a critical time step, ∆tcr, the computed solution might diverge.
The critical time step ∆tcr corresponds to the time that the fastest stress wave
takes to cross the smallest element of the mesh. The critical time step is given by
∆tcr=
lmin
c (2.38)
lmin is the characteristic length of the element, which depends on the element’s
shape. For the tetrahedral element adopted in this thesis, it is the shortest distance from the side of the maximum area to the opposite node. c is the speed of wave propagation defined as
c = s
Ec
ρ (2.39)
where Ec is the constrained modulus of the material defined as
Ec= (1 − ν)E
(1 − 2ν)(1 + v) (2.40)
2.5 Strain smoothing for low order elements
The use of low order triangular elements to solve problems in elasticity dates back to the earliest years of the finite element method. However, there are sometimes difficulties when using low order elements, especially when determining the dis-placement field for a material that is nearly incompressible. For materials with a large bulk modulus, small errors in strain can result in large errors in stress. The mesh may lock when the incompressibility constraints from neighbouring elements are imposed. This locking can eventually propagate through the entire domain.
For illustrating the locking phenomenon, Figure 2.3 shows a discretised domain with 3-noded triangular elements. Elements A and B have the same free node. The displacements of the free node for the two elements are uAand uB if
consid-ering the incompressibility constraint. As shown in Figure 2.3, the two directions are different, so that the free node locks. The same locking happens to all the free nodes attached to the fixed boundaries. Because the material is incompressible, the entire mesh will be locked and yield an unrealistically stiff response.
For high order elements, it is common to prevent mesh locking by reduced in-tegration. For low order elements, a kind of strain smoothing can be applied, which was proposed by Detournay and Dzik (2006). This so-called Nodal Mixed
Figure 2.3: Discretised domain with three-noded triangular elements
Discretisation (NMD) technique overcomes the observed over-stiff behaviour due to the volumetric constraint.
This Nodal Mixed Discretisation technique involves determining the strain rates for each element in the usual manner and then decomposing them into a volu-metric strain rate ˙εv, and a deviatoric strain rate, ˙e:
˙ε = ˙e +1
3ε˙vI (2.41)
with ˙εv= ˙εxx+ ˙εyy+ ˙εzz and I = [1 1 1 0 0 0]T.
The nodal volumetric strain rate for node i is now determined by averaging the volumetric strain rate of the attached elements:
˙¯ εv,i= X e ˙ εv,eVe X e Ve (2.42)
where the sum is over all elements that are attached to node i. Ve is the volume
of element e.
Once the volumetric strain rate is evaluated for node i, the enhanced volumetric strain rate ˙¯εv for element e is determined by averaging the nodal values of the
considered element: ˙¯ εv= 1 nen nen X i=1 ˙¯ εv,i (2.43)
where nen is the number of nodes per element. The working assumption is that
deviatoric strain rates need not be enhanced; only volumetric components. As a result, the final strain rate within an element is redefined as
¯˙ε = ˙e +1
2.6. LOCAL DAMPING
where the deviatoric component is the same as for the uncorrected strain rate field. The new strains are then used to determine the stresses at integration points using the constitutive relation.
It should be clarified here that, for each material of the problem, the enhanced volumetric strain is applied separately.
Furthermore, although the strain smoothing procedure can well solve the locking problem for low order elements, it will also cause energy dissipation due to the numerical process.
2.6 Local damping
The dynamic finite element method can also be used for solving quasi-static problems, where the dynamic terms of the equilibrium equation must be damped out to provide static solutions. This leads to the application of damping in the convergence to a steady state. In geomechanics, a quasi-static analysis will often precede a dynamic analysis. This is also the case for the sea dike wave attack problem which will be discussed later in Chapter 7, as initial stresses caused by soil weight have to be calculated. Soil weight may be applied gradually, but even in this case somewhat oscillating stresses and strains may be obtained. In order to suppress such oscillations damping may be applied.
The use of viscous damping may introduce body forces that retard steady-state collapse and, in extreme cases, may influence the failure mode (Hart et al. 1988). Therefore, Cundall (1982, 1987) described the use of the local damping as an effective method to get a steady state solution.
A damping force term is added to the equation of motion given by Equation (2.28):
MLu = F¨ ext− Fint+ Fdamp (2.45)
where Fdamp is the damping force. For any single degree of freedom i of a considered problem, the damping force fidamp is defined by:
fidamp= −α|fiext− fint
i |sign(vi) (2.46)
where α is a dimensionless damping coefficient and sign(vi) is the sign of a
non-zero value of velocity vi,
sign(vi) =
vi
|vi|
(2.47)
For the considered degree-of-freedom, the damping force fidampis proportional to the out-of-balance force fext
2.7 Absorbing boundary
As mentioned in the previous section, local damping is very beneficial for a quasi-static problem. However, for a truly dynamic analysis, it cannot be applied as it is non-physical. In reality, material bodies under dynamic loading show both hysteretic material damping and radiation damping of the vibrations within the system. Otherwise, the system would oscillate indefinitely when subjected to these dynamic loads. Radiation damping is the process of losing energy through wave radiation in an infinite medium. However, in numerical simulations, when boundaries are introduced the radiation damping is erased. Therefore, artificial absorbing boundaries are prescribed at the boundaries.
Lysmer and Kuhlmeyer (1969) introduced a partial solution to this problem by using dashpots as the absorbing boundaries of the model. The drawback of the dashpots is that the displacements of the boundaries will be continually increasing as long as the dashpots receive stresses from the soil body. In order to limit the displacements of the boundary, a spring may be added parallel to the dashpot to obtain a Kelvin-Voigt type of boundary response (Al-Kafaji 2013) as shown in Figure 2.4, and this is the approach that will be adopted in this thesis.
Figure 2.4: Kelvin-Voigt element
When applying this kind of absorbing boundary condition, one additional term corresponding to the absorbing boundaries is added to the momentum Equation (2.28):
MLu = F¨ ext− Fint− Fvb (2.48)
where Fvbis the force due to the absorbing boundary of the domain.
Because the absorbing boundary force is applied on the surface but not through-out the volume, the force Fvb is only integrated over the surface area where the viscous traction is applied. Even though all nodes of the volume elements are included, the traction force of the interior nodes will be zero.
The Newton-Cotes integration scheme is applied instead of Gaussian integration, because Newton-Cotes integration is not only simple to implement and compu-tationally fast, but is also accurate for springs (Van Langen 1990). According to the Newton-Cotes integration scheme, the integration points of a triangular
2.7. ABSORBING BOUNDARY
element are placed at the same positions as the three nodes, and the weighting coefficient wq for each integration point is 1/3.
The absorbing boundary force Fvb is computed from Fvb= Z S NTτvbdS ≈ NSel X e=1 3 X q=1 wqNT(ξq)τ vb e (ξq)Se (2.49)
where NSel is the number of boundary triangular elements where viscous
bound-aries are applied; τvb is the viscous traction vector and τvb
e (ξq) is the viscous
traction at the integration point q of element e; S is the integrating area of the viscous boundary and Seis the surface area of element e.
Since integration points are located at the element nodes, shape functions N have a value 1 at the corresponding node and a value 0 at other nodes. Equation (2.49) can be simply written as
Fvb≈ NSel X e=1 3 X q=1 1 3τ vb e (ξq)Se (2.50)
Considering the two parts of the Kelvin-Voigt element, the viscous traction τvb can be analysed in two parts; through the spring tractions τspring and dashpot
tractions τdashpot, i.e.
τvb= τspring+ τdashpot (2.51) Formulation for spring
The spring tractions are calculated using spring coefficients and displacements as follows:
τspring= Ku (2.52) where K is the diagonal spring matrix for computational efficiency reasons and contains the spring coefficients. Note that K is the diagonal spring matrix as long as the x, y or z axis is normal to the boundary; otherwise a rotation matrix is applied.
For a node i at the absorbing boundary, the spring matrix yields
ki= ˆ kn 0 0 0 kˆs 0 0 0 ˆks (2.53)
where ˆkn and ˆks are coefficients in the normal and tangential directions,
respec-tively. They are determined by the coefficients of the spring, which is considered as a virtual layer with a thickness of δ that goes beyond the boundaries of the mesh, which has a range of 0 ≤ δ ≤ ∞. Hence,
ˆ kn= Ec δ , ˆ ks= G δ (2.54)
The equations indicate that, when δ = 0, the Kelvin-Voigt element reduces to a rigid boundary; when δ = ∞, the Kelvin-Voigt element reduces to a dashpot boundary.
The shear modulus G given in Equation (2.11) and constrained modulus Ecgiven
in Equation (2.40) are both adopted from the material that is adjacent to the boundary where the Kelvin-Voigt element is prescribed.
Formulation for dashpot
The dashpot tractions are calculated using dashpot coefficients and velocities as follows:
τdashpot= C ˙u (2.55) where C is the diagonal dashpot matrix for computational efficiency reasons and contains the dashpot coefficients. Note that C is the diagonal spring matrix as long as the x, y or z axis is normal to the boundary; otherwise a rotation matrix is applied.
For a node i at the absorbing boundary, the dashpot matrix yields
ci= ˆ cn 0 0 0 ˆcs 0 0 0 ˆcs (2.56)
where ˆcn is the normal coefficient and ˆcsis the tangential coefficient of the
dash-pot. They are defined as ˆ
cn= αnρcp, cˆs= αsρcs (2.57)
where αn and αs are dimensionless parameters in the two directions, and ρ is
the material density which defines the propagation speeds of the primary and secondary waves, cpand cs, respectively, i.e.
cp= s Ec ρ , cs= s G ρ (2.58)
Note that, when the dashpot coefficient αn = αs = 0, the Kelvin-Voigt element
reduces to a single spring; when αn = αs= ∞, the Kelvin-Voigt element reduces
to a rigid boundary.
2.8 Study of the boundary effect
A one dimensional column of unit height L = 1 m is considered here, being instantaneously loaded by a compressive surface traction of q = -1 kPa. The
2.8. STUDY OF THE BOUNDARY EFFECT
column is modelled as a linear elastic material with a Young’s modulus of E = 100 kPa, a Poisson’s ratio of ν = 0 and a density of ρ = 1800 kg/m3. A very fine mesh is required to properly capture the propagating front of the idealized loading, so the column is discretised by 1000 rows of elements in the vertical direction with 6 tetrahedral elements in each row and a single Gaussian point in each element. Roller boundaries are prescribed along the sides so that only longitudinal displacements along the column are possible. A rigid boundary con-dition, i.e. fully fixed, is prescribed at the bottom of the column. For comparison, an absorbing boundary condition with Kelvin-Voigt elements is prescribed at the bottom to reduce the reflection of stress waves, and a parametric study of the absorbing boundary is presented. In order to study the stresses generated by the compressive traction, the body force (i.e. self weight) of the material is not considered.
2.8.1
Rigid bottom boundary
The critical time step is ∆tcr= 7.82 × 10−5 s following Equation (2.38). In order
to obtain a stable solution, a time step size of ∆t = 7.66 × 10−5 s was chosen in this calculation, which is 98% of the critical time step.
Figure 2.5(a) shows the wave front with the stress of -1 kPa at time t = 0.05 s. It has travelled downwards to a depth of 0.37 m which is more or less exactly the analytical result of 0.373 m, as calculated from the wave speed of cp= 7.454 m/s
following from Equation (2.58). At the rigid boundary the wave is reflected back with double magnitude, so that the stress wave propagates upwards with a value of -2 kPa, as shown in Figure 2.5(b). At time t = 0.2 s, the reflected wave front has travelled upwards to a distance of 0.476 m from the bottom, which is close to
(a) wave travels downwards (t = 0.05 s) (b) wave travels upwards (t = 0.2 s)
Figure 2.5: Wave propagation along the column with rigid boundaries prescribed at the bottom
Figure 2.6: Displacements and stresses varying with time for the rigid boundary condition
the analytical value of 0.485 m. When the wave front reaches the top surface, it is reflected and then travels downwards. The reflected wave has the same amplitude as the applied one, but with an opposite sign. As shown in Figure 2.6, the stress in the column is oscillating between -2 kPa and 0 kPa, while reflecting up and down from the bottom and top surfaces.
The steep wave front shown in Figure 2.5 indicates that the mesh is sufficiently fine to properly capture the wave propagation. The observed high frequency oscillations are attributed to the numerical process. This can be suppressed by applying the load gradually over a couple of time steps.
2.8.2
Absorbing bottom boundary
When the rigid boundary is replaced by the Kelvin-Voigt element, no difference occurs before the wave front reaches the bottom of the column. In order to study the choice of parameters of the Kelvin-Voigt element, a parametric study of the absorbing boundary has been carried out with different values of the spring coef-ficient δ and the dashpot coefcoef-ficient α. (Since the absorbing boundary condition is only prescribed in the direction normal to the bottom boundary, the dashpot coefficient α is used in this section instead of αn.)
An investigation of the restrictions of the Kelvin-Voigt element on the critical time step size has been done by Al-Kafaji (2013), and he proposed an additional
2.8. STUDY OF THE BOUNDARY EFFECT
criterion for the critical time step size for a one dimensional problem: ∆tcr=
2αδ cp
(2.59) where cp is the primary propagation speed defined in Equation (2.58).
Absorbing bottom boundary with varying δ/L (α = 1)
Lysmer and Kuhlmeyer (1969) studied the reflection of body waves from a viscous boundary considering only dashpots. They concluded that a value of α = 1 for the dashpot leads to nearly perfect absorption, i.e. non-reflecting conditions. Therefore, this parameter study starts with α = 1 for the dashpot, and different values of δ/L = 0.1 and 1 for the spring.
The critical time step size given by Equation (2.38) is ∆tcr = 7.82 × 10−5 s
for tetrahedral elements. For the absorbing boundary with α = 1, the critical time step size is 0.0268 s for δ = 0.1 and 0.268 s for δ = 1 (following Equation (2.59)). These values are much bigger than the critical time step size given by the first criterion from Equation (2.38). Hence, the smaller value is considered as the critical time step for this problem. Based on a stability study (considering a range of values of ∆t), a time step size of ∆t = 4.69×10−5s provides a stable solution for this problem, which is 60% of the critical time step size. It is therefore concluded that the above criteria are insufficient and that, for Kelvin-Voigt boundaries, an additional critical time step criterion is needed.
Figure 2.7 shows the reflected wave front propagating upwards in the column at different times with δ/L = 0.1 and δ/L = 1 for the spring. The virtual thickness is 0.1 m and 1 m for the two cases and the thicker one indicates a softer spring. When the ratio δ/L is small, i.e. 0.1, as shown in Figure 2.7(a), the absorbing
(a) stiff spring with δ/L = 0.1 (b) soft spring with δ/L = 1
Figure 2.7: Reflected wave propagation along the column with absorbing bound-aries prescribed at the bottom (α = 1)
boundary behaves closer to a rigid boundary and the reflection is almost double the stress. Whereas, when the ratio δ/L is bigger, i.e. 1, as shown in Figure 2.7(b), the reflected stress is considerably reduced, to 38% of the applied pressure. When the wave travels downwards and reaches the bottom, part of the energy transmits into the absorbing boundary, and the spring is compressed, which means the elastic energy is stored in the spring. When the wave travels up-wards, the deformation of the spring returns back, so the spring releases the energy. This results in the increase of the stress at the bottom.
Figure 2.8(a) shows the displacements of the bottom with different values of δ/L for the spring. Analytically, the static displacement of the bottom surface can be calculated from
ubottom=
q
Ecδ (2.60)
So the static displacement of the bottom is found to be 0.01δ. With different δ values of 0.1 m, 1 m, and 3 m, the numerical solutions are 0.001 m, 0.01 m and 0.03 m, respectively, which shows exact agreement with the analytical solution. Figure 2.8(b) shows the relative displacement of the top surface, which is obtained by subtracting the displacement of the bottom surface from the top surface dis-placement. The relative displacement for δ/L = 0.1 keeps oscillating. The cases with δ/L = 1 and δ/L = 3 reach the static state very quickly, even though there is some overshoot at the beginning. The final static results for the 3 cases show the same settlement of 0.01 m, which is exactly the static analytical solution. Figure 2.8(b) shows the trend that the longer the spring, the less the reflection. It indicates the main conclusion that perfect damping is obtained in the limit of δ = ∞, i.e. in the case of a dashpot without a spring.
(a) displacement of the bottom (b) relative displacement
2.8. STUDY OF THE BOUNDARY EFFECT
Absorbing bottom boundary with varying α (δ/L = 1)
A further sensitivity analysis has been done for the dashpot parameter α using δ/L = 1 for the spring.
Following Equation (2.59), for the spring thickness of δ = 1 m, the critical time step is 0.1342 s for α = 0.5 and 0.5366 s for α = 2. They are much bigger than the critical time step of ∆tcr = 7.82 × 10−5 s given by Equation (2.38), so the
smaller value is considered as the critical time step. Based on the stability study, a larger value of α significantly reduces the time step size. For the analysis with α = 0.5, a time step size of ∆t = 6.26 × 10−5 s is necessary to capture a stable solution, which is 80% of the critical time step size; whereas, for the one with α = 2, a time step size of ∆t = 2.35 × 10−5 s is necessary to obtain a stable
solution, which is 30% of the critical time step size.
Figure 2.9 shows the reflected wave front propagating upwards in the column at different times, with α = 0.5 and α = 2 for the dashpot. As shown in Figure 2.9(a), for α = 0.5 the reflected wave front is travelling upwards with a stress of -0.73 kPa, and the stress at the bottom increases up to -1.30 kPa. In Figure 2.9(b) with α = 2, the reflected wave front is -1.32 kPa when travelling upwards. The reflected stress at the bottom increases by up to 53% of the applied pressure. Figure 2.10 shows the relative displacement of the top surface of the column for different values of α. The final static results for the 3 cases show the same settle-ment of 0.01 m, which is exactly the static analytical solution. However, the one with α = 0.5 presents some unrealistic phenomena when the wave front reaches the bottom surface, because of the negative reflective stress at the boundary; moreover, there are more ossifications with α = 0.5. With the α value of 2, similar phenomena occur when the wave front reaches the top and bottom
sur-(a) α = 0.5 (b) α = 2
Figure 2.9: Reflected wave propagation along the column with absorbing bound-aries prescribed at the bottom (δ/L = 1)
Figure 2.10: Relative displacements for different values of α with absorbing boundary condition (δ/L = 1)
faces. When the α value equals to 1, the system behaves more reasonably, and the relative displacement reaches the static solution faster.
Considering the present data, it would seem reasonable to use values of α close to 1.0 for the dashpot.
2.8.3
Conclusions
The introduction of the absorbing boundary with the Kelvin-Voigt element effi-ciently reduces the reflection of the stress at the bottom boundary. The parameter study shows that, the larger value of δ/L, i.e. the softer the spring, the less the reflection at the bottom; the larger the value of α for the dashpot, the more the reflection at the bottom. Furthermore, the stability study of the time step size indicates that using the absorbing boundary reduces the required time step size, and that the critical time step criterion for the absorbing boundary overes-timates the critical time step. Further study of the absorbing boundary effect on the critical time step is required.
