# On CFD Numerical Wave Tank Simulations: Static-Boundary Wave Absorption
Enhancement Using a Geometrical Approach

Muhannad W. Gamaleldin Corresponding Author.
<EMAIL_ADDRESS>
<EMAIL_ADDRESS>Department of Infrastructure Engineering, The
University of Melbourne, Victoria 3010, Australia.  Mechanical Engineering
Department, Alexandria University, Alexandria 21544, Egypt.  Alexander V.
Babanin  Department of Infrastructure Engineering, The University of
Melbourne, Victoria 3010, Australia.  Laboratory for Regional Oceanography and
Numerical Modeling, Qingdao National Laboratory for Marine Science and
Technology, Qingdao 266237, China.  Amin Chabchoub  Centre for Wind, Waves and
Water, School of Civil Engineering, The University of Sydney, NSW 2006,
Australia.  Marine Studies Institute, The University of Sydney, NSW 2006,
Australia.

###### Abstract

The present study aims to extend the applicability of the static-boundary
absorption method in phase-resolving CFD simulations outside the conventional
shallow-water waves limit. Even though this method was originally formulated
for shallow-water waves based on the conventional piston type wavemaker,
extending its use to deeper water conditions provides a more practical and
computationally cost efficient solution compared to other available numerical
wave absorption alternatives. For this sake, absorption of unidirectional
monochromatic waves in a semi-infinite flume by means of a static wall is
investigated theoretically and numerically. Moreover, implementation to a
practical wave-structure interaction application is investigated numerically
and experimentally. A phase-resolving numerical model based on the Reynold-
averaged Navier-Stokes (RANS) equations is implemented using the open source
C++ toolbox OpenFOAM®. The study presents the performance of the static-
boundary method, in a dimensionless manner, by limiting the depth at which the
active-absorption conventional-piston velocity profile is introduced; as a
function of incident wave conditions. Moreover, it is shown that the
performance of the static-boundary method can be significantly enhanced where
wave reflection was reduced to about half of that of the conventional setup in
deep-water conditions. Furthermore, the absorption depth is correlated to the
incident wave conditions; providing an optimization framework for the
selection of the proper dimensions of an absorbing wall. Finally, wave-
structure interaction experimental tests were conducted to validate the
numerical model performance; which shows an acceptable agreement between the
model and the experimental observations. The proposed limiter is straight
forward to be applied in pre-existing wave-structure interaction CFD solvers,
without the need of code modifications.

##  1 Introduction

With the depletion of the relatively easily accessed resources inland and
nearshore, marine industries have been marching further into deeper waters.
This, in turn, led to increasing attention of the scientific community toward
ocean engineering applications taking place in deep-water conditions.
Generally speaking, phase-resolving CFD numerical simulations of ocean
engineering applications are inherently expensive for a multitude of reasons
such as interface tracking methods and the wide scale of physical processes
that need to be resolved. However, one major challenge is the spatial
boundedness of a numerical domain, in which waves need to be artificially
absorbed at the domain bounds to mimic their corresponding real-life
situations. In other words, wave reflection off numerical domain bounds ought
to be prevented to avoid adversely affecting tested subjects. Consequently, a
number of wave absorption techniques have been devised in literature; which
can be classified into _internal_ , _dynamic-boundary_ and _static-boundary_
methods (Windt et al., 2019). First for the internal methods, also referred to
as _passive_ methods, the wave motion is absorbed from the computational
domain by dedicated zones inside the domain rather than using its boundaries;
hence the name _internal_. This can be done, for instance, by geometrically
(Finnegan and Goggins, 2012; Magee et al., 2015) or numerically (Anbarsooz et
al., 2014; Chen et al., 2014) modifying an absorbing zone to imitate a
_dissipative beach_ or _sponge layers_. Another common instantiation is the
use of a _relaxation zone_ (Fuhrman et al., 2006; Jacobsen et al., 2012; Hu et
al., 2016) where the numerical solution is gradually blended to a desired wave
motion and free surface profile. Overall, internal absorption methods are
relatively simple and straight-forward to apply, compared to other methods.
However, they increase the numerical expense due to their need to allocate
relatively large domain zones for optimum wave absorption with minimal
reflection; specially for longer waves. For instance, in (Deng et al., 2015;
Afshar, 2010) it was needed to allocate a three-wavelength long zone to
produce stable absorption.

Second, the dynamic-boundary methods where waves are absorbed using the
numerical domain boundaries (Lara et al., 2010; Troch and De Rouck, 1999);
assimilating experimental _active wave absorption_ systems (Milgram, 1970).
This is done by monitoring incident waves to the domain’s absorbing
termination/wall and moving it correspondingly, hence the name _active_ , in
such a way complying with the wave motion and preventing reflection off the
wall (Schäffer and Klopman, 2000). This in turn gives rise to the need of
using computational techniques to simulate the dynamic walls (Mittal and
Iaccarino, 2005), posing a considerable numerical toll to the problem’s
overall computational cost.

Finally, the static-boundary absorption method, which is the scope of the
present study, where the concept is very close to the previous type except
that the boundary is stagnant. This is simply because, in the numerical realm,
there is no need to physically move a wall to generate or absorb waves;
instead, motion of a wave absorbing wall is _modeled_ by artificially imposing
the desired velocity profiles on the stationary boundary as a Dirichlet
boundary condition. Consequently, this method is distinguished being the one
with the least computational cost amongst other available methods and a
practical choice in numerically demanding simulations (Babaei et al., 2017).
However, this absorption method suffers a severe limitation being applicable
only to shallow water conditions due to the use of the conventional piston
wavemaker in its formulation, which assumes a uniform velocity profile over
the water depth. Clearly, this leads to restricted performance in deeper water
conditions, majorly, due to the mismatch between the shallow and deeper water
wave particle kinematics. For instance in (Higuera et al., 2013), it was
reported that $11.2\%$ of wave amplitude was reflected when this method was
used just about $13\%$ inside the intermediate-water condition range. Some
alleviations have been proposed such as the combination with mesh stretching
or other passive method (Windt et al., 2019), free surface imposition (Torres-
Freyermuth et al., 2010), or even completely neglecting the piston profile and
replacing it with an analytical expression based on a water wave theory
(Higuera, 2020). For more detailed reviews about available techniques, readers
are referred to (Windt et al., 2019; Schmitt and Elsaesser, ; Miquel et al.,
2018).

In the present study, however, revisiting the classical wavemaker theory, one
can see that a wavemaker with a reduced displacement toward the bottom, such
as the flap wavemaker, obviously has better performance in deeper water
conditions compared to the conventional piston. This is because the flap
wavemaker kinematic velocity profile better matches the deep water particle
kinematic velocity profile, compared to the uniform velocity profile posed by
the piston wavemaker (Dean and Dalrymple, 1991). Interestingly, this remark,
however, has not been investigated in the aforementioned implementations of
active wave absorption techniques in CFD numerical simulations. The present
study aims to extend the range of applicability of the static-boundary
absorption method beyond the conventional shallow-water limit. In what
follows, the classical problem of wave absorption by a wall is revisited and
investigated from a hydrodynamical perspective. This is done by modifying the
conventional piston with a _limited absorption depth_ $\alpha h$ to better
match the wave kinematics in deeper water conditions. Moreover, this depth is
correlated to the incident wave conditions, providing an optimization
framework for the selection of the proper piston dimensions. Additionally,
implementation in two and three dimensional phase-resolving CFD numerical
models is presented and validated against theory and physical experiments;
respectively. It is worth highlighting that the proposed limiter is easy and
straight forward to be applied in pre-existing numerical models/packages
without any code modifications. This work is an in-depth analysis and
benchmarking for the author’s presented work at (Gamaleldin and Babanin, 2018,
2019, 2020).

##  2 The wavemaker theory

In what follows, the interaction of a reciprocating wall with a _linear_
unidirectional long-crested monochromatic wave is investigated analytically;
using the potential flow theory. In a numerical context, this wall is merely a
time-dependent velocity profile defined as a Dirichlet boundary condition. As
shown in Fig. 1, the problem at hand is to study the influence of the time-
dependent velocity profile on an incident wave train propagating toward the
wall at which velocity profile is introduced. If the frequency of the velocity
profile is chosen to be equal to that of the incident wave, the problem at
hand can be related to the concept of _mechanical wave-absorbers_ , also known
as “active water-wave absorbers” by (Milgram, 1970); where energy of an
incident wave train may be absorbed by a mechanical paddle. The underlying
concept is explained by the following steps. First, the horizontal velocity
profile $u(0,z,t)$ at the wall can be re-expressed as follows:

$\displaystyle u(0,z,t)=u(z)\cos(\omega t)$ (1)

Where $\omega$ is the frequency of oscillation of the introduced velocity
profile. After that, since the velocity profile is periodic, it might be
replaced by moving the wall instead to generate the same velocity profile as
follows:

$\displaystyle u(z)\cos(\omega t)=\frac{S(z)}{2}\omega\cos(\omega t)$ (2)

Where $S(z)$ is the corresponding stroke profile of the reciprocating wall.
Therefore, if the wall moves in such a way it becomes invisible to the
incident wave; the wave energy is completely absorbed by the wall with no
reflections. In other words, the wall is moved by the incident waves and
complies to their motion (Milgram, 1970). Thereby, what happens here is, more
or less, the opposite of the classical mechanical wavemaker problem; in which
the motion of the wall causes the wave (Dean and Dalrymple, 1991; Ursell et
al., 1960). As a result, in the present context, the interaction of different
subsurface profiles with incident waves is modelled as a mechanical wavemaker
problem; linked together by the kinematic boundary condition Eq. (2).

Figure 1: Schematic diagram of the interaction of a Dirichlet velocity profile
with an incident wave train. Figure 2: Schematic diagram of the proposed step
wavemaker. Figure 3: Wave-height-to-stroke ratio ($H/S$) versus relative depth
($kh$) for various values of $\alpha$. The label “shallow” indicates the
simplified theory.

The proposed wavemaker design is the _step wavemaker_ , shown in Fig. 2. In
comparison to the conventional piston wavemaker, the underlying idea is to
examine the influence of introducing the velocity profile at a portion of the
water depth ( limited absorption depth $\alpha h$) to better match the deep-
water waves particle velocity profile and, as a result, reduce reflections. In
what follows, the wave-height-to-stroke ratio is derived using the simplified
and full wavemaker theories.

###  2.1 Simplified theory

In the present section, the relation between a mechanical wavemaker and
resulted waves is derived based on a simplified assumption. The underlying
concept is to assume that displaced water by a single stroke of a wavemaker is
equal to accumulated water in the resulted wave crest; which is shown to be a
valid approximation in the shallow-water conditions (Galvin Jr, 1964). The
relation between the wavemaker stroke and the resulted wave height ($H/S$) may
be expressed by equating the displaced water by a full stroke to the volume
under the wave crest as shown in Fig. 2, as follows:

$\displaystyle S\alpha h=\int^{\lambda/2}_{0}\frac{H}{2}\sin(kx)dx$ (3)
$\displaystyle S\alpha h=\frac{H}{k}$ (4)
$\displaystyle\therefore\frac{H}{S}\Big{|}_{\text{step}}=\alpha kh$ (5)

Where the expression Eq. (5), plotted in Fig. 3, is valid for relative depths
of the range $kh<\pi/10$; i.e. the shallow-water conditions. Clearly, when
$\alpha=1$, the step wavemaker reduces to Eq. (6); the simplified conventional
piston wavemaker by (Galvin Jr, 1964).

$\displaystyle\frac{H}{S}\Big{|}_{\text{piston}}=kh$ (6)

###  2.2 Complete theory

Assuming that generated waves are small-amplitude, long-crested and
propagating in a potential, incompressible continuum; the governing equation
for the velocity potential is expressed as follows:

$\displaystyle\frac{\partial^{2}\phi}{\partial
x^{2}}+\frac{\partial^{2}\phi}{\partial y^{2}}=0$ (7)

Which is the Laplace equation boundary value problem. The boundary conditions
are:

