Coupling of Viscous and Potential Flow Models with Free Surface: Implementation and Application to Offshore Engineering

In the offshore industry, it is not uncommon that many structures for various purposes were installed at one location for decades. These structures must survive all weather types, including heavy storms or even Tsunami. With the improvements in oil and gas production technology and continuously raising industrial need, the economies of scale are driving the need for ever larger semi-submersibles [1, 2].

Obviously, the design of the offshore structures (e.g. ships, offshore platforms) in coastal and open ocean regions requires very specific tools to evaluate the expected hydro-dynamic forces. Experimental tests on model scale are considered the most reliable way to obtain such data. However, due to developments made within the scientific community in recent years, numerical models are now able to achieve similar levels of reliability. To obtain the most accurate numerical predictions, many non-linear phenomena involved may be taken into account. These include: wave dispersion, wave–wave/current interaction, wave diffraction and interaction with the structure. When dealing with the development of a numerical method, one of the key points to take in to account for practical use is the computational work effort. Indeed, for many existing methods, this may prohibit the desired applications, for instance, of shape optimization or parametric studies. To support cost-efficient analysis with sufficient accuracy there is an important need for efficient computational models.

The motivation for this work stems from the dilemma in modeling the complex flow generated by wave-structure interaction. On the one hand, a large domain (and increased number of mesh elements as a result) is usually needed in numerical simulations so that the structure of interests is far away from the boundary to avoid unphysical effects. On the other hand, there is a super-linear dependence of the computational costs on the number of mesh elements.

The present study is dedicated to the development of an unsteady, heterogeneous, DD method for the simulation of wave–structure interactions. The entire fluid domain is subdivided into near field and far field subdomains in which the viscous and potential flow models are applied, respectively. The general idea of the DD method is illustrated in Figure 1. Between the two regions, the one encompassing the offshore object is designated as the near field, or $\Omega_{NSE}$. In this region, a comprehensive model for viscous flow based on NSE with free- surface capturing is used to reproduce the fine-scale drag, viscosity, vorticity, and turbulence near the structure. Since Navier-Stokes equation (NSE) is computationally expensive, it is not feasible to use the viscous flow model in large flow domains where the waves propagate. On the other hand, the potential flow (PF) model based on Potential Laplacian equation (PLE) is simple and thus used for the far field, or Ω_PF, where the effects of viscosity and vortices can be neglected, in order to take the advantage of relatively low computational cost. In the DD method, the coupling of models usually occurs across the interface between the two subdomains, and the form of the interface equations and their discretization is the crux of the coupled model. More generally, this heterogeneous model coupling is widely used when there is inherently different flow physics across interfaces, e.g., as in Fluid-Structure Interaction (FSI) where the fluid flow equations are coupled with solid mechanics models for structure deformation [3] or coupled aerodynamics and hydrodynamics problems in [4–8]. In this work, however, we pursue an overlapping DD approach and hence the coupling of models occurs within the overlapped regions, which is called buffer zones. More general nomenclature on DD methods can be found in [9].

This work aims to achieve the following goals, (i) allowing the interfaces ($\Gamma_{\text {C,NSE}}^{\mathrm{D}}$ and $\Gamma_{\mathrm{C}, \mathrm{PF}}^{\mathrm{N}}$) to intersect the free surface ($\Gamma_{\mathrm{f}}$) where the subscript letter C stands for Coupling interface, and to have an arbitrary shape; (ii) arbitrary number of regions; (iii) large displacement of the offshore structure. Thus, the heterogeneous coupled NSE-PF model poses several challenges. First, the fact that PLE is a simplification to NSE makes the formulation of a proper coupling model difficult. This is exacerbated by the numerical difficulties associated with the coupling of a traditional multiphase (air–water) formulation of NSE with the single-phase (water only) formulation of PF. Finally, the treatment of the free surface in the buffered zones has to be properly incorporated.

The above challenge (i) is largely due to an inherent incompatibility of PLE and NSE. In particular, since the PLE is based solely on mass conservation, the flow dynamics is missing and should be recovered using Bernoulli’s equation. An NSE-PF coupled model without a free surface was well-posed and reported in [10], but the associated modeling difficulties were also pointed out. A more sophisticated coupled NSE-PF model with free surface was presented in [11] and it assigned PF region to the bottom portion of the flow field, and NSE to the top part, thus the fixed artificial interface between PF and NSE subdomains does not intersect the free surface, which simplifies handling challenges (ii) and (iii). This strategy may not be easy to implement as it appears to be as the interface intersects the wave making and outlet boundaries instead, and more importantly, it is obviously not a good choice in terms of computational cost reduction considering the fact that the majority of mesh elements is clustered near the free surface.

