Latest revision as of 13:59, 14 June 2019

Abstract

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

Keywords: Particle Finite Element Method, Coupled Thermal Analysis, Quasi and Fully Incompressible

1 INTRODUCTION

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

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

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

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

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

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

2 GOVERNING EQUATIONS

We write the governing equations for a quasi-incompressible Newtonian flow problem in the Lagrangian description as follows [3,46].

Momentum equations

In Eq.(1), ${\textstyle \Omega }$ is the analysis domain with boundary ${\textstyle \Gamma }$, ${\textstyle v_{i}}$ and ${\textstyle b_{i}}$ are the velocity and body force components along the ${\textstyle i}$th Cartesian axis, ${\textstyle \rho }$ is the density, ${\textstyle n_{s}}$ is the number of space dimensions (i.e. ${\textstyle n_{s}=3}$ for 3D problems) and ${\textstyle \sigma _{ij}}$ are the Cauchy stresses that are split in the deviatoric (${\textstyle s_{ij}}$) and pressure (${\textstyle p}$) components as

where ${\textstyle \delta _{ij}}$ is the Kronecker delta. Note that the pressure is assumed to be positive for a tension state. Summation of terms with repeated indices is assumed in Eq.(1) and in the following, unless otherwise specified.

The relationship between the deviatoric stresses ${\textstyle s_{ij}}$ and the strain rates ${\textstyle \varepsilon _{ij}}$ has the standard form for a Newtonian fluid,

where ${\textstyle \mu }$ is the viscosity and ${\textstyle \varepsilon _{v}}$ is the volumetric strain rate defined as ${\textstyle \varepsilon _{v}=\varepsilon _{ii}={\partial v_{i} \over \partial x_{i}}}$.

Mass balance equation

The standard mass balance equation for a quasi-incompressible fluid can be written as [3,7,46]

with

In Eq.(4b) ${\textstyle c}$ is the speed of sound in the fluid. For a fully incompressible fluid ${\textstyle c=\infty }$ and Eq.(4a) simplifies to the standard form, ${\textstyle \varepsilon _{v}=0}$. In our work we will retain the quasi-incompressible form of ${\textstyle r_{v}}$ of Eq.(4b) for convenience.

Thermal balance

where ${\textstyle T}$ is the temperature, ${\textstyle c}$ is the thermal capacity, ${\textstyle k}$ is the heat conductivity and ${\textstyle Q}$ is the heat source.

Boundary conditions

Mechanical problem

The boundary conditions at the Dirichlet (${\textstyle \Gamma _{v}}$) and Neumann (${\textstyle \Gamma _{t}}$) boundaries with ${\textstyle \Gamma =\Gamma _{v}\cup \Gamma _{t}}$ are

where ${\textstyle v_{i}^{p}}$ and ${\textstyle t_{i}^{p}}$ are the prescribed velocities and prescribed tractions on ${\textstyle \Gamma _{v}}$ and ${\textstyle \Gamma _{t}}$, respectively and ${\textstyle n_{j}}$ are the components of the unit normal vector to the boundary [3,7,46].

Thermal problem

where ${\textstyle \phi ^{p}}$ and ${\textstyle q_{n}^{p}}$ are the prescribed temperature and the prescribed normal heat flux at the boundaries ${\textstyle \Gamma _{\phi }}$ and ${\textstyle \Gamma _{q}}$, respectively and ${\textstyle n}$ is the direction normal to the boundary.

Remark 1. The term ${\textstyle {\frac {Dv_{i}}{Dt}}}$ in Eq.(1) is the material derivative of the ${\textstyle i}$th velocity component ${\textstyle v_{i}}$. This term is typically computed in a Lagrangian framework as

with

where ${\textstyle {}^{n}v_{i}}$ is the velocity of the material point that has the position ${\textstyle {}^{n}{\bf {x}}}$ at time ${\textstyle t={}^{n}t}$, where ${\textstyle {\bf {x}}}$ is the coordinates vector in a fixed Cartesian system [3,7,46].

3 STABILIZED MASS BALANCE EQUATION

In this work we will use the second order FIC form of the mass balance equation in space for a quasi-incompressible fluid [37,38], as well as the first order FIC form of the mass balance equation in time. These forms have the following expressions:

Second order FIC mass balance equation in space

First order FIC mass balance equation in time

Eq.(12a) is obtained by expressing the balance of mass in a rectangular domain of finite size with dimensions ${\textstyle h_{1}\times h_{2}}$ (for 2D problems), where ${\textstyle h_{i}}$ are arbitrary distances, and retaining up to third order terms in the Taylor series expansions used for expressing the change of mass within the balance domain.