$\displaystyle\eta=\frac{1}{g}\frac{\partial\phi}{\partial t}\rightarrow$
$\displaystyle\text{ at }z=0$ (8) $\displaystyle-\frac{\partial\phi}{\partial
z}=\frac{\partial\eta}{\partial t}\rightarrow$ $\displaystyle\text{ at }z=0$
(9) $\displaystyle-\frac{\partial\phi}{\partial z}=0\rightarrow$
$\displaystyle\text{ at }z=-h$ (10)

Where equations (8), (9) and (10) are the dynamic and kinematic boundary
conditions at the free surface and the impermeable sea bed condition;
respectively. Moreover, the waves considered are propagating in the positive
$x$-direction with no reflection from that end; as shown in Fig. 2.

Now, what is left is to satisfy the kinematic boundary condition on the
wavemaker wall. To address this, first, the horizontal displacement of a
wavemaker of a stroke profile $S(z)$ is described as follows:

$\displaystyle x=\frac{S(z)}{2}\sin(\omega t)$ (11)

Therefore, the wavemaker wall can be described by the following surface
function $F(x,z,t)$:

$\displaystyle F(x,z,t)=x-\frac{S(z)}{2}\sin(\omega t)=0$ (12)

Then, the kinematic boundary condition on the wavemaker wall is expressed as
follows:

$\displaystyle U.\hat{n}=\frac{-\partial F/\partial t}{|\nabla F|}$ (13)
$\displaystyle u-w\frac{\sin(\omega
t)}{2}\frac{dS(z)}{dz}=\frac{S(z)}{2}\omega\cos(\omega t)$ (14) $\displaystyle
u=\frac{S(z)}{2}\omega\cos(\omega t)\rightarrow\text{ at
}x=\frac{S(z)}{2}\sin(\omega t)$ (15)

Finally, if the wavemaker displacement is considered to be relatively small,
the kinematic boundary condition at the wavemaker can be linearized to be:

$u(0,z,t)=\frac{S(z)}{2}\omega\cos(\omega t)\rightarrow\text{ at }x=0$ (16)

Now, after defining the boundary conditions, the general form of Laplace
equation’s solution satisfying the upper and lower boundary conditions is
expressed as follows (Havelock, 1929):

$\phi=A_{p}\cosh(k_{p}(h+z))\sin(k_{p}x-\omega
t)+\sum_{n=1}^{\infty}C_{m}e^{-k_{s(n)}x}\cos(k_{s(n)}(h+z))\cos(\omega t)$
(17)

Where the subscripts $p$ and $s$ stand for _progressive_ and _standing_ waves
generated by the wavemaker, respectively. Those standing waves, also referred
to as _evanescent_ modes (Schäffer and Klopman, 2000), do not propagate and
are rather locked at the vicinity of the wavemaker; decaying exponentially in
the $x$-direction and are negligible two to three water depths of the
wavemaker (Ursell et al., 1960). However, the coefficients $A_{p}$ and $C_{n}$
are found by applying the kinematic boundary condition at the wavemaker wall;
and therefore, are dependent on the type of the wavemaker. Therefore, at
$x=0$:

$\displaystyle u\big{|}_{\text{wavemaker}}$
$\displaystyle=u\big{|}_{\text{flow field}}$ (18)
$\displaystyle\frac{S(z)}{2}\omega\cos(\omega t)$
$\displaystyle=-\frac{\partial\phi}{\partial x}$ (19)

It can be shown that (Dean and Dalrymple, 1991):

$\displaystyle A_{p}$ $\displaystyle=-I_{1}/I_{2}\text{ , where:}$ (20)
$\displaystyle I_{1}$
$\displaystyle=\int^{0}_{-h}\frac{S(z)}{2}\omega\cosh(k_{p}(h+z))dz$ (21)
$\displaystyle I_{2}$
$\displaystyle=k_{p}\int^{0}_{-h}\cosh^{2}(k_{p}(h+z))dz$
$\displaystyle=\frac{1}{4}\sinh(2k_{p}h)+\frac{k_{p}h}{2}$ (22)

Finally, the progressive wave height, relatively away from the wavemaker, can
be linked to the wavemaker type as follows:

$\displaystyle\eta=\frac{H}{2}\cos(k_{p}x-\omega t)$ (23)
$\displaystyle\frac{1}{g}\frac{\partial\phi}{\partial
t}\Big{|}_{z=0}=\frac{H}{2}\cos(k_{p}x-\omega t)$ (24)
$\displaystyle\therefore H=\frac{-2A_{p}}{g}\omega\cosh(k_{p}h)$ (25)

As such, by substituting an expression for $A_{p}$, which is dependent on the
wavemaker shape, into Eq. (25); the _wave-height-to-stroke_ ratio ($H/S$) can
be found. An expression $S(z)$ for the suggested step wavemaker, shown in Fig.
2, is proposed to be as follows:

$\displaystyle S(z)=S(\mathds{1}(z+\alpha h)-\mathds{1}(z))$ (26)

Where $\mathds{1}(z)$ is the Heaviside unit-step function, i.e.:

$\displaystyle\mathds{1}(z)\begin{cases}0&\text{for }z<0\\\ 1&\text{for
}z\geqslant 0\end{cases}$ (27)

Then, substituting into Eq. (21) we get:

$\displaystyle I_{1}=\int_{-h}^{-\alpha h}0dz+\int_{-\alpha
h}^{0}\frac{S}{2}\omega\cosh[k_{p}(h+z)]dz$ (28) $\displaystyle\therefore
I_{1}=\frac{S\omega}{2k_{p}}[\sinh(k_{p}h)-\sinh(k_{p}h(1-\alpha))]$ (29)

Now, substituting into Eq. (20) we get:

$\displaystyle
A_{p}=\frac{-I_{1}}{I_{2}}=\frac{-S\omega[\sinh(k_{p}h)-sinh(k_{p}h(1-\alpha))]}{2k_{p}[\frac{1}{4}\sinh(2k_{p}h)+\frac{k_{p}h}{2}]}$
(30)

Finally, the wave-height-to-stroke ratio is obtained by substitution into the
free surface expression, Eq. (25), as follows:

$\displaystyle\frac{H}{2}=\frac{-A_{p}}{g}\omega\cosh(k_{p}h)$ (31)
$\displaystyle\frac{H}{S}\Big{|}_{\text{step}}=\frac{4\sinh(kh)[\sinh(kh)-\sinh(kh(1-\alpha))]}{\sinh(2kh)+2kh}$
(32)

Equation (32) is plotted in Fig. 3, for several values of $\alpha$. Again, we
can see in the figure that when $\alpha=1$, the step wavemaker reduces to the
conventional piston wavemaker; described by the following expression (Dean and
Dalrymple, 1991):

$\displaystyle\frac{H}{S}\Big{|}_{\text{piston}}=\frac{2(\cosh(2k_{p}h)-1)}{\sinh(2k_{p}h)+2k_{p}h}$
(33)

###  2.3 Static-wall absorption

In what follows, implementation of the previous analysis as a Dirichlet
boundary condition on the wave absorbing wall is presented. Going back to the
original form of the problem, shown earlier in Fig. 1, the step wavemaker can
be transformed into a subsurface velocity profile $u(0,z,t)$ by substituting
the proposed step wavemaker profile Eq. (26) into the kinematic boundary
condition Eq. (16) as follows:

$\displaystyle u(0,z,t)=\frac{S}{2}(\mathds{1}(z+\alpha
h)-\mathds{1}(z))\omega\cos(\omega t)$ (34)

Where $S$ here might be thought of as an amplitude which is dependent on the
incident wave height $H$ and the relative water depth $kh$.
Thereby, the instantaneous velocity profile of the subsurface velocity profile
can be now linked to the instantaneous free surface elevation $\eta(0,t)$,
measured at the wall, by substituting the $H/S$ expression (32) into Eq. (34)
we get:

$u(z,t)=\omega\frac{H}{2}\cos(kx-\omega t)(\mathds{1}(z+\alpha
h)-\mathds{1}(z))\frac{\sinh(2kh)+2kh}{4\sinh(kh)[\sinh(kh)-\sinh(kh(1-\alpha))]}$
(35)

Therefore:

$u(z,t)=\omega(\mathds{1}(z+\alpha
h)-\mathds{1}(z))\eta(0,t)\frac{\sinh(2kh)+2kh}{4\sinh(kh)[\sinh(kh)-\sinh(kh(1-\alpha))]}$
(36)

Where, Eq. (36) is the corresponding subsurface instantaneous velocity profile
of the step profile. As seen in the equation, once the free surface elevation
is monitored, the velocity profile needed to absorb the incident wave is
defined. According to the potential flow theory, if this profile is applied at
a wall, the incident wave will be absorbed; with no reflected components.
Practically speaking, reflected components will be present and need to be
corrected for; which is discussed in the following sections.

Last but not least, it is worth mentioning here that a simplified version of
the previous procedure is what is being used in the C++ library OpenFOAM;
specifically, the IHFOAM and OlaFOAM/OlaFlow packages developed by (Higuera et
al., 2013). This is done by substituting the shallow-water approximation of
H/S for the conventional piston wavemaker Eq. (6), into the kinematic boundary
condition Eq. (16) as follows:

$\displaystyle u(0,z,t)=\frac{\omega}{2}\frac{H}{kh}\cos(\omega t)$ (37)
$\displaystyle u(0,z,t)=u(0,t)=\frac{C}{h}\text{ }\eta(0,t)$ (38)

Where, $C$ is the incident wave celerity; calculated using the shallow-water
approximation $C=\sqrt{gh}$. Moreover, it could have also been calculated
using the dispersion relation. In an experimental context, however, it is also
worth highlighting here that monitoring the celerity of incident waves is a
challenging aspect, specially outside the shallow-water conditions and
dispersive waves, and a number of other techniques have been used in the
context of using discrete-time systems and digital filters (Frigaard and
Christensen, 1994; Andreas, 2006; Wellens and Borsboom, 2020); which is
outside the scope of the present study which investigates the problem within a
hydrodynamic scope.

##  3 Numerical model

The numerical package of choice in the present study is the _Open-source Field
Operations And Manipulations_ (OpenFOAM). Unlike commercial codes, OpenFOAM is
not the type of a black box where the user can modify any step of the solution
process by modifying the source code and it covers a wide range of academic
and industrial problems (liquid sprays, external flows, multiphase
flows…etc.). OpenFOAM is a C++ library utilizing _object-oriented-programming_
to present, discretize and solve computational partial differential equations
(PDEs) (Weller et al., 1998; Ope, 2018). OpenFoam CFD solvers proved strong
potentials and capabilities by numerous case studies in the literature. For
instance, in the present scope, correctly modelling of up to 8th-order
harmonics for regular waves over a submerged bar by (Morgan et al., 2011) and
(Morgan and Zang, 2011). Moreover, (Chen et al., 2014) evaluated the
performance of OpenFOAM, specifically the InterFOAM module, to simulate
interactions of non-linear waves with offshore structures; showing good
agreement with experimental observations. In what follows, details of the
numerical model implemented in the present study is presented; with this font
used as a semantic markup for OpenFOAM specific verbatim entries and
expressions.

Table 1: Finite volume discretization schemes listed in OpenFOAM’s fvSchemes dictionary. Transport Term | Discretization Scheme
---|---
divSchemes | div(phirb,alpha) $\rightarrow$ Gauss interfaceCompression
div(rhoPhi,U) $\rightarrow$ Gauss limitedLinearV 1
div(phi,alpha) $\rightarrow$ Gauss vanLeer
gradSchemes | Gauss linear
snGradSchemes | corrected
interpolationSchemes | linear
timeScheme | Euler
Table 2: Computational domain boundary conditions in OpenFOAM’s conventions.

Boundary |  | Boundary Condition
---|---|---
Field | $p$ | $U$ | $\gamma$
Inlet wall |  | fixedFluxPressure | waveVelocity | waveAlpha
top wall |  | totalPressure | pressureInletOutletVelocity | inletOutlet
Bottom wall |  | fixedFluxPressure | fixedValue | zeroGradient
Outlet wall $\alpha h$ |  | fixedFluxPressure |  waveAbsorption2DVelocity (Eq. (38)) | zeroGradient
Outlet wall $(1-\alpha)h$ |  | fixedFluxPressure | fixedValue | zeroGradient

