Cahn-Hilliard Navier-Stokes simulations for marine free-surface flows

The paper is devoted to the simulation of maritime two-phase flows of air and water. Emphasis is put on an extension of the classical Volume-of-Fluid (VoF) method by a diffusive contribution derived from a Cahn-Hilliard (CH) model and its benefits for simulating immiscible, incompressible two-phase flows. Such flows are predominantly simulated with implicit VoF schemes, which mostly employ heuristic downwind-biased approximations for the concentration transport to mimic a sharp interface. This strategy introduces a severe time step restriction and requires pseudo time-stepping of steady flows. Our overall goal is a sound description of the free-surface region that alleviates artificial time-step restrictions, supports an efficient and robust upwind-based approximation framework, and inherently includes surface tension effects when needed. The Cahn-Hilliard Navier-Stokes (CH-NS) system is verified for an analytical Couette-flow example and the bubble formation under the influence of surface tension forces. 2D validation examples are concerned with laminar standing waves reaching from gravity to capillary scale as well as a submerged hydrofoil flow. The final application refers to the 3D flow around an experimentally investigated container vessel at fixed floatation for Re = 1.4 × 107 and Fn = 0.26. Results are compared with data obtained from VoF approaches, supplemented by analytical solutions and measurements. The study indicates the superior efficiency, resharpening capability, and wider predictive realm of the CH-based extension for free-surface flows with a confined spatial range of interface Courant numbers.


Introduction
Many two-phase flows are characterized by immiscible fluids that feature negligible compressibility. A prominent example refers to maritime free-surface flows. Technical applications of such flows are often subjected to large interface deformations, e.g. breaking waves. The accurate simulation of these flows requires a computational model that conserves the mass of each phase whilst preserving a sharp interface. These requirements still pose a challenge in mesh-based computational fluid dynamics.
Engineering two-phase flow simulations mostly refer to either of two interface-capturing methods Ferziger and Peric (2012): namely the Level-Set Osher and Sethian (1988) and the Volume-of-Fluid (VoF) approach Hirt and Nichols (1981), which both reconstruct the free surface from an indicator function. The Level-Set method introduced by Osher et al. Osher and Sethian (1988) utilizes a signed distance function to characterize the interface by the zero-value iso surface. The continuous distribution of the signed distance simplifies a higher-order discretization of the related transport equation, and the geometry of the interface can be determined with improved accuracy. A drawback of the Level-Set method is that it does not guarantee mass conservation. Two-phase applications of the VoF method suggested by Noh et al. Noh and Woodward (1976) and later refined by Hirt and Nichols Hirt and Nichols (1981) usually employ a scalar volume concentration, of a foreground phase to identify the fluid state of each cell. The method is conservative and capable to predict merging and rupturing of free surfaces. For immiscible fluids, any mixing of both phases is undesired but numerically difficult to avoid. Different strategies are conceivable to improve interface compression: Geometric reconstruction schemes, e.g. SLIC Noh and Woodward (1976), PLIC Hirt and Nichols (1981) or LVIRAPilliod Jr and Puckett (2004), and dedicated downwind-biased advection schemes, e.g. CICSAM Ubbink and Issa (1999), HRIC Muzaferija and Peric (1999), IGDS Jasak et al. (1999) or BRICS Wackers et al. (2011). Geometric reconstruction schemes are afflicted with a considerable algorithmic complexity which reduces their popularity. Dedicated advection schemes are slightly heuristic but fairly simple to implement.
They maintain an approximately sharp interface subject to sufficiently small time steps. On the downside, they require transient implicit simulations even for steady state problems, e.g. the calm-water resistance of steady cruising ships.
To further improve the interface compression, some authors have proposed to add an artificial compression or antidiffusion term, e.g. So et al. (2011); Heyns et al. (2013). These methods rely on heuristic compression factors and improve the compressiveness at the expense of a reduced numerical stability.
If surface tension effects are negligible, VoF models using dedicated advection schemes are deemed a good compromise between efficiency, accuracy and conservation properties. An alternative, much less common approach refers to diffuse interface models, often labelled Cahn-Hilliard (CH) models Cahn and Hilliard (1958). Here, the (ideally sharp) interface is replaced by a (thin) layer where the fluids mix. The approach is able to mimic phase separation and thus promises re-sharpening features which are attractive for engineering simulations. Although, the neglect of surface tension is an acceptable assumption in many engineering problems, it appears that the CH approach incorporates surface tension in a natural way and no additional model, e.g. the Continuum Method Brackbill et al. (1992); Lafaurie et al. (1994), is required.
There exist a variety of different CH approaches for two fluids, e.g. models governed by fluids with matched densities (labelled as Model H Hohenberg and Halperin (1977)), identical viscosities (Boussinesq Fluid Jacqmin (1999)) or so-called thermodynamically consistent systems Lowengrub and Truskinovsky (1998); Abels et al. (2012), just to name a few. In addition to the different CH variants, different strategies for their coupling with the momentum and con-tinuity equations of have been suggested. Further distinctions refer to balancing either the mass or the volume fluxes between both phases Lowengrub and Truskinovsky (1998); Ding et al. (2007), the considered baseline conservation equations Abels et al. (2012); Ding et al. (2007) and the introduction of modifications to ensure thermodynamic consistency Lowengrub and Truskinovsky (1998); Abels et al. (2012). The VoF scheme offers a closed system of PDEs, but entails evolved parametrized approximations. On the contrary, three additional physical parameters occur in the CH method. The first and second correspond to the transition length as well as the surface tension coefficient, the third parameter refers to the mobility that governs the strength of the phase separation process. Combining the CH model with the Navier-Stokes equations, essentially results in an augmented VoF formulation. This inheres a non-linear, diffusive right-hand side of order four, which is zero outside the inter facial region. The non-linear character is beneficial. It supports an accurate computation of surface tension effects when the interface is adequately resolved and the use of stable, upwind-biased advective approximations in under-resolved flow simulations. To this end, a compressive numerical method is suggested for simulations, where the transition length is under-resolved by the numerical grid and surface tension influences cannot be displayed. The latter is based on an automatic adjustment of the mobility parameter. Possible, minimal blurs are bypassed with a non-linear state equation. The resulting system is virtually insensitive for spatial and temporal resolutions aspects.
The remainder of the paper is organized as follows: Sections 2 outlines the mathematical model including a brief introduction into the diffuse interface model. The subsequent third section describes the numerical procedure and outlines implementation aspects. Section 4 covers the verification. The determination of the mobility parameter in under-resolved flows is outlined in Section 5. Sections 6.1 and 6.2 validate the CH-approach and the mobility parameter estimation against frequently studied two-dimensional test cases, i.e. standing waves in the capillary and gravity scale as well as a submerged hydrofoil flow. The comparison of results and computational efforts obtained from a CH-and a VoF-approach for a widely used 3D container ship benchmark case is depicted in Section 6.3. Final conclusions are drawn in Section 7. Within the publication, Einstein's summation convention is used for lower-case Latin subscripts. Vectors and tensors are defined with reference to Cartesian coordinates and dimensionless variables are consistently marked with an asterisk.

