Boundedness-Preserving Implicit Correction of Mesh-Induced Errors for VoF Based Heat and Mass Transfer

TL;DR

This paper introduces a boundedness-preserving correction method for mesh-induced skewness errors in Volume-of-Fluid simulations of heat and mass transfer, improving accuracy and stability on complex unstructured meshes.

Contribution

It develops a novel correction strategy for skewness errors in finite volume discretisation of advection and diffusion in VoF methods, ensuring bounded and accurate transport on complex meshes.

Findings

Effective correction of skewness errors demonstrated

Maintains boundedness and accuracy in heat and mass transfer simulations

Applicable to unstructured meshes of general topology

Abstract

Spatial discretisation of geometrically complex computational domains often entails unstructured meshes of general topology for Computational Fluid Dynamics (CFD). Mesh skewness is then typically encountered causing severe deterioration of the formal order of accuracy of the discretisation, or boundedness of the solution, or both. Particularly methods inherently relying on the accurate and bounded transport of sharp fields suffer from all types of mesh-induced skewness errors, namely both non-orthogonality and non-conjunctionality errors. This work is devoted to a boundedness-preserving strategy to correct for skewness errors arising from discretisation of advection and diffusion terms within the context of interfacial heat and mass transfer based on the Volume-of-Fluid methodology. The implementation has been accomplished using a second-order finite volume method with support for unstructured meshes of general topology. We examine and advance suitable corrections for the finite volume discretisation of a consistent single-field model, where both accurate and bounded transport due to diffusion and advection is crucial. In order to ensure consistency of both the volume fraction and the species concentration transport, i.e. to avoid artificial heat or species transfer, corrections are studied for both cases. The cross interfacial jump and adjacent sharp gradients of species concentration render the correction for skewness-induced diffusion and advection errors additionally demanding and has not so far been addressed in the literature.

Tables

Table 1. Table 1 . Overview of all methods used to discretise the face-interpolated value α f subscript 𝛼 𝑓 \alpha_{f} .

Interpolation Scheme	Correction Approach	Gradient Calculation
UDS	UC	GG
CDS	EC	LSF
CICSAM	SISC
LB

Table 2. Table 2 . Overview of all used methods varied against each other to discretise the diffusive terms of the CST model.

${𝐒}_f\cdot {\left(\nabla \varphi \right)}_f$

UC

NO

NO/NC

Table 3. (a) Laplacian

${𝐒}_f\cdot {\left(\nabla \varphi \right)}_f$

UC

NO

NO/NC

Table 4. (b) Diffusive flux

${⟨\text{D}⟩}_{\text{h,f}}^n$, ${\left(\text{K}c\right)}_f$
Interpolation Scheme	Correction Appr.	Gradient Calc.
CDS	UC	GG
	EC

Equations

$${\partial }_t\varphi +\nabla \cdot \left(\varphi \mathbf{u}\right)=-\nabla \cdot \left(-{\Gamma }_{\varphi }\nabla \varphi \right)+S_{\varphi }.$$

$${\int }_V\nabla \cdot \left(\varphi \mathbf{u}\right)dV$$

$$\approx \sum _f{\mathbf{S}}_f\cdot {\mathbf{u}}_f{\varphi }_f=:\sum _f{F}_f{\varphi }_f,$$

$${\varphi }_f={\varphi }_{f^{\prime }}+\nabla {\varphi }_{f^{\prime }}\cdot \mathbf{m},$$

$${\varphi }_f=\frac{1}{2}\left(\left[{\varphi }_P+\nabla {\varphi }_P\cdot \left({\mathbf{x}}_f-{\mathbf{x}}_P\right)\right]+\left[{\varphi }_N+\nabla {\varphi }_N\cdot \left({\mathbf{x}}_f-{\mathbf{x}}_N\right)\right]\right).$$

$${\oint }_S\mathbf{n}\cdot \left({\Gamma }_{\varphi }\nabla \varphi \right)dS\approx \sum _f{\mathbf{S}}_f\cdot {\left({\Gamma }_{\varphi }\nabla \varphi \right)}_f\approx \sum _f{\left({\Gamma }_{\varphi }\right)}_f{\mathbf{S}}_f\cdot {\left(\nabla \varphi \right)}_f.$$

$${\left({\Gamma }_{\varphi }\right)}_f{\mathbf{S}}_f\cdot {\left(\nabla \varphi \right)}_f={\left({\Gamma }_{\varphi }\right)}_f|{\mathbf{S}}_f|{\nabla }_f^{\perp }\varphi ={\left({\Gamma }_{\varphi }\right)}_f|{\mathbf{S}}_f|\frac{{\varphi }_N-{\varphi }_P}{|\mathbf{d}|}.$$

$${\mathbf{S}}_f=\mathbf{\Delta }+\mathbf{k},$$

$${\mathbf{S}}_f\cdot {\left(\nabla \varphi \right)}_f=\mathbf{\Delta }\cdot {\left(\nabla \varphi \right)}_f+\mathbf{k}\cdot {\left(\nabla \varphi \right)}_f=|\mathbf{\Delta }|\frac{{\varphi }_N-{\varphi }_P}{|\mathbf{d}|}+\mathbf{k}\cdot {\left(\nabla \varphi \right)}_f.$$

$${\varphi }_f={\varphi }_C+{\delta }_f\Psi \left(r\right)\left({\varphi }_D-{\varphi }_C\right),$$

