​

​

The objective of this work is the derivation and implementation of a unified Finite

Element formulation for the solution of fluid and solid mechanics, Fluid-Structure Interaction

(FSI) and coupled thermal problems.

​

The unified procedure is based on a stabilized velocity-pressure Lagrangian formulation. Each time

step increment is solved using a two-step Gauss-Seidel scheme: first the linear momentum equations

are solved for the velocity increments, next the continuity equation

is solved for the pressure in the updated configuration.

​

The Particle Finite Element Method (PFEM) is used for the fluid domains, while the

Finite Element Method (FEM) is employed for the solid ones. As a consequence, the domain is

remeshed only in the parts occupied by the fluid.

​

Linear shape functions are used for both the velocity and the pressure fields. In order to deal

with the incompressibility of the materials, the formulation has been stabilized using an updated

version of the Finite Calculus (FIC) method. The procedure has been derived for

quasi-incompressible Newtonian fluids. In this work, the FIC stabilization procedure has been

extended also to the analysis of quasi-incompressible hypoelastic

solids.

Specific attention has been given to the study of free surface flow problems. In particular,

the mass preservation feature of the PFEM-FIC stabilized procedure has been deeply studied with the

help of several numerical examples. Furthermore, the conditioning of the problem has been analyzed

in detail describing the effect of the bulk modulus on the numerical scheme. A strategy based on

the use of a pseudo bulk modulus for improving the conditioning of the linear system is also

presented.

​

The unified formulation has been validated by comparing its numerical results to experimental

tests and other numerical solutions for fluid and solid mechanics, and FSI problems. The

convergence of the scheme has been also analyzed for most of the prob-

lems presented.

The unified formulation has been coupled with the heat tranfer problem using a staggered

scheme. A simple algorithm for simulating phase change problems is also described. The numerical

solution of several FSI problems involving the temperature is given.

​

The thermal coupled scheme has been used successfully for the solution of an industrial problem.

The objective of study was to analyze the damage of a nuclear power plant pressure vessel induced

by a high viscous fluid at high temperature, the corium.  The

numerical study of this industrial problem has been included in this work.

​

==Resumen==

​

El objectivo de la presente tesis es la derivación e implementación de una formulación unificada con elementos finitos para la solución de problemas de mecánica de fluidos y de sólidos, interacción fluido-estructura (''Fluid-Structure Interaction (FSI)'') y con acoplamiento térmico.  El método unificado está basado en una formulación Lagrangiana estabilizada y las variables incçognitas son las velocidades y la presión. Cada paso de tiempo se soluciona a través de un esquema de dos pasos de tipo Gauss-Seidel. Primero se resuelven  las ecuaciones de momento lineal por los incrementos de velocidad, luego  se calculan las presiones  en la configuración actualizada usando la ecuación de continuidad.  Para los dominios fluidos se utiliza el método de elementos finitos de partículas (''Particle Finite Element Method (PFEM)'') mientras que los sólidos  se solucionan con el método de elementos finitos  (''Finite Element Method (FEM)''). Por lo tanto, se ramalla sólo las partes del dominio ocupadas por el fluido.  Los campos de velocidad y presión se interpolan con funciones de forma lineales. Para poder analizar materiales incompresibles, la formulación ha sido estabilizada con una nueva versión del método ''Finite Calculus (FIC)''. La técnica de estabilización ha sido derivada para fluidos Newtonianos casi-incompresibles. En este trabajo, la estabilización con FIC se usa también para el análisis de sólidos hipoelásticos casi-incompresibles.  En la tesis se dedica particular atención al estudio de flujo con superficie libre. En particular, se analiza en profundidad el tema de las pérdidas de masa y se muestra con varios ejemplos numéricos la capacidad del método de garantizar la conservación de masa en problemas de flujos en superficie libre. Además se estudia con detalle el condicionamiento del  esquema numérico analizando particularmente el efecto del módulo de compresibilidad. Se presenta también una estrategia basada en el uso de un pseudo módulo de compresibilidad para mejorar el condicionamiento del problema.  La formulación unificada ha sido validada comparando sus resultados numéricos con pruebas de laboratorio y resultados numéricos de otras formulaciones. En la mayoría de los ejemplos también se ha estudiado la convergencia del método.  En la tesis también se describe una estrategia segregada para el acoplamiento de la formulación unificada con el problema de transmisión de calor. Además se presenta una simple estrategia para simular el cambio de fase. El esquema acoplado ha sido utilizado para resolver varios problemas de FSI donde se incluye la temperatura y su efecto.  El esquema acoplado con el problema térmico ha sido utilizado con éxito para resolver un problema industrial. El objetivo del estudio era la simulación del daño y la fusión de la vasija de un reactor nuclear provocados por el contacto con un fluido altamente viscoso y a gran temperatura. En la tesis se describe con detalle el estudio numérico realizado para esta aplicación industrial.

​

​

==Acknowledgements==

​

I would like to thank all the people that in some way have helped me to complete this doctoral thesis.   I would like to express my gratitude to Prof. Eugenio Oñate for many reasons. Thanks to him, four years ago I could come to Barcelona and start my PhD at CIMNE. He gave me the balanced doses of freedom and pressure and he encouraged me to test innovative technologies. I would like thank Prof. Oñate also for giving me the opportunity and the material support to spend a three-month period at the College of Engineering of Swansea.  Special thanks go to Dr Josep Maria Carbonell.  Without his support I could not realize this work. He devoted me innumerable hours of his time and he gave me technical and not technical suggestions that helped me to make the right choices during my PhD. Moltes gràcies Josep Maria!  I would like to thank Prof. Javier Bonet for giving me the possibility to work with his group during my stay at Swansea. Special thanks go to Prof. Antonio Gil y Dr. Aurelio Arranz for their very kind treatment, for improving my knowledge of the Immersed Boundary Potential Method and for the useful discussions about FSI strategies. I really hope that this is just the first step for a future fruitful collaboration.  I would like to thank also Mr. Kazuto Yamamura and the Nippon Steel & Sumimoto Metal Corporation for giving us the possibility to publish the PFEM simulation  of the damages of a nuclear power plant pressure vessel.  I also acknowledge the Agencia de Gestió d'Ajuts Universitaris i de Recerca (AGAUR) for the financial support.  I would like to thank also Riccardo, Antonia, Jordi, Pablo and Lorenzo for their contributions for this thesis.   This thesis is dedicated to my parents that still represent to me a fundamental and unique support. Grazie per tutti gli sforzi che avete fatto in questi anni e per le numerose visite a Barcellona. Questa tesi è per voi.  Finally, I would like to thank Lucía simply for staying by my side during these years. Her smile has been the gasoline I used for dealing with all this work.

​

​

==Abbreviations==

​

'''AL'''  '''A'''ugmented '''L'''agrangian 

​

'''ASGS'''  '''A'''lgebraic '''S'''ub'''G'''rid '''S'''cale 

​

'''FCM'''  '''F'''uel '''C'''ontaining '''M'''aterial

​

'''FEM'''   '''F'''inite '''E'''lement  '''M'''ethod 

​

'''FIC'''  '''F'''inite '''I'''ncrement '''C'''alculus 

​

'''FSI'''  '''F'''luid-'''S'''tructure '''I'''nteraction 

​

'''GLS'''  '''G'''alerkin '''L'''east-'''S'''quares 

​

'''IBM'''  '''I'''mmersed  '''B'''oundary '''M'''ethod 

​

'''IFEM'''  '''I'''mmersed '''F'''inite '''E'''lement  '''M'''ethod 

​

'''ISPM'''  '''I'''mmersed '''S'''tructural  '''P'''otential '''M'''ethod 

​

'''LBB'''  '''L'''adyzenskaya-'''B'''abuzka-'''B'''rezzi 

​

'''LFCM'''  '''L'''ava-like '''F'''uel '''C'''ontaining '''M'''aterial

​

'''NPP'''  '''N'''uclear '''P'''ower '''P'''lant

​

'''OSS'''  '''O'''rthogonal '''S'''ub'''S'''cale 

​

