Interface Sharpening Using Temporal Subcycling

Applying temporal sub-cycling to the transport of volume fraction can improve the resolution of the interface between two phases.

The Volume-of-Fluid (VOF) method implemented in Simcenter STAR-CCM+ employs the High-Resolution Interface Capturing (HRIC) scheme to give a second-order accurate approximation of volume fraction $α$ at faces. Similarly, the Mixture Multiphase (MMP) method employs the Adaptive Interface Sharpening (ADIS) discretization scheme for the convective term of the volume fraction, which uses a high-resolution interface capturing (HRIC) scheme in the vicinity of a large scale interface.

To yield well-posed problems for affected transport equations, the HRIC scheme is bound by an upper Courant number (CFL) limit. If the CFL limit (usually < 1) is exceeded, the method falls back to the upwind scheme to sustain stability. As a first-order method, the upwind scheme results in an irreversibly smeared interface.

Introducing a temporal subcycling to the transport of volume fraction allows the poor approximation of the volume fraction at faces to be corrected locally. Hence, a sharp interface solution can be obtained as the effective (local) CFL (denoted as ${C F L}_{s u b}$ ) is reduced. Please see the High-Resolution Interface Capturing (HRIC) for further details.

The transport of volume fraction of phase $i$ is given by Eqn. (2584) or Eqn. (2875).

For standard FV discretization and time-independent control volumes, integrating Eqn. (2584) or Eqn. (2875) over an arbitrary time interval $[t_{n,} t_{n} + Δ t]$ gives:

Figure 1. EQUATION_DISPLAY