$$0\le \Psi \left(r\right)\le \min \left(\frac{1}{{\delta }_f}r\frac{1-{\mathrm{Co}}_f}{{\mathrm{Co}}_f},\frac{1}{{\delta }_f}\right),$$

$${\mathrm{Co}}_f={\left(\frac{\Delta t}{|V|}\sum _f\max \left(F_f,0\right)\right)}_{f,\text{UD}},r=\frac{{\varphi }_C-{\varphi }_U}{{\varphi }_D-{\varphi }_C},$$

$${\varphi }_{U^{\ast }}\approx {\varphi }_D-2\nabla {\varphi }_C\cdot {\mathbf{d}}_{\text{CD}}.$$

$${\varphi }_{U^{\ast }}=\text{min}\left(\text{max}\left({\varphi }_{min}^{nb},{\varphi }_{U^{\ast }}\right),{\varphi }_{max}^{nb}\right),$$

$${\varphi }_f$$

$$\stackrel{!}{=}{\varphi }_C+a_C\left({\varphi }_D-{\varphi }_C\right)$$

$$\Rightarrow a_C$$

$$0\le a_C\le \min \left(r\frac{1-{\mathrm{Co}}_f}{{\mathrm{Co}}_f},1\right)$$

$$\nabla {\varphi }_C,\nabla {\varphi }_D=\frac{1}{|V_P|}\sum _f{\varphi }_f{\mathbf{S}}_f\text{and}\nabla {\varphi }_{f^{\prime }}=f\left(\nabla {\varphi }_C,\nabla {\varphi }_D\right).$$

$${\mathbf{S}}_f\cdot {\left(\nabla \varphi \right)}_f={\mathbf{S}}_f\cdot \left[{\left(\nabla \varphi \right)}_{f^{\prime }}+\nabla {\left(\nabla \varphi \right)}_{f^{\prime }}\cdot \mathbf{m}\right].$$

$${\mathbf{S}}_f\cdot {\left(\nabla \varphi \right)}_f$$

$${\left({\Gamma }_{\varphi }\right)}_f{\mathbf{S}}_f\cdot {\left(\nabla \varphi \right)}_f\stackrel{!}{=}{\left({\Gamma }_{\varphi }\right)}_{f,m}|{\mathbf{S}}_f|\frac{{\varphi }_N-{\varphi }_P}{|\mathbf{d}|}.$$

$${\left({\Gamma }_{\varphi }\right)}_{f,m}=\frac{{\left({\Gamma }_{\varphi }\right)}_f|\mathbf{d}|}{|{\mathbf{S}}_f|}\left[\frac{|\mathbf{\Delta }|}{|\mathbf{d}|}+\frac{\mathbf{k}\cdot {\left(\nabla \varphi \right)}_{f^{\prime }}}{{\varphi }_N-{\varphi }_P}+\frac{{\mathbf{S}}_f\cdot \left[\nabla {\left(\nabla \varphi \right)}_{f^{\prime }}\cdot \mathbf{m}\right]}{{\varphi }_N-{\varphi }_P}\right].$$

$${\partial }_t\alpha +\nabla \cdot \left(\alpha \mathbf{u}\right)=0$$

$${\partial }_{t}c+\nabla \cdot \left(c\mathbf{u}\right)=\nabla \cdot \left(⟨\text{D}{⟩}_{\text{h}}\nabla c\right)+\nabla \cdot \left(⟨\text{D}{⟩}_{\text{h}}\text{K}c\nabla \alpha \right),$$

$$\text{with K}=\frac{H-1}{1+\alpha \left(H-1\right)}\text{ and }H\text{ the Henry coefficient}.$$

$${\alpha }_C^n={\alpha }_C^o-\frac{\Delta t^n}{|V|}\sum _f\frac{1}{2}\left({\alpha }_f^o+{\alpha }_f^n\right)F_f^o,\text{ with }{F}_f^o={\mathbf{u}}_f^o\cdot {\mathbf{S}}_f.$$

$$\sum _f⟨\text{D}{⟩}_{\text{h,f}}^n\left[{\mathbf{S}}_f\cdot {\left(\nabla c\right)}_f^n\right]-\sum _f⟨\text{D}{⟩}_{\text{h,f}}^n{\left(\text{K}c\right)}_f^n\left[{\mathbf{S}}_f\cdot {\left(\nabla \alpha \right)}_f^n\right].$$

Full text

This work is devoted to a boundedness-preserving strategy to correct for skewness errors arising from discretisation of advection and diffusion terms within the context of interfacial heat and mass transfer based on the Volume-of-Fluid methodology. The implementation has been accomplished using a second-order finite volume method with support for unstructured meshes of general topology. We examine and advance suitable corrections for the finite volume discretisation of a consistent single-field model, where both accurate and bounded transport due to diffusion and advection is crucial. In order to ensure consistency of both the volume fraction and the species concentration transport, i.e. to avoid artificial heat or species transfer, corrections are studied for both cases. The cross interfacial jump and adjacent sharp gradients of species concentration render the correction for skewness-induced diffusion and advection errors additionally demanding and has not so far been addressed in the literature.

1. Introduction