###  3.1 Governing equations and discretization schemes

The numerical technique implemented in the present study is classified as a
_CFD phase-resolving_ model (Folley et al., 2012). This is because
displacement of the air-water interface is resolved by means of a sufficiently
fine computational grid in comparison with the wave length. Transport
equations are expressed for two inviscid incompressible immiscible fluids of
different phases: air and water. The governing transport equations are
expressed as follows:

$\displaystyle\nabla\cdot U=0$ (39) $\displaystyle\frac{\partial(\rho
U)}{\partial t}+\nabla\cdot(\rho UU)=-\nabla p+\rho f_{b}$ (40)
$\displaystyle\frac{\partial\gamma}{\partial t}+\nabla\cdot(U\gamma)=0$ (41)

Equations (39), (40) and (41) are the mass conservation, momentum transport,
and phase-fraction transport equations. Where, $U$ is the velocity vector at
any point in the numerical domain, $\gamma$ is the phase-fraction indicator,
$\rho$ is the fluid density, $p$ is the pressure and $f_{b}$ is the term that
includes the effect of body forces per unit fluid mass; specifically, the
gravity and surface-tension forces. The phase-fraction indicator is a piece-
wise function that equals zero and one in the computational domain cells
consisting of gas and liquid, respectively. Thereby, at the interface region
between the two fluids, the phase-fraction indicator can take values between
zero and one. Moreover, this change is, and ought to be, sharp to form a high
phase-fraction gradient. Furthermore, this sudden change/gradient at the
interface region needs to be numerically preserved and propagated to provide
realistic simulations; resembling the coexistence of two immiscible fluids in
the computational domain.

The simulation procedure starts by initializing the computational domain by
allocating each fluid type to its corresponding initial zone; hence, an
interface is introduced. After that, the governing equations are solved to
transport the fluid and preserve a sharp interface. A two-fluid model approach
is implemented (Rusche, 2003) which is based on defining the contributions of
each specific fluid to the velocity field $U$ on a phase-averaged basis, i.e.:

$\displaystyle U=\gamma U_{liquid}+(1-\gamma)U_{gas}$ (42)

The phase-averaged velocity concept facilitates the ability to algebraically
sharpen the interface regions, hence classified as algebraic-VOF, by means of
the so called “compression” term (Rusche, 2003; Berberović et al., 2009). The
compression term is added to the right hand side of the phase-fraction
equation Eq. (41) as follows:

$\displaystyle\frac{\partial\gamma}{\partial
t}+\nabla\cdot(U\gamma)+\nabla\cdot[U_{r}\gamma(1-\gamma)]=0$ (43)

Where, $U_{r}$ is the relative velocity vector between the two phases (i.e.
$=U_{liquid}-U_{gas}$). Even though, the new term in equation Eq. (43)
vanishes analytically111The compression term
$\nabla\cdot[U_{r}\gamma(1-\gamma)]$ vanishes in continuum formulation as
$\gamma$ equals either $0$ or $1$; the term exists only at the interface
region in numerical formulation as $\gamma$ forms a gradient between $0$ and
$1$., it has a fundamental numerical role by compressing the interface region
into a couple-of-cells thickness (Berberović et al., 2009). This method is
broadly used in two-phase simulations in general and in sea waves simulations
in specific; solved using the “multidimensional universal limiters for
explicit solution” (MULES) solution algorithm, guaranteeing boundedness of the
solution (Deshpande et al., 2012). The OpenFOAM numerical solver of choice is
the olaFlow solver which is an updated version of the the older IHFoam package
(Higuera et al., 2013); based on OpenFOAM’s well known interfacial flow solver
interFoam (Ope, 2018).

The selection and understanding of the discretization methods implemented
forms a fundamental role in the numerical simulation; specifically, it
directly influences solution fields conservativeness, boundedness and
Transportiveness (Versteeg and Malalasekera, 2007). Table 1 summarizes the
discretization schemes of choice, in OpenFOAM’s semantics, which are specified
in the finite-volume dictionary fvSchemes. Computational grid cell-center
values are interpolated to the cell’s face-centers using the linear
interpolation scheme, declared in the interpolationSchemes sub-dictionary, as
follows:

$\psi_{f}=0.5\left(\psi_{c}+\psi_{d}\right)$ (44)

Where $\psi$ indicates an arbitrary field variable; and the subscripts $f$,
$c$ and $d$ indicates centers of face, central-node and downstream-node,
respectively. Therefore, interpolating face values using a second-order
central differencing scheme. After that, calculations of the gradient terms
are declared in the gradSchemes to be Gauss linear; indicating that the
standard Gaussian integration of the finite volume method is used to calculate
the face-center values from the cell-center values, using the linear
interpolation scheme described earlier. Then the snGradSchemes sub-dictionary
is used to declare the evaluation method of gradient components in the
direction normal to a cell’s face from the aforementioned calculated cell-
based gradients:

$\nabla^{\perp}_{f}\psi=\hat{n}\cdot(\nabla\psi)_{f}$ (45)

Here, $\hat{n}$ is the unit normal vector to the the cell’s face. The
divSchemes is where the divergence schemes are declared for the convective
terms of the transport equations. The entry div(rhoPhi,U) resembles
$\nabla\cdot(\rho UU)$ term in Eq. (40); set to be Gauss limitedLinearV 1,
where face values are interpolated using the convection-dependent scheme
limitedLinearV belonging to the _Total Variation Diminishing_ (TVD) class
(Harten, 1984). This implies that the interpolation of the face values is
influenced by the flow direction, limiting toward the upwind scheme. This is
chosen instead of a central differencing one to reflect the transportivness
nature of the problem for the convection terms. Then the entry div(phi, alpha)
resembles $\nabla\cdot(U\gamma)$ term in Eq. (43) where the convection
dependent Van Leer discretization scheme is selected; which is a high
resolution nonlinear Fromm-based scheme (Van Leer, 1974)— unbounded second-
order accurate in space. The last entry for divergence terms is div(phirb,
alpha) resembling the fundamental compressions term
$\nabla\cdot[U_{r}\gamma(1-\gamma)]$ of Eq. (43). Here, it is worth mentioning
that the aforementioned MULES algorithm utilizes a delimiter $L_{M}$ which is
equal to $1$ and $0$ at the interface-zone and elsewhere, respectively. This
is done as a mean to switch between using the high resolution scheme at the
interface region and the straight-forward convection scheme elsewhere to
conserve computational cost (Rusche, 2003; Deshpande et al., 2012). Finally,
temporal derivatives are declared at the ddtSchemes subdictionary to be the
implicit first-order Euler scheme:

$\frac{\partial\psi}{\partial t}=\frac{\psi-\psi^{\circ}}{\Delta t}$ (46)

Where the superscript $\circ$ indicates the previous time-level solution.

###  3.2 Geometry and boundary conditions

A computational domain is constructed, as shown in Fig. 4. The domain is
$9~{}\operatorname{m}$, $3.8~{}\operatorname{m}$, and a unity in the x,z and y
directions; respectively. The numerical domain is descretized into a
structured multi-block mesh using the blockMesh(blo, ) utility in OpenFOAM.
Three grid levels where created: coarse, medium and fine; with $0.8$, $1.4$
and $4$ million cells, respectively. This is done to ensure proper grid
resolution where the number of cells per design wave height are $16$, $19$ and
$31$ for the coarse, medium and fines grids, respectively. Moreover, to
maintain a recommended value of grid refinement ratio $r\geq 1.3$ for the sake
of the numerical verification calculations presented later in the succeeding
section.

Figure 4: Schematic diagram for the numerical domain setup. Figure 5: A
graphical representation of the underlying concept of step wavemaker tuning to
meet a desired performance point.

For the temporal resolution, adaptive time step is used to limit the maximum
Courant–Friedrichs–Lewy (CFL) number in the numerical domain to $0.15$.
Moreover, since the model used in the present study is a phase resolving one,
the maximum CFL number of the air-water interface is also limited to a maximum
value of $0.15$. Finally, the simulations where conducted to $60$ seconds of
waves propagation.

Table 2 summarizes the boundary conditions of choice; in OpenFOAM’s
conventions. Since the theoretical analysis in the previous section tackles
linear monochromatic waves, the simulated waves are selected to be so. Waves
are generated by the _inlet wall_ using the Stokes 1st-order theory (Dean and
Dalrymple, 1991). The generated waves at the inlet wall are
$0.01~{}\operatorname{m}$ high, $1.7~{}\operatorname{s}$ period,
$4.5~{}\operatorname{m}$ long and $3~{}\operatorname{m}$ deep. After that,
waves are absorbed by the _outlet wall_ where the subsurface velocity profile
is introduced along the $\alpha h$ portion of the entire water depth $h$;
using the OlaFlow wave generation and absorption library. The _top wall_ and
_bottom wall_ are set to be atmospheric and solid walls; respectively.
Finally, _phase 1_ and _phase 2_ are set to be standard water and air;
respectively. Simulations where conducted on the University of Melbourne’s
high performance computer Spartan (Lafayette et al., 2016).

##  4 Results and discussion

Hereafter, interaction of monochromatic waves with an active absorbing wall is
investigated. Practically, incident waves will not be completely absorbed and
some reflection will take place. This, in turn, leads to a partial standing
wave field in the wave flume. As a result, the wave height will not be
constant; instead, it varies sinusoidally along the channel. Consequently,
wave reflection coefficient forms a fundamental relevant benchmarking quantity
in the present study to evaluate performance of the wave absorbers
investigated. The reflection coefficient can be found using a number of
techniques such as using a carriage-mounted wave-height transducer or using
the three gauges method by (Mansard and Funke, 1980). Nevertheless, since
tests in the present study is conducted in a numerical flume, the wave
envelope is measured along the channel and the wave amplitude reflection
coefficient $\epsilon_{r}$ is found as follows:

$\displaystyle\epsilon_{r}=\frac{H_{r}}{H_{i}}=\frac{|\eta|_{\text{max}}-|\eta|_{\text{min}}}{|\eta|_{\text{max}}+|\eta|_{\text{min}}}$
(47)

Where ${H_{r}}$ and ${H_{i}}$ are the reflected and incident wave heights,
respectively.

Figure 6: Dimensionless absorption depth ($\alpha$) versus relative depth
($kh$). Figure 7: Temporal evolution of the water-air interface during a $10$
seconds duration for various limited absorption depths $\alpha h$.

###  4.1 Active wave absorption enhancement

As illustrated earlier in the theoretical modelling section, the subsurface
velocity profile in the numerical package is calculated using Eq. (38); which
is valid only for shallow-water conditions because it assumes a uniform
particle kinematic velocity profile over the entire depth. In other words, if
a deep-water wave is considered, the value of $H/S$ calculated by the
numerical model, based on the simplified wavemaker theory, will be far from
the actual true value of $H/S$ of the complete wavemaker theory; as can be
seen in Fig. 3. So, the goal now is to correct the value of $H/S$ of the
numerical model at deep-water conditions. In what follows a simple alleviation
for this problem is proposed.

A work around this problem is to set the numerical model to activate the
velocity profile on a portion $\alpha h$ of the outlet wall instead of the
whole depth $h$; assimilating the _step wavemaker_ modelled earlier by
equations (5) and (32). The underlying idea is to make use of the step
wavemaker versatility as the potential flow solution predicts; i.e. for a
given wave condition $kh$, $\alpha$ can be varied to produce a desirable $H/S$
value. Moreover, the step wavemaker profile is expected to better match the
velocity profile of a deep-water wave kinematic compared to the shallow-water
uniform kinematic profile posed by the conventional piston type. This can be
illustrated as follows:

$\displaystyle\frac{H}{S}\Big{|}_{\text{Numerical}}=\frac{H}{S}\Big{|}_{\text{Exact}}$
(48) $\displaystyle\therefore\frac{H}{S}\Big{|}_{\text{Eq.
(\ref{eqn:stepWaveMakerShallow})}}=\frac{H}{S}\Big{|}_{\text{Eq.
(\ref{eqn:stepWaveMakerDeep})}}$ (49) $\displaystyle\alpha
kh=\frac{4\sinh(kh)[\sinh(kh)-\sinh(kh(1-\alpha))]}{\sinh(2kh)+2kh}$ (50)

Therefore, by solving Eq. (50), one can find the optimum value of $\alpha$ so
that the value of $H/S$ calculated by the numerical model (using the shallow-
water approximation Eq. (5)) comply with the corresponding $H/S$ given by the
complete wavemaker theory Eq. (32). Graphically, this can be conceived as
tuning the performance curves (i.e. $kh$ versus $H/S$ curves) of a step
wavemaker until a desired wave condition ($kh$) is met. Figure 5 illustrates
this concept where the value of $\alpha$ is varied in such a way that the
interception point (i.e. solution of Eq. (50)) falls in the desired $kh$
range. Equation (50) is deemed implicit and solved by trial and error as shown
in Fig. 6. Indeed, one can see in the figure that $\alpha$ is, more or less, a
unity in the shallow-water range.

For instance, for the wave conditions considered in the present study, the
relative wave depth is $kh=4.18$. By substituting into Eq. (50) we get
$\alpha=0.38$. This indicates that, for optimum wave absorption of the outlet
wall with minimum reflection, the subsurface velocity profile is to be
introduced along a depth $\alpha h\approx\lambda/4$. In what follows, a
parametric case study is conducted to investigate the influence of varying the
limited absorption depth $\alpha h$ on the performance of an active wave
absorption wall in deep-water conditions.

Table 3: Calculations of the GCI spatial discretization-error estimates. |  $\varphi=\epsilon_{r}$ (-) |  $\varphi=A_{k}$ (cm)
---|---|---
$r_{21}$ | $1.69$ | $1.69$
$r_{32}$ | $1.3$ | $1.3$
$\varphi_{1}$ | $0.136$ | $0.361$
$\varphi_{2}$ | $0.127$ | $0.368$
$\varphi_{3}$ | $0.113$ | $0.367$
$P$ | $2.85$ | $2.87$
$\varphi^{21}_{ext}$ | $0.139$ | $0.359$
$e^{21}_{a}$ | $7.22\%$ | $2.05\%$
$e^{21}_{ext}$ | $2.03\%$ | $0.58\%$
$GCI^{21}_{fine}$ | $2.6\%$ | $0.73\%$

###  4.2 Numerical model verification

Although analytical/theoretical solution is available to validate the
numerical results against, estimating discretization errors is still relevant
as numerical simulations are usually intended to simulate influences of
physical events even when such validation data are absent. Moreover, it is
also deemed important to argue whether the outcome of a numerical model is of
a physical behaviour or, otherwise, random. In the present study, numerical
uncertainty is investigated using the _Grid Convergence Index_ (GCI) method
(Celik et al., 2008; Roache, 1998); which is based on the _Richardson
Extrapolation_ (RE) method (Richardson, 1911; Richardson and Gaunt, 1927). The
GCI represents a measure of how much a computed value is distant off a
numerical model’s asymptotic value.

For this sake, the amplitude reflection coefficient ($\epsilon_{r}$) and the
amplitude spectral peak ($A_{k}$) are chosen as the global variables of issue
for the grid convergence test. Table 3 showcases the GCI calculations. It
shows that the apparent spatial order of convergence of the numerical solution
$P$ is found to be around $2.8$ for both variables. Moreover, values of GCI
for the fine-grid solution for $\epsilon_{r}$ and $A_{k}$ are $2.6\%$ and
$0.73\%$, respectively; indicating that the numerical results are relatively
within acceptable margin of the asymptotic solution.

Table 4: Wave statistical quantities for various limited absorption depths $\alpha h$. $\alpha h$ | $h$ | $\lambda/2$ | $\lambda/4$ | $\lambda/18$
---|---|---|---|---
$\epsilon_{r}(\%)$ | $28.1$ | $22.51$ | $12.66$ | $21.94$
$A_{k}(cm)$ | $0.543$ | $0.484$ | $0.369$ | $0.54$
$\sigma^{2}(cm^{2})$ | $0.178$ | $0.141$ | $0.083$ | $0.173$
$H_{s}(cm)$ | $1.69$ | $1.50$ | $1.15$ | $1.66$

###  4.3 Flow visualization and waves transmission

A parametric study is conducted to investigate the influence of varying the
limited absorption depth $\alpha h$ on the performance of a wave absorbing
wall. Four absorption depths $\alpha h$ are investigated: $h$, $\lambda/2$,
$\lambda/4$ and $\lambda/18$. Figure 7 shows the temporal locus of the air-
water interface during a $10$ second duration, where the interface is defined
numerically by the set of cells where the _phase fraction function_
$\gamma=0.5$.

For the first case $\alpha h=h$, this case resembles the conventional
numerical setting where subsurface velocity profile is introduced along the
whole depth. As have been addressed earlier in the theoretical analysis
section, the velocity profile is introduced in such a way that a conventional
piston wave absorber is used (i.e. equations (6) and (38)). And since this is
based on the shallow-water conditions approximation, high reflection is
anticipated from the absorbing wall. As seen in Fig. 7, reflection takes place
and the wave height varies sinusoidally along the wave flume. Indeed, one can
see that this is a manifestation of the classical _partial standing wave_
system. Bearing that in mind, the amplitude reflection coefficient
$\epsilon_{r}$ can be calculated using Eq. (47) and is found to be
$\epsilon_{r}=28.1\%$.

Figure 8: Reflection coefficient ($\epsilon_{r}$) versus limited absorption
depth ($\alpha h$). Figure 9: Wave amplitude spectra measured by a wave gauge
mounted at mid-flume; at a sampling frequency of $20$ Hz. Figure 10: Schematic
diagram for the numerical domain setup with an added damping zone.

The second case of choice is $\alpha h=\lambda/2$ since this value corresponds
to the effective depth of a linear monochromatic deep-water wave. As can be
seen in Fig. 7, the free surface envelope shows lower variations compared to
the previous case indicating lower reflection off the absorbing wall. The
reflection coefficient for this case is found to be $\epsilon_{r}=22.51\%$.
This is mainly attributed to the improved matching between the absorbing
velocity profile posed on the absorbing wall and the deep-water wave velocity
profile; compared to the previous case.

After that, the value of $\alpha h=\lambda/4$ which is predicted by the
theoretical solution of Eq. (50) to provide optimum absorption. Indeed, one
can see in Fig. 7 that this case provides the least variation in the wave
height along the wave flume. Hence, the best absorbing wall with reflection
coefficient $\epsilon_{r}=12.66\%$, which is more than $50\%$ reduction in
reflection compared to the conventional shallow-water approximation setup of
$\alpha h=h$. For benchmarking purposes, it is worth highlighting here that
the reflection coefficient value ranged between $3.2\%$ and $11.2\%$ in the
shallow-water regime using the conventional setup in (Higuera et al., 2013);
measured using the three gauges method (Mansard and Funke, 1980).

Last case considered is when $\alpha h=\lambda/18$ which seems to be
interesting for two reasons. First, it investigates the scenario when
absorption depth is set to be _lower_ than the optimum value predicted by Eq.
(50). Second, it is the solution of a modified version of Eq. (50) where the
wave number in the left hand side is substituted by:

$k=\frac{\omega}{C}=\frac{\omega}{\sqrt{g\alpha h}}$ (51)

Where $C$ and $g$ are the wave celerity and the gravitational acceleration,
respectively. This is done to take into consideration the method used in the
present numerical packages, such as IHFOAM and OlaFOAM/OlaFlow, to find the
wave celerity; as addressed earlier in the theoretical analysis section.
Substituting into Eq. (50) we get:

$\sqrt{\alpha
kh\cdot\tanh(kh)}=\frac{4\sinh(kh)[\sinh(kh)-\sinh(kh(1-\alpha))]}{\sinh(2kh)+2kh}$
(52)

Where solution to the Eq. (52) is $\alpha h\approx\lambda/18$, by trial and
error, for the wave conditions considered in the present study. One can see in
Fig. 7 that wave height variation has increased compared to the previous case;
the wave amplitude reflection coefficient is found to be
$\epsilon_{r}=21.94\%$. Figure 8 shows the behaviour of $\epsilon_{r}$ versus
$\alpha h$ for of all cases tested in the present work.

Figure 11: Reflection coefficient ($\epsilon_{r}$) versus damping coefficient
($\theta$). Figure 12: Temporal evolution of the water-air interface during a
$10$ seconds duration for two values of the damping coefficient $\theta$.

Finally, a wave gauge is installed mid-flume, i.e. at
$x=4.5~{}\operatorname{m}$, and Fig. 9 shows the amplitude spectra for all
cases. Moreover, table 4 shows the corresponding values of the spectral peaks
$A_{k}$. The table also shows the variance $\sigma^{2}$ and the significant
wave height $H_{S}$ of the wave gauge reading; defined as follows:

$\displaystyle\sigma$
$\displaystyle=\sqrt{\frac{1}{N-1}\sum_{n=1}^{N}\zeta^{2}_{n}}$ (53)
$\displaystyle H_{S}$ $\displaystyle=4\sigma$ (54)

Where, $\zeta_{n}$ is the sampled water surface elevation. Additionally, the
table also shows the values of the reflection coefficient $\epsilon_{r}$ using
Eq. (47).

###  4.4 Absorption performance evaluation

In the preceding sections, active wave absorption was significantly improved
following a hydrodynamical approach. Even though wave reflection has been
dropped significantly compared to the standard shallow-water approximation
setup, wave reflection is still relatively high and will adversely affect test
subjects placed in that flume. The remaining reflection can be attributed to a
number of reasons. First, the inherent nonlinear nature of the problem which
is resembled here by the use of a fully nonlinear phase-resolving CFD model;
in comparison to the linearized potential flow approach followed to derive Eq.
(50) and the assumption of monochromatic linear waves in the present analysis.
Indeed, one can see in Fig. 9 that waves have modulated over the spectrum by
the existence of wave components beside the monochromatic fundamental
frequency.

Second reason is the existing mismatch in the particle kinematic between the
velocity profile of the incident waves and the one posed on the absorbing
wall. For instance, the simulated waves in the present study are of an
exponential velocity profile (deep-water waves) which are absorbed by a step
profile on the absorbing wall Eq. (36).

Third, since wave reflection will inevitably take place in practical
applications, the value of $H$ substituted in the free surface expression
$\eta(x,t)$ used in the previous sections needs to be corrected for the
reflected components. In other words, take Eq. (38) for instance, the
subsurface velocity profile at the wall should be calculated based on the
incident component only rather than the combined incident and reflected ones.
For this sake, the method proposed by (Ursell et al., 1960) to filter out the
incident wave height in mechanical wavemaker setups is proposed as a remedy.
First, the wave height reflection coefficient $\epsilon_{r}$ is found from
measurements as addressed earlier using Eq. (47). After that, the incident
wave height is set to be $H_{i}=2A$ where:

$\displaystyle A=\frac{0.5H_{\text{avg}}}{1+\epsilon_{r}\cos(\delta)}$ (55)
$\displaystyle\text{If $\epsilon_{r}$ is relatively small }\rightarrow
A\approx 0.5H_{\text{avg}}$ (56)

Where: $H_{\text{avg}}$ is the average wave height over half wave length of
the incident wave and $\delta$ is the phase shift between the incident and
reflected components. Moreover, most of the present active wave absorption
numerical packages seem to disregard the existence of the evanescent modes at
the vicinity of the wall, addressed earlier in Eq. (17), where $\eta(x,t)$ is
measured on the absorbing wall itself. Experimentally, this is avoided by
monitoring $\eta(x,t)$ away from the wall where evanescent modes would have
diminished.