'''PFEM'''  '''P'''article '''F'''inite '''E'''lement  '''M'''ethod 

​

'''SPH'''  '''S'''mooth '''P'''article '''H'''ydrodynamics 

​

'''SUPG'''  '''S'''treamline '''U'''pwind '''P'''etrov '''G'''alerkin  

​

'''TL'''  '''T'''otal '''L'''agrangian

​

'''UL'''  '''U'''pdated '''L'''agrangian

​

'''V'''  '''V'''elocity

​

'''VMS'''  '''V'''ariational '''M'''ulti-'''S'''cale 

​

'''VP'''  '''V'''elocity-'''P'''ressure

​

'''VPS'''  '''V'''elocity-'''P'''ressure-'''S'''tabilized

​

​

​

​

''Dai diamanti non nasce niente, <br />

''dal letame nascono i fior.'' 

​

"Nothing grows out of precious diamonds,  <br />

Out of dung the flowers do grow."<br />

Fabrizio De Andrè

​

​

=Chapter 1. Introduction=

​

​

==1.1 Objectives==

​

The objective of this work is to develop a unified formulation for the solution of the fluid and solid mechanics, fluid-structure interaction (FSI) and thermal coupled problems. The aim is to analyze the continuum in a unified manner trying to reduce at minimum the differences between the analysis of fluids and solids. For this purpose the numerical model has been planned in order to meet the specific requirements of  solid and fluid mechanics and  their approximation with the Finite Element Method (FEM) [136], but without limiting excessively the capability of the model. In fact the computational method should be capable to deal with critical problems such as those involving elasto-plastic solids, quasi-incompressible materials, free surface fluids and phase change due to heat transfer solutions.

​

According to these considerations, the computational model has been designed for using a stabilized Velocity-Pressure formulation. The formulation has been particularized for hypoelasto-plastic, compressible and quasi-incompressible solids and quasi-incompressible Newtonian fluids. The algorithm for the FSI problems has been inspired by the analogous unified strategy presented in [63].  For the fluid phase, the Particle Finite Element Method (PFEM) [88] has been used, while for the solid the classical FEM is adopted. The scheme has been stabilized with an updated version of the Finite Calculus (FIC) technique [82]. For dealing with thermal coupled problems, this unified formulation has been extended to treat the heat transfer problem via a staggered scheme. For the modeling of the phase change, a simplified procedure has been developed.  The whole formulation has been implemented in a C++ code.

​

==1.2 State of the art==

​

In this section, an overview of the numerical methods used for simulating a free surface fluid flow interacting with deformable solids is given. For the sake of clarity, the section is divided in three parts representing the main fields involved by this work. First the Eulerian and Lagrangian approaches for free surface fluid dynamics problems are presented. Then an overview of the stabilization for incompressible, or quasi-incompressible, material is given. Finally, the principal algorithms for solving FSI problems are described.

​

===1.2.1 Eulerian and Lagrangian approaches for free surface flow analysis===

​