Mathematical Model
Following the work of Ding et al. Ding et al. (2007), one can distinguish between mass conservative and volume conservative CH strategies. To illustrate this, we define the specie densities (ρ a , ρ b ) by a simple linear equation of state, which connects them to their constant bulk densities (ρ a , ρ b ), i.e. ρ a = cρ a and ρ b = (1 − c)ρ b . The expression c = V a /V represents the volume concentration of the foreground phase, the respective concentration of the background The mass conservation of the species a and b are governed by where the σ a , σ b denote the mass transfer rate into the species a and b, v i refers to the velocity and x i denotes the spatial Cartesian coordinates. Using ρ = ρ a + ρ b , an analogue continuity relation is obtained for the mixture Substituting ρ a and ρ b by cρ a and (1 − c)ρ b , one can reformulate (1) and obtain Summing up expressions (3) yields an alternative continuity equation, which describes the volume change where the volume diffusion fluxes refer to j a k etc.. Two options are conceivable: balanced mass fluxes (σ a = −σ b ) or balanced diffusion fluxes ( j a k = − j b k ). Since the interface is generally thin -particularly in the sharp interface limit, where both options are identical -the preference is related to the employed numerical method.
Most authors opt for a volume conservative approach, i.e. j a k = − j b k , and employ a volume averaged velocity field. In this case, mass is only globally conserved, provided that the inter facial regions don't intersect with the domain boundaries Ding et al. (2007), which might be difficult for travelling waves. The continuity expression (4) conveniently simplifies to a zero velocity divergence, and the conservative momentum equations are augmented by a total mass flux term v i (σ a + σ b ) on the r.h.s., cf. (6). On the contrary, assuming a mass conservative approach, i.e. σ a = −σ b , one has to account for divergence effects of the observed mass-averaged velocity.