Finally, the fourth major reason is the calculation of the incident waves
celerity addressed earlier in the theoretical analysis section. For instance,
the numerical model implemented in the present study is based on shallow-water
approximation and the celerity is found correspondingly. It is worth
highlighting here that this issue might be resolved, as shown in (Higuera,
2020), by solving the dispersion relation iteratively.

Figure 13: Experimental setup photographs for: an Edinburgh Designs wavemaker
(a), a wave-absorbing beach with artificial grass (b), and a vertical rigid
cylinder fitted with a force transducer (c). Table 5: Wave-structure
interaction computational domain boundary conditions in OpenFOAM’s
conventions.

Boundary |  | Boundary Condition
---|---|---
Field | $p$ | $U$ | $\gamma$
Inlet wall |  | fixedFluxPressure | waveVelocity | waveAlpha
top wall |  | totalPressure | pressureInletOutletVelocity | inletOutlet
Bottom wall |  | fixedFluxPressure | slip | zeroGradient
Mid-Flume |  | fixedFluxPressure | slip | zeroGradient
Cylinder wall |  | fixedFluxPressure | slip | zeroGradient
Outlet wall $\alpha h$ |  | fixedFluxPressure | waveAbsorption2DVelocity | zeroGradient
Outlet wall $(1-\alpha)h$ |  | fixedFluxPressure | fixedValue | zeroGradient

Figure 14: Snap shots for the three dimensional numerical domain showing: a
vertical $x-z$ section (a), a horizontal $x-y$ section (c), and close-up views
(b). Figure 15:  A snap shot of the simulated wave passing by the vertical
cylinder, using the CFD numerical model (a rendered mirror view). Free surface
is represented by clipping the value of $\gamma=0.5$. Figure 16: Numerical
versus experimental free surface elevation, measured at a point placed
$\lambda/2$ upstream the cylinder (a), and a close-up view (b). Figure 17:
Numerical versus experimental inline forces on the tested cylinder (a), and a
close-up view (b).

However, even though the reflection coefficient has been dropped significantly
compared to the shallow-water approximation conventional ($\alpha h=h$) setup;
the remainder reflection ought to be further reduced to avoid inversely
affecting tested subjects to be placed in such a flume. Remainder reflection
can be reduced by adding a relatively small passive absorbing zone to amend
performance of the absorbing wall; resulting in a hybrid (i.e. a combination
of active and passive) approach. The underlying concept in such an approach is
that the static wall is intended to absorb the majority of incident waves and
the passive zone to absorb (or dissipate) the remainder wave components that
have been reflected of the static absorbing wall. This simply enables the
utilization of the best-of-both-worlds, where the numerical cost is
significantly reduced in comparison to pure passive absorption and the
reflection coefficient is being within acceptable limits in comparison to pure
active absorption in deep water conditions. Other instantiations of the use of
hybrid absorption approaches have been showcased in (Windt et al., 2019; Chen
et al., 2014; Higuera, 2020; Israeli and Orszag, 1981).

Consequently, a relatively small damping _zone_ of length $\lambda/2$ with a
gradually ramped damping coefficient is added next to the $\alpha h=\lambda/4$
absorbing wall; as shown in Fig. 10. This is done by adding an artificial
gradual vertical damping term to the momentum equation Eq. (40) such that:

$\frac{\partial(\rho U)}{\partial t}+\nabla\cdot(\rho UU)=-\nabla p+\rho
f_{b}-\rho\theta U_{z}\\\ $ (57)

Where $\theta$ is the damping coefficient and $U_{z}$ is the vertical velocity
component. A longer damping zone would result in less reflection but with the
considerable increase in the computational cost. In fact, in pure passive
absorption scenarios, it was recommended that the absorption zone should be at
least $\lambda$ to $2\lambda$ in length for acceptable wave absorption (Windt
et al., 2019; Chen et al., 2014; Miquel et al., 2018). Consequently, the
shorter damping zone was investigated in the present work to investigate the
effectiveness the implemented combination. On the other hand, however, the
damping coefficient seems to have an optimum value where relatively high
damping coefficient $\theta$ values would increase reflection as reported in
(Chen et al., 2014; Romate, 1992).

Figure 11 shows the behaviour of the reflection coefficient ($\epsilon_{r}$)
versus the damping coefficient $\theta$. The figure shows that the optimum
$\theta$ value is 3 where the reflection coefficient was reduced from
$12.66\%$ to $4.44\%$. However, raising the damping coefficient over that
resulted in adverse absorption performance. Moreover, Fig. 12 shows the
temporal locus of the air-water interface for two values of the damping
coefficient $\theta$.

##  5 Application to wave-structure interaction

In what follows, the methodology described in the aforementioned sections is
implemented in a typical engineering application; wave forcing of a rigid
fixed vertical cylinder.

###  5.1 Experimental setup

An experimental setup was constructed in the Wave Flume at the Fluids
Laboratory at the University of Sydney, Australia. The flume is
$30~{}\operatorname{m}$ long, $1~{}\operatorname{m}$ wide and
$1~{}\operatorname{m}$ deep; with a maximum water depth of
$0.75~{}\operatorname{m}$. The experimental setup, depicted in Fig. 13,
entails: an Edinburgh Designs piston-type wavemaker before power-up (a), a
sloped wave-absorbing beach with artificial grass at the opposite end (b), and
a vertical rigid cylinder fitted with a force transducer (c). The cylinder’s
height and diameter are $1.1~{}\operatorname{m}$ and $56~{}\operatorname{mm}$,
respectively. The force transducer is a multi-axis six degrees of freedom load
cell that measures forces and moments in the x, y and z directions. Finally, a
twin-wire resistive wave gauge is installed about $\lambda/2$ upstream of
tested cylinder.

###  5.2 Numerical model

A three-dimensional numerical model is constructed to model the interaction of
regular wave train with a fixed vertical cylinder. The implemented governing
equations and discretizations schemes are the same as those illustrated
earlier in section 3.1. On the other hand, the numerical domain is different
as shown in Fig. 14. As shown in the figure, flow symmetry is assumed around
the “Mid-Flume” boundary and, therefore, half of the domain is simulated.
Moreover, a multi-block structured mesh is implemented at the vicinity of the
cylinder ending up with a high quality mesh as shown in the figure. The
numerical domain is descretized in such a way that $500\text{ cells}/\lambda$
in both $x$ and $y$ directions, and a maximum refinement of $21\text{
cells}/H$ in the $z$ direction at the free surface to ensure acceptable
resolution and reduce numerical diffusion (Finnegan and Goggins, 2012;
Bredmose et al., December 2013; Paulsen et al., 2014). This resulted is a
total cell count of about $7.8$ million cells. Table 5 summarizes the boundary
conditions of choice, using OpenFOAM’s conventions.

###  5.3 Test conditions and results validation

Different wave conditions were considered in this section, to those of section
3.2, to investigate the validity of proposed approach at a different scenario.
The generated waves are $0.03~{}\operatorname{m}$ high,
$1.67~{}\operatorname{s}$ period, $3.7~{}\operatorname{m}$ long and
$0.75~{}\operatorname{m}$ deep. Clearly, this corresponds to a steepness
($ak$) of $0.026$ and an Ursell number ($U_{r}$) of $0.99$, which indicates
that the waves are slightly steeper than the proposed linear theory range of
validity in (Le Méhauté, 1976). Moreover, the relative depth of $kh=1.27$
which corresponds to the intermediate water condition. Furthermore, as seen in
table 5, slip boundary conditions have been implemented which indicate that
viscous effects and boundary layer formation are suppressed. This is still a
valid approximation here since the simulated wave conditions here corresponds
to a Keulegan–Carpenter number ($KC$) of $1.04$, which corresponds to the
inertia-dominated loading regime (Sumer and Fredsøe, 2006). Substituting into
Eq. (50), the proposed optimum limited absorption ratio will be $\alpha=0.83$.
Temporal resolution was varied in such a way that the maximum CFL number is
limited to $0.15$. Finally, the simulation was run for $10.5T$.

Figure 15 shows the wave-structure interaction from the three-dimensional CFD
model. Referring back to phase-fraction transport equation Eq. (41), the free
surface is defined by means of clipping numerical cells where the phase-
fraction $\gamma=0.5$. This illustrates the importance of opting for a high
resolution mesh at the air-water interface so that we end up with relatively a
sharp interface between the two material phases, rather than a smeared
$\gamma$ distribution. Figure 16 shows a comparison between the numerical and
experimental wave gauge measurements at a point placed $\lambda/2$ upstream of
the tested cylinder. It is observed that a very good agreement between both is
achieved, indicating proper absorption of the outlet wall. Finally, Fig. 17
shows a comparison between numerical and experimental inline forces on the
cylinder. Again, a very good agreement is observed between both records,
specially after filtering out noise and structural vibrations (Bredmose et
al., December 2013). This indicate proper effectiveness of the proposed
approach and its validity for engineering applications.

##  6 Concluding remarks

The present study aimed to extend the range of applicability of the static-
boundary absorption method outside the conventional shallow-water waves limit,
providing a computationally cost effective alternative to the other available
methods. To tackle this, absorption of unidirectional monochromatic waves in a
semi-infinite flume by means of a static wall was investigated theoretically
and numerically. This was done by proposing a limited absorption depth $\alpha
h$ which corresponds to the incident wave conditions, to better match the wave
kinematics in deeper water conditions compared to the conventional method of
the shallow-water approximated piston wavemaker.

A theoretical analysis was derived for the interaction of a monochromatic wave
with an absorbing wall using the wavemaker theory. The optimum absorption
depth $\alpha h$ was introduced and linked to the incident wave conditions
using the kinematic boundary condition on the absorbing wall; as shown in Fig.
6 and Eq. (50).

Moreover, a nonlinear wave phase-resolving CFD numerical model was implemented
to validate the proposed theoretical outcomes; using different case scenarios
of the limited absorption depth $\alpha h$. As can be seen in Fig. 7, wave
height variation along the flume was greatly influenced by the variation of
the limited absorption depth $\alpha h$. Compared to the conventional shallow-
water approximation setting ($\alpha h=h$), absorption was significantly
enhanced where wave reflection coefficient dropped from $\epsilon_{r}=28.1\%$
to $12.66\%$. Furthermore, the proposed solution is relatively easy and
straight-forward to be implemented to existing numerical packages without code
modifications. In addition, it can be validated experimentally because the
proposed design is, more or less, a variant of the conventional piston-type
wavemaker.

However, even though the reflection was significantly reduced, still it should
be reduced further more to avoid adversely influencing tested subjects in the
CFD numerical wave flume. This remainder reflection was attributed to a number
of reasons summarized in the following points:

  * •

The use of a linearized theoretical modelling to address the inherent
nonlinear nature of water waves.

  * •

The mismatch between the incident waves and wall absorbing profiles.

  * •

The inclusion of both the reflected waves and evanescent modes in the measured
free surface elevation $\eta(x,t)$.

Aside of the aforementioned caveats/limitations of this method, the remainder
reflection could be further significantly reduced by means of implementing a
relatively small damping zone resulting in a hybrid approach. This has
resulted in reducing wave reflection form $12.66\%$ to $4.44\%$, which is
comparable with typical physical model tests.

Finally, the proposed corrector has been applied and experimentally validated
to a typical engineering application: wave forcing of a vertical rigid
cylinder. As shown in figures 16 and 17, a very good agreement between the
numerical and experimental measurements was observed in both the free surface
profile and the inline forces record.

## Glossary