Consider the  description of the motion of a general continuum  represented in Figure [[#img-1|1]]. The domain <math display="inline">{\Omega _0}</math> represents the body at the initial state at time <math display="inline">{t=t_0}</math> while the domain <math display="inline">{\Omega }</math> represents the same body at time <math display="inline">{t=t_n}</math> after deformation. The domain <math display="inline">{\Omega _0}</math>  is called ''initial configuration'', whereas <math display="inline">{\Omega }</math> is the ''current,'' or ''deformed'', ''configuration''. In order to describe  the kinematics and the deformation of the body, the ''reference configuration'' has to be defined because the motion is defined with respect to this configuration. <div id='img-1'></div>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"

|-

|[[Image:draft_Samper_722607179-motionDescription.png|600px|Description of the motion of a general continuum body]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 1:''' Description of the motion of a general continuum body

|}

​

In solid mechanics, the stresses generally depend on the history of deformation and the undeformed configuration must be specified in order to define the strains. Due to the history dependence, Lagrangian descriptions are prevalent in solid mechanics. In fluids however, the stresses do not depend on the history and it is often unnecessary to describe the motion with respect to a reference configuration. For this reason an Eulerian description represents the most reasonable choice. Furthermore, in problems where the fluid contours are fixed, Eulerian meshes are generally preferred to the Lagrangian ones. This is because Eulerian grids are fixed and they do not deform according to the fluid motion, as  shown in Figure [[#img-2|2]].  <div id='img-2'></div>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"

|-

|[[Image:draft_Samper_722607179-eulerianMesh.png|480px|Motion description using an Eulerian mesh]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 2:''' Motion description using an Eulerian mesh

|}

​

Conversely, in the Lagrangian description, the mesh nodes coincide with the fluid particles and the discretization moves and deforms as the fluid flow (Figure [[#img-3|3]]).  <div id='img-3'></div>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 80%;max-width: 100%;"

|-

|[[Image:draft_Samper_722607179-lagrangianMesh.png|480px|Motion description using a Lagrangian mesh]]

|- style="text-align: center; font-size: 75%;"

| colspan="1" | '''Figure 3:''' Motion description using a Lagrangian mesh

|}

​

Consequently, on the one hand the non-linear convective term disappears from the problem and on the other hand the mesh undergoes large distortions and it requires to be regenerated.

​

In the analysis of free surface flows, the detection of the free surface contours represents a crucial task. Its position is unknown <math display="inline">a priori</math> and it has to be determined at each time increment in order to solve properly the boundary value problem. For these problems, the Lagrangian description may be preferred to the Eulerian one. In fact, with a Lagrangian approach the free surface is detected automatically by the position of the mesh nodes, while an Eulerian approach requires the implementation of a specific technology for this task.

​

Several strategies have been developed in the literature  for tracking the free surface in an Eulerian framework. One of the earliest contributions was given by the so called ''marker and cell'' method [53]. In this approach a set of marker particles that move according the flow are used to detect which regions are occupied by the fluid and which not. An evolution of this technique is the ''volume of fluid'' method [55]. In this case the free surface boundary is detected using a scalar function that assumes the unit value in the fluid cells and the value zero in those ones with no fluid. The cells with an intermedium value are the ones that contain the free surface.  Another possibility is the ''level set'' method [100]. This technique is used in various fields, not only in continuum mechanics, and it allows for detecting shapes or surfaces on a fixed grid without making any parametrization of them. For this reason, this procedure has been also used for matching the free surface contour on an Eulerian mesh [106]. A similar idea was used in [7] where the position of the free surface is detected using a cloud of Lagrangian particles moving over an Eulerian mesh.

​

The free surface flows can be solved also using an hybrid Eulerian-Lagrangian technique. This is the so termed ''Arbitrary Lagrangian-Eulerian (ALE)'' approach [137]. The aim of this method is to exploit the best features of the Eulerian and the Lagrangian procedures and to combine them. The mesh nodes can arbitrary be fixed or can move with the fluid [38]. Generally, far from the moving boundaries a fixed Eulerian grid is used, while near to the interface the mesh moves according to the motion of the boundary  [25]. However if the boundary motion is large or unpredictable, also in the ALE methods the grid may suffer large distortions and may require a proper remeshing procedure [116].

​

In purely Lagrangian approaches, the  mesh needs to be regenerated whenever a  threshold limit for the distortion is exceeded. This is the basis of a particular class of Lagrangian finite element formulation called the  ''Particle Finite Element Method (PFEM)''. The method was initially developed by the group of professors  Idelsohn and Oñate [66,67]. The PFEM treats the mesh nodes of the domain as particles which can freely move and even separate from the rest of the  fluid domain representing, for instance, the effect of water drops. A mesh connects the nodes discretizing the domain where the governing equations are solved using a classical finite element method. These features make the PFEM the ideal numerical procedure to model and simulate free surface flows. In the last years, many scientific publications have shown the efficacy of the PFEM for solving free surface flow problems, see among others [32,68,117]. The PFEM has also been tested successfully in other kind of problems, such as fluid mechanics including thermal convection-diffusion [5,80,95], multi-fluids  [36,62], granular materials [131], FSI [81,132] and  excavation [19].

​

''Meshfree'' methods are other class of Lagrangian techniques. In this strategy the remeshing is not required because the governing equations  are solved over a set of nodes without referring to a mesh. One of the first meshfree techniques is the ''Smooth Particle Hydrodynamics (SPH)'' method. This method was introduced independently by Gingold and Monaghan [50] and Lucy [70] for the simulation of astrophysical problems such as fission of stars. SPH is a particle-based Lagrangian technique where discrete smoothed particles are used to compute approximate values of needed physical quantities and their spatial derivatives. The particles have assigned a characteristic distance, called ‘smoothing length’, over which their properties are ”smoothed” by a kernel function. A typical drawback of the SPH method is that it is hard to reproduce accurately the incompressibility of the materials. The SPH technique has been used successfully for solving fluid-structure interaction problems [3].

​

===1.2.2 Stabilization techniques===

​

A FEM-based procedure may require to be stabilized when incompressible, or quasi- incompressible, materials are analyzed. For example from the finite element solution of the Navier Stokes equations numerical instabilities may arise from two sources. The first one is due to the presence of the convective term in the linear momentum equations. This term introduces a non-linearity in the equations and it needs a proper stabilization for solving high Reynolds number flows with the FEM [38]. Furthermore, the orders of interpolation of pressure and velocity fields cannot be chosen freely but they have to satisfy the so called <math display="inline">{Ladyzenskaya-Babuzka-Brezzi}</math> (LBB), or <math display="inline">{inf-sup}</math>, condition [15]. If the orders of interpolation of the unknown fields do not satisfy this restriction, a stabilization technique is required in order to avoid numerical instabilities, as the spurious oscillations of the pressure field.

​

It is well known that the weak form generated by Galerkin approximation leads to a less diffusive solution than the strong problem. So the main idea of many stabilization techniques consists of adding an artificial diffusion to the problem. The first attempt was made by Von Neumann and Richtmyer [75]. Their solution adds an artificial diffusion to the strong form of the problem. This technique introduces an excessive dissipation to the problem because the diffusivity is added in every direction. An important evolution of this approach was the ''Streamline Upwind Petrov Galerkin (SUPG)'' [58]. In this approach, the artificial diffusion is added by means of the test functions and only on the direction of the streamlines. Furthermore, this is performed in a consistent way: the stabilization terms vanish when the solution is reached. An extension of the SUPG method is the ''Galerking Least-Squares (GLS)'' method [60]. The main difference between the two is that in the GLS method the stabilization terms are applied not only to the convective term but to all the terms of the equation. The ''Variational Multi-Scale (VMS)'' methods [59] split the problem variables in a large-scale and a subscale terms. The large-scale terms represent the part of the solution that can be captured by the finite element mesh, while the subscale part consists of an approximate solution that has to be added to the large-scale term in order to obtain the correct solution. This idea represents the basis of other stabilization methods. Among these, the most largely used are the ''Algebraic Subgrid Scale Formulation (ASGS)'' and the ''Orthogonal Subscales (OSS)'' methods, respectively introduced in [28,29]. Another efficient stabilization technique is the ''Finite Calculus '' (also termed ''Finite Increment Calculus'')  (''FIC'') approach  (see among others [76,77,85,86]). This method  has some analogies with the SUPG technique (for a comparison between these methods see [89]). The FIC approach is based on expressing the equations of balance of mass and momentum in a space-time domain of finite size and retaining higher order terms in the Taylor series expansion typically used for expressing the change in the transported variables within the balance domain. In addition to the standard terms of infinitesimal theory, the FIC form of the momentum and mass balance equations contains derivatives of the classical differential equations in mechanics multiplied by characteristic distances in space and time. In this work an updated version of the FIC method has been derived and tested.

​

In solid mechanics a stabilization procedure may be  required for the solution of problems involving  incompressible, or quasi-incompressible solids. Situations of this type are common in forming processes or in the analysis of the rubber-type materials. Many of the stabilization procedures  for fluid dynamics have been also used also for solid dynamics. For example, the VMS method has been applied in quasi-incompressible solid mechanics in [22,24], the OSS in [26]. In  [48,49] a stabilized multi-field Petrov-Galerkin procedure is used. Finally, the  application of the FIC in solid mechanics is reported [94,103].

​

===1.2.3 Algorithms for FSI problems===

​

Many approaches have been developed for solving FSI problems. Typically the computational techniques for  FSI are distinguished in ''monolithic'' and ''staggered'', or ''partitioned'', approaches. In a monolithic approach the fluid and the solid domains are solved in a single system of equations (see among others [57], or [73]). In this technique the flow of information between the solid and the fluid parts is implicitly performed by the procedure. On the contrary, in  staggered schemes the fluid and the solid dynamics are solved separately and boundary conditions are transferred from a domain to the other at the interface. From the algebraic point of view, the solution is achieved by solving two different linear systems which are coupled by means of the boundary conditions defined along the interface. Within partitioned schemes, depending on the level of coupling between the fluid and the solid dynamics, a further classification can be done. In a ''weakly coupled'' segregated schemes when the transfer of information through the interface is performed only once for each time step,  (for example see [42]), whereas  in ''strongly coupled'' schemes this operation is performed within a convergence loop  (as a  reference see [34]). Clearly, a weakly coupled scheme has a lower computational cost but it can be used only when the interaction between the solid and the fluid domains is not strong or complex. Otherwise the algorithm may not find the correct numerical solution of the problem.

​

Partitioned methods allow the reutilization of existing solvers. Furthermore, the solid and the fluid solvers can be updated independently. Staggered schemes lead to smaller and better conditioned linear systems than the ones obtained with monolithic approaches. On the other hand, monolithic strategies are generally more stable and fast than staggered schemes, and they lead to a more accurate solution of FSI problems, because a stronger coupling is ensured.

​

Another general classification of the FSI algorithms is based upon the treatment of meshes. In ''conforming mesh'' methods,  in order to allow the transfer of information, the fluid and the solid meshes must have in common the nodes along the interface. Consequently, if the position of the interface nodes changes in a domain, also the other domain must modify its grid in order to guarantee the conformity of the two meshes along the interface.  On the contrary, in ''non-conforming mesh'' methods the interface and the related conditions are treated as constraints imposed on the model equations so that the fluid and solid equations can be solved independently from each other with their respective grids  [21]. This represents an important advantage because, typically, the mesh used for the fluid has an average size lower than the one used for the solid and so it is not necessary to refine the solid finite element grid near to the interface. However, non-conforming mesh algorithms are more complex to implement and it is not trivial to guarantee their robustness.

​

In the so-called <math display="inline">{ImmersedBoundaryMethod}</math> (<math display="inline">{IBM}</math>), the fluid is solved using an Eulerian grid and the solids are immersed on top of this mesh [101,104]. The interaction is ensured by penalizing the Navier-Stokes equations with the momentum forcing sources of the immersed structures. An evolution of the IBM is the <math display="inline">{ ImmersedStructural}</math> <math display="inline">{PotentialMethod}</math> (<math display="inline">{ISPM}</math>) where the structure is modeled as a potential energy functional solved over a cloud of integration points that move within the fixed fluid mesh [47,104]. Also in the <math display="inline">{ImmersedFinite}</math> <math display="inline">{ElementMethod}</math> (<math display="inline">{IFEM}</math>) [54,129] the structure acts as a momentum forcing source for the fluid governing equations, but in this case a Lagrangian mesh is employed  for the solid domains.

​

==1.3 Numerical model==

​

The aim of this work is to derive a finite element formulation capable to solve, through a unique set of equations and unknown variables, the mechanics of a general continuum.  The term 'general continuum' refers to a domain that may include compressible and quasi-incompressible solids, fluids or both interacting together. For this reason, the formulation is termed <math display="inline">{Unified}</math>. The Unified formulation is based on a mixed Velocity-Pressure Stabilized procedure and it has been implemented in a sequential <math display="inline">C++</math> code.

​

===1.3.1 Reasons===

​

There are many reasons for undertaking the above objective.

​

The first advantage of the Unified formulation is that it allows for solving fluid and solid dynamics problems by implementing and using a single calculation code.

​

Furthermore, if solids and fluids are solved via the same scheme, it  is simpler to implement the solver for FSI problems because it is not required neither changing the variables, neither implementing the transfer of transmission conditions through the interface. With this formulation solids and fluids represent regions of the same continuum and they differ only by the specific values of the material parameters. As a consequence, the FSI solver requires a small computational effort and it can be implemented by introducing just a few  specific functions. This will be explained in detail in the section dedicated to FSI problems.

​

Additionally, with the Unified formulation the most natural approach for solving  FSI problem is the monolithic one. This brings in the further advantage that the coupling is ensured strongly and an iteration loop is not required, as for staggered procedures.

​

Finally, using the same set of unknowns for the fluid and the solid domains reduces the ill-conditioning of the FSI solver, because the solution system does not include variables of different units of measure.

​

===1.3.2 Essential features===

​

The Unified formulation is a compromise between a fluid and a solid formulation and it should be capable to satisfy the requirements of each problem mechanics and finite element approximations. In other words, in order to solve adequately fluids and  solids with the same formulation, the solution procedure must take into account the constraints that both models impose.

​

Concerning the choice of the unknown variables, these are generally selected depending on the constitutive law of the materials. For Newtonian fluids, the Cauchy stress tensor is related directly to the deformation rate tensor. This implies that the velocities are the most useful unknowns in fluid mechanics. Conversely, in solid mechanics the displacements  are generally preferred as variables, as the stresses are related to the deformations. However, for solids there also exist constitutive laws expressed in rate form, as for hypoelastic models.

​

​

For the analysis of incompressible, or quasi-incompressible, materials a mixed formulation is required in order to overcome the associated numerical instabilities, both in solids and 

fluids. In fluids, the most popular combination of unknown variables is the velocity-pressure 

scheme. In solid mechanics, there are several mixed approaches, as the displacement-stress [16], 

the displacement-strain [23], the displacement-deformation gradient [49], the  

displacement-deformation gradient-jacobian [48], the  displacement-pressure [26] and the 

velocity-pressure formulation [110]. The combination of unknowns is selected depending on the 

desired accuracy for the stress and strain fields and the associated computational cost. Among 

the mentioned mixed approaches, the displacement- pressure and the velocity-pressure formulations 

lead to the lowest computational cost, because the number of unknowns is smaller.  However the 

velocity-pressure formula- tion has the further advantage that is the canonical scheme for fluid 

mechanics. Hence, according to these considerations, the Unified formulation has been based on a 

mixed

velocity-pressure  scheme.

​

Another key decision in the design of the Unified formulation concerns the choice of the description framework. Solids are typically solved  using a (total or updated) Lagrangian description, while for fluids, due to the large deformations they undertake, the Eulerian framework is generally preferred. However, for free surface flows, as the ones treated in this work, a Lagrangian description has the important advantage that the free surface boundaries are detected automatically, because the particles of the fluid coincide with the nodes of the mesh. The price for this is that a remeshing procedure is required in order to avoid the excessive distortion of the mesh. According to these considerations, the Unified formulation has been implemented using an updated Lagrangian description.

​

In this work, the Particle Finite Element Method (PFEM) [88] has been used for solving the fluid domain only. The PFEM is a Langrangian numerical technique that treats the mesh nodes as moving points which can freely move and even separate  from the domain to which they are initially attached representing, for instance, the effect of water drops. The PFEM is based on a remeshing procedure that efficiently combines the Delaunay tessellation and the  Alpha Shape method [39].  In order to reduce the computational cost associated to the remeshing step, all the unknowns are stored in the nodes and linear shape functions are used for the finite element interpolation. Conversely, the solid domain keeps a fixed mesh during the whole analysis. The reason for this is that in non-linear solid mechanics it is required to preserve not only the nodal unknowns but also elemental information as the stress state of the previous time step. A remeshing implies the elimination of the previous discretization and the creation of new elements. In order to keep the elemental information in the remeshing step, a projection procedure from the elements of the previous mesh to its nodes, and from these to the elements of the new mesh is required. These operations may introduce an interpolation error  in the scheme that affects the solution of solid mechanics problems. In the analysis of Newtonian fluids all the information can be stored in the nodes, so the remeshing does not affect the accuracy of the scheme. For these reasons, in this work the PFEM is used only for the fluid, while solids are solved via the classical FEM.

​

In order to improve the conditioning of the solution system, a Gauss-Seidel  segregated scheme is used. This means that the problem is solved iteratively and separating the unkwown variables, the velocities and the pressure. The solution method consists of two steps. In the first one the linear momentum equations are solved for the velocity increments and considering the pressure of the previous non-linear iteration in the residual term. In the second step, the continuity equation is solved for the pressure in the updated configuration  using the velocities computed at the first step. This two-step algorithm leads to a smaller and better conditioned linear system  than using a monolithic approach. A crucial point of this scheme is the derivation of the tangent matrix of the linear momentum equations. The velocity increments are solved by condensing the pressure. According to this procedure, the tangent matrix for the linear momentum equations of the mixed velocity-pressure formulation does not differ from the one of a velocity formulation. Furthermore, the derivation of the tangent matrix for the velocity formulation is easier, as it does not involve the pressure.  For these reasons, the first step towards the unified scheme is the derivation of the Velocity formulation. The mixed Velocity-Pressure formulation will be derived by exploiting the linearization performed for the Velocity formulation. In [79] another procedure for the derivation of the exact tangent matrix for an Updated Lagrangian scheme is described.

​

The condensation of the pressure in the tangent matrix of the linear momentum equations may induce the ill-conditioning of the linear system for the analysis of quasi-incompressible materials. This ill-conditioning emanates from the volumetric counterpart of the tangent matrix that  is governed by  the bulk modulus, that typically has high values that may compromise the conditioning of the linear system. In this work, a thorough study of the conditioning drawbacks associated to this scheme has been performed and a useful and easy to implement technique to overcome these inconveniences has been tested and validated.

​

The numerical method proposed in this work is an improvement versus other numerical schemes for quasi-incompressible materials where an arbitrarily defined pseudo-bulk modulus was used in both the linear momentum and the mass conservation equations [110,109,108].

.

​

The method  can be also related to the so-called Augmented Lagrangian (AL) procedures for solving the Navier-Stokes equations for incompressible [10,11,41,44,124] and weakly compressible [14,123].

​

In order to deal with incompressible (or quasi-incompressible) materials a mixed formulation is required. In this work  both the velocities and the pressure are interpolated using linear shape functions. This combination does not fulfil the <math display="inline">inf-sup</math> condition [15] and the problem needs to be stabilized. The required  stabilization is ensured  using an updated version of the FIC technique [82] applied to the mixed Velocity-Pressure formulation. The FIC stabilization method has a small intrusivity because its terms only affect the continuity equation. In fact for solving the linear momentum equations the same scheme as for the (not stabilized) mixed Velocity-Pressure formulation is used. It will be shown that the PFEM-FIC stabilized procedure guarantees  a good accuracy for the mass conservation of the free surface flows. The stabilization procedure is derived for quasi-incompressible fluids, but it can be easily extended to quasi-incompressible hypoelastic solids.

​

The derivation of the Unified stabilized formulation is carried on trying to maintain as much as possible the generality of the scheme. For this reason the constitutive relations are introduced only as the last step. The Velocity (''V'') and the mixed Velocity-Pressure (''VP'') formulations are derived first for a general material and then particularized for specific constitutive laws. For quasi-incompressible materials only the mixed Velocity-Pressure Stabilized (''VPS'') formulation can be used while for compressible materials the Velocity and the standard (not stabilized) mixed Velocity-Pressure formulations are both suitable. In particular, for  compressible solids,  the hypoelastic model has been chosen.

​

The FSI problem is solved in a monolithic way, hence solids and fluids are solved by the same linear system. For the solid domain it is possible to choose which formulation to use. Depending on the problem, one may chose the V, the VP or the VPS element.

​

The Unified formulation for FSI problems can be easily coupled with the heat transfer problem in order to solve coupled thermal mechanical problem. The coupling is performed via a staggered scheme. The effect of the heat is taken into account by considering the material properties depending on the temperature and including in the strain tensor the deformations induced by the temperature. With the PFEM Unified formulation also phase change problems can be modeled.

​

The Unified stabilized VP formulation with thermal coupling has been used for the analysis of an industrial application. The study concerned the analysis of the damages caused by the dropping of a volume of corium at high temperature on the structure of the pressure vessel in a Nuclear Power Plant (NPP).

​

===1.3.3 Outline===

​

This text is split in the following chapters.

​

In the next chapter two velocity-based FEMs are derived for a general compressible material. In the first section the Velocity formulation is derived and the incremental solution scheme is given. Then, the standard mixed Velocity-Pressure formulation is obtained as an extension of the Velocity formulation. The constitutive laws are not introduced for any scheme. In the third section the formulations are adapted for hypoelastic-plastic compressible solids. The chapter ends with several validation examples for non-linear solid mechanics.

​

In Chapter [[#3 Unified stabilized formulation for quasi-incompressible materials|3]],  the FIC stabilization strategy is introduced in the mixed stabilized VP scheme. This scheme is adapted for   quasi-incompressible Newtonian fluids and hypoelastic solids, in this order. Then the free surface problem is studied in detail. First the PFEM is described highlighting the advantages and the disadvantages of the method. A useful technique for modeling the slip conditions in Lagrangian flows is also explained. Next the mass conservation feature and the conditioning of the scheme are analyzed in detail. Concerning the former point, it will be shown that the PFEM-FIC formulation guarantees a good preservation of the mass in free surface flows. Regarding  the latter point, a practical strategy for improving the conditioning of the linear system is explained and tested. At the end of the chapter several validation examples for quasi-incompressible Newtoninan fluids and hypoelastic solids are given.

​

The FSI solution strategy is described in Chapter [[#4 Unified formulation for FSI problems|4]]. The monolithic scheme for coupling the mechanics of fluids and solids is explained in detail. Finally several validation examples of FSI problems are given.

​

Chapter [[#5 Coupled thermal-mechanical formulation|5]]  is dedicated to couple the Unified formulation for FSI with the heat problem. In the first section the heat problem is introduced and discretized using the FEM. Then the coupling strategy is explained and validated with numerical examples. Next, the procedure for modeling  phase change problems is described and an explicative numerical example is presented.

​

Chapter [[#6 Industrial application: PFEM Analysis Model of NPP Severe Accident|6]]  is fully devoted to the industrial application of the Unified stabilized  and thermally coupled strategy. The damages of the pressure vessel structure caused by the fall of a volume corium at high temperature are studied using two simplified models where the solid structure melts due to the heat transfer from the corium.

​

In Chapter [[#7 Conclusions and future lines of research|7]] the innovative contributions of the work are summarized  and the lines of research opened by this thesis are presented.

​

==1.4 Publications==

​

The following publications have emanated from the work carried out in this doctoral thesis

​

<ol>

​

<li> E. Oñate and A. Franci and J.M. Carbonell, Lagrangian formulation for finite ele-ment analysis of quasi-incompressible fluids with reduced mass losses, ''International Journal for Numerical Methods in Fluids'', 74 (10), 699-731, 2014;</li>

<li> E. Oñate and A. Franci and J.M. Carbonell, A particle finite element method (PFEM) for coupled thermal analysis of quasi and fully incompressible flows and fluid-structure interaction problems, ''Numerical Simulations of Coupled Problems in Engineering''. S.R. Idelsohn (Ed.), 33, 129-164, 2014;</li>

<li>E. Oñate and A. Franci and J.M. Carbonell, A particle finite element method for analysis of industrial forming processes, ''Computational Mechanics'', DOI: 10.1007/s00466-014-1016-2, 2014; </li>

<li> A. Franci and E. Oñate and J.M. Carbonell, On the effect of the bulk tangent ma-trix in partitioned solution schemes for nearly incompressible fluids, ''International Journal for Numerical Methods in Engineering'', DOI: 10.1002/nme.4839, 2014;</li>

<li>A. Franci and E. Oñate and J.M. Carbonell, Unified formulation for solid and fluid mechanics and FSI problems, ''Computer Methods in Applied Mechanics and Engineering'', (in preparation).</li>

​

</ol>

​

=Chapter 2. Velocity-based formulations for compressible materials=

​

In this chapter two velocity-based finite element formulations for compressible materials are presented, namely the Velocity (V) and the mixed Velocity-Pressure (VP) formulations. For both schemes the linear momentum equations are solved iteratively for the velocity increments. The linearization of the governing equations is performed without specifying any constitutive law. The aim of this chapter is  to maintain  as much as possible the generality of the algorithms, leaving the formulations open to different material models. It will be shown that the only requirement demanded to the constitutive laws is that the rate of stress must be linearly related with the rate of deformation.

​

There are several reasons that justify the presentation of the Velocity formulation. First of all, the tangent matrix of the linear momentum equations for the Velocity formulation holds also for the mixed Velocity-Pressure formulation but the linearization procedure is easier because the pressure does not appear. Furthermore, the Velocity formulation is useful for making some interesting and didactic comparisons with the mixed formulations.

​

After deriving the solution scheme  for the Velocity formulation, the mixed Velocity-Pressure method is presented. The governing equations are the linear momentum and the linear pressure-deformation rate equations. The latter is called continuity, or mass balance, equation for its similarity with the incompressibility constraint of the Navier-Stokes problem. As for the velocity scheme, in the mixed formulation the constitutive law is not specified. In Section [[#2.3 Hypoelasticity|2.3]] the hypoelastic-plastic model is presented and this constitutive law is inserted in both the Velocity and the mixed Velocity-Pressure schemes. The incremental solution scheme is explained in detail for both formulations. At the end of the chapter some validation examples for hypoelastic-plastic solids in statics as in dynamics are given.

​

==2.1 Velocity formulation==

​

In this section, the velocity formulation for solving transient problems for a general continuum is derived.  The governing equations are the linear momentum equations and they are derived in the updated Lagrangian (UL) framework. This means that the governing equations are integrated over the unknown configuration <math display="inline">{\Omega }</math> (the so-called updated configuration). As a consequence, the space derivatives for the UL description are computed with respect to the spatial coordinates.

​

===2.1.1 From the local form to the spatial semi-discretization===

​

In this section the spatial semi-discretization of the linear momentum equations is derived.

​

For a general continuum, the local form of the linear momentum equations using the UL description reads

​

<span id="eq-1"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\rho (\hbox{X},t) \frac{ \partial v(\hbox{X},t)}{\partial t}-{\partial \sigma (\hbox{X},t) \over \partial \hbox{x}}-b(\hbox{X},t)=0   \quad \quad in  \Omega \times (0,T)   </math>

|}

| style="width: 5px;text-align: right;" | (1)

|}

​

where <math display="inline">\rho </math> is the density of the material, <math display="inline">v</math> are the velocity vector, <math display="inline">\sigma </math> is the Cauchy stress tensor and <math display="inline">b</math> is the body force vector. The variables within the brackets are the independent variables:  '''X''' are the Lagrangian  or material coordinates vector, '''x''' the Eulerian or spatial coordinates vector and <math display="inline">t</math> is the time. For simplicity, in what follows the independent variables are not specified. The spatial and material coordinates are related through the motion tensor <math display="inline">{\Phi }</math> as

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\hbox{x}=\Phi (\hbox{X},t),  \hbox{X}= \Phi ^{-1}(\hbox{x},t)  </math>

|}

| style="width: 5px;text-align: right;" | (2)

|}

​

The set of governing equations is completed by the following conditions at the Dirichlet (<math display="inline">\Gamma _v</math>) and Neumann (<math display="inline">\Gamma _t</math>) boundaries

​

<span id="eq-3"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>v_i -v_i^p =0 \qquad \hbox{on }\Gamma _v </math>

|}

| style="width: 5px;text-align: right;" | (3)

|}

​

<span id="eq-4"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\sigma _{ij}n_j -t_i^p =0 \qquad \hbox{on }\Gamma _t </math>

|}

| style="width: 5px;text-align: right;" | (4)

|}

​

where <math display="inline">v_i^p</math> and <math display="inline">t_i^p</math> are the prescribed velocities and the prescribed tractions, respectively, and <math display="inline">n</math> is the unit normal vector.

​

In the following summation of terms for repeated indices is assumed, unless otherwise specified.

​

The spaces for the trial and test functions are defined, respectively, as

​

<span id="eq-5"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>v_ i \in {U},  \quad \quad {U}=\{ v_ i |v_ i  \in C^0,  v_ i=v_i^p   on \Gamma _v \}   </math>

|}

| style="width: 5px;text-align: right;" | (5)

|}

​

<span id="eq-6"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>w_ i \in {U_0}, \quad \quad {U_0}=\{ w_ i | w_ i  \in C^0,  w_ i=0   on \Gamma _v \}   </math>

|}

| style="width: 5px;text-align: right;" | (6)

|}

​

Multiplying  Eqs.([[#eq-1|1]]) by the test functions and integrating over the updated configuration domain, the following global integral form is obtained

​

<span id="eq-7"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\int _\Omega w_i \left(\rho \dot{v}_i- {\partial \sigma _{ij} \over \partial x_j}-b_i\right)d\Omega =0 </math>

|}

| style="width: 5px;text-align: right;" | (7)

|}

​

where the symbol <math display="inline">\dot{(\cdot )}</math> represents the material time derivative.

​

Integrating by parts  the term involving <math display="inline">\sigma _{ij}</math> in Eq.([[#eq-7|7]]) and using the Neumann boundary conditions ([[#eq-4|4]]) yields the weak variational form of the momentum equations as

​

<span id="eq-8"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\int _\Omega w_i \rho \dot{v}_i d\Omega + \int _\Omega \frac{\partial w_i}{\partial x_j} \sigma _{ij} d\Omega - \int _\Omega w_i b_i d\Omega - \int _{\Gamma _t} w_i t_i^p d\Gamma =0  </math>

|}

| style="width: 5px;text-align: right;" | (8)

|}

​

Eq.([[#eq-8|8]]) is the standard form of the Principle of Virtual Power [9].

​

The spatial discretization is introduced using the classical FEM-Galerkin procedure  [136]. Hence both the trial and the test functions are interpolated in space in terms of their nodal values by means of the same shape functions <math display="inline">{N}</math>

​

<span id="eq-9"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>v_i = \sum \limits _{I=1}^n N_I(X) \bar v_{iI} \quad ,\quad w_i = \sum \limits _{I=1}^n N_I(X) \bar w_{iI}     </math>

|}

| style="width: 5px;text-align: right;" | (9)

|}

​

where, assuming the use of  simplicial elements, <math display="inline">n=3/4</math> for 2D/3D problems is the number of the nodes of the element, <math display="inline">\bar{(\cdot )}</math> denotes a nodal value, the capital subscript specifies the node and the lower case subscript represents the cartesian direction. In this work, linear shape functions have been used for <math display="inline">N_I</math>.

​

Since Eq.([[#eq-9|9]]) must hold for all the test functions in the interpolation space, introducing the spatial discretization ([[#eq-9|9]]) into  Eq.([[#eq-8|8]]), the spatial semi-discretized form of the momentum equations in the UL framework for the node <math display="inline">{I}</math> reads

​

<span id="eq-10"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\underbrace{  \int _\Omega N_I \rho d\Omega \dot{v_i}  }_{\displaystyle{{{f}^{dyn}_{Ii}}}}  +  \underbrace{  \int _\Omega \frac{\partial N_I}{\partial x_j} \sigma _{ij} d\Omega  }_{\displaystyle{{{f}^{int}_{Ii}}}}= \underbrace{  \int _\Omega N_I b_i d\Omega + \int _{\Gamma _t} N_I t_i^p d\Gamma  }_{\displaystyle{{{f}^{ext}_{Ii}}}}  </math>

|}

| style="width: 5px;text-align: right;" | (10)

|}

​

where <math display="inline">{f^{dyn}}</math>, <math display="inline">{f^{int}}</math> and <math display="inline">{f^{ext}}</math> are the dynamic, internal and external forces, respectively, expressed in the UL framework.

​

For convenience, the semi-discretized form of the momentum equations in the total Lagrangian (TL) framework is also presented here. This is written as [9]

​

<span id="eq-11"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\underbrace{  \int _{\Omega _0} N_I \rho _0 d\Omega \dot{v_i} }_{\displaystyle{{^{TL}{f}^{dyn}_{Ii}}}}+ \underbrace{  \int _{\Omega _0} \frac{\partial N_I}{\partial X_j} P_{ij} d\Omega }_{\displaystyle{{^{TL}{f}^{int}_{Ii}}}} = \underbrace{  \int _{\Omega _0} N_I b_i d\Omega + \int _{\Gamma _0} N_I t^p_{0i} d\Gamma  }_{\displaystyle{{^{TL}{f}^{ext}_{Ii}}}} </math>

|}

| style="width: 5px;text-align: right;" | (11)

|}

​

where <math display="inline">{P}</math> is the first Piola-Kirchhoff stress tensor, or the nominal stress tensor, and <math display="inline">{^{TL}f^{dyn}}</math>, <math display="inline">{^{TL} f^{int}}</math> and <math display="inline">{^{TL} f^{ext}}</math> are the dynamic, internal and external forces, respectively, expressed in the TL framework.  All the variables with vectors subscript <math display="inline">{(\cdot )}_0</math> refer to the last known configuration. Note that Eq.([[#eq-11|11]]) can be obtained from Eq.([[#eq-10|10]]) by pull back transformations on all its terms [9].

​

For the sake of clarity in the notation, the terms referred to the TL description are denoted with the left index <math display="inline">{^{TL}(\cdot )}</math>. Unless otherwise specified, the variables belong to the UL description.

​

===2.1.2 Time integration===

​

In this work, the kinematic variables  have been integrated in time using a second order scheme. In particular, the implicit Newmark's integration rule has been adopted. For the general case, this rule states that accelerations and displacements are computed as

​

<span id="eq-12"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>^{n+1}{\dot{v}}</math>

| style="text-align: right;" | <math>= \frac{1}{\gamma \Delta t} \left({^{n+1}v} -{^{n}v}\right)- \frac{1-\gamma }{\gamma }{^{n}\dot{v}}  </math>

|-

| style="text-align: center;" | <math> ^{n+1}u</math>

| style="text-align: right;" | <math>={^{n} u}+ {\Delta t}\frac{\gamma - \beta }{\gamma }{^{n}v} + \Delta t\frac{\beta }{\gamma } {^{n+1}v}+ {\Delta t^2}\frac{\gamma{-} 2\beta }{2\gamma }{^{n} \dot{v}} </math>

|}

| style="width: 5px;text-align: right;" | (12)

|}

​

Where <math display="inline">{\beta }</math> and <math display="inline">{\gamma }</math> are the so-called Newmark's parameters [9]. This time integration scheme is unconditionally stable if the following relation holds

​

<span id="eq-13"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\gamma \ge 2 \beta \ge \frac{1}{2} </math>

|}

| style="width: 5px;text-align: right;" | (13)

|}

​

In the present work, the Newmark's parameters chosen are <math display="inline">{\beta = \frac{1}{4}}</math> and <math display="inline">{\gamma = \frac{1}{2}}</math>.

​

Replacing the numerical values of the constants in Eq.([[#eq-12|12]]) yields

​

<span id="eq-14"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>^{n+1}{\dot{v}}= \frac{2}{ \Delta t} \left({^{n+1}v} - {^{n}v}\right)- {^{n}\dot{v}}  </math>

|}

| style="width: 5px;text-align: right;" | (14)

|}

​

<span id="eq-15"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>^{n+1}u={^{n}u} + \frac{\Delta t}{2}\left({^{n+1}v}+ {^{n}v}\right)  </math>

|}

| style="width: 5px;text-align: right;" | (15)

|}

​

===2.1.3 Linearization===

​

Although the problem is set out in the UL framework, the linearization for the velocities increment of the momentum equations is performed first on the TL  semi-discretized form ([[#eq-11|11]]). The UL linearized form is  obtained by performing a push-forward transformation on the TL form. This is justified by the easier derivation of the tangent matrix in the TL framework. In fact, in Eq.([[#eq-11|11]]) the only variable that depends on time is the nominal stress <math display="inline">{P}</math>, while in the UL form ([[#eq-10|10]]) the time-dependent variables are the updated domain <math display="inline">{\Omega }</math>, the Cauchy stress tensor <math display="inline">{\sigma }</math> and the spatial derivatives <math display="inline">{   \partial N/ \partial x}</math>. For the sake of clarity, the linearization of the internal and dynamic forces will be performed separately.  

​

''Internal component of the tangent matrix''

​

From  Eq.([[#eq-11|11]]) the internal forces in the TL description are defined as

​

<span id="eq-16"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>^{TL}f^{int}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} P_{ij} d\Omega _0  </math>

|}

| style="width: 5px;text-align: right;" | (16)

|}

​

The constitutive relation is expressed in rate form. Hence it is more convenient to perform the linearization of the material derivative of the internal forces and then integrate for the time step increment <math display="inline">\Delta t</math>. The material time derivative of ([[#eq-16|16]]) is

​

<span id="eq-17"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>^{TL}\dot{f}^{int}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}\dot{P}_{ij} d\Omega _0  </math>

|}

| style="width: 5px;text-align: right;" | (17)

|}

​

The infinitesimal increment of Eq.([[#eq-17|17]]) is

​

<span id="eq-18"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>^{TL}\delta \dot{f}^{int}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} \delta  \dot{P}_{ij} d\Omega _0  </math>

|}

| style="width: 5px;text-align: right;" | (18)

|}

​

The first Piola-Kirchhoff stress tensor <math display="inline">{P}</math> is not typically used because it is not symmetric and its rate is a non-objective measure. For these reasons, in the TL framework it is more convenient to work with the second Piola-Kirchhoff stress tensor <math display="inline">{S}</math> and its rate. These stress rate measures are related each other through the following relation

​

<span id="eq-19"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\dot{P}_{ij}=\dot{S}_{ir} F^T_{rj} + S_{ir} \dot{F}^T_{rj}  </math>

|}

| style="width: 5px;text-align: right;" | (19)

|}

​

where <math display="inline">{F}</math> is the  deformation gradient tensor defined as

​

<span id="eq-20"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{F}_{ij}= \frac{\partial x_i}{\partial X_j}  </math>

|}

| style="width: 5px;text-align: right;" | (20)

|}

​

Substituting Eq.([[#eq-19|19]]) into ([[#eq-18|18]]) yields

​

<span id="eq-21"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>^{TL}\delta \dot{f}^{int}_{Ii}= \underbrace{\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} F_{ir} \delta \dot{S}_{jr}d\Omega _0} _{\displaystyle{{^{TL}\delta \dot{f}^{m}_{Ii}}}} +  \underbrace{\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}  S_{ir}  \delta \dot{F}^T_{rj} d\Omega _0} _{\displaystyle{{^{TL}\delta \dot{f}^{g}_{Ii}}}}     </math>

|}

| style="width: 5px;text-align: right;" | (21)

|}

​

Eq.([[#eq-21|21]]) shows that the increment of the material time derivative of the  internal forces can be split into material and   geometric parts, <math display="inline">{^{TL} \delta \dot{f}^{m}}</math> and <math display="inline">{^{TL}\delta \dot{f}^{g}}</math>, respectively. The former accounts for the material response through the rate of the second Piola-Kirchhoff stress tensor and the second term is the initial stress term that contains the information of the updated stress field. Note that, up to now, no constitutive relationships have been introduced and the above derivation holds for a general continuum. 

​

''Material tangent  matrix''

​

From Eq.([[#eq-21|21]]), the material part of the material time derivative of the internal forces reads

​

<span id="eq-22"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}\delta \dot{f}^{m}_{Ii}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} F_{ir} \delta \dot{S}_{jr}d\Omega _0}   </math>

|}

| style="width: 5px;text-align: right;" | (22)

|}

​

For the derivation of the material tangent matrix, the constitutive relation between the stress and the strain measures is required. In order to maintain the formulation as general as possible, the stress rate measure is related to the deformation rate through a tangent constitutive tensor, as

​

<span id="eq-23"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\dot{S}_{ij}=C_{ijkl} \dot{E}_{kl}  </math>

|}

| style="width: 5px;text-align: right;" | (23)

|}

​

where <math display="inline">{C}</math> is a fourth-order  tensor and <math display="inline">{E}</math> is the Green-Lagrange strain tensor. Eq.([[#eq-23|23]]) can also be expressed in Voigt notation as

​

<span id="eq-24"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\{ \dot{S} \}=\left[C\right]\{ \dot{E}\}  </math>

|}

| style="width: 5px;text-align: right;" | (24)

|}

​

where <math display="inline">{(\cdot )}</math> denotes a vector with components <math display="inline">\left[(\cdot )_{11}, (\cdot )_{22}, (\cdot )_{33}, (\cdot )_{12}, (\cdot )_{13},  (\cdot )_{23} \right]</math>.  As it will be shown in the following sections, Eq.([[#eq-23|23]]) can represent  both a Kirchhoff solid material and a Newtonian fluid.  If different constitutive laws are used, Eq.([[#eq-23|23]]) should be modified accordingly in order to derive the material part of the tangent matrix.

​

The Green-Lagrange strain tensor can be expressed in terms of the nodal velocities as

​

<span id="eq-25"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\dot{E}_{kl}= \frac{\partial N_J}{\partial X_l}  F_{sk}  \bar  v_{Js}  </math>

|}

| style="width: 5px;text-align: right;" | (25)

|}

​

In Voigt notation

​

<span id="eq-26"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\{ \dot{E} \}= B_0  \bar{v}  </math>

|}

| style="width: 5px;text-align: right;" | (26)

|}

​

where for a plane strain problem

​

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: right;" | <math>\displaystyle{B_{0}}= \left[                    \begin{matrix}\displaystyle                   \frac{\partial N_I}{\partial {X}}\frac{\partial x}{\partial {X}} & \displaystyle                    \frac{\partial N_I}{\partial {X}}\frac{\partial y}{\partial {X}} \\[.4cm] \displaystyle   \frac{\partial N_I}{\partial {Y}}\frac{\partial x}{\partial {Y}} & \displaystyle                    \frac{\partial N_I}{\partial {Y}}\frac{\partial y}{\partial {Y}} \\[.4cm] \displaystyle    \frac{\partial N_I}{\partial {X}}\frac{\partial x}{\partial {Y}}+\frac{\partial N_I}{\partial {Y}}\frac{\partial x}{\partial {X}} &   \displaystyle   \frac{\partial N_I}{\partial {X}}\frac{\partial y}{\partial {Y}}+\frac{\partial N_I}{\partial {Y}}\frac{\partial y}{\partial {X}} \\                      \end{matrix} \right] </math>

|}

