2.2 Discrete Primal Two-Phase Flow
2.2.6 Algebraic Equation System
φ∗dt∗ discretize→ ∑
T
∫
∆t∗(T)
φ∗dt∗. (2.70)
2.2.5 Temporal Approximation
Temporal Integrals
Analogous to the spatial integration, time integrals employ the midpoint rule, viz.
∫
∆t∗(T)
φ∗dt∗ = [φ∗∆t∗]T+O(∆t∗2). (2.71) The evaluation of time integrals solely occurs during the theoretical derivation of the adjoint system.
Temporal Derivatives
Temporal changes are approximated either by a first-order implicit backward Euler (BE) method or a second-order accurate implicit three time level approach (ITTL). The former discretizes a continuous quantity φ∗ based on temporal information from the actual time instance (T) as well as the previous (T-1), viz.
∂ φ∗
∂ t∗ = φ∗T−φ∗T−1
∆t∗ +O(∆t∗) (2.72)
and refers to the method of choice when the intermediate solution of a pseudo-transient sim-ulation is not of interest, e.g. the total ship resistance in calm water where∂ φ∗/∂ t∗|t∗end → 0. Such quasi-steady simulations frequently employ a non-constant, adaptive time step pro-cedure that hooks up on the local Courant number distribution where typical values are between 0.2 ≤ Co ≤ 0.4. If the temporal behavior of the simulated system is relevant, an additional previous time step (T-2) is considered to the approximation. The resulting second-order accurate ITTL stencil reads
∂ φ∗
∂ t∗ = 3φ∗T−4φ∗T−1+φ∗T−2
2 ∆t∗ +O(∆t∗2) (2.73)
based on a constant time step size.
2.2.6 Algebraic Equation System
The employed FV approximation of (2.63) yields a discrete system of size NP×NP, where each line corresponds to a particular CV-balance
A∗φ,Pφ∗P− ∑
NB(P)
A∗φ,NBφ∗NB=S∗φ,P. (2.74)
Each discrete quantity also varies in time as well as per outer-iteration, e.g. strictly speaking the general quantity and source term reads φ∗P,M,T and S∗φ,P,M,T at the discrete
cell center P, the outer iteration M as well as the discrete time step T. For reasons of readability, the last two incremented indices (M,T) are suppressed if not needed, viz.
φ∗P ≡ φ∗P,M,T and S∗φ,P ≡ S∗φ,P,M,T. The preceding part of this section presents the discrete translation of the individual terms in (2.63) which contribute to (2.74).
The resulting system matrix should render two properties to support a stable solution process and physically realistic solutions. To maintain bounded solutions, one requires positive coefficients A∗P, A∗N B > 0. At the same time, the diagonal dominance of the matrix is required to support the iterative solution procedure, e.g. A∗P ≥∑
A∗N B. Once the discrete relation is assembled, an under-relaxation is employed
A∗φ,P→ A∗φ,P
ωφ and S∗φ,P← (1−ωφ)A∗φ,P
ωφ φ˜∗P, (2.75) where ˜φ∗P follows from the previous outer iteration. The relaxation supports the diagonal dominance. Related parameters are case specific, but typical values refer to, e.g.,ωvi = 0.6, ωk = 0.4 and ωc = 0.4. The approximation of the fluid dynamic pressure is tailored to employωp= 1.0. Relaxation of the effective pressure is realized within the actual velocity-pressure coupling routine and usually employs ˜ωp = 0.2.
Transient Contributions using the implicit BE approximation (2.72) together with the integration rule (2.64) result in the following contributions to the main diagonal and the r.h.s.
A∗P,T ← ∆ Ω∗P
∆t∗ and S∗φ,P,T ← ∆ Ω∗P
∆t∗ φ∗P,T−1. (2.76) Similarly, the second-order accurate ITTL (2.73) and the integration rule (2.64) return the following contributions to (2.74)
A∗φ,P,T ← 3 2
∆ Ω∗P
∆t∗ and S∗φ,P,T ← ∆Ω∗P
∆t∗ [
2φ∗P,T−1− 1
2φ∗P,T−2 ]
. (2.77) Convective Contributions are characterized by their directional information transport.
Following the approximation of surface integrals (2.67), the convective flux is fractioned
into ∑
F(P)
[v∗kφ∗∆Γ∗k]F =∑
F(P)
V˙∗Fφ∗F (2.78)
where ˙V∗F =vk∗F∆Γ∗kF represents the volume flux across the face F, and the velocity vk∗F is interpolated linearly in line with Eqn. (2.60). To account for the directional information transport, upwind biased formulae are employed to approximateφ∗F. To this end, φ∗F is split into a sum of an implicit and an explicit term, e.g.
∑
F(P)
[vk∗φ∗∆Γ∗k]F =∑
F(P)
V˙∗Fφ∗F →∑
F(P)
V˙∗F[
φ∗U+ DC∗C,φ,F]
. (2.79)
A simple, first-order accurate upwind differencing scheme (UDS) is implicitly considered to retain the desired properties of the matrix coefficients, viz.
V˙∗Fφ∗U:= max(
V˙∗F,0)
φ∗P+ min(
V˙∗F,0)
φ∗NB. (2.80)
The second term of (2.79) refers to a convective (upper index C) Deferred Correction (DC) contribution, which is explicitly considered on the r.h.s and computed from values of the previous iteration. The latter optionally accounts for high-order contributions as a supplement to the UDS part, viz.
V˙∗FDC∗C,φ,F = V˙∗F 4
[(1 +κ)(
φ∗D−φ∗U)
+ (1−κ)(
φ∗U−φ∗UU)]
. (2.81)
Equation (2.81) offers a generic way to account for various baseline convection schemes based on the constant κ ∈ [−1,1]. The preferred choice in this thesis is the third-order Quadratic Upwind Interpolation of Convective Kinematics (QUICK) scheme of Leonard [1979] which refers to κ= 0.5. Other popular schemes supported by (2.81) are the central differencing (CDS, κ = 1) and linear upwind differencing (LUDS, κ = −1) schemes.
On unstructured grids, the φ∗UU value is not explicitly available and reconstructed with second-order accurate central differences, viz. φ∗UU ≈ φ∗D −2d∗Fi (∂ φ∗U/∂ x∗i), cf. Fig.
2.4. Though this restricts the attainable reconstruction accuracy, it does not additionally harm the spatial order of accuracy for this term due to the upstream second-order accurate numerical integration.
In line with Godunov’s Theorem, only first-order schemes can be monotone (bounded) and avoid the generation of new extreme values, cf. Godunov [1954, 1959]. Therefore nonlinear elements are introduced to (2.81) to avoid the generation and amplification of extremes. Van Leer [1979] presented a popular Monotonic Upstream Scheme for Conser-vation Laws (MUSCL) that scrambles the nonlinearities into a damping function ψ based on a local sensor r, e.g.
V˙∗FDC∗C,φ,F = V˙∗F 4
[
(1 +κ)(
φ∗D−φ∗U)
ψ(r) + (1−κ)(
φ∗U−φ∗UU) ψ
(1 r
)]
and r = φ∗U−φ∗UU
φ∗D−φ∗U . (2.82)
Harten [1997] motivated the so-called total variation diminishing (TVD) scheme which ac-counts for monotonicity preservation by locally decreasing the approximation order. Many other approaches exist that, e.g., restrict ψ, c.f. Lien and Leschziner [1994a,b]. Related discussions frequently employ a Normalized Variable Diagramm (NVD), and the imple-mentation of dedicated NVD schemes follows the same strategy compared to the TVD methods. A major contribution of this work is concerned with adjoint NVD schemes dedicated to the convective transport of the concentration c in conjunction with VoF ap-proaches and are considered separately in Sec. 2.2.7. A detailed overview of primal / adjoint TVD schemes can be found in St¨uck [2012]. The nonlinear TVD/NVD functions are particularly unpleasant w.r.t. their corresponding adjoint operators as they produce unfavorable derivatives. Hence, the function ψ is usually kept frozen and thus not con-sidered in the adjoint calculus, cf. St¨uck [2012], St¨uck and Rung [2013], Kr¨oger [2016], Kr¨oger et al. [2018], Manzke [2018].
The convective contributions to (2.74) read:
A∗φ,P←∑
F(P)
−min(
V˙∗F,0)
, A∗φ,NB← −min(
V˙∗F,0)
, S∗φ,P← −∑
F(P)
V˙∗FDC∗C,φ,F. (2.83)
A word of caution should be noted regarding the main diagonal contribution. The latter employs∑
F(P)−min( ˙V∗F,0) =∑
F(P)max( ˙V∗F,0) which is strictly speaking only true for a converged (mass conservative) state, but supports the dominance of the matrix diagonal and the iterative procedure. Mind that the volume flux is a function of the velocity, and a Picard linearization is employed within the iterative solution procedure, cf. Alg. 1. For the transport of density-weighted field quantities, e.g. the linear momentum φ∗i = ρ∗vi∗, the volume flux is replaced by a mass flux, viz. ˙V∗F →m˙∗F = [ρ∗vk∗∆Γ∗k]F.
Diffusive Contributions follow from an unbiased information transport. The discrete representation of the diffusive flux in (2.67) is also fractioned into an implicit part and an explicit deferred correction
where DC∗D,φ,F represents a diffusive (superscript D) deferred correction contribution. The coefficient α = [dini/||d∗k||]F optionally accounts for unfavorable face-to-normal arrange-ments and employed exponents refer to n = 1 for the assembly of the pressure correction equation and n = 0 in all other cases. A CDS resembling, implicit approximation accounts for the unbiased information transport, viz. ∂ φ∗F/∂ x∗i ≈(φ∗NP−φ∗P)/||dFk∗|| , supports the diagonal dominance and maintains positive coefficients. The explicit DC contribution accounts for possible non-orthogonality
where the face gradient is determined by (2.60). The explicit DC∗D,φ,F contribution tends to zero if the computational grid is perfectly orthogonal ([d∗k/||d∗i|| − nk]F → 0). The diffusive contribution to (2.74) reads
A∗φ,P← ∑ Sources and Sinks are considered as explicit contributions to the r.h.s. using the mid-point rule (2.64)
S∗φ,P←[q∗φ∆Ω∗]P . (2.87)
The explicit contributions to the r.h.s. differ strongly from the respective conserved quan-tities, e.g. q∗vi,P = fi∗ST,P if surface tension contributions are added to the momentum balance.
Pressure Determination
The pressure-velocity coupling follows a modified version of the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) pressure correction algorithm, originally proposed
by Patankar and Spalding [1972]. At first, a velocity distribution ˜vi∗ is determined from the discretized momentum equation using an estimated pressure ˜p∗. Usually, these velocities do not satisfy the discrete continuity equation and the resulting discrete continuity defect of each CV, e.g. ∑
F(P)[˜vi∗Γ∗i]F−[s∗p∆Ω∗]P, is used to compute a correction of the estimated pressure in an iterative procedure, e.g. ˜p∗P ← p′∗P. A combination of the discretized mo-mentum (in reduced form) and continuity yields the inherently discrete pressure correction equation, cf. Ferziger and Peric [2012], which shares features with the pressure Poisson equation, cf. Chorin [1968]. Analogous to Yakubov et al. [2015], the pressure correction scheme employed herein accounts for velocity divergence sources that are of relevance for the CH-VoF model, cf. (2.43). The algebraic pressure correction scheme features only discretized diffusive fluxes and r.h.s contributions. The former are assembled in line with Eqn. (2.84) using µ∗p,P = [∆Ω∗/A∗vi]P and respective face values follow from a linear in-terpolation. The r.h.s. inheres contributions from the non-solenoidal velocity field as well as further CH-related sources, viz.
S∗p′,P← [dc
dt∗∆Ω∗ ]P
+∑
F(P)
[
˜ v∗i∆Γ∗i
]F
. (2.88)
Mind that the first r.h.s. contribution a priori vanishes for the VoF scheme and the converged pressure field should globally support zero velocity divergence. The pressure update involves an under-relaxation ˜p∗P ← ω˜p′p′∗P that usually employs ˜ωp′ = 0.2. The odd-even decoupling problem of the co-located variable arrangement is suppressed with a fourth-order artificial dissipation term on the r.h.s of the continuity equation along a route outlined by Rhie and Chow [1983]. The latter also considers contributions from gravity and surface tension in order to strengthen the iterative procedure.