## References

  * Windt et al. [2019] Christian Windt, Josh Davidson, Pál Schmitt, and John V Ringwood.  On the assessment of numerical wave makers in cfd simulations.  _Journal of Marine Science and Engineering_ , 7(2):47, February 2019.  Multidisciplinary Digital Publishing Institute, doi:10.3390/jmse7020047.
  * Finnegan and Goggins [2012] William Finnegan and Jamie Goggins.  Numerical simulation of linear water waves and wave–structure interaction.  _Ocean Engineering_ , 43:23–31, January 2012.  Elsevier, doi:10.1016/j.oceaneng.2012.01.002.
  * Magee et al. [2015] Allan R Magee, Varjola Nelko, Kian Yew Lim, and Lup Wai Chew.  Benchmarking experiments for wave absorption modeling.  Newfoundland, Canada, May-June 2015.  American Society of Mechanical Engineers (ASME) 2015 34th International Conference on Ocean, Offshore and Arctic Engineering (OMAE), doi:10.1115/OMAE2015-42300.
  * Anbarsooz et al. [2014] Morteza Anbarsooz, Mohammad Passandideh-Fard, and Mohammad Moghiman.  Numerical simulation of a submerged cylindrical wave energy converter.  _Renewable Energy_ , 64:132–143, 2014.  Elsevier, doi:10.1016/j.renene.2013.11.008.
  * Chen et al. [2014] LF Chen, J Zang, AJ Hillis, GCJ Morgan, and AR Plummer.  Numerical investigation of wave–structure interaction using openfoam.  _Ocean Engineering_ , 88:91–109, 2014.  Elsevier, doi:10.1016/j.oceaneng.2014.06.003, ISSN 0029-8018.
  * Fuhrman et al. [2006] David R Fuhrman, Per A Madsen, and Harry B Bingham.  Numerical simulation of lowest-order short-crested wave instabilities.  _Journal of Fluid Mechanics_ , 563:415–441, September 2006\.  Cambridge University Press, doi:10.1017/S0022112006001236.
  * Jacobsen et al. [2012] Niels G Jacobsen, David R Fuhrman, and Jørgen Fredsøe.  A wave generation toolbox for the open-source cfd library: Openfoam®.  _International Journal for Numerical Methods in Fluids_ , 70(9):1073–1088, November 2012.  Wiley Online Library, doi:10.1002/fld.2726.
  * Hu et al. [2016] Zheng Zheng Hu, Deborah Greaves, and Alison Raby.  Numerical wave tank study of extreme waves and wave-structure interaction using openfoam®.  _Ocean Engineering_ , 126:329–342, 2016.  Elsevier, doi:10.1016/j.oceaneng.2016.09.017.
  * Deng et al. [2015] Yanfei Deng, Jianming Yang, Wenhua Zhao, Longfei Xiao, and Xin Li.  An efficient focusing model of freak wave generation considering wave reflection effects.  _Ocean Engineering_ , 105:125–135, 2015.  Elsevier, doi:10.1016/j.oceaneng.2015.04.058.
  * Afshar [2010] Mostafa Amini Afshar.  Numerical wave generation in openfoam®.  Master’s thesis, Technical University of Denmark, 2010.  Chalmers tekniska högskola.
  * Lara et al. [2010] Javier L Lara, Andrea Ruju, and Inigo J Losada.  Reynolds averaged navier–stokes modelling of long waves induced by a transient wave group on a beach.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 467(2129):1215–1242, November 2010\.  The Royal Society Publishing, doi:10.1098/rspa.2010.0331.
  * Troch and De Rouck [1999] Peter Troch and Julien De Rouck.  An active wave generating–absorbing boundary condition for vof type numerical model.  _Coastal Engineering_ , 38(4):223–247, December 1999.  Elsevier, doi:10.1016/S0378-3839(99)00051-4.
  * Milgram [1970] Jerome H Milgram.  Active water-wave absorbers.  _Journal of Fluid Mechanics_ , 42(4):845–859, March 1970.  Cambridge University Press, doi:10.1017/S0022112070001635.
  * Schäffer and Klopman [2000] Hemming A Schäffer and Gert Klopman.  Review of multidirectional active wave absorption methods.  _Journal of waterway, port, coastal, and ocean engineering_ , 126(2):88–97, March 2000.
  * Mittal and Iaccarino [2005] Rajat Mittal and Gianluca Iaccarino.  Immersed boundary methods.  _Annu. Rev. Fluid Mech._ , 37:239–261, January 2005.  Annual Reviews, doi:10.1146/annurev.fluid.37.061903.175743.
  * Babaei et al. [2017] Hossein Babaei, Scott Baker, and Andrew Cornett.  Validation of an open-source cfd tool to support efficient design of offshore gravity-based structures exposed to extreme waves.  In _The 27th International Ocean and Polar Engineering Conference_ , pages 662–669, San Francisco, CA, USA, June 2017.  International Society of Offshore and Polar Engineers (ISOPE), ISBN: 9781880653975\.
  * Higuera et al. [2013] Pablo Higuera, Javier L Lara, and Inigo J Losada.  Realistic wave generation and active wave absorption for navier–stokes models: Application to openfoam®.  _Coastal Engineering_ , 71:102–118, January 2013.  Elsevier, doi:10.1016/j.coastaleng.2012.07.002.
  * Torres-Freyermuth et al. [2010] Alec Torres-Freyermuth, Javier L Lara, and Inigo J Losada.  Numerical modelling of short-and long-wave transformation on a barred beach.  _Coastal Engineering_ , 57(3):317–330, 2010.  Elsevier, doi:10.1016/j.coastaleng.2009.10.013.
  * Higuera [2020] Pablo Higuera.  Enhancing active wave absorption in rans models.  _Applied Ocean Research_ , 94:102000, 2020.  ISSN:0141-1187, doi:10.1016/j.apor.2019.102000.
  * [20] Pál Schmitt and Björn Elsaesser.  A review of wave makers for 3d numerical simulations.
  * Miquel et al. [2018] Adria Miquel, Arun Kamath, Mayilvahanan Alagan Chella, Renata Archetti, and Hans Bihs.  Analysis of different methods for wave generation and absorption in a cfd-based numerical wave tank.  _Journal of Marine Science and Engineering_ , 6(2):73, June 2018.  Multidisciplinary Digital Publishing Institute, doi:10.3390/jmse6020073.
  * Dean and Dalrymple [1991] RG Dean and RA Dalrymple.  _Water wave mechanics for engineers and scientists_ , volume 2.  Word Scientific, 1991.  ISBN: 9810204205.
  * Gamaleldin and Babanin [2018] M.W. Gamaleldin and A.V. Babanin.  On the numerical simulation of extreme non-breaking waves interaction with piled structures.  In _Australia-China Centre for Maritime Engineering (ACCME) $3^{rd}$ Workshop on Ocean Surface Waves and Wave-Coupled Processes_, Wuyang International hotel, Hangzhou, China, Oct 2018.
  * Gamaleldin and Babanin [2019] M.W. Gamaleldin and A.V. Babanin.  On phase-resolving simulations of deep-water waves: Wave absorption using a static-boundary method.  In _2nd International Workshop on Waves, Storm Surges, and Coastal Hazards 16th International Workshop on Wave Hindcasting and Forecasting_ , Sheraton Hotel, Melbourne, Australia, 10-15 November 2019.  P26.
  * Gamaleldin and Babanin [2020] M.W. Gamaleldin and A.V. Babanin.  On wave-absorbing boundaries for numerical wave flumes: Geometrical optimization of the static-boundary method for deep-water conditions.  In _The 15th OpenFOAM Workshop_ , Arlington, VA, USA, 22-26 June 2020\.  Technical Session II-A.
  * Ursell et al. [1960] Fritz Ursell, Robert George Dean, and YS Yu.  Forced small-amplitude water waves: a comparison of theory and experiment.  _Journal of Fluid Mechanics_ , 7(1):33–52, January 1960.  Cambridge University Press, doi:10.1017/S0022112060000037.
  * Galvin Jr [1964] Cyril J Galvin Jr.  Wave-height prediction for wave generators in shallow water.  Technical report, March 1964.  pages 21, COASTAL ENGINEERING RESEARCH CENTER VICKSBURG MS, Accession Number: AD0440880.
  * Havelock [1929] TH Havelock.  Lix. forced surface-waves on water.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 8(51):569–576, April 1929.  Taylor & Francis, doi:10.1080/14786441008564913.
  * Frigaard and Christensen [1994] Peter Frigaard and Morten Christensen.  An absorbing wave-maker based on digital filters.  In _24th International Conference on Coastal Engineering_ , volume 1, pages 168–168, Kobe, Japan, October 1994.  doi:10.1061/9780784400890.014.
  * Andreas [2006] Antonio Andreas.  _Digital signal processing: Signals, systems, and filters_.  McGraw-Hill, New York, 2006.  ISBN: 0071454241, doi:10.1036/0071454241.
  * Wellens and Borsboom [2020] Peter Wellens and Mart Borsboom.  A generating and absorbing boundary condition for dispersive waves in detailed simulations of free-surface flow interaction with marine structures.  _Computers & Fluids_, 200:104387, 2020.  Elsevier, doi:10.1016/j.compfluid.2019.104387.
  * Weller et al. [1998] Henry G Weller, Gavin Tabor, Hrvoje Jasak, and Christer Fureby.  A tensorial approach to computational continuum mechanics using object-oriented techniques.  _Computers in physics_ , 12(6):620–631, 1998\.  AIP, doi:10.1063/1.168744.
  * Ope [2018] _OpenFOAM The Open Source CFD Toolbox User Guide_.  OpenCFD Limited, v1812 edition, 12 2018.  https://www.openfoam.com/documentation/user-guide/.
  * Morgan et al. [2011] Gerald Charles John Morgan, Jun Zang, Deborah Greaves, Andrew Heath, Chris Whitlow, and John Young.  Using the rasinterfoam cfd model for wave transformation and coastal modelling.  _Coastal Engineering Proceedings_ , 1(32):23, 2011\.
  * Morgan and Zang [2011] Gerald Morgan and Jun Zang.  Application of openfoam to coastal and offshore modelling.  _The 26th IWWWFB. Athens, Greece_ , 2011.
  * Folley et al. [2012] Matt Folley, Aurélien Babarit, Ben Child, David Forehand, Louise O’Boyle, Katie Silverthorne, Johannes Spinneken, Vasiliki Stratigaki, and Peter Troch.  A review of numerical modelling of wave energy converter arrays.  In _ASME 2012 31st International Conference on Ocean, Offshore and Arctic Engineering_ , pages 535–545, Rio de Janeiro, Brazil, July 2012.  American Society of Mechanical Engineers (ASME), OMAE2012-83807, doi:10.1115/OMAE2012-83807.
  * Rusche [2003] Henrik Rusche.  _Computational fluid dynamics of dispersed two-phase flows at high phase fractions_.  PhD thesis, Imperial College London (University of London), 2003.  1–343, ISNI 0000 0001 3539 7210.
  * Berberović et al. [2009] Edin Berberović, Nils P van Hinsberg, Suad Jakirlić, Ilia V Roisman, and Cameron Tropea.  Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution.  _Physical Review E_ , 79(3):036306(1–15), March 2009.  doi:10.1103/PhysRevE.79.036306.
  * Deshpande et al. [2012] Suraj S Deshpande, Lakshman Anumolu, and Mario F Trujillo.  Evaluating the performance of the two-phase flow solver interfoam.  _Computational science & discovery_, 5(1):1–36, November 2012.  IOP Publishing, doi:10.1088/1749-4699/5/1/014016.
  * Versteeg and Malalasekera [2007] Henk Kaarle Versteeg and Weeratunge Malalasekera.  _An introduction to computational fluid dynamics: the finite volume method_.  Pearson education, 2007.  ISBN: 9780131274983.
  * Harten [1984] Ami Harten.  On a class of high resolution total-variation-stable finite-difference schemes.  _SIAM Journal on Numerical Analysis_ , 21(1):1–23, July 1984.  SIAM, doi:10.1137/0721001.
  * Van Leer [1974] Bram Van Leer.  Towards the ultimate conservative difference scheme. ii. monotonicity and conservation combined in a second-order scheme.  _Journal of computational physics_ , 14(4):361–370, March 1974.  Elsevier, doi:10.1016/0021-9991(74)90019-9.
  * [43] Mesh generation with the blockmesh utility.  https://cfd.direct/openfoam/user-guide/v6-blockmesh/.  CFD Direct, Accessed: 2019-10-03.
  * Lafayette et al. [2016] Lev Lafayette, Greg Sauter, Linh Vu, and Bernard Meade.  Spartan performance and flexibility: An hpc-cloud chimera.  Barcelona, OpenStack Summit, October 27 2016.  doi: org/10.4225/49/58ead90dceaaa.
  * Mansard and Funke [1980] Etienne PD Mansard and ER Funke.  The measurement of incident and reflected spectra using a least squares method.  _Coastal Engineering Proceedings_ , 1(17):154–172, December 1980.  doi:10.1061/9780872622647.008.
  * Celik et al. [2008] Ishmail B Celik, Urmila Ghia, Patrick J Roache, et al.  Procedure for estimation and reporting of uncertainty due to discretization in cfd applications.  _Journal of fluids Engineering-Transactions of the ASME_ , 130(7), 2008.  American Society of Mechanical Engineers (ASME).
  * Roache [1998] Patrick J Roache.  _Verification and validation in computational science and engineering_.  Hermosa Publishers, ISBN: 0913478083, August 1998.
  * Richardson [1911] Lewis Fry Richardson.  The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam.  _Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character_ , 210:307–357, January 1911.  doi:10.1098/rsta.1911.0009.
  * Richardson and Gaunt [1927] Lewis F Richardson and J Arthur Gaunt.  The deferred approach to the limit. part i. single lattice. part ii. interpenetrating lattices.  _Philosophical Transactions of the Royal Society of London. Series A, containing papers of a mathematical or physical character_ , 226:299–361, 1927.  JSTOR, ISSN 02643952.
  * Israeli and Orszag [1981] Moshe Israeli and Steven A Orszag.  Approximation of radiation boundary conditions.  _Journal of computational physics_ , 41(1):115–135, 1981.  Elsevier, doi:10.1016/0021-9991(81)90082-6.
  * Romate [1992] JE Romate.  Absorbing boundary conditions for free surface waves.  _Journal of computational Physics_ , 99(1):135–145, 1992.  Elsevier, doi:10.1016/0021-9991(92)90281-3.
  * Bredmose et al. [December 2013] Henrik Bredmose, Jesper Mariegaard, Bo Terp Paulsen, Bjarne Jensen, Signe Schløer, Torben J Larsen, Taeseong Kim, and Anders Melchior Hansen.  The wave loads project.  December 2013.  DTU Wind Energy Report E–0045.
  * Paulsen et al. [2014] Bo Terp Paulsen, Henrik Bredmose, and Harry B Bingham.  An efficient domain decomposition strategy for wave loads on surface piercing circular cylinders.  _Coastal Engineering_ , 86:57–76, 2014.  Elsevier, doi:10.1016/j.coastaleng.2014.01.006.
  * Le Méhauté [1976] Bernard Le Méhauté.  _An introduction to hydrodynamics and water waves_.  Springer Science & Business Media, 1976.  ISBN: 9783642855696, doi:10.1007/978-3-642-85567-2.
  * Sumer and Fredsøe [2006] B Mutlu Sumer and Jørgen Fredsøe.  _Hydrodynamics around cylindrical strucures_ , volume 26.  World scientific, 2006.  ISBN: 9812700390.