Algebraic Volume-of-Fluid (VoF) methods [1, 2, 3, 4] are widely used for numerical simulations of flows consisting of two incompressible fluid phases with a sharp deformable interface separating them. Recently a model for interfacial species transfer within the algebraic VoF framework, namely the Continuous Species Transfer (CST) model, was presented in [5, 6]. Within this modeling framework, heat transfer is to be seen as a special case without jump and can thus be described by a simplified CST model equation. Both methods, VoF and CST, are derived by applying the Conditional Volume-Averaging (CVA) technique to the local instantaneous sharp interface model equations and the corresponding interfacial jump conditions so as to transform them into volume-averaged single-field equations. In this form the governing equations are readily suitable for consistent discretisation using a Finite Volume Method (FVM).

When considering computational domains entailing complex geometries, unstructured meshes of general topology are required. Then, depending on the mesh quality, discretisation errors related to mesh skewness arise and can severely detoriate the numerical fidelity of the method resulting in loss of boundedness, accuracy and order of convergence. Typically, two types of skewness errors are to be distinguished: errors due to non-conjunctionality and errors due to non-orthogonality at the three-point stencil around a cell face (see Figure 1). In order to avoid detoriation of numerical fidelity on such meshes, skewness correction approaches have been devised for FV discretisation and are nowadays well established practice. Classical approaches are explicit and correct for non-orthogonality during discretisation of diffusion terms and non-conjunctionaly during FV discretisation of advection terms, respectively. For an overview of state-of-the-art classical approaches to skewness correction the interested reader is referred to [7, 8, 9]. A non-classical explicit approach is the ghost-point method set out by Perić in [8], correcting for both non-conjunctionality and non-orthogonality simultaneously.

Applying classical (explicit) skewness correction approaches to the VoF advection equation on distorted meshes results directly in severe unboundedness of the volume fraction and as a consequence in substantial problems with stability and/or accuracy. The reason for this resides in the fact that due to the typically high density ratios even small errors in the volume fraction cause large errors in mass fraction and consequently in the momentum. It has been found that in such cases implicit correction approaches are needed, see [10, 11, 12] and more recently [13, 14]. Taking interfacial heat and mass transfer into consideration renders the situation even more challenging, since there are now both the discretised advection and diffusion terms to be skewness corrected. While it is obvious that classical explicit approaches remain unsuitable (likewise causing unboundedness errors as for the VoF transport), up to now it is unclear whether implicit correction approaches can be used to correct for mesh skewness when applied to the discretised form of the CST model for interfacial heat and mass transfer.

This contribution is concerned with enhancements of implicit correction approaches based on a detailed investigation into their suitability to maintain numerical fidelity such as boundedness, accuracy and order of convergence in the context of VoF/CST-based simulations of interfacial heat and mass transfer on distorted meshes of general topology. In contrast to previous work, we introduce implicit skewness corrections for finite volume discretisations of both advection and diffusion terms. Moreover, we use finite volume meshes which have been distorted in a targeted manner rather than randomly (as is common practice in the existing literature) in order to allow a systematic study on mesh-induced errors and avoid error cancellation. The objectives of this work are two-fold, viz.

• •

to devise a suitable method for mesh skewness correction for the discretised advection and diffusion terms of the single-field CST model for VoF-based interfacial heat and mass transfer. Notably, one of the CST model terms is discretised in a non-standard way [6], that is, despite originating from a diffusive transport term, it has been found that it is best to treat one term implicitly by means of a FV divergence rather than a FV laplacian operator.

• •

to identify and clarify the influence of a jump in the scalar transport variable on the order of convergence when refining a systematically distorted finite volume mesh and applying implicit skewness corrections.

This study has been accomplished using and extending OpenFOAM (version 4.x).

2. Standard finite volume discretisation and Mesh-Induced Errors

In order to better understand the proposed approach to implicit skewness correction for VoF/CST-based interfacial heat and mass transfer on distorted meshes of general topology, the errors on skewed meshes in finite volume discretisation of a prototypical scalar advection-diffusion equation are briefly recapped and a review of prominent publications addressing these errors is given.

The starting point of our investigations is the advection-diffusion equation for an arbitrary scalar quantity $\phi$, i.e.

[TABLE]

In the finite volume context, the standard discretisation procedure is based on the application of Gauss’ rule, transferring the divergence terms into face integrals, and on the mid-point rule. This yields a formal second-order accurate approximation. In collocated finite volume methods, as employed in this work, the values of quantity $\phi$ are stored in the cell centre locations, which necessitates the interpolation of cell-centred data onto the face centres. Below cell-face interpolation will be denoted by $\left(\cdot\right)_{{\!f}}$. In order to keep the discretisation second order accurate, the interpolation schemes have to be at least of the same order. However, in general, even if higher order interpolation schemes are selected, the convergence order as well as the overall accuracy can still deteriorate, particularly in cases of complex computational domain geometries and meshes of general topology (polyhedral meshes). The errors stem directly from local mesh skewness, namely non-orthogonality and non-conjunctionality. They can be corrected during discretisation of the respective transport terms in (1).

2.1. Advection term

The standard finite volume discretisation of advective terms reads

[TABLE]

with the volumetric face flux $F_{{\!f}}:=\mathbf{S}_{{\!f}}\cdot\mathbf{u}_{{\!f}}$. The face-centred values of $\phi$ are obtained by using interpolation schemes. In general, second order accuracy of the discretisation formally relies on the assumption that the face-centre position corresponds to the intersection point of the face and the connecting line between the neighbouring cell centres. If this assumption does not hold, as the constellation depicted in Figure 1(a) shows by way of example, so called non-conjunctionality errors occur.

