Under-Resolved Interfacial Flow

This section briefly discusses means to model the mobility parameter M^∗ in flows where the interface thicknessγ^∗ is hardly resolved by the numerical grid and the surface tension influence is neglected due to the coarse resolution of the interface. As outlined in Sec.

2.1.4, C^∗b = 0 in the absence of surface tension and the concentration equation (2.27) simplifies to b from Eqn. (2.17) yields

∂

Depending on the concentration value, (2.54) acts locally diffusive (ν^∗c≥0) or compressive (ν^∗c < 0). As outlined in Fig. 2.2, the apparent viscosity in (2.54) vanishes at c = 0.5(1±1/√

3) and is negative over approximately 58% of the transition region. Aiming at a closure for the mobility parameter in under-resolved simulations, the mobility is separated into a physical and a modeled part, i.e.

M^∗ =M^∗phys+M^∗mod. (2.55)

0.0 0.5 1.0

−1

−0.5 0 0.5 1

c[-]

(∂b/∂n)/max(|(∂b/∂n|) [-]

n= 0 n= 1 n= 2 n= 3

(a)

0.0 3−√ 3

6 0.5 3+√

3 6 1.0 0

2 4

c [-]

νc /(M∗ C∗a [-]

(b)

Figure 2.2: Double well potential: (a) Normalized fourth-order polynomialb as well as its first three (normalized) derivatives and (b) the apparent viscosity basedb^′′with indicated roots.

The physical part is usually assigned to fairly small values, e.g. M^∗phys << 1×10⁻¹⁵m³s/kg, cf. Jacqmin [2000] or Magaletti et al. [2013]. Jacqmin [1999] reports that the mobility typi-cally scales with the transition lengthM^∗ ∝γ^∗n, wherenvaries between 1≤n ≤2. More-over, a recent publication of Magaletti et al. [2013] suggestsn= 2 and thus Pe∼Ca⁻¹ or γ^∗ ∼√

σ^∗a,bM^∗/V^∗. Most engineering settings are therefore unable to sufficiently resolve the transition length and the numerical contributionM^∗mod is considered to be dominant.

Sec. 2.2.9 proposes an approach to estimate the homogeneous mobility parameter based on leading errors of the approximation of the convective concentration transport.

Under-Resolved Equation System

A divergence-free primal velocity field is achieved if a) the VoF approach is used, b) both participating fluids feature the same bulk density or c) a nonlinear EoS is employed, cf.

Eqn. (2.12). Hence, under-resolved CH-VoF studies consistently employ f^ρ or W e → 0 that further simplifies the primal conservation of mass, momentum and concentration. The resulting subsystem of equations forms the relevant set of PDEs which serves as the basis for all adjoint studies, viz.

R^∗p = ∂v_k^∗

∂x^∗_k = 0 (2.56)

R^∗c = ∂ c

∂ t^∗ + ∂ v_k^∗c

∂x^∗_k − ∂

∂x^∗_k [

ν^∗c ∂ c

∂ x^∗_k ]

= 0 (2.57)

R^∗_i^vi =ρ^∗ [∂ v_i^∗

∂ t^∗ +v^∗_k∂ v^∗_i

∂ x^∗_k ]

+ ∂

∂ x^∗_k [

p^∗effδij−2µ^∗effS_ij^∗ ]

−ρ^∗g^∗_i = 0. (2.58)

Their boundary conditions in Tab. 2.2 close (2.56)-(2.58). Constant values of the pres-sure gradient typically follow from hydrostatic contributions of an extended prespres-sure, viz.

p^∗eff → p^∗ +ρ^∗x^∗_kg^∗_k and thus ∂p^∗eff/∂x^∗_i → ∂p^∗/∂x^∗_i +ρ^∗g_i^∗. In all cases, inlet, symme-try and slip-wall boundaries are perpendicular to the gravity vector and no hydrostatic boundary contributions occur.

boundary type v^∗_i p^∗ c