# Attraction Between Electron Pairs in High Temperature Superconductors

Phil Attard<EMAIL_ADDRESS>
4 Mar., 2022

###### Abstract

It is proposed that in high temperature superconductors Cooper pairs form and
condense due to the monotonic-oscillatory transition in the pair potential of
mean force, which occurs quite generally at high coupling in charge systems.
It is shown that the predicted transition temperatures are broadly in line
with measured superconducting transition temperatures for reasonable values of
the total electron density and the residual dielectric permittivity arising
from the immobile electrons. The predicted transition is independent of the
isotopic masses of the solid. Consequent design principles for high
temperature superconductors are discussed.

##  I Introduction

The Bardeen-Cooper-Schrieffer (BCS)BCS57  theory of superconductivity is based
on electron pairs forming zero-momentum, zero-spin bosons bound together by an
attractive interaction. For normal low temperature superconductors the
accepted mechanism for the attraction is the dynamic interaction of the
electron pair with phonons, which are the quantized vibrations of the solid
lattice. The dependence of the transition temperature on the isotopic masses
of the solidMaxwell50 ; Reynolds50  confirms BCS theory.

In the late 1980’s high temperature superconductors with transition
temperatures above 30 K were discovered. Bednorz86 ; Wu87  These are
independent of the isotopic masses, which rules out the phonon exchange
mechanism. New materials with ever higher superconducting transition
temperatures have been discovered in the ensuing years, but despite many
proposals, Anderson87 ; Bickers87 ; Inui88 ; Gros88 ; Kotliar88 ; Mann11 ;
Monthoux91  no consensus has emerged for the nature of the attractive
potential that forms Cooper pairs in these materials.

In this paper I propose that the electron attraction responsible for Cooper
pairs in high temperature superconductors is due to the oscillatory pair
static correlation function that occurs at high coupling. This has long been
established for like-charged particles in the one component plasma and in
primitive model electrolytes. Brush66 ; Stillinger68 ; Fisher69 ; Outhwaite78
; Stell76 ; Parrinello79 ; Attard93 ; Ennis95  Here the high-temperature
superconductor is modeled as a one-component plasma _in media_ , with the
relatively few electrons in the Fermi foam comprising the fluid charges, and
the fixed nuclei and majority immobile electrons in the Fermi sea forming the
neutralizing background and static relative permittivity (dielectric
constant). The predicted temperatures for the monotonic-oscillatory transition
in this model encompass the measured transition temperatures for high
temperature superconductors for physically reasonable values of the total
electron density and static relative permittivity.

##  II Model and Analysis

The solid conductor is modeled as a one component plasma (mobile electrons in
a uniform counter-charge background), together with a finite relative
permittivity, $\epsilon_{\mathrm{r}}={\cal O}(10^{2})$, that results from the
remaining immobile but polarizable electrons (i.e. those deep in the Fermi
sea). The mobile electrons at the Fermi foam have number density
$\rho_{\mathrm{F}}(T)$ derived below. Arguments concerning this model are
canvassed in the conclusion.

I begin with the restricted primitive model electrolyte for three reasons:
First, there are a wealth of analytic, numeric, and experimental results known
for electrolytes. Second, it shows the generality of the monotonic-oscillatory
transition in charge systems. And third, it gives a specific value for the
width of the accessible energy states, which is required to determine the
electron density of the Fermi foam.

In the restricted primitive model electrolyte (ions of equal hard sphere
diameter), the pair distribution function undergoes an oscillatory transition
whenAttard93

$\kappa_{\mathrm{D}}d\geq\sqrt{2},$ (2.1)

where $d$ is the hard core diameter of the ions. The inverse Debye screening
length for the binary symmetric electrolyte is
$\kappa_{\mathrm{D}}=\sqrt{(4\pi\beta/\epsilon)2\rho_{\mathrm{F}}q^{2}}$,
where $q$ is the ionic charge (in this case the electron charge), and
$\rho_{\mathrm{F}}$ is the number density of each type of ion. Here
$\beta=1/k_{\mathrm{B}}T$ is the inverse temperature, and
$\epsilon=4\pi\epsilon_{0}\epsilon_{\mathrm{r}}$ is the total permittivity of
the medium, $\epsilon_{0}$ being the permittivity of free space (SI units).
This result is based on the Debye-Hückel form for the pair distribution
function combined with the exact Stillinger-Lovett second moment condition.
More accurate analytic and numeric approximations exist,Attard93  but this is
sufficient for the present purposes.

To make the connection with the one component plasma, which does not impose a
hard core diameter, the distance of closest approach of the electrons can be
set as the point at which the Coulomb potential _in media_ reaches several
times the thermal energy, $u(d)=\alpha k_{\mathrm{B}}T$, or $d=\beta
e^{2}/\epsilon\alpha$. With these the oscillatory transition in the symmetric
electrolyte occurs when

$2\leq\frac{4\pi\beta
2\rho_{\mathrm{F}}e^{2}}{\epsilon}\frac{\beta^{2}e^{4}}{\epsilon^{2}\alpha^{2}}=\frac{6}{\alpha^{2}}\Gamma^{3}.$
(2.2)

Here the plasma coupling parameter with finite relative permittivity is
$\Gamma\equiv\beta e^{2}/[\epsilon(3/4\pi\rho_{\mathrm{F}})^{1/3}]$.

As mentioned, the parameter $\alpha$ is the multiple of the thermal energy
which bounds the accessible states. Choosing $\alpha=\surd 24\approx 4.9$, the
transition criterion becomes

$\Gamma\geq 2.$ (2.3)

With this value of $\alpha$, the value of the coupling constant at the
transition given for the restricted primitive model electrolyte agrees with
that found by Monte Carlo simulations of the one component plasma.Brush66

I now estimate the density of the electrons in the Fermi foam modeling them as
a non-interacting ideal gas, also known as the free electron model. The Fermi
momentum and the Fermi energy for ideal fermions arePathria72

$p_{\mathrm{F}}=2\pi\hbar\left(\frac{3\rho}{8\pi}\right)^{1/3},\mbox{ and
}\epsilon_{\mathrm{F}}=\frac{(2\pi\hbar)^{2}}{2m}\left(\frac{3\rho}{8\pi}\right)^{2/3},$
(2.4)

where $\rho=N/V$ is the total electron number density. The thermal wavelength
is $\Lambda=\sqrt{2\pi\beta\hbar^{2}/m}$, and
$\beta\epsilon_{\mathrm{F}}=2\pi(3\rho\Lambda^{3}/8\pi)^{2/3}/2$, which is
much larger than unity.

With momentum state spacing being $\Delta_{p}=2\pi\hbar/L$, Messiah61 ;
Merzbacher70  where the volume is $V=L^{3}$, the number in the Fermi foam is

$\displaystyle N_{\mathrm{F}}$ $\displaystyle=$ $\displaystyle
2\Delta_{p}^{-3}\int_{\epsilon_{\mathrm{F}}-\alpha/\beta}^{\epsilon_{\mathrm{F}}+\alpha/\beta}\mathrm{d}\epsilon\,4\pi
m\sqrt{2m\epsilon}\frac{e^{-\beta(\epsilon-\epsilon_{\mathrm{F}})}}{1+e^{-\beta(\epsilon-\epsilon_{\mathrm{F}})}}$
(2.5) $\displaystyle\approx$ $\displaystyle 4\alpha
V\Lambda^{-3}(3\rho\Lambda^{3}/8\pi)^{1/3}.$

An expansion to leading order for large $\beta\epsilon_{\mathrm{F}}$ has been
made to obtain the final equality. That is

$\rho_{\mathrm{F}}\Lambda^{3}=4\alpha\left(\frac{3\rho\Lambda^{3}}{8\pi}\right)^{1/3}.$
(2.6)

The total excitable electron density, $\rho_{\mathrm{F}}$, is significantly
less than the total electron density, $\rho$. It is proportional to the number
of accessible energy states at the Fermi energy, which is fixed by equating
the results of the restricted primitive model to those of the one component
plasma, $\alpha=\surd 24$

The idea that only the electrons at the Fermi surface contribute to screening
also underlies the Thomas-Fermi model of the electron gas. Kittel76  This idea
is taken a little further here by modeling the remaining immobile electrons as
being polarisable and contributing to the residual dielectric constant.

##  III Numerical Results

The predicted oscillatory-monotonic transition temperature is now explored for
a range of the two free parameters: the total electron density $\rho$ and the
residual static relative permittivity $\epsilon_{\mathrm{r}}$. For these a
guide is provided by values for ceramic materials. The total electron density
of zirconia ZrO2 is $\rho=1.65\times 10^{30}$ m-3. The relative permittivity
of typical ceramic insulators is on the order of $\epsilon_{\mathrm{r}}=$
$10^{1}$–$10^{2}$. ceram

Figure 1:  Coion pair potential of mean force as a function of separation in
the symmetric binary monovalent electrolyte ($d=3.41$ Å,
$\epsilon_{\mathrm{r}}=100$, 1 M concentration, HNC approximation). The solid
curve is for $T=130$ K ($\kappa_{\mathrm{D}}^{2}d^{2}=2.3$, $\Gamma=1.8$), the
short-dash curve is for $T=100$ K ($\kappa_{\mathrm{D}}^{2}d^{2}=2.9$,
$\Gamma=2.3$), the long-dash curve is for $T=70$ K
($\kappa_{\mathrm{D}}^{2}d^{2}=4.2$, $\Gamma=3.3$), and the dot-dash curve is
for $T=50$ K ($\kappa_{\mathrm{D}}^{2}d^{2}=5.9$, $\Gamma=4.5$). The solid
line is an eye guide.