The standard approach to account for this error is to correct the face-interpolated value $\phi_{{\!f}^{\prime}}=f(\phi_{{\!N}},\phi_{{\!P}})$ by an explicit term reading

[TABLE]

where $\mathbf{m}$ denotes the vector from the face intersection point to the face centre.

The above time-explicit correction of the mesh-induced non-conjunctionality error is very simple and computationally inexpensive but may lead to unboundedness of the transported quantity [7], e.g. when high-resolution schemes (based on TVD, NVD face limiters etc.) are utilised for the discretisation as is common practice in FVM. Magalhaes et al. [15] introduce a similar approach for adaptive hexahedral meshes to account for non-conjunctionality errors, where the face value is reconstructed from the face-neighbouring cell-centred values and the respective gradients, and then averaged:

[TABLE]

An entirely different correction approach, the so-called ghost point-based interpolation scheme, which corrects the face interpolation by introducing a virtual non-skew stencil, was proposed by Ferziger and Peric [8] and successfully employed e.g. in [16]. However, if boundedness of the transported quantity is of utmost importance, as e.g. in the VoF-based transport of the volumetric phase fraction or species concentration fields, the above explicit methods are not usable as they do not ensure this property.

Therefore, this work makes use of a modified implicit version of the standard correction approach described by Equation (3), which makes use of an improved method for preserving the boundedness of the transported quantity and is described in detail in Section 3.

2.2. Diffusive term

Analogously to Equation (2), the discretised form of the diffusive term is given by

[TABLE]

The term $\left(\Gamma_{\!\phi}\right)_{\!f}\mathbf{S}_{\!f}\cdot\left({\nabla\!}\phi\right)_{\!f}$ is preferably to be discretised in a compact way, that is only information of the two cells neighbouring the respective face is considered:

[TABLE]

This approximation of the surface normal gradient $\nabla_{\!f}^{\perp}\phi$, however, only holds true for orthogonal meshes. As a consequence, a numerical error is introduced in case of local non-orthogonality, i.e. if $\mathbf{S}_{\!f}\nparallel\mathbf{d}$ (see Figure 1(b)).

One way to account for this error was introduced in [17] and [18], where cell vertex-interpolated values, obtained by either inverse distance weighting or least-squares interpolation, are utilized to correct for mesh non-orthogonality, which in effect significantly increases the numerical stencil. This, however, counteracts the accurate discretization of discontinuities, which are inherently smoothed by a large discretization stencil.

A different approach, proposed by [19, 20, 21], is to split up the vector $\mathbf{S}_{\!f}$ in such a way that

[TABLE]

with the restriction of $\mathbf{\Delta}$ being parallel to $\mathbf{d}$. Employing the above splitting, the discretisation can be performed as

[TABLE]

For the computation of the magnitude of $\mathbf{\Delta}$ and $\mathbf{k}$ different options are available (Jasak [7], Demirdzic [22]). For further reading regarding non-orthogonal corrections the interested reader is also referred to the works of Demirdzic [19, 20, 22], as well as Mathur and Murthy [21]. These pionieering works represent the main contributions to non-orthogonality corrections and their development history. An overview of alternative but similar correction approaches can also be found in [23]. All of the above mentioned approaches eventually result in a formulation similar to or exactly the same as that given in (8).

If in addition the mesh exhibits local non-conjunctionality, i.e. is generally skewed (see Figure 2), further corrections may have to be applied. One available explicit method for correcting both skewness errors simultaneously is the ghost-point interpolation proposed by [8]. Generally, non-conjunctionality corrections in the discretisation of diffusion terms are simply omitted without comment. In this work, however, their influence on accuracy will be investigated.

3. Proposed Bounded Implicit Correction of Mesh-Induced Errors

As the explicit correction of mesh-induced errors can lead to unboundedness, suitable implicit corrections must be devised. For this purpose, we first briefly recap boundedness-preserving criteria from the literature and then adapt the formulation of the skewness correction in a way such that they can be applied.

3.1. Boundedness Criteria of Advection Term on Unstructured Meshes

To enforce boundedness during discretisation of the advection term, different boundedness criteria based on positivity [24], are readily available in literature. Most prominent (inter alia) are the Total Variation Diminishing (TVD) [25] and Local Extrema Diminishing (LED) [26] criteria as well as the transient Convection Boundedness Criterion (CBC) [27].

Defining an interpolation scheme based on a flux limiter formulation as [28, 25]

[TABLE]

with the mesh weight $\delta_{\!f}=\frac{\mathbf{d}_{\!\text{Cf}}\cdot\mathbf{S}_{\!f}}{\mathbf{d}_{\!\text{CD}}\cdot\mathbf{S}_{\!f}}$ and the limiter function $\Psi\!\left(r\right)$, the TVD conditions can be transferred into

[TABLE]

with $\operatorname{Co}_{\!f}$ denoting the face Courant number and $r$ the gradient ratio, being defined as

[TABLE]

where $F_{\!f}$ denotes the volumetric face flux and $\left(\cdot\right)_{{\!f},\text{UD}}$ denotes upwind interpolation. Contrary to the condition Sweby [28] derived for a flux limiter function $\Psi$, the condition in Equation (10) does not require the mesh to be uniform. Additionally, the above limiter condition makes use of the whole boundedness region, while the original limiter condition introduced in [28] only offers a very conservative estimation of the boundedness region.