Eq.(12b), on the other hand, is obtained by expressing the balance of mass in a space-time domain of infinitesimal length in space and finite dimension ${\textstyle \delta }$ in time [20]. The derivation of Eqs.(12) for a 1D problem are shown in [41].

The FIC terms in Eqs.(12) play the role of space and time stabilization terms respectively. In the discretized problem, the space dimensions ${\textstyle h_{i}}$ and the time dimension ${\textstyle \delta }$ are related to characteristic element dimensions and the time step increment, respectively as it will be explained later. Note that for ${\textstyle h_{i}=0}$ and ${\textstyle \delta =0}$ the standard infinitesimal form of the mass balance equation, ${\textstyle r_{v}=0}$, is obtained.

After some transformations the stabilized mass balance equation (12a) is written as [41]

where ${\textstyle \tau }$ is a stabilization parameter given by [41]

and ${\textstyle {\hat {r}}_{m_{i}}}$ is a static momentum term defined as

Eq.(13) is used as the starting point for deriving the stabilized FEM formulation as explained in the following sections.

4 VARIATIONAL EQUATIONS

4.1 Momentum equations

Multiplying Eq.(1) by arbitrary test functions ${\textstyle w_{i}}$ with dimensions of velocity and integrating over the analysis domain ${\textstyle \Omega }$ gives the weighted residual form of the momentum equations as [3,7,46]

Integrating by parts the term involving ${\textstyle \sigma _{ij}}$ and using the Neumann boundary conditions (7) yields the weak variational form of the momentum equations as

where ${\textstyle \delta \varepsilon _{ij}={\partial w_{i} \over \partial x_{j}}+{\partial w_{j} \over \partial x_{i}}}$ is an arbitrary (virtual) strain rate field. Eq.(17) is the standard form of the principle of virtual power [3,7,46].

Substituting the expression of the stresses from Eq.(2) into (17) gives

Eq.(18) can be written in matrix form as

In Eq.(19) ${\textstyle {\bf {w}},{\bf {v}},{\boldsymbol {\varepsilon }}}$ and ${\textstyle \delta {\boldsymbol {\varepsilon }}}$ are vectors containing the test functions, the velocities, the strain rates and the virtual strain rates respectively; ${\textstyle {\bf {b}}}$ and ${\textstyle {\bf {t}}^{p}}$ are body force and surface traction vectors, respectively; ${\textstyle {\bf {D}}}$ is the viscous constitutive matrix and ${\textstyle {\bf {m}}}$ is an auxiliary vector. These vectors are defined as (for 3D problems)

4.2 Mass balance equations

We multiply Eq.(13) by arbitrary (continuous) test functions ${\textstyle q}$ (with dimensions of pressure) defined over the analysis domain ${\textstyle \Omega }$. Integrating over ${\textstyle \Omega }$ gives

Integrating by parts the last integral in Eq.(21) and using (15) gives after some transformations [39,40]

Expression (24) holds for 2D and 3D problems. The terms involving the first and second material time derivative of the pressure and the boundary term in Eq.(24) are important for preserving the conservation of mass in free-surface flow problems [10,41].

4.3 Thermal balance equation

Application of the standard weighted residual method to the heat balance equations (5) and (9) leads, after standard operations, to [7,44]

where ${\textstyle {\hat {w}}}$ are the space weighting functions for the temperature.

5 FEM DISCRETIZATION

We discretize the analysis domain into finite elements with ${\textstyle n}$ nodes in the standard manner leading to a mesh with a total number of ${\textstyle N_{e}}$ elements and ${\textstyle N}$ nodes. In our work we will choose simple 3-noded linear triangles (${\textstyle n=3}$) for 2D problems and 4-noded tetrahedra (${\textstyle n=4}$) for 3D problems with local linear shape functions ${\textstyle N_{i}^{e}}$ defined for each node ${\textstyle i}$ (${\textstyle i=1,n}$) of element ${\textstyle e}$ [34,44]. The velocity components, the pressure and the temperature are interpolated over the mesh in terms of their nodal values in the same manner using the global linear shape functions ${\textstyle N_{j}}$ spanning over the elements sharing node ${\textstyle j}$ (${\textstyle j=1,N}$) [34,44,46]. In matrix form

where

with ${\textstyle {\bf {N}}_{j}=N_{j}{\bf {I}}_{n}}$ where ${\textstyle {\bf {I}}_{n}}$ is the ${\textstyle n\times n}$ unit matrix.