Momentum and Continuity
The present research opts for a mass-averaged velocity field v i governed by the Navier-Stokes (NS) equations supplemented by a surface tension force f ST Assuming σ a = −σ b , the continuity equation follows from (2), and yields a conservative formulation for mass and where ρ, µ e and p refer to the density, dynamic viscosity and pressure of the mixture. The unit coordinates and the strain-rate tensor are marked by δ ij and S ij = 1 2 ( ∂v i ∂x j + ∂v j ∂x i ). The four right-hand side terms of the momentum balance (6) denote viscous, pressure, body and surface tension forces. The pressure is a numerical property which often inheres further trace terms related to the adopted phase field model and surface tension force expression. The framework supports laminar and Reynolds-averaged (modelled) turbulent flows (RANS). In the latter case, v i and p correspond to Reynolds-averaged properties and p is additionally augmented by a turbulent kinetic energy (k) term, i.e. 2ρk/3. Along with the Bousinesq hypothesis, the dynamic viscosity µ e = µ + µ t of turbulent flows consists of a molecular and a turbulent contribution (µ t ), and the system is closed using a two-equation turbulence model to determine µ t and k. Details of the turbulence modelling practice are omitted to safe space and can be found in textbooks, e.g. (Wilcox, 1998).

Equation of State
We expect both fluid phases to be incompressible and virtually immiscible. Moreover, we assume no-slip between the fluid phases along the interface and model the flow as a mixture between fluids which share the velocity field governed by equation (6). The continuity equation (7) serves to determine the pressure and the local fluid properties ρ and µ follow from an equation of state. A more general equation of state (EoS) refers to a weighted sum of the bulk properties of the participating phases, viz.
An alternative approach employed herein reads Here γ m is a non-dimensional model parameter to adjust the transition regime. For 0.1 ≤ γ m ≤ 0.5 the difference between (9) and (10) is limited to 0.009-3.597% at the transition points c = 0 and c = 1, resulting in a slight offset of fluid properties. The formulation (10) serves the regularization of unbounded concentration values and helps to sharpen partly blurred interfaces.

Diffusive Interface Model
The subsection briefly summarizes the diffusive interface model for isothermal two-phase flows as suggested in a landmark paper by Cahn and Hilliard Cahn and Hilliard (1958) and later elucidated by Jacqmin Jacqmin (1999). The present approach refers to a classical CH-model and is based upon the free energy E of the interface Γ between two isothermal fluid phases The coefficients C 1 [Pa] and C 2 [N] can be determined from the interface thickness γ[m] and surface tension σ a,b [N/m] between the two fluids, as outlined below. The foreground phase concentration c represents a measure of phase and ranges from a foreground state (c a ) to a background state (c b ), i.e. c a = 1 and c b = 0. Mind that alternative CHformulations exist, which employ the mass concentration or other energetic contributions.
The first term of e refers to the bulk energy density and aims to separate the fluids. The second term represents the gradient energy which widens the interface. To model separated (immiscible) fluids, a fourth-order polynomial, labelled "double well potential", is often used to describe the bulk energy density, In equilibrium conditions, E is minimized with respect to c. Using variational calculus, this relates to the root of a chemical potential ψ for the equilibrium state of plane interfaces (i.e. ψ = 0) where x n represents the interface normal direction. Substituting (12) into (13), one obtains a hyperbolic tangent concentration profile. This also renders a relation between the coefficients C 1 , C 2 and an interface thickness γ, viz.
Similarly, surface tension forces can be related to the concentration. Jacqmin Jacqmin (1999) outlined, that the convective rate of change of the free energy widens or compresses the interface. The term reads δE/δc This immediately reveals the surface tension force used herein Some authors rearrange this definition into an apparent pressure term and a term involving the chemical potential to express the surface tension force by the divergence of a surface tension stress −C 2 (∇ k c)(∇ i c) and the gradient of a related apparent pressure C 1 b + 0.5 C 2 (∇ k c) 2 Lowengrub and Truskinovsky (1998); Abels et al. (2012). Expression (16) associates vanishing surface tension effects with C 2 = 0. Jacqmin Jacqmin (1999) also deduced a link between the surface tension and C 1 , C 2 for a plane interface Substituting (14) into (17) yields Once γ and σ a,b have been chosen, both CH coefficients can be determined from the plane interface relations (14, 18) In the remainder of the paper, the neglect of the surface tension in under-resolved flows is modelled by C 2 = 0 N and, for the sake of simplicity, C 1 = 1 Pa.