On unstructured meshes of general topology the computation of $r$ leads to the problem that $\phi_{\!U}$ is required, the definition of the upwind node, however, is not clear.

In literature, a variety of different possible remedies have been introduced: [10] chooses the upwind cell from all face-neighbouring cells of the centre node as the one whose centre position is closest to the line through centre and downwind nodes (see Figure 4(a)), whereas [29] presents alternatives based on an extended stencil (see Figure 4(b)) and compact stencil (see Figure 5(a)). The latter is the most widely utilized approach and introduces a virtual upwind node as the backward projection of the downwind to centre node.

The value in the virtual upwind node is then computed by a second-order Taylor series approximation (see Figure 5(b))

[TABLE]

On uniform hexahedral meshes the exact value for the upwind cell value is recovered, if the gradient in Equation (12) is calculated using Gauss’ approximation. Different strategies for calculating the $r$-factor and/or virtual upwind value reported in literature can be found in [30, 31, 32, 29, 33]. To maintain a bounded solution, different limiting strategies for the virtual upwind have been proposed in literature [31, 32]. The simplest strategy is to bound the virtual upwind value by the physical bounds of the transported quantity (e.g. between 0 and 1 for the volumetric phase fraction). [31] introduced a local limiting strategy based on the values in the neighbouring (nb) cells

[TABLE]

which on unstructured meshes requires an efficient search algorithm to find the computational cell containing the virtual upwind node. This approach is adopted in the presented work, employing the local vicinity search algorithm of [34, 35].

3.2. Implicit Mesh-induced Error Correction for Advection

In order to apply the previously described boundedness criteria to advection schemes accounting for non-conjunctionality errors, the formulation of the correction (Equation 3) has to be rearranged. To this end, an approach has recently been introduced by Denner and vanWachem [13, 14], which is based upon including the explicit correction into the interpolation weights. This is done by extending the general interpolation scheme already defined in Equation (9) by the correction term and solving the resulting equation for a newly defined interpolation weight $a_{\!C}$:

[TABLE]

To obtain a bounded solution the TVD condition (Equation 10) now no longer has to be enforced on the flux limiter $\Psi$ of the base scheme, but on the resulting interpolation weight $a_{\!C}$:

[TABLE]

Thus, in this approach the correction is still explicit but is incorporated directly into the interpolation weights which allows for a bounded and formally time-implicit discretisation.

Another problem which must be addressed when aiming at boundedness-preserving non-conjunctionality corrections, is the coupling between face-interpolation and gradient computation, which has been noted by numerous researchers [10, 11, 12]. A face interpolation with non-conjunctionality correction (cf. Equation 14) depends on the face gradient that is interpolated from gradients in the face-neighbouring cell centres, which in return may require the face interpolated values to be known (e.g. using Gauss’ theorem)

[TABLE]

To overcome this problem, [11] and [10] introduce different iterative methods, requiring the computation of the face interpolated value and face gradient multiple times in a correction loop. [12] introduces an improved approach, based on an implicit calculation of the cell-centred gradient using a least-squares method.

It is shown in this work (Section 5.1), that the method utilized to compute the gradient in the correction term has a major influence on the accuracy of the non-conjunctionality correction for advection terms. In this work, the least-squares gradient computation introduced in [12] and the gradient computation using standard Gauss’ theorem are utilized and compared in respect of the accuracy of non-conjunctionality corrections. Both gradient computations are explicit without any correction loops.

3.3. Implicit Mesh-induced Error Correction for Diffusion

A straight forward way to account for both mesh-induced errors simultaneously is to consecutively apply the non-conjunctionality and non-orthogonality corrections already described. In a first step the non-conjunctional corrected gradient of the transported variable $\phi$ analogous to Equation (3) is obtained:

[TABLE]

Using Equation (8) to discretise $\mathbf{S}_{\!f}\cdot\left({\nabla\!}\phi\right)_{{\!f}^{\prime}}$ finally yields:

[TABLE]

Analogously to the approach previously introduced for advection terms, the above discretisation of the diffusive term can also be rearranged in a way that allows for implicit bounding. This is done by rearranging Equation (17) to formally correspond to the uncorrected discretisation of the diffusive term (cf. Equation 6):

[TABLE]

Comparing Equation (18) with Equation (6) reveals that the diffusive coefficient $\left(\Gamma_{\!\phi}\right)_{{\!f},m}$ has to be defined as

[TABLE]

Theoretically this formulation now enables the implicit application of limiting criteria to the diffusive coefficient $\left(\Gamma_{\!\phi}\right)_{{\!f},m}$, which in turn contains the explicit correction. However, until now it remains uncertain as to whether limiting is needed at all and if so, the criterion required is missing. Both issues will be addressed in the results section.

4. Utilized Finite Volume Discretisation

The developed correction strategies are applied to the discretised phase fraction advection equation, i.e. the VoF transport eqn.

[TABLE]

and the discretised species transport equation underlying the CST model [6]. Employing a harmonic mean diffusion coefficient, the simplest version of the latter can be written as (cf. [6])

[TABLE]

The above equations can be discretised in many different ways to account for mesh-induced errors. Suitable candidates – based on the selection in the previous Sections – are identified in the remainder of this Section and tested in Section 5. The discretised form of the phase fraction transport equation, using time-centered Crank-Nicolson time discretisation, reads

[TABLE]