inlet v_i^∗ =v_i^{∗in ∂p}_∂n^∗_∗ = 0 c=cⁱⁿ outlet _∂n^∂vi∗_∗ = 0 p^∗ =ρ^∗g_k^∗x^∗_{k ∂n}^∂c_∗ = 0 symmetry v_i^∗n_i = 0, ^{∂ v}_{∂ n}^∗i_∗t_i ^∂p_∂n^∗_∗ = 0 _∂n^∂c_∗ = 0 wall (slip) v_i^∗n_i = 0, ^{∂ v}_{∂ n}^∗i_∗t_i ^∂p_∂n^∗_∗ = 0 _∂n^∂c_∗ = 0 wall (no-slip) v_i^∗ =v_i^∗w _∂n^∂p∗_∗ =ρ^∗g^∗_kn^∗_{k ∂n}^∂c_∗ = 0

Table 2.2: Boundary conditions for the primal equations, where ti [ni] refer to the local boundary tangential [normal] vector.

The closely coupled Eqns. (2.43)-(2.48) are discretized and approximated in space and time to obtain a numerical solution. This section describes the baseline procedure as introduced in Rung et al. [2009], briefly summarizes the discretization and outlines the numerical approximation based upon a second-order accurate Finite-Volume method (FVM).

All primal flow discretization concepts will be reused in Sec. 4.2 for the discretization of the adjoint flow. Mind that the primal/dual system shares the same Euclidean space and time horizon in possibly reversed mode. This motivates a more detailed insight to the discretization of individual expressions in Sec. 2.2 and 4.2 based on a generic transport equation.

2.2.1 Baseline Procedure

After spatial discretization of (2.43)-(2.48), a coupled system of NP equations is obtained after a suitable finite approximation for each field quantity n^φ = [1, ...,N^φ] with N^φ = 5 [N^φ = 7] if a laminar [turbulent] flow is considered. The discrete field variables are algorithmically linked via a sequential procedure, in which the individual equation systems are solved in a fixed sequence, cf. Alg. 1. The latter is implemented by means of outer-iterationsthat are performed until all –or only a selected subset of equations– are converged below a prescribed tolerance R^∗φ,min or a maximum number of outer iterations N^out,max is reached.

Potential cross-coupling lags and employs values of the last outer iteration for the respec-tive system variable(s), which also applies to the application of boundary conditions, cf.

Sec. 2.2.8. The iterative solution procedure of the algebraic system of each outer iteration employs the Portable, Extensible Toolkit for Scientific Computation (PETSc, cf. Balay et al. [2019]) library that distinguishes between solvers for symmetric (e.g. Conjugate Gradient, CG) and asymmetric (e.g. Biconjugate Gradient, BiCG) matrices which in turn follow from the basic physical characteristics of the underlying equation. Their iterative counter refers to so-called inner-iterations and a termination criterion employs either a prescribed convergence tolerance or a maximum number of inner-iterations. The described sequential iteration technique is modularized and the algorithm changes only slightly if additional field equations are added. Several modules are implemented generically and invoked in a modular fashion in line with a generic transport equation, cf. Sec. 2.2.6.

The outer iterations are embedded in a time loop n^T = [1, ...,N^T] if the system is advanced in time, cf. Alg. 1. Hence, the simulated time horizon t^∗ = [t^start∗, t^end∗] is discretized into NT equally distributed time slices ∆t^∗ = t^∗(n^T) where n^T = [1, ...,N^T] represents the actual time step. The discrete time step size ∆t^∗ usually follows from an a priori known, user defined value that fixes the total number of time steps NT. For steady applications, the time progresses in a pseudo-time and resembles an additional relaxation.

Despite the steady nature of all applications in this thesis, the integration in pseudo-time is unavoidable for VoF procedures. The adjoint-based studied cost functionals are free of temporal expressions and their adjoint systems are linearized around the temporally converged primal state.

Algorithm 1:Schematic overview of the sequential primal solution procedure based on outer-iterations embedded in an optional temporal loop.

while n^T ≤N^T do

while (n^out ≤N^out,max) or (R^∗φ ≤R^∗φ,min) do approximate linearized momentum equations approximate pressure correction equation approximate linearized turbulence equations approximate linearized concentration equation

2.2.2 Spatial Discretization