| style="width: 5px;text-align: right;" | (30)

|}

​

Substituting Eq.([[#eq-25|25]]) in  Eq.([[#eq-23|23]]), yields

​

<span id="eq-31"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\dot{S}_{ij}=C_{ijkl}  \frac{\partial N_J}{\partial X_l}  F_{sk} \bar   v_{Js}  </math>

|}

| style="width: 5px;text-align: right;" | (31)

|}

​

Similarly, substituting  Eq.([[#eq-26|26]]) in  Eq.([[#eq-24|24]])

​

<span id="eq-32"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\{ \dot{S} \}=\left[C\right] B_0 v  </math>

|}

| style="width: 5px;text-align: right;" | (32)

|}

​

Substituting Eq.([[#eq-31|31]]) into Eq.([[#eq-22|22]]) yields

​

<span id="eq-33"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}\delta \dot{f}^{m}_{Ii}}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}  F_{ir} C_{jrkl}  \frac{\partial N_J}{\partial X_l}  F_{sk}  d\Omega _0  \delta \bar   v_{Js}  </math>

|}

| style="width: 5px;text-align: right;" | (33)

|}

​

where <math display="inline">{ \bar{v}_{Js}}</math> is the <math display="inline">{s}</math>-component of the velocity of node <math display="inline">{J}</math>. In Voigt notation Eq.([[#eq-33|33]]) reads

​

<span id="eq-34"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}\delta \dot{f}^{m}}=\int _{\Omega _0} B_0^T \left[C \right]B_0 d\Omega _0 \delta \bar{v}   </math>

|}

