Casals 2018d

Error creating thumbnail: File missing

Eq. 6Eq. 6Eq. 6Estabilized

1. Introduction

Tdyn is a numerical simulation solver for multi-physics problem that uses the stabilized (FIC) finite element method. Using this stabilized numerical scheme, Tdyn is able to solve for instance heat transfer problems in both, fluids and solids, turbulent flows, advection of species and free surface problems. This document gives a short overview of the theoretical principles that Tdyn is based on.

2. Navier Stokes solver (ransol module)

2.1 Governing equations

The incompressible Navier-Stokes-Equations in a given three-dimensional domain Ω and time interval (0, t) can be written as:

$\left.{\begin{matrix}\rho \left({\frac {\partial {\mathit {\boldsymbol {u}}}}{\partial t}}+\left({\mathit {\boldsymbol {u}}}{\mathit {\boldsymbol {\cdot }}}\nabla \right){\mathit {\boldsymbol {u}}}\right)+\nabla p-\nabla \cdot \left(\mu \nabla {\mathit {\boldsymbol {u}}}\right)=\rho f\,&\,in\,\Omega \times \left(0,t\right)\\\nabla {\mathit {\boldsymbol {u}}}=0\,&\,in\,\Omega \times \left(0,t\right)\end{matrix}}\right\}$

(2-1)

where u = u (x, t) denotes the velocity vector, p = p (x, t) the pressure field, ρ the (constant) density, μ the dynamic viscosity of the fluid and f the volumetric acceleration. The above equations need to be combined with the following boundary conditions:

${\begin{matrix}{\mathit {\boldsymbol {u}}}={\mathit {\boldsymbol {u}}}_{c}&in\,{\Gamma }_{D}\times \left(0,t\right)\\p={p}_{c}&in\,{\Gamma }_{P}\times \left(0,t\right)\\{\mathit {\boldsymbol {n}}}\cdot \sigma \cdot {\mathit {\boldsymbol {g}}}_{1}=0,\,\,{\mathit {\boldsymbol {n}}}\cdot \sigma \cdot {\mathit {\boldsymbol {g}}}_{1}=0,\,\,{\mathit {\boldsymbol {n}}}\cdot {\mathit {\boldsymbol {u}}}={\mathit {\boldsymbol {u}}}_{M}&in\,{\Gamma }_{M}\times \left(0,T\right)\\{\mathit {\boldsymbol {u}}}\left({\mathit {\boldsymbol {x}}},0\right)={\mathit {\boldsymbol {u}}}_{0}({\mathit {\boldsymbol {x}}})&in\,{\Omega }_{D}\times \left(0\right)\\p\left({\mathit {\boldsymbol {x}}},0\right)={p}_{0}({\mathit {\boldsymbol {x}}})&in\,{\Omega }_{D}\times \left(0\right)\end{matrix}}$

(2-2)

In the above equations, ${\textstyle \Gamma :=\partial \Omega }$ , denotes the boundary of the domain ${\textstyle \Omega }$ , with n the normal unit vector, and g₁, g₂ the tangent vectors of the boundary surface ∂Ω. u_c is the velocity field on Γ_D (the part of the boundary of Dirichlet type, or prescribed velocity type), p_c the prescribed pressure on Γ_P (prescribed pressure boundary). σ is the stress field, ${\textstyle {\mathit {\boldsymbol {u}}}_{M}}$ the value of the normal velocity and u₀, p₀ the initial velocity and pressure fields. The union of Γ_D, Γ_P and Γ_M must be Γ; their intersection must be empty, as a point of the boundary can only be part of one of the boundary types, unless it is part of the border between two of them.

The spatial discretization of the Navier-Stokes equations has been done by means of the finite element method, while for the time discretization an iterative algorithm that can be consider as an implicit two steps "Fractional Step Method" has been used. Problems with dominating convection are stabilized by the so-called "Finite Increment Calculus" method, presented below.

2.2 CFD algorithms and stability

Using the standard Galerkin method to discretize the incompressible Navier-Stokes equations leads to numerical instabilities coming from two sources. Firstly, owing to the advective-diffusive character of the governing equations, oscillations appear in the solution at high Reynolds numbers, when the convection terms become dominant. Secondly, the mixed character of the equations limits the choice of velocity and pressure interpolations such that equal order interpolations cannot be used.

In recent years, a lot of effort has been put into looking for ways to stabilize the governing equations, many of which involve artificially adding terms to the equations to balance the convection, for example, the artificial diffusion method [7].

A new stabilization method, known as finite increment calculus, has recently been developed [5, 6]. By considering the balance of flux over a finite sized domain, higher order terms naturally appear in the governing equations, which supply the necessary stability for a classical Galerkin finite element discretization to be used with equal order velocity and pressure interpolations. This method has been Christianized “Finite Calculus” (FIC).

2.2.1 Finite calculus (FIC)formulation

To demonstrate this method, consider the flux problem associated with the incompressible conservation of mass in a 2D source less finite sized domain defined by four nodes, as shown in Figure 1.

Error creating thumbnail: File missing

Figure 1 Mass balance in 2 dimensions

Consider the flux problem (mass balance) through the domain, taking the average of the nodal velocities for each surface.

${\frac {{h}_{y}}{2}}\left(u\left(x-{h}_{x}\,,y\right)+u\left(x-{h}_{x},y-{h}_{y}\right)-\right.$$\left.u\left(x,y\right)+u\left(x,y-{h}_{y}\right)\right)+{\frac {{h}_{x}}{2}}\left(u\left(x,y-\right.\right.$$\left.\left.{h}_{y}\right)+u\left(x-{h}_{x},y-{h}_{y}\right)-u\left(x-{h}_{x},y\right)+\right.$$\left.u\left(x,y\right)\right)=0$

(2-3)

Now expand these velocities using a Taylor expansion, retaining up to second order terms. Write u (x,y) = u:

$u\left(x-{h}_{x},y\right)=u-{h}_{x}{\frac {\partial u}{\partial x}}+{\frac {{h}_{x}^{2}}{2}}{\frac {{\partial }^{2}u}{\partial {x}^{2}}}-$$O({h}_{x}^{3})\,$

(2-4)

and similarly for the other two components of the velocity.

Substituting this back into the original flux balance equation gives, after simplification, the following stabilized form of the 2D mass balance equation:

${\frac {\partial u}{\partial x}}+{\frac {\partial y}{\partial y}}-{\overline {{\frac {{h}_{x}}{2}}\left({\frac {{\partial }^{2}u}{\partial {x}^{2}}}-{\frac {{\partial }^{2}v}{\partial x\partial y}}\right)}}-$${\overline {{\frac {{h}_{y}}{2}}\left({\frac {{\partial }^{2}u}{\partial x\partial y}}-{\frac {{\partial }^{2}v}{\partial {y}^{2}}}\right)}}=$$0$

(2-5)

Note that as the domain size tends to zero, i.e. 72px , the standard form of the incompressible mass conservation equation in 2D is recovered. The underlined terms in the equation provide the necessary stabilization to allow a standard Galerkin finite element discretization to be applied. They come from admitting that the standard form of the equations is an unreachable limit in the finite element discretization, i.e. by admitting that the element size cannot tend to zero, which is the basis for the finite element method. It also allows equal order interpolations of velocity and pressure to be used.

2.2.2 Stabilized Navier-Stokes equations

The Finite increment Calculus methodology presented above is used to formulate stabilized forms of the momentum balance and mass balance equations and the Neumann boundary conditions. The velocity and pressure fields of an incompressible fluid moving in a domain Ω R^d(d=2,3) can be described by the incompressible Navier Stokes equations:

$\left.{\begin{matrix}\rho {\frac {\partial {u}_{i}}{\partial t}}+\rho {\frac {\partial {u}_{i}{u}_{j}}{\partial {x}_{j}}}+{\frac {\partial p}{\partial {x}_{i}}}-{\frac {\partial {s}_{ij}}{\partial {x}_{j}}}=\rho {f}_{i}\\{\frac {\partial {u}_{i}}{\partial {x}_{i}}}=0\end{matrix}}\right\}\,i,j=$$1,N$

(2-6)

Where 1 ≤ i, j ≤ d, ρ is the fluid density field, u_i is the ith component of the velocity field u in the global reference system x_i, p is the pressure field and s_ij is the viscous stress tensor defined by:

${s}_{ij}=2\nu \left({\epsilon }_{ij}-{\frac {1}{3}}{\frac {\partial {u}_{k}}{\partial {x}_{k}}}{\delta }_{ij}\right),\,{\quad \epsilon }_{ij}=$${\frac {1}{2}}\left({\frac {\partial {u}_{i}}{\partial {x}_{j}}}+{\frac {\partial {u}_{j}}{\partial {x}_{i}}}\right)$

(2-7)

The stabilized FIC form of the governing differential equation (2-6) can be written as

${r}_{{m}_{i}}-{\frac {1}{2}}{h}_{ij}^{2}{\frac {\partial {r}_{{m}_{i}}}{\partial {x}_{j}}}=$$0\quad in\,\Omega$

(2-8)

${r}_{d}-{\frac {1}{2}}{h}_{j}^{d}{\frac {\partial {r}_{d}}{\partial {x}_{j}}}=0\quad in\,\Omega$

(2-9)

Summation convention for repeated indices in products and derivatives is used unless otherwise specified. In above equations, terms $r_{m_{i}}$ and $r_{d}$ denote the residual of Eq. (2-6) , and $h_{ij}^{m}$ , $h_{j}^{d}$ are the characteristic length distances, representing the dimensions of the finite domain where balance of mass and momentum is enforced. Details on obtaining the FIC stabilized equations and recommendation for the calculation of the stabilization terms can be found in Oñate 1998.

Let n be the unit outward normal to the boundary ∂Ω, split into two sets of disjoint components Γ_t, Γ_u where the Neumann and Dirichlet boundary conditions for the velocity are prescribed respectively. The boundary conditions for the stabilized problem to be considered are (Oñate 2004):

$\left.{\begin{matrix}{n}_{j}{\sigma }_{ij}-{t}_{i}+{\frac {1}{2}}{h}_{ij}^{m}{n}_{j}{r}_{{m}_{i}}=0&\,on\,{\Gamma }_{t}\\{u}_{j}={u}_{j}^{p}&\,on\,{\Gamma }_{u}\end{matrix}}\right\}$

(2-10)

Where $t_{i},u_{j}^{p}$ are the prescribed surface tractions and ${\sigma }_{ij}$ is the total stress tensor, defined as

${\sigma }_{ij}={s}_{ij}-p{\delta }_{ij}$

(2-11)

Equations (2-8) - (2-10) constitute the starting point for deriving stabilized FEM for solving the incompressible Navier-Stokes equations. An interesting feature of the FIC formulation is that it allows using equal order interpolation for the velocity and pressure variables (García Espinosa 2003, Oñate 2004).

2.2.3 Stabilized integral forms

From the momentum equations, the following relations can be obtained (Oñate 2004)

$\left.{\begin{matrix}{\frac {\partial {r}_{d}}{\partial {x}_{i}}}={\frac {{h}_{ii}^{m}}{2{a}_{i}}}{\frac {\partial {r}_{{m}_{i}}}{\partial {x}_{j}}}\\{a}_{i}={\frac {2\mu }{3}}+{\frac {{u}_{i}{h}_{i}^{d}}{2}}\,\end{matrix}}\right\}\quad no\,sum\,in\,i$

(2-12)

Substituting the equation above into Eq. (2-9) and retaining only the terms involving the derivatives of ${\textstyle r_{m_{i}}}$ with respect to ${\textstyle x_{i}}$ , leads to the following alternative expression for the stabilized mass balance equation

${r}_{d}-\sum _{i=1}^{{n}_{d}}{\tau }_{i}{\frac {\partial {r}_{{m}_{i}}}{\partial {x}_{i}}}=$$0,\quad {\tau }_{i}={\left({\frac {8\mu }{3{h}_{jj}^{m}{h}_{j}^{d}}}+{\frac {2\rho {u}_{i}}{{h}_{ii}^{m}}}\right)}^{-1}\quad \,no\,sum\,in\,i$

(2-13)

The ${\tau }_{i}$ ’s in Eq. (2-13) when multiplied by the density are equivalent to the intrinsic time parameters, seen extensively in the stabilization literature. The interest of Eq. (2-13) is that it introduces the first space derivatives of the momentum equations into the mass balance equation. These terms have intrinsic good stability properties as explained next.

The weighted residual form of the momentum and mass balance equations is written as

$\left.{\begin{matrix}\int _{\Omega }^{\,}{\nu }_{i}\left[{r}_{{m}_{i}}-{\frac {1}{2}}{h}_{j}{\frac {\partial {r}_{{m}_{i}}}{\partial {x}_{j}}}\right]d\Omega +\int _{{\Gamma }_{i}}^{\,}{\nu }_{i}\left[{n}_{j}{\sigma }_{ij}-{t}_{i}+{\frac {1}{2}}{h}_{j}{n}_{j}{r}_{{m}_{i}}\right]d\Gamma =0\\\int _{\Omega }^{\,}q\left[{r}_{d}-\sum _{i=1}^{N}{\tau }_{i}{\frac {\partial {r}_{{m}_{i}}}{\partial {x}_{i}}}\right]d\Omega =0\end{matrix}}\right\}$

(2-14)

where ${\textstyle v_{i},q}$ are generic weighting functions. Integrating by parts the residual terms in the above equations leads to the following weighted residual form of the momentum and mass balance equations:

$\left.{\begin{matrix}\int _{\Omega }^{\,}{\nu }_{i}{r}_{{m}_{i}}d\Omega +\int _{{\Gamma }_{i}}^{\,}{\nu }_{i}\left({n}_{j}{\sigma }_{ij}-{t}_{i}\right)d\Gamma +\int _{\Omega }^{\,}{\frac {1}{2}}{h}_{j}{\frac {\partial {\nu }_{i}}{\partial {x}_{j}}}{r}_{{m}_{i}}d\Omega =0\\\int _{\Omega }^{\,}q{r}_{d}d\Omega +\int _{\Omega }^{\,}\left[\sum _{i=1}^{N}{\tau }_{i}{\frac {\partial q}{\partial {x}_{i}}}{r}_{{m}_{i}}\right]d\Omega -\int _{\Gamma }^{\,}\left[\sum _{i=1}^{N}q{r}_{i}{n}_{i}{r}_{{m}_{i}}\right]d\Gamma =0\end{matrix}}\right\}$

(2-15)

We will neglect here onwards the third integral in Eq. (2-15) by assuming that ${\textstyle r_{m_{i}}}$ is negligible on the boundaries. The deviatoric stresses and the pressure terms in the second integral of Eq. (2-15) are integrated by parts in the usual manner. The resulting momentum and mass balance equations are:

$\int _{\Omega }^{\,}{\nu }_{i}\rho \left({\frac {\partial {u}_{i}}{\partial t}}+{u}_{j}{\frac {\partial {u}_{i}}{\partial {x}_{j}}}\right)+{\frac {\partial {\nu }_{i}}{\partial {x}_{j}}}\left(\mu {\frac {\partial {u}_{i}}{\partial {x}_{j}}}-{\delta }_{ij}p\right)d\Omega -\int _{\Omega }^{\,}{\nu }_{i}\rho {f}_{i}d\Omega -\int _{\Gamma }^{\,}{\nu }_{i}{t}_{i}d\Gamma +$$\int _{\Omega }^{\,}{\frac {{h}_{j}}{2}}{\frac {\partial {\nu }_{i}}{\partial {x}_{j}}}\,{r}_{{m}_{i}}d\Omega =$$0$