Concentration Transport and Velocity Divergence
The mixture fraction is computed from an additional transport equation that models the mass transfer between the phases. Two options will be discussed, referring to the classical VoF and the CH approach outlined in Section 2.3.
Using the classical VoF approach, we assume that the material property of the fluid must not change, viz.
On the contrary, the CH-approach of Lowengrub and Truskinovsky Lowengrub and Truskinovsky (1998) refers to the mass concentration c m = cρ 1 /ρ of the foreground phase and involves a diffusive phase transfer term Using the continuity relation (7) leads to The desired conservative volume-flux based transport follows from (22), with ρc m = ρ 1 c and M m = M/ρ 1 Here the mobility parameter M of dimension m 3 s/kg [= ν/Pa] is a free parameter that controls the strength of the phase separation process and will be explored for under-resolved flows in Section 5.
VoF methods for immiscible incompressible fluids are usually based upon a pressure correction/projection scheme and prefer to observe volume instead of the mass fluxes to avoid the density jump. The latter can be extended for hydrodynamic flows featuring ∂v i /∂x i = 0 along the route outlined in (Yakubov et al., 2015). Since the bulk densities are deemed incompressible, the density solely depends on the transition function m, which in turn depends on the concentration (9), (10). An alternative continuity equation is derived from (7) Substituting (23) into (24), one finally arrives at In conjunction with VoF, a solenoidal velocity field is recovered due to the neglect of diffusive mass transfer. Note that f * is virtually zero in combination with the non-linear EoS (10) in the sharp interface limit, which returns a divergence-free velocity field.

Non-dimensional Governing Equations
The non-dimensional equations support the discussion of influences and assist the verification. Assuming a spatially constant mobility M, the non-conservative concentration, continuity and momentum equations read: The left-hand sides of the balance equations represent the classical, incompressible VoF-method. In the case of a nonzero right-hand side of (26), the velocity field is no longer divergence free in (27), which in turn introduces additional terms to the momentum balance (28). An exemplary relationship between a dimensional quantity, a reference value and a non-dimensional quantity marked with an asterisk reads v i = Vv * i . The dimensionless parameters are defined by The quantities utilized for the non-dimensionalisation are given in Table 1. Local discrete similarity parameters employ L = δx i , T = δt and V = v i etc.. It should be pointed out that the transition length γ can be small compared to a local grid spacing δx i , resulting in small (discrete) Cahn-numbers, which supports the neglect of the second term on the r.h.s. of (26) in an under-resolved sharp interface limit.
field quantity