The spatial domain is discretized by NP contiguous control volumes (CV) or cells of in-dividual size ∆Ω^∗P, where P ∈ [1, N P] denotes to the cell center. The control volumes are of arbitrary polyhedral shape and form an unstructured numerical grid. In principle, the control volumes could move and deform over time along the route of an Arbitrary Langrangian-Eulerian approach (V¨olkner et al. [2017], Schubert and Rung [2019], Schu-bert [2019], Luo-Theilen and Rung [2017, 2019], Theilen [2020]), which is however not meaningful for steady or quasi steady problems in the focus of this study. A collocated variable arrangement is employed, hence, all transported quantities are stored in the cell center. Variable values at other locations have to be reconstructed using the geometry and topology information outlined below.

Geometrical and Topological Relations

Next to the cell centered values of the dependent variables, the employed FV approximation requires the values at the center F of the faces between two control volumes. The numerical grid and the face-based data structure involves NF faces separating the cells from each other and the exterior domain. The definition of the face oriented vectors connecting the relevant positions reads

d^∗_i^F = ˜d^∗_i^F+ ¯d^∗_i^F with d˜^∗_i^F =x^∗_i^Fx^∗_i^P and d¯^∗_i^F =x^∗_i^NBx^∗_i^F. (2.59) The considered domain is enclosed by NΓ degenerated boundary cells as depicted Fig. 2.3 (b), which have no volume and a midpoint coordinate that coincides with the center of the adjacent face, viz. ¯d^∗_i^F = 0 → d^∗_i^B = ˜d^∗_i^F. Figure 2.3 (a) displays an arrangement of two adjacent interior control volumes and (b) the analogue boundary situation. The local approximations are constructed from values extracted at the cell center P, the face centers F separating the control volume around P from the neighboring control volumes and the centers of the neighboring control volume denoted by NB (or B). They employ the cell size ∆Ω^∗P as well as the area ∆Γ^∗F and the respective area vector ∆Γ^∗_i^F = ∆Γ^∗Fni of all cell adjacent faces where ni represents the local face normal. Some crucial aspects for the approximation are: a) The connecting vector d^∗F_i between two face adjacent cell centers does not necessarily pierce through the face center F and sometimes does not even pierce through the face at all, b) the connecting vector d^∗_i^F is not necessarily co-linear to the face normal ni and c) faces are not necessarily centered between the center of the neighboring cell.

P

F

∆Γ^∗_i^F d^∗F_i

d˜^∗_i^F d¯^∗F_i NB F˜

·

(a)

P B

d^∗_i^B

∆Γ^∗_k^B t^B_k

(b)

Figure 2.3: Schematic representation of a Finite-Volume arrangement (a) in the field and (b) along the boundary.

Reconstruction of Face Values

Face values are obtained by interpolating adjacent cell-centered values linearly to an inter-mediate location ˜F along the connecting line. As displayed in Fig. 2.3 (a), the intermediate point is located perpendicular above the location of the face center F and the interpolation is second-order accurate, viz.

φ^∗F = (

1−λ^F)

φ^∗P+λ^Fφ^∗NP

  

interp.

+(

x^∗_i^F−x^∗_i^F˜)∂ φ^∗

∂ x^∗_i

⏐⏐

⏐⏐

F˜

  

extrap.

+O(h^∗2)

with λ^F = ∆Γ^∗_i^F(

x^∗_i^F−x^∗_i^P)

[∆Γ^∗_kd^∗_k]^F . (2.60) Since, this intermediate point ˜F does not necessarily coincide with the face center F, a sub-sequent explicit non-orthogonality correction is optionally added. The explicit correction can be regarded as an extrapolation that ensures second-order accuracy even on unstruc-tured grids. The required gradient in the intermediate point ˜F follows from a second-order accurate linear interpolation

∂ φ^∗

∂ x^∗_i

⏐⏐

⏐⏐

F˜

=(

1−λ^F)∂ φ^∗

∂ x^∗_i

⏐⏐

⏐⏐

P

+λ^F∂ φ^∗

∂ x^∗_i

⏐⏐

⏐⏐

NB

+O(h^∗2). (2.61)

The majority of field variables and gradients are interpolated based on Eqn. (2.60).

Many flows considered in this work –and most engineering flows– are characterized by large Reynolds numbers and thus convection dominates. The interpolation of face values to approximate the convective fluxes is therefore crucial and will be discussed in greater detail below. Moreover, the interpolation of face values is especially relevant from a hybrid continuous/discrete adjoint perspective.

2.2.3 Spatial Approximation