This paper is arranged as follows. In Section 2 we present the mathematical formulations of the heterogeneous flow model together with the matching conditions at the interface. Numerical results from two examples compared to some physical experiments are reported in Section 3, which is followed by conclusions given in Section 4.

Throughout the paper, we use the following notation and assumptions for the flow models. We assume incompressibility of the fluid and, for simplicity, set density to one for air and one thousand for water. We assume also that the domain of flow $\Omega$, is an open-bounded Lipschitz domain with boundary $\partial \Omega$. The boundary $\partial \Omega$ of the domain $\Omega$ consists of the free-surface boundary $\Gamma_f$ at the air-water interface, the wavemaker boundary ($\Gamma_{I, P F}^N$ in Figure 1), the outlet boundary ($\tilde{\Gamma}_{O, P F}^N$), the atmosphere boundary ($\Gamma_{\text {atm, NSE}}^N$) and lastly no-slip boundary at the walls and bottom of the flow domain ($\Gamma_{W, N S E}^D$ and $\Gamma_{W, P F}^N$). We recall that $\Gamma_f$ is not a priori known that makes the exact geometry of time dependent.

The domain $\Omega$ is decomposed into overlapping subdomains $\Omega_{\text {NSE}}$ and $\Omega_{P F}$, and accordingly the free-surface boundary $\Gamma_f$ is further decomposed into $\Gamma_{f, P F}^D$ and $\Gamma_{f, N S E}$. The domain $\Omega_{\text {NSE}}$ is the viscous flow subdomain where flow is governed by NSE, and $\Omega_{P F}$ is governed by PLE. In addition, we put a letter $N$ in the boundary superscription to denote that a Neumann type boundary condition is employed at this particular boundary and a letter $D$ to denote a Dirichlet type boundary condition. For a variable $a$, we denote by $a_i=a \mid \Omega_i, i=N S E, P F$, and $a_{\Gamma}=\left.a\right|_{\Gamma}$, their restrictions to the subdomains and boundaries, respectively. In this section, the respective flow models will be discussed, which is followed by the coupling strategy. Next, we will then address the free-surface matching as well as the interpolation between nonconforming grids.

In $\Omega_{\text {NSE}}$, we conduct a two-phase flow simulation based on volume-of-fluid (VOF) model to solve the important variables velocity vector $\vec{U}$ and pressure $p$ and other fluid parameters include densities $\rho_l$ and $\rho_g$, viscosities $\mu_l$ and $\mu_g$ and surface tension coefficient. Here the subscript $l$ denotes the liquid phase and $g$ the gaseous phase. The flow is governed by the following equations [12]:

And

where $\Gamma_{f, N S E}$ is the air-water interface, $\delta\left(x-x_s\right)$ is the three-dimensional (3D). In this work, the continuum surface force (CSF) model of Brackbill et al. [13] is employed and as such the dynamic and kinematic boundary conditions on $\Gamma_{f, N S E}$ are incorporated into the system as a source term. Readers are referred to [13] for a detailed discussion of the integration of surface tension term and approximation of the local interfacial curvature $\kappa$. In the above expressions, the fluid density $(\rho)$ and viscosity $(\mu)$ are obtained by

and volume fraction function $\gamma$ is solved by a separate equation and is known at time level $t^n$ and $t^{n+1}$ when solving momentum equations, the density and viscosity fields in Eq. (3) are also known. The solver interDyMFoam in the framework of the finite volume method (FVM) is employed as a modified approach similar to the one proposed in [14] to solve the NSE flow problem.

The PF model which approximates NSE in $\Omega_{P F}$ makes the assumption that the viscosity, surface tension and compressibility are neglected. The fluid flow is thus irrotational so that a velocity potential, $\Phi(x, y, z, t)$, exists and the fluid velocity, $\vec{U}$, is given by its gradient, or $\vec{U}=\nabla \Phi$. This assumption reduces the mass conservation Eq. (2) to the Laplace equation:

Accordingly, the momentum conservation Eq. (1), after appropriate manipulations and integration in space, becomes Bernoulli's equation as described in [15]. It is noted that all the viscous terms disappear from Bernoulli's equation by virtue of Eq. (4), and this equation is only evaluated on some of the boundaries as the potential inside of $\Omega_{P F}$ is not needed.

The prescribed velocities are imposed on the external part of the boundaries to represent the wavemaker and wall boundaries [16]. This can be translated to a Neumann condition in $\Omega_{P F}$ on boundary $\Gamma_{I, P F}^N$.

where $\vec{U}^*$ is the boundary data which may be acquired from a wave theory [16] and $\vec{n}$ the unit normal pointing outward to the domain.

Furthermore, the following dynamic and kinematic free-surface boundary conditions are imposed on $\Gamma_{f, P F}^D$ [17]

Here $p^*$ is the atmospheric pressure (taken as constant zero) and $\frac{D}{D t}$ denotes the material derivative given as $\frac{\partial}{\partial t}+\vec{U} \cdot \nabla$. The source term $\Phi_b=g \eta$, where $g$ is the gravity acceleration as appeared in Eq. (1) and $\eta$ the surface elevation. In particular, Eq. (6) is the aforementioned Bernoulli's equation applied on the interface and needs to be solved for $\Phi$ which in turn will be assigned on $\Gamma_{f, P F}^D$ as the boundary condition. In this work, a second order accuracy backward Euler scheme is employed to discretize Eq. (6). The physical meaning behind Eq. (7) is that the net mass flux across free surface, $\Gamma_{f, P F}^D$, should be exact zero. This will be served as the criterion in adjusting the shape of the free-surface boundary $\Gamma_{f, P F}^D$. The entire mesh is then adapted to digest this boundary adjustment by using a vertex-based automatic mesh motion solver [18].

The model is finally complemented by the interface matching conditions, as well as the special treatment introduced for the buffer zone. Assume each of the NSE and PF solvers are provided with an entire set of external boundary conditions, and that each solver has its way of handling the free-surface boundary. The remaining condition needed is that on the coupling interfaces ($\Gamma_{C, N S E}$ and $\Gamma_{C, P F}$).

In the non-overlapping DD methods, $\Gamma_{C, N S E} \equiv \Gamma_{C, P F}=\Gamma_C$, and it is well known [9] that one needs two interface conditions on the coupling interface to account for the proper matching of primary variables and of their fluxes across an interface. These two conditions must involve three variables $\vec{U}, p, \Phi$, but only two of these have a physical meaning on each side of the interface, i.e., on each of the subdomains. These conditions can be derived from the matching of the velocities (thus mass) and by approximating the matching of stresses (thus momentum), respectively. The coupling procedure at the time step $t \rightarrow t+\Delta t$ can be described as the following:

1. The NSE is advanced in time from $t \rightarrow t+\Delta t$ by using $\left.\vec{U}_{P F}(t)\right|_{\Gamma}$ as a boundary condition on the interface $\Gamma_{C, N S E}^D$;

2. The VOF equation is solved in $\Omega_{\text {NSE}}$ to obtain a tentative volume fraction ( $\gamma$ ) field at the new time step;

3. The tentative velocity field at the new time step in $\Omega_{N S E}$ is used to get the normal velocity component $\left.\vec{U}_{N S E}\right|_{\Gamma} \cdot \vec{n}$ at $\Gamma_{C, P F}^N$;

4. The PLE is solved in $\Omega_{P F}$;

5. Adjust the free-surface boundary $\Gamma_{f, P F}^D$ such that Eq. (7) is satisfied, and morph the $\Omega_{P F}$ domain mesh accordingly;

6. The tentative velocity field at the new time step in $\Omega_{P F}$ is evaluated on $\Gamma_{C, N S E}^D$ and used as boundary conditions for the NSE solver at the next iteration;

7. Iterations are repeated till a convergence criteria is satisfied.

It shall noted that a free-surface matching procedure is usually required in Step (5) to enforce the continuity of surface elevations between PF and NSE regions, and this special procedure will be introduced next.

