​

​

==Abstract==

​

We present a Lagrangian formulation for coupled thermal analysis of quasi and fully incompressible flows and fluid-structure interaction (FSI) problems that has excellent mass preservation features. The success of the formulation lays on  a  residual-based stabilized expression of the mass balance equation obtained using the Finite Calculus (FIC) method. The  governing equations are discretized with the FEM using simplicial elements with equal linear interpolation for the velocities, the pressure and the temperature. The merits of the formulation in terms of reduced mass loss and overall accuracy are verified in the solution of 2D and 3D adiabatic and thermally-coupled quasi-incompressible free-surface flow problems using the Particle Finite Element Method (PFEM). Examples include the sloshing of water in a tank and the falling  of a water sphere and a cylinder into a  tank containing water.

​

'''keywords''' Particle Finite Element Method, Coupled Thermal Analysis, Quasi and Fully Incompressible

​

==1 INTRODUCTION==

​

The analysis of thermally coupled flows and their interaction with structures is relevant in many fields of engineering. In this work we present a Lagrangian numerical technique for solving this kind of problems for quasi and fully incompressible fluids using the Particle Finite Element Method (PFEM, www.cimne.com/pfem).

​

The PFEM treats the mesh nodes in the analysis domain as particles which can freely move and even separate from the  domain representing, for instance, the effect of water drops or cutting particles in drilling problems. A mesh connects the nodes discretizing the domain where the governing equations are solved using a stabilized FEM. Examples of application of PFEM to problems in fluid and solid mechanics including fluid-structure interaction (FSI) situations can be found in <span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-28'></span><span id='citeF-29'></span><span id='citeF-35'></span><span id='citeF-36'></span><span id='citeF-39'></span><span id='citeF-40'></span><span id='citeF-41'></span><span id='citeF-42'></span><span id='citeF-43'></span>[[#cite-4|[4,5,6,8,9,10,11,12,13,14,15,16,17,18,19,28,29,35,36,39,40,41,42,43]]]. Early attempts of the PFEM for solving thermally coupled flows were reported in <span id='citeF-1'></span><span id='citeF-2'></span>[[#cite-1|[1,2]]].

​

In Lagrangian analysis procedures (such as  PFEM) the motion of the fluid particles is tracked during the transient solution. Hence, the convective terms vanish in the momentum and heat transfer equations and no numerical stabilization is needed for treating those terms. Two other sources of mass loss, however, remain in the numerical solution of Lagrangian flows, i.e. that due to the treatment of the incompressibility constraint by a stabilized numerical method, and that induced by the inaccuracies in tracking the flow particles and, in particular, the free surface.

​

In this work the PFEM equations for analysis of thermally coupled flows and FSI problems are derived using the stabilized formulation based in the Finite Calculus (FIC) method proposed by Oñate et al. <span id='citeF-20'></span><span id='citeF-21'></span><span id='citeF-22'></span><span id='citeF-23'></span><span id='citeF-24'></span><span id='citeF-25'></span><span id='citeF-26'></span><span id='citeF-27'></span><span id='citeF-30'></span><span id='citeF-31'></span><span id='citeF-32'></span><span id='citeF-37'></span><span id='citeF-38'></span><span id='citeF-39'></span>[[#cite-20|[20,21,22,23,24,25,26,27,30,31,32,37,38,39]]] that has excellent mass preservation features.

​

The lay-out of the paper is the following. In the next section we present the basic equations for conservation of linear momentum, mass and heat transfer for a quasi-incompressible fluid in a Lagrangian framework. A full incompressible fluid can be considered as a particular limit case of the former. Next we derive the stabilized FIC form of the mass balance equation. Then the finite element discretization using simplicial element with equal order approximation for the velocity,  the pressure and the temperature is presented and the relevant matrices and vectors of the discretized problem are given.  Details of the implicit solution of the Lagrangian FEM equations in time using a Newton iterative scheme are presented. The relevance of the bulk stiffness terms in the tangent matrix  for enhancing the convergence and accuracy of the iterative solution scheme is discussed. The basic steps of the PFEM for solving coupled free-surface FSI problems are described.

​

The efficiency and accuracy of the PFEM  technique are verified by solving a set of adiabatic and thermally coupled quasi-incompressible free surface flow problems in two (2D) and three (3D) dimensions with the PFEM. The adiabatic problems are the sloshing of water in a tank and the penetration of a water sphere into a cylindrical tank containing water. The thermally coupled problems considered are the extended 2D version of the adiabatic cases. The excellent performance of the numerical method proposed in terms of mass conservation and general accuracy is highlighted.

​

==2 GOVERNING EQUATIONS==

​

We write the governing equations for a quasi-incompressible Newtonian flow problem in the Lagrangian description as follows <span id='citeF-3'></span><span id='citeF-46'></span>[[#cite-3|[3,46]]].

​

===Momentum equations===

​

<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 \frac{Dv_i}{Dt}-{\partial \sigma _{ij} \over \partial x_j}-b_i=0\quad , \quad i,j=1,n_s \quad  \hbox{in } \Omega  </math>

|}

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

|}

​

In Eq.([[#eq-1|1]]), <math display="inline">\Omega </math> is the analysis domain with boundary <math display="inline">\Gamma </math>, <math display="inline">v_i</math> and <math display="inline">b_i</math> are the velocity and body force components  along the <math display="inline">i</math>th Cartesian axis, <math display="inline">\rho </math> is the density, <math display="inline">n_s</math> is the number of space dimensions (i.e. <math display="inline">n_s=3</math> for 3D problems) and <math display="inline">\sigma _{ij}</math> are the Cauchy stresses that are split in the deviatoric (<math display="inline">s_{ij}</math>) and pressure (<math display="inline">p</math>) components as

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\sigma _{ij} = s_{ij}+ p \delta _{ij} </math>

|}

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

|}

​

where <math display="inline">\delta _{ij}</math> is the Kronecker delta. Note that the pressure is assumed to be positive for a tension state. Summation of terms with repeated indices  is assumed in Eq.([[#eq-1|1]]) and in  the following, unless otherwise specified.

​

The relationship between the deviatoric stresses <math display="inline">s_{ij}</math> and the strain rates <math display="inline">\varepsilon _{ij}</math> has the standard form for a Newtonian fluid,

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>s_{ij} = 2\mu \left(\varepsilon _{ij} - {1 \over 3}\varepsilon _v \delta _{ij}\right)\quad \hbox{with}\quad  \varepsilon _{ij} = {1 \over 2}\left({\partial v_i \over \partial x_j} +  {\partial v_j \over \partial x_i}  \right) </math>

|}

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

|}

​

where <math display="inline">\mu </math> is the viscosity  and <math display="inline">\varepsilon _v</math> is the volumetric strain rate defined as <math display="inline">\varepsilon _v= \varepsilon _{ii} = {\partial v_i \over \partial x_i} </math>.

​

===Mass balance equation===

​

The standard mass balance equation for a quasi-incompressible fluid can be written as <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]]

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>r_v =0 </math>

|}

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

|}