| style="width: 5px;text-align: right;" | (34)

|}

​

In order to obtain the increment of the internal forces, the material time derivative of  the internal forces increment is integrated over a time step increment <math display="inline">\Delta t</math> as

​

<span id="eq-35"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}\delta {f}^{m}}= {^{TL}\delta \dot{f}^{m}} \Delta t  </math>

|}

| style="width: 5px;text-align: right;" | (35)

|}

​

From Eq.([[#eq-35|35]]) and Eq.([[#eq-33|33]]), yields

​

<span id="eq-36"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}\delta{f}^{m}_{Ii}}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}  F_{ir} \Delta t C_{jrkl}  \frac{\partial N_J}{\partial X_l}  F_{sk}  d\Omega _0  \delta \bar   v_{Js}  </math>

|}

| style="width: 5px;text-align: right;" | (36)

|}

​

and, similarly,  from Eq.([[#eq-35|35]]) and Eq.([[#eq-34|34]])

​

<span id="eq-37"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}\delta {f}^{m}}=\int _{\Omega _0} B_0^T \Delta t\left[C \right]B_0 d\Omega _0 \delta \bar{v}   </math>

|}

| style="width: 5px;text-align: right;" | (37)

|}

​

From Eq.([[#eq-36|36]]) and Eq.([[#eq-37|37]]), the material tangent matrix in the TL description can be computed as

​

<span id="eq-38"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>^{TL}{K}^{m}_{IJis}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}  F_{rj} \Delta t C_{ijkl}  \frac{\partial N_J}{\partial X_s}  F_{kl}  d\Omega  </math>

|}

| style="width: 5px;text-align: right;" | (38)

|}