(2-16)

and

$\int _{\Omega }^{\,}q{\frac {\partial {u}_{i}}{\partial {x}_{i}}}d\Omega +\int _{\Omega }^{\,}\left[\sum _{i=1}^{N}{\tau }_{i}{\frac {\partial q}{\partial {x}_{i}}}{r}_{{m}_{i}}\right]d\Omega \,=$$0$

(2-17)

In the derivation of the viscous term in Eq. (2-16) we have used the following identity holding for incompressible fluids (prior to the integration by parts)

${\frac {\partial {s}_{ij}}{\partial {x}_{j}}}=2\mu {\frac {\partial {\epsilon }_{ij}}{\partial {x}_{j}}}=$$\mu \,{\partial }^{2}{u}_{i}/\partial {x}_{j}\partial {x}_{j}$

(2-18)

2.2.4 Convective and pressure gradient projections

The computation of the residual terms are simplified if we introduce the convective and pressure gradient projections ${\textstyle c_{i}}$ and ${\textstyle {\pi }_{i}}$ , respectively defined as

${\begin{matrix}{c}_{i}={r}_{{m}_{i}}-{u}_{j}{\frac {\partial {u}_{i}}{\partial {x}_{j}}}\\{\pi }_{i}={r}_{{m}_{i}}-{\frac {\partial p}{\partial {x}_{i}}}\end{matrix}}$

(2-19)

We can express the residual terms in Eq. (2-16) and Eq. (2-17) in terms of ${\textstyle c_{i}}$ and ${\textstyle {\pi }_{i}}$ , respectively which then become additional variables. The system of integral equations is now augmented in the necessary number by imposing that the residuals vanish (in average sense). This gives the final system of governing equations as:

	$\int _{\Omega }^{\,}{\nu }_{i}\rho \left({\frac {\partial {u}_{i}}{\partial t}}+\,{u}_{j}{\frac {\partial {u}_{i}}{\partial {x}_{j}}}\right)+{\frac {\partial {\nu }_{i}}{\partial {x}_{j}}}\left(\mu {\frac {\partial {u}_{i}}{\partial {x}_{j}}}-{\delta }_{ij}p\right)d\Omega -\,\int _{\Omega }^{\,}{\nu }_{i}\rho {f}_{i}d\Omega -\,-\,\int _{{\Gamma }_{t}}^{\,}{\nu }_{i}{t}_{i}d\Gamma +\,\int _{\Omega }^{\,}{\frac {{h}_{k}}{2}}\rho {\frac {\partial {\nu }_{i}}{\partial {x}_{k}}}\left({u}_{j}{\frac {\partial {u}_{i}}{\partial {x}_{j}}}+{c}_{i}\right)d\Omega =0\,$
	$\int _{\Omega }^{\,}q{\frac {\partial {u}_{i}}{\partial {x}_{i}}}d\Omega +\,\int _{\Omega }^{\,}\left[\sum _{i=1}^{N}{\tau }_{i}{\frac {\partial q}{\partial {x}_{i}}}\left({\frac {\partial p}{\partial {x}_{i}}}+{\pi }_{i}\right)\right]d\Omega =$$0$	(2-20)
	$\int _{\Omega }^{\,}{b}_{i}\,\left({u}_{j}{\frac {\partial {u}_{i}}{\partial {x}_{j}}}+{c}_{i}\right)d\Omega =$$0\quad \quad \quad \,\,no\,sum\,in\,i$
	$\int _{\Omega }^{\,}{w}_{i}\,\left({\frac {\partial {p}_{i}}{\partial {x}_{j}}}+{\pi }_{i}\right)d\Omega =$$0\quad \quad \quad \quad \,no\,sum\,in\,i$

Being i,j,k = 1,N and b_i, w_i the appropriate weighting functions.

2.2.5 Monolithic time integration scheme

In this section an implicit monolithic time integration scheme, based on a predictor corrector scheme for the integration of Eq. (2-20) is presented.

Let us first discretize in time the stabilized momentum Eq. (2-8), using the trapezoidal rule (or θ method) as (see Zienkiewicz 1995):

$\rho \left({\frac {{{u}_{i}^{n+1}-u}_{i}^{n}}{\delta t}}+{\frac {\partial }{\partial {x}_{i}}}\,{\left({u}_{i}{u}_{j}\right)}^{n+\theta }\right)+$${\frac {\partial {p}^{n+1}}{\partial {x}_{i}}}-\,{\frac {\partial {s}_{ij}^{n+\theta }}{\partial {x}_{j}}}-$$\,\rho {f}_{i}^{n+\theta }-{\frac {1}{2}}\rho {h}_{mj}{\frac {\partial }{\partial {x}_{j}}}\,\left({u}_{j}^{n+\theta }{\frac {\partial {u}_{i}^{n+\theta }}{\partial {x}_{j}}}+\right.$$\left.\,{c}_{i}^{n+\theta }\right)=0$

(2-21)

where superscripts n and θ refer to the time step and to the trapezoidal rule discretization parameter, respectively. For θ = 1 the standard backward Euler scheme is obtained, which has a temporal error of 0(t). The value θ = 0.5 gives a standard Crank Nicholson scheme, which is second order accurate in time 0(t²).

An implicit fractional step method can be simply derived by splitting Eq. (2-21) as follows:

	$\rho \left({\frac {{u}_{i}^{\ast ,n+1}-{u}_{i}^{n}}{\delta t}}+{\frac {\partial }{\partial {x}_{i}}}\,{\left({u}_{i}{u}_{j}\right)}^{n+\theta }\right)+$${\frac {\partial {p}^{n+1}}{\partial {x}_{i}}}-{\frac {\partial {s}_{ij}^{n+\theta }}{\partial {x}_{j}}}-$$\,\rho {f}_{i}^{n+\theta }-{\frac {1}{2}}\,\rho {h}_{mj}{\frac {\partial }{\partial {x}_{j}}}\,\left({u}_{j}^{n+\theta }{\frac {\partial {u}_{i}^{n+\theta }}{\partial {x}_{j}}}+\right.$$\left.\,{c}_{i}^{n+\theta }\right)=0$	(2-22)
	${u}_{i}^{n+1}=\,{u}_{i}^{\ast ,n+1}-{\frac {\delta t}{\rho }}\left({\frac {\partial {p}^{n+1}}{\partial {x}_{i}}}-\right.$$\left.{\frac {\partial {p}^{n}}{\partial {x}_{i}}}\right)$

On the other hand, substituting the last equation above into Eq. (2-9) and after some algebra leads to the following alternative mass balance equation