​

with

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>r_v:=-\frac{1}{c^2}\frac{Dp}{Dt}+\rho \varepsilon _v </math>

|}

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

|}

​

In Eq.([[#eq-5|4b]]) <math display="inline">c</math> is the speed of sound in the fluid. For a fully incompressible fluid <math display="inline">c=\infty </math> and Eq.([[#eq-4|4a]]) simplifies to the standard form, <math display="inline">\varepsilon _v =0</math>. In our work we will retain the quasi-incompressible form of <math display="inline">r_v</math> of Eq.([[#eq-5|4b]]) for convenience.

​

===Thermal balance===

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\rho c \frac{DT}{Dt} - {\partial  \over \partial x_i} \left(k {\partial T \over \partial x_i}\right)+ Q =0  \quad i=1,n_s  \quad \hbox{in }\Omega   </math>

|}

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

|}

​

where <math display="inline">T</math> is the temperature, <math display="inline">c</math> is the thermal capacity, <math display="inline">k</math> is the heat conductivity and <math display="inline">Q</math> is the heat source.

​

===Boundary conditions===

​

===''Mechanical problem''===

​

The boundary conditions at the Dirichlet (<math display="inline">\Gamma _v</math>) and Neumann (<math display="inline">\Gamma _t</math>) boundaries with <math display="inline">\Gamma  = \Gamma _v \cup \Gamma _t</math> are

​

<span id="eq-6"></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;" | (6)

|}

​

<span id="eq-7"></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 \quad \hbox{on }\Gamma _t \quad i,j=1,n_s  </math>

|}

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

|}

​

where <math display="inline">v_i^p</math> and <math display="inline">t_i^p</math> are the prescribed velocities and prescribed tractions on <math display="inline">\Gamma _v</math> and <math display="inline">\Gamma _t</math>, respectively and <math display="inline">n_j</math> are the components of the unit normal vector to the boundary <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]].

​

===''Thermal problem''===

​

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

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\phi - \phi ^p =0 \quad \hbox{on }\Gamma _\phi </math>

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

|-

| style="text-align: center;" | <math> k {\partial \phi  \over \partial n} + q_n^p =0 \quad \hbox{on }\Gamma _q  </math>

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

|}

|}

​

where <math display="inline">\phi ^p</math> and <math display="inline">q_n^p</math> are the prescribed temperature and the prescribed normal heat flux at the boundaries <math display="inline">\Gamma _\phi </math> and <math display="inline">\Gamma _q</math>, respectively and <math display="inline">n</math> is the direction normal to the boundary.

​

'''Remark 1.''' The term <math display="inline">\frac{Dv_i}{Dt}</math> in Eq.([[#eq-1|1]]) is the ''material derivative'' of the <math display="inline">i</math>th velocity component <math display="inline">v_i</math>. This term is typically computed in a Lagrangian framework as

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\frac{Dv_i}{Dt} = \frac{{}^{n+1} v_i - {}^n v_i }{\Delta t}     </math>

|}

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

|}

​

with

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>{}^{n+1} v_i:=v_i ({}^{n+1} {\bf x},{}^{n+1} t)\quad ,\quad {}^n v_i:=v_i ({}^n{\bf x},{}^n t)       </math>

|}

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

|}

​

where <math display="inline">{}^n v_i</math> is the velocity of the material point that has the position <math display="inline">{}^n{\bf x}</math> at time <math display="inline">t={}^nt</math>, where <math display="inline">{\bf x}</math> is the coordinates vector in a fixed Cartesian  system <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]].

​

​

==3 STABILIZED MASS BALANCE EQUATION==

​

In this work we will use  the ''second order FIC form'' of the mass balance equation in space for a quasi-incompressible fluid <span id='citeF-37'></span><span id='citeF-38'></span>[[#cite-37|[37,38]]], as well as the first order FIC form of the mass balance equation in time. These forms have the following expressions:

​

===''Second order FIC mass balance equation in space''===

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>r_v + \frac{h_i^2}{12} \frac{\partial ^2 r_v}{\partial x^2_i}=0\qquad \hbox{in }\Omega \qquad i=1,n_s </math>

|}

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

|}

​

===''First order FIC mass balance equation in time''===

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>r_v + \frac{\delta }{2} \frac{D r_v}{D t}=0 \qquad \hbox{in }\Omega  </math>

|}

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

|}

​