In this section, we apply the overlapping DD method to simulate the propagation of a solitary wave. This kind of shallow water wave is frequently used in ocean or coastal engineering research, particularly pertaining to the study of tsunamis. The computational domain is comprised of two overlapping sub-regions. As shown in Figure 2, the left half domain is $\Omega_{P F}$, and the right half domain is $\Omega_{N S E}$. The buffer zone $\Omega_{\text {Buff}}$ is located at the middle. At the beginning of the simulation, a soliton is generated in $\Omega_{P F}$ using a wave generation tool waves2Foam [19]. Such tool implements a large range of wave theories and among which, the solitary wave theory is implemented based on the limit solution for infinite wave length in the KdV equation, which is also the one describing the solution for regular cnoidal waves. The outputs of waves2Foam are the time history of velocity components and free-surface elevation. The resulting velocity data ($\vec{U}^*$) is then applied at the wavemaker boundary ($\Gamma_{I, P F}^N$) using eqn (5) and the free-surface elevation is used to move the free-surface boundary $\Gamma_{f, P F}^D$ near $\Gamma_{I, P F}^N$.

As the soliton propagates, it enters $\Omega_{\text {NSE}}$ through $\Omega_{\text {Buff }}$. Let $n_{1, h}$ and $n_{2, h}$ denote the number of discretization at $\Gamma_{C, N S E}^D$ from $\Omega_{N S E}$ and $\Gamma_{C, P F}^N$ from $\Omega_{P F}$, respectively. Thus, $n_{1, h}$ indicates the number of finite volume faces at $\Gamma_{C, N S E}^D$, whereas $n_{2, h}$ gives the number of finite volume faces on $\Gamma_{C, P F}^N$. In this work, we have $n_{2, h}=35$ and $n_{1, h}=120$ and this translates to $\Delta Z_{P F} \approx 3.43 \Delta Z_{N S E}$. In the horizontal direction, discretization size ratio is $\Delta X_{P F} \approx 2 \Delta X_{N S E}$. A PF simulation and an NSE simulation for the whole domain were also performed for comparison, results of which are hereafter referred to as "PF results" and "NSE results" respectively.

The computing CPU time for different methods is listed in Table 1 to investigate the efficiency of the Decomposed Domain method coupling PF and NSE. It can be found that coupling PF and NSE can reduce the computing CPU time to one third of that for the pure NSE method. The reduction of CPU time is associated with the fact that part of the domain is solved using fast tool PF instead of NSE.

The time history of free-surface elevation at four wave gauges is plotted in Figure 3, and it can be seen that although all the solvers produced almost the same accurate results at gauge 1, which is the closest one to the wavemaker, for all other gauges further away, the NSE solver suffers the most from the numerical diffusion as its predicted wave profile decreased a lot compared to PF and coupled DD method.

Snapshots of wave profile given in Figure 4 show a general agreement between coupled DD results and PF benchmark results. It is interesting to see that comparing to the analytic solution, the soliton travels faster in the NSE method while a bit slower in PF/coupled DD method. This somehow explains why the energy dissipates faster and thus wave height is lower from the NSE results.

The consistency of the coupled DD method can also be examined by looking at the nodal velocity values in/outside the buffer zone. Four fixed nodal probes, A-D, as shown in Figure 2, have been configured for this purpose. Figure 5 shows velocity time history at these probes. In particular, note the deviation of the predicted velocity values to the analytic solution near the end of the simulations is due to the numerical reflection and the fact that probes C and D

are close to the right wall of ΩNSE. Other than examining the wave profile and nodal velocity, we also compare the crest propagation celerity from the various results. As shown in Figure 5, the propagation of the wave generally matches very well. A slight phase shift occurs after tC λ =/.085 as the wave travels in slightly different speed in the various sub-domains as we have already pointed out.

In this paper, we presented our recent CFD work on coupling potential flow solver and vis-cous flow solver to study wave dynamics using open source code OpenFOAM. It helps to save computational resources especially for the multi-scale problems where large domain is needed. In this method, the computational domain is decomposed into two parts, the near-field subdomain and the far-field subdomain. In particular, only the field near offshore structures is calculated using NSE, while the field far from offshore structures are handled using PLE which can be solved faster than using NSE. The capability of the tool based on coupling models was validated by investigating the propagation of solitary wave in a numerical wave flume. The accuracy of the numerical tool in terms of wave profile, wave height and wave crest travelling speed was found to be in acceptable level. With this numerical tool, the offshore structure with big displacement can be studied with much less computation cost with the help of coupling strategy and further development of overset mesh to be developed in the future.