In Eq.(25) vectors ${\textstyle {\bar {v}}}$, ${\textstyle {\bar {p}}}$ and ${\textstyle {\bar {T}}}$ contain the nodal velocities, the nodal pressures and the nodal temperatures for the whole mesh, respectively and the upperindex denotes the nodal value for each vector or scalar magnitude.

Substituting Eqs.(24) into Eqs.(17), (22) and (25) and choosing a Galerkin formulation with ${\textstyle w_{i}=q={\hat {w}}_{i}=N_{i}}$ leads to the following system of algebraic equations

where ${\textstyle {\dot {\bar {\bf {a}}}}}$ and ${\textstyle {\ddot {\bar {\bf {a}}}}}$ denote the first and second material time derivatives of the components of a vector ${\textstyle {a}}$. The different matrices and vectors in Eqs.(26) are assembled from the element contributions given in Box 1.

Remark 2. The boundary terms of vector ${\textstyle {\bf {f}}_{p}}$ can be incorporated in the matrices of Eq.(26b). This leads to a non symmetrical set of equations. These boundary terms are computed here iteratively within the incremental solution scheme.

Remark 3. The presence of matrix ${\textstyle {\bf {M}}_{b}}$ in Eq.(26b) allows us to compute the pressure without the need of prescribing its value at the free surface. This eliminates the error introduced when the pressure is prescribed to zero in free boundaries, which leads to considerable mass losses in viscous flows [15].

Remark 4. For transient problems the stabilization parameter ${\textstyle \tau }$ of Eq.(14) is computed for each element ${\textstyle e}$ using ${\textstyle h=l^{e}}$ and ${\textstyle \delta =\Delta t}$ as

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

For steady state problems the stabilization parameter is computed with Eq.(27) substituting ${\textstyle \Delta t}$ by ${\textstyle {\frac {l^{e}}{\vert {\bf {v}}^{e}\vert }}}$. The characteristic boundary length ${\textstyle h_{n}}$ in the expression of ${\textstyle {\bf {f}}_{p}}$ (Box 1) has been taken equal to ${\textstyle l^{e}}$ in our computations.

6 TRANSIENT SOLUTION OF THE DISCRETIZED EQUATIONS

Eqs.(26) are solved in time with an implicit Newton-Raphson type iterative scheme [3,7,44,46]. The basic steps within a time increment ${\textstyle [n,n+1]}$ are:

• Initialize variables: ${\textstyle ({}^{n+1}{\bf {x}}^{1},{}^{n+1}{\bar {\bf {v}}}^{1},{}^{n+1}{\bar {\bf {p}}}^{1},{}^{n+1}{\bar {\bf {T}}}^{1},{}^{n+1}{\bar {\bf {r}}}_{m}^{1})\equiv \left\{{}^{n}{\bf {x}},{}^{n}{\bar {\bf {v}}},{}^{n}{\bar {\bf {p}}},{}^{n}{\bar {\bf {T}}},{}^{n}{\bar {\bf {r}}}_{m}\right\}}$.
• Iteration loop: ${\textstyle i=1,NITER}$.

  For each iteration.

Step 1. Compute the nodal velocity increments ${\displaystyle \Delta {\bar {v}}}$Step 1. Compute the nodal velocity increments Δ v ¯ {\displaystyle \Delta {\bar {v}}}

From Eq.(26a), we deduce

with the momentum residual ${\textstyle {\bar {\bf {r}}}_{m}}$ and the iteration matrix ${\textstyle {\bf {H}}_{v}}$ given by

Step 2. Update the nodal velocities

Step 3. Compute the nodal pressures ${\displaystyle {}^{n+1}{\bar {\bf {p}}}^{i+1}}$Step 3. Compute the nodal pressures n + 1 p ¯ i + 1 {\displaystyle {}^{n+1}{\bar {\bf {p}}}^{i+1}}

From Eq.(26b) we obtain

with

Step 4. Update the nodal coordinates

A more accurate expression for computing ${\textstyle {}^{n+1}{\bf {x}}^{i+1}}$ can be used involving the nodal accelerations [40].

Step 5.Compute the nodal temperatures

From Eq.(26c) we obtain

with

Step 6. Check convergence

Verify the following conditions:

where ${\textstyle e_{v}}$, ${\textstyle e_{p}}$ and ${\textstyle e_{T}}$ are prescribed error norms. In our examples we have set ${\textstyle e_{v}=e_{p}=e_{T}=10^{-3}}$.