DISCRETE MODEL OF STRUCTURE BUILT OF SATURATED PORO-PLASTIC MEDIUM

In this Chapter, we describe the plane strain coupled discrete beam lattice model of structure built of saturated poro-plastic cohesive material. The coupled discrete beam lattice model is based on Voronoi cell representation of the domain with cohesive links as inelastic Timoshenko beam finite elements enhanced with additional kinematics in terms of embedded strong discontinuities. The coupling between the solid phase and the internal (pore) fluid is governed by Biot's porous media theory and Darcy's law for pore fluid flow. The internal fluid pore pressure field is approximated with CST (Constant Strain Triangle) finite elements, which coincide with Delaunay triangles. By exploiting Hammer quadrature rule for numerical integration on CST finite elements and duality property between Voronoi diagram and Delaunay triangle representations, the proposed formulation results with an additional pressure-type degree of freedom placed at each node of Timoshenko beam finite element. The application of the model is illustrated through several numerical simulations.

ACKNOWLEDGEMENTS

The work on this thesis has been supported by the University of Sarajevo, Faculty of Civil Engineering and the French Government scholarship for doctoral level. These supports are gratefully acknowledged.

I would like to express my warmest thanks

To my supervisors:

• professor Adnan Ibrahimbegovic, for the patient guidance, immense knowledge generously shared, continuous support, encouragement and motivation throughout this journey, which helped me to move forward and to grow, not just as a researcher but as a person. I would also like to thank his wife Nita, who is always very kind and welcoming, and with whom I shared many laughs.

• professor Samir Dolarević, for persistent help and support in every situation, whether it was about research or personal struggles.

To the reviewers of the thesis, professor Željana Nikolić and professor José Luis Pérez Aparicio, and the members of the jury, emeritus professor Muhamed Zlatar, professor Abdellatif Ouashine, professor Hermann G. Matthies To all my friends and colleagues at the Faculty of Civil Engineering, University of Sarajevo for their help, many constructive discussions, and for all funny talks we shared throughout this period.

To the wonderful new friends I met during my PhD studies at UTC. Thanks to Tea and Ivan for always being there for me, my Sarita and Milan for all their support and encouragement, Pablito -FEAP guru for his friendship, patience and selfless help with FEAP and also thanks to Abir, Adela, Xuan Nam, Andjelka and Mijo for all the fun we had on this journey.

To all my friends and especially my best friend Lejla, for all supporting words and encouragement.

And finally, to my family. Thank you for everything. Without you, this would not be possible.

xii

ABSTRACT Title: Analysis of pore pressure influence on failure mechanisms in structural systems

This thesis studies the issue of the overall safety of structures built of heterogeneous and pore-saturated materials under extreme loads in application to fluid-structure interaction problems, such as the dam-reservoir interaction. We propose a numerical model of interaction capable of predicting main tendencies and overall behavior of pore-saturated dam structure interacting with the reservoir in failure analyses of practical interest. The proposed numerical model is first presented in two-dimensional (2D) framework and later extended to three-dimensional (3D) framework.

We consider the structure built of porous cohesive material. We assume that the external fluid in interaction with the structure acts as a source of pore saturation.

We model the response of the pore-saturated structure with the coupled discrete beam lattice model based on Voronoi cell representation of domain with inelastic Timoshenko beam finite elements enhanced with additional kinematics in terms of embedded strong discontinuities acting as cohesive links. The coupling between the solid phase and the pore fluid is handled with Biot's porous media theory, and Darcy's law governing the pore fluid flow. The numerical consideration of internal coupling results with an additional pressure-type degree of freedom placed at each node of the Timoshenko beam finite element, which is later used at the fluidstructure interface.

The confined conditions met for external fluid placed in the reservoir enable the modeling of external fluid motion with the acoustic wave theory. For the numerical representation of the external fluid limited to small (irrotational) motion, we choose a Lagrangian formulation and the mixed displacement/pressure based finite element approximation. The end result are the displacement and pressure degrees of freedom per node of external fluid finite elements, which allows for the issue of the fluid-structure interface to be solved in an efficient and straightforward manner by directly connecting the structure and external fluid finite elements at common nodes. As a result, all computations can be performed in a fully monolithic manner.

All numerical implementations and computations are performed with the research version of the computer code FEAP (Finite Element Analysis Program).

The proposed numerical models of structure, external fluid and ultimately numerical model of interaction are validated in the linear elastic regime of structure response by comparing computed results against reference values obtained either with analytical solutions or continuum models. The numerical simulations in the nonlinear regime of structure response are performed with the aim to demonstrate the proposed coupled discrete beam lattice model capabilities to capture complete macro-scale response and failure mechanisms in pore-saturated structures. Finally, the proposed numerical model of interaction ability to deal with the progressive localized failure of a dam structure built of porous cohesive material under damreservoir interaction for a particular loading program was tested. To account for the temperature effects, the thermal coupling is introduced in the numerical model of the structure. La structure est considérée comme un milieu poreux saturé constitué d'un matériau cohésif. On suppose que le fluide externe en interaction avec la structure agit comme une source de saturation des pores. La réponse de la structure en matériau saturé est décrite avec un modèle lattice discrete couplé de type poutre, basé sur la discrétisation du domaine avec la tessellation de Voronoï, où les liens cohésifs sont représentés par des poutres de Timoshenko non linéaires avec un champ de déplacements enrichi en termes de discontinuités fortes. Le couplage entre la phase solide et le fluide dans les pores est traité avec la théorie de Biot et la loi de Darcy décrivant l'écoulement d'un fluide à travers d'un milieu poreux.

La prise en compte numérique du couplage interne ajoute un degré de liberté supplémentaire du type pression à chaque noeud de l'élément fini de Timoshenko, qui est ensuite utilisé pour résoudre les problèmes d'interface entre la structure et le fluide.

On considère que le fluide externe dans le réservoir est limité à des petits mouvements, ce qui nous permet de le modéliser avec la théorie des ondes acoustiques.

Pour cela, la formulation lagrangienne avec l'approximation mixte déplacementpression est choisie. Le traitement de l'interface fluide-structure dans le modèle numérique d'interaction est résolu d'une manière simple et efficace. Notamment, les éléments finis de la structure et du fluide externe partagent les mêmes degrés de liberté dans les noeuds communs, permettant ainsi la résolution du système d'équations avec une approche de calcul monolithique.

Toutes les implémentations et les simulations numériques sont effectués avec la version recherche du code informatique FEAP (Finite Element Analysis Program).

Les modèles numériques proposés pour la structure, le fluide externe et le modèle d'interaction sont validés dans le régime élastique linéaire en comparant les résultats calculés avec les valeurs de référence obtenues soit avec des solutions analytiques, soit avec des modèles continus. Les simulations numériques dans le régime non linéaire ont comme but de démontrer les capacités du modèle proposé de capturer la réponse complète à l'échelle macro et les mécanismes de rupture des structures en matériaux saturés. Enfin, la capacité du modèle d'interaction proposé de traiter la défaillance localisée progressive d'un barrage construit en matériau cohésif poreux sous l'interaction barrage-réservoir a été testé pour un programme de chargement spécifique. Pour prendre en compte les effets de la température, le couplage thermique est introduit dans le modèle numérique de la structure.

Mots clés: modèle lattice discrete cohésif de type poutre, poutre de Timoshenko, discontinuité intégrée, couplage, pression, milieu poreux, théorie de Biot, théorie des ondes acoustiques, approximation mixte déplacement-pression, interaction fluide-structure, analyse de défaillance, sécurité globale, température xvi SAŽETAK Naslov: Analiza uticaja pornog pritiska na mehanizme loma konstruktivnih sistema Ova doktorska teza analizira problem globalne sigurnosti konstrukcija od heterogenih i zasićenih materijala izloženih djelovanju ekstremnih opterećenja, s primjenom na probleme interakcije fluida i konstrukcije kao što je slučaj interakcije brana-rezervoar. U tezi je predložen numerički model s kojim je moguće predvidjeti glavne tendencije i globalno ponašanje konstrukcije brane izgra dene od zasićenog materijala u interakciji sa rezervoarom u analizama sloma od praktičnog interesa.

Predloženi numerički model je najprije predstavljen kao dvodimenzionalan, a zatim proširen dodavanjem treće dimenzije. Konstrukcija se posmatra kao zasićena porozna sredina izgra dena od kohezivnog materijala. Pretpostavlja se da vanjski fluid u interakciji sa konstrukcijom djeluje kao izvor poro-zasićenja. Odgovor konstrukcije od zasićenog materijala je modeliran diskretnim rešetkastim modelom koji je baziran na idealizaciji domene Voronoi ćelijama, sa nelinearnim Timoshenko-vim grednim konačnim elementima s obogaćenim poljem pomjeranja u vidu ugra denih jakih diskontinuiteta kao kohezivnim vezama. Interakcija izme du čvrste faze i fluida u porama materijala je u modelu definisana sa Biot-ovom teorijom i Darcy-evim zakonom tečenja fluida. Numerička implementacija problema unutrašnje interakcije rezultuje sa dodatnim stepenom slobode kretanja -porni pritisak po svakom čvoru konačnog elementa Timoshenkove grede, koji se kasnije koristi za rješavanje problema sučelja izme du konstrukcije i fluida.

Malo pomjeranje vanjskog fluida u rezervoaru omogućuje modeliranje vanjskog fluida sa akustičnom teorijom. Za numeričko modeliranje malih pomjeranja vanjskog fluida, izabrana je Lagrange-ova formulacija i mješovita pomjeranje/pritisak aproksimacija. Problem zajedničkog sučelja konstrukcije i fluida u numeričkom modelu interakcije je jednostavno riješen direktnim povezivnjem konačnih elemenata konstrukcije i konačnih elemenata vanjskog fluida koji dijele iste stepene slobode kretanja u zajedničkim čvorovima, što omogućava monolitan pristup rješavanju rezultujućeg sistema jednačina. Sva numerička implementacija i numerički proračuni su ura deni u istraživačkoj verziji kompjuterskog koda FEAP (Finite Element Analysis Program). Predloženi numerički modeli konstrukcije, vanjskog fluida kao i numerički model interakcije su validirani u linearno elastičnom području odgovora konstrukcije pore denjem izračunatih rezultata sa referentnim vrijednostima dobivenim analitičkim rješenjima ili na modelima kontinuuma. Numeričke simulacije u nelinearnom području odgovora konstrukcije su sprovedene s ciljem demonistriranja mogućnosti predloženog modela konstrukcije da predvidi kompletan odgovor i mehanizme loma u konstrukciji izgra denoj od zasićenog materijala. Konačno, testirana je mogućnost predloženog modela interakcije da opiše progresivni lokalizirani lom u konstrukciji brane izgra dene od kohezivnog zasićenog materijala u uslovima interakcije sa rezervoarom za odre deni program opterećenja. S ciljem analiziranja uticaja temperature, u numerički model konstrukcije je uvedeno temperaturno polje.

Ključne riječi: kohezivni diskretni gredni rešetkasti model, Timoshenko-va greda, ugra deni diskontinuitet, kuplovanje, pritisak, porozna sredina, Biot-ova teorija, akustična teorija, mješovita pomjeranje/pritisak aproksimacija, interakcija fluidkonstrukcija, analiza sloma, globalna sigurnost, temperatura xviii Introduction

Motivation

This research is motivated by the practical considerations of localized failure in structures built of heterogeneous and pore-saturated materials under extreme loads. Namely, the presence of the pore fluid can determine the final response and in many cases increase the risk of the ultimate failure of the pore-saturated structure. Examples of this happening can be found widely in everyday engineering practice, whether it is a landslide triggered due to heavy rain in a slope that was once stable, liquefaction phenomena occurring as a result of the pore fluid presence in combination with extreme dynamic loading, fracturing and collapse of the dams resulting from the combined influences of poor monitoring and accidental loads, or issue of sound design of floating wind turbines (Figure 1.1).

(a) Landslide on Cecil Lake Road The majority of dams in Bosnia and Herzegovina are built in the pre-war period, at least 30 years ago with most of them still operating at full capacity. Due to the poor maintenance and monitoring, the issue of the overall safety and durability of these structures is raised. Even in the developed countries where the maintenance and monitoring of dam structures are obligatory, the accidents that can pose a threat to structural integrity can occur. One recent example is the dam in Oroville, California, where in the spillway of the dam one large macro-crack formed suddenly due to severe erosion, raising a fear that a whole structure might ultimately collapse. In underdeveloped countries, this issue is even more pronounced. The recent failure of the Patel dam in Kenya, resulting with a death toll of 41 and hundreds of people forced to leave a completely devastated area, raised the question about the current state of the infrastructure in the country. For inspecting the current state and potential risk of failure of pore-saturated structures, aside from experimental results and in-situ data measurements, numerical models and numerical simulations can serve as a powerful tool. To achieve the goal of providing an overall safety assessment of pore-saturated structures under extreme loads, the numerical model has to provide a realistic description of the nonlinear response of the structure. Here, one of the most important features of the model is its ability to predict the true cracking pattern for a particular loading program, where cracks coalescence results with the ultimate failure mode of the structure.

Aims and scopes

The main scientific aim of this doctoral thesis is to provide a numerical model capable of performing an overall safety assessment of structures built of pore-saturated materials under extreme loads, in application to fluid-structure interaction problems such as dam-reservoir interaction. Namely, our goal is to provide a numerical model for investigating the current state of dam structures, in order to be able to determine the potential risk of failure or severe cracking that can endanger structural integrity, durability, and functionality.

When modeling fluid-structure interaction problems such as dam-reservoir interaction, we have to consider three parts where numerical modeling of each part can be identified as one research aim. The first is the numerical modeling of the structure. The majority of dam structures are built of porous cohesive materials, such as soil, rock, or concrete. The external fluid in interaction with structure acts as a source of external loading and also as a source of pore saturation. The presence of the pore fluid in heterogeneous and porous materials in many cases increases the risk of the final failure. The latter suggests that the numerical model of structure ought to account for the internal fluid-structure interaction. In other words, the structure in fluid-structure interaction ought to be represented in terms of a saturated porous medium. In order to perform the safety analysis, the chosen model of the pore-saturated structure has to be able to capture the true macro-scale response, along with the fracture process zone, and localized failure. In addition, a numerical model has to give an efficient performance in terms of fast convergence rates and efficient computational time.

The second part in numerical modeling of fluid-structure interaction is the numerical modeling of the external fluid. Numerical representation of the external fluid motion has to take into account the confined conditions typically met for the external fluid placed in a reservoir. The third and final part is the issue of the

Literature review

fluid-structure interface. In other words, how to connect the structure and the external fluid finite elements at the common boundary, which affects the manner (monolithic or partitioned) in which the resulting equations of the fluid-structure system can be solved.

Literature review

The fluid-structure interaction problems are frequently encountered in engineering practice. These problems range from one extreme where the external fluid is in large overall motion, to another extreme where fluid displacements can be regarded as very small, and the fluid behavior can be described with the acoustic wave theory.

In the latter group of acoustic fluid-structure interaction problems fall engineering structures like dams, reservoirs, containers or storage tanks. Concrete, rock or soil (for earth dams) are typically used for constructing many such structures. Both concrete, and especially soil are porous materials in which the presence of the pore fluid influences the final response and in many cases increases the risk of ultimate failure. For example, heavy rain can trigger a landslide in the slope that was once stable due to the reduction in the shear strength of the soil. For rocks, fluid flow can occur through the network of initial discontinuities and thus contribute to the rocks collapse or rock-falls.

In fluid-structure interaction problems, the external fluid acts as a source of loading on the structure and also as a source of pore saturation, keeping the material of the structure fully saturated in every time step. In already damage structure, fluid flow can also occur through cracks which have formed in the structure subjected to extreme loading. In this case, to fully describe the structure part in fluidstructure interaction problems, we ought to have a better numerical formulation of the structure response in terms of a saturated porous medium. In other words, we ought to account for the interaction between the solid phase and pore (internal) fluid and resulting coupling effects.

The pioneering works in the domain of solid phase-pore fluid interaction are Terzaghi's theory of one-dimensional consolidation [1], and Biot's theory of threedimensional consolidation [2]. Both theories are limited to the case of linear elastic behavior of the material. Biot's porous media theory is further extended to take into account nonlinearities of the material in the pre-peak part of the response [3][4][5]. Biot's porous media theory has also been extensively used as one of the main ingredients in numerical modeling of response, fracturing, and localized failure in saturated and partially saturated porous media. The successful numerical implementations can be found in works which exploit Extended Finite Element Method (XFEM) [START_REF] Réthoré | A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks[END_REF][START_REF] Mohammadnejad | Hydro-mechanical modeling of cohesive crack propagation in multiphase porous media using the extended finite element method[END_REF][START_REF] Cao | Interaction between crack tip advancement and fluid flow in fracturing saturated porous media[END_REF][START_REF] Remij | An investigation of the step-wise propagation of a mode-ii fracture in a poroelastic medium[END_REF][START_REF] Vahab | X-fem modeling of multizone hydraulic fracturing treatments within saturated porous media[END_REF], Partition-of-Unity Finite Element Method (PUFEM) [START_REF] Borst | A numerical approach for arbitrary cracks in a fluid-saturated medium[END_REF][START_REF] Kraaijeveld | Twodimensional mode i crack propagation in saturated ionized porous media using partition of unity finite elements[END_REF],

Embedded Discontinuity Finite Element Method (EDFEM) [START_REF] Armero | An analysis of strong discontinuities in a saturated poro-plastic solid[END_REF][START_REF] Callari | Finite element methods for the analysis of strong discontinuities in coupled poro-plastic media[END_REF][START_REF] Callari | Analysis and numerical simulation of strong discontinuities in finite strain poroplasticity[END_REF][START_REF] Callari | Strong discontinuities in partially saturated poroplastic solids[END_REF][START_REF] Nguyen | Modelling hydraulic fractures in porous media using flow cohesive interface elements[END_REF], Finite element methods with adaptive remeshing techniques [START_REF] Cao | Interaction between crack tip advancement and fluid flow in fracturing saturated porous media[END_REF][START_REF] Schrefler | On adaptive refinement techniques in multi-field problems including cohesive fracture[END_REF][START_REF] Secchi | Mesh adaptation and transfer schemes for discrete fracture propagation in porous materials[END_REF][START_REF] Secchi | A method for 3-d hydraulic fracturing simulation[END_REF][START_REF] Cao | Porous media fracturing dynamics: stepwise crack advancement and fluid pressure oscillations[END_REF], Central force lattice model [START_REF] Cao | Interaction between crack tip advancement and fluid flow in fracturing saturated porous media[END_REF][START_REF] Milanese | Avalanches in dry and saturated disordered media at fracture[END_REF][START_REF] Milanese | Avalanches in dry and saturated disordered media at fracture in shear and mixed mode scenarios[END_REF], Phase field model [START_REF] Wheeler | An augmented-lagrangian method for the phase-field approach for pressurized fractures[END_REF][START_REF] Mikelic | A phase-field method for propagating fluid-filled fractures coupled to a surrounding porous medium[END_REF][START_REF] Mikelić | Phase-field modeling of a fluid-driven fracture in a poroelastic medium[END_REF][START_REF] Lee | Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model[END_REF] or Discrete lattice model [START_REF] Bolander | Simulation of shrinkage induced cracking in cement composite overlays[END_REF][START_REF] Grassl | A lattice approach to model flow in cracked concrete[END_REF][START_REF] Šavija | Lattice modeling of chloride diffusion in sound and cracked concrete[END_REF][START_REF] Šavija | Lattice modeling of rapid chloride migration in concrete[END_REF][START_REF] Nikolic | Discrete element model for the analysis of fluid-saturated fractured poro-plastic medium based on sharp crack representation with embedded strong discontinuities[END_REF]. The solid phase-pore fluid interaction is further extended to account for the temperature effects in the saturated porous media. The pioneering work in the finite element analysis of thermo-hydro-mechanical coupling problem is the work of Aboustit et al. [START_REF] Aboustit | Variational principles and finite element simulations for thermo-elastic consolidation[END_REF][START_REF] Aboustit | Finite element invesitagtions of thermo-elastic and thermoplastic consolidation[END_REF], followed by the works of [START_REF] Britto | Finite element analysis of the coupled heat flow and consolidation around hot buried objects[END_REF][START_REF] Lewis | A coupled finite element model for the consolidation of nonisothermal elastoplastic porous media[END_REF][START_REF] Lewis | Finite element modelling of twophase heat and fluid flow in deforming porous media[END_REF][START_REF] Noorishad | Coupled thermohydroelasticity phenomena in variably saturated fractured porous rocks-formulation and numerical solution[END_REF][START_REF] Gatmiri | A formulation of fully coupled thermal-hydraulicmechanical behaviour of saturated porous media-numerical approach[END_REF][START_REF] Cui | An alternative coupled thermo-hydro-mechanical finite element formulation for fully saturated soils[END_REF].