​

<span id="eq-39"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}{K}^{m}}=\int _{\Omega _0} B_0^T \Delta t\left[C \right]B_0 d\Omega _0   </math>

|}

| style="width: 5px;text-align: right;" | (39)

|}

​

The material tangent matrix for the UL framework is obtained by applying a push-forward transformation on each term of Eq.([[#eq-38|38]]) and integrating over the updated domain <math display="inline">\Omega </math>. The following relations hold

​

<span id="eq-40"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>d{\Omega _0}=\frac{ d \Omega }{J}   </math>

|}

| style="width: 5px;text-align: right;" | (40)

|}

​

<span id="eq-41"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\frac{\partial N_I}{\partial X_j}=\frac{\partial N_I}{\partial x_k} F_{kj}  </math>

|}

| style="width: 5px;text-align: right;" | (41)

|}

​

<span id="eq-42"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>C_{ijkl}=F^{-1}_{mi} F^{-1}_{nj} F^{-1}_{ok} F^{-1}_{pl} c^{\tau }_{mnop} = F^{-1}_{mi} F^{-1}_{nj} F^{-1}_{ok} F^{-1}_{pl} c^{\sigma }_{mnop} J  </math>

|}

| style="width: 5px;text-align: right;" | (42)

|}

​

where <math display="inline">{c^{\tau }}</math> is the tangent moduli corresponding to the material time derivative of the Kirchhoff stress tensor <math display="inline">{\tau ^{\bigtriangledown }}</math> and <math display="inline">{c^{ \sigma }}</math> is the tangent moduli for the rate of the Cauchy stress <math display="inline">{\sigma ^{\bigtriangledown }}</math>. The rate of the Cauchy stress tensor is related to the rate of deformation <math display="inline">{d}</math> through the fourth-order tensor <math display="inline">{c^{\sigma }}</math> by the following expression