α^{n + Δ t} V - α^{n} V + \int_{t^{n}}^{t^{n} + Δ t} (\sum_{A} α_{f} (s) v \cdot d a) d s = \int_{t^{n}}^{t^{n} + Δ t} (S (α (s), ...) d s

(2622)

where

α_{f}

is the face value volume fraction,

Δ t

is the current time-step size and

t^{n}

is the current time level. The right hand-side contributions from Eqn. (2584) or Eqn. (2875) have been combined into a single source term

S (α (s), ...)

Single-Step

For the Single-Step solution strategy the solver approximates Eqn. (2584) or Eqn. (2875) as in the following relation:

Figure 2. EQUATION_DISPLAY

α^{n + Δ t} V - α^{n} V + Δ t (\sum_{A} α_{f}^{n + Δ t} (Δ t) v \cdot d a) = Δ t S (α^{n + Δ t}, ...)

(2623)

The face value $α_{f}^{n + Δ t}$ is calculated by the HRIC scheme, and treated implicitly with respect to time.

Multi-stepping resolves the time-step restrictions on the global time-step size

Δ t

by sub-dividing the time interval

[{t^{n}, t}^{n} + Δ t

into N uniform sub-intervals

[t_{i}, t_{i + 1}], i = 0, ..., N - 1

the $t_{i}$ sub-intervals are defined as $t_{i} = t^{n} + i τ$
$τ$ is the corresponding sub-stepping time-step size, defined as $τ = \frac{Δ t}{N}$
The intermediate solutions $α_{i + 1}$ are then computed for each of the sub-intervals i.

There are two multi-step solver options in Simcenter STAR-CCM+ :

Explicit Multi-Step: Available only with Volume-of-Fluid (VOF) multiphase model.
Figure 3. EQUATION_DISPLAY

$α_{i + 1} V - α_{i} V + τ (\sum_{A} α_{f, i} (τ) v \cdot d a) = τ S (α^{n + Δ t}, ...), f o r i = 0, 1, ..., N - 1$

(2624)

Implicit Multi-Step: Figure 4. EQUATION_DISPLAY

$α_{i + 1} V - α_{i} V + τ (\sum_{A} α_{f, i + 1} (τ) v \cdot d a) = τ S (α^{n + Δ t}, ...), f o r i = 0, 1, ..., N - 1$

(2625)

In multi-stepping, the calculation of the HRIC face value

α_{f}

is subject to the reduced time-step size

τ

, and thus the CFL dependent upwind-blending is reduced. The sum across all sub-steps then gives:

Figure 5. EQUATION_DISPLAY

α^{n + Δ t} V - α^{n} V + \sum_{i = 0}^{N - 1} τ (\sum_{A} α_{f, i (i + 1)} (τ) v \cdot d a) = Δ t S (α^{n + Δ t}, ...)

(2626)

The following table provides a brief summary of the first-order time integrators in Simcenter STAR-CCM+:


	Single-Step	Implicit Multi-Step	Explicit Multi-Step (VOF only)
$\int_{t^{n}}^{t^{n} + Δ t} (\sum_{A} α_{f} (s) v \cdot d a) d s$	Figure 6. EQUATION DISPLAY $Δ t (\sum_{A} α_{f}^{n + Δ t} (Δ t) v \cdot d a)$ (2627)	Figure 7. EQUATION DISPLAY $\sum_{i = 0}^{N_{i m p - 1}} τ (\sum_{A} α_{f, i + 1} (τ) v \cdot d a)$ (2628) where $N_{i m p}$ is the number of implicit sub-steps	Figure 8. EQUATION DISPLAY $\sum_{i = 0}^{N_{e x p - 1}} τ (\sum_{A} α_{f, i} (τ) v \cdot d a)$ (2629) where $N_{e x p}$ is the number of explicit sub-steps
Stability	Unconditionally Stable	Unconditionally Stable	Conditionally Stable: $\frac{C F L}{N_{\exp}} \leq 1$
Accuracy (sharp interface)	$C F L \leq 0.5$	$\frac{C F L}{N_{i m p}} \leq 0.5$	$\frac{C F L}{N_{\exp}} \leq 0.5$
Computational Cost	Low	$N_{i m p}$ determined by user input	$N_{\exp}$ determined by stability constraint

The explicit integration imposes a CFL condition where

\frac{C F L}{N_{\exp}} \leq 1

and the VOF Explicit Multi-Step solver determines the number of required explicit sub-steps

N_{\exp}

automatically.

Implicit multi-stepping for MMP and VOF is not bound to a time-step size restriction, and the number of implicit sub-steps $N_{i m p}$ is chosen by the user. If the reduced time-step size $τ$ exceeds the HRIC scheme's CFL limitations, the volume fraction solution may become increasingly diffusive.

At each inner iteration, multi-stepping uses the volume fraction solution of the previous time-step as initial data. After sub-stepping, the VOF Explicit Multi-Step Solver applies under-relaxed corrections only to the final solution at ${t = t}^{n} + Δ t$ of the previous inner iteration, whereas the Implicit Multi-Step Solver applies under-relaxed corrections to all intermediate solutions.

Implicit multi-stepping solves the transport equation Eqn. (2625) at each sub-interval, therefore the residual

r_{i}

of this equation is calculated for each sub-interval. Simcenter STAR-CCM+ defines the residual

r

for the entire time-step as the average of these sub-interval residuals:

Figure 9. EQUATION_DISPLAY

r = \frac{{\sum r_{i}}_{i = 0}^{N - 1}}{N}

(2630)

Automatic Determination of Sub-step Size in the VOF Explicit Multi-Step (Deprecated) Solver

The process for determining the sub-step size is as follows:

1. Identify the cells that belong to the VOF interface

All cells within which the magnitude of the gradient of any of the VOF phases exceeds a minimum threshold value are considered to belong to the VOF interface:

Figure 10. EQUATION_DISPLAY

\nabla α_{i} \geq \nabla α_{\min}

(2631)

The minimum threshold for the magnitude of the volume fraction gradient is equal to:

Figure 11. EQUATION_DISPLAY

\begin{matrix} \nabla α_{\min} = \frac{Δ α_{\max}}{l_{r e f}} = \frac{1}{f_{c} V^{\frac{1}{3}}} & f_{c} : = 20 \end{matrix}

(2632)

where:

$Δ α_{\max}$ is the maximum possible change in volume fraction (which is 1.0)
$l_{r e f}$ is the cell reference length scale (the length scale of the cell related to its volume)
$V$ is the cell volume
$f_{c}$ is an arbitrarily chosen and hard-coded constant (20) that is used in defining the free surface.

Cells that have a volume fraction gradient at least equal to the maximum possible change in volume fraction over a distance equal to $f_{c}$ times $l_{r e f}$ are considered to belong to the free surface. That is, only interfaces that are smeared over up to 20 cells or less are considered.

2. Scale the CFL number in the VOF interface cells

The CFL number in the cells that were identified in the previous step is scaled to account only for the convective velocity normal to the VOF interface. If the velocity is tangential to the VOF interface (such as at the side of a liquid jet in an air domain) the correction for the liquid and gas volume fractions is zero (or very close to zero) although the CFL number may be large. These cells can be omitted. The scaled CFL number ${C F L}^{*}$ is:

Figure 12. EQUATION_DISPLAY

{C F L}^{*} = C F L 〈 \frac{v}{| v |}, \frac{\nabla α}{| \nabla α |} 〉

(2633)

Where $C F L$ is the "regular" Courant number and $v$ is the velocity vector.

3. Compute the sub-step size

Figure 13. EQUATION_DISPLAY

\begin{matrix} Δ t_{k} = \frac{Δ t}{k} & k = \frac{\max ({C F L}^{*})}{{C F L}_{\max}^{*}} \end{matrix}

(2634)

where ${C F L}_{\max}^{*}$ is the user-specified target Courant number. See Explicit Multi-Step (Deprecated) Solver Properties.