The motion of the fluid is, in a general case, described with Navier-Stokes equations [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. However, the small motion of the fluid in contained conditions, such as reservoirs, water containers or storage tanks, allows us to derive the equations governing the external fluid motion from the acoustic wave theory [START_REF] Ohayon | Structural Acoustics and Vibration[END_REF][START_REF] Fahy | Sound and Structural Vibration: Radiation, Transmission and Response[END_REF]. In numerical modeling of acoustic fluids by the finite element method, two common approaches are Eulerian and Lagrangian.

In the Eulerian approach, the fluid motion is formulated in terms of fluid pressures, displacement potentials, or velocity potentials as state variables. The fluid finite element formulations falling into this category are the pressure formulation [START_REF] Zienkiewicz | Fluid-structure dynamic interaction and wave forces. an introduction to numerical treatment[END_REF], the displacement potential-pressure formulation [START_REF] Morand | Variational formulations for the elasto acoustic vibration problem-finite element results[END_REF][START_REF] Morand | Substructure variational analysis of the vibrations of coupled fluid-structure systems. finite element results[END_REF][START_REF] Morand | Fluid-Structure Interaction: Applied Numerical Methods[END_REF], the velocity potential formulation [START_REF] Everstine | A symmetric potential formulation for fluid-structure interaction[END_REF][START_REF] Olson | Analysis of fluid-structure interactions. a direct symmetric coupled formulation based on the fluid velocity potential[END_REF] and the velocity potential-density formulation [START_REF] Kock | Fluid-structure interaction analysis by the finite element method-a variational approach[END_REF]. The Eulerian approach is commonly used in solving fluid-structure systems, in which the behavior of the structure is assumed to be linear or nonlinear [START_REF] Saini | Coupled hydrodynamic response of concrete gravity dams using finite and infinite elements[END_REF][START_REF] Cervera | Seismic evaluation of concrete dams via continuum damage models[END_REF][START_REF] Calayir | Earthquake analysis of gravity dam-reservoir systems using the eulerian and lagrangian approaches[END_REF][START_REF] Everstine | Finite element formulations of structural acoustics problems[END_REF][START_REF] Küçükarslan | Transient analysis of damreservoir interaction including the reservoir bottom effects[END_REF][START_REF] Gogoi | Influence of sediment layers on dynamic behavior of aged concrete dams[END_REF][START_REF] Wang | Dynamic analysis of a water-soil-pore water coupling system[END_REF][START_REF] Bouaanani | Assessment of potential-based fluid finite elements for seismic analysis of dam-reservoir systems[END_REF][START_REF] Mitra | 2d simulation of fluid-structure interaction using finite element method[END_REF][START_REF] Mandal | 2d finite element analysis of rectangular water tank with separator wall using direct coupling[END_REF]. However, both the equilibrium condition and the compatibility condition have to be enforced at the fluid-structure interface, because of different state variables in the fluid and the structure.

In the Lagrangian approach, the fluid motion is formulated in terms of nodal displacements [START_REF] Zienkiewicz | Fluid-structure dynamic interaction and wave forces. an introduction to numerical treatment[END_REF][START_REF] Chopra | Earthquake analysis of reservoir-dam systems[END_REF][START_REF] Shugar | Development of finite element head injury model[END_REF][START_REF] Kiefling | Fluid-structure finite element vibrational analysis[END_REF][START_REF] Shantaram | Dynamic transient behaviour of two-and three-dimensional structures including plasticity, large deformation effects and fluid interaction[END_REF][START_REF] Belytschko | A fluid-structure finite element method for the analysis of reactor safety problems[END_REF][START_REF] Khalil | Parametric study of head response by finite element modeling[END_REF][START_REF] Wilson | Finite elements for foundations, joints and fluids[END_REF][START_REF] Bathe | On transient analysis of fluid-structure systems[END_REF]. The advantage of displacement based formulation is in its finite element implementation because fluid finite elements share the same degrees of freedom as the standard finite elements commonly used in the numerical modeling of the structure. Hence, equilibrium and compatibility conditions are automatically satisfied through element assembly procedure. However, it has been reported that the pure displacement based fluid finite element formulation suffers from the existence of spurious zero-energy deformation modes [START_REF] Kiefling | Fluid-structure finite element vibrational analysis[END_REF][START_REF] Hamdi | A displacement method for the analysis of vibrations of coupled fluid-structure systems[END_REF][START_REF] Akkas | Applicability of general-purpose finite element programs in solid-fluid interaction problems[END_REF][START_REF] Olson | A study of displacement-based fluid finite elements for calculating frequencies of fluid and fluid-structure systems[END_REF][START_REF] Chen | Vibration analysis of fluid-solid systems using a finite element displacement formulation[END_REF].

Various numerical treatments have been proposed with the aim to improve the proposed formulation, such as the penalty method [START_REF] Hamdi | A displacement method for the analysis of vibrations of coupled fluid-structure systems[END_REF], the finite element formulation with rotational constraints and a reduced integration technique [START_REF] Wilson | Finite elements for the dynamic analysis of fluidsolid systems[END_REF], and the finite element formulation with a reduced integration technique and element mass projection [START_REF] Chen | Vibration analysis of fluid-solid systems using a finite element displacement formulation[END_REF]. With the aim to solve the issue of the existence of spurious modes, the mixed displacement/pressure based fluid finite element formulation is proposed, with nodal displacements and pressures as unknown variables [START_REF] Bathe | A mixed displacement-based finite element formulation for acoustic fluid-structure interaction[END_REF][START_REF] Wang | Displacement/pressure based mixed finite element formulations for acoustic fluid-structure interaction problems[END_REF][START_REF] Wang | On mixed elements for acoustic fluid-structure interactions[END_REF].

The Lagrangian approach with displacement variables in both the structure and the fluid is employed in many numerical models of fluid-structure interaction with linear or nonlinear response of structure, such as [START_REF] Calayir | Earthquake analysis of gravity dam-reservoir systems using the eulerian and lagrangian approaches[END_REF][START_REF] Calayir | Static and dynamic analysis of fluid and fluid-structure systems by the lagrangian method[END_REF][START_REF] Calayir | A continuum damage concrete model for earthquake analysis of concrete gravity dam-reservoir systems[END_REF][START_REF] Calayir | Seismic fracture analysis of concrete gravity dams including dam-reservoir interaction[END_REF][START_REF] Akköse | Elasto-plastic earthquake response of arch dams including fluid-structure interaction by the lagrangian approach[END_REF][START_REF] Bilici | Stochastic dynamic response of dam-reservoir-foundation systems to spatially varying earthquake ground motions[END_REF][START_REF] Pelecanos | Numerical modelling of hydrodynamic pressures on dams[END_REF][START_REF] Kartal | Non-linear earthquake response of cfr damreservoir-foundation systems[END_REF][START_REF] Pelecanos | The effects of dam-reservoir interaction on the nonlinear seismic response of earth dams[END_REF].

Beside the Eulerian and Lagrangian approach, an added mass or Westergaard approach is also commonly used for solving fluid-structure systems. The added mass approach, which follows the work of Westergaard [START_REF] Westergaard | Water pressures on dams during earthquakes[END_REF], is the first and the most simple approach for solving fluid-structure interaction. In this approach, the fluid-structure interaction is solved by adding the fraction of the fluid mass to the structural model. The added mass is typically evaluated based on the assumption of the incompressible fluid. The numerical models of fluid-structure interaction, which exploits an added mass approach can be found in [START_REF] Kuo | Fluid-structure interactions: Added mass computations for incompressible fluid[END_REF][START_REF] Bhattacharjee | Seismic cracking and energy dissipation in concrete gravity dams[END_REF][START_REF] Ghrib | An application of damage mechanics for seismic analysis of concrete gravity dams[END_REF][START_REF] Lee | A plastic-damage concrete model for earthquake analysis of dams[END_REF][START_REF] Ghaemian | Nonlinear seismic response of concrete gravity dams with dam-reservoir interaction[END_REF][START_REF] Calayir | A continuum damage concrete model for earthquake analysis of concrete gravity dam-reservoir systems[END_REF].

Thesis overview

In this thesis, we propose a numerical model capable of dealing with nonlinear fluid-structure interaction problems, such as progressive localized failure of a dam structure built of porous cohesive material under dam-reservoir interaction. The numerical model of the structure, the numerical model of external fluid and ultimately the numerical model of interaction are first presented in two-dimensional (2D) framework and later extended to three-dimensional (3D) framework.

We first start with the numerical model of the structure. Here, we assume that the structure is built of porous cohesive material and that the external fluid in interaction with structure acts as a source of pore saturation. developed by R.L. Taylor [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. The finite element mesh in numerical examples is generated in GMSH [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF], with cross-sectional dimensions of beam finite elements computed by using MATLAB [START_REF]version 8.5.0.197613[END_REF].

Introduction

In this Chapter, we describe the plane strain coupled discrete beam lattice model used in this thesis for numerical modeling of response and localized failure in structures built of saturated porous cohesive materials. The theoretical formulation for the pore-saturated structure is based upon Biot's porous media theory extended to localized poro-plasticity. The internal fluid-structure interaction, i.e., the interaction between the solid phase and the pore fluid is handled by using a judicious combination of lattice models for structure cracking and finite element approximation for internal (pore) fluid flow.

The outline of the Chapter is as follows: In Section 2.2, we briefly describe Biot's porous media theory. In Section 2.3, we provide a detailed description of discrete lattice models and their application in numerical modeling of response and localized failure in structures. We also explain in detail the manner in which the internal fluid-structure interaction is numerically treated. In Section 2.4, we present a finite element formulation of the proposed coupled discrete lattice model, together with the computational procedure. The finite element formulation is given in fully dynamics framework, which can be easily adjusted for quasi-static simplification of structure built of a saturated porous medium, or dry medium. In Section 2.5, we present results of several numerical simulations with the aim to validate the proposed coupled discrete beam lattice model and to demonstrate its capabilities to capture the response and localized failure in pore-saturated structures, where the cracks coalescence results with the ultimate failure modes. In Section 2.6, we give concluding remarks.

Biot's porous media theory

In Biot's porous media theory, the equations governing the coupled problem are derived by combining equilibrium and continuity equations [2]. In a dynamic setting, we consider the equation of motion instead of the equilibrium equation.

In deriving the strong form of the equations of motion, we use the d'Alembert principle. The main idea behind the d'Alembert principle is that a motion of a solid body in each time step can be described by the equilibrium equations, which includes an additional external load in terms of the inertia force. Thus, the strong form of equations of motion is written as

∇ • σ + b -ρa = 0 (2.1)
where σ is the total stress tensor, b are the body forces vector, a is the solid phase acceleration vector and ρ is the mass density of the mixture assumed to be constant, which is defined as ρ = (1-n)ρ s +nρ f , with n as the porosity, ρ s as the mass density of solid phase and ρ f as the mass density of pore fluid.

In the formulation of the coupled problem, Biot's theory exploits Terzaghi's principle of effective stresses, which states that the total normal stress is equal to the sum of the effective stress carried by the solid phase, and the pore pressure carried by the fluid in pores. Terzaghi's principle of effective stresses is written as

σ = σ ′ -Ibp (2.2)
where σ ′ is the effective stress tensor, I is the second order identity tensor, p is the pore pressure assumed positive in compression, and b is Biot's constant defined as b = 1 -K t /K s , with K t and K s as the bulk moduli of the porous skeleton and the solid phase, respectively. If the material grains are assumed to be incompressible, then Biot's constant is equal to b = 1.

The continuity equation for pore fluid flow is written as

ζ + ∇ • q = 0 (2.3)
where ζ is the variation of the fluid content, and q is the fluid flux.

The variation of the fluid content ζ is written as

ζ = 1 M p + b∇ • u (2.4)
where u is the displacement of the solid phase, and M is Biot's modulus, defined as 1/M = (bn)/K s + n/K f , with K f as the bulk modulus of the pore fluid. If the material grains and the pore fluid are assumed to be incompressible, then Biot's

modulus is M → ∞ (1/M → 0).
According to Darcy's law, the fluid flux q is equal to

q = - k γ f ∇p (2.5)
where k is the coefficient of permeability of the isotropic porous medium, and γ f is the specific weight of the fluid.

By combining Equations (2.3-2.5), a new form of continuity equation can be written as 1

M ṗ + b∇ • u -∇ • k γ f ∇p = 0 (2.6)
In everyday engineering practice, the response of the structure built of heterogeneous materials is usually computed with homogenized, continuum type of numerical models with strength and deformability defined at the macro-scale, in which the domain is meshed with two-dimensional (2D) or three-dimensional (3D) finite elements. The main advantage of these models is in a relatively simple and efficient way of predicting the main tendencies and overall behavior of heterogeneous materials in failure analyses of interest for engineering practice. However, the use of 2D or 3D finite elements can be challenging when trying to numerically represent the inelastic behavior, and especially cracking and localized failure in the structure built of heterogeneous materials. Successful numerical implementations of strain localization phenomenon can be found in [START_REF] Armero | An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids[END_REF][START_REF] Oliver | Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. part 2: Numerical simulation[END_REF][START_REF] Ibrahimbegovic | Combined hardening and softening constitutive model of plasticity: precursor to shear slip line failure[END_REF][START_REF] Linder | Finite elements with embedded strong discontinuities for the modeling of failure in solids[END_REF][START_REF] Brancherie | Novel anisotropic continuum-discrete damage model capable of representing localized failure of massive structures: Part i: theoretical formulation and numerical implementation[END_REF][START_REF] Hautefeuille | Failure model for heterogeneous structures using structured meshes and accounting for probability aspects[END_REF][START_REF] Dujc | Quadrilateral finite element with embedded strong discontinuity for failure analysis of solids[END_REF][START_REF] Dujc | Stress-hybrid quadrilateral finite element with embedded strong discontinuity for failure analysis of plane stress solids[END_REF][START_REF] Saksala | Anisotropic viscodamage-viscoplastic consistency constitutive model with a parabolic cap for rocks with brittle and ductile behaviour[END_REF][START_REF] Saksala | Combined continuum damage-embedded discontinuity model for explicit dynamic fracture analyses of quasi-brittle materials[END_REF][START_REF] Do | Dynamics framework for 2d anisotropic continuum-discrete damage model for progressive localized failure of massive structures[END_REF][START_REF] Benedetti | 3d numerical modelling of twisting cracks under bending and torsion of skew notched beams[END_REF].

Contrary to the continuum type of numerical models, in discrete lattice models the domain is meshed with one-dimensional (1D) finite elements. The use of 1D, instead of 2D or 3D elements, has its advantages in terms of numerical implementation of different kind of phenomena commonly observed in heterogeneous materials, improving rates of convergence, and decreasing computational time. whether it is based on spring or beam lattice network. More details can be found in works [START_REF] Nikolić | Lattice element models and their peculiarities[END_REF][START_REF] Schlangen | Fracture simulations of concrete using lattice models: computational aspects[END_REF][START_REF] Ostoja-Starzewski | Lattice models in micromechanics[END_REF][START_REF] Karihaloo | Lattice modelling of the failure of particle composites[END_REF][START_REF] Bolander | Irregular lattice model for quasistatic crack propagation[END_REF][START_REF] Berton | Crack band model of fracture in irregular lattices[END_REF].

The successful application of discrete lattice models in numerical modeling of inelastic response and localized failure in materials rests on the idea of multiscale modeling. Namely, by appropriate numerical modeling of the micro-scale represented by cohesive links, we can obtain the true macro-scale response of the structure. The discrete lattice models with geometrically exact shear deformable beams acting as cohesive links have been used for modeling dynamic fracture of structures built of brittle material [START_REF] Ibrahimbegovic | Microscale and mesoscale discrete models for dynamic fracture of structures built of brittle material[END_REF]. In this discrete model, fracture criterion is stated in terms of strains with two possible modes of fracture, the separation between two adjacent cells and the relative bending deformation. The fracture limits for both modes are assumed to follow Gaussian random distribution. The discrete lattice models with truss bar elements capable of representing the response of a two-phase material [START_REF] Ibrahimbegovic | Embedded discontinuity finite element method for modeling of localized failure in heterogeneous materials with structured mesh: an alternative to extended finite element method[END_REF], have been used in numerical modeling of the response and fracturing in concrete [START_REF] Benkemoun | Failure of heterogeneous materials: 3d meso-scale fe models with embedded discontinuities[END_REF][START_REF] Benkemoun | Anisotropic constitutive model of plasticity capable of accounting for details of meso-structure of two-phase composite material[END_REF]. The truss bar elements enhanced with additional kinematics in terms of embedded strong discontinuity were able to model the crack formation and propagation in mode I, that relates to crack opening.

Beside mode I, the crack can also propagate in mode II that relates to crack sliding.

To be able to capture both mode I and mode II of failure, Timoshenko beam finite elements with an embedded strong discontinuity in both axial and transverse direction have been proposed instead of truss bar elements. The discrete lattice models with enhanced Timoshenko beam finite elements have been successfully used to describe the response and fracturing in rocks, both in 2D and 3D setting [START_REF] Nikolic | Brittle and ductile failure of rocks: embedded discontinuity approach for representing mode i and mode ii failure mechanisms[END_REF][START_REF] Cvitanovic | Influence of specimen shape deviations on uniaxial compressive strength of limestone and similar rocks[END_REF][START_REF] Nikolic | Rock mechanics model capable of representing initial heterogeneities and full set of 3d failure mechanisms[END_REF][START_REF] Nikolic | Rock mechanics, failure phenomena with pre-existing cracks and internal fluid flow through cracks[END_REF].

In the discrete lattice model proposed in [START_REF] Nikolic | Brittle and ductile failure of rocks: embedded discontinuity approach for representing mode i and mode ii failure mechanisms[END_REF], a distinction between two phases has been made. Phase I is an intact rock which is not likely to fail. This phase is represented with standard linear elastic Timoshenko beam finite element. Phase II represents a weak phase or initial discontinuities in rocks. Inelastic Timoshenko beam finite elements with enhanced kinematics are used to model this phase.

Our final goal in this thesis is to tackle the acoustic fluid-structure interaction problems. In the acoustic fluid-structure interaction problems on hand, the external fluid acts both as a source of the loading on the structure, and also as a source of the pore saturation. For a better numerical description of the underlying phenomena, we ought to model the structure as a saturated porous medium. In this thesis, for the numerical representation of the pore-saturated structures, we generalize the development of the coupled discrete beam lattice model presented in [START_REF] Nikolic | Discrete element model for the analysis of fluid-saturated fractured poro-plastic medium based on sharp crack representation with embedded strong discontinuities[END_REF][START_REF] Nikolic | Rock mechanics, failure phenomena with pre-existing cracks and internal fluid flow through cracks[END_REF],

extending its ability to deal with the nonlinear response and localized failure in pore-saturated structures in acoustic fluid-structure interaction problems of our interest. The proposed coupled discrete beam lattice model is an extension of the discrete beam lattice model with enhanced Timoshenko beam finite elements as cohesive links presented in [START_REF] Nikolic | Brittle and ductile failure of rocks: embedded discontinuity approach for representing mode i and mode ii failure mechanisms[END_REF][START_REF] Nikolic | Rock mechanics, failure phenomena with pre-existing cracks and internal fluid flow through cracks[END_REF], with Biot's porous media theory governing the coupling between the solid phase and pore fluid as a new ingredient. Different types of failure mechanisms in different materials are, among others, the result of material heterogeneities, loading, and boundary conditions. Because of this, in the coupled discrete beam lattice model used in this thesis, we allow for all Timoshenko beams to fail [START_REF] Hadzalic | Failure mechanisms in coupled soil-foundation systems[END_REF][START_REF] Hadzalic | Failure mechanisms in coupled poro-plastic medium[END_REF][START_REF] Hadzalic | Macroscale response of dry and saturated soil subjected to footing load[END_REF].

To the best of author's knowledge, the literature search did not reveal any confirmed relations between elastic properties of the continuum model and Timoshenko beam lattice in discrete models based on Voronoi cell discretization. However, from the condition that no stiffness is gained or lost compared to the contin-

Finite element formulation

The starting point for our work on coupled discrete beam lattice model of poresaturated structure is the Timoshenko beam finite element presented in [START_REF] Nikolic | Brittle and ductile failure of rocks: embedded discontinuity approach for representing mode i and mode ii failure mechanisms[END_REF][START_REF] Nikolic | Discrete element model for the analysis of fluid-saturated fractured poro-plastic medium based on sharp crack representation with embedded strong discontinuities[END_REF][START_REF] Nikolic | Rock mechanics, failure phenomena with pre-existing cracks and internal fluid flow through cracks[END_REF]. The detailed descriptions of various numerical procedures and numerical solutions used in the finite element formulation of the Timoshenko beam finite element on hand can be found in [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF].

Kinematics

Consider a straight plane Timoshenko beam finite element of length L e and crosssectional area A e . The element has two nodes, and three degrees of freedom per node: axial displacement, transverse displacement, and rotation of cross-section.

Displacement fields are interpolated with standard linear interpolation functions, written as

u = N s u ū (2.7)
where

u T = u, v, θ ; ūT = u 1 , v 1 , θ 1 , u 2 , v 2 , θ 2 ; N s u =     N 1 0 0 N 2 0 0 0 N 1 0 0 N 2 0 0 0 N 1 0 0 N 2     ; N 1 , N 2 = 1 -x L e , x L e
(2.8)

Finite element formulation

For simplicity, we considered the shape functions for a Timoshenko beam placed along the global coordinate x axis, which can be easily adapted to any arbitrary element orientation by using a local element frame.

For geometrically linear Timoshenko beam, strain fields are defined as

ε xx = d u d x ; γ x y = d v d x -θ; κ z = d θ d x (2.9)
where ε xx is the axial strain, γ x y is the shear strain and κ z is the curvature. In what is to follow, we will denote these strains with the following symbols ε, γ and κ.

Following the standard finite element procedure, strain fields are interpolated as

ϵ = B s u ū (2.10)
where

ϵ T = ε, γ, κ ; B s u =     B 1 0 0 B 2 0 0 0 B 1 -N 1 0 B 2 -N 2 0 0 B 1 0 0 B 2     ; B 1 , B 2 = d N 1 d x , d N 2 d x = -1 L e , 1 L e (2.11)
For modeling crack formation in mode I and mode II, i.e., the softening behavior, localization limiter based on embedded strong discontinuity is employed [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF][START_REF] Jukić | Embedded discontinuity finite element formulation for failure analysis of planar reinforced concrete beams and frames[END_REF][START_REF] Jukić | Failure analysis of reinforced concrete frames by beam finite element that combines damage, plasticity and embedded discontinuity[END_REF].

Namely, strong discontinuities are introduced in displacement fields in the axial and transverse direction (Figure 2.5), which are now represented as a sum of a regular part and the discontinuous part. The discontinuous part is represented as a product of an interpolation function M and parameter α, which represents the displacement jump located at the center of the element.

The enhanced displacement fields are interpolated as

u = N s u ū + Mα (2.12)
where

α T = α u , α v , 0 ; M =     M 0 0 0 M 0 0 0 0     ; M = H x -N 2 ; H x =    0, x ≤ x 1, x > x (2.13)
Here, α u and α v represent displacement jumps in axial and transverse direction. The finite element implementation procedure for the localization limiter based on embedded strong discontinuity fits well within the framework of incompatible mode method [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]. Namely, the parameter α can be interpreted as an incompatible mode parameter, which can be statically condensed on the element level, so that on the global level only remaining unknowns are nodal displacements.

δ x N 1 ( ) N 2 ( ) M( ) G(x) H x (x) u 1 v 1 θ 1 1 v 2 θ 2 u 2 2 α v α u L e x y x
The enhanced strain fields are interpolated as

ϵ = B s u ū + Gα (2.14)
where

G =     G 0 0 0 G 0 0 0 0     ; G = d M d x =    Ḡ, x ∈ [0, x) ∪ x ∈ ( x, L e ] Ḡ + δ x , x = x (2.15)
where Ḡ = -1 L e , and δ x is the Dirac function

δ x =    0, x ∈ [0, x) ∪ x ∈ ( x, L e ] ∞, x = x (2.16)
The pore pressure field is approximated with CST finite elements (Figure 2.6).

The finite element approximation for the pore pressure field is written as p = N s p p (2.17) where pT = p 1 , p 2 , p 3 ;

N s p = N p 1 , N p 2 , N p 3 ; N p 1 = 1 2A (x 2 y 3 -x 3 y 2 ) + (y 2 -y 3 )x + (x 3 -x 2 )y] N p 2 = 1 2A (x 3 y 1 -x 1 y 3 ) + (y 3 -x 1 )x + (x 1 -x 3 )y] N p 3 = 1 2A (x 1 y 2 -x 2 y 1 ) + (y 1 -x 2 )x + (x 2 -y 1 )y] (2.18)
where A is the area of the triangular element, x,y, are global coordinates and x i , y i are nodal coordinates of CST element.

Continuity equation and equations of motion 2.4.2.1 Continuity equation

The coupling of the mechanics and the pore fluid flow occurs through the axial direction of the Timoshenko beam finite element. The continuity equation for fluid flow through the porous body is written as

1 M ṗ + b ε - k γ f ∇ • (∇p) = 0 (2.20)
The weak from of the continuity equation is obtained through the principle of virtual works, written as

Ω e C ST δp 1 M ṗ - k γ f ∇ • (∇p) d Ω + L e 0 δp ‵ b εd x = 0 (2.21)
where δp is the virtual pore pressure field interpolated in the same manner as the real pore pressure field δp = N s p δ p; δ pT = δp 1 , δp 2 , δp 3 ;

p ‵ = N s up p‵ ; N s up = N 1 , N 2 ; p‵T = p 1 , p 2 ; δp ‵ = N s up δ p‵ ; δ p‵T = δp 1 , δp 2 (2.22)

Equations of motion