$\rho {\frac {\partial {u}_{i}^{\ast }}{\partial {x}_{i}}}-\,\delta t{\frac {{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}}\left({p}^{n+1}-\right.$$\left.{p}^{n}\right)+\,\sum _{i=1}^{N}{\tau }_{i}{\frac {\partial }{\partial {x}_{i}}}\,\left({\frac {\partial {p}^{n+1}}{\partial {x}_{i}}}+\right.$$\left.\,{\pi }_{i}^{n+1}\right)=0$

(2-23)

The weighted residual form of the above equations can be written as follows:

	$\int _{\Omega }^{\,}{\nu }_{i}\rho \left({\frac {{u}_{i}^{\ast ,n+1}-{u}_{i}^{n}}{\Delta t}}+\right.$$\left.\,{u}_{j}^{\ast ,n+\theta }{\frac {\partial }{\partial {x}_{j}}}\,{u}_{i}^{\ast ,n+\theta }\right)d\Omega +$$\,\int _{\Omega }^{\,}\rho {\frac {\partial {\nu }_{i}}{\partial {x}_{j}}}\left({\frac {{\partial u}_{i}^{\ast ,n+\theta }}{\partial {x}_{j}}}-\right.$$\left.\,{\delta }_{ij}{p}^{n}\right)d\Omega -\,\int _{\Omega }^{\,}{\nu }_{i}\rho {f}_{i}^{n+\theta }d\Omega -$$\,\int _{{\Gamma }_{t}}^{\,}{\nu }_{i}{t}_{i}^{n+\theta }d\Gamma +\,\int _{\Omega }^{\,}{\frac {{h}_{k}}{2}}\rho {\frac {\partial {\nu }_{i}}{\partial {x}_{k}}}\left({u}_{j}^{\ast ,n+\theta }{\frac {{\partial u}_{i}^{\ast ,n+\theta }}{\partial {x}_{j}}}+\right.$$\left.{c}_{i}^{n}\right)d\Omega =0$	(2-24)
	$\int _{\Omega }^{\,}q\left(\rho {\frac {{\partial u}_{i}^{\ast }}{\partial {x}_{i}}}-\right.$$\left.\,\Delta t{\frac {{\partial }^{2}}{\partial {x}_{i}\partial {x}_{j}}}\,({p}^{n+1}-\right.$$\left.{p}^{n})\right)d\Omega +\,\int _{\Omega }^{\,}{\tau }_{i}{\frac {\partial q}{\partial {x}_{i}}}\left({\frac {{\partial p}_{\,}^{n+1}}{\partial {x}_{i}}}+\right.$$\left.\,{\pi }_{i}^{n+1}\right)d\Omega =0$
	$\int _{\Omega }^{\,}{\nu }_{i}\left({u}_{i}^{n+1}-{u}_{i}^{\ast }+{\frac {\Delta t}{\rho }}\,\left({\frac {\partial {p}^{n+1}}{\partial {x}_{i}}}-\right.\right.$$\left.\left.{\frac {\partial {p}^{n}}{\partial {x}_{i}}}\right)\right)d\Omega =$$0\quad \quad \,no\,sum\,in\,i$

At this point, it is important to introduce the associated matrix structure corresponding to the variational discrete FEM form of Eq. (2-24):

	$M{\frac {1}{\delta t}}{\overline {U}}^{n+1}-\,M{\frac {1}{\delta t}}\,{U}^{n}+\,K\left({\overline {U}}^{n+\theta }\right){\overline {U}}^{n+\theta }+$$\,{S}_{1}{\overline {U}}^{n+\theta }+\,{S}_{2}C-G{P}^{n}+\,{M}_{{\Gamma }_{t}}T=$$F$
	$\left(\delta tL+{L}^{{\tau }_{2}}\right){P}^{n+1}+\,{D}^{{\tau }_{2}}\Pi +{G}^{T}{\overline {U}}^{n+1}=$$\delta t\,L{P}^{n}$
	$M{\frac {1}{\delta t}}{U}^{n+1}-\,M{\frac {1}{\delta t}}\,{\overline {U}}^{n+1}-G{P}^{n+1}-$$\,G{P}^{n}=0$	(2-25)
	$N{U}^{n+1}+MC=0$
	$G{P}^{n+1}+M\Pi =0$

Where U, P are the vectors of the nodal velocity and pressure fields, T is the vector of prescribed tractions and Π, C the vectors of convective and pressure gradient projections. Terms denoted by over-bar identify the intermediate velocity obtained from the fractional momentum equation.

The system of equations above includes an error due to the splitting of the momentum equation. This error can be eliminated by considering the analogous system of equations

	$M{\frac {1}{\delta t}}{U}_{i+1}^{n+1}-\,M{\frac {1}{\delta t}}\,{U}^{n}+\,K\left({U}_{i}^{n+\theta }\right){U}_{i+1}^{n+\theta }+$$\,{S}_{1}{U}_{i+1}^{n+\theta }+\,{S}_{2}{C}_{i}-G{P}_{i}^{n+1}+\,{M}_{{\Gamma }_{t}}T=$$F$
	$\delta tL\left({P}_{i+1}^{n+1}-{P}_{i}^{n+1}\right)+\,{L}^{{\tau }_{2}}{P}_{i+1}^{n+1}+$${D}^{{\tau }_{2}}\Pi +{G}^{T}{U}_{i+1}^{n+1}=0$	(2-26)
	$N{U}_{i+1}^{n+1}+M{C}_{i+1}=0$
	$G{P}_{i+1}^{n+1}+M{\Pi }_{i+1}=0$

Where i is the iteration counter of the monolithic scheme. Basically, in this final formulation, the convergence of the resulting monolithic uncoupled scheme is enforced by the first term of the second equation of Eq. (2-26) (see Soto 2001).

2.2.6 Compressible flows

In the previous sections the basics of the incompressible flow solver implemented in Tdyn has been introduced. The extension of that algorithm to compressible flows is straightforward.

Eq. (2-23) can be easily modified to take into account the compressibility of the flow. The resulting equation is (see M. Vázquez, 1999)

$\rho {\frac {\partial {u}_{i}^{\ast }}{\partial {x}_{i}}}+{\frac {\alpha }{\delta t}}\,\left({p}^{n+1}-\right.$$\left.\,{p}^{n}\right)-\,\delta t{\frac {{\partial }^{2}}{\partial {x}_{i}\partial {x}_{i}}}\left({p}^{n+1}-\right.$$\left.{p}^{n}\right)+\,\sum _{i=1}^{N}{\tau }_{i}{\frac {\partial }{\partial {x}_{i}}}\left({\frac {\partial {p}^{n+1}}{\partial {x}_{i}}}+\right.$$\left.\,{\pi }_{i}^{n+1}\right)=f$

(2-27)

In the above equation the terms α and f depend on the compressibility law used.

2.3 Turbulence solvers

At high Reynolds numbers the flow will certainly be turbulent, and the resulting fluctuations in velocities need to be taken into account in the calculations. A process known as Reynolds averaging [3] is applied to the governing equations whereby the velocities, u_i are split into mean and a fluctuating component, where the fluctuating component, u_i’, is defined by:

${u}_{i}=\,{\overline {u}}_{i}+\,{u}_{i}^{'}$

(2-28)

In the expression above:

${\overline {u}}_{i}\left(x,t\right)={\frac {1}{T}}\,\int _{-T/2}^{T/2}{u}_{i}\left(x,t+\tau \right)d\tau \,$

(2-29)

This leads to extra terms in the governing equations that can be written as a function of ‘Reynolds stresses’ tensor, defined in Cartesian coordinates as:

${\tau }_{ij}^{R}=\,-\,{\overline {\rho {u}_{i}^{'}{u}_{j}^{'}}}\,$

(2-30)

Since the process of Reynolds averaging has introduced extra terms into the governing equations, we need extra information to solve the system of equations. To relate these terms to the other flow variables, a model of the turbulence is required. A large number of turbulence models exist of varying complexity, from the simple algebraic models to those based on two partial differential equations. The more complex models have increased accuracy at the expense of longer computational time. It is therefore important to use the simplest model that gives satisfactory results.

Several turbulence models as Smagorinsky, k, k-ε, k-ω, k-kt and Spalart Allmaras have been implemented in Tdyn. The final form of the so-called Reynolds Averaged Navier Stokes Equations (RANSE) using these models is:

${r}_{mi}-{\frac {{h}_{mj}}{2}}{\frac {\partial {r}_{mi}}{\partial {x}_{j}}}=$$0\quad \quad \quad on\,\Omega ,\,i,j=1,2,3.No\,sum\,in\,i$

(2-31)

${r}_{d}-{\frac {{h}_{dj}}{2}}{\frac {\partial {r}_{d}}{\partial {x}_{j}}}=$$0\quad \quad \quad on\,\Omega ,\,j=1,2,3.$

(2-32)

Where:

	${r}_{mi}=\rho \,\left[{\frac {\partial {u}_{i}}{\partial t}}+{\frac {\partial \left({u}_{i}{u}_{j}\right)}{\partial {x}_{j}}}\right]+$${\frac {\partial p}{\partial {x}_{i}}}-{\frac {\partial }{\partial {x}_{j}}}\left[\left(\mu +\right.\right.$$\left.\left.{\mu }_{T}\right){\frac {\partial {u}_{i}}{\partial {x}_{j}}}\right]-$$\,\rho {f}_{i}$	(2-33)
	${r}_{d}=\rho {\frac {\partial {u}_{i}}{\partial {x}_{i}}}$

being μ_T the so-called eddy viscosity.

The theory behind the above mentioned models may be found in the references, but a basic outline of the three most representative models is given below.

2.3.1 Zero equation (algebraic) model: the Smagorinsky model

Zero equation models are the most basic turbulence models. They assume that the turbulence is generated and dissipated in the same place, and so neglect the diffusion and convection of the turbulence. They are based on the idea of turbulent viscosity, and prescribe a turbulence viscosity, ν_t, either by means of empirical equations or experiment.

The Smagorinsky turbulence model is based on the idea of large eddy simulations (LES) in which the coherent large-scale structures are modelled directly within the computational mesh, whilst the small scales are modelled with the concept of eddy viscosity [4]. The Reynolds stress is written in terms of the turbulence kinetic energy, k, and eddy viscosity as:

${\overline {{u}_{i}^{'}{u}_{j}^{'}}}={\frac {2}{3}}k{\delta }_{ij}-\,2{\nu }_{t}{\epsilon }_{ij}$

(2-34)

where

$\epsilon ={\frac {1}{2}}{\left({\overline {\nabla }}\cdot {\overline {{u}^{'}}}+\,\left({\overline {\nabla }}\cdot {\overline {{u}^{'}}}\right)\right)}^{T}$

(2-35)

Smagorinsky proposed that the eddy viscosity depends on the mesh density and velocity gradients, and suggested the following expression:

${\nu }_{t}={\frac {{\mu }_{t}}{\rho }}=C{{h}^{e}}^{2}{\sqrt {{\epsilon }_{ij}{\epsilon }_{ij}}}$

(2-36)

where ν_T is the kinematic eddy viscosity, C is a constant of the order of 0.001, and h^e is the element size.

2.3.2 One-equation models: the kinetic energy model (k-model)

One-equation models attempt to model turbulent transport, by developing a differential equation for the transport of one of turbulent quantities. In the kinetic energy model, the velocity scale of the turbulence is taken as the square root of the turbulence kinetic energy, k. Kolmogorov and Prandtl independently suggested the following relationship between the eddy viscosity, this velocity scale, 24px , and length scale, L:

${\nu }_{t}={C}_{D}{\sqrt {kL}}$

(2-37)

where

$k={\frac {1}{2}}\,{\overline {{\overline {{u}^{'}}}^{2}}}={\frac {1}{2}}\,{\overline {{{u}^{'}}^{2}+\,{{v}^{'}}^{2}+\,{{w}^{'}}^{2}}}$

(2-38)

and C_D is an empirical constant, usually taken as 1.0.

Manipulation of the Reynolds and momentum equations leads to a differential equation for k [1]. The only gap in this analysis is that left by the length scale, L. This has to be specified either from experiment, or empirical equation. In the analysis used here, it has been specified as 1% of the smallest mesh size. Clearly this is the largest downfall of the analysis, and many consider it insufficiently accurate [3]. However, it is more accurate than the zero-equation models and less computationally expensive than the more accurate 2-equation models, and so it can be a good compromise between computational expense and accuracy, depending on the problem.

2.3.3 Two-equation models: the k-ε model

In order to increase the accuracy of our turbulence modelling, it is therefore necessary to develop a second differential transport equation to provide a complete system of closed equations without the need for empirical relationships. The k-ε turbulence model develops two such differential transport equations: one for the turbulence kinetic energy, k, and a second for the turbulent dissipation, ε. As for the k-model, the k-ε model relies on the Prandtl-Kolmogorov expression for the eddy viscosity above. From dimensional arguments, the turbulent dissipation, can be written in terms of the turbulence kinetic energy and the turbulence length scale, and the eddy viscosity written in terms of k and ε

	$\epsilon ={C}_{\epsilon }{\frac {{k}^{3/2}}{L}}$	(2-39)
	${\nu }_{t}={C}_{\nu }{\frac {{k}^{2}}{\epsilon }}$

Where C_ε and C_ν are empirical constants, and are both taken to have a value of 0.09.

In a similar way to that of the k-model, differential transport equations can be formulated for k and ε [3], thus closing the system of equations, without having to empirically define any of the turbulent quantities.

Since this method directly models the transport of all of the turbulent quantities, it is the most accurate. However, it involves solving 2 differential equations and is so the most computationally expensive. It is therefore necessary to evaluate the importance of turbulence in the problem before a model is selected. In many cases the easiest way to do this may be by trying several methods and using the simplest which gives the most accurate results. If turbulence is a relatively small feature of the flow then it is not necessary to waste computational expense by modelling it with a 2-equation model. At the same time, a zero-equation model would inadequately model a problem with large areas of turbulence, and in many cases this will lead to instabilities in the simulation.

2.3.4 Implicit LES turbulence model

In recent years a significant progress has been carried out in the development of new turbulence models based on the fact that not the entire range of scales of the flow is interesting for the majority of engineering applications. In this type of applications information contained in "the large scales" of the flow is enough to analyze magnitudes of interest as velocity, temperature ... Therefore, the idea that the global flow behaviour can be correctly approximated without the necessity to approximate the smaller scales correctly, is seen by many authors as a possible great advance in the modelling of turbulence. This fact has originated the design of turbulence models that describe the interaction of small scales with large scales. These models are commonly known as Large Eddy Simulation models (LES).

Numerous applications have shown that an extension of the FIC method allows to model low and high Reynolds number flows. In the standard large eddy simulation, a filtering process is applied to the Navier Stokes Eq. (2-6). After the filtering, a new set of equations are obtained, variables of this new set of equations are the filtered velocities or also called large scale velocities. As a consequence of the filtering process a new term called the subgrid scale tensor appears into the momentum equation. The subgrid scale tensor is defined in function of the large scale velocities and the subscale velocities too. Then it is necessary to model this term in function only of the large scale velocities. Several models exist for the subgrid scale tensor, but basically all of them propose an explicit description of the subgrid scale, an analytic expression in term of large scale velocities. In conclusion, all of these models define the subgrid scale tensor as a new nonlinear viscous term.

It is our understanding that the stabilized FIC formulation, with an adequate evaluation of the characteristic lengths introduces the same effects of a LES model without giving an explicit expression of the subgrid scale tensor. All the effects related to turbulence modelling are included into the non-linear stabilization parameters.

2.3.5 Computation of the characteristic lengths

The computation of the stabilization parameters is a crucial issue as they affect both the stability and accuracy of the numerical solution. The different procedures to compute the stabilization parameters are typically based on the study of simplified forms of the stabilized equations. Contributions to this topic are reported in (Oñate 1998 and García Espinosa 2005). Despite the relevance of the problem there still lacks a general method to compute the stabilization parameters for all the range of flow situations.

Error creating thumbnail: File missing

Figure 2 Decomposition of the velocity for characteristic length evaluation

The application of the FIC/FEM formulation to convection-diffusion problems with sharp arbitrary gradients has shown that the stabilizing FIC terms can accurately capture the high gradient zones in the vicinity of the domain edges (boundary layers) as well as the sharp gradients appearing randomly in the interior of the domain (Oñate 1998). The FIC/FEM thus reproduces the best features of both, the so called transverse (cross-wind) dissipation or shock capturing methods.

The approach proposed in this work is based on a standard decomposition of the characteristic lengths for every component of the momentum equation as

${h}_{i}=\alpha {\frac {\mathit {\boldsymbol {u}}}{\left|{\mathit {\boldsymbol {u}}}\right|}}+$$\,\beta {\frac {\nabla {u}_{i}}{\left|\nabla {u}_{i}\right|}}$

(2-40)

Being α the projection of the characteristic length component in the direction of the velocity and β the projection in the direction of the gradient of every velocity component. The decomposition defined by Eq. 40 can be demonstrated to be unique if the characteristic lengths are understood as tensors written in a system of coordinates aligned with the principal curvature directions of the solution. Obviously, Eq. (2-40) is a simplification of that approach. However, it can be easily demonstrated that the variational discrete FEM form of the resulting momentum equations are equivalent using one or the other approach.

Then, the characteristic lengths are calculated as follows

$\alpha =\left(\mathrm {coth} \,{\gamma }^{u}-{\frac {1}{{\gamma }^{u}}}\right){u}_{j}{l}_{j}\,,\quad {\gamma }^{u}=$${\frac {{u}_{j}{l}_{j}}{2\mu }}$

(2-41)

Where l_j are the maximum projections of the element on every coordinate axis (see Figure 3).

Error creating thumbnail: File missing

Figure 3 Calculation of l_i in a triangular element

and

	$\beta =\,\left(\mathrm {coth} \,{\gamma }_{ij}-{\frac {1}{{\gamma }_{ij}}}\right){l}^{\nabla {u}_{i}}\,,\quad {\gamma }_{ij}=$${u}_{i}{\frac {\nabla {u}_{i}}{\left|\nabla {u}_{i}\right|}}{\frac {{l}^{\nabla {u}_{i}}}{2\mu }}$	(2-42)
	${l}^{\nabla {u}_{i}}=\,{\frac {\nabla {u}_{i}}{\left|\nabla {u}_{i}\right|}}{l}_{i}\,,\,j=$$1,2$

As for the length parameters ${\textstyle h_{i}^{d}}$ in the mass conservation equation, the simplest assumption ${\textstyle h_{i}^{d}=h^{d}}$ has been taken.

The overall stabilization terms introduced by the FIC formulation above presented have the intrinsic capacity to ensure physically sound numerical solutions for a wide spectrum of Reynolds numbers without the need of introducing additional turbulence modelling terms.

2.4 Flow through porous media

The incompressible Navier-Stokes quations, in a given domain Ω and time interval (0,t) is given by Eq.(2-1). In the case of solids, small velocities can be considered and the term ${\textstyle \left(u\cdot \nabla \right)u}$ can be neglected. Therefore, assuming and incompressible flow (constant density) in a certain domain Ω and considering the mass conservation equation, the Navier-Stokes equations can be written as follows:

	$\rho {\frac {\partial u}{\partial t}}-\mu {\nabla }^{2}u+\nabla p=f$	$on{\mathit {\boldsymbol {\,}}}{\boldsymbol {\Omega }}x\left({\mathit {\boldsymbol {0,t}}}\right)$	(2-43)
	$\nabla \cdot u=0$	$on{\mathit {\boldsymbol {\,}}}{\boldsymbol {\Omega }}x\left({\mathit {\boldsymbol {0,t}}}\right)$

The general form of the Navier-Stokes equation is valid for the flow inside pores of a porous media, but its solution cannot be generalized to describe the flow in porous region. Therefore, the general form of the Navier-Stokes equation must be modified to describe the flow through porous media. To this aim, Darcy's law is used to describe the linear relation between the velocity ${\textstyle u}$ and the gradient of pressure ${\textstyle p}$ within the porous media. It defines the permeability resistance of the flow in a porous media as:

$\nabla p=-\mu Du$ $on{\mathit {\boldsymbol {\,}}}{\Omega }_{p}x\left({\mathit {\boldsymbol {0,t}}}\right)$

(2-44)

where ${\textstyle {\boldsymbol {\Omega }}_{\boldsymbol {p}}}$ is the porous domain, ${\textstyle D}$is the Darcy's law resistance matrix and ${\textstyle u}$ the velocity vector. In the case of considering an homogeneous porous media, ${\textstyle D}$ is a diagonal matrix with coefficients ${\textstyle 1/\alpha }$ , being ${\textstyle \alpha }$ the permeability coefficient.

As usual, the Reynolds number is defined as:

$Re={\frac {\rho UL}{\mu }}$

(2-45)

being ${\textstyle U}$and ${\textstyle L}$ a characteristic velocity and a characteristic length scale, respectively. In order to characterize the inertial effects, it is possible to define the Reynolds number associated to the pores as:

${Re}_{p}={\frac {\rho U\delta }{\mu }}$

(2-46)

where ${\textstyle \delta }$ is the characteristic pore size.

Whereas Darcy's law is reliable for values of ${\textstyle {Re}_{p}<1}$, for larger values of ${\textstyle {Re}_{p}}$ it is necessary to consider a more general model, which accounts also for the inertial effects, in the form:

$\nabla p=-\left(\mu Du+{\frac {1}{2}}\rho Cu\left|u\right|\right)$ $on{\mathit {\boldsymbol {\,}}}{\Omega }_{p}x\left({\mathit {\boldsymbol {0,t}}}\right)$

(2-47)

where ${\textstyle C}$ is the inertial resistance matrix.

The modified Navier-Stokes equation in the whole domain results from considering the two source terms associated to the resistance induced by the porous medium (linear Darcy and inertial loss term). Hence, these source terms are added to the standard fluid flow momentum equation as follows:

$\rho {\frac {\partial u}{\partial t}}-\mu {\nabla }^{2}u+\nabla p-\mu Du-{\frac {1}{2}}\rho Cu\left|u\right|=$$0$ $on{\mathit {\boldsymbol {\,}}}\Omega x\left({\mathit {\boldsymbol {0,t}}}\right)$

(2-48)

Again, considering an homogeneous porous media, ${\textstyle D}$ is a diagonal matrix with coefficients ${\textstyle 1/\alpha }$ and ${\textstyle C}$ is also a diagonal matrix. Therefore, a modified Darcy's resistance matrix should be used in Tdyn as follows:

	$\rho {\frac {\partial u}{\partial t}}-\mu {\nabla }^{2}u+\nabla p=\mu {D}^{\ast }u$	$on{\mathit {\boldsymbol {\,}}}\Omega x\left({\mathit {\boldsymbol {0,t}}}\right)$	(2-49)
	${D}^{\ast }=D+{\frac {1}{2\mu }}\rho C\left|u\right|I$

being ${\textstyle I}$ the identity matrix.

It should be noted that in laminar flows through porous media, the pressure ${\textstyle p}$ is proportional to the velocity ${\textstyle u}$, and ${\textstyle C}$ can be considered zero ${\textstyle \left({D}^{\ast }=\right.}$$\left.D\right)$ . Therefore, the Navier-Stokes momentum equation can be rewritten as:

$\rho {\frac {\partial u}{\partial t}}-\mu {\nabla }^{2}u+\nabla p=-\mu Du$ $on{\mathit {\boldsymbol {\,}}}\Omega x\left({\mathit {\boldsymbol {0,t}}}\right)$

(2-50)

Note that this momentum equation is solved by Tdyn in the case of solids, where the ${\textstyle \left(u\cdot \nabla \right)u}$ term vanishes due to small velocities. In the case that ${\textstyle \left(u\cdot \nabla \right)u}$ cannot be neglected in the modelization (i.e. high velocities), then Tdyn should solve the followimg momentum equation within a fluid instead of a solid:

$\rho \left({\frac {\partial u}{\partial t}}+\left(u\cdot \nabla \right)u\right)-$$\mu {\nabla }^{2}u+\nabla p=-\mu Du$ $on{\mathit {\boldsymbol {\,}}}\Omega x\left({\mathit {\boldsymbol {0,t}}}\right)$

(2-51)

3. Overlapping domain decomposition level set (fsurf module)

This section introduces one of the algorithms implemented in Tdyn for the analysis of free surface flows. The main innovation of this method is the application of domain decomposition concept in the statement of the problem, in order to increase accuracy in the capture of free surface as well as in the resolution of governing equations in the interface between the two fluids. Free surface capturing is based on the solution of a level set equation, while Navier Stokes equations are solved using the iterative monolithic predictor-corrector algorithm presented above, where the correction step is based on the imposition of the divergence free condition in the velocity field by means of the solution of a scalar equation for the pressure.

3.1 Governing equations

The velocity and pressure fields of two incompressible and immiscible fluids moving in the domain Ω R^d(d=2,3) can be described by the incompressible Navier Stokes equations for multiphase flows, also known as non-homogeneous incompressible Navier Stokes equations:

	${\frac {\partial \rho }{\partial t}}+{\frac {\partial }{\partial {x}_{j}}}\left(\rho u\right)=$$\,0$
	${\frac {\partial \rho {u}_{i}}{\partial t}}+{\frac {\partial }{\partial {x}_{j}}}\left(\rho {u}_{i}{u}_{j}\right)\,+$${\frac {\partial p}{\partial {x}_{i}}}-{\frac {\partial {\tau }_{ij}}{\partial {x}_{j}}}=$$\rho {f}_{i}$	(3-1)
	${\frac {\partial {u}_{i}}{\partial {x}_{j}}}=0$

Where ${\textstyle 1\leq i,j\leq d}$, ${\textstyle \rho }$ is the fluid density field, ${\textstyle {u}_{i}}$ is the ith component of the velocity field ${\textstyle u}$ in the global reference system ${\textstyle {x}_{i}}$, ${\textstyle p}$ is the pressure field and ${\textstyle \tau }$ is the viscous stress tensor defined by

${\tau }_{ij}=\mu ({\partial }_{i}{u}_{j}+\,{\partial }_{j}{u}_{i})$

(3-2)

where ${\textstyle \mu }$ is the dynamic viscosity.

Let ${\textstyle {\Omega }_{1}=\left\{x\in \Omega \,\right|x\in Fluid1\rbrace }$ be the part of the domain ${\textstyle \Omega }$ occupied by the fluid number 1 and let ${\textstyle {\Omega }_{2}=}$$\left\{x\in \Omega \,\right|x\in Fluid2\rbrace \,$ be the part of the domain ${\textstyle \Omega }$ occupied by fluid number 2. Therefore ${\textstyle {\Omega }_{1},\,{\Omega }_{2}}$ are two disjoint subdomains of ${\textstyle \Omega }$ . Therefore (see Figure 4)

$\Omega =int({\overline {{\Omega }_{1}\cap \,{\Omega }_{2}}})$

(3-3)

where ‘int’ denotes the topological interior and the over bar indicates the topological adherence of a given set, for more details see [18]. The system in Eq. (3-1) must be completed with the necessary initial and boundary conditions, as shown below.

It is usual in the literature to consider that the first equation in Eq. (3-1) is equivalent to impose a divergence free velocity field, since the density is taken as a constant. However, in the case of multiphase incompressible flows, density cannot be consider constant in ${\textstyle \Omega \,\times (0,\,T)}$. Actually, it is possible to define ${\textstyle \rho ,\,\mu }$ fields as follows:

$\rho ,\,\mu =\,\left\{{\begin{matrix}{\rho }_{1},\,{\mu }_{1}\quad \,\,x\,\in \,{\Omega }_{1}\\{\rho }_{2},\,{\mu }_{2}\quad \,\,x\,\in \,{\Omega }_{2}\end{matrix}}\right.$

(3-4)

Let ${\textstyle \Psi \,:\Omega \,\times \left(0,\,T\right)\rightarrow R}$ be a function (named Level Set function) defined as follows:

$\Psi \,\left(x,\,t\right)=\,\left\{{\begin{matrix}d\left(x,t\right)\quad \quad x\,\in \,{\Omega }_{1}\\0\quad \quad \quad \,\,x\,\in \Gamma \\-d\left(x,t\right)\quad x\,\in \,{\Omega }_{2}\end{matrix}}\right.$

(3-5)

Where ${\textstyle d(x,t)}$ is the distance to the interface between the two fluids, denoted by ${\textstyle \Gamma }$ , of the point ${\textstyle x}$ in the time instant ${\textstyle t}$. From the above definition it is trivially obtained that:

$\Gamma =\left\{\,x\,\in \Omega \,\right|\Psi \left(x,\cdot \right)=$$\,0\,\rbrace$

(3-6)

Since the level set 0 identifies the free surface between the two fluids, the following relations can be obtained:

$n\left(x,t\right)=\nabla \,\Psi {\left.\right|}_{(x,t)}\,;\kappa \left(x,t\right)=$$\nabla \cdot (n\left(x,t\right))$

(3-7)

Where ${\textstyle n}$ is the normal vector to the interface ${\textstyle \Gamma }$ , oriented from fluid 1 to fluid 2 and ${\textstyle \kappa }$ is the curvature of the free surface. In order to obtain the above relations it has been assumed that function ${\textstyle \Psi }$ fulfils:

$\left|\nabla \Psi \right|=1\quad \quad \quad \forall \,\left(x,t\right)\in \Omega \,\times (0,T)$

(3-8)

Therefore, it is possible to redefine density and viscosity as follows:

$\rho ,\,\mu =\,\left\{{\begin{matrix}{\rho }_{1},\,{\mu }_{1}\,\,\Psi >0\\{\rho }_{2},\,{\mu }_{2}\,\,\Psi <0\end{matrix}}\right.$

(3-9)

Let us write the density fields in terms of the level set function ${\textstyle \Psi }$ as

$\rho \left(x,t\right)=\rho \left(\Psi \left(x,t\right)\right)\,\,\forall \,\left(x,t\right)\in \Omega \times (0,\,T)$

(3-10)

Then, density derivatives can be written as

${\frac {\partial \rho }{\partial t}}={\frac {\partial \rho }{\partial \Psi }}{\frac {\partial \Psi }{\partial t}},\quad {\frac {\partial \rho }{\partial {x}_{i}}}=$${\frac {\partial \rho }{\partial \Psi }}{\frac {\partial \Psi }{\partial {x}_{i}}}$

(3-11)

Inserting the above relation into the first line of Eq. (3-1) gives

${\frac {\partial \rho }{\partial t}}+{\frac {\partial }{\partial {x}_{i}}}\,\left(\rho {u}_{i}\right)=$$\,{\frac {\partial \rho }{\partial t}}+\,{u}_{i}{\frac {\partial \rho }{\partial {x}_{i}}}=$$\,{\frac {\partial \rho }{\partial \Psi }}{\frac {\partial \Psi }{\partial t}}+\,{\frac {\partial \rho }{\partial \Psi }}{u}_{i}{\frac {\partial \Psi }{\partial {x}_{i}}}=$$\,{\frac {\partial \rho }{\partial \Psi }}\left[{\frac {\partial \Psi }{\partial t}}+\right.$$\left.\,{u}_{i}{\frac {\partial \Psi }{\partial {x}_{i}}}\right]=0$

(3-12)

which gives as a result that the multiphase Navier Stokes problem is equivalent to solve the following system of equations:

	${\frac {\partial \rho {u}_{i}}{\partial t}}+{\frac {\partial }{\partial {x}_{j}}}\left(\rho {u}_{i}{u}_{j}\right)+$${\frac {\partial p}{\partial {x}_{i}}}-{\frac {\partial {\tau }_{ij}}{\partial {x}_{j}}}=$$\rho {f}_{i}$	(3-13)
	${\frac {\partial {u}_{i}}{\partial {x}_{i}}}=0$

coupled with the equation

${\frac {\partial \Psi }{\partial t}}+\,{u}_{i}{\frac {\partial \Psi }{\partial {x}_{i}}}=$$0$

(3-14)

Eq. (3-14) defines the transport of the level set function due to the velocity field obtained by solving Eq. (3-13).

As a conclusion, the free surface capturing problem can be described by the above equations. In this formulation, the interface between the two fluids is defined by the level 0 of ${\textstyle \Psi }$ .

Denoting prescribed values by an over-bar, the boundary conditions to be considered for the above presented problem are

	$u=\,{\overline {u}}\quad \,\,on\,{\Gamma }_{u}$
	$p=\,{\overline {p}},\,\,{n}_{j}{\tau }_{ij}={\overline {t}}_{i}\quad \quad on\,{\Gamma }_{p}$	(3-15)
	${u}_{j}{n}_{j}=\,{\overline {u}}_{n}\,,\quad \left.{\begin{matrix}{n}_{j}{\tau }_{ij}{g}_{i}=\,{\overline {t}}_{1}\\{n}_{j}{\tau }_{ij}{s}_{i}=\,{\overline {t}}_{2}\end{matrix}}\right\}\quad on\,{\Gamma }_{\tau }$

Where the boundary ${\textstyle \partial \Omega }$ of the domain ${\textstyle \Omega }$ has been split in three disjoint sets: ${\textstyle {\Gamma }_{u},\,{\Gamma }_{p}}$ where the Dirichlet and Neumann boundary conditions are imposed and ${\textstyle {\Gamma }_{\tau }}$ where the Robin conditions for the velocity are set. In the above, vectors ${\textstyle g,\,s}$ span the space tangent to ${\textstyle {\Gamma }_{\tau }}$. In a similar way, the boundary conditions for Eq. (3-14) are defined

$\Psi =\,{\overline {\Psi }}\quad on\,{\Gamma }_{u}$

(3-16)

Finally, the initial conditions of the problem are

$u={u}_{0}\quad on\,\Omega ,\quad \Psi ={\Psi }_{0}\quad on\,\Omega$

(3-17)

Where ${\textstyle {\Gamma }_{0}=\left\{x\in \Omega \,\right|{\Psi }_{0}\,\left(x\right)=}$$0\rbrace$ defines the initial position of the free surface between the two fluids.

3.2 How to calculate a signed distance to the interface

Based on flow properties as density or viscosity, we can calculate an initial guess for ${\textstyle \Psi }$ as follows

$\Psi \left(x,t\right)=\,\left\{{\begin{matrix}1\quad \,\,if\,\rho \left(x,t\right)=\,{\rho }_{1}\\0\quad \,if\,x\,\in \Gamma (t)\\-1\quad if\,\rho \left(x,t\right)=\,{\rho }_{2}\end{matrix}}\right.$

(3-18)

Clearly, the level set function obtained from Eq. (3-18) is not a signed distance to ${\textstyle \Gamma }$ . We need to reinitiate ${\textstyle \Psi }$ in order to guarantee that ${\textstyle \Psi }$ fulfils the definition. We use a technique based on the following property to reinitiate the level set function as a signed distance to ${\textstyle \Gamma }$ . If ${\textstyle \Psi }$ is a signed distance function to ${\textstyle \Gamma }$ then,

$\left\|\nabla \Psi \right\|=1\quad \forall \,\left(x,t\right)\in \Omega \,\times (0,T]$

(3-19)

where ${\textstyle \left\|\cdot \right\|}$ denotes the Euclidian norm in ${\textstyle {R}^{d}}$.

We use Eq. (3-19) to recalculate ${\textstyle \Psi }$ . In case that Eq. (3-19) is not satisfied, then

$1-\,\left\|\nabla \Psi \right\|=r$

(3-20)

The residual r in Eq. (3-20) can be understood as the time derivate of ${\textstyle \Psi }$ with respect to a pseudo-time ${\textstyle \tau }$ , i.e.

$r\,:={\partial }_{\tau }\Psi$

(3-21)

Then, the stationary solutions of Eq. (3-21) (i.e. when ${\textstyle {\partial }_{\tau }\Psi =}$$0$) satisfy Eq. (3-19). In summary, for a given time ${\textstyle t\,\in [0,\,T]}$, we calculate ${\textstyle \Psi }$ as the stationary solution of the following hyperbolic problem.

	${\partial }_{\tau }\Psi +\left\|\nabla \Psi \right\|=1\quad on\,\Omega \,$
	$\Psi \left(\tau =0\right)=\Psi (\cdot ,t)$	(3-22)
	$\Psi \left(x,\tau \right)=\,0\quad \forall \,x\,\in \,\Gamma (t)$

The problem stated in Eq. (3-22) can be rewritten in a more convenient manner as

${\partial }_{\tau }\Psi +\left({\boldsymbol {w}}\cdot \nabla \right)\Psi =sign\left(\Psi \right);\quad {\boldsymbol {w}}=$$sign\left(\Psi \right){\frac {\nabla \Psi }{\left\|\nabla \Psi \right\|}}$

(3-23)

Since the level set values identify the free surface between the two fluids, the following relations can be obtained:

${\mathit {\boldsymbol {n}}}\left({\mathit {\boldsymbol {x}}},t\right)=\nabla \Psi {\left.\right|}_{x,t}\,;\,\kappa \left({\mathit {\boldsymbol {x}}},t\right)=$$\nabla \cdot ({\mathit {\boldsymbol {n}}}\left({\mathit {\boldsymbol {x}}},t\right))$

(3-24)

where ${\textstyle {\mathit {\boldsymbol {n}}}}$ is the normal vector to the interface ${\textstyle \Gamma }$ , oriented from fluid 2 to fluid 1 and ${\textstyle \kappa }$ is the local curvature of the interface.

From equation Eq. (3-23) it can be easily shown that the reinitialization process begins from the interface and it propagates to the whole domain. This property is attractive from the computational point of view, since we are only interested in accurate values of ${\textstyle \Psi }$ in the region close to the interface. Then we only need to reach the steady state solution of Eq. (3-23) in a narrow band close to the interface as proposed in [13]. The reinitialization algorithm used in this work has been optimized using the concept of nodal levels described next.

Given a level set distribution ${\textstyle \Psi }$ , for every time step ${\textstyle t}$ one can assign a level ${\textstyle {l}_{i}}$ to the node ${\textstyle {\mathit {\boldsymbol {x}}}_{i}}$ based on the following rules:

• A nodal point i of coordinates ${\textstyle {\mathit {\boldsymbol {x}}}_{i}}$ belongs to the level ${\textstyle {l}_{i}=}$$1$, if ${\textstyle \Psi \left({\mathit {\boldsymbol {x}}}_{i},t\right)\geq 0}$, and it is connected to at least one node j with ${\textstyle \Psi \left({\mathit {\boldsymbol {x}}}_{j},t\right)<0}$

• A nodal point i of coordinates ${\textstyle {\mathit {\boldsymbol {x}}}_{i}}$ belongs to the level ${\textstyle {l}_{i}=\,-}$$1$, if ${\textstyle \Psi \left({\mathit {\boldsymbol {x}}}_{i},t\right)<0}$, and it is connected to at least one node j with ${\textstyle \Psi \left({\mathit {\boldsymbol {x}}}_{j},t\right)\geq 0}$

• A nodal point i of coordinates ${\textstyle {\mathit {\boldsymbol {x}}}_{i}}$ belongs to the level ${\textstyle {l}_{i}>1}$, if ${\textstyle \Psi \left({\mathit {\boldsymbol {x}}}_{i},t\right)>0}$, and it is connected to at least one node of level ${\textstyle {l}_{i-1}}$

• A nodal point i of coordinates ${\textstyle {\mathit {\boldsymbol {x}}}_{i}}$ belongs to the level ${\textstyle {l}_{i}<-}$$1$, if ${\textstyle \Psi \left({\mathit {\boldsymbol {x}}}_{i},t\right)<0}$, and it is connected to at least one node of level ${\textstyle {l}_{i+1}}$

Once the levels have been assigned to the mesh points, the reinitialization is performed in a narrow band of n levels around the interface. Eq. (3-23) is solved using a forward Euler scheme until the steady state in the n levels is achieved. k iterations with k > n are performed in such a way that for every iteration i = 1,…,k. The iterations are carried out only for those nodes j of the mesh with |l_j| ≤ i.

Due to the evolution of the interface defined by Eq. (3-14), the level set function does not fulfil Eq. (3-19) after some time. It is therefore recommended to apply the reinitialization scheme described above at every time step.

3.3 Interfacial boundary conditions

At a moving interface the jump condition applies. For immiscible fluids we can write:

${\begin{matrix}\left[{\mathit {\boldsymbol {n}}}\cdot \left({\mathit {\boldsymbol {\sigma }}}\cdot {\mathit {\boldsymbol {n}}}\right)\right]=\gamma \kappa \quad \quad on\,\Gamma \\\left[{\mathit {\boldsymbol {s}}}\cdot \left({\mathit {\boldsymbol {\sigma }}}\cdot {\mathit {\boldsymbol {n}}}\right)\right]=0\quad \quad \,\,on\,\Gamma \end{matrix}}$

(3-25)

where ${\textstyle \gamma \,}$ is the coefficient of surface tension (a constant of the problem), ${\textstyle {\mathit {\boldsymbol {s}}}}$ is any unit vector tangent to the interface and ${\textstyle [\cdot ]}$ defines the jump across the interface.

$\left[A\right]={A}_{1}-\,{A}_{2}$

(3-26)

Taking into account that ${\textstyle {\mathit {\boldsymbol {\sigma }}}_{i}=\,-}$$p{\boldsymbol {I}}+\,{\mu }_{i}\left(\nabla {\mathit {\boldsymbol {u}}}+\nabla {\mathit {\boldsymbol {u}}}^{T}\right)\quad i=$$1,2$ and using taking into account that ${\textstyle \cdot {\mathit {\boldsymbol {n}}}=}$$0$ , the second equation in Eq. 67 yields

${\boldsymbol {s}}\cdot \left[\left({\mu }_{1}-\,{\mu }_{2}\right)\left(\nabla {\mathit {\boldsymbol {u}}}+\right.\right.$$\left.\left.\nabla {\mathit {\boldsymbol {u}}}^{T}\right)\cdot {\mathit {\boldsymbol {n}}}\right]=$$\,0$

(3-27)

From Eq. (3-27) we can conclude that ${\textstyle \left({\mu }_{1}-\right.}$$\left.\,{\mu }_{2}\right)\left(\nabla {\mathit {\boldsymbol {u}}}+\nabla {\mathit {\boldsymbol {u}}}^{T}\right)\cdot {\mathit {\boldsymbol {n}}}\in {\left\langle s\right\rangle }^{\perp }$, where ${\textstyle {\left\langle {\mathit {\boldsymbol {v}}}\right\rangle }^{\perp }}$, denotes the orthogonal subspace to ${\textstyle {\mathit {\boldsymbol {v}}}}$. Then it is clear that

$\left({\mu }_{1}-\,{\mu }_{2}\right)\left(\nabla {\mathit {\boldsymbol {u}}}+\nabla {\mathit {\boldsymbol {u}}}^{T}\right)\cdot {\mathit {\boldsymbol {n}}}=$$\lambda {\mathit {\boldsymbol {n}}},\quad \,\lambda \in Spec({\mathit {\boldsymbol {\tau }}}_{1}-$$\,{\mathit {\boldsymbol {\tau }}}_{2})$

(3-28)

Integrating over any closed surface ${\textstyle \partial M}$ containing ${\textstyle \Gamma }$ in both sides of Eq. (3-28) yields

$\int _{\partial M}^{\,}\left({\mu }_{1}-\,{\mu }_{2}\right){\mathit {\boldsymbol {n}}}\cdot (\nabla {\mathit {\boldsymbol {u}}}+\nabla {\mathit {\boldsymbol {u}}}^{T})=$$\,\int _{\partial M}^{\,}\lambda {\mathit {\boldsymbol {n}}}$

(3-29)

Applying the Gauss theorem to the right hand side of Eq. (3-29) we can deduce that ${\textstyle \lambda =}$$0$ or equivalently

${\mathit {\boldsymbol {n}}}\cdot \left(\nabla {\mathit {\boldsymbol {u}}}+\nabla {\mathit {\boldsymbol {u}}}^{T}\right)\quad \,on\,\Gamma$

(3-30)

Introducing Eq. (3-30) into the first condition in Eq. (3-25) the following equivalent condition is obtained:

${p}_{1}{\mathit {\boldsymbol {n}}}={p}_{2}{\mathit {\boldsymbol {n}}}+\gamma \kappa {\mathit {\boldsymbol {n}}}$

(3-31)

This form of the jump condition is very useful in the development of the method presented in the following sections.

In the classical level set formulation for multiphase flow problems [15] the jump condition Eq. (3-25) can be introduced into the momentum equation by adding an additional source term, i.e.

${\begin{matrix}{\partial }_{t}\left(\rho {\mathit {\boldsymbol {u}}}\right)+\,\nabla \left(\rho {\mathit {\boldsymbol {u}}}\otimes {\mathit {\boldsymbol {u}}}\right)-\,\nabla {\mathit {\boldsymbol {\sigma }}}=\rho {\tilde {\mathit {\boldsymbol {f}}}}\\{\tilde {\mathit {\boldsymbol {f\,}}}}{\mathit {\boldsymbol {=f+}}}\gamma \kappa {\mathit {\boldsymbol {n}}}\delta (\Gamma )\end{matrix}}$

(3-32)

Where ${\textstyle \delta \left(\Gamma \right)}$ is a Dirac’s delta function on ${\textstyle \Gamma }$ . Since the properties of the fluids change discontinuously across the interface, the direct solution of Eq. (3-32) has important numerical drawbacks: inaccurate resolution of the pressure at the interface, fluid mass losses, inaccurate interface location, etc. A common approach to overcome these difficulties is to smooth the discontinuous transition of the flow properties and the Dirac’s delta functions across the interface using of a smoothed Heaviside function. This function is usually expressed as:

${H}_{\epsilon }\left(\Psi \right)=\,\left\{{\begin{matrix}0&\Psi <\,-\epsilon \\{\frac {1}{2}}\left[1+{\frac {\Psi }{\epsilon }}+{\frac {1}{\pi }}\mathrm {sin} \,\left({\frac {\pi \Psi }{\epsilon }}\right)\right]&\left|\Psi \right|\leq \epsilon \\1&\Psi >\,\epsilon \end{matrix}}\right.$

(3-33)

where ${\textstyle 2\epsilon }$ is the thickness of the transition area. Then a property ${\textstyle \eta }$ is approximated as

$\eta \left(\phi \right)=\,{\eta }_{1}+\,\left({\eta }_{1}-\,{\eta }_{2}\right){H}_{\epsilon }(\phi )$

(3-34)

and a Dirac’s delta function as

$\delta \left(\Gamma \right)=\,{d}_{\Psi }{H}_{\epsilon }$

(3-35)

This approach also presents several deficiencies: the accuracy of the results depends on the ad hoc thickness of the transition area ${\textstyle 2\epsilon }$ , the jump condition is not properly captured, the interface propagation velocity is not correctly calculated, etc.

Other schemes, such as the Ghost Cell method or lagrangian methods as the Particle FEM Method, attempt to overcome the problem of properties smoothing at the interface region. However, these methods have several drawbacks, such us limitations of conservative and convergence properties in some cases or an increase of the computational effort. As a conclusion, a method able to deal with a discontinuous transition of properties across the interface, the jump interface condition and computationally affordable for large models is required.

The overlapping domain decomposition method described in the following sections has been developed to overcome these difficulties.

3.4 Stabilized FIC-FEM formulation

The stabilized FIC form of the governing differential equations Eq. (3-13) and Eq. (3-14) is written as

${r}_{m}{\overline {-{\frac {1}{2}}\left({h}_{m}\cdot \nabla \right){r}_{m}}}=$$0,\quad {r}_{d}{\overline {-{\frac {1}{2}}\left({h}_{d}\cdot \nabla \right){r}_{d}}}=$$0,\quad {r}_{\Psi }{\overline {-{\frac {1}{2}}\left({h}_{\Psi }\cdot \nabla \right)}}=$$0$

(3-36)

The boundary conditions for the stabilized FIC problem are written as

	${\mathit {\boldsymbol {u}}}={\overline {\mathit {\boldsymbol {u}}}}\quad \,on\,{\Gamma }_{D}$
	${\mathit {\boldsymbol {n}}}\cdot {\mathit {\boldsymbol {\sigma }}}-{\overline {{\frac {1}{2}}\,{\mathit {\boldsymbol {n}}}\cdot {\mathit {\boldsymbol {h}}}_{m}\cdot {\mathit {\boldsymbol {r}}}_{m}}}=$${\overline {\mathit {\boldsymbol {t}}}}\quad \,\,on\,{\Gamma }_{N}$
	${\mathit {\boldsymbol {n}}}\cdot {\mathit {\boldsymbol {u}}}={\overline {\mathit {\boldsymbol {u}}}}_{n},\quad \,{\mathit {\boldsymbol {n}}}\cdot {\mathit {\boldsymbol {\sigma }}}\cdot {\mathit {\boldsymbol {g}}}{\overline {-{\frac {1}{2}}{\mathit {\boldsymbol {n}}}\cdot {\mathit {\boldsymbol {h}}}_{m}\cdot {\mathit {\boldsymbol {r}}}_{m}\cdot {\mathit {\boldsymbol {g}}}}}=$${\overline {\mathit {\boldsymbol {t}}}}_{1}\quad \,on\,{\Gamma }_{M}$	(3-37)
	${\mathit {\boldsymbol {n}}}\cdot {\mathit {\boldsymbol {\sigma }}}\cdot {\mathit {\boldsymbol {s}}}{\overline {-{\frac {1}{2}}{\mathit {\boldsymbol {n}}}\cdot {\mathit {\boldsymbol {h}}}_{m}\cdot {\mathit {\boldsymbol {r}}}_{m}\cdot {\mathit {\boldsymbol {s}}}}}=$${\overline {\mathit {\boldsymbol {t}}}}_{2}$
	$\Psi =\,{\overline {\Psi }}\quad \,on\,{\Gamma }_{inlet}$

The underlined terms in Eq. (3-36) and Eq. (3-37) introduce the necessary stabilization for the numerical solution of the Navier Stokes problem [2], [4], [6].

Note that terms ${\textstyle {\mathit {\boldsymbol {r}}}_{m}}$, ${\textstyle {\mathit {\boldsymbol {r}}}_{d}}$ and ${\textstyle {\mathit {\boldsymbol {r}}}_{\Psi }}$ denote the residuals of Eq. (3-13) and Eq. (3-14), i.e.

	${\mathit {\boldsymbol {r}}}_{m}\,:=\,{\partial }_{t}\left(\rho {\mathit {\boldsymbol {u}}}\right)+$$\,\nabla \cdot \left(\rho {\mathit {\boldsymbol {u}}}\otimes {\mathit {\boldsymbol {u}}}\right)+$$\,\nabla p-\nabla \cdot {\mathit {\boldsymbol {\sigma }}}-\rho {\mathit {\boldsymbol {f}}}$
	${\mathit {\boldsymbol {r}}}_{d}=\nabla \cdot {\mathit {\boldsymbol {u}}}$	(3-38)
	${\mathit {\boldsymbol {r}}}_{\Psi }\,:=\,{\partial }_{t}\Psi +\left({\mathit {\boldsymbol {u}}}\cdot \nabla \right)\Psi$

The characteristic lengths vectors ${\textstyle {\boldsymbol {h}}_{m},{\boldsymbol {h}}_{d}}$ and ${\textstyle {\boldsymbol {h}}_{\psi }}$ contain the dimensions of the finite space domain where the balance laws of mechanics are enforced. At the discrete level these dimensions are of the order of the element or grid dimension used for the computation. Details on obtaining the FIC stabilized equations and the recommendation for computing the characteristic length vectors can be found in [4]-[6].

3.5 Overlapping domain decomposition method

To overcome the drawbacks related to the discontinuity of fluid properties across the interface and to impose the interfacial boundary condition we split the domain ${\textstyle \Omega }$ into two overlapping subdomains. Based on this subdivision of the domain, we apply a Dirichlet-Neumann overlapping domain decomposition technique with appropriate boundary conditions in order to satisfying the interfacial condition.

Let ${\textstyle K=\bigcup _{e=1}^{\#K}\left({N}^{e},{K}^{e},{\Theta }^{e}\right)\,}$ be a finite element partition of domain the ${\textstyle \Omega }$ , where ${\textstyle N^{e}}$ are element nodes, ${\textstyle K^{e}}$ is the element spatial domain, ${\textstyle {\Theta }^{e}}$ are element shape functions and ${\textstyle \mathrm {K} }$ is the total number of elements in the partition.

We assume that ${\textstyle \mathrm {K} }$ satisfies the following approximation property:

$dist\left(\partial \Omega ,\,\partial \left(\bigcup _{e=1}^{\#K}{K}^{e}\right)\right)\leq \mathrm {} \,\lbrace diam\left({K}^{e}\right)\rbrace \,$

(3-39)

for a fixed instant ${\textstyle t,\,t\,\in (0,T]}$. Let us consider a domain decomposition of domain ${\textstyle \Omega }$ into three disjoint sub domains ${\textstyle {\Omega }_{3}\left(t\right),{\mbox{ }}{\Omega }_{4}\left(t\right)}$ and ${\textstyle {\Omega }_{5}\left(t\right)}$ in such a way that ${\Omega }_{\mbox{3}}\left(t\right)={\underset {e}{\cup }}K_{3}^{e}$ , ${\Omega }_{\mbox{5}}={\underset {e}{\cup }}K_{5}^{e}$ , where ${\textstyle K_{e}^{3}}$ are the elements of the finite element partition ${\textstyle \mathrm {K} }$ , such as ${\textstyle \forall {\boldsymbol {x}}\in K_{3}^{e}\vert \psi \left(x,t\right)>0}$ and ${\textstyle K_{e}^{5}}$ are the elements of the finite element partition such as ${\textstyle \forall {\boldsymbol {x}}\in K_{5}^{e}\vert \psi \left(x,t\right)<0}$ . The geometrical domain decomposition is completed with

${\Omega }_{4}\left(t\right)=\Omega \backslash \left({\Omega }_{3}\left(t\right)\cup {\Omega }_{5}\left(t\right)\right)$

(3-40)

From this partition we define the two overlapping domains:

${\overline {\Omega }}_{1}\left(t\right):=int\left({\overline {{\Omega }_{3}\left(t\right)\cup {\Omega }_{4}\,\left(t\right)}}\right),\quad \,{\overline {\Omega }}_{2}\left(t\right):=int\left({\overline {{\Omega }_{4}\left(t\right)\cup {\Omega }_{5}\,\left(t\right)}}\right)$

(3-41)

Some comments about the ${\textstyle {\tilde {\Omega }}_{1}}$ - ${\textstyle {\tilde {\Omega }}_{2}}$ partition are given next.

In order to simplify the notation we will omit the time dependency of domains in the rest of the section.

Let ${\textstyle \partial {\tilde {\Omega }}_{i},{\mbox{ }}i=1,2}$ be the boundary of ${\textstyle {\tilde {\Omega }}_{i}}$ , then:

$\partial {\tilde {\Omega }}_{i}\,\cap \partial \Omega =\cup \cup$

(3-42)

It has to be noted that ${\textstyle \partial {\tilde {\Omega }}_{i}}$ contains an additional term ${\textstyle {\tilde {\Gamma }}_{i}}$ that is not included into ${\textstyle \partial \Omega }$ . This term comes from the presence of an interface inside ${\textstyle \Omega }$ . Since we are using a capturing technique, the position of the interface will not lay on mesh nodes, as usually happens in a tracking technique. Therefore some elements can be intersected by the interface. Then ${\textstyle {\tilde {\Gamma }}_{i}}$ is defined as follows:

${\tilde {\Gamma }}_{i}=\partial {\Omega }_{4}\cap {\tilde {\Omega }}_{i+{\left(-1\right)}^{i-1}}$

(3-43)

Note that ${\textstyle {\tilde {\Gamma }}_{i}}$ is not coincident with ${\textstyle \Gamma }$ , but the following condition is satisfied:

$dist\left(\Gamma ,{\tilde {\Gamma }}_{i}\right)\leq \mathrm {max} \,\lbrace diam({K}^{e})\vert {K}^{e}\subset {\Omega }_{4}$

(3-44)

Summarizing, we have

$\partial {\tilde {\Omega }}_{i}={\Gamma }_{iD}\cup {\Gamma }_{iN}\cup {\Gamma }_{iM}\cup {\tilde {\Gamma }}_{i}$

Figure 5 is a sketch of the domain partition described above.

Error creating thumbnail: File missing Figure 5 Geometrical domain decomposition

We have two restrictions for the definition of the boundary condition on ${\textstyle {\tilde {\Gamma }}_{i}}$ : fluid velocities must be compatible at the interface and the jump boundary condition Eq. (3-25) must be satisfied on ${\textstyle \Gamma }$ . The domain decomposition technique chosen is of Dirichlet-Neumann type and it allows us to fulfil both restrictions.

Let us introduce a standard notation. We denote by ${\textstyle x_{i}}$ a variable related to ${\textstyle {\tilde {\Omega }}_{i}}$ . For example u₁ is the velocity field solution of the problem in Eq. (3-37), where the domain ${\textstyle \Omega }$ has been replaced by ${\textstyle {\tilde {\Omega }}_{1}}$ .

We apply the Dirichlet conditions on ${\textstyle {\tilde {\Gamma }}_{1}}$ and the Neumann conditions on ${\textstyle {\tilde {\Gamma }}_{2}}$ . For the boundary ${\textstyle {\tilde {\Gamma }}_{1}}$ , we make use of the compatibility of velocities at the interface:

${\boldsymbol {u}}_{1}={\boldsymbol {u}}_{2}\quad \,on\,{\tilde {\Gamma }}_{1}$

(3-45)

For the boundary ${\textstyle {\tilde {\Gamma }}_{2}}$ we make use of the jump boundary condition. We are looking for a traction vector ${\textstyle {\boldsymbol {\tilde {t}}}}$ such that:

${\boldsymbol {n}}_{2}\cdot {\boldsymbol {\sigma }}_{2}={\tilde {\boldsymbol {t}}}{\boldsymbol {\,}}\quad on\,{\tilde {\Gamma }}_{2}$

(3-46)

and ${\textstyle {\boldsymbol {\tilde {t}}}}$ must guarantee that the jump boundary condition Eq. (3-25) holds, that is to say:

${p}_{1}{\mathit {\boldsymbol {n}}}={p}_{2}{\mathit {\boldsymbol {n}}}+\gamma \kappa {\mathit {\boldsymbol {n}}}\quad \,on\,\Gamma$

(3-47)

Let us write the variational problem considering the domain decomposition above described.

Domain 1

	${\left({\rho }_{1}{\partial }_{t}{\boldsymbol {u}}_{1},{\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{1}}+$${\left\langle {\rho }_{1}\left({\boldsymbol {u}}_{1}\cdot \nabla \right){\boldsymbol {u}}_{1},{\boldsymbol {v}}\right\rangle }_{{\tilde {\Omega }}_{1}}+$$\,{\left({\mathit {\boldsymbol {\sigma }}}_{1},\nabla {\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{1}}+$${\frac {1}{2}}{\left({\boldsymbol {r}}_{1m},\left({\boldsymbol {h}}_{1m}\cdot \nabla \right){\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{1}}+$$\,{\left({\overline {\boldsymbol {t}}},{\boldsymbol {v}}\right)}_{{\Gamma }_{1N}}+$$\,{\left({\overline {\boldsymbol {t}}}_{\boldsymbol {1}}{\boldsymbol {g+}}{\overline {\boldsymbol {t}}}_{\boldsymbol {2}}{\boldsymbol {s,v}}\right)}_{{\Gamma }_{1M}}=$${\left\langle {\rho }_{1}{\boldsymbol {f}},{\boldsymbol {v}}\right\rangle }_{{\tilde {\Omega }}_{1}}\,$
	${\left(q,\nabla \cdot {\boldsymbol {u}}_{1}\right)}_{{\tilde {\Omega }}_{1}}+{\frac {1}{2}}{\left({r}_{1d},\left({\boldsymbol {h}}_{1d}\cdot \nabla \right)q\right)}_{{\tilde {\Omega }}_{1}}=$$0$	(3-48)
	${\boldsymbol {u}}_{1}={\boldsymbol {u}}_{2}\quad \,on\,{\tilde {\Gamma }}_{1}$

Domain 2

	${\left({\rho }_{2}{\partial }_{t}{\boldsymbol {u}}_{2},{\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{2}}+$${\left\langle {\rho }_{2}\left({\boldsymbol {u}}_{2}\cdot \nabla \right){\boldsymbol {u}}_{2},{\boldsymbol {v}}\right\rangle }_{{\tilde {\Omega }}_{2}}+$$\,{\left({\mathit {\boldsymbol {\sigma }}}_{2},\nabla {\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{2}}+$${\frac {1}{2}}{\left({\boldsymbol {r}}_{2m},\left({\boldsymbol {h}}_{2m}\cdot \nabla \right){\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{2}}+$$\,{\left({\overline {\boldsymbol {t}}},{\boldsymbol {v}}\right)}_{{\Gamma }_{2N}}+$$\,{\left({\overline {\boldsymbol {t}}}_{\boldsymbol {1}}{\boldsymbol {g+}}{\overline {\boldsymbol {t}}}_{\boldsymbol {2}}{\boldsymbol {s,v}}\right)}_{{\Gamma }_{2M}}=$${\left\langle {\rho }_{2}{\boldsymbol {f}},{\boldsymbol {v}}\right\rangle }_{{\tilde {\Omega }}_{2}}$	(3-49)
	${\left(q,\nabla \cdot {\boldsymbol {u}}_{2}\right)}_{{\tilde {\Omega }}_{2}}+{\frac {1}{2}}{\left({r}_{2d},\left({\boldsymbol {h}}_{2d}\cdot \nabla \right)q\right)}_{{\tilde {\Omega }}_{2}}=$$0$

An iterative strategy between the two domains is used to reach a converged global solution in the whole domain ${\textstyle \Omega }$ . It is expected that this global solution will satisfy both restrictions presented above. For this reason we define the following stop criteria:

${\begin{matrix}\left\|{u}_{1}-\,{u}_{2}\right\|\leq {tol}_{u}\quad \,on\,\Omega \\\left|{\boldsymbol {n}}\cdot \left({\mathit {\boldsymbol {\sigma }}}_{1}-{\mathit {\boldsymbol {\sigma }}}_{2}\right)\cdot {\boldsymbol {n}}-\gamma \kappa \right|\leq {tol}_{\sigma }\quad \,on\,\Gamma \end{matrix}}$

(3-50)

The global velocity solution obtained from Eq. (3-49) is used as the convective velocity for the transport of the level set function, i.e.

${\left({\partial }_{t}\Psi ,\zeta \right)}_{\Omega }+{\left\langle \left({\boldsymbol {u}}\cdot \nabla \right)\Psi ,\zeta \right\rangle }_{\Omega }+$${\frac {1}{2}}{\left({r}_{\Psi },\left({\boldsymbol {h}}_{\Psi }\cdot \nabla \right)\zeta \right)}_{\Omega }=$$0\,$

(3-51)

Unfortunately there is not theoretical result that can ensure the convergence of this method. However our experience in applying this method to a wide range of problems involving moving interfaces has showed that the method is stable and robust.

An expected property of the method is its dependency with the mesh size around the interface. This issue is directly related with how good is ${\textstyle \Gamma }$ approximated by ${\textstyle {\tilde {\Gamma }}_{i}}$ , as it has been point out in property Eq. (3-44). A fine mesh close to the interface reduces significantly the amount of iterations needed to reach a converged global solution, as well as the accuracy of the velocities and the pressure at the interface.

This methodology can be viewed as a combination of domain decomposition and level set techniques. This justifies the name given to the new method: ODDLS (by Overlapping Domain Decomposition Level Set).

3.6 Discretization by the finite element method (FEM)

In this section we present the final discrete system of equations associated to problem Eq. (3-48) and Eq. (3-49) using the FEM method. Let us consider a uniform partition of the time interval of analysis ${\textstyle \left[0,T\right]}$ , with time step size ${\textstyle \delta t}$ . We will denote by a superscript the time step at which the algorithmic solution is computed. We assume ${\textstyle {\boldsymbol {u}}_{1}^{n},p_{1}^{n},{\boldsymbol {u}}_{2}^{n}}$ and ${\textstyle p_{2}^{n}}$ are known. If ${\textstyle \theta \in \left[0,1\right]}$ , the trapezoidal rule applied to Eq. (3-48) and Eq. (3-49) consists of finding ${\textstyle {\boldsymbol {u}}_{1}^{n+1},p_{1}^{n+1},{\boldsymbol {u}}_{2}^{n+1},p_{2}^{n+1}}$ as the solution of the problem

	${\left({\rho }_{i}^{n}{\delta }_{t}{\boldsymbol {u}}_{i}^{n},{\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{i}}+$${\left\langle {\rho }_{i}^{n}\left({\boldsymbol {u}}_{i}^{n+\theta }\cdot \nabla \right){\boldsymbol {u}}_{i}^{n+\theta },{\boldsymbol {v}}\right\rangle }_{{\tilde {\Omega }}_{i}}+$$\,{\left({\mathit {\boldsymbol {\sigma }}}_{i}^{n+\theta },\nabla {\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{i}}+$${\frac {1}{2}}{\left({\boldsymbol {r}}_{im}^{n+\theta },\left({\boldsymbol {h}}_{im}^{n+\theta }\cdot \nabla \right){\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{i}}+$$\,{\left({\tilde {\boldsymbol {t}}}{\boldsymbol {,v}}\right)}_{{\Gamma }_{iN}}+\,{\left({\tilde {\boldsymbol {t}}}_{\boldsymbol {1}}{\boldsymbol {g+}}{\tilde {\boldsymbol {t}}}_{\boldsymbol {2}}{\boldsymbol {s,v}}\right)}_{{\Gamma }_{iM}}+$$\,\left(i-1\right){\left({\tilde {\boldsymbol {t}}}{\boldsymbol {,v}}\right)}_{{\tilde {\Gamma }}_{i}}=$$\,{\left\langle {\rho }_{i}^{n}{\boldsymbol {f}},{\boldsymbol {v}}\right\rangle }_{{\tilde {\Omega }}_{i}}$	(3-52)
	${\left(q,\nabla \cdot {\boldsymbol {u}}_{i}^{n+\theta }\right)}_{{\tilde {\Omega }}_{i}}+$${\frac {1}{2}}{\left({r}_{id}^{n+\theta },\left({\boldsymbol {h}}_{im}^{n+\theta }\cdot \nabla \right)q\right)}_{{\tilde {\Omega }}_{i}}=0\quad \quad \quad \quad \quad \quad \quad i=$$1,2$

where ${\textstyle x^{n+\theta }:=\theta x^{n+1}+\left(1-\theta \right)x^{n}}$ and ${\textstyle {\delta }_{t}x^{n}:={\frac {1}{\delta t}}\left(x^{n+1}-x^{n}\right)}$ .

This problem is nonlinear due to the convective term and the evolution of the free surface. Prior to the finite element discretization we linearize it using a Picard method, which leads to an Oseen problem for each iteration step. Several strategies can be adopted to deal with the two iterative algorithms involved in the ODDLS method: the overlapping domain decomposition iterative scheme and the linearization scheme (Picard). In our case, for each iteration step of the linearization scheme the domain decomposition scheme is also updated. This simultaneous iterative strategy requires to complete Eq. (3-50) with a stop criteria for the Picard iteration. The same iteration index is used for both schemes.

The adopted strategy reduces noteworthy the computational effort versus other nested iterative schemes. We denote by a superscript n the time step at which the algorithmic solution is computed and by a superscript j the iteration step of the domain decomposition scheme (which for a simultaneous iterative strategy will be the same as for the Picard iteration step). The form of the iterative scheme is as follows

	${\left({\rho }_{i}^{n,j-1}{\delta }_{t}{\boldsymbol {u}}_{i}^{n,j},{\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{i}}+$$\,{\left\langle {\rho }_{i}^{n}\left({\boldsymbol {u}}_{i}^{n+\theta ,j-1}\cdot \nabla \right){\boldsymbol {u}}_{i}^{n+\theta ,j},{\boldsymbol {v}}\right\rangle }_{{\tilde {\Omega }}_{i}}+$$\,{\left({\mathit {\boldsymbol {\sigma }}}_{i}^{n+\theta ,j},\nabla {\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{i}}+$${\frac {1}{2}}{\left({\boldsymbol {r}}_{im}^{n+\theta ,j},\left({\boldsymbol {h}}_{im}^{n+\theta ,j}\cdot \nabla \right){\boldsymbol {v}}\right)}_{{\tilde {\Omega }}_{i}}+$$\,{\left({\tilde {\boldsymbol {t}}}{\boldsymbol {,v}}\right)}_{{\Gamma }_{iN}}+\,{\left({\tilde {\boldsymbol {t}}}_{\boldsymbol {1}}{\boldsymbol {g+}}{\tilde {\boldsymbol {t}}}_{\boldsymbol {2}}{\boldsymbol {s,v}}\right)}_{{\Gamma }_{iM}}+$$\,\left(i-1\right){\left({\tilde {\mathit {\boldsymbol {t}}}}{\mathit {\boldsymbol {,v}}}\right)}_{{\tilde {\Gamma }}_{i}}=$$\,{\left\langle {\rho }_{i}^{n,j-1}{\boldsymbol {f}},{\boldsymbol {v}}\right\rangle }_{{\tilde {\Omega }}_{i}}$	(3-53)
	${\left(q,\nabla \cdot {\boldsymbol {u}}_{i}^{n+\theta ,j}\right)}_{{\tilde {\Omega }}_{i}}+$${\frac {1}{2}}{\left({r}_{id}^{n+\theta ,j},\left({\boldsymbol {h}}_{im}^{n+\theta ,j}\cdot \nabla \right)q\right)}_{{\tilde {\Omega }}_{i}}=0\quad \quad \quad \quad \quad \,i=$$1,2$

The Galerkin approximation of Eq. (3-53) is straightforward and stable, thanks to the FIC stabilized formulation [2-7]. Equal polynomial order can be used for the approximation of the velocities and the pressure. Based on the finite element partition introduced in Eq. (3-39), let $V_{h}\subset V$ and $Q_{h}\subset Q$ be two finite dimensional conforming finite element spaces, being ${\textstyle h}$ the maximum of the elements partition diameters.

Let ${\tilde {V}}_{h,i}^{0}:=\left\{{\boldsymbol {v}}\in V_{h}\vert {{\boldsymbol {v}}\vert }_{\partial {\tilde {\Omega }}_{i}}=\right.$$\left.{\boldsymbol {0}}\right\}{\begin{array}{cc},&Q_{h,i}^{0}:=\left\{q\in Q_{h}\vert {q\vert }_{\partial {\tilde {\Omega }}_{i}}=0\right\}\end{array}}$ be the finite element spaces for the test functions. The discretized form of Eq. 94 is

	${\left({\rho }_{h,i}^{n,j-1}{\delta }_{t}{\boldsymbol {u}}_{h,i}^{n,j},{\boldsymbol {v}}_{h}\right)}_{{\tilde {\Omega }}_{i}}+$$\,{\left\langle {\rho }_{h,i}^{n}\left({\boldsymbol {u}}_{h,i}^{n+\theta ,j-1}\cdot \nabla \right){\boldsymbol {u}}_{h,i}^{n+\theta ,j},{\boldsymbol {v}}_{h}\right\rangle }_{{\tilde {\Omega }}_{i}}+$$\,{\left({\mathit {\boldsymbol {\sigma }}}_{h,i}^{n+\theta ,j},\nabla {\boldsymbol {v}}_{h}\right)}_{{\tilde {\Omega }}_{i}}+$${\frac {1}{2}}{\left({\boldsymbol {r}}_{h,im}^{n+\theta ,j},\left({\boldsymbol {h}}_{h,\,\,\,ega\,omain\,decompositionim}^{n+\theta ,j}\cdot \nabla \right){\boldsymbol {v}}_{h}\right)}_{{\tilde {\Omega }}_{i}}+$$\,{\left({\tilde {\boldsymbol {t}}}{\boldsymbol {,}}{\boldsymbol {v}}_{h}\right)}_{{\Gamma }_{iN}}+$$\,{\left({\tilde {\boldsymbol {t}}}_{\boldsymbol {1}}{\boldsymbol {g+}}{\tilde {\boldsymbol {t}}}_{\boldsymbol {2}}{\boldsymbol {s,}}{\boldsymbol {v}}_{h}\right)}_{{\Gamma }_{iM}}+$$\,\left(i-1\right){\left({\tilde {\boldsymbol {t}}}{\boldsymbol {,}}{\boldsymbol {v}}_{h}\right)}_{{\tilde {\Gamma }}_{i}}=$$\,{\left\langle {\rho }_{h,i}^{n,j-1}{\boldsymbol {f}},{\boldsymbol {v}}_{h}\right\rangle }_{{\tilde {\Omega }}_{i}}$
	${\left({q}_{h},\nabla \cdot {\boldsymbol {u}}_{h,i}^{n+\theta ,j}\right)}_{{\tilde {\Omega }}_{i}}+$${\frac {1}{2}}{\left({r}_{h,id}^{n+\theta ,j},\left({\boldsymbol {h}}_{h,im}^{n+\theta ,j}\cdot \nabla \right){q}_{h}\right)}_{{\tilde {\Omega }}_{i}}=0\quad \quad \quad \quad \quad i=$$1,2$	(3-54)
	$\forall ({\boldsymbol {v}}_{h},{q}_{h})\in {\boldsymbol {V}}_{h}\times {Q}_{h}$

Remark

Updating of the level set equation is done once the convergence of the Navier-Stokes equations is obtained. Therefore the complete iteration loop is as follows:

Step 1. Solve equations (54) in ${\tilde {\Omega }}_{1}$.

Step 2. Solve equations (54) in ${\tilde {\Omega }}_{2}$.

Step 3. Check convergence on pressure and velocity fields. If it is satisfied, go to Step 4, otherwise go to Step 1.

Step 4. Solve equation (51).

Step 5. Check convergence on level set equation. If it is satisfied, advance to the next time step, otherwise go to Step 1.

4. Heat Transfer Solver (heatrans module)

4.1 Governing equations

Tdyn solves the transient Heat Transfer Equations in a given domain Ω:={Ω_F U Ω_S}(being Ω_F the fluid domain and Ω_Sthe solid domain) and time interval (0, t):

	$\rho {C}_{P}\left({\frac {\partial \theta }{\partial t}}+\,\left({\boldsymbol {u}}\cdot \nabla \right)\theta \right)-$$\,\nabla \cdot \left({\kappa }_{F}\nabla \theta \right)={\rho }_{F}{q}_{F}\quad \,in\,{\Omega }_{F}\times (0,t)$	(4-1)
	${\rho }_{S}C{\frac {\partial \theta }{\partial t}}-\,\nabla \cdot \left({\kappa }_{S}\nabla \theta \right)=$$\,{\rho }_{S}{q}_{S}\quad \,in\,{\Omega }_{S}\times (0,t)$

where θ = θ (x, t) denotes the temperature field, ρ_S the (constant) solid density, k_S the solid thermal conduction tensor, k_F the fluid thermal conduction constant, C_P the fluid specific heat constant, C the solid specific heat constant and q_S and q_F the solid and fluid, volumetric heat source distributions respectively. The above equations need to be combined with the following boundary conditions:

	$\theta ={\theta }_{c}\quad \quad \quad \quad \quad \quad in\,{\Gamma }_{\theta }\times (0,t)$
	$n{\kappa }_{F}\nabla \theta ={f}_{F}\quad \quad \quad \,in\,{\Gamma }_{F}\times (0,T)$	(4-2)
	$n{\kappa }_{S}\nabla \theta ={f}_{S}\quad \quad \quad \,\,in\,{\Gamma }_{S}\times (0,T)$
	$\theta \left(x,0\right)=\,{\theta }_{0}\left(x\right)\quad \,in\,\Omega \times \lbrace 0\rbrace$

In the above equations, Γ := ∂Ω denotes the boundary of the domain Ω, with n the normal unit vector, θ_c is the temperature field on Γ_θ (the part of the boundary of Dirichlet type, or prescribed temperature type), f_F the prescribed heat flux on Γ_F (prescribed heat flux fluid boundary), f_S the prescribed heat flux on Γ_S (prescribed heat flux solid boundary) and θ₀ the initial temperature field.

The spatial discretization of the heat transfer equations is done by means of the finite element method, while for the time discretization implicit first and second order schemes have been implemented.

4.2 Stabilized heat transfer equations

The Finite Increment Calculus theory presented above is also valid for the heat transfer balance equations Eq. (4-1). The stabilized form of the governing differential equations in 3 dimensions is written as follows:

	${r}_{F}-{\frac {{h}_{Fj}}{2}}{\frac {\partial {r}_{F}}{\partial {x}_{j}}}=$$0\quad \quad \quad \quad \quad \quad \,in\,{\Omega }_{F},\,j=1,2,3$	(4-3)
	${r}_{S}-{\frac {{h}_{Sj}}{2}}{\frac {\partial {r}_{S}}{\partial {x}_{j}}}=$$0\quad \quad \quad \quad \quad \quad \,\,in\,{\Omega }_{S},\,j=1,2,3$

Where

	${\textstyle {r}_{F}=\rho {C}_{P}\left({\frac {\partial \theta }{\partial t}}\right)-}$$\,{\kappa }_{F}\Delta \theta -\rho {q}_{F}$	(4-4)
	${r}_{S}=\,{\rho }_{S}C{\frac {\partial \theta }{\partial t}}-\nabla \cdot \left({\kappa }_{S}\nabla \theta \right)-$${\rho }_{S}{q}_{S}$

and the new boundary conditions are:

	$\theta =\,{\theta }_{c}\quad \,in\,{\Gamma }_{\theta }\times (0,t)$
	$n{\kappa }_{F}\nabla \theta +n{\frac {{h}_{Fj}}{2}}{r}_{F}={f}_{F}\quad \,in\,{\Gamma }_{F}\times (0,T)$	(4-5)
	$n{\kappa }_{S}\nabla \theta +n{\frac {{h}_{Sj}}{2}}{r}_{S}={f}_{S}\quad \,\,in\,{\Gamma }_{S}\times (0,T)$
	$\theta \left(x,0\right)=\,{\theta }_{0}\left(x\right)\quad \quad \quad \quad \,in\,\Omega \times \lbrace 0\rbrace$

The stabilized Heat transfer balance equations are discretized in time in the standard way, to give the following equations:

	${\theta }^{j+1,n+1}={\theta }^{n}-{\frac {\Delta t}{\rho {C}_{P}}}{\left[\rho {C}_{P}\left({\boldsymbol {u}}\cdot \nabla \right)\theta -\nabla \cdot \left({\kappa }_{F}\nabla \theta \right)-\,\rho {q}_{F}\right]}^{j,n+\theta }$	(4-6)
	${\theta }^{j+1,n+1}={\theta }^{n}-{\frac {\Delta t}{\rho C}}{\left[\nabla \cdot \left({\kappa }_{S}\nabla \theta \right)-\,\rho {q}_{S}\right]}^{j,n+\theta }$

4.3 Radiation models

Heat transfer by radiation can be taken into account within Tdyn by using one of the following radiation models.

4.3.1 P-1 radiation model

The energy conservation equation in terms of the temperature ${\textstyle T}$, can be written as:

$\rho c\left[{\frac {\partial T}{\partial t}}+{u}_{j}{\frac {\partial T}{\partial {x}_{j}}}\right]-$${\frac {\partial }{\partial {x}_{j}}}\left(k\cdot {\frac {\partial T}{\partial {x}_{j}}}\right)+$$rT+p{\frac {\partial {u}_{j}}{\partial {x}_{j}}}-\Phi ={s}_{T}$

(4-7)

Being ${\textstyle c}$ the heat capacity, ${\textstyle k}$ the thermal conductivity, ${\textstyle r}$ the heat reaction term, ${\textstyle s}$ the heat source term, and ${\textstyle \Phi }$ the viscous dissipation function given by:

$\Phi ={\mathit {\boldsymbol {\tau }}}:\nabla u$

(4-8)

Where ${\textstyle {\mathit {\boldsymbol {\tau }}}}$ is the viscous stress tensor of the fluid.

The FEM formulation for the temperature equation can be written as:

	$\int _{\Omega }^{}W{\mathit {\boldsymbol {\rho c}}}{\frac {{\mathit {\boldsymbol {T}}}^{n+\theta ,k}-{\mathit {\boldsymbol {T}}}^{n}}{\theta \Delta t}}d\Omega +$$\int _{\Omega }^{}W{\mathit {\boldsymbol {\rho c}}}\left[{\mathit {\boldsymbol {u}}}_{j}^{n+\theta ,k-1}{\frac {\partial }{\partial {x}_{j}}}{\mathit {\boldsymbol {T}}}^{n+\theta ,k}\right]d\Omega +$$\int _{\Omega }^{}\tau {\mathit {\boldsymbol {\rho c}}}\left[{\mathit {\boldsymbol {u}}}_{j}^{n+\theta ,k-1}{\frac {\partial W}{\partial {x}_{j}}}\right]\left[{\mathit {\boldsymbol {u}}}_{j}^{n+\theta ,k-1}{\frac {\partial }{\partial {x}_{j}}}{\mathit {\boldsymbol {T}}}^{n+\theta ,k}\right]d\Omega +$$\int _{\Omega }^{}k{\frac {\partial W}{\partial {x}_{j}}}{\frac {\partial }{\partial {x}_{j}}}{\mathit {\boldsymbol {T}}}^{n+\theta ,k}d\Omega +$$\int _{\Omega }^{}W{\mathit {\boldsymbol {r}}}{\mathit {\boldsymbol {T}}}^{n+\theta ,k}d\Omega =$$-\int _{\Omega }^{}W{s}_{T+\Phi +p}d\Omega +\int _{\Gamma }^{}{n}_{j}W{\mathit {\boldsymbol {q}}}_{\mathit {\boldsymbol {j}}}^{\mathit {\boldsymbol {\ast }}}{\mathit {\boldsymbol {\,}}}d\Gamma +$$\int _{\Omega }^{}\tau {\mathit {\boldsymbol {\rho c}}}\left[{\mathit {\boldsymbol {u}}}_{j}^{n+\theta ,k-1}{\frac {\partial W}{\partial {x}_{j}}}\right]{\mathit {\boldsymbol {\varsigma }}}_{i}^{n+\theta ,k-1}d\Omega$	(4-9)
	$\int _{\Omega }^{}W{\mathit {\boldsymbol {\varsigma }}}_{i}^{n+1,k}d\Omega =\int _{\Omega }^{}W\left[{\mathit {\boldsymbol {u}}}_{j}^{n+\theta ,k}{\frac {\partial }{\partial {x}_{j}}}{\mathit {\boldsymbol {T}}}^{n+\theta ,k}\right]d\Omega$	(4-10)

Where ${\textstyle {s}_{T+\Phi +p}}$ is the total heat/sink source term, including viscous dissipation and rate of work for volume change, and ${\textstyle {\mathit {\boldsymbol {q}}}^{\mathit {\boldsymbol {\ast }}}}$ is the heat flux at the walls.

The so called P-N approximation reduces the integral equations of radiative transfer in media to differential equations by approximating the transfer relations with a finite set of moment equations. The starting point in the following equation for the variation of intensity of radiation ${\textstyle {i}^{'}}$ (W/m2·sr) at a position ${\textstyle {\mathit {\boldsymbol {x}}}}$, along the ${\textstyle {\boldsymbol {s}}}$ direction [20]:

${\frac {d{i}_{\lambda }^{'}}{d{\boldsymbol {s}}}}={a}_{\lambda }{i}_{\lambda b}^{'}-\left({a}_{\lambda }+\right.$$\left.{\sigma }_{s\lambda }\right){i}_{\lambda }^{'}+{\frac {{\sigma }_{s\lambda }}{4\pi }}\int _{0}^{4\pi }{i}_{\lambda }^{'}\cdot \phi \cdot d\omega$

(4-11)

Where ${\textstyle {a}_{\lambda }}$ is the absorption coefficient and ${\textstyle {\sigma }_{s\lambda }}$ is the dispersion coefficient for the given wave-length ${\textstyle \lambda }$ , and the sub-index ${\textstyle b}$ refers to the black body. If the medium is assumed gray, with uniform scattering and absorption coefficients, and to scatter isotropically, above equation can be integrated over all ${\textstyle \lambda }$ , to give:

${\frac {d{i}^{'}}{d{\boldsymbol {s}}}}=a{i}_{b}^{'}-\left(a+{\sigma }_{s}\right){i}^{'}+$${\frac {{\sigma }_{s}}{4\pi }}\int _{0}^{4\pi }{i}^{'}\cdot d\omega$

(4-12)

Where the first term at the right hand side represents the gain of intensity of radiation by emission, the second term the losses by absorption and scattering, and the third term the gain by scattering into ${\textstyle s}$ direction.

To develop the P-N method, the intensity ${\textstyle {i}^{'}}$ is expanded in an orthogonal series of spherical harmonics. The P-1 approximation is obtained by retaining the first two terms of the series, and gives [20]:

$\sum _{i=1}^{3}{\frac {\partial {i}^{'(i)}}{\partial {x}_{i}}}=a\left(4\pi {i}_{b}^{'}-\right.$$\left.{i}^{'(0)}\right)$

(4-13)

In the above equation ${\textstyle {i}^{'(0)}}$ and ${\textstyle {i}^{'(i)}}$ represents the zeroth and first moments of the radiation intensity:

	${i}^{'(0)}=G=\int _{0}^{4\pi }{i}^{'}{\boldsymbol {\cdot }}d\omega$	(4-14)
	${i}^{'(i)}={q}_{r,i}=\int _{0}^{4\pi }{i}^{'}{\boldsymbol {\cdot }}{s}_{i}{\boldsymbol {\cdot }}d\omega ,\quad {\mathit {\boldsymbol {q}}}_{r}=$$\int _{0}^{4\pi }{i}^{'}{\boldsymbol {\cdot s\cdot }}d\omega$	(4-15)

Where ${\textstyle {\mathit {\boldsymbol {q}}}_{r}}$ is the radiative heat flux vector, and ${\textstyle G}$ is the incident radiation. Therefore, the P-1 approximation yields:

$-\nabla \cdot {\mathit {\boldsymbol {q}}}_{r}=a\left(G-4\pi {i}_{b}^{'}\right)$

(4-16)

The first term in the above equation, represents the net radiative energy supplied per unit volume.

Furthermore, P-1 model gives a following additional relation between ${\textstyle G}$ and ${\textstyle {\mathit {\boldsymbol {q}}}_{r}}$:

${\mathit {\boldsymbol {q}}}_{r}=-{\frac {1}{3\left(a+{\sigma }_{s}\right)}}\nabla G=$$-\gamma \nabla G$

(4-17)

Combining above relations, the resulting transport equation for ${\textstyle G}$ is:

$\nabla \cdot \left(\gamma \nabla G\right)-aG=-4e\sigma {T}^{4}$

(4-18)

Where ${\textstyle \sigma }$ is the Stefan-Boltzmann constant (5.672 ·10^-8  W/m²·K⁴), ${\textstyle \epsilon }$ is the emissivity of the material, and the black body radiant energy has been evaluated using:

${i}_{b}^{'}={\frac {\sigma {T}^{4}}{\pi }}$

(4-19)

The FEM formulation for the incident radiation equation can be written as:

$\int _{\Omega }^{}\gamma {\frac {\partial W}{\partial {x}_{j}}}{\frac {\partial }{\partial {x}_{j}}}{\mathit {\boldsymbol {G}}}^{k}d\Omega +$$\int _{\Omega }^{}Wa{\mathit {\boldsymbol {G}}}^{k}d\Omega =\int _{\Omega }^{}W4\epsilon \sigma {\left({\mathit {\boldsymbol {T}}}^{n+\theta ,k}\right)}^{4}d\Omega +$$\int _{\Gamma }^{}W{q}_{r,w}{\mathit {\boldsymbol {\,}}}d\Gamma$

(4-20)

Where ${\textstyle {q}_{r,w}}$ is the flux of the incident radiation al the walls.

From the solution of the previous equation, the heat source term that must be added to the energy equation (temperature) is calculated as follows:

${s}_{T}=-\nabla \cdot {\mathit {\boldsymbol {q}}}_{r}=aG-4\epsilon \sigma {T}^{4}$

(4-21)

In the implemented model the emissivity is assumed to be equal to the absorptivity ${\textstyle a}$, and therefore above equations can be written as:

	$\int _{\Omega }^{}\gamma {\frac {\partial W}{\partial {x}_{j}}}{\frac {\partial }{\partial {x}_{j}}}{\mathit {\boldsymbol {G}}}^{k}d\Omega +$$\int _{\Omega }^{}Wa{\mathit {\boldsymbol {G}}}^{k}d\Omega =\int _{\Omega }^{}W4a\sigma {\left({\mathit {\boldsymbol {T}}}^{n+\theta ,k}\right)}^{4}d\Omega +$$\int _{\Gamma }^{}W{g}_{w}{\mathit {\boldsymbol {\,}}}d\Gamma$	(4-22)
	${s}_{T}=-\nabla \cdot {\mathit {\boldsymbol {q}}}_{r}=a\left(G-4\sigma {T}^{4}\right)$	(4-23)

The boundary condition for the incident radiation equation at the walls is given by (assuming an inwards normal vector ${\textstyle {\boldsymbol {n}}}$):

${g}_{w}=\gamma \nabla G\cdot {\boldsymbol {n}}=\gamma {\frac {\partial G}{\partial n}}$

(4-24)

Where ${\textstyle {g}_{w}}$ can be calculated, assuming that walls are diffuse gray surfaces, as [21]:

${g}_{w}={\frac {{\epsilon }_{w}}{2\left(2-{\epsilon }_{w}\right)}}\left(4\sigma {T}_{w}^{4}-\right.$$\left.{G}_{w}\right)$

(4-25)

Where ${\textstyle {\epsilon }_{w}}$ is the wall emissivity, and ${\textstyle {T}_{w}}$ is the wall temperature.

Remark: It is usually assumed that the emissivity at the inlet and outlets is 1 (black body).

The final FEM formulation for the incident radiation equation can be written as:

$\int _{\Omega }^{}\Gamma {\frac {\partial W}{\partial {x}_{j}}}{\frac {\partial }{\partial {x}_{j}}}{\mathit {\boldsymbol {G}}}^{k}d\Omega +$$\int _{\Omega }^{}Wa{\mathit {\boldsymbol {G}}}^{k}d\Omega +{\frac {{\epsilon }_{w}}{2\left(2-{\epsilon }_{w}\right)}}\int _{\Gamma }^{}W{\mathit {\boldsymbol {G}}}^{k}{\mathit {\boldsymbol {\,}}}d\Gamma =$$\int _{\Omega }^{}W4a\sigma {\left({\mathit {\boldsymbol {T}}}^{n+\theta ,k}\right)}^{4}d\Omega +$${\frac {{\epsilon }_{w}}{2\left(2-{\epsilon }_{w}\right)}}\int _{\Gamma }^{}W4\sigma {\mathit {\boldsymbol {T}}}_{w}^{4}{\mathit {\boldsymbol {\,}}}d\Gamma$

(4-26)

4.3.2 Surface to surface (S2S) radiation model

The surface-to-surface (S2S) radiation model can be used to account for the radiation exchange in an enclosure of gray-diffuse surfaces. This energy exchange depends on:

• Surfaces size

• Separation distance

• Surfaces orientation

All these factors are casted together into a geometric function called the view factor ${\textstyle {F}_{ij}}$ that can be computed as follows:

${F}_{ij}={\frac {1}{{A}_{i}}}\int _{{A}_{i}}^{}\int _{{A}_{j}}^{}{\frac {cos{\theta }_{i}\cdot cos{\theta }_{j}}{\pi {r}^{2}}}{\delta }_{ij}d{A}_{i}d{A}_{j}$

(4-27)

where ${\textstyle {\delta }_{ii}=0}$, ${\textstyle {\delta }_{ij}=}$$0$ when i,j are blocking surfaces and ${\textstyle {\delta }_{ij}=}$$1$ when i,j are non-blocking surfaces.

The S2S model main assumptions are:

• All surfaces are gray and diffuse. According to this, if a certain amount of radiation is incident on a surface, then a fraction is reflected, a fraction is absorbed, and a fraction is transmitted.

• Absorption, emission and scattering of radiation by the medium can be ignored (non-participating media)

• Emissivity and absorptivity are independent of the radiation wave length

• Emissivity is taken to be identical to absorptivity ${\textstyle \left(\epsilon \equiv \alpha \right)}$

• Reflectivity is independent of the outgoing or the incoming directions

In most applications of the S2S model the surfaces are opaque to thermal radiation (in the infrared spectrum). In such a situation the following relations are fulfilled:

• ${\textstyle \tau =0}$

• ${\textstyle \alpha =\epsilon }$

• ${\textstyle \rho =1-\epsilon }$

and the only independent parameter characteristic of the material is the emissivity ${\textstyle \epsilon }$ .

When the S2S model is used, it is also possible to define a “partial enclosure” that allows disabling the view factor calculation for walls with negligible emission/absorption or walls that have uniform temperature.

The energy flux leaving a given surface is composed of directly emitted and reflected energy. Hence, the total energy flux leaving the surface ${\textstyle k}$ ( ${\textstyle {q}_{out,k}}$) has two contributions, one depending on its own temperature and emissivity, and another one depending on the energy flux incident on the surface from the surroundings.

${q}_{out,k}=\,{\epsilon }_{k}\cdot \sigma \cdot {T}_{k}^{4}+{\rho }_{k}\cdot {q}_{in,k}$

(4-28)

The energy flux incident from the surroundings is given by:

${q}_{in,k}=\sum _{j=1}^{N}{F}_{kj}\cdot {q}_{out,j}$

(4-29)

where N is the number of surface clusters (i.e. element patches) considered. Hence, the total energy flux leaving a given surface can be expressed as follows:

${q}_{out,k}=\,{\epsilon }_{k}\cdot \sigma \cdot {T}_{k}^{4}+{\rho }_{k}\sum _{j=1}^{N}{F}_{kj}\cdot {q}_{out,j}$

(4-30)

${q}_{out,k}-{\rho }_{k}\sum _{j=1}^{N}{F}_{kj}\cdot {q}_{out,j}=\,{\epsilon }_{k}\cdot \sigma \cdot {T}_{k}^{4}$

(4-31)

This can be rewritten in matrix form as:

${\boldsymbol {K}}\cdot {\boldsymbol {J}}={\boldsymbol {E}}$

(4-32)

where ${\textstyle {\boldsymbol {K}}}$ is a NxN matrix, ${\textstyle {\boldsymbol {J}}}$ is the radiosity vector, and ${\textstyle {\boldsymbol {E}}}$ is the emissive power.

In the above equations, ${\textstyle {T}_{k}^{}}$and the material dependent property ${\textstyle {\epsilon }_{k}}$ are calculated for each patch by averaging the element’s values as follows:

${T}_{k}={\frac {\sum _{e}^{}{A}_{e}\cdot {T}_{e}}{\sum _{e}^{}{A}_{e}}},\,{\quad \quad \epsilon }_{k}=$${\frac {\sum _{e}^{}{A}_{e}\cdot {\epsilon }_{e}}{\sum _{e}^{}{A}_{e}}}$

(4-33)

Finally, equation 5 can be solved to obtain the outgoing fluxes ${\textstyle {q}_{out,k}}$ . And the net heat fluxes can be obtained as follows:

${q}_{net,k}={{q}_{out,k}-{q}_{in,k}=q}_{out,k}-\sum _{j=1}^{N}{F}_{kj}\cdot {q}_{out,j}$

(4-34)

where Eq. (3) has been used to express the incident fluxes ${\textstyle {q}_{in,k}}$.

The solution process outlined above provides the net heat flux at the element/patch level. This can be applied as boundary condition of the boundary value problem to be solved by the finite elements method after proper transformation of element/patch values to nodal values.

5. Species Advection Solver (advect module)

5.1 Governing equations

Tdyn solves the transient Species Advection Equations in a given fluid domain Ω_F for a number of different species, and time interval (0, t):

${\frac {\partial \varphi }{\partial t}}+\,\left({\boldsymbol {u}}\cdot \nabla \right)\varphi +$$\left(c\cdot \varphi \cdot g\cdot \nabla \right)\varphi -\nabla \cdot \left({\kappa }_{P}\nabla \varphi \right)=$$\,{q}_{F}\quad \,in\,{\Omega }_{F}\times (0,t)$

(5-1)

where φ = φ (x, t) denotes the concentration of species field, c the decantation coefficient, k_S the total diffusion coefficient (including turbulent effects) and q_P a volumetric source. The above equations need to be combined with the standard boundary conditions.

As in the previous examples, the spatial discretization of the Species Advection equations has been done by means of the finite element method, while for the time discretization implicit first and second order schemes have been implemented. Problems with dominating convection are stabilized by the Finite Calculus (FIC) method, presented above.

6. References

[1] R. Codina "A discontinuity capturing crosswind dissipation for the finite element solution of the convection diffusion equation", Comp. Meth. Appl. Mech. Eng., pages 110:325-342, 1993.

[2] J. García, "A Finite Element Method for the Hydrodynamic Analysis of Naval Structures" (in Spanish), Ph.D. Thesis, Univ. Politécnica de Cataluña, Barcelona, December 1999.

[3] C. Hirsch "Numerical computation of internal and external flows", John Wiley and Sons Ltd., Chichester, UK, 1991.

[4] C. Kleinstreuer "Engineering Fluid Dynamics, An Interdisciplinary Systems Approach", Cambridge University Press, Cambridge, UK, 1997.

[5] E. Oñate "A stabilized finite element method for incompressible viscous flows using a finite increment calculus formulation. Research report No. 150", CIMNE, Barcelona, January 1999.

[6] García, J., Oñate, E., Sierra, H., Sacco, C. and Idelsohn, S. “A Stabilized Numerical Method for Analysis of Ship Hydrodynamics”. Proceedings Eccomas Conference on CFD, 7-11 September 1998, Athens, John Wiley, 1998.

[7] E. Oñate, M. Manzan, "Stabilization techniques for finite element analysis of convection-diffusion problems", Publication CIMNE No. 183, Barcelona, February 2000.

[8] O.C. Zienkiewicz, T.L.Taylor "The finite element method". McGrawHill Book Company, Volume 1, Maidenhead, UK, 1989.

[9] E. Oñate, J. García and S. Idelsohn, “On the stabilisation of numerical solutions of advective-diffusive transport and fluid flow problems”. Comp.Meth. Appl. Mech. Engng, 5, 233-267, 1998.

[10] J.E. Bardina, P.G. Huang and T.J. Coakley, “Turbulence Modelling Validation, Testing, and Development”. NASA technical Memorandum 110446, Moffett Field, California, April 1997.

[11] D.C. Wilcox "Turbulence Modelling for CFD". DCW industries, Inc. 5354 Palm Drive, La Cañada, California, 1993.

[12] M. Vázquez, R. Codina and O.C. Zienkiewicz "Numerical modelling of Compressible Laminar and Turbulent Flow: The CBS Algorithm". Monograph CIMNE Nº 50, Barcelona, April 1999.

[13] J. García-Espinosa and E. Oñate. “An Unstructured Finite Element Solver for Ship Hydrodynamics”. Jnl. Appl. Mech. Vol. 70: 18-26. (2003).

[14] E. Oñate, J. García-Espinosa y S. Idelsohn. “Ship Hydrodynamics”. Encyclopedia of Computational Mechanics. Eds. E. Stein, R. De Borst, T.J.R. Hughes. John Wiley & Sons (2004).

[15] O.C. Zienkiewicz and Codina, R., “A general algorithm for compressible and incompressible flow. Part I, The Split, Characteristic Based Scheme”, Int. Journ. Num. Meth. Fluids, 20, 869-885 (1995).

[16] O. Soto, R. Lohner, J Cebral. An Implicit Monolithic Time Accurate Finite Element Scheme for Incompressible Flow Problems. 15th AIAA Computational Fluid Dynamics Conference. Anaheim (USA) 2001.

[17] J. García-Espinosa, E. Oñate and J. Bloch Helmers. “Advances in the Finite Element Formulation for Naval Hydrodynamics Problems”. International Conference on Computational Methods in Marine Engineering (Marine 2005). June 2005. Oslo (Norway).

[18] J. García-Espinosa, A. Valls and E. Oñate. “ODDLS: A new unstructured mesh finite element method for the analysis of free surface flow problems”. Int. Jnl. Num. Meth. Engng. 76(9): 1297-1327 (2008).

[20] R. Siegel and J.R. Howell. “Thermal radiation heat transfer”. Hemisphere Publishing Corporation, Washington DC, 1992.

[21] M.F. Modest. “Radiative heat transfer”. McGraw Hill, 2003.

[22] J.H. Lienhard. “A heat transfer text book”. Phlogiston Press, Cambridge 2008.