​

<span id="eq-43"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\sigma ^{\bigtriangledown }=c^{\sigma } : d  </math>

|}

| style="width: 5px;text-align: right;" | (43)

|}

​

Substituting Eqs.([[#eq-40|40]]-[[#eq-42|42]]) into ([[#eq-38|38]]) and using the minor symmetries,  the material tangent matrix  for the UL reads

​

<span id="eq-44"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>K^{m}_{IJrs}=\int _{\Omega } \frac{\partial N_I}{\partial x_k} \delta _{ri} \Delta t c^{\sigma }_{kijl}  \frac{\partial N_J}{\partial x_l} \delta _{sj}    d\Omega  </math>

|}

| style="width: 5px;text-align: right;" | (44)

|}

​

Using Voigt notation, the same matrix reads

​

<span id="eq-45"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>K^{m}_{IJ}=\int _{\Omega ^e} B^{T}_I  \Delta t \left[c^{\sigma } \right] B_J    d\Omega   </math>

|}

| style="width: 5px;text-align: right;" | (45)

|}

​

For the node <math display="inline">I</math> of a 3D element,  matrix <math display="inline">B</math> is

​

<span id="eq-51"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: right;" | <math>\displaystyle{B}_I = \left[\begin{matrix}\displaystyle {\partial N_I \over \partial x} &0&0\\ \displaystyle{0}& \displaystyle {\partial N_I \over \partial y}&0\\ \displaystyle{0}&0&\displaystyle {\partial N_I \over \partial z}\\ \displaystyle {\partial N_I \over \partial y}&\displaystyle {\partial N_I \over \partial x}&0\\[.25cm] \displaystyle {\partial N_I \over \partial z}&0&\displaystyle {\partial N_I \over \partial x}\\[.25cm] \displaystyle{0}&\displaystyle {\partial N_I \over \partial z}&\displaystyle {\partial N_I \over \partial y}          \end{matrix}  \right]  </math>