In deriving the weak form of the equations of motion, we use the d'Alembert principle, which is a dynamic analog to the principle of virtual work in the quasistatic setting [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]. The main idea behind the d'Alembert principle is that a motion of a solid body in each time step can be described by the equilibrium equation, which includes an additional external load in terms of the inertia force. Here, we assume that the coupling of the mechanics and the internal fluid flow occurs through the axial direction of the Timoshenko beam finite element. Thus, the strong form of the equations of motion for a Timoshenko beam is written as

-ρ A e ∂ 2 u ∂t 2 + ∂N ∂x + n(x, t ) = 0 (2.23) -ρ A e ∂ 2 v ∂t 2 + ∂V ∂x + q(x, t ) = 0 (2.24) -ρI e ∂ 2 θ ∂t 2 + ∂M ∂x + V + m(x, t ) = 0 (2.25)
where ρ is the mass density, A e and I e are the cross-sectional properties of Timoshenko beam: area and second moment of inertia of a cross-section, n, q and m are the distributed external loads, and N = N ′bp ‵ A e , V = V ′ and M = M ′ are the total axial force, the shear force and the bending moment, respectively. The superscript ′ denotes effective force.

The weak form of equations of motion for Timoshenko beam finite element is derived from the d'Alembert principle, which states that the virtual work of internal forces is equal to the virtual work of external forces, which include the Finite element formulation inertia force. The virtual displacements are imaginary displacements which are kinematically admissible and small enough so that virtual strains can be regarded as infinitesimal. The virtual displacement fields and virtual pore pressure field refer to a particular deformed configuration and thus are not time-dependent. The virtual displacement and virtual strain fields are interpolated in the same manner as the real displacement and real strain fields, written as

δu = N s u δ ū + Mδα (2.26) δϵ = B s u δ ū + Gδα (2.27)
where δu T = δu, δv, δθ ;

δα T = δα u , δα v , 0 ; δ ūT = δu 1 , δv 1 , δθ 1 , δu 2 , δv 2 , δθ 2 ;
δϵ T = δε, δγ, δκ

(2.28)
The virtual work of external forces is computed as The virtual work of internal forces for Timoshenko finite beam finite element is equal to

G ext ,e = δ ūT f ext ,e -
G i nt ,e = L e 0 B s u δ ū T σd x + L e 0 Gδα T σ u d x (2.31)
where

σ T = N ,V, M = N ′ -bp ‵ A e ,V ′ , M ′ ; σ T = σ T u -σ T p ; σ T u = N ′ ,V ′ , M ′ ; σ T p = bp ‵ A e , 0, 0 (2.32)

Discrete model of structure built of saturated poro-plastic medium

The virtual work of internal forces consists of two parts. The first is the virtual work of internal forces in the bulk part of the element (part of the element outside the discontinuity), and the second is the virtual work of internal forces acting at the discontinuity, written as

G i nt ,e = δ ūT f i nt ,e + δα T h e (2.33)
where f i nt ,e is the internal force vector, and h e is the residual vector due to discontinuity, computed as

f i nt ,e = L e 0 B s,T u σd x = f i nt ,e (σ u ) -f i nt ,e (σ p )
;

h e = L e 0 G T σ u d x = L e 0 ḠT σ u d x + t; Ḡ =     Ḡ 0 0 0 Ḡ 0 0 0 0     ; t T = t u , t v , 0 (2.34)
Here, t is the internal force vector acting at the discontinuity.

By exploiting the principle of virtual work we obtain 0 = G i nt ,e -G ext ,e = δ ūT (f i nt ,ef ext ,e + f acc,e ) + δα T h e (2.35)

The previous equation is satisfied if

f acc,e + f i nt ,e -f ext ,e = 0 (2.
)

h e = 0 (2.37)
The Equation (2.36) relates to the bulk part of the element, and the Equation (2.37)

relates to the discontinuity.

The internal force vector t acting at the discontinuity is computed from the following conditions

h e = 0 L e 0 Ḡd x = L e 0 - 1
L e d x = -1        ⇒ t =     1 0 0 0 1 0 0 0 0     L e 0 σ u d x (2.38)

Constitutive model

The pre-peak response of the Timoshenko beam finite element in both axial and transverse direction is described with the elasto-viscoplastic constitutive model with implemented linear hardening [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], and Fredrick-Armstrong nonlinear kinematic hardening law [START_REF] Armstrong | A mathematical representation of the multiaxial bauschinger effect[END_REF]. The post-peak response of the element is described with exponential softening. The behavior of the element in bending is purely linear elastic.

Plasticity model

The proposed model of plasticity is described with [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF] 1. Additive decomposition of the total strain into elastic and viscoplastic part, with only the elastic part influencing the value of the stresses

ε = ε e + ε v p γ = γ e + γ v p κ = κ e
(2.39)

2. Strain energy function in terms of strains and internal variables, plastic deformations ε v p , γ v p and strain-like hardening variables ξu , ξv

ψu ε, ε v p , ξu = 1 2 ε -ε v p E A e ε -ε v p + 1 2 ξu H u l h A e ξu ψv γ, γ v p , ξv = 1 2 γ -γ v p k c G A e γ -γ v p + 1 2 ξv H v l h A e ξv (2.40)
where E is Young's modulus, G = E 2(1+ν) is the shear modulus with ν as Poisson's ratio, H u l h , H v l h are the linear hardening moduli and k c is the shear correction factor. For rectangular cross-sections shear correction factor is equal to k c = 5/6.

3. Yield function in terms of stresses (or stress resultants), stress-like hardening variables qu , qv , and back-stress variables χu , χv

φu N ′ , χu , qu = |N ′ -χu A e | -N y -qu A e ≤ 0 φv V ′ , χv , qv = |V ′ -χv A e | -V y -qv A e ≤ 0 (2.41)
where N y , V y are the elastic (yield) limits.

The stress resultant values N ′ , V ′ and M ′ are computed from the elastic part of total deformation, as

N ′ = E A e ε -ε v p V ′ = k c G A e γ -γ v p M ′ = E I e κ (2.42)
where the area A e and second moment of inertia I e for rectangular crosssections with cross-sectional width equal to one and cross-sectional height equal to h e are A e = h e , I e = (h e ) 3 /12.

Remark: We note that, in a general case, the element local coordinate axis x does not pass through the center of gravity of the cross-section (Figure 2.

7).

In what is to follow, we will neglect the additional coupling between axial and bending strains as a result of the eccentricity e c . In experiments with cyclic loading, it has been observed that in a typical cycle with a load reversal, the plasticity threshold limit is reduced from the previous value. This phenomenon is called the Bauschinger effect. To track the change in the value of the plasticity threshold limit with respect to the previous plastic state, back-stress variable χ is introduced. We employ Fredrick-Armstrong nonlinear kinematic hardening law [START_REF] Armstrong | A mathematical representation of the multiaxial bauschinger effect[END_REF], which we can write as

C C e c

χu = H u l k εvp -H u nl k ξu χu χv = H v l k γvp -H v nl k ξv χv (2.44)
where H u l k , H v l k are the kinematic hardening moduli, and H u nl k , H v nl k are the non-dimensional material parameters that characterize the nonlinear kinematic hardening behavior.

3. The evolution equations for internal variables ϵ v p , γ v p and ξu , ξv obtained through the principle of maximum plastic dissipation and Kuhn-Tucker optimality conditions εvp = γu si g n(N ′ -χu A e ); ξu = γu ; γu ≥ 0 γvp = γv si g n(V ′ -χv A e ); ξv = γv ; γv ≥ 0 (2.45) where γu , γv are the plastic multipliers whose values are equal to

γu = 〈 φu (N ′ , χu , qu )〉 η u γv = 〈 φv (V ′ , χv , qv )〉 η v (2.46)
where η u , η v are the viscosity parameters.

Exponential softening

Once the element enters the softening phase, the jump in displacement is activated.

All plastic deformation from that point on remains localized at the discontinuity, whereas the bulk part of the element elastically unloads. The exponential softening is described with 1. The yield function defined in terms of stresses (or stress resultants), and dual variables qu , qv

φu t u , qu = |t u | -N f -qu A e ≤ 0 φv t v , qv = |t v | -V f -qv A e ≤ 0 (2.47)
where N f , V f are the fracture limits, and qu , qv are the stress-like softening variables that for exponential softening take the following form

qu = N f A e 1 -exp -ξu (N f /A e ) G u f qv = V f A e 1 -exp -ξv (V f /A e ) G v f (2.48)
where G u f , G v f are the fracture energies, and ξu , ξv are the strain-like softening variables. The stress resultant values in the bulk part of an element in the softening phase, which determine the value of internal forces at the discontinuity t u , t v are computed as

N ′ = E A e ε -ε v p + Ḡα u V ′ = k c G A e γ -γ v p + Ḡα v (2.49)
2. The evolution equations for internal variables α u , α v and ξu , ξv with the loading/unloading conditions obtained through the principle of maximum plastic dissipation and Kuhn-Tucker optimality conditions [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF] αu = γu si g n(t u ); ξu = γu ; γu ≥ 0; φu ≤ 0; γu φu = 0 αv = γv si g n(t v ); ξv = γv ; γv ≥ 0; φv ≤ 0; γv φv = 0

Computational procedure

The end result of the finite element discretization procedure is two sets of equations.

The first is the first order differential evolution equations for internal variables defined locally (on the element level) at each Gauss quadrature point, and the second is the second order differential equations in time governing the coupled problem defined globally. The solution of these equations is computed at discrete pseudotime steps t 1 , t 2 , ..., t n by using the operator split solution procedure [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]. Here, the solution procedure is divided into the local phase and the global phase, which are treated separately. In the local phase, the solution of evolution equations is computed by using implicit backward Euler time integration scheme. The unknown values of displacement jumps are statically condensed on the element level from the condition that the residual at the discontinuity is equal to zero. In the global phase, the solution in terms of the unknown nodal displacements and pore pressures is computed in a fully monolithic manner using Newmark time-integration scheme and Newton's iterative method.

Finite element formulation

Following the standard finite element discretization procedure, the system of equations governing the coupled problem for single Timoshenko beam finite element is written as

M e uu ü + f i nt ,e (σ u ) -K e up p‵ = f ext ,e (2.52
)

K e,T up u + D e pp ṗ‵ + K e pp p‵ = q ext ,e (2.53)
where M e uu is the mass matrix, f i nt ,e (σ u ) is the internal load vector resulting from displacements (i.e. effective stresses), K e up is the coupling matrix, D e pp is the compressibility matrix, K e pp is the permeability matrix, and f e,ext and q e,ext are the load vectors. The matrices K e up , D e pp and K e pp are computed as

K e up = L e 0 B s,T up bN s up d x; B s up = B 1 0 0 B 2 0 0 D e pp = Ω e C ST N s,T p 1 M N s p d Ω; K e pp = Ω e C ST (∇N s p ) T k γ f ∇N s p d Ω (2.54)
Here we note that in the global system of equations (Equations 2.52 and 2.53)

we take the parts of the D e pp , and K e pp matrices (Equation 2.54) that correspond to the nodes of Timoshenko beam finite element.

The Timoshenko beam finite element mass matrix M e in Equation (2.52) is obtained by distributing the total mass of an element to nodes (Figure 2.8), resulting in a diagonally lumped mass matrix, written as [START_REF] Ibrahimbegovic | Microscale and mesoscale discrete models for dynamic fracture of structures built of brittle material[END_REF]]

M e uu = 1
⇒            N ′ t r i al n+1 = E A e ε (i ) n+1 -ε v p n qu,tr i al n+1 = -H u l h ξu n χu,tr i al n+1 = H u l k ε v p n -H u nl k ξu n χu n ; ⇒            V ′ t r i al n+1 = k c G A e γ (i ) n+1 -γ v p n qv,tr i al n+1 = -H v l h ξv n χv,tr i al n+1 = H v l k γ v p n -H v
χu n + η u ∆t ⇒                  ε v p n+1 = ε v p n + γu n+1 si g n N ′t r i
A e k c G + H v l h + H v l k -H v nl k χv n + η v ∆t ⇒                  γ v p n+1 = ϵ v p n + γv n+1 si g n V ′t r i al n+1 -χv,tr i al n+1 A e ξv n+1 = ξv n + γv n+1 qv n+1 = qv n -H v l h γv n+1 χv n+1 = χv n + H v l k γv n+1 si g n V ′t r i al n+1 -χv,tr i al n+1 A e -H v nl k γv n+1 χv n
The value of elasto-plastic tangent modulus C ep,v at time step t n+1 is

C ep,v n+1 = k c G H v l h + H v l k -H v nl k χv n+1 + η v ∆t k c G + H v l h + H v l k -H v nl k χv n+1 + η v ∆t

Local phase: discontinuity

Computational procedure for exponential softening is as follows

Start with the best iterative values of nodal displacements, then

• Assume elastic trial step γu,tr i al n+1 = 0 ⇒    α u,t r i al n+1 = α u n ξu,tr i al n+1 = ξu n ; γv,tr i al n+1 = 0 ⇒    α v,t r i al n+1 = α v n ξv,tr i al n+1 = ξv n
• Calculate trial value of internal forces acting at discontinuity and stress-like softening variables

⇒          t u,t r i al n+1 = E A e ε (i ) n+1 -ε v p + Ḡα u n qu,tr i al n+1 = N f A e 1 -exp -ξu n (N f /A e ) G u f ⇒          t v,t r i al n+1 = k c G A e γ (i ) n+1 -γ v p + Ḡα v n qv,tr i al n+1 = V f A e 1 -exp -ξv n (V f /A e ) G v f • Calculate trial value of yield functions φu,tr i al n+1 = t u,t r i al n+1 -N f -qu,tr i al n+1 A e φv,tr i al n+1 = t v,t r i al n+1 -V f -qv,tr i al n+1 A e → if φu,tr i al n+1
≤ 0 the current step in axial direction is indeed elastic, and the values of internal variables at time step t n+1 remain the same as at time

step t n γu n+1 = 0 ⇒            α u n+1 = α u n ξu n+1 = ξu n qu n+1 = qu,tr i al n+1 if φu,tr i al n+1
> 0 the current step is in axial direction is plastic, and the values of internal variables needs to be updated in order to ensure the plastic admissibility of stress. The values of internal variables for exponential softening at time step t n+1 are computed iteratively from the condition φu n+1 ≤ t ol

The values of internal variables at time step t n+1 and iteration ( j ) are

∆ γu,(j) n+1 = φu,(j) n+1 A e -ḠE + K u,( j ) α ⇒      α u,( j +1) n+1 = α u,( j ) n+1 + ∆ γu,(j) n+1 si g n t u,t r i al n+1 ξu,(j+1) n+1 = ξu,(j) n+1 + ∆ γu,(j) n+1 where K u,( j ) α = - (N f /A e ) 2 G u f exp -ξu,(j) (N f /A e ) G u f t u,( j ) n+1 = E A e ϵ (i ) n+1 -ϵ v p + Ḡα u,( j ) n qu,(j) n+1 = N f A e 1 -exp - ξu,(j) n+1 (N f /A e ) G u f → if φv,tr i al n+1
≤ 0 the current step in transverse direction is indeed elastic, and the values of internal variables at time step t n+1 are

γv n+1 = 0 ⇒              α v n+1 = α v n ξv n+1 = ξv n qv n+1 = qv,tr i al n+1 if φv,tr i al n+1
> 0 the current step in transverse direction is plastic. The values of internal variables for exponential softening at time step t n+1 are computed iteratively from the condition

φv n+1 ≤ t ol
The values of internal variables at time step t n+1 and iteration ( j ) are

∆ γv,(j) n+1 = φv,(j) n+1 A e -ḠGk c + K v,( j ) α ⇒      α v,( j +1) n+1 = α v,( j ) n+1 + ∆ γv,(j) n+1 si g n t v,t r i al n+1 ξv,(j+1) n+1 = ξv,(j) n+1 + ∆ γv,(j) n+1 where K v,( j ) α = - (V f /A e ) 2 G v f exp -ξv,(j) (V f /A e ) G v f t v,( j ) n+1 = k c G A e γ (i ) n+1 -γ v p + Ḡα v,( j ) n qv,(j) n+1 = V f A e 1 -exp - ξv,(j) n+1 (V f /A e ) G u f 2.

Global phase

The computation of internal variables in the local phase is followed by the global phase. The global phase of the computation is performed by means of Newton's iterative method. Based on the converged values of internal variables, the values of internal forces are computed and compared with the values of external forces to check whether the equilibrium is satisfied. The equilibrium condition is written as

n el em A e=1 f ext -f i nt ,(i ) -f acc,(i ) e n+1 < tol (2.56)
If the previous condition is satisfied, the computation for new time-step is performed. If not, new iterative values of unknown nodal displacements and pore pressures are to be computed.

The global phase of the computation is solved using Newmark time-integration scheme, which solves the second-order transient problem in the following manner

[41] R n+1 = F n+1 -P u n+1 , un+1 , ün+1 = 0 (2.57)
where

u n+1 = u n + ∆t un + ∆t 2 (0.5 -β) ün + β ün+1 un+1 = un + ∆t (1 -γ) ün + γ ün+1 (2.58)
in which β and γ are parameters controlling stability and numerical dissipation, and ∆t is the time increment.

For a time step t n+1 and iteration i , the global system of equations to be solved is written as

n el em A e=1   1 β∆t 2 M uu + Kuu -K up γ β∆t K T up γ β∆t D pp + K pp   e,(i ) n+1    ∆ ū ∆ p‵    e,(i ) n+1 =    r u r p    e,(i ) n+1 (2.59)
where Ke uu is the element tangent stiffness matrix, and r e,(i ) u,n+1 and r e,(i ) p,n+1 are residuals pertaining to the solid and the pore fluid part.

After solving the global system of equations, the new iterative values of unknown fields are updated as

ū(i+1) n+1 = ū(i) n+1 + ∆ ū(i) n+1 p(i+1) n+1 = p(i) n+1 + ∆ p(i) n+1 (2.60)
In each iteration, the tangent stiffness matrix for every element is constructed.

The tangent stiffness matrix of an element is obtained by linearization of Equations (2.36) and (2.37) with respect to nodal displacements and displacement jumps.

If the element is elasto-viscoplastic then only Equation (2.36) with respect to nodal displacements is linearized since the jump in displacement is not yet activated. The element tangent stiffness matrix is defined as

Ke,(i) uu,n+1 = K e,(i ) uu,n+1 = ∂f i nt (σ u ) δ ū e,(i ) n+1 = L e 0 B s,T u C ep,(i ) n+1 B s u d x (2.61)
Elasto-viscoplastic tangent matrix is written as

C ep,(i ) n+1 =     C ep,u n+1 A e 0 0 0 C ep,v n+1 A e 0 0 0 E I e     (2.62)
where C ep,u n+1 and C ep,v n+1 are elasto-plastic tangent moduli for axial and transverse direction.

If the element is in the softening then both Equations (2.36) and (2.37) are linearized. The element tangent stiffness matrix is obtained by performing static condensation procedure in which the unknown values of displacement jumps are eliminated from the condition that the residual at the discontinuity is equal to zero.

The statically condensed element tangent stiffness matrix is written as

Ke,(i) uu,n+1 = Ke,(i) uu,n+1 = K e,(i ) uu,n+1 -F e,(i ) n+1 H e,(i ) n+1 + K α -1 F e,(i ),T n+1 + K d (2.63)
where

K e,(i ) uu,n+1 = ∂f i nt (σ u ) δ ū e,(i ) n+1 = L e 0 B s,T u C ep,(i ) n+1 B s u d x F e,(i ) n+1 = ∂f i nt (σ u ) δα e,(i ) n+1 = L e 0 B s,T u C ep,(i ) n+1 Ḡd x F T + K d e,(i ) n+1 = ∂h δ ū e,(i ) n+1 = L e 0 ḠT C ep,(i ) n+1 B s u d x + K d H + K α e,(i ) n+1 = ∂h δα e,(i ) n+1 = L e 0 ḠT C ep,(i ) n+1 Ḡd x + K α (2.64)
If the current step in the softening is elastic, then

K d = C * B s u ; K α = 0; C * =     E A e 0 0 0 k c G A e 0 0 0 0     (2.65)
Else, if the current step in the softening is plastic, then

K d = 0; K α =     K u α A e 0 0 0 K v α A e 0 0 0 0     (2.66)
Remark 1: When computing the element tangent stiffness matrix all different combinations that can occur, have to be taken into account. For example, in axial direction softening can occur while in transverse direction element is elastoviscoplastic. After the appropriate tangent stiffness matrix for every element is constructed, an assembly procedure to take into account the contribution of every element to global equilibrium is performed. We note that in the assembly procedure denoted by operator A, the local-global transformation procedure is included.

Remark 2:

We note that the equal order of finite element interpolation (linear interpolation) is used for both the displacement and pore pressure fields. If the undrained limit state is considered, with permeability and compressibility matrix equal to zero, then this kind of approximation (unless stabilization techniques are implemented [START_REF] Pastor | Stabilized low-order finite elements for failure and localization problems in undrained soils and foundations[END_REF]) can cause stability issues and finite elements which satisfy Babuska-Brezzi condition [START_REF] Babuška | The finite element method with lagrangian multipliers[END_REF][START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers[END_REF], or Zienkiewicz-Taylor mixed patch test [START_REF] Zienkiewicz | The patch test for mixed formulations[END_REF] should be used (e.g., see [START_REF] Sanavia | A formulation for an unsaturated porous medium undergoing large inelastic strains[END_REF]). However, if the undrained limit state is not considered (as is in our case), then equal order of interpolation for unknown fields is justified [5,[START_REF] Zienkiewicz | Static and dynamic behaviour of soils: a rational approach to quantitative solutions. i. fully saturated problems[END_REF].

R n+1 = F n+1 -P u n+1 , un+1 = 0 (2.67) where un+1 = 1 ∆t u n+1 -u n (2.68)
For a time step t n+1 and iteration i , the global system of equations to be solved under the quasi-static assumption is written as

n A e=1   Kuu -K up 1 ∆t K T up 1 ∆t D pp + K pp   e,(i ) n+1    ∆ ū ∆ p‵    e,(i ) n+1 =    r u r p    e,(i ) n+1
(2.69)

Remark 4:

In numerical simulations of structures built of dry material in quasistatic setting, the global system of equations is reduced to the system of algebraic equilibrium equations. The global phase of computation is solved using Newton's iterative method. For a time step t n+1 and iteration i , the global system of equations to be solved under the quasi-static assumption is written as

n A e=1 Ke,(i) uu,n+1 ∆ ūe,(i) n+1 = {r u } e,( i ) n+1 (2.70)
For solving the global phase, the arc-length method can also be used in which both the current values of external load and displacements are sought iteratively.

Remark 5:

We note that the boundary conditions in the numerical test for identification of lattice Young's modulus shown in Figure 2.4a, correspond to the conditions of constrained lateral deformations. The conditions of constrained later deformations are typically met in the oedometer test, which is commonly used in soil mechanics to determine the consolidation characteristics of soils. In the oedometer test the movement of the sample is allowed only in the vertical direction, which makes this test very useful for obtaining the value of constrained (oedometer) elastic modulus.

The relation between the oedometer elastic modulus E oed and Young's modulus E for continuum can be derived by exploiting Hooke's law. Hooke's law for isotropic material in three dimensions with ν as Poisson's ratio is written as

            σ 11 σ 22 σ 33 σ 23 σ 13 σ 12             = E (1 + ν)(1 -2ν)             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                         ε 11 ε 22 ε 33 ε 23 ε 13 ε 12             (2.71)
where σ is the stress vector and ε is the strain vector.

Hence, the oedometer elastic modulus E oed and Young's modulus E for continuum are related by

E oed = E (1 -ν) (1 + ν)(1 -2ν) (2.72)
The ratio between the oedometer modulus E oed and Young's modulus E for different values of Poisson's ratio are shown in Table 2.1.

Numerical results

In this section, we present the results of several numerical simulations, which serve to illustrate the coupled discrete beam lattice model performance. All numerical implementations and computations are performed with the research version of the computer code FEAP, developed by R.L. Taylor [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. In all numerical simulations, the finite element mesh is generated by using Delaunay triangulation and GMSH software [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The cross-sectional properties of the Timoshenko beam finite elements are computed from the Voronoi diagram by using MATLAB software [START_REF]version 8.5.0.197613[END_REF].

Illustrative examples

In this section, we demonstrate the capabilities of the proposed discrete beam lattice model on a set of illustrative examples for a structure with very small water content (i.e., structure built of dry, non-porous material). We describe in details the phases of the macro-scale response of the structure computed with the proposed model. Namely, we aim to demonstrate the idea of multi-scale modeling. In other words, we aim to show that with the appropriate modeling of the micro-scale, we are able to capture the true response of the structure on the macro-scale level.

Hardening response of discrete beam lattice model

In our first example, we observe a cantilever beam modeled with a single Timoshenko beam finite element, shown in Figure 2.9a. The material parameters of the Timoshenko beam finite element are: Young's modulus E = 10 5 kPa, the yield limit σ y = 10 kPa, the linear hardening modulus H l h = 10 3 kPa, the linear hardening modulus H l k = 10 3 kPa, the nonlinear hardening parameter H nl k = 10 2 , and the viscosity parameter η = 10 3 kPa s.

We subject the beam to the increasing triangular pulse load program (Figure The computed force-displacement curve for the case when the ratio between yield stress values in compression and tension is equal to 1 is shown in Figure 2.13a.

We can see that the computed hysteresis curve is of a similar shape as the one obtained in the cyclic loading experiments with materials such as steel. In Figure 2.13b, we can see that for the case when the same ratio is equal to 10, the shape of the computed hysteresis curve corresponds to the one observed in the experiments with materials such as concrete.

Hardening and softening response of discrete beam lattice model

The macro-scale response of structure built of cohesive heterogeneous material, such as rock or concrete, can be divided into three phases. The first phase is the linear elastic phase with no cracks occurring in the structure. The second phase represents the fracture process zone where cracks start to form and propagate, resulting in a loss of stiffness, but the load carrying capacity of the structure still increases. In the third phase, cracks coalescence results in an ultimate failure mode.

Namely, after reaching the ultimate load value, one large macro-crack is starting to form and propagate, leading to a complete failure of the structure.

The constitutive model chosen for the Timoshenko beam finite element in combination with the enhanced kinematics in terms of embedded strong discontinuities in the axial and transverse direction can indeed capture all three phases of the macro-scale response. With the aim of illustrating these phases, we perform a numerical simulation of uniaxial tension test on a rock specimen in the quasi-static framework. We perform uniaxial tension test by imposing uniform vertical displacements on the top base of the specimen while keeping the bottom base fixed in the horizontal and vertical direction (Figure 2.14a). The global phase of computation under the quasi-static assumption is solved by using Newton's iterative method in combination with the line search algorithm [START_REF] Matthies | The solution of nonlinear finite element equations[END_REF] (see Remark 4).

E = 10000 ν = 0.2 σ y,t = 18 H l h,t = 1000 µ G,t = 22 G f ,t = 15 σ G,t = 1 σ y,c = 180 H l h,c = 1000 µ G,c = 220 G f ,c = 300 σ G,c = 10 σ y,s = 22 H l h,s = 1000 µ G,s = 26 G f ,s = 150 σ G,s = 1
The geometry of the specimen and boundary conditions are shown in Figure 2.2. To take into account material heterogeneities, the ultimate values of stresses are randomly assigned to every finite element using Gaussian random distribution defined with mean µ G and standard deviation σ G (with σ 2 G as the variance). The probability density function (PDF) of Gaussian random distribution of Gaussian random variable X is written as (Figure 2.15)

f (x |µ G , σ G ) = 1 σ G 2π e -1 2 (x-µ G ) 2 σ 2 G (2.

Validation examples

In this section, we aim to validate the coupled discrete beam lattice model of structure built of a saturated porous medium. For the quasi-static case, we compare our results against those computed with a continuum model using commercial software PLAXIS [START_REF] Swolfs | PLAXIS 2D Manuals[END_REF]. In PLAXIS, the coupling between the solid phase and the pore fluid is governed by Biot's porous media theory, and the equal order of interpolation for both the displacement and pressure fields is available. For such comparison, we choose 15-node triangular elements with fourth-order polynomial interpolations. For the dynamic case, we compare the computed results against the reference values provided in [START_REF] De Boer | One-dimensional transient wave propagation in fluid-saturated incompressible porous media[END_REF].

Quasi-static case

First, we observe a saturated poro-elastic column, and we perform a numerical simulation of a one-dimensional plane strain consolidation test. The stress state in this example can be regarded as homogeneous. Second, we observe a saturated

Discrete model of structure built of saturated poro-plastic medium

We repeat the computation with a discrete model with the ramp loading reaching its maximum value at t = 1 s. The computed results are shown in Figures 2.23a and 2.23b. We compare computed results against those shown in Figures 2.21a and 2.21b. From Figure 2.23b, we can conclude that the maximum value of excess pore pressure is, as expected, greater in the case of higher loading rate. This observation is in accordance with Terzaghi's principle of effective stresses. In the case of undrained conditions, where water is not allowed to drain, the total applied force will be carried by the water in the pores because the water is assumed to be much stiffer than the porous skeleton. If the water is allowed to drain, then the part of the force will be carried by the solid phase and part by the water in the pores. In the case of a higher loading rate, the water has less time to drain; hence, the maximum value of excess pore pressure is greater than in the case of lower loading rate. We observe a good match between the results computed with the discrete and continuum model, with the results being practically mesh independent. We can conclude that nonlinear behavior results in an increase in the values of computed vertical displacements and excess pore pressures, as previously observed in [5,[START_REF] Siriwardane | Two numerical schemes for nonlinear consolidation[END_REF].

Dynamic case

In this example, we perform a validation computation of the proposed discrete model of a saturated porous medium in a dynamic setting, where we compare the computed results against the reference values provided in [START_REF] De Boer | One-dimensional transient wave propagation in fluid-saturated incompressible porous media[END_REF]. We observe a saturated poro-elastic column (Figure 2.30a) subjected to step and sinusoidal loading (Figures 2.31a and 2.31b), defined with following expressions

F = 3 [kN/m 2 ] F = 3[1 -cos(ωt )] [kN/m 2 ]; ω = 75 s -1
The computation in a continuum model is performed for Lame's parameters We can conclude that a good match between the computed results and reference values is obtained, with the results being practically mesh independent.

Plane strain compression test

With the aim to investigate the influence of coupling effects on the failure mechanisms in a saturated porous medium, we simulate a plane strain compression test on a water-saturated specimen. The geometry of the specimen and boundary conditions are shown in Figure 2

E = 20000 ν = 0.2 σ y,t = 10 H l h,t = 2000 µ G,t = 12 G f ,t = 20 σ G,t = 1 σ y,c = 100 H l h,c = 2000 µ G,c = 120 G f ,c = 100 σ G,c = 10 σ y,s = 12 H l h,s = 2000 µ G,s = 15 G f ,s = 100 σ G,s = 1
The material parameters of the Timoshenko beam finite element are shown in Table 2.3. The ultimate values of stresses are randomly assigned to every finite element using Gaussian random distribution. The ultimate shear stress for each Timoshenko beam finite element is defined with Mohr-Coulomb law, written as

σ * f ,s =    σ f ,s + |σ ′ | • tan φ, σ ′ < 0 σ f ,s , σ ′ ≥ 0 (2.74)
where σ f ,s is the ultimate shear stress independent of the normal effective stress (i.e. cohesion, c), φ is the angle of internal friction of material and σ ′ is the effective normal stress in the element, assumed negative in compression.

The angle of internal friction is selected as φ = 11 • . The coefficient of permeability is k = 10 -7 m/s, the specific weight of the water is γ w = 10 kN/m 3 , Biot's constant is b = 1 and Biot's modulus is 1/M → 0.

Our goal is to capture the post-peak response and the failure mechanisms in a water-saturated specimen. For this reason, the test is performed with vertical displacements imposed on the top base of the specimen with a constant rate. We

E = 160 ν = 0.3 σ y,t = 0.03 H l h,t = 60 σ f ,t = 0.09 G f ,t = 0.6 σ y,c = 0.30 H l h,c = 60 σ f ,c = 0.90 G f ,c = 60 σ y,s = 0.07 H l h,s = 60 σ f ,s = 0.13 G f ,s = 30
We simulate the response of the rigid footing by imposing uniform vertical displacements with a constant rate along the length of the footing. We assume that the inertial effects are negligible. To investigate the coupling effects, we perform the The proposed coupled discrete beam lattice model shows an efficient performance in terms of fast convergence rates and efficient computational time. In Table 2.5), we give the residual and energy convergence rates for a typical time step in the post-peak part of the response of the soil subjected to footing load. The computational time depends on the density of the mesh, the loading program, the loading rate, and the selected time step. For the case of coupled soil-foundation system (k = 1 • 10 -6 m/s ), the computational time for loading rate v = 5

Introduction

In this Chapter, we propose a novel approach for numerical modeling of fluidstructure interaction nonlinear problems [START_REF] Hadzalic | Fluid-structure interaction system predicting both internal pore pressure and outside hydrodynamic pressure[END_REF][START_REF] Hadzalic | Theoretical formulation and seamless discrete approximation for localized failure of saturated poroplastic structure interacting with reservoir[END_REF][START_REF] Hadzalic | Novel method for acoustic fluid-structure interaction in application to overall safety of structures in quasi-static setting[END_REF]. The main original aspects are

External fluid model in mixed Lagrangian formulation 3.2.1 Assumptions and governing equations

We assume a small motion of inviscid, isentropic, and homogeneous fluid with constant density. The fluid motion starts from the rest, the fluid velocities remain small, and the fluid flow is considered to be irrotational as the conditions typically met for fluid placed in a reservoir. The governing momentum and continuity equations derived from acoustic wave theory are written as

ρ v + ∇p = 0 (3.1) ∇ • v + ṗ β = 0 (3.2)
where v is the velocity vector, p is the pressure, ρ is the mass density, and β is the bulk modulus of external fluid.

For irrotational flows, the vorticity is constrained to zero

∇ × v = 0 (3.3)

Mixed displacement/pressure based formulation

Mixed displacement/pressure based finite element formulation for acoustic fluids is based on the analogy with the mixed formulations for nearly incompressible solids [START_REF] Zienkiewicz | The Finite Element Method[END_REF][START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF][START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. Here, we assume that the fluid motion starts from the rest with the fluid displacements remaining small. Moreover, we follow [START_REF] Hamdi | A displacement method for the analysis of vibrations of coupled fluid-structure systems[END_REF][START_REF] Wilson | Finite elements for the dynamic analysis of fluidsolid systems[END_REF] by assuming that the fluid motion remains irrotational. However, contrary to the penalty-type approach for enforcing irrotational constraint, here we use the mixed variational formulation to achieve the same goal. The main advantage of such an approach is in resulting set of degrees of freedom, featuring both fluid displacements and pressure, and providing a seamless connection to the poro-plasticity formulation described in Chapter 2. The resulting variational formulation is equivalent to [START_REF] Bathe | A mixed displacement-based finite element formulation for acoustic fluid-structure interaction[END_REF][START_REF] Wang | Displacement/pressure based mixed finite element formulations for acoustic fluid-structure interaction problems[END_REF][START_REF] Wang | On mixed elements for acoustic fluid-structure interactions[END_REF], and can be written as

Π e f = Ω e f - p 2 2β -p ∇ • u - Λ • Λ 2ϑ + Λ • ∇ × u -u • f b d Ω (3.4)
where p is the pressure, u is the displacement vector, and Λ is the 'vorticity moment'

or Lagrangian multiplier enforcing zero vorticity in terms of perturbed Lagrangian with ϑ as the penalty parameter. The parameter β is the bulk modulus of the external fluid, and f b is the external load vector, that next to the body forces also includes inertia force -ρ ü.

The corresponding weak form is obtained from the first variation of Equation (3.4) and is written as

Ω e f δp - 1 β p -∇ • u + -p ∇ • δu + Λ ∇ × δu -δu • f b + δΛ - 1 ϑ Λ + ∇ × u d Ω = 0 (3.5)
We further ought to introduce the finite element approximations for the displacements, pressure, and 'vorticity moment'. For a typical finite element, we have

u = N f u ū; p = N f p p; Λ = N f λ λ u = N f u ü; p = N f p p; Λ = N f λ λ ∇ • u = (∇ • N f u ) ū = V f ū ∇ × u = (∇ × N f u ) ū = D f ū (3.6)
where N f u , N f p , and N f λ are interpolation matrices, ū is the vector of unknown nodal displacements, p is the vector of unknown nodal pressures, and λ is the vector of unknown nodal 'vorticity moments'.

The end result of the finite element discretization procedure is the following system of equations governing the discrete problem

n el em A e=1            A uu 0 0 0 0 0 0 0 0       e            ü p λ           e +       0 L up L uλ L T up L pp 0 L T uλ 0 L λλ       e            ū p λ           e =            f f 0 0            e      (3.7)
where

A e uu = Ω e f ρN f ,T u N f u d Ω; L e up = - Ω e f V f ,T N f p d Ω; L e uλ = Ω e f D f ,T N f λ d Ω L e pp = - Ω e f 1 β N f ,T p N f p d Ω; L e λλ = - Ω e f 1 ϑ N f ,T λ N f λ d Ω (3.8)
The 'vorticity moment' degree of freedom can be statically condensed on the element level so that the only unknown variables remaining on the global level are displacements and pressures.

The strong form of the governing equations follows from Equation (3.5)

∇p + ∇ × Λ -f b = 0 (3.9) ∇ • u + p β = 0 (3.10) ∇ × u - Λ ϑ = 0 (3.11)
Based on previously made assumptions, we note that Equations (3.9-3.10) are penalized version of the strong form described in (3.1-3.3).

Finite element approximation post-processing

Mixed displacement/pressure based finite element formulation, in combination with the proper choice of finite elements that satisfy the inf-sup condition, has

shown to be a very efficient tool when dealing with nearly incompressible behavior of acoustic fluids in the frequency or dynamic analysis. The range of elements satisfying the inf-sup condition can be found in literature [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF][START_REF] Bathe | Finite Element Procedures[END_REF][START_REF] Chapelle | The inf-sup test[END_REF][START_REF] Bathe | The inf-sup condition and its evaluation for mixed finite element methods[END_REF].

In our numerical simulations of fluid-structure interaction in the dynamic setting for the external fluid domain we choose linear finite element approximations for the displacements, with constant approximations for the pressure and the 'vor- For the purpose of solving the global system of equations, we can statically condense (e.g., [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF]) the pressure and the 'vorticity moment' unknowns at the element level. We can then obtain the values of the pressures and the 'vorticity moments' from the computed values of displacements. However, our goal is to directly connect the external fluid finite elements with structure finite elements at the common nodes in order to ensure the direct exchange of both the motion and the pressure at the fluid-structure interface. This is achieved by reconstructing the pressure field for the external fluid finite elements, by extrapolating the pressure computed inside an element to the nodes of a Q4 finite element used for the displacement approximation. In particular, the pressure at each node is evaluated as an average value of the pressures computed in the finite elements that share that node, as shown in Figure 3.2. For the mesh of regular elements, the error in the value of the pressure computed in this manner is only limited to the finite elements located at the bottom (and/or top, depending on the boundary conditions) of the mesh. The error in the computed pressure increases if the size of the finite elements varies along the y-coordinate. However, by increasing the density of the mesh, the error can be decreased to better control the accuracy. The proposed approach enables us to compute the displacements, pressure, and 'vorticity moment' on a Q4-P1-Λ1 finite element by performing the static condensation process. With such a postprocessing procedure, we obtain the 'Q4-P4' finite element, which we can then use for direct transfer at the fluid-structure interface. we directly obtain 'Q4-P4' finite element, without any post-processing procedures needed. The finite element approximation of this kind is predictive enough for delivering the solution in terms of proper treatment of the fluid-structure interface, as it will be shown in the numerical simulations.

Remark 1:

We note that practically same results for the quasi-static damreservoir interaction presented in this thesis (section 3.3.1.2) can be obtained by using Q4-P1-Λ1 instead of Q4-P4-Λ4 finite elements. However, in this case, a small mass has to be added in the external fluid model.

Remark 2:

In general case, the penalty parameter ϑ in mixed up -Λ formulation is a numerically large value ranging from 10 2 β ≤ ϑ ≤ 10 6 β. By assigning a numerically small value to a ϑ, the up formulation is obtained [START_REF] Wang | Displacement/pressure based mixed finite element formulations for acoustic fluid-structure interaction problems[END_REF].

Numerical results

In this section, we present the results of several numerical simulations, which serve to illustrate the capabilities of the proposed numerical approach for dealing with acoustic fluid-structure interaction problems. All numerical implementations and computations are performed with the research version of the computer code FEAP, developed by R.L. Taylor [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. In all numerical simulations, the finite element mesh is generated by using Delaunay triangulation and GMSH software [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The cross-sectional properties of the Timoshenko beam finite elements are computed from the Voronoi diagram by using MATLAB software [START_REF]version 8.5.0.197613[END_REF].

Quasi-static case

Rigid and flexible water container

The first example presents a simple validation computation of the proposed numerical model of acoustic fluid-structure interaction in the quasi-static setting.

Namely, we study a concrete water container with the geometry shown in Figure 3.4a.

Dam-reservoir system

In this numerical example, we observe a small-sized dam subjected to the selfweight and the hydrostatic loading. shown in Figure 3. We first apply an increase of the self-weight to simulate the construction of the dam, followed by the hydrostatic loading of the external reservoir. In the final part of this analysis, we would like to illustrate the capability of the proposed approach to represent the localized failure, where the cracks coalescence results with the ultimate failure mode. In this manner, we seek to investigate the safety factor of the complete dam structure, which (contrary to usual factors for different materials) ought to be defined for particular load pattern. This is done by computing the maximum overload that can be applied on top of (already) acting self-weight and hydrostatic pressure of the external reservoir. The computation is carried out by imposed displacement, bringing the dam structure each time to the ultimate state of localized failure. We present both results computed for either horizontal or vertical overload, quantifying the localized failure brought by increasing horizontal or vertical load, respectively.

We start by assuming that the behavior of the dam remains linear elastic, and we apply first the self-weight and then the hydrostatic loading. Here, we analyze in detail the results of the fluid-structure interaction obtained with the mixed displacement/pressure based external fluid finite element formulation. Next, we assume that the behavior of the dam is inelastic. We introduce yield and fracture limits for Timoshenko beam finite elements to take into account the possibility of plastic zones and cracks formation as a result of combined body and hydrostatic loading. Finally, we determine the admissible horizontal and vertical overloads of the dam by imposing corresponding displacements on the top base of the dam until ultimate failure. In the first two computation phases, application of the self-weight and hydrostatic loading, the time step is set to ∆t = 0.1 day. In the final stage of localized failure, the time step is set to ∆t = 1 s.

→ Linear elastic behavior

Young's modulus and Poisson's ratio of the Timoshenko beam finite element are E = 10 4 MPa and ν = 0. The coefficient of permeability is k = 10 -7 m/s, Biot's constant is b = 1 and Biot's modulus is M = 10 6 MPa. The bulk modulus of the outside water is β = 10 3 MPa. The value of penalty parameter is ϑ = 10 -5 β. We compare the results obtained with the mixed displacement/pressure based finite element formulation of fluid motion against analytical values. Figure 3.9a show the contours for the computed pressure field in the reservoir and in a porous medium from which the dam is built, by using a numerical model of acoustic fluid-structure interaction. The results are practically identical with the analytical solution for the hydrostatic pressure distribution (Figure 3.9b). From the computed excess pore pressure distribution inside the dam, we can conclude that the water in the dam reservoir can keep the porous material from which the dam is built fully saturated in every time step, which is confirmed by the results in Figure 3.10a showing the time evolution of total excess pore pressure at the bottom of the dam with respect to applied loading. The horizontal displacement of the tip of the dam → point A (3,[START_REF] Réthoré | A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks[END_REF][m] gives additional insight into the results of fluid-structure interaction. From Figure 3.10d, we can see that the total horizontal displacement of the tip of the dam, after hydrostatic loading has been applied is a sum of displacements due to self-weight, hydrostatic forces and hydrostatic pressures. We can verify this if we perform the uncoupled computation of the structure response (Biot's constant b = 0 → matrix K e up is a null matrix) for the hydrostatic loading case. This way, we can obtain the displacements of the dam only from hydrostatic forces exerted on the structure. To validate our results, we perform an uncoupled computation on a model of a dam alone, where we apply hydrostatic loading in terms of equivalent nodal horizontal forces. We can conclude that the value of horizontal displacement for the uncoupled fluid-structure interaction case and the case of the dam alone are practically the same (Figure 3.10e). Next, we subtract the value of displacements for the cou-

→ Inelastic behavior

In this case, we admit that due to the self-weight and the hydrostatic loading, plastic zones and cracks can form in the dam. Thus, we introduce yield and fracture limits for the Timoshenko beam finite element in tension, compression and shear: → Localized failure analysis of dam structure for either horizontal or vertical overload The self-weight and the hydrostatic loading are considered as the fundamental load cases acting upon the dam, which must be used for the structure design. During the dam life cycle, one can imagine occasional extreme loads that can threaten structure integrity and lead to localized failure. Here, we want to quantify the remaining resistance to possible overload (on top of dam dead-load and reservoir pressure) that can be applied on the dam leading to the complete failure. The overload of this kind can be interpreted as the structure safety factor, which has to be properly taken into account in the design process. In particular, we carry on such a safety factor computation by computing the final failure mode for either applied horizontal or vertical overloads.

Dynamic case

Analytical solutions for hydrodynamic pressures

When the fluid-structure systems, such as dam-reservoir are subjected to extreme dynamic loads brought by for example seismic ground movement, next to the hydrostatic pressures and forces, additional hydrodynamic pressures and forces are exerted on the upstream face of the structure. Their values and distribution have to be quantified either with analytical solutions proposed in the literature or with numerical methods, in order to provide the sound design of the structure.

The problem of evaluating the value and the distribution of the hydrodynamic pressure was first examined in the work of Westergaard who provided a solution for the case of a rigid dam with vertical upstream face subjected to horizontal harmonic ground motion [START_REF] Westergaard | Water pressures on dams during earthquakes[END_REF], followed by von Kármán who provided the solution very close to Westergard's using linear momentum balance principle [START_REF] Kármán | Discussion of water pressures on dams during earthquakes[END_REF]. Later, Chwang and Housner, also using the momentum method, derived the analytical solution for the general case of an inclined upstream face of the dam subjected to constant horizontal acceleration [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 1. momentum method[END_REF]. In the second part of his work, Chwang derived the analytical solution using two-dimensional potential flow theory [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 2. exact theory[END_REF].

In addition to the inclined upstream face of the dam, Liu proposed an analytical solution for the sloped reservoir base [START_REF] Liu | Hydrodynamic pressures on rigid dams during earthquakes[END_REF].

→ Westergaard's solution

Westergaard [START_REF] Westergaard | Water pressures on dams during earthquakes[END_REF] is the first who studied the problem of evaluating the value and distribution of hydrodynamic pressure exerted on the upstream face of the dam. His work was focused on the simple two-dimensional dam-reservoir system subjected to the horizontal harmonic ground motion (Figure 3.19). The dam was assumed to be rigid with vertical upstream face, and the length of the reservoir was assumed to be infinite. Westergaard derived an analytical solution for the hydrodynamic pressure distribution in terms of series of sine functions. Maximum hydrodynamic pressure distribution on the vertical upstream face of the rigid dam, according to Westergaard, is described with the following expression

p = 8a 0 ρH π 2 n 1,3,5,... 1 n 2 c n sin nπy 2H c n = 1 - 16ρH 2 n 2 g K T 2 (3.12)
where a 0 is the maximum horizontal acceleration of the foundation, ρ is the density of the retained water, H is the depth of the reservoir, K is the bulk modulus of water, and T is the period of the horizontal acceleration of the foundation. This expression, even though represents a conservative approximation, gives satisfactory results and is widely used in everyday engineering practice.

According to Westergaard, the total horizontal force exerted on the vertical, upstream face of the dam is equal to

F = F x = 0.543ρa 0 H 2 (3.14)

→ Von Kármán's solution

An analytical solution very close to Westergaard's was derived by von Kármán.

In his work, von Kármán [START_REF] Kármán | Discussion of water pressures on dams during earthquakes[END_REF] exploited linear momentum-balance principle and derived an expression for hydrodynamic pressure distribution which reads as

p = 0.707ρa 0 H y (3.15)
According to von Kármán, the total horizontal force exerted on the vertical, upstream face of the dam is equal to .20 Dam-reservoir system studied by [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 1. momentum method[END_REF] and [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 2. exact theory[END_REF] In their analytical solution, as in the Westergaard's and von Kármán 's solutions, the value of hydrodynamic pressures ranges from zero at the top of the reservoir to maximum at the bottom of the reservoir. Chwang, in the second part of his work [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 2. exact theory[END_REF], derived an exact solution for hydrodynamic pressure distribution using two-dimensional potential flow theory. Here, except for the case of the vertical upstream face, the maximum value of hydrodynamic pressure does not occur at the bottom of the reservoir but is moved up to a certain distance.

F = F x = 0.
The hydrodynamic pressure distribution on the inclined upstream face is described with the following expression

p = C p ρa 0 H (3.17)
where C p is the pressure coefficient. The values of the pressure coefficient for different angles of inclination are shown in Figure 3.21, for both methods. The expression for computing the values of pressure coefficient can be found in [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 1. momentum method[END_REF] and [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 2. exact theory[END_REF].

The total horizontal and vertical force exerted on the inclined upstream face of the dam are equal to

F x = C x ρa 0 H 2 F y = C y ρa 0 H 2 (3.18)
where C x and C y are force coefficients. The values of force coefficients for different angles of inclination are shown in Figure 3.22, for both methods. The expressions for computing the values of force coefficient can be found in [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 1. momentum method[END_REF] and [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 2. exact theory[END_REF]. where n, m are integers, and c is the acoustic wave speed.

The computed values of first four frequencies are shown in Table 3.1. We can conclude that by increasing the mesh density, the results obtained with Q4-P1-Λ1

approach the analytical solution.

Linear elastic separator wall

In this example, we observe a linear elastic separator wall 12 m high, and 1. We subject the separator wall and the contained fluid to a horizontal ground acceleration reaching its maximum value of a 0 = 1 m/s 2 at t = 1s, after which is kept constant (Figure 3.26). According to [START_REF] Chopra | Dynamics of Structures: Theory and Applications to Earthquake Engineering[END_REF], we model the horizontal ground movement in terms of equivalent horizontal forces. Here, we consider a system with a fixed base, and we apply at every node of the finite element mesh an equivalent force with a value equal to the product of the corresponding mass and the imposed base acceleration. We compare the computed results against those provided by Chwang [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 2. exact theory[END_REF], with the aim to validate the proposed model for predicting hydrodynamic pressures and hydrodynamic forces exerted on the upstream face of the structure.

Dam-reservoir system

In this example, we observe a small-size gravity dam, shown in Figure 3.29. The length of the reservoir is chosen as L = 20.40 m in order to eliminate the influence of the boundary effects. First, we assume that the response of dam remains linear elastic, and we subject the dam-reservoir system to a horizontal ground acceleration reaching its maximum value of a 0 = 1 m/s 2 at t = 1s, which is then kept constant (Figure 3.30a). The time step is set to ∆t = 0.01 s. The horizontal ground movement is modeled in terms of equivalent horizontal forces acting on the system with fixed base [START_REF] Chopra | Dynamics of Structures: Theory and Applications to Earthquake Engineering[END_REF]. We compare computed results against analytical solutions provided by Chwang [START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 2. exact theory[END_REF].

Concluding remarks

In this Chapter, we presented a novel approach to numerical modeling of the acoustic fluid-structure interaction nonlinear problems. The structure is represented in terms of a saturated porous medium. The structure response is modeled with coupled discrete beam lattice model, which is based on Voronoi cell representation with inelastic Timoshenko beam finite elements enhanced with additional kinematics as cohesive links. The internal coupling is handled with Biot's porous media theory. The motion of the external fluid in interaction with the structure is described with the Lagrangian formulation and mixed displacement/pressure based finite element approximation. The main advantage of the proposed numerical model for acoustic fluid-structure interaction is in the resulting set of degrees of freedom per node of structure and external fluid finite elements, which permits for the elements to be connected directly at the common nodes without any need for special numerical consideration of the fluid-structure interface. As a result, the direct exchange of both motion and the pressure at the common boundary is ensured, and all unknown fields can be computed in a fully monolithic manner. We have confirmed through numerical simulations that the proposed numerical model of interaction is capable of predicting both the hydrostatic and hydrodynamic pressures and forces exerted on the upstream face of the structure in close agreement with analytical solutions. In addition, we are able to model the formation of cracks and plastic zones in the structure subjected to fundamental load cases, such as the self-weight and hydrostatic loading.

Finally, we have illustrated the proposed model ability to represent localized failure modes, and thus quantify the overall safety of the dam structure to potential overload with respect to the fundamental load cases. The computation of this kind in a quasi-static setting is performed by imposing either horizontal or vertical displacement upon the dam deformed configuration, which is obtained under the dead-load and reservoir pressure. In a dynamic setting, the localized failure of dam structure is brought by subjecting the dam-reservoir system to the linear increasing horizontal ground acceleration, which is modeled in terms of equivalent horizontal forces. We have shown qualitatively how different types of failure modes affect the computed value of admissible overload compared to the uncoupled case of structure response. Namely, for different types of loading programs, the presence of the pore fluid can increase or decrease the computed value of admissible overload.

Thus, the ability of the proposed formulation to quantify the overall structure safety for a particular loading program is very important from the standpoint of engineering practice dealing with issues of the structure design, both in terms of the structural integrity and structural durability.

For a more detailed explanation of computed admissible overload for a particular loading program, we go back to Terzaghi's principle of effective stresses (see Chapter 2), which states that the total normal stress is equal to the sum of effective stress carried by the solid phase, and pore pressure carried by the pore fluid. We can rewrite this principle in terms of total reaction, which we can interpret as the ultimate load value structure can withstand until failure. Namely, the total reaction in the coupled model R c at each time step is equal to the sum of contributions of effective reaction carried by the solid phase R e , and the reaction carried by the fluid in pores R p , written as

R c = R e + R p
For uncoupled model, the total reaction R uc is equal to

R uc = R e
In our case, we can state that the reaction R e obtains approximately the same value regardless of the structure response being coupled or uncoupled. To explain this, we return to Mohr-Coulomb's law, which states that the strength of the material depends on the cohesion and the friction, and that friction is directly proportional to the effective normal stress. The strength of cohesive materials, such as rocks or concrete, dominantly depends on the value of cohesion and contribution of the friction is not as significant. This is the case in our proposed coupled discrete beam lattice model, which makes this model suitable for cohesive types of materials.

The fracture limits in tension, compression, and shear define the ultimate load level structure can withstand until ultimate failure for a particular loading program, and this part is constant regardless of the structure response being coupled or uncoupled. An increase in the ultimate load value can be obtained if the Mohr-Coulomb's law is implemented for defining the fracture limit in shear (see Section 2.5.3). Thus, if the latter is not implemented, we can state that the contribution of the solid phase to the total value of the total reaction is constant. Hence, the value of the total reaction changes due to the contribution of the pore fluid. Depending on the time evolution of the pore pressures in the body of the dam for different types of failure modes, the contribution of the pore fluid can increase or decrease the total reaction compared to the uncoupled case.

The proposed model of interaction, even when pushed to the extreme while computing the final phase of localized failure, exhibits good performance in terms of fast convergence rates and efficient computational time. In Tables 3.2 and 3.3, we give the residual and energy convergence rates for a typical time step per each computation phase of localized failure analysis of pore-saturated dam structure under horizontal overload in a quasi-static setting. The computational times for horizontal and vertical overloads in a quasi-static setting with included first two computation phases are 720.83 s, and 815.52 s, respectively. The computational time for horizontal overload in a dynamic setting is 1050.07 s (Processor: Intel Core i7-3632QM/2.2GHz, RAM 8.00GB, FEAP: single-core).

Introduction

In this Chapter, we extend the two-dimensional (2D) numerical models of poresaturated structure, external fluid and ultimately acoustic fluid-structure interaction previously presented in Chapters 2 and 3 to three-dimensional (3D) setting.

The outline of the Chapter is as follows: In Section 4.2, we describe the main features of the coupled discrete beam lattice model of the structure in a 3D setting.

In Section 4.3, we present the finite element formulation of the proposed coupled discrete beam lattice model, together with the computational procedure. The finite element formulation is given in fully dynamics framework, which can be easily adjusted for quasi-static simplification. In Section 4.4, we provide finite element approximation for external fluid. In Section 4.5, we present the results of several numerical simulations. First, we simulate the response of the structure alone, with the aim to validate the proposed model of the structure and to demonstrate its capabilities to deal with nonlinear behavior and progressive localized failure in structures built of saturated poro-plastic cohesive material. Second, we simulate the response of dam structure under reservoir interaction. In Section 4.6, we give concluding remarks.

Discrete model of structure

For the numerical representation of the structure built of a saturated poro-plastic medium, we extend the proposed 2D plane strain coupled discrete beam lattice model to the 3D setting. For the construction of the discrete lattice model in the 3D framework, we use the same procedure described previously in Chapter 2, Section Voronoi cells and is perpendicular to the polygon shared between these two cells.

We place along each edge of tetrahedra a cohesive link, whose behavior we model with the 1D finite element. The cross-section of each finite element is the polygon shared between two adjacent cells. In general case, the polygonal cross-section for every finite element is of a different shape. To simplify the matter, we replace the polygonal cross-section for every finite element with an equivalent circular cross-section (Figure 4.1). Namely, we compute the diameter of an equivalent circular cross-section from the condition that the area of the polygon obtained from Voronoi cell representation is equal to the area of a circular cross-section. In a 3D numerical model of structure built of a saturated poro-plastic medium, we spread the pore fluid flow across the mesh of tetrahedra (Tet4 -linear tetrahedral finite elements) that coincides with the mesh of tetrahedra obtained by Delaunay tetrahedralization of the domain. For numerical integration on tetrahedral finite elements, we choose nodal point rule which positions the integration points at every node of tetrahedra [START_REF] Gellert | Moderate degree cubature formulas for 3-d tetrahedral finite-element approximations[END_REF]. The integration rule of this kind eliminates the contribution of two nodes and leaves the contribution of the two nodes that correspond to the Timoshenko beam finite element. This allows us to treat the pore pressure as an additional degree of freedom placed at each node of the Timoshenko beam finite element, which we later use for the fluid-structure interface (Figure 4.2). where

1 v 2 ¡ z,2 u 2 2 L e x L y L z L w 2 ¡ x,2 ¡ y,2 v 1 ¡ z,1 u 1 w 1 ¡ x,1 ¡ y,1 α v α u x α w
u T = u, v, w, θ x , θ y , θ z ; ūT = u 1 , v 1 , w 1 , θ x,1 , θ y,1 , θ z,1 , u 2 , v 2 , w 2 , θ y,2 , θ z,2 ; N s u =                 N 1 0 0 0 0 0 N 2 0 0 0 0 0 0 N 1 0 0 0 0 0 N 2 0 0 0 0 0 0 N 1 0 0 0 0 0 N 2 0 0 0 0 0 0 N 1 0 0 0 0 0 N 2 0 0 0 0 0 0 N 1 0 0 0 0 0 N 2 0 0 0 0 0 0 N 1 0 0 0 0 0 N 2                 ; N 1 , N 2 = 1 -x L e , x L e ; M =                 M 0 0 0 0 0 0 M 0 0 0 0 0 0 M 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 ; α T = α u , α v , α w , 0, 0, 0 (4.2)
Here, α u , α v and α w represent displacement jumps in axial, in-plane transverse and out-of-plane transverse direction. The enhanced strain fields are interpolated as

ϵ = B s u ū + Gα (4.3)
where

ϵ T = ε xx , γ x y , γ xz , κ x , κ y , κ z ⇒                    ε xx = d u d x κ x = d θ x d x γ x y = d v d x -θ z κ y = d θ y d x γ xz = d w d x + θ y κ z = d θ z d x ; B s u =                 B 1 0 0 0 0 0 B 2 0 0 0 0 0 0 B 1 0 0 0 -N 1 0 B 2 0 0 0 -N 2 0 0 B 1 0 N 1 0 0 0 B 2 0 N 2 0 0 0 0 B 1 0 0 0 0 0 B 2 0 0 0 0 0 0 B 1 0 0 0 0 0 B 2 0 0 0 0 0 0 B 1 0 0 0 0 0 B 2                 ; B 1 , B 2 = d N 1 d x , d N 2 d x = -1 L e , 1 L e ; G =                 G 0 0 0 0 0 0 G 0 0 0 0 0 0 G 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 (4.4)
where ε xx is the axial strain, γ x y , γ xz are the shear strains, and κ x , κ y , κ z are the curvatures. In what is to follow, we will denote these strains with following symbols ε, γ y , γ z , κ x , κ y and κ z .

The pore pressure field is approximated with tetrahedral finite elements with four nodes (Figure 4.6). The finite element approximation for the pore pressure field is written as

p = N s p p (4.5)
where pT = p 1 , p 2 , p 3 , p 4 ;

N s p = N p 1 , N p 2 , N p 3 , N p 4 ; N p i = a i + b i x + c i y + d i z 6V , i = 1, 4; 6V = d et 1 x 1 y 1 z 1 1 x 2 y 2 z 2 1 x 3 y 3 z 3 1 x 4 y 4 z 4 ; a 1 = d et x 2 y 2 z 2 x 3 y 3 z 3 x 4 y 4 z 4 ; b 1 = -d et 1 y 2 z 2 1 y 3 z 3 1 y 4 z 4 ; c 1 = -d et x 2 1 z 2 x 3 1 z 3 x 4 1 z 4 ; d 1 = -d et x 2 y 2 1
x 3 y 3 1

x 4 y 4 1 ; et c. (4.6)
where V is the volume of the tetrahedral element, x,y,z are global coordinates, x i , y i , z i are nodal coordinates of tetrahedral element, and constants a 2-4 , b 2-4 , c 2-4 , d 2-4 are defined by cyclic interchange of the subscripts in the order 1,

Continuity equation and equations of motion

Continuity equation

The coupling of the mechanics and the pore fluid flow occurs through the axial direction of the Timoshenko beam finite element. The continuity equation for fluid flow through the porous body is written as

1 M ṗ + +b ε - k γ f ∇ • (∇p) = 0 (4.8)
The weak from of the continuity equation is obtained through the principle of virtual works, written as

Ω e Tet 4 δp 1 M ṗ - k γ f ∇ • (∇p) d Ω + L e 0 δp ‵ b εd x = 0 (4.9)
where δp is the virtual pore pressure field interpolated in the same manner as the real pore pressure field

δp = N s p δ p; δ pT = δp 1 , δp 2 , δp 3 , δp 4 ; p ‵ = N s up p‵ ; N s up = N 1 , N 2 ; p‵T = p 1 , p 2 ; δp ‵ = N s up δ p‵ ; δ p‵T = δp 1 , δp 2 (4.10)

Equations of motion

The strong form of the equations of motion for a 3D Timoshenko beam with a circular cross-section is written as

-ρ A e ∂ 2 u ∂t 2 + ∂N ∂x + n(x, t ) = 0 (4.11) -ρ A e ∂ 2 v ∂t 2 + ∂V y ∂x + q y (x, t ) = 0 (4.12) -ρ A e ∂ 2 w ∂t 2 + ∂V z ∂x + q z (x, t ) = 0 (4.13) -ρI e pol ∂ 2 θ x ∂t 2 + ∂M x ∂x + m x (x, t ) = 0 (4.14) -ρI e y y ∂ 2 θ y ∂t 2 + ∂M y ∂x -V z + m y (x, t ) = 0 (4.15) -ρI e zz ∂ 2 θ z ∂t 2 + ∂M z ∂x + V y + m z (x, t ) = 0 (4.16)
where ρ is the mass density, A e is the are of a cross-section, I e y y , I e zz are the second moments of inertia of a cross-section, I e pol is the polar moment of inertia of a cross-section, n, q y , q z , m x , m y and m z are the distributed external loads, and

N = N ′ -bp ‵ A e is the total axial force, V y = V ′ y , V z = V ′ z are the shear forces, M x = M ′ x is the torsion moment, M y = M ′ y , M z = M ′
z are the bending moments (Figure 4.7). The superscript ′ denotes effective force. Here, we assume that the coupling of the mechanics and the internal fluid flow occurs through the axial direction of the Timoshenko beam finite element. The weak form of equations of motion for Timoshenko beam finite element is derived from the d'Alembert principle. The virtual displacement and virtual strain fields are interpolated in the same manner as the real displacement and real strain fields, written as

M z N x L y L z L M x M y V y V z
δu = N s u δ ū + Mδα (4.17) δϵ = B s u δ ū + Gδα (4.18)
where δu T = δu, δv, δw, δθ x , θ y , δθ z ;

δ ūT = δu 1 , δv 1 , δw 1 , δθ x,1 , δθ y,1 , δθ z,1 , δu 2 , δv 2 , δw 2 , δθ x,2
, δθ y,2 , δθ z,2 ;

δα T = δα u , δα v , δα w , 0, 0, 0 ;

δϵ T = δε, δγ y , δγ z , δκ x , δκ y , δκ z (4.19)
The virtual work of external forces is computed as

G ext ,e = δ ūT f ext ,e -δ ūT f acc,e (4.20)
where The virtual work of internal forces for Timoshenko finite beam finite element is equal to

f acc,e = L e 0 N s,T u σd x; σT = ρ A e ü, ρ A e v, ρ A e ẅ,
G i nt ,e = L e 0 B s u δ ū T σd x + L e 0 Gδα T σ u d x (4.22)
where

σ T = N ,V y ,V z , M x , M y , M z = N ′ -bp ‵ A e ,V ′ y ,V ′ z , M ′ x , M ′ y , M ′ z ; σ T = σ T u -σ T p ; σ T u = N ′ ,V ′ y ,V ′ z , M ′ x , M ′ y , M ′ z ; σ T p = bp ‵ A e , 0, 0, 0, 0, 0 (4.23)
The virtual work of internal forces consists of two parts. The first is the virtual work of internal forces in the bulk part of the element (part of the element outside the discontinuity), and the second is the virtual work of internal forces acting at the discontinuity, written as

G i nt ,e = δ ūT f i nt ,e + δα T h e ( 4.24)
where f i nt ,e is the internal force vector, and h e is the residual vector due to discontinuity, computed as

f i nt ,e = L e 0 B s,T u σd x = f i nt ,e (σ u ) -f i nt ,e (σ p ); h e = L e 0 G T σ u d x = L e 0 ḠT σ u d x + t; Ḡ =                 Ḡ 0 0 0 0 0 0 Ḡ 0 0 0 0 0 0 Ḡ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 ; t T = t u , t v , t w , 0 , 0 , 0 (4.25)
Here, t is the internal force vector acting at the discontinuity.

By exploiting the principle of virtual work we obtain 0 = G i nt ,e -G ext ,e = δ ūT (f i nt ,ef ext ,e + f acc,e ) + δα T h e (4.26)

Finite element formulation

Previous equation is satisfied if f acc,e + f i nt ,ef ext ,e = 0 (4.27)

h e = 0 (4.28)
The Equation (4.27) relates to the bulk part of the element, and the Equation (4.28)

relates to the discontinuity.

The internal force vector t acting at the discontinuity is computed from the following conditions

h e = 0 L e 0 Ḡd x = L e 0 - 1
L e d x = -1        ⇒ t =                 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 L e 0 σ u d x (4.29)

Constitutive model

The pre-peak response of the 3D Timoshenko beam finite element in axial, and both transverse directions is described with the elasto-viscoplastic constitutive model with implemented linear hardening [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF], and Fredrick-Armstrong nonlinear kinematic hardening law [START_REF] Armstrong | A mathematical representation of the multiaxial bauschinger effect[END_REF]. The post-peak response of the element is described with exponential softening. The behavior of the element in bending and torsion is purely linear elastic.

Plasticity model

The proposed model of plasticity is described with [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF] 1. Additive decomposition of the total strain into elastic and viscoplastic part, with only the elastic part influencing the value of the stresses

ε = ε e + ε v p κ x = κ e x γ y = γ e y + γ v p y κ y = κ e y γ z = γ e z + γ v p z κ z = κ e z ( 4.30)
For Fredrick-Armstrong nonlinear kinematic hardening law [START_REF] Armstrong | A mathematical representation of the multiaxial bauschinger effect[END_REF], the backstress variables χu , χv , χw are computed as where γu , γv , γw are the plastic multipliers whose values are equal to

χu = H u l k εvp -H u nl k ξu χu χv = H v l k γvp y -H v nl k ξv χv χw = H w l k γvp z -H w nl k ξw χw
γu = 〈 φu (N ′ , χu , qu )〉 η u γv = 〈 φv (V ′ y , χv , qv )〉 η v γw = 〈 φw (V ′ z , χw , qw )〉 η w (4.37)
where η u , η v , η w are the viscosity parameters.

Exponential softening

The exponential softening is described with [START_REF] Ibrahimbegovic | Nonlinear solid mechanics: theoretical formulations and finite element solution methods[END_REF] 1. The yield function defined in terms of stresses (or stress resultants), and dual variables qu , qv , qw

φu t u , qu = |t u | -N f -qu A e ≤ 0 φv t v , qv = |t v | -V y, f -qv A e ≤ 0 φw t v , qw = |t w | -V z, f -qw A e ≤ 0 (4.38)
where N f , V y, f , V z, f are the fracture limits, and qu , qv , qw are the stress-like softening variables.

The stress-like softening variables for exponential softening take the following

form qu = N f A e 1 -exp -ξu (N f /A e ) G u f qv = V y, f A e 1 -exp -ξv (V y, f /A e ) G v f qw = V z, f A e 1 -exp -ξw (V z, f /A e ) G w f (4.39)
where G u f , G v f , G w f are the fracture energies, and ξu , ξv , ξw are the strain-like softening variables. The stress resultant values in the bulk part of an element in the softening phase, which determine the value of internal forces at the discontinuity t u , t v , t w are computed as The Timoshenko beam finite element mass matrix M e uu in Equation (4.43) is obtained by distributing the total mass of an element to nodes (Figure 4.8), resulting in a diagonally lumped mass matrix, written as where V t ot is obtained as one-third of the total volume of the cylinder with a base diameter equal to d e and height equal to L e (Figure 4.8).

N ′ = E A e ε -ε v p + Ḡα u V ′ y = k c G A e γ y -γ v p y + Ḡα v V ′ z = k c G A e γ z -γ v p z + Ḡα w
M e uu = 1 2 ρ d i ag V t ot , V t ot , V t
3 The solution in terms of unknown nodal displacements and pore pressures is computed at discrete pseudo-time steps t 1 , t 2 , ..., t n using Newmark time-integration scheme. For a time step t n+1 and iteration i , the global system of equations to be solved is written as

V tot = A e •L e •1 L e d e
n el em A e=1    1 β∆t 2 M uu + Kuu -K up γ β∆t K T up γ β∆t D pp + K pp    e,(i ) n+1      ∆ ū ∆ p‵      e,(i ) n+1 =      r u r p      e,(i ) n+1 (4.47)
where Ke uu is the element tangent stiffness matrix, and r e,(i ) u,n+1 and r e,(i ) p,n+1 are residuals pertaining to the solid and the pore fluid part.

After solving the global system of equations, the new iterative values of unknown fields are updated as

ū(i+1) n+1 = ū(i) n+1 + ∆ ū(i) n+1 p(i+1) n+1 = p(i) n+1 + ∆ p(i) n+1 (4.48)
The element tangent stiffness matrix Ke depends whether the element is in elasto-viscoplastic or softening part of the response. If the element is elastoviscoplastic, the tangent stiffness matrix is defined as Ke,(i)

uu,n+1 = K e,(i ) uu,n+1 = L e 0 B s,T u C ep,(i ) n+1 B s u d x (4.49)
Elasto-plastic tangent matrix is written as If the element is in the softening, the element tangent stiffness matrix is obtained by performing static condensation procedure. The statically condensed element tangent stiffness matrix is written as Ke,(i) uu,n+1 = Ke,(i) uu,n+1 = K e,(i ) uu,n+1 -F e,(i ) n+1 H e,(i The matrices K d and K α depend on the current step in softening being elastic or plastic. If the current step in the softening is elastic, then

C ep,(i ) n+1 =                 C ep,
) n+1 + K α -1 F e,(i ),T n+1 + K d (4.
K d = C * B s u ; K α = 0; C * =                 E A e 0 0 0 0 0 0 k c G A e 0 0 0 0 0 0 k c G A e 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 (4.53)
Else, if the current step in the softening is plastic, then

K d = 0; K α =                 K u α A e 0 0 0 0 0 0 K v α A e 0 0 0 0 0 0 K w α A e 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 (4.54)
The local phase procedure for computing matrices K u α , K v α and K w α is given in Chapter 2, Section 2.4.4.2.

External fluid model

The external fluid model is extended from 2D to the 3D setting. The strong form of governing equations and variational formulation is equivalent to that presented in Chapter 3, Section 3.2.2. The end result of the finite element discretization procedure is the following system of equations governing the small and irrotational motion of the external fluid, written as

n el em A e=1            A uu 0 0 0 0 0 0 0 0       e            ü p λ           e +       0 L up L uλ L T up L pp 0 L T uλ 0 L λλ       e            ū p λ           e =            f f 0 0            e      (4.55)
where

A e uu = Ω e f ρN f ,T u N f u d Ω; L e up = - Ω e f V f ,T N f p d Ω; L e uλ = Ω e f D f ,T N f λ d Ω L e pp = - Ω e f 1 β N f ,T p N f p d Ω; L e λλ = - Ω e f 1 ϑ N f ,T λ N f λ d Ω (4.56)
The 'vorticity moment' degree of freedom can again be statically condensed on the element level so that the only unknown variables remaining on the global level are displacements and pressures.

The finite element approximations for displacements, pressure and 'vorticity moment' are given as

u = N f u ū; p = N f p p; Λ = N f λ λ u = N f u ü; p = N f p p; Λ = N f λ λ ∇ • u = (∇ • N f u ) ū = V f ū ∇ × u = (∇ × N f u ) ū = D f ū (4.57)
In the 3D numerical model of acoustic fluid-structure interaction, for the external fluid domain, we choose linear finite element approximations (hexahedral finite element with eight nodes) for the displacements, with constant approximations for the pressure and the 'vorticity moment'. In other words, we choose the Hex8-P1-Λ1 finite element depicted in Figure 4.9. In order to solve the issue of fluid-structure interface, we use the same finite element approximation post-processing procedure previously described in Chapter 3, Section 3.2.3 for the Q4-P1-Λ1 finite element.

Numerical results

In this section, we present the results of several numerical simulations, which serve to illustrate the coupled discrete beam lattice model performance and the capabilities of proposed numerical approach for dealing with acoustic fluid-structure interaction problems in the 3D framework. With the aim of validating and comparing results, the numerical simulations are performed for the numerical examples presented in Chapters 2 and 3, which are extended by adding the third dimension to the model. All numerical implementations and computations are performed with the research version of the computer code FEAP, developed by R.L. Taylor [START_REF] Zienkiewicz | The Finite Element Method[END_REF].

In all numerical simulations, the finite element mesh is generated by using Delaunay tetrahedralization and GMSH software [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The cross-sectional properties of the Timoshenko beam finite elements are computed from the Voronoi diagram by using MATLAB software [START_REF]version 8.5.0.197613[END_REF], which uses Qhull code [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF].

Saturated poro-elastic column

In this section, we aim to validate the coupled discrete beam lattice model of structure built of a saturated porous medium in the 3D framework. We observe saturated poro-elastic column shown in Figure 4.10a, and we perform a quasi-static onedimensional consolidation test. We compare our results against those computed with a continuum model using commercial software PLAXIS [START_REF] Swolfs | PLAXIS 2D Manuals[END_REF]. In PLAXIS, the coupling between the solid phase and the pore fluid is governed by Biot's porous media theory, and the equal order of interpolation for both the displacement and pressure fields is available. For such comparison, we choose 10-node tetrahedral elements with second-order polynomial interpolations. Hence, it can be easily identified from standard experimental tests.

Compression test

In this section, we perform a compression test on a water-saturated specimen shown in

E = 20000 ν = 0.2 σ y,t = 10 µ G,t = 12 G f ,t = 20 σ G,t = 1 σ y,c = 100 µ G,c = 120 G f ,c = 100 σ G,c = 10 σ y,s = 12 µ G,s = 15 G f ,s = 100 σ G,s = 1
The linear hardening modulus is H l h = 2 • 10 3 kPa, the linear kinematic hardening modulus is H l k = 2 • 10 3 kPa, nonlinear hardening parameter is

H nl k = 2 • 10 2 ,
and the viscosity parameter is η = 20 kPa s.

We first perform uncoupled computation (Biot's constant b = 0 → matrix K e up is a null matrix), to simulate the response of dry specimen. We compare the results computed for a 3D model of specimen against those computed for the 2D model of the specimen. For both cases, we present results in terms of average stress, which is obtained as total vertical reaction divided by the area of the specimen (100 mm x 1 mm for 2D case, and 100 mm x 100 mm for 3D case). The comparison of results is shown in Figure 4.15a. We can conclude that the 3D model results in a higher ultimate value of average stress, and higher value of fracture energy needed to drive the specimen to a complete failure.

Dam-reservoir system

In this numerical example, we observe the 3D problem of dam-reservoir interaction. The geometry of the dam-reservoir system is shown in Figure 4.17 Next, we perform localized failure analysis with the aim to compute admissible horizontal overload (in x direction) that can be applied on top of (already) acting self-weight and hydrostatic loading. The computation of admissible overload in the quasi-static setting is performed by imposing horizontal displacements on the top base of the dam. The admissible overload in a dynamic setting is computed by subjecting the dam-reservoir system to the linear increasing horizontal ground acceleration, which is modeled in terms of equivalent horizontal forces. Throughout this numerical example, the vorticity constraint is completely omitted.

→ Hydrodynamic pressure distribution

We subject the dam-reservoir system to a horizontal ground acceleration increasing linearly with time and reaching its maximum value of a 0 = 1 m/s 2 at t = 1s, after which is kept constant (

Concluding remarks

In this Chapter, we presented a 3D numerical model of acoustic fluid-structure interaction. The proposed model is an extension of the previously presented 2D numerical model, with modifications that arise due to the added dimension of the model.

First, we have simulated the response of the pore-saturated structure. The results of numerical simulations in the linear elastic regime of structure response confirm that the coefficient of permeability of coupled discrete beam lattice model of pore-saturated structure closely matches that of an equivalent continuum model.

The results of numerical simulations in the nonlinear regime of structure response

Introduction

In this Chapter, we account for the thermal effects in the structure built of a saturated porous medium. Temperature changes can induce additional stresses in the structure as a result of restrained movement, which may contribute to cracking.

Furthermore, the thermal actions can affect the stiffness and strength properties of the material, which has to be properly accounted for in order to provide a sound design of the structure. To be able to model these phenomena, which influence the overall safety and durability of the structure, we extend the coupled discrete beam lattice model presented in Chapter 4, for the isothermal case to nonisothermal case by introducing the thermal coupling into the model.

The outline of the Chapter is as follows: In Section 5.2, we describe the equations governing the response of the nonisothermal saturated porous medium. In Section 5.3, we present the finite element formulation of the proposed thermo-hydromechanical coupled discrete beam lattice model, together with the computational procedure. The finite element formulation is given in fully dynamics framework, which can be easily adjusted for quasi-static simplification. In Section 5.4, we present the results of several numerical simulations with the aim to validate the proposed coupled discrete beam lattice model of structure built of a nonisothermal saturated porous medium, and to demonstrate its capabilities to predict response and cracking in the structure subjected to combined thermal and mechanical loads.

In Section 5.5, we give concluding remarks.

Governing equations of nonisothermal saturated porous medium

The equations governing the response of nonisothermal saturated porous medium are derived by combining equilibrium equation imposed on a porous medium, continuity equation imposed on a fluid flow and energy equation imposed on heat flow through such a porous medium [START_REF] Aboustit | Variational principles and finite element simulations for thermo-elastic consolidation[END_REF][START_REF] Aboustit | Finite element invesitagtions of thermo-elastic and thermoplastic consolidation[END_REF][START_REF] Lewis | A coupled finite element model for the consolidation of nonisothermal elastoplastic porous media[END_REF]5]. In a dynamic setting, we refer to the equation of motion instead of the equilibrium equation. The isothermal equations of motion and continuity equation are explained in detail in Chapter 2, Section 2.2.

Next, we present the equations governing the response of the nonisothermal saturated porous medium. The governing equations are derived under the assumptions that no phase change occurs and that the thermal equilibrium between the solid phase and the pore fluid is achieved.

Energy equation

The energy equation for porous medium under the assumption that both the solid phase and pore fluid have the same temperature at the coincident point and with convection ignored is written as

ρC T Ṫ + ∇ • q T -s = 0 (5.6)
where ρC T is the effective heat capacity defined as

ρC T = (1 -n)ρ s C s T + nρ f C f T , with C s
T as the heat capacity of the solid phase and C f T as the heat capacity of the pore fluid, q T is the heat flux and s is the heat source.

The heat flux q T is defined with Fourier law for heat conduction, written as

q T = -k T ∇T (5.7)
where k T is the coefficient of thermal conductivity of the isotropic porous medium.

In the formulation of energy equation the convection, the pressure and temperature dependence of the solid and fluid densities, and mechanical contributions to energy balance are neglected. These simplifications are in many cases justified and result in an uncoupled form of the energy equation [5,[START_REF] Lewis | A coupled finite element model for the consolidation of nonisothermal elastoplastic porous media[END_REF][START_REF] Booker | Consolidation around a point heat source[END_REF].

Discrete model of structure

The starting point for our work on coupled discrete beam lattice model of the nonisothermal pore-saturated structure is the model presented in Chapter 4, for the isothermal case. Next, we present the modifications in the finite element formulation that arise due to the introduction of thermal coupling in the model.

Finite element formulation

Kinematics

The kinematics for inelastic Timoshenko beam finite element, and finite element approximation for pore pressure field are given in Chapter 4, Section 4.3.1. The temperature field is approximated in the same manner as the pore pressure field, using tetrahedral finite elements with four nodes (Figure 5.1).

The finite element approximation for temperature field is written as

T = N s T T (5.8)
where TT = T 1 , T 2 , T 3 , T 4 ;

N s T = N t 1 , N t 2 , N t 3 , N t 4 ; N t i = a i + b i x + c i y + d i z 6V , i = 1, 4; 6V = d et 1 x 1 y 1 z 1 1 x 2 y 2 z 2 1 x 3 y 3 z 3 1 x 4 y 4 z 4 ; a 1 = d et x 2 y 2 z 2 x 3 y 3 z 3 x 4 y 4 z 4 ; b 1 = -d et 1 y 2 z 2 1 y 3 z 3 1 y 4 z 4 ; c 1 = -d et x 2 1 z 2 x 3 1 z 3 x 4 1 z 4 ; d 1 = -d et x 2 y 2 1
x 3 y 3 1

x 4 y 4 1

; et c.

(5.9)

where V is the volume of the tetrahedral element, x,y,z are global coordinates and

M ṗ + b ε -βs f Ṫ -∇ • k γ f ∇p = 0 (5.11)
By introducing finite element approximations into the weak form of continuity equation and performing standard finite element discretization procedure (Chapter 4, Section 4.3), we obtain

Ω e Tet 4 N s,T p 1 M N s p ṗ -N s,T p βs f N s T Ṫ + (∇N s p ) T k γ f ∇N s p p d Ω + L e 0 N s,T up bB s up ud x = q ext ,e
(5.12)

Energy equation

The energy equation is written as

ρC T Ṫ -∇ • k T ∇T -s = 0 (5.13)
The weak form of energy equation is written as

Ω e Tet 4 δT ρC T Ṫ -k T ∇ • (∇T ) -s d Ω = 0 (5.14)
where δT is the virtual temperature field interpolated in the same manner as the real temperature field δT = N s T δ T; δ TT = δT 1 , δT 2 , δT 3 , δT 4 ;

p ‵ = N s uT T‵ ; N s uT = N 1 , N 2 ; T‵T = T 1 , T 2 ; δT ‵ = N s uT δ T‵ ; δ T‵T = δT 1 , δT 2 (5.15)
By introducing finite element approximations into the weak form of continuity equation and performing standard finite element discretization procedure, we

obtain

Ω e Tet 4 N s,T T ρC T N s T Ṫ + (∇N s T ) T k T ∇N s T T d Ω = s ext ,e
(5.16)

Equations of motion

The derivation of the weak form of equations of motion for the isothermal porous medium is presented in detail in Chapter 4, Section 4.3.2.2. The resulting weak form of equations of motion for the nonisothermal porous medium is written as f acc,e + f i nt ,ef ext ,e = 0 (5.17

)

h e = 0 (5.18)
where f i nt ,e is the vector of inertial forces, f i nt ,e is the internal force vector, f acc,e is the external load vector, and h e is the residual vector due to discontinuity, computed as

f acc,e = L e 0 N s,T u σd x; f i nt ,e = L e 0 B s,T u σd x = f i nt ,e (σ uT ) -f i nt ,e (σ p ); h e = L e 0 G T σ uT d x = L e 0 ḠT σ uT d x + t; Ḡ =                 Ḡ 0 0 0 0 0 0 Ḡ 0 0 0 0 0 0 Ḡ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 ; t T = t u , t v , t w , 0 , 0 , 0 (5.19)
By taking into account that the coupling between the mechanics with the pore fluid flow and heat flow occurs through the axial direction of the Timoshenko beam finite element, the vectors in (5.19) are computed as σT = ρ A e ü, ρ A e v, ρ A e ẅ, ρI e pol θx , ρI e y y θy , ρI e zz θz ;

σ T = N ,V y ,V z , M x , M y , M z ; σ T = σ T uT -σ T p ; σ T uT = N ′ uT ,V ′ uT,y ,V ′ uT,z , M ′ uT,x , M ′ uT,y , M ′ uT,z ; σ T p = bp ‵ A e , 0, 0, 0, 0, 0 ; σ T uT = σ T u -σ T T ; σ T u = N ′ ,V ′ y ,V ′ z , M ′ x , M ′ y , M ′ z ; σ T T = β uT T (T ‵ -T 0 )A e , 0, 0, 0, 0, 0 (5.20)
where β uT T is the thermal stress for thermo-mechanical coupling in the axial direction. The stress resultant vector σ T u results from displacements, and is computed in the same manner as described in Chapter 4, Section 4.3. The internal force vector t acting at the discontinuity is computed as

t =                 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 L e 0 σ uT d x (5.21)
For linear elasticity, the thermal stress β uT T is a function of lattice Young's modulus E and thermal expansion coefficient of the solid phase βs , written as

β uT T = Ek β uT T βs (5.22)
The coefficient k β uT T takes into account the difference between 1D thermomechanical coupling and the 3D thermo-mechanical coupling. Namely, in 3D thermo-mechanical coupling for isotropic case the thermal stress β T can be expressed in terms of thermal expansion coefficient βs and Lame's parameters for continuum model λ L and µ L , written as

β T = (3λ L + 2µ L ) βs (5.23)
The last expression can be rewritten in terms of oedometer modulus E oed and Poisson's ratio ν for continuum model, as

β T = 1 + ν 1 -ν E oed βs (5.24)
Hence, the coefficient k β uT T can then be computed as

k β uT T = 1 + ν 1 -ν E oed /E (5.25)

Computational procedure

The system of equations governing the thermo-hydro-mechanical coupled problem for single Timoshenko beam finite element is written as where M e uu is the mass matrix, f i nt ,e (σ uT ) is the internal load vector resulting from displacements and temperature changes, K e up and K e pT are the coupling matrices, D e pp is the compressibility matrix, K e pp is the permeability matrix, D e T T is the heat capacity matrix, K e T T is the conductivity matrix, and f e,ext , q e,ext , s

M e uu ü + f i nt ,e (σ uT ) -K e up p‵ = f ext ,e (5
Here we note that in the global system of equations (Equations 5.27 and 5.28)

we take the parts of the D e pp , K e pp , K e pT , D e T T and K e T T matrices (Equation 5.29) that correspond to the nodes of Timoshenko beam finite element.

The solution in terms of unknown nodal displacements, pore pressures and temperatures is computed at discrete pseudo-time steps t 1 , t 2 , ..., t n using Newmark time-integration scheme. For a time step t n+1 and iteration i , the global system of equations to be solved is written as

n el em A e=1       1 β∆t 2 M uu + Kuu -K up -KuT γ β∆t K T up γ β∆t D pp + K pp γ β∆t K pT 0 0 γ β∆t D T T + K T T       e,(i ) n+1            ∆ ū ∆ p‵ ∆ T‵            e,(i ) n+1 =            r u r p r T            e,(i ) n+1 (5.30)
where Ke uu is the tangent stiffness matrix, Ke uT is the tangent coupling matrix , and r e,(i ) u,n+1 , r e,(i ) p,n+1 and r e,(i ) T,n+1 are residuals pertaining to the solid, the pore fluid and the temperature part.

After solving the global system of equations, the new iterative values of unknown fields are updated as

ū(i+1) n+1 = ū(i) n+1 + ∆ ū(i) n+1 p(i+1) n+1 = p(i) n+1 + ∆ p(i) n+1 T(i+1) n+1 = T(i) n+1 + ∆ T(i) n+1
(5.31)

The tangent stiffness matrix Ke uu and tangent coupling matrix Ke uT depend whether the element is in elasto-viscoplastic or softening part of the response. If the element is elasto-viscoplastic, the tangent stiffness matrix and tangent coupling matrix are defined as Ke,(i)

uu,n+1 = K e,(i ) uu,n+1 = ∂f i nt (σ uT ) δ ū e,(i ) n+1 = L e 0 B s,T u C ep,(i ) n+1 B s u d x Ke,(i) uT,n+1 = K e,(i ) uT,n+1 = ∂f i nt (σ uT ) δ T e,(i ) n+1 = L e 0 B s,T uT C ep,u,(i ) n+1 k β uT T βs A e N s uT d x;
B s uT = B 1 0 0 0 0 0 B 2 0 0 0 0 0 (5.32)

Elasto-plastic tangent matrix is written as If the element is in the softening, the tangent stiffness matrix and tangent coupling matrix are obtained by performing static condensation procedure. The statically condensed tangent stiffness matrix is written as Ke,(i) uu,n+1 = Ke,(i) uu,n+1 = K e,(i ) uu,n+1 -F e,(i ) n+1 H e,(i

C ep,(i ) n+1 =                 C ep,u n+1 A e 0 0 0 0 0 0 C ep,v n+1 A e 0 0 0 0 0 0 C ep,w n+1 A e 0 0 0 0 0 0 G I e pol 0 0 0 0 0 0 E I e 0 0 0 0 0 0 E I e                 (
) n+1 + K α -1
F e,(i ),T n+1 + K d (5.34) where

K e,(i ) uu,n+1 = ∂f i nt (σ uT ) δ ū e,(i ) n+1 = L e 0 B s,T u C ep,(i ) n+1 B s u d x F e,(i ) n+1 = ∂f i nt (σ uT ) δα e,(i ) n+1 = L e 0 B s,T u C ep,(i ) n+1 Ḡd x F T + K d e,(i ) n+1 = ∂h δ ū e,(i ) n+1 = L e 0 ḠT C ep,(i ) n+1 B s u d x + K d H + K α e,(i ) n+1 = ∂h δα e,(i ) n+1 = L e 0 ḠT C ep,(i ) n+1 Ḡd x + K α (5.35)
The matrices K d and K α depend on the current step in softening being elastic or plastic. If the current step in the softening is elastic, then

K d = C * B s u ; K α = 0; C * =                 E A e 0 0 0 0 0 0 k c G A e 0 0 0 0 0 0 k c G A e 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 (5.36)
Else, if the current step in the softening is plastic, then

K d = 0; K α =                 K u α A e 0 0 0 0 0 0 K v α A e 0 0 0 0 0 0 K w α A e 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                
(5.37)

The local phase procedure for computing matrices K u α , K v α and K w α is given in Section 2.4.4.2.

By taking into account that the coupling of solid phase with fluid flow occurs through the axial direction of the Timoshenko beam finite element, the statically condensed tangent coupling matrix is written as Ke,(i)

uT,n+1 = Ke,(i)

uT,n+1 = K e,(i ) uT,n+1 -Q e,(i ) n+1 Y e,(i ) n+1 + K αT -1 W e,(i ) n+1 + K T (5.38)
where

K e,(i ) uT,n+1 = ∂f i nt (σ uT ) δ T e,(i ) n+1 = L e 0 B s,T uT C ep,u,(i ) n+1 k β uT T βs A e N s uT d x Q e,(i ) n+1 = ∂f i nt (σ uT ) δα e,(i ) n+1 = L e 0 B s,T uT C ep,u,(i ) n+1 k β uT T βs A e Ḡd x W + K T e,(i ) n+1 = ∂h δT e,(i ) n+1 = L e 0 ḠC ep,u,(i ) n+1 k β uT T βs A e N s uT d x + K T Y + K αT e,(i ) n+1 = ∂h δα e,(i ) n+1 = L e 0 ḠC ep,u,(i ) n+1 A e Ḡd x + K αT (5.39)
The matrices K T and K αT depend on the current step in softening being elastic or plastic. If the current step in the softening is elastic, then

K T = C ep,u,(i ) n+1
A e N s uT ; K αT = 0 (5.40)

Else, if the current step in the softening is plastic, then

K T = 0; K αT = K u α A e
(5.41)

Numerical results

In this section, we present the results of several numerical simulations, which serve to illustrate the proposed approach ability to deal with nonisothermal problems. All numerical implementations and computations are performed with the research version of the computer code FEAP, developed by R.L. Taylor [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. In all numerical simulations, the finite element mesh is generated by using Delaunay tetrahedralization and GMSH software [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The cross-sectional properties of the Timoshenko beam finite elements are computed from Voronoi diagram by using MATLAB software [START_REF]version 8.5.0.197613[END_REF], which uses Qhull code [START_REF] Barber | The quickhull algorithm for convex hulls[END_REF].

Nonisothermal saturated poro-elastic column

In this section, we perform a quasi-static thermo-elastic one-dimensional consolidation analysis of saturated column, with the aim to validate the coupled discrete beam lattice model of structure built of saturated nonisothermal porous medium.

The problem of one-dimensional thermo-elastic consolidation was first solved by Aboustit et al. [START_REF] Aboustit | Variational principles and finite element simulations for thermo-elastic consolidation[END_REF][START_REF] Aboustit | Finite element invesitagtions of thermo-elastic and thermoplastic consolidation[END_REF], and is later used as a benchmark by Lewis et al. [START_REF] Lewis | A coupled finite element model for the consolidation of nonisothermal elastoplastic porous media[END_REF],

Noorishad and Tsang [START_REF] Noorishad | Coupled thermohydroelasticity phenomena in variably saturated fractured porous rocks-formulation and numerical solution[END_REF], and Cui et al. [START_REF] Cui | An alternative coupled thermo-hydro-mechanical finite element formulation for fully saturated soils[END_REF] to test their finite element formulations in 2D framework. To test our model, we extend the problem analyzed in the literature for 2D case (Figure 5.2a) to the 3D case (Figure 5.2b) and we compare our results against reference solutions provided by Lewis et al. [START_REF] Lewis | A coupled finite element model for the consolidation of nonisothermal elastoplastic porous media[END_REF] and Cui et al. [START_REF] Cui | An alternative coupled thermo-hydro-mechanical finite element formulation for fully saturated soils[END_REF].

Because the problem analyzed is the problem of one-dimensional consolidation, the results computed for the 3D case should match those obtained for the 2D case.

The problem analyzed is unit-less. 5.1 [START_REF] Lewis | A coupled finite element model for the consolidation of nonisothermal elastoplastic porous media[END_REF][START_REF] Cui | An alternative coupled thermo-hydro-mechanical finite element formulation for fully saturated soils[END_REF]. We perform three types of computations. The first is the isothermal consolidation with the applied surface pressure only, the second is the thermal consolidation with the applied surface temperature only, and the third is the thermal consolidation with applied both the surface temperature and the surface pressure.

The computed results in terms of vertical displacements of column top for all three cases are shown in Figures 5.4a-5.4c. We can conclude that a good match between the computed results and reference values is obtained. The computed results are practically mesh independent. We note that the results for the thermal consolidation with the applied surface temperature only are not given for comparison in [START_REF] Cui | An alternative coupled thermo-hydro-mechanical finite element formulation for fully saturated soils[END_REF]. and reference solutions, with results being practically mesh independent. Next, we perform a thermo-mechanical computation (the pore pressure field is excluded from the computation) for the case of applied surface temperature only, and we compare the computed results with those obtained with the 3D numerical model of a column constructed in FEAP by using FEAP built-in solid elements [START_REF] Zienkiewicz | The Finite Element Method[END_REF]. The temperature evolution computed with discrete model is shown in Figures 5.7a

Thermal mechanical compression test

In this section, we perform a thermal mechanical compression test. We first subject the specimen to high temperatures, after which we impose vertical displacements on the top base of the specimen. We assume that the inertial effects are negligible.

The geometry of the specimen and the boundary conditions are shown in Figure To investigate the temperature effects on the stiffness and load carrying capacity of the specimen, we assume that the material parameters: Young's modulus, yield limits, fracture limits and fracture energies of Timoshenko beam finite element are temperature dependent. We assume linear temperature dependence, written as m = m 1 -ω T (T -T 0 ) , T > T 0 (5.42) where m is the material parameter, T and T 0 are the temperature and the reference temperature in the element and ω T is the parameter controlling the temperature dependence. The value of ω T depends on the type of material and has to be identified through more elaborate parameter identification procedures. For example, the exposure of the rock or concrete specimen to the elevated temperatures results in a decrease in the compressive strength, whereas for clay specimen the compressive strength increases.

Dam-reservoir system

In concrete dams, the temperature gradient near the exposed surface can cause surface cracks, which can affect structural durability. These surface cracks can also result in the long-term degradation of the load carrying capacity and stiffness of the dam, especially due to the large seasonal ambient temperature variations.

Hence, the thermal effects in the structure have to be evaluated in order to set the initial conditions for the further failure analyses and to ensure the durability of the exposed surfaces and structure in general. We note that in this numerical example, we assumed the temperature dependence of the material parameters of dam material. However, this assumption can be omitted without significant influence on the computed results, because the temperature in the dam varies in a relatively narrow range.

The computed temperature distribution in the body of the dam is shown in Figure 5.13. We can conclude that the temperature variation (Figure 5.12b) affect the temperature less than 0.5 m into a cross-section, which is due to the high heat capacity and relatively low thermal conductivity of the dam material. The proposed numerical analysis can be further extended to account for the loading history throughout the dam operating period, which allows for a more precise investigation of the current state and remaining resistance of the existing dams.

The safety analysis of this kind involves more complex loading program, which includes ambient temperature variations over the dam operating period. Namely, the large seasonal ambient temperature changes contribute to the long-term degradation in dam structure. Hence, aside from the self-weight and hydrostatic loading of the external reservoir, which are considered as the fundamental load cases, the effects of the large seasonal ambient temperature changes have to be computed in order to set the initial conditions for the computation of remaining resistance. The ambient temperatures over a period of months or years, depending on the dam operating period, are obtained from in-situ data measurements, and can include daily temperature variations, the weekly averages or can be approximated by sinusoidal function if the data follow sinusoidal pattern (Figure 5.16, for illustrative purposes not result of real data measurements). To more precisely model the effects of ambient temperature changes on the response of dam structure, convective heat transfer as a result of temperature differences between exposed surfaces and ambient temperature should be added in the numerical model, as well as radiation.

The self-weight, hydrostatic loading, and long-term temperature effects define initial conditions for further safety analysis, which can also include accidental loads such as earthquake or fire that can occur during a dam lifespan. From the standpoint of engineering practice, it is essential to be able to quantify the remaining resistance of the structure, which previously experienced considerable damage from extreme loads of this kind. To decrease computational cost, this kind of numerical analysis would require, operator split solution procedure with multi-scale analysis in time can be utilized [START_REF] Ibrahimbegovic | Thermomechanical coupling in folded plates and non-smooth shells[END_REF][START_REF] Niekamp | Formulation, solution and ctl software for coupled thermomechanics systems[END_REF]. Namely, a large time step can be selected for the thermal part followed by a number of smaller time steps for the mechanical part, where the time step during mechanical computation can be controlled with the damage growth.

Concluding remarks

In this Chapter, we presented a coupled discrete beam lattice model of structure built of a nonisothermal saturated porous medium, which is an extension of the coupled discrete beam lattice model previously presented in Chapter 4, for the isothermal case.

We have confirmed through the results of validation computations in the linear elastic regime of structure response that the coefficient of thermal conductivity, the effective heat capacity and the thermal expansion coefficient of coupled discrete beam lattice model closely match those of an equivalent continuum model. Hence, they can be easily identified from standard experimental tests.

Next, we have illustrated the temperature effects on the stiffness and load carrying capacity of the structure by implementing the linear temperature dependence of the material parameters. The parameter controlling the temperature dependence has to be identified through more elaborate parameter identification procedures.

Finally, we have simulated the response and cracking of the dam structure subjected to combined mechanical and thermal loads. This type of numerical analysis can be further extended to failure analyses, that can also account for the long-term effects of the large seasonal ambient temperature changes on the integrity and durability of the dam structure.

CHAPTER 6 CONCLUSIONS AND PERSPECTIVES Conclusions

In this thesis, we studied the issue of the overall safety of structures built of saturated poro-plastic cohesive material in application to acoustic fluid-structure interaction problems, such as dam-reservoir interaction. We presented a novel approach to numerical modeling of dam-reservoir interaction, with the emphasis on the numerical representations of the structure and the external fluid.

The first feature of the proposed approach is the numerical representation of the structure in terms of a saturated porous medium. The macro-scale response and ultimate failure modes in the pore-saturated structure are captured with a coupled discrete beam lattice model based on Voronoi cell representation of the domain with cohesive links as inelastic Timoshenko beam finite elements enhanced with additional kinematics in terms of embedded strong discontinuities. The internal interaction of solid phase and pore fluid is handled with Biot's porous media theory and Darcy's law for pore fluid flow, resulting with an additional pressure-type degree of freedom placed at each node of a Timoshenko beam finite element.

The second feature of the proposed approach is the representation of the small motion of external fluid in confined conditions with the acoustic wave theory, Lagrangian formulation and mixed displacement/pressure finite element approximation, resulting with the displacement and pressure degrees of freedom. The finite element representations of the pore-saturated structure and the external fluid allows for the structure and the external fluid finite elements to share the displacement and pressure degrees of freedom. Thus, the issue of the fluid-structure interface is solved in an efficient and straightforward manner by directly connecting the structure and external fluid finite elements at common nodes.

  and professor Ivica Kožar for careful reading and constructive comments. I am deeply grateful to emeritus professor Muhamed Zlatar for his constant support from the very beginning of my academic and research career. To madam Catherine Palpant from the Department of Culture and Scientific Cooperation, French Embassy in Bosnia and Herzegovina and madam Željka Divković from Campus France office in Bosnia and Herzegovina, for their immense support.

Keywords:

  cohesive discrete beam lattice model, Timoshenko beam, embedded discontinuity, coupling, pressure, porous medium, Biot's theory, acoustic wave theory, mixed displacement/pressure approximation, fluid-structure interaction, failure analysis, overall safety, temperature xiv RÉSUMÉ Titre: Analyse des mécanismes de rupture des systèmes structurels en poroplasticité Cette thèse porte sur la sécurité globale des structures en matériaux hétérogènes saturés soumis à des charges extrêmes, et est appliquée à des problèmes d'interaction fluide-structure, tels que l'interaction barrage-réservoir. Un modèle numérique d'interaction est proposé pour prédire les principales tendances et le comportement général d'un barrage en matériau saturé en interaction avec le réservoir dans des analyses de défaillance d'intérêt pratique. Le modèle numérique proposé est d'abord présenté dans un cadre bidimensionnel (2D), puis étendu à un cadre tridimensionnel (3D).

(Figure 1 . 1 Figure 1 . 2

 1112 Photo by R. Couture (CGC), 2017)

Figure 1 . 3

 13 Figure 1.3 Failure in dam structure

  For the numerical representation of the pore-saturated structure, we propose the use of coupled discrete beam lattice, which is based on Voronoi cell representation of the domain with inelastic Timoshenko beam finite elements enhanced with additional kinematics in terms of embedded strong discontinuities acting as cohesive links. The internal fluid-structure interaction, i.e., the coupling between the solid phase and pore fluid is introduced through Biot's porous media theory and Darcy's law governing the pore fluid flow, ultimately resulting with an additional pressure-type degree of freedom placed at each node of Timoshenko beam finite element. Next, we focus on numerical modeling of external fluid and finally fluid-structure interaction. We assume a small motion of the external fluid in confined conditions typically met for fluid placed in a reservoir. We describe the small motion of the external fluid with the acoustic wave theory. We model the external fluid motion with the Lagrangian formulation and mixed displacement/pressure approximation. The finite element approximation for external fluid in combination with coupled Introduction discrete beam lattice model for pore-saturated structure allows for external fluid and structure finite elements to share both the displacement and pressure degrees of freedom at the fluid-structure interface. Thus, the issue of the fluid-structure interface in the numerical model of interaction is solved in a straightforward manner through the element assembly procedure. The thesis is divided into six chapters. In Chapter 1, we outline the aim of the thesis, together with the motivation for research work and literature review. In Chapter 2, we describe in detail the plane strain coupled discrete beam lattice model of the pore-saturated structure, and we discuss the results of several numerical examples. In Chapter 3, we present the numerical model of external fluid and ultimately numerical model of fluid-structure interaction, and its application through a set of numerical simulations. In Chapter 4, we extend the proposed numerical model of the structure, external fluid, and acoustic fluid-structure interaction to the 3D setting. In Chapter 5, we account for the temperature effects by introducing thermal coupling in the proposed 3D numerical model of the poresaturated structure. Finally, in Chapter 6, we summarize all the main findings and contributions of the thesis together with the suggestions for future works. All numerical implementations and computations are performed with the research version of the computer code FEAP (Finite Element Analysis Program),

For the numerical representationFigure 2 . 1

 21 Figure 2.1 Voronoi diagram and Delaunay triangulation

Figure 2 . 2 Figure 2 . 3

 2223 Figure 2.2 Discretization of displacement and pore pressure fields

Figure 2 . 5

 25 Figure 2.5 Timoshenko beam finite element: interpolation functions

Figure 2 . 6

 26 Figure 2.6 CST finite element

L e 0 δuρ( 2

 02 A e ü + δvρ A e v + δθρI e θ d x = δ ūT f ext ,eδ ūT f acc,e u σd x; σT = ρ A e ü, ρ A e v, ρI e θ (2.30)

Figure 2 . 7

 27 Figure 2.7 Eccentricity of the center of the gravity

(2. 50 )

 50 where γu , γv are the plastic multipliers whose values are obtained from

Figure 2 . 8

 28 Figure 2.8 Timoshenko beam finite element lumped mass matrix computation

2. 4 . 4 . 1 •

 441 Local phase: bulk Computational procedure for proposed plasticity model is as follows Start with the best iterative values of nodal displacements, then • Assume elastic trial step γu,tr i al n+1 Calculate trial values of resultant forces, stress-like hardening variables and back-stress variables

Remark 3 :

 3 In quasi-static setting, with inertial effects being negligible, the resulting set of equations on the global level are the first order differential equations in time. The global phase of computation under quasi-static assumptions is solved using Newton's iterative method and backward Euler scheme which solves the first order transient problem in the following manner[START_REF] Zienkiewicz | The Finite Element Method[END_REF]

2 .Figure 2 .

 22 Figure 2.9 Nonlinear cantilever beam

Figure 2 .

 2 Figure 2.11 Nonlinear column

1 - 10 Figure 2 .

 1102 Figure 2.13 Nonlinear column: computed hysteresis loops

Figure 2 .

 2 Figure 2.14 Uniaxial tension test

2 .

 2 14a. The finite element mesh and Voronoi diagram are shown in Figures 2.14b

and 2 .

 2 14c. The material parameters of Timoshenko beam finite elements are shown in Table

Figure 2 .

 2 Figure 2.15 PDF of Gaussian random distribution of fracture limit in compression

Figure 2 .Figure 2 .

 22 Figure 2.16 Uniaxial tension test: computed results

Figure 2 .

 2 Figure 2.19 Saturated poro-elastic column

and 2 .

 2 21b. We can conclude that the results obtained with the discrete model show an excellent agreement with the results obtained with a continuum model.The computed results are practically mesh independent. What is important to emphasize is that the coefficient of permeability of a coupled discrete beam lattice model matches that of an equivalent continuum model. Hence, it can be easily identified from standard experimental tests.

  Figure 2.21 Saturated poro-elastic column: computed results

  Figure 2.23 Comparison of loading rates

Figure 2 .Figure 2 .

 22 Figure 2.24 Footing on soil stratum

Figure 2 .

 2 Figure 2.28 Computed excess pore pressure at point D(0,6)

  , 2.27b and 2.28b.

  The excess pore pressure fields at t = 5 days computed with the continuum and discrete model are shown in Figures 2.29a

Figure 2 .

 2 Figure 2.29 Nonlinear behavior, computed excess pore pressure field

λ L = 5 .Time 2 ]

 52 Figure 2.30 Saturated poro-elastic column, dynamic case

Figure 2 .

 2 Figure 2.32 Saturated poro-elastic column, dynamic case: computed results

  .33a. The finite element mesh and Voronoi discretization are shown in Figures 2.33b

  and 2.33c.

  Figure 2.33 Plane strain compression test

Figure 2 .Figure 2 .Figure 2 .

 222 Figure 2.34 Computed results, different displacement rates

Figure 2 .Figure 2 .Figure 2 . 2 . 6

 22226 Figure 2.39 Computed results, different displacement rates

  the proposed approach capabilities to carry on with computations of the external fluid interacting with structure undergoing localized failure with cracks and to represent internal fluid (pore) pressure influence on ultimate failure modes with fully saturated cracks. Namely, we generalize the development of the coupled discrete beam lattice model described in Chapter 2, extending its ability to deal with fluid-structure interaction with the emphasis on the numerical representation of the structure response in terms of a saturated porous medium. We limit ourselves to fluid-structure problems under conditions that enable the modeling of the external fluid with the acoustic wave theory. Namely, the latter is justified for the small motion of the external fluid in confined conditions typically met for fluid placed in a reservoir. For the numerical representation of the external fluid motion, we choose a Lagrangian formulation and mixed displacement/pressure based finite element approximation. The coupled discrete model for the structure built of a saturated porous medium in combination with the mixed displacement/pressure based formulation for the external fluid allows for the structure and fluid finite elements to share both the displacement and pressure degrees of freedom at the fluid-structure interface. This permits for the elements to be connected directly at the common nodes without any special numerical treatment of the fluid-structure interface. As a result, the exchange of both the motion and the pressure at the fluid-structure interface is ensured, and all computations can be performed in a fully monolithic manner.The main benefit of the proposed model is its ability to provide a better numerical representation of the structure in interaction with confined fluids. In previous attempts to model the fluid-structure interaction problems of this kind, the use of pressure degree of freedom at the fluid-structure interface led to a cumbersome exchange of information in terms of added mass. Here, the numerical model of the structure contains pore pressure as an additional degree of freedom per node, which enables direct exchange of the pressures at the fluid-structure interface. With such a numerical model of the structure, we are able to simulate the full saturation of the material of the structure at every time step, with the external fluid acting as the source of pore saturation. In addition, we are able to take into account the inelastic behavior, to predict the crack formation and propagation in saturated materials and quantify overall safety of pore-saturated structures subjected to the extreme loads.The outline of the Chapter is as follows: In Section 3.2, we present equations governing the external fluid motion derived from acoustic wave theory under small motion hypothesis. We describe mixed displacement/pressure based finite element formulation, and we provide two different finite elements, for the quasistatic case and dynamic case. In Section 3.3, we present the results of several numerical simulations. We simulate acoustic fluid-structure interaction problems in a quasi-static and dynamic setting, respectively. First, we perform validation computations of the proposed external fluid model. Second, we investigate the proposed numerical model of acoustic fluid-structure interaction capabilities to predict hydrostatic and hydrodynamic pressure distributions on the upstream face of the structure, which are in close agreement with analytical solutions. Finally, we focus on the proposed numerical model ability to deal with the loss of structural integrity, ultimately providing the overall safety assessment of the dam structure to potential overload with respect to the fundamental load cases. In Section 3.4, we give concluding remarks.

  ticity moment'. In other words, for the discretization of the external fluid domain, we use the Q4-P1-Λ1 finite element depicted in Figure 3.1.

Figure 3 . 1

 31 Figure 3.1 Finite element approximation for external fluid, dynamic case: Q4-P1-Λ1

p c = p 2 Figure 3 . 2

 232 Figure 3.2 Extrapolation of the pressures

Figure 3 . 3

 33 Figure 3.3 Finite element approximation for external fluid, quasi-static case: Q4-P4-Λ4

Figure 3 . 4 Figure 3 . 5 Figure 3 . 6

 343536 Figure 3.4 Water container

7 .

 7 The finite element mesh and Voronoi discretization of the structure domain are shown in Figures 3.8a and 3.8b. The loading program is shown in Figure 3.8c. The specific weight of the dam material is γ s = 20 kN/m 3 , and the specific weight of the water is γ w = 10 kN/m 3 .

Figure 3 . 7 Figure 3 . 8

 3738 Figure 3.7 Geometry of dam-reservoir system, quasi-static case

Figure 3 . 9

 39 Figure 3.9 Linear elastic behavior: computed results

Figure 3 .

 3 Figure 3.10 Linear elastic behavior: computed results

  σ y,t = 0.015 MPa; σ y,c = 0.20 MPa; σ y,s = 0.015 MPa; σ f ,t = 0.02 MPa; σ f ,c = 0.30 MPa; σ f ,s = 0.02 MPa. The fracture energies in tension, compression and shear are: G f ,t = 0.01 GN/m; G f ,c = 0.1 GN/m; G f ,s = 0.01 GN/m. The linear hardening modulus is H l h = 10 3 MPa.

Figure 3 .Figure 3 .

 33 Figure 3.11 Inelastic behavior: broken cohesive links and plastic zones

and 3 .

 3 11c, broken cohesive links in the dam at the end of the loading program are shown. These broken cohesive links (marked red) represent Timoshenko beam finite elements that had entered softening phase of the response in mode I, and mode II. The plastic zones and broken cohesive links result in the less stiff response of the dam. This leads to an increase in the value of horizontal displacement of the tip of the dam, compared to the linear elastic case (Figure3.12).

First, we compute0Figure 3 .Figure 3 .Figure 3 .Figure 3 .

 3333 Figure 3.13 Horizontal overload: computed results

Figure 3 .

 3 Figure3.19 Dam-reservoir system studied by Westergaard[START_REF] Westergaard | Water pressures on dams during earthquakes[END_REF]

555ρa 0 H 2 ( 3 . 16 )→Figure 3

 23163 Figure3.20 Dam-reservoir system studied by[START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 1. momentum method[END_REF] and[START_REF] Chwang | Hydrodynamic pressures on sloping dams during earthquakes. part 2. exact theory[END_REF]

Figure 3 .Figure 3 .Figure 3 . 3 . 3 . 2 . 2 Figure 3 .

 33333223 Figure 3.21 Pressure coefficient provided by [149] and [150], Adapted from "Hydrodynamic pressures on sloping dams during earthquakes. Part 2. Exact theory", by A.T. Chwang, 1978, The Journal of Fluid Mechanics, 87(2), pp. 347. Copyright 1978 Cambridge University Press, Adapted with permission

Figure 3 .

 3 Figure 3.25 Linear elastic separator wall: problem statement

Figure 3 . 8 Figure 3 .

 383 Figure 3.26 Linear elastic separator wall: loading program

Figure 3 .

 3 Figure 3.28 Pressure and pore pressure distribution: L/H = 4

Figure 3 . 10 kN/m 3 .Figure 3 .

 31033 Figure 3.29 Geometry of dam-reservoir system, dynamic case

and 3 .Figure 3 .Figure 3 .Figure 3 .Figure 3 .

 33333 Figure 3.31 Linear elastic behavior: computed results, Θ = 90 •

Figure 3 .Figure 3 .

 33 Figure 3.35 Horizontal overload, dynamic case: computed results

Figure 3 .

 3 Figure 3.37 Comparison of horizontal overloads, dynamic case

2. 3 .

 3 Namely, the construction of the model is performed by exploiting the duality property between the Voronoi cell representation and Delaunay tetrahedralization of the domain (Figure 4.1). The end result of Delaunay tetrahedralization is the mesh of tetrahedra. Every edge of tetrahedra connects the centers of two adjacent

Figure 4 . 1

 41 Figure 4.1 Voronoi diagram and Delaunay tetrahedralization

Figure 4 . 2 Figure 4 . 3

 4243 Figure 4.2 Displacement and pore pressure fields finite element approximations, 3D case

Figure 4 . 4

 44 Figure 4.4 Timoshenko beam finite element in 3D framework

Figure 4 . 5

 45 Figure 4.5 Three modes of crack formation

Figure 4 . 6

 46 Figure 4.6 Linear tetrahedral finite element

Figure 4 . 7

 47 Figure 4.7 Stress resultants

(4. 35 ) 3 .z

 353 The evolution equations for internal variables ϵ v p , γ and ξu , ξv , ξw εvp = γu si g n(N ′ -χu A e ); ξu = γu ; γu ≥ 0 γvp y = γv si g n(V ′ y -χv A e ); ξv = γv ; γv ≥ 0 γvp z = γw si g n(V ′ z -χw A e ); ξw = γw ; γw ≥ 0 (4.36)

(4. 40 ) 2 .

 402 The evolution equations for internal variables α u , α v , α w and ξu , ξv , ξw with the loading/unloading conditions αu = γu si g n(t u ); ξu = γu ; γu ≥ 0; φu ≤ 0; γu φu = 0 αv = γv si g n(t v ); ξv = γv ; γv ≥ 0; φv ≤ 0; γv φv = 0 αw = γw si g n(t w ); ξw = γw ; γw ≥ 0; φw ≤ 0; γw φw = 0(4.41) where γu , γv , γw are the plastic multipliers whose values are obtained from

Figure 4 . 8

 48 Figure 4.8 Timoshenko beam finite element lumped mass matrix computation in 3D setting

  are elasto-plastic tangent moduli for axial and transverse directions. The local phase procedure for computing elasto-plastic tangent moduli is given in Section 2.4.4.1.

Figure 4 . 9

 49 Figure 4.9 Finite element approximation for external fluid in 3D numerical model Hex8-P1-Λ1

Figure 4 .

 4 Figure 4.10 3D Saturated poro-elastic column

Figure 4 .

 4 Figure 4.11 Finite element mesh densities

and 4 .

 4 12b. We can conclude that the results obtained with the discrete model show an excellent agreement with the results obtained with a continuum model, with the results being practically mesh independent.

Figure 4 . 12 3D

 412 Figure 4.12 3D Saturated poro-elastic column: computed results

and 4 .

 4 13b. We compare computed results against those shown in Figures 4.12a and 4.12b. From Figure4.13b, we can conclude that a higher loading rate result in an increase in the value of excess pore pressure. As the value of excess pore pressure approaches zero, the value of vertical displacement approaches constant value marking the end of the consolidation.

Figure 4 .

 4 Figure 4.13 Comparison of loading rates

Figure 4 .Figure 4 .

 44 Figure 4.14 3D Compression test

Figure 4 .

 4 Figure 4.15 Comparison of computed results

Figure 4 .

 4 Figure 4.17 Geometry of 3D dam-reservoir system

Remark 1 :

 1 The finite element mesh of the structure domain is obtained by performing Delaunay tetrahedralization of the domain, which results with the irregular triangular mesh of Timoshenko beam finite elements on the domain surfaces. With the goal of connecting the structure finite elements with the external fluid finite elements, we impose regular triangular mesh on the fluid-structure boundary. By doing so, we are able to connect structure and external fluid finite elements at the common nodes(Figure 4.18). It is important to note that since we are connecting line elements with volume elements at the fluid-structure boundary, we ensure the exchange of motion and pressure through the common nodes, whereas the compatibility across the surface is not achieved.

Figure 4 . 2 : 4 . 5 . 3 . 1 →

 424531 Figure 4.18 Mesh of Timoshenko beam finite elements at fluid-structure boundary

Figure 4 .

 4 Figure 4.19 Loading program for validation analysis, quasi-static case

Figure 4 .

 4 Figure 4.20 Hydrostatic pressure and pore pressure distribution

Figure 4 .

 4 Figure 4.21 Linear elastic behavior, dynamic case: computed results

Figure 4 .Figure 4 .Figure 4 .

 444 Figure 4.22 Loading program for validation analysis, dynamic case

Figure 4 .

 4 Figure 4.24 Hydrodynamic pressure and pore pressure distribution

Figure 4 .

 4 Figure 4.25 Localized failure analysis, loading programs

  and 4.25b. We first apply an increase of the self-weight, followed by the hydrostatic loading. The time step in the first two computation phases is set to ∆t = 0.1 day. In a quasi-static setting, we compute the admissible horizontal overload by imposing horizontal displacements with a constant rate of 1 • 10 -3 mm/s on the top base of the dam. The time step is set to ∆t = 1 s. In a dynamic setting, we compute the admissible horizontal overload by subjecting the dam-reservoir system to the linear increasing horizontal ground acceleration with a constant rate of 2 m/s 2 /1s. The initial time step is set to ∆t = 0.001 s. Prior to reaching the ultimate value of horizontal overload, the time step is reduced to ∆t = 0.0001 s.→ Quasi-static caseThe computed admissible horizontal overload force is 101.2393 kN(16.8732 kN/m') (Figure4.26). The total horizontal reaction resulting from the self-weight and hydrostatic loading is 749.58 kN (124.93 kN/m'). Hence, we can state that for this particular geometry and the mechanical properties of the dam, the factor of safety of the dam against failure is 1.14. The broken cohesive links in increasing softening at the end of the loading program are shown inFigures 4.27a

  -4.27c.

Figure 4 .

 4 Figure 4.26 Horizontal overload (x direction), quasi-static case

  -4.29c.

Figure 4 .Figure 4 .

 44 Figure 4.28 Horizontal overload (x direction), dynamic case

  x i , y i , z i are nodal coordinates of tetrahedral element, and constants a 2-4 , b 2-4 , c 2-4 , d 2-4 are defined by cyclic interchange of the subscripts in the order 1, 2, 3, 4.

Figure 5 . 1 3 . 2 . 1

 51321 Figure 5.1 Linear tetrahedral finite element

  are elasto-plastic tangent moduli for axial and transverse directions. The local phase procedure for computing elasto-plastic tangent moduli is given in Section 2.4.4.1.

Figure 5 . 2

 52 Figure 5.2 Thermo-elastic consolidation of saturated column, problem statement

Figure 5 . 3

 53 Figure 5.3 Finite element mesh densities

Figure 5 . 4

 54 Figure 5.4 Saturated poro-elastic column: vertical displacement of column top

Figure 5 . 5

 55 Figure 5.5 Thermal case, applied surface temperature and pressure: computed results

  -5.7c.

Figure 5 . 6 Figure 5 . 7

 5657 Figure 5.6 Thermo-mechanical computation, applied surface temperature:comparison with 3D numerical model in FEAP

5 .

 5 8a. The finite element mesh is shown in Figure 5.8b. The material parameters of the Timoshenko beam finite element are shown in Table 5.2. The linear hardening modulus is H l h = 2 • 10 3 kPa, the linear kinematic hardening modulus is H l k = 2 • 10 3 kPa, nonlinear hardening parameter is H nl k = 2 • 10 2 , and the viscosity parameter is η = 20 kPa s. The coefficient of permeability is k = 10 -8 m/s, the specific weight of the water is γ w = 10 kN/m 3 , Biot's constant is b = 1 and Biot's modulus is 1/M → 0. The coefficient of thermal conductivity is k T = 2 W/mK, the effective heat capacity is ρC T = 1850 kJ/Km 3 , the thermal expansion coefficient of the solid phase is βs = 0.00001 ( • C ) -1 and the reference temperature is T 0 = 0 • . The temperature term in the continuity equation is omitted, i.e. βs f is set to 0. The coefficient k β ut T is selected as 1.0.

Figure 5 . 8 2 σ

 582 Figure 5.8 Thermal mechanical compression test

Figure 5 . 9 Figure 5 .

 595 Figure 5.9 Thermal mechanical compression test: computed results

Figure 5 .Figure 5 .

 55 Figure 5.11 Geometry of dam-reservoir system

Figure 5 .Figure 5 .

 55 Figure 5.13 Temperature evolution in the dam

Figure 5 .

 5 Figure 5.15 Horizontal displacements, x direction

Figure 5 .

 5 Figure 5.16 Sinusoidal ambient temperature change, illustrative example

Discrete model of structure built of saturated poro-plastic medium

	(a) Young's modulus	(b) Poisson's ratio
	Figure 2.4 Numerical tests for identification of linear elastic parameters of discrete
	lattice model: boundary conditions
	Next, we give an overview of the finite element formulation of the plane strain
	coupled discrete beam lattice model.
	uum, linear elastic parameters of the discrete lattice model (lattice Young's modulus
	and Poisson's ratio) can be easily identified from standard experimental tests in
	combination with numerical tests, shown in Figures 2.4a and 2.4b.

Table 2 .

 2 1 Oedometer vs. Young's modulus for different values of Poisson's ratio

	Poisson's ratio Oedometer modulus/Young's modulus
	ν	E oed /E
	0	1
	0.1	1.02
	0.2	1.11
	0.3	1.35
	0.4	2.14
	0.45	3.79

Table 2 .

 2

	Young's	Poisson's	Yield limit	Hardening	Fracture	Fracture
	modulus	ratio		modulus	limit	energy
	[MPa]	[/]	[MPa]	[MPa]	[MPa]	[N/mm]

2 Uniaxial tension test: material parameters of the finite element

Table 2 .

 2 3 Compression test: material parameters of the finite element

	Young's	Poisson's	Yield limit	Hardening	Fracture	Fracture
	modulus	ratio		modulus	limit	energy
	[kPa]	[/]	[kPa]	[kPa]	[kPa]	[kN/m]

Table 2 .

 2 4 Coupled soil-foundation system: material parameters of the finite element

	Young's	Poisson's	Yield limit	Hardening	Fracture	Fracture
	modulus	ratio		modulus	limit	energy
	[MPa]	[/]	[MPa]	[MPa]	[MPa]	[GN/m]

  • 10 -5 m/s is 236.15s, and for loading rate v = 2 • 10 -5 m/s is 556.34s (Processor: Intel Core i7-3632QM/2.2GHz, RAM 8.00GB, FEAP: single-core).

	Table 2.5 Coupled soil foundation-system (k = 1 • 10 -6 m/s, v = 5 • 10 -5 m/s ),
	convergence rates for time step: 4000
	Iteration Residual norm Energy norm
	1	6.7516 • 10 -03	4.9926 • 10 -07
	2	3.9441 • 10 -03	1.6061 • 10 -07
	3	6.3604 • 10 -04	4.0685 • 10 -09
	4	4.0008 • 10 -05	9.1335 • 10 -11
	5	4.3019 • 10 -06	1.0264 • 10 -13
	6	2.8904 • 10 -11	3.5076 • 10 -23

OUTLINE 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 External fluid model in mixed Lagrangian formulation . . . . . . . . . 3.2.1 Assumptions and governing equations . . . . . . . . . . . . . . 3.2.2 Mixed displacement/pressure based formulation . . . . . . . . 3.2.3 Finite element approximation post-processing . . . . . . . . . 3.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Quasi-static case . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1.1 Rigid and flexible water container . . . . . . . . . . . . 3.3.1.2 Dam-reservoir system . . . . . . . . . . . . . . . . . . . 3.3.2 Dynamic case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2.1 Analytical solutions for hydrodynamic pressures . . . 3.3.2.2 Modal analysis of rigid cavity problem . . . . . . . . . 3.3.2.3 Linear elastic separator wall . . . . . . . . . . . . . . . 3.3.2.4 Dam-reservoir system . . . . . . . . . . . . . . . . . . . 3.4 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 3 .

 3

	1 Rigid cavity problem: computed frequencies
	Mesh density		Frequency [Hz]
	4x3	174.4	374.9	444.6	463.0
	8x6	171.1	348.8	429.9	459.3
	32x24	170.1	340.5	425.3	457.8
	64x48	170.0	340.1	425.1	457.8
	Analytical solution	170.0	340.0	425.0	457.7

Table 3 .

 3 2 Convergence rates: residual norm for selected time steps

	Iteration	Self-weight phase Hydrostatic loading phase Localized failure Time step: 50 Time step: 180 Time step: 1610
	1	2.4177 • 10 -04	3.4008 • 10 -04	1.2834 • 10 -02
	2	5.4542 • 10 -04	4.6168 • 10 -04	8.8185 • 10 -04
	3	2.0887 • 10 -05	1.7170 • 10 -05	4.4117 • 10 -04
	4	8.7336 • 10 -16	1.6796 • 10 -13	1.3242 • 10 -04
	5			2.5074 • 10 -05
	6			1.0504 • 10 -04
	7			3.9459 • 10 -07
	8			2.0083 • 10 -13

Table 3 .

 3

		3 Convergence rates: energy norm for selected time steps
	Iteration	Self-weight phase Hydrostatic loading phase Localized failure Time step: 50 Time step: 180 Time step: 1610
	1	9.7579 • 10 -10	8.2718 • 10 -09	2.6056 • 10 -08
	2	7.6616 • 10 -11	9.7831 • 10 -11	5.8175 • 10 -10
	3	8.2985 • 10 -14	2.7183 • 10 -14	1.4944 • 10 -10
	4	3.5316 • 10 -35	1.0109 • 10 -29	5.8688 • 10 -12
	5			7.2201 • 10 -13
	6			6.4927 • 10 -13
	7			2.3475 • 10 -17
	8			6.3199 • 10 -29

  ot , I e pol , I e , I e , V t ot , V t ot , V t ot , I e pol , I e , I e ) (4.46)

Table 4 .

 4 1 3D Compression test: material parameters of the finite element

	Young's	Poisson's ratio	Yield limit	Fracture limit	Fracture
	modulus				energy
	[kPa]	[/]	[kPa]	[kPa]	[kN/m]

Table 5 .

 5

		1 Time stepping scheme
	Time interval Number of time steps
	0.01	10
	0.1	10
	10	10
	100	10
	1000	20

CHAPTER 3 ACOUSTIC FLUID-STRUCTURE INTERACTION ABSTRACT

In this Chapter, we propose a novel approach for dealing with acoustic fluidstructure interaction nonlinear problems, such as progressive localized failure of a dam structure built of porous cohesive material in interaction with reservoir under extreme static and/or dynamic loads. The main original aspects of the proposed approach concern ability to deal with the loss of structural integrity in fluid-structure interaction problems, as well as a seamless discretization between structure and external fluid achieved by using a judicious combination of a lattice model for the pore-saturated structure and finite element approximation for external fluid. Namely, the interaction of structure built of a poro-plastic medium is here restricted to acoustic fluids, which allows to use Lagrangian formulation for the small motion of the external fluid and mixed discrete approximation resulting with the same displacement pressure degrees of freedom and seamless approximation between structure and fluid. The latter is justified for the small motion of the external fluid confined in the reservoir and simplification of the acoustic wave theory for irrotational flows and quasi-incompressible fluids. The application of the model is illustrated through several numerical simulations.

CHAPTER 4 EXTENSION TO 3D CASE ABSTRACT

In this Chapter, we extend the proposed two-dimensional (2D) numerical model of acoustic fluid-structure interaction to the three-dimensional (3D) setting. The 3D numerical model of interaction inherits all the main features and benefits of the 2D model, with modifications that arise due to the added dimension of the model. The structure response is described with coupled discrete beam lattice model of structure built of saturated poro-plastic cohesive material. The small motion of external fluid motion in contained conditions is described with acoustic wave theory and modeled with the Lagrangian formulation and mixed displacement/pressure based finite element approximation. The finite element representations of the structure and the external fluid allow for the issue of the fluid-structure interface to be solved through element assembly procedure. The application of the model is illustrated through several numerical simulations.

2. Strain energy function in terms of strains and internal variables, plastic deformations ϵ v p , γ v p y , γ v p z and strain-like hardening variables ξu , ξv , ξw

where for circular cross-sections shear correction factor is equal to k c = 9/10.

3. Yield function in terms of stresses (or stress resultants), stress-like hardening variables qu , qv , qw and back-stress variables χu , χv , χw

where N y , V y,y , V z,y are the elastic (yield) limits.

The stress resultant values N ′ , V ′ y , V ′ z , M ′ x , M ′ y and M ′ z are computed from the elastic part of total deformation. For a circular cross-section for which the element local coordinate axis x coincides with the center of gravity, the stress resultant values are computed as

Computational procedure

The finite element discretization procedure reduces the nonlinear coupled problem to two sets of equations. The first is the first order differential evolution equations of internal variables defined locally at each Gauss point, and the second is the second order differential equations governing the coupled problem defined globally. The solution of these equations is again computed by using the operator split solution procedure, in the same manner as described in Chapter 2, Section 2.4.4. The solution procedure is split into local and global phase, which are treated separately.

In the local phase, the solution of evolution equations is computed by using implicit backward Euler time integration scheme (see Sections 2.4.4.1 and 2.4.4.2). In the global phase, the solution in terms of the unknown nodal displacements and pore pressures is computed in a fully monolithic manner using Newmark timeintegration scheme and Newton's iterative method.

The system of equations governing the coupled problem for single Timoshenko beam finite element is written as

)

where M e uu is the mass matrix, f i nt ,e (σ u ) is the internal load vector resulting from displacements (i.e. effective stresses), K e up is the coupling matrix, D e pp is the compressibility matrix, K e pp is the permeability matrix, and f e,ext and q e,ext are the load vectors. The matrices K e up , D e pp and K e pp are computed as

Here we note that in the global system of equations (Equation 4.44) we take the parts of the D e pp , and K e pp matrices (Equation 4.45) that correspond to the nodes of Timoshenko beam finite element.

For numerical integration on tetrahedral finite elements, we choose the nodal point integration rule. The nodal point integration rule is a linear order rule, which results with a diagonal form of compressibility matrix D e pp . However, this does not affect significantly computed results. Namely, the Biot's modulus M in the problems of our interest is always a large number or infinity, which results with the values of entries of a matrix D e pp close to zero or zero.

demonstrate the capabilities of the proposed numerical model to capture inelastic response and localized failure in the pore-saturated structure.

Second, we have simulated dam-reservoir interaction. We have confirmed that the proposed numerical model of acoustic fluid-structure interaction is capable of predicting the hydrostatic and hydrodynamic pressure distributions exerted on the upstream face of the dam structure in close agreement with analytical solutions.

Furthermore, we have illustrated the proposed numerical model ability to deal with the progressive localized failure of dam structure by computing the admissible horizontal overload that dam structure can withstand on top of already acting self-weight and hydrostatic pressure.

The proposed 3D numerical model of the structure and numerical model of acoustic fluid-structure interaction exhibit good performance in terms of fast convergence rates. However, the computational time is significantly increased when compared to the 2D case. For instance, the computational time for horizontal overload with included first two computation phases for the 2D case is 720.83 s, and for the 3D case 232566.9 s.

CHAPTER 5 EXTENSION TO NONISOTHERMAL CASE ABSTRACT

In this Chapter, we account for the thermal effects in the structure built of a saturated porous medium. Namely, we present a coupled discrete beam lattice model of structure built of a nonisothermal saturated porous medium, which is an extension of the coupled discrete beam lattice model proposed in Chapter 4, for the isothermal case. The modifications pertain to the introduction of thermal coupling in the model. The temperature field is approximated with linear tetrahedral finite elements, resulting with an additional temperature degree of freedom placed at each node of the Timoshenko beam finite element. The heat flow in the model is governed by Fourier law for heat conduction. The application of the model is illustrated through several numerical simulations.

Equations of motion

The strong form of equations of motion is written as

where σ is the total stress tensor, b are the body forces vector, a is the solid phase acceleration vector and ρ is the mass density of the mixture assumed to be constant. Following Terzaghi's principle of effective stresses, the total stress is decomposed into effective stress and pore pressure, written as

where σ ′ is the effective stress tensor, I is the second order identity tensor, p is the pore pressure assumed positive in compression, and b is Biot's constant.

Furthermore, the effective stress tensor can be decomposed into mechanical part σ u resulting from displacements and thermal part σ T resulting from changes in temperature, written as

with σ T computed as

where β T is the thermal stress tensor for isotropic case defined as β T = β T I, and T 0 is the reference temperature.

Continuity equation

The continuity equation for fluid flow through a nonisothermal porous medium is written as 1

where M is Biot's modulus, k is the coefficient of permeability of the isotropic porous medium, γ f is the specific weight of the fluid and βs f is the thermal expansion coefficient of the mixture defined as βs f = (bn) βs + n βf , with βs as the thermal expansion coefficient of the solid phase and βf as the thermal expansion coefficient of the pore fluid. In the formulation of continuity equation, the temperature dependence of the solid density is neglected.

Conclusions and Perspectives

The proposed numerical models of pore-saturated structure, external fluid and ultimately numerical model of acoustic fluid-structure interaction are first presented in 2D framework, and later extended to the 3D framework. The proposed numerical model of the pore-saturated structure in the 3D framework is further extended to the nonisothermal case.

We performed several numerical simulations with the aim to test the potential applications of the proposed approach. The successful application of discrete lattice models depends on the appropriate modeling of the micro-scale, which is in

Perspectives

The proposed approach for acoustic fluid-structure interaction can have practical importance and can be used in various situations in engineering practice. The proposed numerical model of fluid-structure interaction can be used for predicting the safety and durability of existing dams in Bosnia and Herzegovina while taking into account existing defects or weak spots in the dam structure. It could also serve as a tool for numerical computations in the structural analysis and design of new dams and reservoirs.

Future works on the presented research include

• Testing the application of the 3D model of acoustic fluid-structure interaction in failure analyses of dam structures with complex geometries.

• Failure analyses of dam structure with complex loading programs involving large seasonal ambient temperature variations and extreme mechanical loads.

• Implementation of probability based identification procedures for determining the material parameters of Timoshenko beam finite element, such as yield limits, fracture limits and fracture energies, which would enable the comparison of numerical results with in-situ data measurements of real dam structures.

• Development of the coupled discrete beam lattice model of the structure built of a partially saturated porous medium.

LIST OF FIGURES