Figure 1 shows the coion pair potential of mean force above and below the
oscillatory transition close to contact. The results were obtained with the
hypernetted chain approximation for the restricted primitive model binary
electrolyte. Attard93  The qualitative difference that the transition makes is
apparent. At higher temperatures than those shown, the potential of mean force
is monotonic repulsive and exponentially decaying. At the highest temperature
shown, which is below the one component plasma transition value but above that
predicted for the restricted primitive model, the potential of mean force
turns negative at about $r=12.5$ Å and has very long oscillatory period
(${\cal O}(10^{3})$ Å). The primary minimum in the clearly oscillatory case
following the transition is rather shallow. It becomes deeper with
approximately constant width as the temperature is decreased (i.e. the
coupling is increased). For the present parameters, the location of the
transition predicted by the Debye-Hückel–Stillinger-Lovett approximation for
the primitive model electrolyte, Eq. (2.1), is more or less equal to that
observed in the one component plasma, Eq. (2.3).

The molecular interpretation of the attraction is that it arises at high
coupling from over-charging by counterions (or the background charge) combined
with packing constraints.Attard93  As mentioned, it also occurs in the one
component plasma.Brush66

Figure 2:  Plasma coupling parameter _in media_ using electron density
$\rho_{\mathrm{F}}(T)$, Eq. (2.6), and, from bottom to top, relative
permittivity $\epsilon_{\mathrm{r}}=$ 200, 150, 100, and 75. The dotted line
marks the transition to oscillatory behavior. In all cases $\alpha=\surd 24$
and $\rho=10^{30}$ m-3.

Figure 2 shows the value of the plasma coupling parameter as a function of the
relative permittivity for fixed total electron density. Values $\Gamma>2$ mark
oscillatory pair correlation functions. One sees that this occurs for low
temperatures. For the lowest relative permittivities shown,
$\epsilon_{\mathrm{r}}=$ 100 and 75, the transition temperatures are greater
than those measured for high temperature superconductors. Bednorz86 ; Wu87

Figure 3 shows the transition temperature as a function of relative
permittivity for several values of the electron density. The transition
temperature increases with decreasing relative permittivity and with
increasing electron density. It is more sensitive to changes in the
permittivity than the electron density. The range of calculated transition
temperatures encompasses those measured for high temperature superconductors.
From the figure one can conclude that a high superconducting transition
temperature requires both a high electron density and a low static relative
permittivity. Since the latter can be expected to increase with increasing
electron density, there is an obvious competition between these two
requirements. This undoubtedly restricts the structure and composition of
suitable candidate materials for high temperature superconductivity.

Figure 3:  Transition temperature as a function of the relative permittivity.
From left to right the curves are for a total electron density of $\rho=$ 0.1,
1, and 10 $\times 10^{30}$ m-3.

##  IV Conclusion

Perhaps the most controversial aspect of the proposed model is the invocation
of a finite relative permittivity for the solid. The perceived wisdom is that
a conductor has infinite permittivity. Also, the one component plasma is
usually modeled _in vacuo_ , $\epsilon_{\mathrm{r}}=1$. In response to these
anticipated objections I make the following points. First, one has to
distinguish between an experimental challenge and a fundamental limitation of
nature. There is no doubt that a macroscopic measurement of the dielectric
constant of a conductor yields an infinite dielectric constant in the zero
frequency limit. In my opinion, all this says is that the conductivity
dominates the measurement; it does not say that the residual relative
permittivity is either unity or infinity. Second, fixed atoms surrounded by
the electrons deep in the Fermi sea, which include the immobile inner shell
electrons, remain polarizable, and therefore they must contribute to a finite
relative permittivity. Third, if the relative permittivity were truly infinite
at the molecular level, then the immobile electrons could not interact via the
Coulomb repulsion, and they would be utterly transparent to each other, which
is obvious nonsense. Fourth, and finally, insulators that are close in
chemical composition and physical structure to specific high temperature
superconductors, have a measured finite relative permittivity greater than
unity. For example, cuprate superconductors are insulators if the doping
fraction is less than 0.1 holes per CuO2.CuO2  It seems plausible that the
small changes in composition that turn these into conductors do not much
change the residual static relative permittivity. Instead the simpler
interpretation is that the infinite conduction permittivity swamps any attempt
to measure macroscopically the finite residual static relative permittivity.

The most reliable way to measure the residual static dielectric permittivity
of a high temperature superconductor may turn out to be by extrapolation from
an insulator of similar chemical and structural composition.

As mentioned in connection with Fig. 3, a high monotonic-oscillatory
transition temperature relies upon the competing requirements of high electron
density and low relative permittivity. One might speculate that the reason for
the prevalence of layered structures amongst high temperature superconductors
is that they combine a high electron density within the conducting planes
together with a low electron density and hence low polarisability in the
interlayer space. At the level of dielectric continuua, the consequent image
charges increase the coupling between the electrons in the dense layerAttard88
and increase the temperature of the oscillatory transition compared to a
uniform dielectric medium.

Of course one has to be a little skeptical about the quantitative accuracy of
modeling a crystalline layered conductor as a homogeneous charge fluid. On the
one hand the universality of the monotonic-oscillatory transition for the
like-charge pair potential of mean force at high coupling is robust and
undeniable. On the other hand, the quantitative prediction of the transition
temperature and the depth of the primary minimum cannot be taken too literally
for electrons in the layered crystalline solids that are of interest in high
temperature superconductivity.

This paper proposes that the oscillatory potential of mean force that occurs
at high coupling is responsible for Cooper pair formation and
superconductivity. The BCS theory is predicated upon an attractive potential,
not an attractive potential of mean force. The difference between the pair
potential and the pair potential of mean force is fundamentally the difference
between quantum mechanics and quantum statistical mechanics. It is the latter,
not the former, that is the appropriate theory for condensed matter. In the
accompanying paper I give a new quantum statistical mechanical theory for
superconductivity that shows explicitly how Cooper pairs form and condense
depending upon the potential of mean force. The theory uses the classical
phase space formulation of quantum statistical mechanics,Attard18 ; Attard21
together with techniques recently developed for superfluidity. Attard21 ;
Attard22

## References

  * (1) J. Bardeen, L. N. Cooper, and J. R. Schrieffer, “Theory of Superconductivity”, _Phys. Rev._ 108, 1175 (1957).
  * (2) E. Maxwell, “Isotope Effect in the Superconductivity of Mercury”, _Phys. Rev._ 78, 477 (1950).
  * (3) C. A. Reynolds, B. Serin, W. H. Wright, and L. B. Nesbitt, “Superconductivity of Isotopes of Mercury”, _Phys. Rev._ 78, 487 (1950).
  * (4) J. G. Bednorz and K. A. Möller, “Possible high $T_{\mathrm{c}}$ superconductivity in the Ba-La-Cu-O system”, _Z. Phys. B_ 64, 189 (1986).
  * (5) M. K. Wu, J. R. Ashburn, C. J. Torng, P. H. Hor, R. L. Meng, L. Gao, Z. J. Huang, Y. Q. Wang, and C. W. Chu, “Superconductivity at 93 K in a New Mixed-Phase Y-Ba-Cu-O Compound System at Ambient Pressure”, _Phys. Rev. Lett._ 58, 908 (1987).
  * (6) P. Anderson, “The resonating valence bond state in La2CuO4 and superconductivity”, _Science_ 235, 1196 (1987).
  * (7) N. E. Bickers, D. J. Scalapino, and R. T. Scalettar, “CDW and SDW mediated pairing interactions”, _Int. J. Mod. Phys. B_ 1, 687 (1987).
  * (8) M. Inui, S. Doniach, P. J. Hirschfeld, A. E. Ruckenstein, Z. Zhao, Q. Yang, Y. Ni, and G. Liu, “Coexistence of antiferromagnetism and superconductivity in a mean-field theory of high-Tc superconductors”, _Phys. Rev. B_ 37, 5182 (1988).
  * (9) C. Gros, D. Poilblanc, T. M. Rice, and F. C. Zhang, “Superconductivity in correlated wavefunctions”, _Physica C_ , 153–155, 543 (1988).
  * (10) G. Kotliar, and J. Liu, “Superexchange mechanism and d-wave superconductivity”, _Phys. Rev. B_ 38, 5142 (1988).
  * (11) A. Mann, “High-temperature superconductivity at 25: Still in suspense”, _Nature_ 475 280 (2011).
  * (12) P. Monthoux, A. V. Balatsky, and D. Pines, “Toward a theory of high-temperature superconductivity in the antiferromagnetically correlated cuprate oxides”, _Phys. Rev. Lett._ 67, 3448 (1991).
  * (13) S. G. Brush, H. L. Sahlin and E. Teller, “Monte Carlo Study of a One-Component Plasma. I. _J. Chem. Phys._ 45, 2102 (1966).
  * (14) F. H. Stillinger and R. Lovett, “Ion-Pair Theory of Concentrated Electrolytes. I. Basic Concepts”, _J. Chem. Phys._ 48, 3858 (1968).
  * (15) M. Parrinello and M. P. Tosi, “Structure and dynamics of simple ionic liquids”, _Rev. Nuovo Cimento_ 2 1 (1979).
  * (16) C. W. Outhwaite, “Modified Poisson-Boltzmann equation in electric double layer theory based on the Bogoliubov-Born-Green-Yvon integral equations”, _J. Chem. Soc. Faraday Trans. 2_ 74, 1214 (1978).
  * (17) M. E. Fisher and B. Widom, “Decay of Correlations in Linear Systems”, _J. Chem. Phys._ 50, 3756 (1969).
  * (18) G. Stell, K. C. Wu, and B. Larsen, “Critical Point in a Fluid of Charged Hard Spheres”, _Phys. Rev. Lett._ 37, 1369 (1976).
  * (19) P. Attard, “Asymptotic Analysis of Primitive Model Electrolytes and the Electrical Double Layer”, _Phys. Rev. E_ 48, 3604 (1993).
  * (20) J. Ennis, R. Kjellander, and D. J. Mitchell, “Dressed Ion Theory for Bulk Symmetric Electrolytes in the Restricted Primitive Model”, _J. Chem. Phys._ 102, 975 (1995).
  * (21) R. K. Pathria, _Statistical Mechanics_ (Pergamon Press, Oxford, 1972).
  * (22) A. Messiah _Quantum Mechanics_ (Vol 1 and 2) (North-Holland, Amsterdam, 1961).
  * (23) E. Merzbacher, _Quantum Mechanics_ 2nd edn (Wiley, New York, 1970).
  * (24) C. Kittel, _Introduction to Solid State Physics_ 5th edn (Wiley, New York, 1976).
  * (25) /www.ceramtec.com/ceramic-materials/dielectric/ (Accessed 1 Feb., 2022).
  * (26) http://hyperphysics.phy-astr.gsu.edu/hbase/Solids
/hitc.html (Accessed 1 Feb., 2022).

  * (27) P. Attard, D. J. Mitchell, and B. W. Ninham, “Beyond Poisson-Boltzmann: Images and Correlations in the Electric Double Layer. II. Symmetric Electrolyte.” J. Chem. Phys. 89, 4358 (1988).
  * (28) P. Attard, “Quantum Statistical Mechanics in Classical Phase Space. Expressions for the Multi-Particle Density, the Average Energy, and the Virial Pressure”, arXiv:1811.00730 (2018).
  * (29) P. Attard, _Quantum Statistical Mechanics in Classical Phase Space_ , (IOP Publishing, Bristol, 2021).
  * (30) P. Attard, “Bose-Einstein Condensation, the Lambda Transition, and Superfluidity for Interacting Bosons”, arXiv:2201.07382 (2022).