The most crucial decision in discretising the advection term is how to compute the face-interpolated value $\alpha_{{\!f}}$. This work considers four different schemes, namely the Upwind UDS, Central Differences (CDS), CICSAM [1] and a specific local blended scheme (LB). The employed LB scheme is similar to the CICSAM scheme in that it uses the CICSAM scheme’s blending function to blend between the central differencing and downwind scheme, resulting in a formally unbounded interpolation scheme. In the LB case, boundedness will be enforced by implicit correction once, while in the CICSAM case, boundedness is enforced twice: by the non-linear limiter of the scheme and by limiting the skewness corrected face interpolation.

Each of the schemes is tested without correction (UC), with standard explicit correction (EC) and with the newly proposed correction, from now on referred to as Semi-Implicit Skewness Correction (SISC) scheme. Both corrections make use of the gradient operator which is discretised by both, Gaussian Gradient (GG) and Least Squares Fit (LSF). An overview of all interpolation schemes, correction approaches and gradient calculation methods varied against each other for discretisation of the advection term is given in Table 1.

In this work, we test different gradient schemes for calculation of the correction part and utilize the least-squares gradient to compute the virtual upwind value on non-uniform meshes. It is emphasised that on uniform Cartesian meshes the Gaussian gradient should be utilized for the virtual upwind calculation to maintain accuracy.

It is crucial to discretise the left hand side of the species transport equation (21) identical to the phase fraction transport equation in order to avoid artificial species transfer, i.e. identical discretisation schemes for the respective terms and corrections need to be employed. A suitable discretisation of the right hand side as utilized in the presented work reads

[TABLE]

Since the phase fraction transport equation is solved first, the phase fraction field $\alpha^{n}$ at the new time level is already explicitly available, allowing for the fully time implicit discretisation shown above. The gradient in face normal direction in both terms need to be discretised identically, whereby in the scope of this work three different discretisation methods are tested. Firstly, without any correction, which therefore corresponds to Equation (6). Secondly, by applying the explicit non-orthogonal correction as proposed by Jasak [7], cf. Equation (8), and thirdly, by additionally accounting for non-conjunctionality errors as formulated in (19). These three variants are from now on referred to as UC, NO and NO/NC, respectively.

The diffusion coefficient $\langle\text{D}\rangle_{\!\text{h,f}}^{n}$ and the term $\left(\text{K}c\right)_{f}$ are linearly interpolated, once with and once without non-conjunctionality correction. Table 2 summarises all used discretisations for the different terms of the CST model.

5. Results and Discussion

To evaluate the performance of all discretisation methods presented in Section 4, including the newly developed correction strategies SISC and NO/NC, they are verified using a series of test cases involving steep gradients and discontinuities. Two aspects are considered to judge the discretisation performance: the obtained solution should accurately approximate the analytical solution and at the same time boundedness must be preserved (as required for both phase fraction and species concentration).

The evaluation comprises of two steps. First, the discretisation of the phase fraction transport equation is varied according to Table 1 and verified. Then secondly, the diffusive terms in the CST-model are discretised with different variants as given in Table 2. In both steps,the respective discretization performance is verified via different test cases.

In all test cases 2D meshes are used. To provoke mesh induced errors the mesh is systematically distorted in such a way that non-orthogonality as well as non-conjunctionality errors are present. Figure 6 shows the mesh structure which is identical for all cases.

5.1. Correction of VoF-based phase fraction advection term

5.1.1. Test Cases

First we consider a simple 2D plug flow, whereby the transported field is advected from left to right by a constant velocity field (see Figure 7(a)). The second case examines the shape-preserving properties of different advection schemes when a circular shape is translated through a distorted mesh over a distance equal to 12 times its diameter (see Figure 7(b)). Here, the utilised average mesh resolution is approximately 15 cells per diameter of the circle. The numerical time step in both simulations is set in such a way that the Courant number limit $\text{Co}=|\mathbf{u}|\Delta t/\Delta x\leq 0.1$ is maintained.

5.1.2. Results

Results are visualized by plotting the phase fraction profile over the path ‘L1’ and ‘L2’ respectively (see Figure 7). Without any correction all interpolation schemes show poor results in both test cases. As might be expected the usage of CDS or LB generates strongly unbounded, inaccurate solutions in all cases. While boundedness is ensured by the use of UDS or CICSAM interpolation, the accuracy is insufficient. More specifically, high diffusion is observed when applying UDS and severe deformation is ovserved with CICSAM. The respective plots are shown in Figure 8, where it should be noted that some lines cannot be detected as the unboundedness of the solution is too severe.

Adding explicit correction worsens the results even further. Except for the UDS scheme all generated solutions are unbounded and all of them suffer heavily from a lack of accuracy. This is the case irrespective of whether the gradient is calculated using GG or LSF, where applying the latter still produces the better results.

Now, evaluating the SISC scheme, we can state that for all tested combinations boundedness is preserved up to $10^{-6}$ (see Figure 9). Another interesting finding is that the LSF of the gradient outperforms the Gaussian gradient in all cases. Except for the UDS scheme, where the strong diffusion overrides all other effects, the difference is quite significant. Using the SISC-bounded CDS, diffusion has a huge effect and smears out the initial step profiles. Applying CICSAM yields a significant improvement because it is better able to preserve steep gradients. The best result, however, is achieved by the LB scheme, where almost the exact shape is preserved when using LSF for gradient computation.