The governing Eqns. (2.43)-(2.48) can be expressed in a generic residual expression for a conserved quantity φ^∗, viz.

R^∗φ = ∂ φ^∗ The individual terms of (2.62) refer to temporal,convective anddiffusive transport as well as external sources from left to right, respectively. Mind that the last two contributions refer to s^∗φ in Eqn. (2.1). Using the FV method, the equation is treated in a weak form for each CV, e.g. where the Gauss-theorem has been applied to the integral over the convective and diffusive terms. Spatial approximations are required to determine volume and surface integrals, as well as the spatial gradients at the cell faces. All spatial integrals are approximated based on the second-order accurate mid-point rule as described below. Spatial gradients refer to second-order central differencing schemes (CDS).

Volume Integrals

The approximation of volume integrals by the midpoint rule yields:

∫

∆Ω^∗(P)

ψ^∗dΩ^∗ = [ψ^∗∆Ω^∗]^P+O(h^∗2). (2.64) For the example of the generic transport equation (2.63), one obtains

∫

The closed surface integral is fractioned into integrals over the individual faces F between P and its neighbors NB. The face integrals are again approximated using the midpoint rule, viz. For the example of the generic transport equation (2.63), one obtains

∮

The convective and diffusive fluxes have to be evaluated at the face center. The evaluation of the face centered viscosities (diffusivities)µ^∗φ follows the procedure outlined in (2.60).

The same procedure is employed to predict the face velocity v^∗_k of the convective term.

The face value of the transported variableφ^∗ is usually reconstructed with a specific rule on a case by case basis, cf. Sec. 2.2.6. The approximation of spatial gradients at the cell centers is outlined below and subsequently interpolated to the face center also using (2.60).

Spatial Gradients

The numerical procedure requires approximations of spatial gradients at various positions of the code. The discrete approximation of a gradient at the cell center is either computed from a fully conservative approach, using integration by parts –a.k.a. Gauss’s Theorem–, or from a non-conservative –but more robust– weighted Least-Square method.

Gaussian Approaches follow from an integral formulation of spatial gradients using the second-order formulae (2.64) and (2.66), viz.

(1) ∫

Face values are linearly interpolated from the adjacent cell-centered values, cf. Sec. 2.2.2.

Least-Squares Approaches follow from a first-order Taylor series expansion of the face-centered value for each face per cell, viz. φ^∗NB ≈ φ^∗P+d^∗F_i [∂ φ^∗/∂ x^∗_i]^P+O(h^∗) that can be rearranged to compute for the desired gradient. The relation is over-determined since the number of adjacent cells always exceeds the dimension. Hence, a regression is applied based on the symmetric matrix D^∗P_ik that regularizes the expression in the Least-Squares sense and agrees with second-order central differences on regular grids

φ^∗F ≈φ^∗P+d^∗_i^F∂ φ^∗ Equation (2.69) is preconditioned by the reciprocal cell center distance ||d^∗_m|| to counter-act issues caused by distorted grids. The approach is not conservative but improves the accuracy and stability for unfavorable –e.g. highly skewed– cells.

2.2.4 Temporal Discretization

In general the temporal discretization involves variable time step sizes which are fed into implicit two and three time layer approximations. The applications of this thesis refer to steady problems. Time stepping is thus restricted to converge the primal (and adjoint)

flow field in pseudo-time to steady state using constant time step sizes. Necessary for the study of the discrete adjoint system, a continuous time integral is discretized as follows

∫

φ^∗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 =v_k^∗F∆Γ^∗_k^F represents the volume flux across the face F, and the velocity v_k^∗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)

[v_k^∗φ^∗∆Γ^∗_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^∗F_i (∂ φ^∗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 = ρ^∗v_i^∗, the volume flux is replaced by a mass flux, viz. ˙V^∗F →m˙^∗F = [ρ^∗v_k^∗∆Γ^∗_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 α = [d_in_i/||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)/||d^F_k^∗|| , 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 = f_i^∗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 ˜v_i^∗ 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. ∑

by Patankar and Spalding [1972]. At first, a velocity distribution ˜v_i^∗ 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. ∑

Im Dokument Adjoint-Based Shape Optimization Constraint by Turbulent Two-Phase Navier-Stokes Systems (Seite 43-0)