|}

| style="width: 5px;text-align: right;" | (51)

|}

​

​

The geometric tangent matrix for the UL framework is  derived next using the same procedure adopted for the material components. Hence, first the linearization is performed using the TL form and then the UL tangent matrix is obtained by performing the required transformation over the TL terms.

​

From Eq.([[#eq-21|21]])

​

<span id="eq-52"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}\delta \dot{f}^{g}_{Ii}}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j}  S_{ir}  \delta \dot{F}^T_{rj} d\Omega _0  </math>

|}

| style="width: 5px;text-align: right;" | (52)

|}

​

where the rate of the deformation gradient is defined as

​

<span id="eq-53"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>\dot{F}_{ij}= \frac{\partial N_J}{\partial X_i}  \bar v_{Jj}  </math>

|}

| style="width: 5px;text-align: right;" | (53)

|}

​

Substituting Eq.([[#eq-53|53]]) into Eq.([[#eq-52|52]]), the geometric components of the internal power in the TL description can be written as

​

<span id="eq-54"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}\delta \dot{f}^{g}_{Ii}}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} S_{ir} \frac{\partial N_J}{\partial X_r}   d\Omega _0  \delta \bar  v_{Jj}  </math>

|}

| style="width: 5px;text-align: right;" | (54)

|}

​

Integrating Eq.([[#eq-54|54]]) on time for a time step increment <math display="inline">\Delta t</math> yields

​

<span id="eq-55"></span>

{| class="formulaSCP" style="width: 100%; text-align: left;" 

|-

| 

{| style="text-align: left; margin:auto;" 

|-

| style="text-align: center;" | <math>{^{TL}\delta {f}^{g}_{Ii}}=\int _{\Omega _0} \frac{\partial N_I}{\partial X_j} \Delta t S_{ir} \frac{\partial N_J}{\partial X_r}   d\Omega  \delta \bar  v_{Jj}  </math>

|}

| style="width: 5px;text-align: right;" | (55)

|}

​