To finally assess all relevant scheme combinations the phase fraction field is also visualized (see Figures 10 and 11). Set-ups not considered as relevant are those which lead to unbounded results and those using the Gaussian gradient approximation because LSF performs equally or better in all cases. Other than the results gained with the SISC scheme this leaves only the non-corrected UDS and CICSAM schemes, and for the explicitly corrected set-ups solely the UDS scheme.

Figure 10 depicts the fields excluding SISC. The UDS scheme without and also with explicit correction is dominated by diffusive effects which leads to a smeared out $\alpha$-distribution. CICSAM results in a slightly sharper solution, which, however, still can not be considered as satisfactory. It can therefore be concluded that none of the commonly used advection schemes meet the requirements concerning accuracy and boundedness.

Looking at the distribution of the phase fraction obtained with the SISC scheme (depicted in Figure 11) allows a more extensive evaluation. The UDS scheme clearly shows a diffusive effect which is so pronounced that neither the shape of the plug flow nor the circular shape are preserved. To a lesser degree, this effect is also visible when CDS is used. CICSAM and the LB scheme show good results in both test cases. However, although CICSAM preserves sharp gradients quite well, the LB scheme is able to keep the initial steepness of the profile. Addressing the ability to preserve the circular shape it can be seen that both schemes tend to transform the circle into a quadratic shape, for the mesh topology under consideration. Hence, for algebraic VoF methods the mesh topology has a major influence on accuracy (see Figure 12). Flux/cell-face alignment is important for the accuracy of algebraic VoF methods and thus polyhedral mesh topologies are to be favoured over to quadrilateral (or hexahedral) meshes. Numerical simulations using standard discretisation on polyhedral meshes, however, suffer from mesh-induced skewness errors and hence the present work becomes generally important for algebraic VoF methods.

To conclude it can be stated that the proposed correction scheme, in contrast with all other tested schemes, is able to keep the boundedness of the solution and at the same time, if combined with a suitable interpolation scheme, is capable of producing accurate results.

5.2. Correction of CST-based species transfer term

Different variants for discretising the terms of the CST-model are tested, with special attention devoted to mesh induced error correction (cf. Table 2). The three main questions are:

• •

How severe is the impact of skewed meshes on the accuracy of the CST-model and which of the proposed correction strategies is needed to restore the expected order of accuracy on non-distorted meshes?

• •

Does explicit correction lead to unboundedness of the solution? In other words, is implicit limiting necessary?

• •

Is it sufficient to account for non-orthogonality (NO) only, or must also non-conjunctional error correction (NO/NC) be taken into account when discretising a gradient in face normal direction?

Even though the last question is commonly answered in the negative it shall nevertheless be discussed here, since species concentration fields commonly exhibit large discontinuities at the interface in both value and gradient. Consequently, errors introduced by skewed meshes are more pronounced and may necessitate a full skewness correction to obtain accurate results.

5.2.1. Test Case

In order to find the scheme which best meets the requirements, the simulation results of a planar diffusion case are evaluated. To avoid effects of the convective term the velocity field is set to zero. To evaluate the different correction strategies, species transfer from a surrounding stagnant gas into a stagnant liquid film is considered (see Figure 13). Five different Henry coefficients ($\operatorname{H}=0.033,0.2,1,5,30$) are tested, where the diffusion coefficient in the gas phase is $D_{g}=1\cdot 10^{-1}\,\text{m}^{2}/\text{s}$ and in the liquid phase $D_{l}=1\cdot 10^{-5}\,\text{m}^{2}/\text{s}$. The initial gas concentration is $c=1$ and the liquid phase is initialized with $c=0$.

In order to obtain the time-dependent numerical exact reference solution, the diffusion equation is discretised in 1D using the Finite Difference Method (FDM). Here the same numerical set-up (domain size, boundary and initial conditions) is employed and the domain is sufficiently resolved. The solution which is obtained by implicitly solving the system of equations with Octave version 3.2.4 is referred to later as Exact.

5.2.2. Results

Applying a specific correction when discretising the face-interpolated values (Table 2(b)), the results in some cases slightly improve, in others, however, deteriorate. A possible explanation is that to effectively correct the calculation of the face-interpolated values within a steep profile, a small gradient stencil and, at the same time, an accurate gradient computation are required. To maintain a small stencil, we utilize the Gaussian gradient for correction of the face-interpolation $\left(\text{K}c\right)_{f}$. However, as severely distorted meshes are used, the accuracy of the approximation suffers which is probably the cause for the observed poor performance. Consequently in the following only the variation of the gradient discretisation (Table 2(a)) is discussed while the face-interpolated values are discretised by simply using the uncorrected CDS.

Results are visualized by plotting the profile of the concentration $c$ over the path ‘L3’ (see Figure 13). As shown in Figure 14, this is done for all considered Henry coefficients $\operatorname{H}$ and correction strategies. Despite the heavily distorted mesh, all approaches lead to reasonable results. For all Henry coefficients, the uncorrected approach tends to overestimate the diffusion and therefore does not correspond to the exact result. However, NO correction and NO/NC correction, result in a good approximation to the exact solution indicating that both approaches effectively correct the mesh-skewness induced errors (see Figure 15).