Numerical Procedure
Results of the present study were obtained from the Navier-Stokes procedure FreSCo + (Rung et al., 2009). The implicit finite volume procedure uses a segregated algorithm based on the strong conservation form and employs a cell-centred, co-located storage arrangement for all transport properties. Unstructured grids, based on arbitrary polyhedral cells or hanging nodes, can be used. The solution is iterated to convergence using a modified pressurecorrection scheme (Yakubov et al., 2015). Various turbulence-closure models are available with respect to statistical (RANS) or scale-resolving (LES, DES) approaches. Time derivatives are approximated by an implicit Euler scheme.
The numerical integration employs a second-order mid-point rule, diffusive fluxes are determined from second-order accurate central differencing and convective fluxes use higher-order upwind biased interpolation formulae. Since the data structure is generally unstructured, suitable preconditioned iterative sparse-matrix solvers for symmetric and nonsymmetric systems, e.g. GMRES, BiCG, QMR, CGS or BiCGStab, are used. The procedure is parallelized for several thousand processes using a domain decomposition methods and the MPI communications protocol (Yakubov et al., 2013). It supports local mesh refinement, overset grids (Völkner et al., 2017), node-based adjoint shape-optimization Kröger and Rung (2015); Kröger et al. (2018) and fluid-structure interactions between mechanically coupled floating bodies (Luo-Theilen and . Cell-centred fluid properties are determined from (8). For face-based properties, a linear interpolation between adjacent cell-centre values is used. The discretization of the additional terms that originate from the CH-NS approach are discussed in the upcoming lines for equations (23 -concentration), (25 -continuity) and (6 -momentum). The discussion refers to the symbolic finite-volume approximation of a variable φ located in the center P of a control volume of size δΩ P , with neighbour control columns NB, i.e. A P φ P − NB A NB φ NB = S φ δΩ P . Here, the right-hand side source term S φ is treated explicitly during the (time-implicit) iteration of a segregated solution procedure Ferziger and Peric (2012). All presented own VoF studies employ the compressive HRIC scheme for convective concentration transport Muzaferija and Peric (1999).

Concentration Conservation
The numerical solution of (23) follows a deferred correction approach. Hereto the right-hand side is notionally split into a bulk density contribution and a gradient term, i.e. Using , an inherently positive contribution to the bulk density term is identified. This leads to an implicit contribution and an explicit source S c which uses values of the previous iteration The integration over a control volume yields a discretized surface integral over all sub surfaces, i.e. δΓ ( f ) i of δΩ P for the implicit part of (30) f (δΩ P ) (2M C 1 )∇ i c δΓ i , which is discretized using central differences and the mid-point integration rule. The explicit source term follows from a mid-point integration over the control volume. Upwind biased schemes are used to approximate the convective fluxes of (23) for the CH approach and a compressive downwind biased approach Muzaferija and Peric (1999) is used for the VoF approach.

Momentum and Mass Conservation
A simple mid point integration is employed to account for the additional explicit CH-NS sources. This involves f ∇ i (M∇ i ψ)δΩ P for the conservation of mass (25) in combination with a pressure correction scheme (Yakubov et al., 2015), and f ST i δΩ P (15) for the momentum equation (6).

Planar Couette Flow
The implementation is verified for a planar Couette flow under the influence of vertical gravity, wherefore a non-dimensional analytical solution is constructed and compared with the numerical results. Figure 1  The velocity is assumed to be unidirectional, i.e. v * 1 (x * 2 ), and in a fully developed, laminar, steady state. Moreover, the concentration field is also considered steady and homogeneous in the primary direction (x * 1 ), viz.
C : Using the linear EoS (9) allows an integration of (31-33) and results in the following analytical solution This solution follows from Dirichlet conditions for the concentration and velocities along the top as well as the bottom wall (v * 1 (1) = c(0) = 1, v * 1 (0) = c(1) = 0). Additionally, a prescribed top-wall pressure (p * (1) = 0) was employed. Interestingly, the solution is independent from the Peclet number und thus also from the mobility parameter M.
Computational results for the CH-NS system (6,23,25) were obtained on a 2D grid in conjunction with periodic boundary conditions in stream wise direction. Convective fluxes for momentum and concentration were discretized using first-order upwind differencing (UDS) and the mobility was assigned to a value that results in Pe = 1 × 10 5 .

Stationary Bubble
The influence of the surface tension model is verified by computing the transition from an initial non-equilibrium (rectangular) bubble into an equilibrium (circular) bubble. The example is restricted to advancing a 2D flow field without gravitational effects (Fr = ∞) in pseudo time. As outlined in Figure 3, a lighter foreground phase (ρ a /ρ b =1/100) rectangle with an edge length of L = 0.005 m is initially embedded into a heavier phase, such that the surface ten- pressure approximation, the same situation is simulated with an increased transition length for three different surface tension values Ca L = 0.1, Pe L = 200/σ a,b , We = 4/σ a,b . Fig. 3 (b) shows the resulting pressure distributions over a radial coordinate which reveal a fair agreement with theoretical results.

Mobility Parameter in Under-Resolved Flows
The section discusses means to model the mobility parameter in under-resolved flow simulations, where the surface tension influence is neglected due to the coarse resolution of the interface thickness γ. As outlined in Sect. 2.3, the neglect of surface tension yields C 2 = 0 and the concentration equation (25) simplifies towards Using ∇ i b = (∇ i c) ∂b/∂c together with the definition of the double well potential b (12), the r.h.s. of (37) reads Depending on the concentration value, (38) acts locally diffusive (ν c ≥ 0) or compressive (ν c < 0). The sign-change of the apparent viscosity ν c = 2C 1 M(6c 2 − 6c + 1) resembles compressive approximations of the convective term, which switch from upwind to downwind approximations along the interface, and thus from positive to negative (apparent) viscosities, to keep the interface sharp Ubbink and Issa (1999); Muzaferija and Peric (1999). The apparent viscosity in (38) zeroes 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, we separate the mobility into a physical and a modelled part, i.e. M = M phys + M mod . The physical part is usually assigned to fairly small values, e.g M phys << 1 × 10 −15 m 3 s/kg Jacqmin (2000); Magaletti et al. (2013). Jacqmin Jacqmin (1999) reports that the mobility typically scales with the transition length M ∝ γ n , where n varies between 1 ≤ n ≤ 2. Moreover, a recent publication of Magaletti et al. Magaletti et al. (2013) suggests n = 2 and thus Pe ∼ Ca −1 or γ ∼ σ a,b M/V. Most engineering settings are therefore unable to sufficiently resolve the transition length and we consider the numerical contribution M mod to be dominant.

Homogeneous Mobility Model
The formulation of M mod is based on the interface blurring introduced by upwind-biased schemes. An estimation of the tensorial numerical diffusion at a cell face returned by a first-order upwind scheme might read Here v ( f ) i denotes the velocity at the face center, δx ( f ) j refers to the connecting vector from the upstream to the downstream adjacent cell center and λδx ( f ) j approximates the distance between the face and the upstream cell. Depending on the time discretization scheme, the related error might be included into the estimate of the mobility parameter from a modified equation analysis. For the example of a first-order implicit time discretization, the modified equation analysis suggests a simple supplement of a Courant number term, i.e. ν UDS (39) is spatially and temporally variable. Spatially volatile mobility distributions are deemed to obstruct the robustness of the procedure. Hence, we confine our interest to homogenized approaches and estimate the mobility based on the maximum norm of the matrix valued numerical diffusion Note that the field is filtered to extract the interface region, i.e. only faces with a projected concentration gradient above δ c = 10 −3 are considered. A von-Neumann stability analysis of the first-order discretized 1D equation at the interface location (c = 0.5) yields the following estimates for stable solutions (cf. appendix 9.1) where ϕ represents the phase angle. Approaching the steady state limit, the analysis only excludesM = 1. In case of ϕ → 0, the branchM ≤ 1 might be a saver recommendation. ThereforeM = 0.1 was used for all applications displayed in section 6, for which no stability problems were observed. If alternative convection schemes are used, the 1st-order upwind based analysis still provides reliable estimates.
Further inspection reveals, that under no circumstances the predicted phase transition spans less than five control volumes. Improved sharpening is obtained from the non-linear equation of state (10). An illustrative 1D example is used to demonstrate this. In this example, a free surface is transported by a prescribed flow field on the grid depicted in Fig. 4 (a). The horizontal flow field is directed from left to right with a constant velocity of v 1 = 1 m/s, and only the concentration equation (23) is computed. The employed grid is homogeneous (λ = 0.5) and features δx 1 = 10 −3 m and δx 2 = 1m. The simulation is initialised with a sharp interface along a vertical line at the centre location x 1 = 1m. Fig. 4 (b) displays a partly blurred interface from one CH simulation withM = 0.1 and Co = 1 after t = 5 s. The corresponding density field obtained from the non-linear equation of state using γ m = 0.05 is displayed in Fig. 4 (c).
Although the concentration field is slightly blurred, the resulting density and viscosity fields are sharp.

Gravity and Capillary Wave
The first example deals with the decay of standing waves which are initialised according to Fig. 5(a). We aim to assess, if the shear driven energy exchange between the two fluids, labelled a & b, is correctly captured in both, the capillary and the gravity regime under the influence of vertical gravitational acceleration g 2 . Numerical results are compared with analytical solutions of Prosperetti (1981), which exist for identical kinematic viscosities (µ a /ρ a = µ b /ρ b ) in the linear (laminar) flow regime. The initial wave length complies with a unit wave number (k = 2π/λ = 1) and the initial wave amplitude corresponds to a = 0.01λ. The reference velocities refer to V = |g 2 | λ and V = σ a,b /ρ b for the gravity and the capillary case, respectively. The extent of the 2D computational domain depicted in Fig. 5 (b) reads λ × λ. A locally refined grid with approximately 250.000 isotropic control volumes was employed. The resolution of the free surface region refers to δx 1 = δx 2 = λ/4000. The time step was assigned to δt = λ/(4000V) which was sufficient to ensure Courant numbers below Co < 0.1. Symmetry (no-slip) conditions were used along constant x 1 (x 2 ) boundaries of the domain.
For the gravity wave, surface tension influences were neglected and the density ratio and Reynolds-number read For the capillary case, the density and viscosity ratios read ρ a /ρ b = 1/100 and µ a /µ b = 1/10 respectively, the interface thickness is resolved by 10 vertical cells and follows from a Cahn-number of Ca λ = γ/λ = 1/400. The is prescribed in accordance with a Peclet-number of Pe λ = 2 · 10 10 together with the linear material law. Conclusions drawn for the capillary case are similar to the gravity case, as indicated by Figure

Submerged Hydrofoil
The second example refers to the wave pattern downstream of a submerged NACA0012 hydrofoil at 5 • incidence in accord with experimental data of Duncan Duncan (1981, 1983, cf. Fig. 7 The utilized unstructured numerical grid is displayed in Fig. 7

Flow around a KCS container vessel
The final application refers to the fully turbulent flow around an unappended Kriso container ship hull (KCS).
Experimental resistance data and wave fields were published by Kim et al. (2001)  pattern. Based on the current Froude-number a dimensionless wavelength of λ/L pp = 2 π Fn 2 = 0.4247 is expected, which is approximated with roughly 100 cells. Fig. 9 indicates the different refinement levels for the near and the far field.
At the inlet, outlet, outer and lower boundaries, Dirichlet values for velocity and concentration are specified, while the pressure is extrapolated. A reverse situation is given at the top face which corresponds to a pressure boundary.
Symmetry and wall boundary conditions are declared along the midship plane as well a the hull. Turbulence is modelled by a high-Re k − model (Wilcox, 1998). Convective momentum transport is realised by a monotonicitypreserving QUICK scheme. Similar to the hydrofoil case, data obtained from CH-NS simulation is compared with   Fig. 12. The predictive discrepancy is generally small and the non-linear CH-NS tends to provide slightly larger amplitudes. Mind that the non-linear EoS 10 leads to a significant sharpening of the density field, as illustrated by Fig. 13.

Discussion and Outlook
The paper presents an alternative approach for the simulation of marine free surface flows at engineering scale. The method (labelled CH-NS) can be displayed as a generalized Volume-of-Fluid (VoF) approach which is extended by a diffusive right hand side of order four obtained from a Cahn-Hilliard (CH) system. While the CH-framework is often used to describe phase separation processes along diffuse interfaces, the phase transition region usually falls below the typical grid resolution in engineering settings. Therefore these settings closely resemble a numerical sharp interface limit. The phase separating characteristics of the CH-NS is supported by negative diffusivity of the concentration equation in the central transition regime. This involves ≈ 60% negative diffusion for the double-well potential used in this paper, and scales with a mobility. A homogeneous but temporally/iteratively variable mobility, adapted to a spatially averaged error expression, is suggested which leads to robust results with a fair predictive accuracy. The model involves a free parameter which is assigned to a value well below the stability limit. An alternative non-linear material law is supplementary used to improve minimal blurring.
The implementation of the CH-NS system was verified for an analytical laminar Couette flow. Subsequently, different laminar and turbulent two-phase flows were used to validate the results against experimental data, reaching from the capillary to the gravity scale. The CH-NS approach naturally includes surface tension effects but the more relevant advantage refers to the efficiency of the approach. Unlike a VoF approach, accurate steady simulations can be performed much faster without any CFL-constraint.
Efficiency benefits over the traditional VoF method were demonstrated for a fully turbulent two phase flow around a container ship hull at realistic Froude-and large Reynolds-numbers in steady state. For all cases presented in this paper, the computational time is reduced by at least one order of magnitude and minor predictive differences were Fig. 11: Comparison of predicted wave field around the Kriso container vessel at Re=1.4×10 7 and Fn=0.26 obtained by the VoF (bottom) and CH-NS (top) approach. White horizontal lines indicate evaluation planes used for the wave cuts displayed in Fig. 12 . observed in comparison with the VoF method. x 1 /L pp [-] x 3,FS /L pp [-] ξ ≤ 1. With attention restricted to the under-resolved situation (C 1 = 1 Pa, C 2 = 0 N), we obtain ξ = ξ m+1 ξ m = 1 1 + Co + 2Co diff b 1 − cos(ϕ) + i Co sin(ϕ) = 1 + Co + 2Co diff b 1 − cos(ϕ) 1 + Co + 2Co diff b 1 − cos(ϕ) 2 + (Co sin(ϕ)) 2 − i Co sin(ϕ) 1 + Co + 2Co diff b 1 − cos(ϕ) 2 + (Co sin(ϕ)) 2