Eq.([[#eq-13|12a]]) is obtained by expressing the balance of mass in a rectangular domain of finite size with dimensions <math display="inline">h_1\times h_2</math> (for 2D problems), where <math display="inline">h_i</math> are arbitrary distances, and retaining up to third order terms in the Taylor series expansions used for expressing the change of mass within the balance domain.

​

Eq.([[#eq-14|12b]]), on the other hand, is obtained by expressing the balance of mass in a space-time domain of infinitesimal length in space and finite dimension <math display="inline">\delta </math> in time <span id='citeF-20'></span>[[#cite-20|[20]]]. The derivation of Eqs.([[#''First order FIC mass balance equation in time''|12]]) for a 1D problem are shown in <span id='citeF-41'></span>[[#cite-41|[41]]].

​

The FIC terms in Eqs.([[#''First order FIC mass balance equation in time''|12]]) play the role of space and time stabilization terms respectively. In the discretized problem, the space dimensions <math display="inline">h_i</math> and the time dimension <math display="inline">\delta </math> are related to characteristic element dimensions and the time step increment, respectively as it will be explained later. Note that for <math display="inline">h_i =0</math> and <math display="inline">\delta =0</math> the standard infinitesimal form of the mass balance equation, <math display="inline">r_v =0</math>, is obtained.

​

After some transformations the stabilized mass balance equation ([[#eq-13|12a]]) is written as <span id='citeF-41'></span>[[#cite-41|[41]]]

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>-\frac{1}{\kappa }\frac{Dp}{Dt}+\varepsilon _v  - \frac{\tau }{c^2} \frac{D^2p}{Dt^2} +\tau {\partial \hat r_{m_i} \over \partial x_i} =0  </math>

|}

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

|}

​

where <math display="inline">\tau </math> is a ''stabilization parameter'' given by <span id='citeF-41'></span>[[#cite-41|[41]]]

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\tau = \left(\frac{8\mu }{h^2}+ \frac{2\rho }{\delta } \right)^{-1}  </math>

|}

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

|}

​

and <math display="inline">\hat r_{m_i}</math> is a ''static momentum term'' defined as

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\hat r_{m_i}= {\partial  \over \partial x_j} (2\mu \varepsilon _{ij}) + \frac{\partial p}{\partial x_i} + b_i   </math>

|}

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

|}

​

Eq.([[#eq-13|13]]) is used as the starting point for deriving the stabilized FEM formulation as explained in the following sections.

​

==4 VARIATIONAL EQUATIONS==

​

===4.1 Momentum equations===

​

Multiplying Eq.(1) by arbitrary test functions <math display="inline">w_i</math> with dimensions of velocity and integrating over the analysis domain <math display="inline">\Omega </math> gives the weighted residual form of the momentum equations as <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]]

​

<span id="eq-16"></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 \frac{Dv_i}{Dt}- {\partial \sigma _{ij} \over \partial x_j}-b_i\right)d\Omega =0 </math>

|}

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

|}

​

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

​

<span id="eq-17"></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 \frac{Dv_i}{Dt} d\Omega + \int _\Omega \delta \varepsilon _{ij} \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;" | (17)

|}

​

where <math display="inline">\delta \varepsilon _{ij}= {\partial w_i \over \partial x_j}+{\partial w_j \over \partial x_i}</math> is an arbitrary (virtual) strain rate field. Eq.([[#eq-17|17]]) is the standard form of the ''principle of virtual power''  <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-46'></span>[[#cite-3|[3,7,46]]].

​

Substituting the expression of the stresses from Eq.([[#eq-2|2]]) into ([[#eq-17|17]]) gives

​

<span id="eq-18"></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 \frac{Dv_i}{Dt} d\Omega + \!\int _\Omega \left[\delta \varepsilon _{ij} 2\mu \left(\!\varepsilon _{ij} - \frac{1}{3}\varepsilon _{v} \delta _{ij}\!\right)+ \delta \varepsilon _{v}p\right]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;" | (18)

|}

​

Eq.([[#eq-18|18]]) can be written in matrix form as

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\displaystyle  \int _\Omega {\bf w}^T \rho \frac{D{\bf v}}{Dt} d\Omega + \!\int _\Omega \delta {\boldsymbol \varepsilon }^T {\bf D} {\boldsymbol \varepsilon } d\Omega  + \!\int _\Omega \delta {\boldsymbol \varepsilon }^T {\bf m} p d\Omega - \int _\Omega {\bf w}^T {\bf b}d\Omega - \!\int _{\Gamma _t} {\bf w}^T {\bf t}^p d\Gamma =0  </math>

|}

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

|}

​

In Eq.([[#eq-19|19]]) <math display="inline">{\bf w},{\bf v},{\boldsymbol \varepsilon }</math> and <math display="inline">\delta {\boldsymbol \varepsilon }</math> are vectors containing the test functions, the velocities, the strain rates and the virtual strain rates respectively; <math display="inline">{\bf b}</math> and <math display="inline">{\bf t}^p</math> are body force and surface traction vectors, respectively; <math display="inline">{\bf D}</math> is the viscous constitutive matrix and <math display="inline">{\bf m}</math> is an auxiliary vector. These vectors are defined as (for 3D problems)

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{l}\textbf{w} =[w_1,w_2,w_3]^T\quad ,\quad \textbf{v} =[v_1,v_2,v_3]^T \quad ,\quad \textbf{b} =[b_1,b_2,b_3]^T \quad ,\quad  \textbf{t}^p =[t_1^p,t_2^p,t_3^p]^T\\[.5cm] {\boldsymbol \varepsilon }= [\varepsilon _{11},\varepsilon _{22}, \varepsilon _{33},\varepsilon _{12},\varepsilon _{13}, \varepsilon _{23}]^T \quad ,\quad  \delta {\boldsymbol \varepsilon }= [\delta \varepsilon _{11},\delta \varepsilon _{22}, \delta \varepsilon _{33},\delta \varepsilon _{12},\delta \varepsilon _{13}, \delta \varepsilon _{23}]^T\\[.5cm]  {\bf D}=\mu \left[                      \begin{array}{cccccc}4/3 & -2/3 & -2/3 & 0 & 0 & 0 \\                         &  4/3 & -2/3 &0  &0  & 0 \\                         &  & 4/3 & 0 & 0 & 0 \\                         &  &  & 2 & 0 & 0 \\                        {Sym.} &  &  &  & 2 & 0 \\                         &  &  &  &  & 2 \\                      \end{array} \right]\qquad ,\qquad {\bf m}=[1,1,1,0,0,0]^T \end{array} </math>

|}

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

|}

​

===4.2 Mass balance equations===

​

We multiply Eq.([[#eq-13|13]]) by arbitrary (continuous) test functions <math display="inline">q</math> (with dimensions of pressure) defined over the analysis domain <math display="inline">\Omega </math>. Integrating over <math display="inline">\Omega </math> gives

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\int _\Omega -\frac{q}{\kappa } \frac{Dp}{Dt}d\Omega -\int _\Omega q \frac{\tau }{c^2} \frac{D^2p}{Dt^2}d\Omega + \int _\Omega q \varepsilon _v d\Omega  + \int _\Omega q \tau {\partial \hat r_{m_i} \over \partial x_i} d\Omega =0  </math>

|}

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

|}

​

Integrating by parts the last integral in Eq.([[#eq-21|21]]) and using ([[#eq-15|15]]) gives after some transformations <span id='citeF-39'></span><span id='citeF-40'></span>[[#cite-39|[39,40]]]

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{r}\displaystyle \int _\Omega \frac{q}{\kappa } \frac{Dp}{Dt}d\Omega + \int _\Omega q \frac{\tau }{c^2} \frac{D^2p}{Dt^2}d\Omega - \int _\Omega q \varepsilon _v d\Omega + \int _\Omega \tau {\partial q \over \partial x_i} \left({\partial  \over \partial x_i} (2\mu \varepsilon _{ij})+ {\partial p \over \partial x_i}+b_i\right)d\Omega \\ \displaystyle - \int _{\Gamma _t} q \tau \left[\rho \frac{Dv_n}{Dt}-\frac{2}{h_n} \left(2\mu {\partial v_n \over \partial n} + p -t_n\right)\right]d\Gamma =0 \end{array}   </math>

|}

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

|}

​

Expression ([[#eq-24|24]]) holds for 2D and 3D problems. The terms involving the first and second material time derivative of the pressure and the boundary term in Eq.([[#eq-24|24]]) are important for preserving the conservation of mass in  free-surface flow problems <span id='citeF-10'></span><span id='citeF-41'></span>[[#cite-10|[10,41]]].

​

===4.3 Thermal balance equation===

​

Application of the standard weighted residual method to the heat balance equations ([[#eq-5|5]]) and ([[#eq-9|9]]) leads, after standard operations, to <span id='citeF-7'></span><span id='citeF-44'></span>[[#cite-7|[7,44]]]

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\int _\Omega \hat w \rho c {\partial T \over \partial t} d\Omega + \int _\Omega {\partial \hat w  \over \partial x_i}  k {\partial T \over \partial x_i}d\Omega - \int _\Omega \hat w Q d\Omega + \int _{\Gamma _q} \hat w q_n^p d\Gamma =0 </math>

|}

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

|}

​

where <math display="inline">\hat w</math> are the space weighting functions for the temperature.

​

==5 FEM DISCRETIZATION==

​

We discretize  the analysis domain into finite elements with <math display="inline">n</math> nodes in the standard manner leading to a mesh with a total number of <math display="inline">N_e</math> elements and <math display="inline">N</math> nodes. In our work we will choose simple 3-noded linear triangles (<math display="inline">n=3</math>) for 2D problems and 4-noded tetrahedra (<math display="inline">n=4</math>) for 3D problems with local linear shape functions <math display="inline">N_i^e</math> defined for each node <math display="inline">i</math> (<math display="inline">i=1,n</math>) of element <math display="inline">e</math> <span id='citeF-34'></span><span id='citeF-44'></span>[[#cite-34|[34,44]]]. The velocity components, the  pressure and the temperature are interpolated over the mesh in terms of their nodal values in the same manner using the  global linear shape functions <math display="inline">N_j</math> spanning over the elements sharing node <math display="inline">j</math> (<math display="inline">j=1,N</math>)  <span id='citeF-34'></span><span id='citeF-44'></span><span id='citeF-46'></span>[[#cite-34|[34,44,46]]]. In matrix form

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>{\bf v} = {\bf N}_v\bar {\bf v}  ~,~p = {\bf N}_p\bar {\bf p}\quad  , \quad T = {\bf N}_T \bar {\bf T}  </math>

|}

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

|}

​

where

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{c}\displaystyle \bar {\bf v} = \left\{\begin{matrix}\bar{\bf v}^1\\\bar{\bf v}^2\\\vdots \\ \bar{\bf v}^N\end{matrix}  \right\}\quad \hbox{with}~ \bar{\bf v}^i = \left\{\begin{matrix}\bar{v}^i_1\\\bar{v}^i_2\\\bar{v}^i_3\end{matrix}  \right\}~, ~ \bar {\bf p} =\left\{\begin{matrix}\bar{p}^1\\\bar{p}^2\\\vdots \\ \bar{p}^{N}\end{matrix}  \right\}~, ~ {\bf T}= \left\{\begin{matrix}\bar{T}_1\\\bar{T}_2\\\vdots \\ \bar{T}^N\end{matrix}  \right\}\\ \displaystyle {\bf N}_v = [{\bf N}_1, {\bf N}_2,\cdots , {\bf N}_N ]^T ~,~ {\bf N}_p = {\bf N}_T = [{N}_1, { N}_2,\cdots , {N}_N ] \end{array} </math>

|}

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

|}

​

with <math display="inline">{\bf N}_j = N_j {\bf I}_n</math> where <math display="inline">{\bf I}_n</math> is the <math display="inline">n\times n</math> unit matrix.

​

In Eq.([[#eq-25|25]]) vectors <math display="inline">\bar{v}</math>, <math display="inline">\bar {p}</math> and <math display="inline">\bar {T}</math>  contain the nodal velocities, the nodal pressures and the nodal temperatures for the whole mesh, respectively and the upperindex denotes the nodal value for each vector or scalar magnitude.

​

Substituting Eqs.([[#eq-24|24]]) into Eqs.([[#eq-17|17]]), ([[#eq-22|22]]) and ([[#eq-25|25]]) and choosing a Galerkin formulation with <math display="inline">w_i = q =  \hat w_i =N_i</math> leads to the following system of algebraic equations

​

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

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>{\bf M}_0 {\dot{\bar{\bf v}}} + {\bf K}\bar{\bf v}+{\bf Q}\bar {\bf p}- {\bf f}_v={\bf 0} </math>

|}

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

|}

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\textbf{M}_1 {\dot{\bar{\bf p}}}+\textbf{M}_2 {\ddot{\bar{\bf p}}}-{\bf Q}^T \bar{\bf v} + ({\bf L}+{\bf M}_b)\bar {\bf p}- {\bf f}_p={\bf 0} </math>

|}

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

|}

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>{\bf C} {\dot{\bar{\bf T}}}+\hat {\bf L} \bar{\bf T} -  {\bf f}_T={\bf 0}  </math>

|}

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

|}

​

where <math display="inline">{\dot{\bar{\bf a}}}</math> and <math display="inline">{\ddot{\bar{\bf a}}}</math> denote the first and second material time derivatives of the components of  a vector <math display="inline">{a}</math>. The different matrices and vectors in Eqs.([[#eq-26|26]]) are assembled from the element contributions given in Box 1.

​

<br/>

​

'''Remark 2.''' The boundary terms of vector <math display="inline">{\bf f}_p</math> can be incorporated in the matrices of Eq.([[#eq-26b|26b]]). This leads to a non symmetrical set of equations. These boundary terms are computed here iteratively within the incremental solution scheme.

​

'''Remark 3.''' The presence of matrix <math display="inline">{\bf M}_b</math> in Eq.([[#eq-26b|26b]]) allows us to compute the pressure without the need of prescribing its value at the free surface. This eliminates the error introduced when the pressure is prescribed to zero in free boundaries, which leads to considerable mass losses in viscous flows <span id='citeF-15'></span>[[#cite-15|[15]]].

​

'''Remark 4.''' For ''transient problems'' the stabilization parameter <math display="inline">\tau </math> of Eq.([[#eq-14|14]]) is computed for each element <math display="inline">e</math> using <math display="inline">h=l^e</math> and <math display="inline">\delta = \Delta t</math> as

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\tau = \left(\frac{8\mu }{(l^e)^2} + \frac{2\rho }{\Delta t}  \right)^{-1} </math>

|}

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

|}

​

where <math display="inline">\Delta t</math> is the time step used for the transient solution and <math display="inline">l^e</math> is a characteristic element length computed as <math display="inline">l^e = 2(\Omega ^e)^{1/n_s}</math> where <math display="inline">\Omega ^e</math> is the element area (for 3-noded triangles) or volume (for 4-noded tetrahedra). For fluids with heterogeneous material the values of <math display="inline">\mu </math> and <math display="inline">\rho </math>  are computed at the element center.

​

For ''steady state'' problems the stabilization parameter is computed with Eq.([[#eq-27|27]]) substituting <math display="inline">\Delta t</math> by <math display="inline">\frac{l^e}{\vert {\bf v}^e\vert }</math>.  The characteristic boundary length <math display="inline">h_n</math> in the expression of <math display="inline">{\bf f}_p</math> (Box 1) has been taken equal to <math display="inline">l^e</math> in our computations.

​

{|  class="floating_tableSCP wikitable" style="text-align: left; margin: 1em auto;min-width:50%;"

|-

| style="border-left: 2px solid;border-right: 2px solid;border-top: 2px solid;" | <math>\begin{array}{l}  \displaystyle \mathbf{M}_{0_{ij}}^e  =\int _{\Omega ^e}\rho {N}_{i}^e {N}_{j} {\bf I}_3 d\Omega ~,~\mathbf{K}^e_{ij} =\int _{\Omega ^e} \mathbf{B}_i^{eT} \mathbf{D}\mathbf{B}_j^e d\Omega ~,~\mathbf{Q}^e_{ij}  =\int _{\Omega ^e} \mathbf{B}_i^{eT} \mathbf{m} {N}_j^e d\Omega \\[.4cm] \displaystyle {M}_{1_{ij}}^e =\int _{\Omega ^e} \frac{1}{\kappa }{N}_i^e {N}_j^e d\Omega ~,~ {M}_{2_{ij}}^e =\int _{\Omega ^e} \frac{\tau }{ c^2 }{N}_i^e{N}_j^e d\Omega ~,~ {M}_{b_{ij}}^e = \int _{\Gamma _t} \frac{2\tau }{h_n} {N}_i^e{N}_j^e  d\Gamma \\[.4cm] 

\displaystyle {L}^e_{ij}= \int _{\Omega ^e} \tau ({\boldsymbol \nabla }^T {N}_i^e) {\boldsymbol \nabla }  N_j^e d\Omega ~,~ \mathbf{f}^e_{v_i}=  \int _{\Omega ^e}\mathbf{N}_{i}^e \mathbf{b} d\Omega + \int _{\Gamma _t}  \mathbf{N}_{i}^e {\bf t}^p d\Gamma \\[.4cm] 

\displaystyle {f}^e_{p_i}=\int _{\Gamma _t}\tau N_i^e \left[\rho \frac{Dv_n}{Dt}-\frac{2}{h_n} (2\mu \varepsilon _n -t_n)\right]d\Gamma - \int _{\Omega ^e} \tau  {\boldsymbol \nabla }^T {N}_i^e {b} d\Omega \\[.4cm] 

\displaystyle {C}^e_{ij}= \int _{\Omega ^e} \rho c N_i^e N_j^e d\Omega ~,~ \hat L^e_{ij}= \int _{\Omega ^e} k ({\boldsymbol \nabla }^T {N}_i^e) {\boldsymbol \nabla } N_j^e d\Omega ~,~ \displaystyle {f}^e_{T_i}= \int _{\Omega ^e} N_i^e Q d\Omega - \int _{\Gamma _q} N_i^e q^p_n d\Gamma \\[.4cm] \hbox{with } i,j=1,n.\\ \hbox{For 3D problems}\\[.4cm] 

\displaystyle{\bf B}_i^e = \left[\begin{matrix}  \displaystyle {\partial N_i^e \over \partial x_1} &0&0\\ \displaystyle{0}& \displaystyle {\partial N_i^e \over \partial x_2}&0\\ \displaystyle{0}&0&\displaystyle {\partial N_i^e \over \partial x_3}\\ \displaystyle {\partial N_i^e \over \partial x_2}&\displaystyle {\partial N_i^e \over \partial x_1}&0\\[.25cm] \displaystyle {\partial N_i^e \over \partial x_3}&0&\displaystyle {\partial N_i^e \over \partial x_1}\\[.25cm] \displaystyle{0}&\displaystyle {\partial N_i^e \over \partial x_3}&\displaystyle {\partial N_i^e \over \partial x_2}          \end{matrix}  \right]~,~ \mathbf{\bf N}_{i}^e = N_i^e {\bf I}_3 \quad \hbox{and} \quad  {\boldsymbol \nabla } = \left\{\begin{matrix}\displaystyle {\partial  \over \partial x_1} \\ \displaystyle {\partial  \over \partial x_2}\\\displaystyle {\partial  \over \partial x_3}  \end{matrix}  \right\} \end{array}</math>

|- style="border-bottom: 2px solid;"

| style="border-left: 2px solid;border-right: 2px solid;" | <math>N_i^e</math>: Local shape function of node <math display="inline">i</math> of element <math display="inline">e</math> [34,44]

|}

​

<div class="center" style="font-size: 75%;">'''Box 1'''. Element form of the matrices and vectors in Eqs.(26)</div>

​

==6 TRANSIENT SOLUTION OF THE DISCRETIZED EQUATIONS==

​

Eqs.([[#eq-26|26]]) are solved in time with an implicit Newton-Raphson type iterative scheme <span id='citeF-3'></span><span id='citeF-7'></span><span id='citeF-44'></span><span id='citeF-46'></span>[[#cite-3|[3,7,44,46]]]. The basic steps within a time increment <math display="inline">[n,n+1]</math> are:<br/>

​

* Initialize variables: <math display="inline">({}^{n+1}\hbox{\bf x}^1,{}^{n+1}\bar {\bf v}^1,{}^{n+1}\bar{\bf p}^1,  {}^{n+1}\bar{\bf T}^1, {}^{n+1}\bar{\bf r}^1_m)\equiv  \left\{{}^{n}{\bf x},{}^{n}\bar{\bf v},{}^{n}\bar{\bf p},{}^{n}\bar{\bf T}, {}^{n}\bar{\bf r}_m \right\}</math>.

* Iteration loop: <math display="inline">i=1,NITER</math>.<p> For each iteration. 

​

</p>

​

===''Step 1. Compute the nodal velocity increments <math>\Delta \bar {v}</math>''===

​

From Eq.([[#eq-26a|26a]]), we deduce

​

<span id="eq-28a"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>{}^{n+1}{\bf H}_v^i \Delta \bar {\bf v} = - {}^{n+1}\bar{{\bf r}}_m^i \rightarrow  \Delta \bar{\bf v} </math>

|}

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

|}

​

with the momentum residual <math display="inline">\bar{\bf r}_m</math> and the iteration matrix <math display="inline">{\bf H}_v</math>  given by

​

<span id="eq-28b"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>\bar{\textbf{r}}_m = \textbf{M}_0 {\dot{\bar{\bf v}}} + \textbf{K}\bar{\bf v}+\textbf{Q}\bar {\bf p}-\textbf{f}_v\quad , \quad \textbf{H}_v = \frac{1}{\Delta t} \textbf{M}_0 + \textbf{K} + \textbf{K}_v  </math>

|}

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

|}

​

===''Step 2. Update the nodal velocities''===

​

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

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

|-

| 

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

|-

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

|}

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

|}

​

===''Step 3. Compute the nodal pressures'' <math> {}^{n+1}\bar{\bf p}^{i+1}</math>===

​

From Eq.([[#eq-26b|26b]]) we obtain

​

<span id="eq-30a"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>{}^{n+1} \textbf{H}_p^i  {}^{n+1}\bar{\bf p}^{i+1}= \frac{1}{\Delta t} \textbf{M}_1 {}^{n+1}\bar{\bf p}^i +  \frac{1}{\Delta t^2} \textbf{M}_2 (2 ^n \bar{\bf p} - ^{n-1} \bar{\bf p}) + {\bf Q}^T {}^{n+1}\bar{\bf v}^{i+1}+ {}^{n+1}\bar{\bf f}^{i}_p  \rightarrow  {}^{n+1}\bar{\bf p}^{i+1}   </math>

|}

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

|}

​

with

​

<span id="eq-30b"></span>

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

|-

| 

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

|-

| style="text-align: center;" | <math>\textbf{H}_p = \frac{1}{\Delta t} \textbf{M}_1+ \frac{1}{\Delta t^2} \textbf{M}_2+\textbf{L} + \textbf{M}_b </math>

|}

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

|}

​

===''Step 4. Update the nodal coordinates''===

​

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

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

|-

| 

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

|-

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

|}

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

|}

​

A more accurate expression for computing <math display="inline">{}^{n+1}{\bf x}^{i+1}</math> can be used involving the  nodal accelerations <span id='citeF-40'></span>[[#cite-40|[40]]].

​

===''Step 5.Compute the nodal temperatures''===

​

From Eq.([[#eq-26c|26c]]) we obtain

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\left[\frac{1}{\Delta t}{\bf C} + \hat{\bf L}   \right]\Delta \bar{\bf T}=- {}^{n+1} \bar{\bf r}^i_T \quad ,\quad  {}^{n+1}\bar{\bf T}^{i+1}= {}^{n}\bar{\bf T}^i + \Delta \bar{\bf T} </math>

|}

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

|}

​

with

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>{\bf r}_T= {\bf C} {\dot{\bar{\bf T}}} + \hat{\bf L} \bar{\bf T} - {\bf f}_T </math>

|}

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

|}

​

===''Step 6. Check convergence''===

​

Verify the following conditions:

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\begin{array}{c}\displaystyle \Vert {}^{n+1}\bar {\bf v}^{i+1}- {}^{n+1}\bar{\bf v}^{i}\Vert \le e_v \Vert {}^{n}\bar {\bf v}\Vert \\ \displaystyle \Vert {}^{n+1}\bar {\bf p}^{i+1}-{}^{n+1}\bar {\bf p}^i\Vert \le e_p  \Vert {}^{n}\bar{\bf p}\Vert \\ \displaystyle \Vert {}^{n+1}\bar {\bf T}^{i+1} - {}^{n+1}\bar {\bf T}^i\Vert \le e_T  \Vert {}^{n}\bar{\bf T}\Vert  \end{array}   </math>

|}

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

|}

​

where <math display="inline">e_v</math>, <math display="inline">e_p</math> and <math display="inline">e_T</math> are prescribed error norms. In our examples we have set <math display="inline">e_v = e_p=e_T = 10^{-3}</math>.

​

If  conditions ([[#eq-34|34]]) are satisfied then make <math display="inline">n \leftarrow n+1 </math> and proceed to the next  time step. Otherwise, make the iteration counter <math display="inline">i \leftarrow i+1 </math> and repeat Steps 1&#8211;5.

​

'''Remark 5.''' In Eqs.([[#eq-28|28]])&#8211;([[#eq-34|34]]) <math display="inline">{}^{n+1}(\cdot )</math> denotes the values of a matrix or a vector computed using the nodal unknowns at time <math display="inline">n+1</math>. In this work the derivatives and integrals in all the matrices and the residual vectors <math display="inline">\bar{r}_m</math>  and <math display="inline">\bar{\bf r}_T</math> are computed on the ''discretized geometry at time'' <math display="inline">n</math>  while the nodal force vectors <math display="inline">{\bf f}_v</math>, <math display="inline">{\bf f}_p</math> and <math display="inline">{\bf f}_v</math> are computed on the current configuration at time <math display="inline">n+1</math>. This is equivalent to using an ''updated Lagrangian formulation'' <span id='citeF-3'></span><span id='citeF-45'></span><span id='citeF-46'></span>[[#cite-3|[3,45,46]]].

​

'''Remark 6.'''  Including the  bulk stiffness matrix <math display="inline">\textbf{K}_v</math> in <math display="inline">\textbf{H}_v</math> has proven to be essential for the fast convergence, mass preservation and overall accuracy of the iterative solution <span id='citeF-10'></span><span id='citeF-41'></span>[[#cite-10|[10,41]]].     The element expression of <math display="inline">\textbf{K}_v</math> can be obtained as <span id='citeF-41'></span>[[#cite-41|[41]]]

​

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

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

|-

| 

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

|-

| style="text-align: center;" | <math>\hbox{K}_v^e = \int _{\Omega ^e} \textbf{B}^T \texbf{m} \theta \Delta t\kappa \textbf{m}^T \textbf{B}d\Omega  </math>

|}

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

|}

​

where <math display="inline">\theta </math> is a positive number such that <math display="inline">0<   \theta \le 1</math> that has the role of preventing the ill-conditioning of the iteration matrix <math display="inline">{\bf H}_v</math> for highly incompressible fluids. An adequate selection of <math display="inline">\theta </math> also improves the overall accuracy of the numerical solution and the preservation of mass for large time steps <span id='citeF-10'></span>[[#cite-10|[10]]]. For fully incompressible fluids (<math display="inline">c</math> and <math display="inline">\kappa =\infty </math>), a finite  value of <math display="inline">\kappa </math> is used in practice in <math display="inline">{\bf K}_v</math> as this helps to obtaining an accurate solution for velocities and pressure  with reduced mass loss in few iterations per time step <span id='citeF-10'></span>[[#cite-10|[10]]]. These considerations, however, do not affect the value of <math display="inline">\kappa </math> within matrix <math display="inline">{\bf M}_1</math> in Eq.([[#eq-26b|26b]]) that vanishes for a fully incompressible fluid. Clearly, the value of the terms of <math display="inline">{\bf K}_v^e</math> can also be limited by reducing the time step size. This, however, leads to an increase in the cost of the computations. A similar approach for improving mass conservation in incompressible flows was proposed  in <span id='citeF-42'></span>[[#cite-42|[42]]].

​

'''Remark 7.''' The iteration matrix <math display="inline">\textbf{H}_v</math> in Eq.([[#eq-28a|28a]]) is an ''approximation of the exact tangent matrix'' in the updated Lagrangian formulation for a quasi-incompressible fluid <span id='citeF-40'></span>[[#cite-40|[40]]]. The simplified form of <math display="inline">\textbf{H}_v</math>  used in this work has yielded good results  with convergence achieved for the nodal velocities,  the pressure and the temperature in 3&#8211;4 iterations in all the problems analyzed.

​

'''Remark 8.''' The time step within a time interval <math display="inline">[n,n+1]</math> is chosen as <math display="inline">\Delta t =\min \left(\frac{^n l_{\min }^e}{\vert {}^n{\bf v}\vert _{\max }},\Delta t_b\right)</math> where <math display="inline">^n l_{\min }^e</math> is the minimum characteristic distance  of all elements in the mesh, with <math display="inline">l^e</math> computed as explained in Remark 4, <math display="inline">\vert ^n{\bf v}\vert _{\max }</math> is the maximum value of the modulus of the velocity of all nodes in the mesh  and <math display="inline">\Delta t_b</math> is the critical time step of all nodes approaching a solid boundary defined as <math display="inline">\Delta t_b = \min \left(\frac{^n l_b}{\vert ^n{\bf v}_b\vert _{\max }}\right)</math> where <math display="inline">^n l_b</math> is the distance from the node to the boundary and <math display="inline">^n{\bf v}_b</math> is the velocity of the node. This definition of <math display="inline">\Delta t</math> intends that no node crosses a solid boundary during a time step.

​

A method that allows using large time steps in the integration of the PFEM equations can be found in <span id='citeF-16'></span>[[#cite-16|[16]]].

​

==7 ABOUT THE PARTICLE FINITE ELEMENT METHOD (PFEM)==

​

===7.1 The basis of the PFEM===

​

Let us consider a domain <math display="inline">V</math> containing  fluid and solid subdomains. Each subdomain is characterized by a set of points, hereafter termed ''particles''. The particles contain all the information  for defining the geometry and the material and mechanical properties of the underlying subdomain. In the PFEM both subdomains are modelled using an ''updated'' ''Lagrangian formulation'' <span id='citeF-3'></span><span id='citeF-45'></span>[[#cite-3|[3,45]]].

​

The solution steps within a time step in the PFEM are as follows:

​

<div id='img-1'></div>

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

|-

|[[Image:draft_Samper_264539605-Figure2a.png|600px|Sequence of steps to update a “cloud” of nodes representing a domain containing a fluid and a solid part from time n   (t=ⁿt)  to   time n+2 (t=ⁿt +2∆t)  ]]

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

| colspan="1" | '''Figure 1:''' Sequence of steps to update a “cloud” of nodes representing a domain containing a fluid and a solid part from time <math>n</math>   (<math>t=^n t</math>)  to   time <math>n+2</math> (<math>t=^n t +2\Delta t</math>)  

|}

​

<ol>

​

<li>The starting point at each time step is the cloud of points <math display="inline">C</math> in the fluid and solid domains. For instance <math display="inline">^{n} C</math> denotes the cloud at time <math display="inline">t={}^n t </math> (Figure [[#img-1|1]]). </li>

​

<li>Identify the boundaries defining the analysis domain <math display="inline">^{n} V</math>, as well as the subdomains in the fluid and the solid. This is an essential step as some boundaries (such as the free surface in fluids) may be severely distorted during the solution, including separation and re-entering of nodes. The Alpha Shape method <span id='citeF-8'></span>[[#cite-8|[8]]] is used for the boundary definition. Clearly, the accuracy in the reconstruction of the boundaries depends on the number of points in the vicinity of each boundary and on the Alpha Shape parameter. In the problems solved in this work the Alpha Shape method has been implementation as described in <span id='citeF-12'></span><span id='citeF-28'></span>[[#cite-12|[12,28]]]. </li>

​

<li> Discretize the the analysis domain <math display="inline">^{n} V</math> with a finite element mesh <math display="inline">{}^{n} M.</math>We use an efficient mesh generation scheme based on an enhanced Delaunay tesselation <span id='citeF-11'></span><span id='citeF-12'></span>[[#cite-11|[11,12]]]. </li>

​

<li> Solve the Lagrangian equations of motion for the overall continuum using the standard FEM. Compute the state variables in at the next (updated) configuration for <math display="inline">^n t+\Delta t</math>: velocities, pressure and viscous stresses in the fluid and displacements, stresses and strains in the solid. </li>

​

<li>  Move the mesh nodes to a new position <math display="inline">^{n+1} C</math> where ''n''+1 denotes the time <math display="inline">{}^n t +\Delta t</math>, in terms of the time increment size. </li>

​

<li> Go back to step 1 and repeat the solution for the next time step to obtain <math display="inline">{}^{n+2} C</math>. </li>

​

</ol>

​

Note that the key differences between the PFEM and the classical FEM are the remeshing technique and the identification of the domain boundary at each time step.

​

The CPU time required for meshing grows linearly with the number of nodes. As a general rule,  meshing consumes for 3D problems around 15% of the total CPU time per time step, while the solution of the equations (with typically 3 iterations per time step) and the system assembly consume approximately 70% and 15% of the CPU time per time step, respectively. These figures refer to analyses in a single processor Pentium IV PC <span id='citeF-36'></span>[[#cite-36|[36]]]. Considerable speed can be gained using parallel computing techniques.

​

In this work we will apply the PFEM to problems involving a rigid domain containing fluid particles only. Application of the PFEM in fluid and solid mechanics and in fluid-structure interaction problems can be found in <span id='citeF-4'></span><span id='citeF-5'></span><span id='citeF-6'></span><span id='citeF-8'></span><span id='citeF-9'></span><span id='citeF-10'></span><span id='citeF-11'></span><span id='citeF-12'></span><span id='citeF-13'></span><span id='citeF-14'></span><span id='citeF-15'></span><span id='citeF-16'></span><span id='citeF-17'></span><span id='citeF-18'></span><span id='citeF-19'></span><span id='citeF-28'></span><span id='citeF-29'></span><span id='citeF-35'></span><span id='citeF-36'></span><span id='citeF-39'></span><span id='citeF-40'></span><span id='citeF-41'></span><span id='citeF-42'></span><span id='citeF-43'></span>[[#cite-4|[4,5,6,8,9,10,11,12,13,14,15,16,17,18,19,28,29,35,36,39,40,41,42,43]]], as well in www.cimne.com/pfem.

​

==8 EXAMPLES==

​

===8.1 Sloshing of water in prismatic tank===

​

The problem has been solved first in 2D. Figure [[#img-2|2]] shows the analysis data. The fluid oscillates due to the hydrostatic forces induced by its original position.

​

<div id='img-2'></div>

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

|-

|[[Image:draft_Samper_264539605-Figure2_SL2D_official_input.png|570px|2D analysis of sloshing of water in rectangular tank. Initial geometry, analysis data and mesh of 5064 3-noded triangles discretizing the water in the tank]]

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

| colspan="1" | '''Figure 2:''' 2D analysis of sloshing of water in rectangular tank. Initial geometry, analysis data and mesh of 5064 3-noded triangles discretizing the water in the tank

|}

​

The problem has been run using different values of the parameter <math display="inline">\theta </math> in the tangent bulk stiffness matrix <math display="inline">{K}_v^e</math> (Eq.([[#eq-41|41]])). The first set of results (Figures [[#img-3|3]] and [[#img-4|4]]) were obtained with <math display="inline">\theta =1</math>. The problem was then solved for <math display="inline">\theta = 0.08</math>, thereby, reducing in one order the magnitude the diagonal terms in <math display="inline">{K}_v^e</math>.