Enlarged sections of two representative concentration profiles are shwon in Figure 15. It is shown that even for the lowest and highest Henry coefficients, no unboundedness occured, which indicates that no limiting of the correction terms is required. However, one might intuitively try limiting the diffusive coefficient in the NO/NC formulation (cf. Equation 19), which is also tested in the present work. Due to the lack of a mathematically sound limiter criterion, the most intuitive one is applied, which is to restrict the modified diffusion coefficient to positive values. When applied, the results deteriorated and the solution showed flickering in time. Based on this finding and the fact that the limiting does not seem to be crucial no further investigations were carried out and all the results presented below do not include any limiting.

To conclusively assess the convergence performance of the different schemes a mesh convergence study was conducted. The average molar concentration in the gas phase after $t=0.1$ s and $t=0.5$ s was chosen as error measure . The latter is almost constant as the diffusive coefficient in the gas phase exceeds that in the liquid phase by orders of magnitude. By this error definition the focus of the evaluation is not on the representation of the exact concentration profile (which we have already investigated, see Figure 14 and 15) but on the correct prediction of the time integrated species flux over the interface. Figures 16 and 17 show the $L_{1}$ error plotted over the average cell size for each Henry coefficient at different times $t$. It can be seen that the skewness-correction approaches employed generally lead to a substantial improvement of accuracy and convergence order. In all cases at least first order, for $\operatorname{H}=1$ (heat transfer) even second order convergence is achieved. In contrast the uncorrected discretization converges with a much lower order or exhibits no convergence at all.

The observed convergence characterisic in the corrected case corresponds exactly to the formal order of the CST formulation, which is second order for $\operatorname{H}=1$ and tends to first order for $\operatorname{H}\ll 1$ and $\operatorname{H}\gg 1$ even on uniform cartesian meshes. It therefore can be stated that applying the corrections restores the formal order of convergence on skewed meshes.

In the uncorrected case if $\operatorname{H}\leq 1$ and in the corrected case if $\operatorname{H}=30$, unusual convergence behaviour is detected, that is the plot shows a local minimum of error. This can be explained by the fact that in these cases the concentration profile crosses the exact one when beeing evaluated using different mesh resolutions. In consequence a minimum error is found at some random intermediate mesh resolution and seems to converge for lower and to diverge for higher resolutions. In addition, these cases obviously converge to a wrong solution. The magnitude of the absolute error thus corresponds to the $L_{1}$ error value it converges to for further refined meshes.

By considering both the convergence study and also the concentration profiles depicted in Figure 14, the question about the impact of non-conjunctionality errors when discretising diffusive terms can finally be answered: it is indeed negligible, since neither the NO nor the NO/NC profiles differ significantly, even for cases where sharp gradients are obviously present.

Considering the heavily distorted mesh, the achieved solution accuracy might be surprising. Hence, we also raise the question as to whether errors may possibly be cancelled out because of the use of a systematically distorted mesh (see Figure 6). To eliminate this possibility all test cases are additionally solved on a randomly distorted hexahedral and tetrahedral mesh, both depicted in Figure 18, as these are often used in literature.

Figures 19 and 20 show the results analogous to the already discussed set of graphs in Figure 14. In both cases, independent of the selected discretisation method or Henry coefficient, the results almost perfectly match the solution, which leads to the conclusion that randomly distorted meshes (unlike the systematically distorted meshes studied here) result in error cancelation and are thus inappropriate for the systematic investigation of skewness correction approaches.

6. Summary and Conclusions

This work is aimed at the targeted investigation and enhancement of mesh-skewness correction approaches with special focus on implicit and bounded corrections in the context of VoF-based interfacial heat and mass transfer. The presence of a jump and/or steep gradients substantiate a sound reconsideration of the skewness correction approaches, which are commonly recognised as accepted knowledge in the literature. Discretisation procedures for the advection and diffusion terms are verified separately using different test cases on systematically distorted meshes. Randomly distorted meshes are shown to lead to error cancelation and thus are proven to be inappropriate for a systematic study of mesh-induced errors.

For the phase-fraction advection term we investigate the boundedness and accuracy achieved with different discretisation strategies, i.e. the shape-preservation and sharpness properties of different approaches to skewness correction. The main finding is that implicit correction using the novel semi-implicit skewness correction (SISC) approach is imperative. Further, it is shown that the gradient scheme used for computation of the correction term has a strong influence on accuracy. In our study, the performance of the least square gradient computation is found superior over Gauss gradient approximation.

For the diffusive heat/species transfer terms in the CST method, again we investigate the boundedness and accuracy achieved with different discretisation strategies. Additionally, focus is on convergence on strongly distorted meshes. It is shown that the approach for semi-implicit skewness correction (SISC) used for the advection term can be adopted for diffusion terms as well. However, for the test case studied here applying a relatively small time step, boundedness of the diffusion term is also attained using explicit correction approaches, although not a-priori guaranteed as in the implicit approach. Non-orthogonality correction (NO) has been found to perform reasonably well even for the transport of fields with discontinuities. Simultaneuos implicit correction of non-orthogonality and non-conjunctionality errors (NO/NC) has shown no further benefit even for quantities exhibiting strong solution curvature, i.e. strong spatial changes in the gradient. Finally, the authors have demonstrated that the convergence order on non-distorted meshes can be retained even on strongly distorted meshes only when suitable correction approaches such as NO or NO/NC are used.

To conclude, it is shown that the CST-method is able to predict accurate solutions on moderatly distorted meshes even without corrections. For severe distortions the accuracy can still be maintained by applying suitable, i.e. bounded implicit, mesh-skewness corrections.