# Numerical study of natural oscillations of supported drops with free and
pinned contact lines

Jordan Sakakeeny  Yue Ling<EMAIL_ADDRESS>Department of Mechanical
Engineering, Baylor University, Waco, TX 76798, USA

###### Abstract

The oscillation of droplets supported by solid surfaces is important to a wide
variety of applications such as dropwise condensation. In the present study,
the axisymmetric natural oscillations of a liquid drop supported by a flat
surface is investigated by direct numerical simulation. The liquid-gas
interface is captured using a geometric volume-of-fluid (VOF) method. A
parametric study is carried out by varying the equilibrium contact angle and
the gravitational Bond number (Bo). Both positive and negative gravities are
considered, and thus the results cover both pendant and sessile drops. To
incorporate the effect of contact line mobility, the two asymptotic limits,
namely the pinned contact line (PCL) and free contact line (FCL) conditions,
are considered and their effects on the drop oscillation features are
characterized. The predicted oscillation frequencies for PCL and FCL serve as
the upper and lower bounds for general situations. The drop oscillation is
initiated by increasing the gravity magnitude for a short time. The first mode
due to the drop centroid translation dominates the excited oscillation. The
oscillation frequency scales with the capillary frequency, and the normalized
frequency monotonically decreases with the equilibrium contact angle. For zero
gravity, the computed frequencies for all contact angles agree remarkably well
with the inviscid theory for both the PCL and FCL conditions. The kinetic
energy correction factor is introduced to account for the additional
contribution of the oscillation-induced internal flow to the overall kinetic
energy of the drop. Both the frequency and the kinetic energy correction
factor increase with Bo, decrease with the contact angle, and increase when
the contact line condition changed from FCL to PCL. The variation of
oscillation frequency due to the change of Bo is particularly significant when
the contact angle is large, suggesting that the gravity effect must be
incorporated to accurately predict the oscillation frequency for drops
supported by hydrophobic or superhydrophobic surfaces.

##  I Introduction

The oscillation of droplets supported by solid surfaces is important to a wide
variety of applications, such as drop shedding on condensation surfaces Yao
_et al._ (2017) and water harvesting Dai  _et al._ (2018). When the surface
normal vector is aligned and opposite of the gravity direction, the supported
drops are also referred to as pendant and sessile drops, respectively. Drop
oscillation induced by mechanical surface vibrations or acoustics has been
shown to enhance drop mobility on the supported surface Yao  _et al._ (2017),
or even to cause drops to detach from the surface Boreyko and Chen (2009). Due
to the resonance effect, when the frequencies for the external forcing match
the natural frequencies of the supported drops, the excited oscillation
amplitude will be maximized for a given energy input Boreyko and Chen (2009);
Yao  _et al._ (2017). Therefore, it is advantageous to accurately predict the
natural oscillation frequencies for supported drops on surfaces of different
material properties.

The natural oscillation of a liquid drop, when there is no external forcing,
is a classic fluid mechanics problem and has been extensively studied in the
past. For a free isolated drop, Rayleigh provided the explicit expression for
the oscillation frequency for a given mode $n$ in the inviscid, free-surface,
and small-amplitude limit Rayleigh (1879). The Rayleigh frequency
$\omega_{\text{Ra}}$ scales with the capillary frequency $\omega_{c}$, and the
ratio $\omega/\omega_{c}$ is a function of the mode number $n$,

$\frac{\omega_{\text{Ra}}^{2}}{\omega_{c}^{2}}=(n-1)n(n+2)$ (1)

for $n\geq 2$. The capillary frequency is defined as
$\omega_{c}=\sqrt{\sigma/(\rho_{l}R_{0}^{3})}$, where $\sigma$ and $\rho_{l}$
are the surface tension and the liquid density, respectively, and $R_{0}$ is
the radius of the spherical drop. The effect of the surrounding fluid has been
incorporated by Lamb Lamb (1932) and the Lamb frequencies can be expressed as

$\frac{\omega_{\text{Lamb}}^{2}}{\omega_{c}^{2}}=\frac{(n-1)n(n+1)(n+2)}{(n+1)+n\rho_{g}/\rho_{l}}\,,$
(2)

where $\rho_{g}$ is the density of the surrounding gas. For cases with a small
density ratio $\rho_{g}/\rho_{l}\ll 1$, the difference between the Lamb and
Rayleigh frequencies is small.

The effect of liquid viscosity on drop oscillation is generally characterized
by the Ohnesorge number (Oh). For drops with finite Oh, the oscillation
amplitude will decrease over time due to viscous dissipation. For small-
amplitude oscillations, the decay of oscillation amplitude follows an
exponential function $A(t)\sim\exp(-\beta t)$, where $\beta$ is the damping
rate, which scales with the viscous frequency $\omega_{v}=\nu_{l}/R_{0}^{2}$,
where $\nu_{l}$ is the kinematic viscosity of the drop liquid. The normalized
damping rate $\beta/\omega_{v}$ is also a function of the mode number as shown
by Lamb Lamb (1932),

$\frac{\beta_{\text{Lamb}}}{\omega_{v}}=(n-1)(2n+1)\,.$ (3)

The oscillation frequency generally decreases with Oh, though quite slowly.
The leading order correction to the oscillation frequency is quadratic,
$\omega^{2}=\omega_{\text{Ra}}^{2}-\beta^{2}$, which can be expanded as

$\frac{\omega^{2}}{\omega_{c}^{2}}=(n-1)n(n+2)-(n-1)^{2}(2n+1)^{2}\
\text{Oh}^{2}+O(\text{Oh}^{3})\,,$ (4)

It is observed that, for drops with low Oh, the viscous effect on oscillation
frequency is small unless the mode number $n$ is very large.

The aforementioned studies all assume the drop oscillation amplitude is small.
As a result, the oscillation is linear and a superposition of different
oscillation modes is allowed. When the oscillation amplitude is finite, the
oscillation becomes nonlinear Trinh and Wang (1982); Tsamopoulos and Brown
(1983) and the additional effects such as inter-mode coupling arise Basaran
(1992). Furthermore, when a drop is moving, such as falling under the action
of gravity, the surrounding gas flows can also influence the drop oscillation
Helenbrook and Edwards (2002); Bergeles  _et al._ (2018); Zhang and Ling
(2019). In the present study, the focus is on low Ohnesorge-number drops
(water droplets of millimeter sizes) and small-amplitude oscillations
(oscillation amplitude lower than 10% of drop radius). Furthermore, the drop
liquid density and viscosity are significantly larger than those of the
surrounding gas. As a result, if the drop is not supported by the solid
surface, the oscillation frequency and the damping rate are expected to follow
the Rayleigh frequencies and Lamb’s damping rate.

When the drop is in contact with and supported by a solid surface, additional
complexities arise due to the interaction between the drop and the surface.
First of all, the supported drop exhibits a first-mode ($n=1$) oscillation,
which is associated with the drop centroid translation Strani and Sabetta
(1984); Sakakeeny and Ling (2020). For a free drop, the first mode and the
corresponding centroid motion does not trigger a shape deformation, if the
effect of ambient fluid is ignored. For a supported drop, however, when the
distance between the drop centroid and the surface varies, there must be a
corresponding deformation of the drop surface Sakakeeny and Ling (2020).

Furthermore, the natural oscillations of the supported drop will also be
influenced by the surface material properties, such as the equilibrium contact
angle and contact-line dynamics Strani and Sabetta (1984); Noblin  _et al._
(2004); Chang  _et al._ (2013); Bostwick and Steen (2014). Oscillation of
supported drops can induce motion of the contact line. When contact-line
hysteresis is present, the contact angle varies as the contact line moves. The
angle for an advancing contact line is typically larger than that for the
receding counterpart. The effect of hysteresis is typically characterized by
the difference between the advancing and receding contact angles. Modeling
moving contact lines in continuum mechanics remains an unresolved challenge
Bonn  _et al._ (2009); He and Yao (2019). The present study is focused only on
the two asymptotic limits for the contact-line mobility: the pinned contact
line (PCL) and the free contact line (FCL), see Figs. 1 (a) and (b),
respectively, where the typical singularity behaviors for moving contact lines
Snoeijer and Andreotti (2013) are alleviated. For PCL, it is considered that
the hysteresis effect is strong and the oscillation amplitude is small, so the
contact angle always lies between the receding and advancing contact angles.
As a result, the contact line is fixed/pinned while the contact angle can
vary. For FCL, it is considered that there is no hysteresis effect, so the
contact line can move freely while the contact angle is fixed at its
equilibrium value. The oscillation frequencies for supported drops with
general contact-line conditions will be bounded by these two limits. The
previous inviscid theoretical models Strani and Sabetta (1984); Bostwick and
Steen (2014); Sakakeeny and Ling (2020) all indicated that the oscillation
frequencies decrease with the contact angle for all oscillation modes.
Bostwick and Steen Bostwick and Steen (2014) further indicated that the
oscillation frequency for a given mode number and contact angle increases
significantly if the contact-line condition changes from FCL to PCL.

Finally, the gravity also affects the oscillation of a supported drop. In the
present study, it is considered that the surface is flat and the gravity is
normal to the surface. It is taken that the gravity $g$ is positive when it is
opposite to the surface normal. Therefore, positive and negative values of $g$
represent sessile drops and pendant drops, respectively. The effect of gravity
can be characterized by the gravitational Bond number
$\text{Bo}=\rho_{l}gR_{d}^{2}/\sigma$, where $R_{d}=(3V_{d}/4\pi)^{1/3}$ is
the volume-based radius and $V_{d}$ is the volume of the drop. The Bond number
serves as a measure for the ratio between the gravity and surface tension
contributions. Previous studies of supported drop oscillations often ignore
the effect of Bo Strani and Sabetta (1984); Bostwick and Steen (2014).
Nevertheless, numerical studies have shown that the oscillation frequency
increases with Bo for sessile drops Sakakeeny and Ling (2020) and decreases
with the magnitude of Bo for pendant drops Basaran and DePaoli (1994). With
the present definitions for $g$, Bo varies from negative to positive values,
and the normalized oscillation frequency will then increase monotonically with
Bo, though a more comprehensive investigation is required.

Figure 1: Schematics of supported drop oscillations with (a) free contact line
(FCL) and (b) pinned contact line (PCL), where the solid lines represent the
equilibrium state.

The goal of the present study is to investigate the natural oscillation of a
liquid drop supported by a flat solid surface through direct numerical
simulation. The oscillation is excited by increase the gravity magnitude for a
short time. Since the normal of the surface is taken to be aligned with the
gravity, the induced oscillation is axisymmetric. Particular attention is paid
to the first oscillation mode, because it generally dominates the excited
oscillations. The effects of contact angle, contact line mobility, and
gravitational Bond number on important oscillation features, including the
oscillation frequency, the damping rates, and the oscillation-induced internal
flow, will be characterized through parametric simulations. As an extension of
our former study on oscillation of sessile drops with FCL condition Sakakeeny
and Ling (2020), the present study is focused on the PCL condition, so that a
comprehensive understanding of the effect of contact-line mobility on the
oscillation of supported drops can be established. We will vary the Bond
number for both negative and positive values, to cover both the pendant and
sessile drops regimes.

The rest of the paper will be organized as follows. The simulation approaches,
including the governing equations, the numerical methods, and the simulation
setup, will be presented in Section II. The simulation results will be shown
and discussed in Section III. Finally, we will conclude the key findings in
Section IV.

##  II Simulation methods

###  II.1 Governing equations

The liquid-gas two-phase flow is resolved using the one-fluid approach,
wherein the two phases, liquid and gas, are treated as one fluid with material
properties that change abruptly across the gas-liquid interface. The Navier-
Stokes equations for incompressible flow with surface tension are given as

$\rho(\partial_{t}\mathbf{u}+\mathbf{u}\cdot\nabla\mathbf{u})=-\nabla
p+\nabla\cdot(2\mu\mathbf{D})+\sigma\kappa\delta_{s}\mathbf{n},$ (5)
$\nabla\cdot\mathbf{u}=0,$ (6)

where $\rho$, $\mathbf{u}$, $p$, and $\mu$, represent density, velocity,
pressure, and viscosity, respectively. The strain-rate tensor is denoted by
$\mathbf{D}$. The surface tension term on the right-hand side of Eq. (5) is a
singular term, with the Dirac distribution function $\delta_{s}$ localized on
the interface. The surface tension coefficient is represented by $\sigma$, and
$\kappa$ and $\mathbf{n}$ are the local curvature and unit normal of the
interface, respectively.

The two different phases are distinguished by the liquid volume fraction $C$.
While $C=0$ indicates that the cell is full of gas, $C=1$ indicates that the
cell is full of liquid. For cells with interfaces, $0<C<1$. The temporal
evolution of $C$ satisfies the advection equation,

$\partial_{t}C+\mathbf{u}\cdot\nabla C=0.$ (7)

The fluid density and viscosity are determined by

$\displaystyle\rho$ $\displaystyle=C\rho_{l}+(1-C)\rho_{g}\,,$ (8)
$\displaystyle\mu$ $\displaystyle=C\mu_{l}+(1-C)\mu_{g}\,,$ (9)

where the subscripts $g$ and $l$ correspond to the gas and the liquid phases,
respectively.

###  II.2 Numerical methods

The governing equations (Eqs. (5), (6), and (7)) are solved using the open-
source, multiphase flow solver Basilisk Popinet . The Basilisk solver uses a
finite-volume approach based on a projection method. An adaptive quadtree
spatial discretization is used, which allows for adaptive mesh refinement
(AMR) in user-defined regions. The advection equation (Eq. (7)) is solved via
the piecewise-linear geometrical volume-of-fluid (VOF) method Scardovelli and
Zaleski (1999); Popinet (2009). Compared to other popular interface-capturing
methods, such as the front-tracking Unverdi and Tryggvason (1992), level-set
Sussman  _et al._ (1994), the VOF method has the important advantage of
conserving mass/volume, which is crucial to predicting oscillation frequency,
as desired in the present study, since the frequency is a function of the drop
volume. The surface tension calculation in VOF framework can induce numerical
parasitic currents near the interface Renardy and Renardy (2002). This
numerical issue is solved by combining the balanced-force continuum-surface-
force method for surface tension discretization and the height-function (HF)
method for curvature calculation Popinet (2009). The HF method is additionally
used to specify the contact angle at the surface. The Basilisk solver utilizes
a staggered-in-time discretization of the volume fraction/density and
pressure, leading to a formally second-order-accurate time discretization
Popinet (2009). Numerous validation studies for the numerical methods, as well
as examples of a wide variety of interfacial multiphase flows, can be found on
the Basilisk website and in previous studies e.g., Zhang  _et al._ (2019);
Marcotte and Zaleski (2019); Zhang  _et al._ (2020); Mostert and Deike (2020);
Sakakeeny and Ling (2020).

###  II.3 Physical parameters

In the present study, we consider the axisymmetric natural oscillations of a
viscous liquid drop supported by a flat surface, as shown in Fig. 2. The
physical properties for the liquid and gas phases are taken to be similar to
water and air, respectively. The volume of the drop is kept constant,
$V_{d}=65.45$ µL, across all cases, for which the volume-based radius
$R_{d}=2.5$ mm. The wettability of the surface is characterized by the
equilibrium contact angle, $\theta_{0}$, which is varied from 50 to $None$.
The range of contact angles considered here is sufficient to cover common
hydrophilic, hydrophobic, and superhydrophobic surfaces Yao  _et al._ (2017).
The values of the key physical parameters are listed in Table 1.

While the equilibrium shape for a free drop is a sphere, the equilibrium shape
for a supported drop is a spherical cap, when gravity is absent. The the
radius of the spherical cap is also denoted by $R_{0}$, which varies with
$\theta_{0}$ for a given volume $V_{d}$ as

$\displaystyle R_{0}$
$\displaystyle=\left(\frac{3V_{d}}{\pi(2+cos\theta_{0})(1-cos\theta_{0})^{2}}\right)^{1/3}\,.$
(10)

Since $R_{0}$ better represents the surface curvature, the capillary frequency
is defined based on $R_{0}$ as $\omega_{c}=\sqrt{\sigma/(\rho_{l}R_{0}^{3})}$.

The key dimensionless parameters, defined based on scaling variables $R_{d}$,
$\rho_{l}$, and $\sigma$, are listed in Table 2. It can be seen that the gas-
to-liquid ratios for density and viscosity are quite small, thus the effect of
the surrounding gas on the liquid drop is minimal. The Ohnesorge number
$\text{Oh}=0.0024$ indicates that the effect of viscosity is weak.
Furthermore, variation of Oh due to moderate change of drop volume will have
little effect on the normalized oscillation frequency, $\omega/\omega_{c}$,
where $\omega$ and $\omega_{c}$ are the drop oscillation and capillary
frequencies, respectively. For this reason, we have considered only one drop
volume.

The effect of gravity is characterized by the Bond number, Bo. The value of Bo
can be varied by changing $V_{d}$ or $g$. In the present study, we keep
$V_{d}$ fixed and vary $g$ from -0.98 to 9.8 m/s2. Negative Bo and $g$
represent the cases for pendant drops. The resulting range of Bo is -0.088 to
0.88. For pendant drops with large $\theta_{0}$, the drop can be unstable and
detach from the surface if $|\text{Bo}|$ is large. Therefore, a smaller range
of Bo is considered for the pendant drop than the sessile drop. It is
confirmed that for the range of Bo considered, the equilibrium state of the
supported drop is stable for all $\theta_{0}$ considered. We have also
considered only small-amplitude oscillations, so that the drop will not detach
from the wall.

While Oh and Bo are defined based on $R_{d}$, the Bond and Ohnesorge numbers
can be alternatively defined based on $R_{0}$ as
$\text{Bo}_{0}=\rho_{l}gR_{0}^{2}/\sigma$ and
$\text{Oh}_{0}=\mu_{l}/(\rho_{l}\sigma R_{0})$, which will then vary with
$\theta_{0}$.

| $\rho_{l}$
---
$(kg/m^{3})$
| $\rho_{g}$
---
$(kg/m^{3})$
| $\mu_{l}$
---
$(Pa\cdot s)$
| $\mu_{g}$
---
$(Pa\cdot s)$
| $\sigma$
---
$(N/m)$
| $V_{d}$
---
$(\mu L)$
| $\theta_{0}$
---
(°)
| $g$
---
(m/s2)
$1000$ | $1.2$ | $1\times 10^{-3}$ | $1\times 10^{-5}$ | $0.07$ | $65.45$ | 50 to 150 | -0.98 to 9.8
Table 1: Physical parameters. r | $m$ | Oh | $\theta_{0}$ | Bo
---|---|---|---|---
$\rho_{g}/\rho_{l}$ | $\mu_{g}/\mu_{l}$ | $\mu_{l}/\sqrt{\rho_{l}\sigma R_{d}}$ | (°) | $\rho_{l}gR_{d}^{2}/\sigma$
$0.0012$ | $0.01$ | 0.0024 | 50 to 150 | -0.088 to 0.88
Table 2: Key dimensionless parameters. Figure 2: Simulation setup for (a) free
contact line (FCL) and (b) pinned contact line (PCL) cases.

###  II.4 Simulation setup

####  II.4.1 Computational domain and boundary conditions

The computational domain is the same for all cases, see in Fig. 2. The length
of the square domain edge is $H=4R_{d}$. The axisymmetric boundary condition
is applied on the left surface, while the top and right surfaces are slip
walls. The drop is in contact with the bottom surface. For the FCL condition,
the contact angle on the bottom surface is fixed as the equilibrium contact
angle, namely $\theta=\theta_{0}$, see Fig. 2(a), and the contact line can
move freely when the drop oscillates. The interface normal for a given contact
angle is specified using the height function method Afkhami and Bussmann
(2009). To alleviate the singular behavior at the contact line, such as the
diverging viscous stress Snoeijer and Andreotti (2013), the bottom surface for
FCL cases is taken to be a slip wall.

####  II.4.2 Pinned-contact-line boundary conditions

For the PCL condition, the contact line is pinned at its equilibrium position
$x_{cl,0}$ for a given Bo, and the contact angle can vary freely when the drop
oscillates. To be consistent with the contact line condition, we treat the
bottom surface as a no-slip wall, see Fig. 2(b). To pin the contact line at a
given location $x_{cl,0}$, the contact angle on the bottom surface is
specified as $\theta=\theta_{\min}$ for $x<x_{cl,0}$ and
$\theta=\theta_{\max}$ $x>x_{cl,0}$, see Fig. 2(b). The contact angle needs to
reach $\theta_{\min}$ and $\theta_{\max}$ for the contact line to move to the
the left and right, respectively. If $\theta$ varies between $\theta_{\min}$
and $\theta_{\max}$ when the drop oscillates, then the contact line will not
move. For the present study, the equilibrium contact angle $\theta_{0}$ is
varied from $None$ to $None$. For small-amplitude oscillations, the contact
angle $\theta$ only varies in a small extent around $\theta_{0}$. As long as
$\theta_{\min}$ and $\theta_{\max}$ are sufficiently small and large, the
specific values are immaterial. In the present study, we have set
$\theta_{\min}=$$$ and $\theta_{\max}=$$$, which are shown to be sufficient to
statisfy $$$\ll\theta\ll$$$ for all cases considered and to pin the contact
line effectively when the drop oscillates. Though the abrupt change of contact
angle on the bottom surface can effectively pin the contact line, the boundary
condition introduces small-amplitude velocity fluctuations in the cell where
the contact line is located, which leads to non-physical kinetic energy. To
eliminate this numerical artifact, the fluid velocity in the cells that are
$l_{0}$ away from the contact line location is manually set to be zero in
every time step. For all the simulations, $l_{0}/\Delta_{\min}=4$, which have
been verified to effectively eliminate the numerical velocity oscillation
without influencing the oscillation dynamics.

####  II.4.3 Initial conditions

For a given combination of $\theta_{0}$ and Bo, the initial shape of the
supported drop is taken to be the equilibrium shape. The geometry of the
equilibrium supported drop and the contact line location $x_{cl,0}$ can be
obtained by the equilibrium drop theory and numerically solving a system of
ODE. The details can be found in our previous study and thus are not repeated
here Sakakeeny and Ling (2020). To initiate the shape oscillation, the gravity
magnitude is increased for a short duration. The magnitude of gravity $|g|$ is
increased by $g_{pert}$ for $t\leq t_{pert}$. The perturbation gravity
$g_{pert}=4.9$ m/s2 and the perturbation time
$t_{pert}=0.14\sqrt{\rho_{l}R_{d}^{3}/\sigma}$ for all cases. Due to the
change of gravity, the drop will be pushed down ($g>0$) or pulled up ($g<0$)
and deviate from the equilibrium shape. Once the gravity returns to the
original value for $t>t_{pert}$, the drop deforms toward the equilibrium shape
and starts to oscillate. Since the surface normal is aligned with the gravity,
only the axisymmetric zonal modes will be excited. Though all oscillation
modes will be excited to some extent by this method, the first mode ($n=1$)
dominates other high-order modes.

####  II.4.4 Mesh resolution

A quadtree adaptive mesh is used to discretize the domain. The local cell size
is adapted based on the estimated discretization errors of the liquid volume
fraction and the velocity components. The assessment of discretization error
for each variable is made through a wavelet transform van Hooft  _et al._
(2018). If the estimated error is larger than the specified threshold, the
mesh will be locally refined, or vice versa. For the present simulation, the
normalized error thresholds for the volume fraction and the velocity are set
as 0.001 and 0.0001, respectively. Tests have been made to verify these
thresholds are sufficiently small. A representative snapshot of the mesh close
to to the contact line is shown in Fig. 2(c). The minimum cell size in the
quadtree mesh is controlled by the maximum refinement level, $L$, i.e.,
$\Delta x_{\min}=H/2^{L}$. The mesh for $L=11$ is used in the present
simulation, which corresponds to $R_{0}/\Delta x_{\min}\approx 512$, namely
512 cells across the drop radius.

###  II.5 Summary of simulation cases

To systematically investigate the effects of the equilibrium contact angle
($\theta_{0}$), the Bond number (Bo), and the contact line mobility on the
oscillation of a sessile drop, 11 different values of $\theta_{0}$ (from
$None$ to $None$ with an increment of $None$) and 9 different values of Bo
(from -0.088 to 0.88) have been used for both the FCL and PCL conditions.
Therefore, a total of 198 cases are simulated in the parametric study.

The simulations were performed on the Baylor University cluster _Kodiak_ using
4 to 18 CPU cores (Intel E5-2695 V4). Each simulation case takes about 133 to
195 hours of CPU time to reach the time $t\omega_{c}\approx 105$ (51 s). The
simulation time has been verified to be sufficiently long to measure the
frequencies for the first oscillation mode.

##  III Results

###  III.1 Grid refinement and validation

A grid refinement study varying $L=9$ to 12 has been performed for
$\theta_{0}=$$$ and $\text{Bo}=0$. The results for the temporal evolution of
the drop centroid height $z_{c}$ and the corresponding frequency spectra
obtained by the Fourier transform of the temporal signals are shown in Fig. 3.
The difference between the results for $L=11$ and 10 are almost invisible,
demonstrating that the refinement level $L=11$ is sufficient to fully resolve
the oscillations of supported drops.

Figure 3: (a) Temporal evolution and (b) frequency spectrum of the drop
centroid height $z_{c}$ for $\theta_{0}=$$$ and $\text{Bo}=0$, for three
different mesh refinement levels $L=$9, 10, and 11.

Validation for the present simulation setup for the FCL condition can be found
in our previous study Sakakeeny and Ling (2020). For additional validation of
the PCL condition, we examine whether the contact line is effectively pinned
when the drop oscillates. Representative drop surfaces corresponding to the
maximum, equilibrium, and minimum centroid heights in one first-mode
oscillation cycle for $\theta_{0}=$$,$$,$$$ are shown in Figs. 4(a)-(c),
respectively. The temporal evolution of the contact line $x$-location $x_{cl}$
for $\theta_{0}=$$,$$,$$$ is shown in Fig. 4(d), and it is clearly shown that
the contact line is successfully pinned for all cases shown. The equilibrium
contact line locations for $\theta_{0}=$$$ and $None$ are the same.

Figure 4: Validation study for the PCL condition. (a)-(c) Representative drop
surfaces corresponding to the maximum, equilibrium, and minimum drop heights
in one first-mode oscillation cycle for $\theta_{0}=$$$, $None$, and $None$.
(d) Temporal evolution of the contact line location.

###  III.2 Oscillation frequency

The frequency of the shape oscillation depends on the mode number. In the
present study we focus on the dominant first mode $n=1$. The frequency of the
$n=1$ mode can be measured through the temporal evolution of the centroid
height $z_{c}$. The temporal evolution of $z_{c}$ for $\theta_{0}=$$$ and
$None$ and the PCL condition are shown in Fig. 5(a). The drop for
$\theta_{0}=$$$ and FCL is also shown for comparison. For all three cases
shown here, $\text{Bo}=0$. It can be seen that $z_{c}$ oscillates with respect
to the equilibrium value $z_{c,0}$. The oscillation amplitude is generally
small compared to $z_{c,0}$ so the oscillation is expected to be linear.

Fourier transform is performed to generate frequency spectra, which are used
to identify oscillation frequencies (shown as peaks in the spectra). Higher-
order modes $n>1$ are also observed in the spectra, though the first mode is
clearly the dominant one. As addressed in previous studies (Bostwick and
Steen, 2014; Sakakeeny and Ling, 2020), the flat surface for $\theta_{0}=$$$
and FCL is identical to the symmetric boundary condition. Therefore, a
supported drop with $\theta_{0}=$$$ is equivalent to the top half of a free
drop with twice the size. The oscillation frequency for the $n^{\text{th}}$
mode for a supported drop with $\theta_{0}=$$$ and FCL is identical to that
for the $(2n)^{\text{th}}$ mode for a free drop. Since the oscillation
frequencies for the free drop for the present fluid properties are well
predicted by the Rayleigh frequencies, the values of $\omega_{\text{Ra}}$ for
the $n=2$, 4, and 6 are plotted in Fig. 5(b) for comparison. It is clearly
shown that $\omega_{1}$, $\omega_{2}$, and $\omega_{3}$ for the supported drop
with $None$ and FCL agree very well with $\omega_{2,\text{Ra}}$,
$\omega_{4,\text{Ra}}$, and $\omega_{6,\text{Ra}}$.

The oscillation frequency depends on both the contact angle and the contact
line mobility. It is shown that $\omega/\omega_{c}$ decreases from about 4.5
to 2.1 when $\theta_{0}$ increases from $None$ to $None$. For the same contact
angle, $\theta_{0}=$$$, $\omega/\omega_{c}$ decreases from 4.5 to 2.8 when the
contact line mobility changes from PCL to FCL. For the same initial shape and
perturbation method, the oscillation amplitude of $z_{c}$ for FCL is
significantly larger than that for PCL, since the constraint of the latter
condition on the drop is stronger. Correspondingly, the amplitude for FCL in
the spectrum is also higher than that for PCL.

Figure 5: (a) Temporal evolution and (b) frequency spectrum of the drop
centroid height for different $\theta_{0}$. (c) Variations of normalized
first-mode oscillation frequency $\omega/\omega_{c}$ and the difference
between the two $(\omega_{P}-\omega_{F})/\omega_{c}$ as a function of contact
angle $\theta_{0}$ for $\text{Bo}=0$. In (b) the Rayleigh frequencies (Eq.
(1)) for the $n=2$, 4, and 6 modes are shown for comparison. In (c) the
simulation results (symbols) are compared with the inviscid theory of Bostwick
and Steen Bostwick and Steen (2014) (solid lines). The dashed lines represent
fitting correlations for the simulation results.

The first-mode frequency for a supported drop with both FCL and PCL for
$\text{Bo}=0$ are shown as a function of the equilibrium contact angle
$\theta_{0}$ in Fig. 5(c). For both FCL and PCL conditions,
$\omega/\omega_{c}$ monotonically decreases with $\theta_{0}$. The decrease is
more profound for small $\theta_{0}$ (hydrophilic surfaces), and is more
gradual for large $\theta_{0}$ (hydrophobic or superhydrophobic surfaces).
When $\theta_{0}\to$$$, the supported drop approaches a free drop and the
constraint from the surface disappears. As a result, the oscillation frequency
will reduce to zero. When $\theta_{0}$ is close to $None$, the decrease of
$\omega$ over $\theta_{0}$ becomes very rapid Strani and Sabetta (1984). Yet,
a detailed analysis of the asymptotic behavior of the frequency near the limit
of $\theta_{0}\to$$$ is out of the scope of the present paper.

For all $\theta_{0}$, $\omega/\omega_{c}$ for PCL is higher than that for FCL.
The difference between the two, $(\omega_{P}-\omega_{F})/\omega_{c}$, also
decreases with $\theta_{0}$, where $\omega_{P}$ and $\omega_{F}$ represent the
oscillation frequencies for the PCL and FCL, respectively. As $\theta_{0}$
increases, the contact area decreases. As a result, the constraint of the
surface on the drop reduces, and the effect of contact line mobility
conditions will also become less important. Since PCL and FCL represent the
two asymptotic limiting conditions for the contact line mobility, the
predicted frequencies for PCL and FCL shown here represent the upper and lower
bounds for the first-mode oscillation frequencies for general situations. The
results are useful to estimate the natural frequency of a supported drop on
arbitrary material surfaces.

Bostwick and Steen Bostwick and Steen (2014) have established an inviscid
theoretical model to predict the oscillation frequency for supported drop at
$\text{Bo}=0$. Since Oh in the present case is small, the inviscid theory of
Bostwick and Steen is expected to be a good approximation. Their theoretical
predictions are available for $$$<\theta_{0}<$$$ and are plotted in Fig. 5(c)
for comparison. The agreement between the simulation and theoretical model is
excellent for both the FCL and PCL conditions. The good agreement observed
further validates the simulation results.

For convenience of using the present results, correlations for first-mode
oscillation frequencies for the PCL and FCL conditions as a function of
$\theta_{0}$ are fitted in the following form:

$\displaystyle\log(\omega/\omega_{c})=c_{0}+c_{1}(1+\cos\theta_{0})+\left[\exp\left(\frac{(1+\cos\theta_{0})^{c_{2}}}{c_{3}}\right)-1\right]\,.$
(11)

The fitted constants are
$[c_{0},c_{1},c_{2},c_{3}]=[-0.0901,\,1.15,\,20.3,\,2.16\times 10^{5}]$ for
FCL and $[0.287,\,1.22,\,20.7,\,2.81\times 10^{5}]$ for PCL. The fitting
correlations are plotted in Fig. 5(c) and are shown to well represent the
simulation results. For large $\theta_{0}$, $1+\cos(\theta_{0})$ is small, and
the expression above reduces to a linear function, e.g.,
$\log(\omega/\omega_{c})\approx c_{0}+c_{1}(1+\cos\theta_{0})$. The linear
relation between $\log(\omega/\omega_{c})$ holds for all hydrophobic cases
$\theta>$$$. The correction term
$[\exp\left({(1+\cos\theta_{0})^{c_{2}}}/{c_{3}}\right)-1]$ is mainly used to
account for the deviation of the hydrophilic cases from the linear function.
It is also worth mentioning that the correlation Eq. (11) is strictly valid
for the range of $\theta_{0}$ studied. It is not intended to capture the
asymptotic behavior at $\theta_{0}=$$$, i.e., where the frequency drops
rapidly over $\theta_{0}$ near singularity location $\theta_{0}=$$$ Strani and
Sabetta (1984).

###  III.3 Kinetic energy correction factor

The kinetic energy of the liquid drop can be expressed as

$E_{k}=\rho_{l}\int_{V_{d}}\frac{|\boldsymbol{u}|^{2}}{2}dV\,,$ (12)

which will vary over time as the drop oscillates. The temporal variation of
$E_{k}$ is due to two contributions. The first contribution is related to the
bulk motion of the drop following the velocity of the drop centroid,

$E_{kc}=\frac{m_{d}|\boldsymbol{u_{c}}|^{2}}{2}=\frac{m_{d}w_{c}^{2}}{2}\,.$
(13)

where $w_{c}$ is the $z$-component of centroid velocity, $\boldsymbol{u_{c}}$.
For a free drop, $E_{k}=E_{kc}$ since the translation of the drop does not
induce shape deformation. However, $E_{k}>E_{kc}$ for a supported drop due to
the additional contribution of the internal flow induced by the shape
oscillation.

Figure 6: Temporal evolutions of (a) the flow field in the drop reference
frame, (b) the centroid height $z_{c}$ and velocity $u_{c}$ and (c) the
droplet kinetic energy $E_{k}$ for $\text{Bo}=0$ and PCL condition. The
snapshots of the flow fields in (a) are for $\theta_{0}=$$$ and $None$ and
four critical phases of an oscillation cycle, as indicated in (b) and (c).

The velocity fields around the drop for different contact angles are shown in
Fig. 6(a). The snapshots shown in Fig. 6(a) correspond to the valley, peak,
and two equilibrium positions of the centroid in an oscillation cycle of the
dominant $n=1$ mode, which are also indicated in the time evolutions of
$z_{c}$ and $E_{k}$ in Figs. 6(b) and (c). The velocity here is in the drop
reference frame and thus the contribution of the bulk motion has been
subtracted. When the drop moves upward, see column (ii), the typical straining
flow pattern can be recognized inside the drop. Furthermore, a clockwise
circulation is generated on the left top corner of the drop. When the drop
moves downward, see column (iv), the directions of the circulation and
internal straining flow reverse. The flow pattern for $\theta=$$$ is quite
similar to that for a free drop undergoing a $n=2$ mode oscillation. This
similarity in the drop shape for supported drops with large $\theta_{0}$ has
been observed by Strani and Sabetta Strani and Sabetta (1984). When
$\theta_{0}$ decreases, such as $\theta=$$$, the flow pattern will become less
similar to the free drop $n=2$ mode.

When the drop centroid passes the equilibrium positions (ii) and (iv), the
centroid velocity reaches local maximum, and in the mean time the internal
flow is also intense. When the drop centroid reaches the local minimum (i) and
maximum (iii), $u_{c}$ becomes zero and the flow around the drop is minimal.
This indicates that the two contributions to $E_{k}$, i.e., the one from the
centroid motion ($E_{kc}$) and the one due to the shape-oscillation-induced
internal flow ($E_{k}-E_{kc}$) are in phase. Correspondingly, the temporal
evolutions of $E_{k}$ and $E_{kc}$ are also in phase. This conclusion is
further confirmed in Fig. 7(a)), where $E_{k}$ is plotted as a function of
$E_{kc}$ for $\theta_{0}=$$$ and both FCL and PCL conditions. It is observed
that, for both cases, $E_{k}$ varies approximately linearly with $E_{kc}$.
This interesting feature allows us to make an approximation of $E_{k}$ as

$E_{k}\approx\zeta E_{k,c}\,,$ (14)

where $\zeta$ is the kinetic energy correction factor, which is time
independent. The value of $\zeta$ can be obtained by fitting the simulation
results of $E_{k}$ vs $E_{kc}$ (see Fig. 7(a)). It can been seen that the PCL
case exhibits a steeper slope than the FCL case, thus the $\zeta$ value is
greater. The temporal evolutions of $E_{k}$ and the approximation $\zeta
E_{kc}$ are plotted in Fig. 7(b), which affirms that $\zeta E_{kc}$ agrees
well with $E_{k}$ for all time for both FCL and PCL condtions. Here we only
show the results for $\theta_{0}=$$$ as an example. The approximation Eq. (14)
is valid for all $\theta_{0}$. The kinetic energy correction factor and the
approximation Eq. (14) are useful to develop theoretical model to predict the
oscillation frequency as shown in previous study Sakakeeny and Ling (2020).

Figure 7: (a) Variation of $E_{k}$ as a function of $E_{kc}$ and (b) the
temporal evolutions for $E_{k}$ and $\zeta E_{kc}$ for $\theta_{0}=$$$,
$\text{Bo}=0$, and both FCL and PCL conditions. (c) Variations of $\zeta$ for
PCL and FCL conditions and the difference between the two as a function of
contact angle $\theta_{0}$ for $\text{Bo}=0$. The dashed lines in (c) are
fitting correlations.

The variation of $\zeta$ as a function of $\theta_{0}$ for $\text{Bo}=0$ is
shown in Fig. 7(c). Similar to $\omega/\omega_{c}$, $\zeta$ also decrease with
$\theta_{0}$. For all $\theta_{0}$, $\zeta$ for PCL is larger that the FCL
counterpart. This is again due to the stronger constraint from the surface for
PCL. The difference between the values of $\zeta$ for PCL and FCL conditions,
i.e., $\zeta_{P}-\zeta_{F}$, also decreases with $\theta_{0}$. As $\theta_{0}$
approaches $None$, both $\zeta_{P}$ and $\zeta_{F}$ approaches one since there
is neither shape deformation for the first mode nor the additional kinetic
energy contribution from the oscillation-induced flow. As a result,
$\zeta_{P}-\zeta_{F}$ will reach zero.

Similar to the oscillation frequency, correlations are also made for $\zeta$
for the PCL and FCL conditions as a function of $\theta_{0}$ in a similar
form:

$\displaystyle\log(\zeta(\theta_{0}))=e_{0}+e_{1}(1+\cos\theta_{0})+\left[\exp\left(\frac{(1+\cos\theta_{0})^{e_{2}}}{e_{3}}\right)-1\right]\,.$
(15)

The fitted constants are $[e_{0},e_{1},e_{2},e_{3}]=[0,\,0.753,\,5.87,\,54.6]$
for FCL and $[0.14,\,1.00,\,9.54,\,599]$ for PCL. The fitting correlations are
plotted in Fig. 7(c) are found to well represent the simulation results.

###  III.4 Viscous damping of oscillation

Due to the viscous effect, fluid motion induced by shape oscillation will
dissipate the energy provided by the initial excitation. As a result, the
oscillation amplitude will decay over time. In the linear regime, the
oscillation amplitude $A$ follows the exponential function in time,

$A(t)=A_{0}e^{-\beta t},$ (16)

where $A_{0}$ is the initial amplitude. For a free drop, the damping rate
normalized by the viscous frequency, i.e., $\beta/\omega_{v}$ is a function of
the mode number, as indicated in Eq. (3). For a supported drop, the damping
rate will also be influenced by the contact angle and the contact line
mobility. The decay of the oscillation amplitude for the present problem
mainly reflects the damping rate of the dominant $n=1$ mode.

The temporal evolution of $|z_{c}-z_{c,0}|/z_{c}$, for $\theta_{0}=$$$, is
plotted in Fig. 8(a) for both FCL and PCL. As discussed above, the $n=1$ mode
of the supported drop with FCL and $\theta_{0}=$$$ is similar to the $n=2$
mode of the free drop with twice the size. Therefore, the damping rate $\beta$
is expected to be the same as the $\beta_{\text{Lamb}}$ for $n=2$, as given in
Eq. (3). In Fig. 8(a), it is can be observed that the oscillating amplitude
decay for FCL and $\theta_{0}=$$$ agrees very well with the Lamb’s prediction.
When the contact line is pinned, the damping rate increases slightly. The
oscillation amplitude damping is due to viscous dissipation of kinetic energy.
For FCL, since the slip boundary condition is invoked on the bottom surface,
dissipation is only caused by viscous fluid motion inside the drop. For PCL,
additional dissipation is induced by the no-slip boundary condition on the
surface and the pinned contact line. Therefore, the oscillation damping rate
is higher for the PCL than the FCL cases.

Figure 8: Decay of oscillation amplitude for (a) different contact line
mobility and (b) different contact angles for $\text{Bo}=0$. The Lamb damping
rate (Eq. (3)) for the $n=2$ mode is shown for comparison.

The viscous damping of oscillation amplitude for different $\theta_{0}$ are
shown in Fig. 8(b). The damping rate generally decreases with increasing
$\theta_{0}$. As $\theta_{0}$ increases, the contact area decreases and the
constraint to the drop shape deformation is reduced. In the limit of
$\theta_{0}\to$$$, the drop will not deform due to the $n=1$ mode. As a
result, there will be no viscous dissipation due to the shape oscillation and
$\beta\to 0$.

###  III.5 Effect of gravitational Bond number

The results discussed so far are only for $\text{Bo}=0$, which represents the
oscillation dynamics of supported drops in a zero-gravity environment. The
results can serve as approximations for tiny drops with very small Bo.
However, for the droplet considered ($R_{d}=5$ mm), when gravity is present,
the equilibrium shape of the drop will significantly deviate from the
spherical cap. In the present study, we allow $g$ and Bo to vary from negative
to positive values. For $\text{Bo}>0$, the drop will be flattened, while for
$\text{Bo}<0$, the drop will be elongated. The addition of the hydrostatic
pressure will also change the pressure balance at the drop surface. The radius
of curvature of the drop at equilibrium state will not be constant as for
Bo=0. The gravity effect is shown to influence both the oscillation frequency
and the kinetic energy of the supported drop, though its effect on the viscous
damping rate seems to be very minor.

Figure 9: Variation of the supported drop first-mode oscillation frequency as
functions of (a) $\theta_{0}$ and (b) Bo for the PCL condition.

The simulation results of the first-mode frequencies for supported drops for
PCL and different $\theta_{0}$ and Bo are shown in Fig. 9. It is observed from
Fig. 9(a) that for all Bo, the decreasing trends of $\omega/\omega_{c}$ over
$\theta_{0}$ are similar. Furthermore, the oscillation frequencies increases
monotonically with Bo from negative to positive values for all $\theta_{0}$.
In other words, the oscillation frequency increases with gravity for sessile
drops, while for pendant drops, the frequency will decrease due to the gravity
effect. The trend of variation of the oscillation frequency is consistent with
former observations for both pendant and sessile drops Basaran and DePaoli
(1994); Sakakeeny and Ling (2020). Comparing the cases with zero and full
gravity, i.e., $\text{Bo}=0$ and 0.88, it is seen that the rise of frequency
due to the gravity effect increases with $\theta_{0}$.

To better illustrate the change in frequency due to the gravity effect, the
oscillation frequency is normalized by that for $\text{Bo}=0$ and is plotted
as a function of Bo in Fig. 9(b). Again, the monotonic increase of
$\omega/\omega_{c}$ over Bo from -0.088 to 0.88 can be clearly seen. The rate
of increase, indicated by the slopes of the curves, generally decreases with
Bo.

It is further observed that, the rate of change of $\omega/\omega_{Bo=0}$ over
Bo is more significant for large $\theta_{0}$. For $\theta_{0}=$$$,
$\omega/\omega_{Bo}$ increases about 72% when Bo increases from 0 to 0.088. In
contrast, with the same increase of Bo, $\omega/\omega_{Bo=0}$ for
$\theta_{0}=$$$ only increases less than 2%. It is also worth noting that for
$\theta_{0}=$$$, the frequency decreases quite rapidly with Bo when
$\text{Bo}<0$. The results indicate that the oscillation frequency for drops
supported by hydrophobic/super-hydrophobic surfaces can be very sensitive to
the change of Bo.

The variation of the oscillation frequency with $\theta_{0}$ and Bo can be
explained by the inviscid theoretical model developed in our previous study
Sakakeeny and Ling (2020), in which the first-mode oscillation can be modeled
as a mass-spring harmonic oscillator,

$k(z_{c}-z_{c,1})+m\frac{d^{2}(z_{c}-z_{c,1})}{dt^{2}}=0,$ (17)

where $k$ and $m$ are the effective spring constant and drop mass, while
$z_{c,1}$ represents the equilibrium centroid location for finite Bo. For
Bo=0, $z_{c,1}=z_{c,0}$. Here, the restoring force, $k(z_{c}-z_{c,1})$, is
mainly due to surface tension. It was shown that $k\sim\eta$, where $\eta$ is
a parameter that characterizes the increase of the drop surface area as the
centroid deviates from the equilibrium position, namely
$(S-S_{1})/S_{0}=\eta((z_{c}-z_{c,1})/R_{0})^{2}$, where $S_{0}$ and $S_{1}$
are the equilibrium drop surface area for zero and finite Bo. As a result, the
oscillation frequency

$\omega^{2}=k/m\sim\eta\,.$ (18)

The $S$-$z_{c}$ curve and $\eta$ for a given $\theta_{0}$ and Bo can be
estimated by the equilibrium drop theory, see Ref. Sakakeeny and Ling (2020)
for details.

The equilibrium drop theory indicates that $\eta$ monotonically increases with
Bo for all $\theta_{0}$. In other words, when Bo increases, the drop
equilibrium shape deviates from the spherical cap, and the increase of surface
area $(S-S_{1})$ for a given centroid deviation $(z_{c}-z_{c,1})$ becomes
higher. As a consequence, the restoring force increases, and thus too the
oscillation frequency increase. Similarly, the rate of increase of $\eta$ with
Bo also increases with $\theta_{0}$. For supported drops with large
$\theta_{0}$, the surface area increase is more “responsive” to the centroid
deviation and the change of Bo, due to the smaller contact area and constraint
from the surface. Therefore, the difference between the surface area
increments for the same $(z_{c}-z_{c,1})$ for zero and finite Bo, namely
$(S-S_{1})-(S-S_{1})_{\text{Bo=1}}$, is higher for larger $\theta_{0}$.
Correspondingly, the increases in both the restoring force and the frequency
are also magnified as $\theta_{0}$ increases.

Figure 10: Variation of the kinetic energy correction factor $\zeta$ as
functions of (a) $\theta_{0}$ and (b) Bo for the PCL condition.

The gravity effect also modifies the fluid flow induced by shape oscillation
and the drop kinetic energy correction factor $\zeta$. The simulation results
of $\zeta$ for different $\theta_{0}$ and Bo and for the PCL condition are
shown in Fig. 10. Similar to the oscillation frequency, $\zeta$ decreases with
$\theta_{0}$ for all Bo and increases with Bo for all $\theta_{0}$. Yet unlike
$\omega$, the increase of $\zeta$ due to the rise of Bo is more profound for
smaller $\theta_{0}$. Since the value of $\zeta$ also decreases as
$\theta_{0}$ decreases, it ends up that the normalized results, i.e.,
$\zeta/\zeta_{Bo=0}$, for different $\theta_{0}$ collapse approximately, see
Fig. 10(b), and can be fit via a linear function of Bo as

$\zeta(\theta_{0},\text{Bo})=\zeta_{Bo=0}(\theta_{0})(1+\alpha\text{Bo})\,.$
(19)

For the PCL results for all $\theta_{0}$ and Bo, the fitting yields
$\alpha=0.432$. Similar scaling behavior has been observed for the FCL
condition Sakakeeny and Ling (2020), where $\alpha=0.358$. Therefore, the
increase of the kinetic energy correction factor over Bo is faster when the
contact line changes from the FCL to the PCL condition.

For all $\theta_{0}$ and Bo, the oscillation frequency $\omega/\omega_{c}$ and
the kinetic energy correction factor $\zeta$ for PCL are always larger than
their FCL counterparts. The difference between the PCL and FCL values of
$\omega/\omega_{c}$ and $\zeta$, normalized by the difference at Bo=0, i.e.,
$(\omega_{P}-\omega_{F})/(\omega_{P}-\omega_{F})_{\text{Bo}=0}$ and
$(\zeta_{P}-\zeta_{F})/(\zeta_{P}-\zeta_{F})_{\text{Bo}=0}$, are plotted in
Fig. 11. For both variables, the results for different $\theta_{0}$ collapse
approximately. It is further observed that
$(\omega_{P}-\omega_{F})/(\omega_{P}-\omega_{F})_{\text{Bo}=0}$ varies little
with Bo. The collapsed results approximately follow a linear function passing
through the point (0,1) with a very small slope (about 0.10). This implies
that, though the frequencies $\omega_{P}$ and $\omega_{F}$ increase over Bo,
the difference between the two actually changes little. As a result, the
difference between oscillation frequencies for the PCL and FCL at Bo=0,
$(\omega_{P}-\omega_{F})_{\text{Bo}}=0$, see Fig. 5(c), is a good
approximation for non-zero Bo cases. On the other hand, the results for
$(\zeta_{P}-\zeta_{F})/(\zeta_{P}-\zeta_{F})_{\text{Bo}=0}$ for different
$\theta_{0}$ also approximately collapse and agree with a linear function, but
the slope is bigger than that for the frequency, i.e., about 0.690. As a
result, the difference between the PCL and FCL values of $\zeta$, i.e.,
$(\zeta_{P}-\zeta_{F})_{\text{Bo}=0}$ as shown in Fig. 7(c), needs to be
corrected using the results shown in Fig. 11(b) to represent the non-zero Bo
cases.

Figure 11: The differences between the values corresponding to the PCL and FCL
conditions for (a) $\omega$ and (b) $\zeta$ as functions of $\theta_{0}$ and
Bo.

##  IV Conclusions

The axisymmetric natural oscillations of a liquid drop supported by a flat
surface have been studied by direct numerical simulation. The parameters,
including the equilibrium contact angle ($\theta_{0}$) and the gravitational
Bond number (Bo), are varied to systematically investigate their effects on
the oscillation frequency and the induced flow around the drop. The two
asymptotic limits of contact line hysteresis and mobility, i.e., the pinned
and free contact line conditions, are considered to investigate the effect of
contact line mobility on oscillation. For the pinned contact line (PCL)
condition, the drop contact angle can vary freely, but the contact line cannot
move. For the free contact line (FCL) condition, the drop contact angle is
fixed, while the contact line is allowed to move freely. The results of
oscillation frequencies for these two limiting cases also serve as the upper
and lower bounds for general contact line conditions. In total, over 198
simulation cases were performed to study a wide range of equilibrium contact
angles ($$$\leq\theta_{0}\leq$$$) and Bond numbers ($-0.088\leq\text{Bo}\leq
0.88$) for both the FCL and PCL conditions. The negative and positive Bo
represent the pendant and sessile drops, respectively.

The drop oscillation is initiated by changing the gravity for a short period
of time. The first oscillation mode due to the drop centroid translation is
observed to dominate the excited oscillations. The oscillation frequency
$\omega$ scales with the capillary frequency $\omega_{c}$, and the normalized
frequency $\omega/\omega_{c}$ decreases with $\theta_{0}$. Remarkable
agreement between the simulation results with the inviscid theory of Bostwick
and Steen Bostwick and Steen (2014) is achieved, which validates the present
simulations. The shape oscillations induce flows within the drop that
contributes to the kinetic energy of the drop. The kinetic energy correction
factor $\zeta$ is defined as the ratio between the total kinetic energy of the
drop and that for the bulk motion. Similar to $\omega$, $\zeta$ also decreases
with $\theta_{0}$. The viscous damping rate $\beta$ of the oscillation
amplitude is also observed to decrease with $\theta_{0}$.

When Bo increases from -0.088 to 0.88, both $\omega/\omega_{c}$ and $\zeta$
increase. The increase in $\omega/\omega_{c}$ due to the rise of gravity
becomes more profound for larger $\theta_{0}$, indicating that the drop
oscillation frequency for hydrophobic/superhydrophobic surfaces can be quite
sensitive to the gravity effect. In contrast, the increase of $\zeta$ due to
gravity is more significant for small $\theta_{0}$. Furthermore, the results
of $\zeta$ for different $\theta_{0}$ collapse to a linear function if they
are normalized by the values at zero Bo. For all $\theta_{0}$ and Bo, the
values of $\omega/\omega_{c}$ and $\zeta$ for PCL are always greater than
their respective FCL values. The difference between the frequencies for FCL
and PCL, $\omega_{P}-\omega_{F}$ for different $\theta_{0}$ scales with the
counterpart for $\text{Bo}=0$, and the normalized difference,
$(\omega_{P}-\omega_{F})/(\omega_{P}-\omega_{F})_{\text{Bo}=0}$, varies little
with Bo.

## Acknowledgement

This work was supported by the startup fund at Baylor University and the
National Science Foundation (1853193, 1942324). The Baylor High Performance
and Research Computing Services (HPRCS) have provided the computational
resources that have contributed to the research results reported in this
paper. The authors also acknowledge Dr. Stéphane Popinet for the contribution
to the development of the Basilisk code.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Yao  _et al._ (2017) C.-W. Yao, C.-L. Lai, J. L. Alvarado, J. Zhou, K. T. Aung, and J. E. Mejia, “Experimental study on effect of surface vibration on micro textured surfaces with hydrophobic and hydrophilic materials,” Appl. Surf. Sci. 412, 45–51 (2017).
  * Dai  _et al._ (2018) X. Dai, N. Sun, S. O. Nielsen, B. B. Stogin, J. Wang, S. Yang, and T.-S. Wong, “Hydrophilic directional slippery rough surfaces for water harvesting,” Sci. Adv. 4, eaaq0919 (2018).
  * Boreyko and Chen (2009) J. B. Boreyko and C.-H. Chen, “Restoring superhydrophobicity of lotus leaves with vibration-induced dewetting,” Phys. Rev. Lett. 103, 174502 (2009).
  * Rayleigh (1879) L. Rayleigh, “On the capillary phenomena of jets,” Proc. R. Soc. Lond. 29, 71–97 (1879).
  * Lamb (1932) H. Lamb,  _Hydrodynamics_ (Cambridge University Press, 1932).
  * Trinh and Wang (1982) E. Trinh and T. G. Wang, “Large-amplitude free and driven drop-shape oscillations: experimental observations,” J. Fluid Mech. 122, 315–338 (1982).
  * Tsamopoulos and Brown (1983) J. A. Tsamopoulos and R. A. Brown, “Nonlinear oscillations of inviscid drops and bubbles,” J. Fluid Mech. 127, 519–537 (1983).
  * Basaran (1992) O. A. Basaran, “Nonlinear oscillations of viscous liquid drops,” J. Fluid Mech. 241, 169–198 (1992).
  * Helenbrook and Edwards (2002) B. T. Helenbrook and C. F. Edwards, “Quasi-steady deformation and drag of uncontaminated liquid drops,” Int. J. Multiphase Flow 28, 1631–1657 (2002).
  * Bergeles  _et al._ (2018) K. Bergeles, Y. Hardalupas, and A. M. K. P. Taylor, “On the transient flow inside and around a deforming millimetre class oil droplet falling under the action of gravity in stagnant air,” Phys. Fluids 30, 013305 (2018).
  * Zhang and Ling (2019) B. Zhang and Y. Ling, “High-fidelity modeling and simulation of primary breakup fo a gasoline surrogate jet,” in  _Proc. ILASS-Americas 30th Annual Conference on Liquid Atomization and Spray Systems 30th Annual Conference on Liquid Atomization and Spray Systems_ (2019).
  * Strani and Sabetta (1984) M. Strani and F. Sabetta, “Free vibrations of a drop in partial contact with a solid support,” J. Fluid Mech. 141, 233–247 (1984).
  * Sakakeeny and Ling (2020) J. Sakakeeny and Y. Ling, “Natural oscillations of a sessile drop on flat surfaces with mobile contact lines,” Phys. Rev. Fluids 5, 123604 (2020).
  * Noblin  _et al._ (2004) X. Noblin, A. Buguin, and F. Brochard-Wyart, “Vibrated sessile drops: Transition between pinned and mobile contact line oscillations,” Eur. Phys. J. E 14, 395–404 (2004).
  * Chang  _et al._ (2013) C.-T. Chang, J. B. Bostwick, P. H. Steen, and S. Daniel, “Substrate constraint modifies the rayleigh spectrum of vibrating sessile drops,” Phys. Rev. E 88, 023015 (2013).
  * Bostwick and Steen (2014) J. B. Bostwick and P. H. Steen, “Dynamics of sessile drops. part 1. inviscid theory,” J. Fluid Mech. 760, 5–38 (2014).
  * Bonn  _et al._ (2009) D. Bonn, J. Eggers, J. Indekeu, J. Meunier, and E. Rolley, “Wetting and spreading,” Rev. Mod. Phys 81, 739 (2009).
  * He and Yao (2019) P. He and C.-W. Yao, “Simulating contact angle hysteresis using pseudo-line tensions,” MRS Commun. 9, 1060–1066 (2019).
  * Snoeijer and Andreotti (2013) J. H. Snoeijer and B. Andreotti, “Moving contact lines: scales, regimes, and dynamical transitions,” Annu. Rev. Fluid Mech. 45, 269–292 (2013).
  * Basaran and DePaoli (1994) O. A. Basaran and D. W. DePaoli, “Nonlinear oscillations of pendant drops,” Phys. Fluids 6, 2923–2943 (1994).
  * (21) S. Popinet, “The basilisk code.” Available from http://basilisk.fr/.
  * Scardovelli and Zaleski (1999) R. Scardovelli and S. Zaleski, “Direct numerical simulation of free-surface and interfacial flow,” Annu. Rev. Fluid Mech. 31, 567–603 (1999).
  * Popinet (2009) S. Popinet, “An accurate adaptive solver for surface-tension-driven interfacial flows,” J. Comput. Phys. 228, 5838–5866 (2009).
  * Unverdi and Tryggvason (1992) S. O. Unverdi and G. Tryggvason, “A front-tracking method for viscous, incompressible, multi-fluid flows,” J. Comput. Phys. 100, 25–37 (1992).
  * Sussman  _et al._ (1994) M. Sussman, P. Smereka, and S. Osher, “A level set approach for computing solutions to incompressible two-phase flow,” J. Comput. Phys. 114, 146–159 (1994).
  * Renardy and Renardy (2002) Y. Renardy and M. Renardy, “PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method,” J. Comput. Phys. 183, 400–421 (2002).
  * Zhang  _et al._ (2019) B. Zhang, Y. Ling, P.-H. Tsai, A.-B. Wang, S. Popinet, and S. Zaleski, “Short-term oscillation and falling dynamics for a water drop dripping in quiescent air,” Phys. Rev. Fluids 4, 123604 (2019).
  * Marcotte and Zaleski (2019) F. Marcotte and S. Zaleski, “Density contrast matters for drop fragmentation thresholds at low ohnesorge number,” Phys. Rev. Fluids 4, 103604 (2019).
  * Zhang  _et al._ (2020) B. Zhang, S. Popinet, and Y. Ling, “Modeling and detailed numerical simulation of the primary breakup of a gasoline surrogate jet under non-evaporative operating conditions,” Int. J. Multiphase Flow 130, 103362 (2020).
  * Mostert and Deike (2020) W. Mostert and L. Deike, “Inertial energy dissipation in shallow-water breaking waves,” J. Fluid Mech. 890, A12 (2020).
  * Afkhami and Bussmann (2009) S. Afkhami and M. Bussmann, “Height functions for applying contact angles to 3d VOF simulations,” Int. J. Numer. Meth. Fluids 61, 827–847 (2009).
  * van Hooft  _et al._ (2018) J. A. van Hooft, S. Popinet, C. C. van Heerwaarden, S. J. A. van der Linden, S. R. de Roode, and B. J. H. van de Wiel, “Towards adaptive grids for atmospheric boundary-layer simulations,” Bound.-Layer Meteor. 167, 421–443 (2018).


# An Adaptive Receiver for Underwater Acoustic Full-Duplex Communication with
Joint Tracking of the Remote and Self-Interference Channels
††thanks:  This work was supported primarily by the Engineering Research
Centers Program of the National Science Foundation under NSF Cooperative
Agreement No. (CNS-1704097, CNS-1704076). Any opinions, findings and
conclusions or recommendations expressed in this material are those of the
authors and do not necessarily reflect those of the National Science
Foundation.

Mohammad Towliat1, Zheng Guo2, Leonard J. Cimini3, Xiang-Gen Xia4, and Aijun
Song5

1,3,4Department of Electrical and Computer Engineering
University of Delaware
Newark, DE, USA
{mtowliat, cimini<EMAIL_ADDRESS>2,5Department of Electrical and
Computer Engineering
University of Alabama
Tuscaloosa, AL, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Full-duplex (FD) communication is a promising candidate to address the data
rate limitations in underwater acoustic (UWA) channels. Because of
transmission at the same time and on the same frequency band, the signal from
the local transmitter creates self-interference (SI) that contaminates the the
signal from the remote transmitter. At the local receiver, channel state
information for both the SI and remote channels is required to remove the SI
and equalize the SI-free signal, respectively. However, because of the rapid
time-variations of the UWA environment, real-time tracking of the channels is
necessary. In this paper, we propose a receiver for UWA-FD communication in
which the variations of the SI and remote channels are jointly tracked by
using a recursive least squares (RLS) algorithm fed by feedback from the
previously detected data symbols. Because of the joint channel estimation, SI
cancellation is more successful compared to UWA-FD receivers with separate
channel estimators. In addition, due to providing a real-time channel tracking
without the need for frequent training sequences, the bandwidth efficiency is
preserved in the proposed receiver.

###### Index Terms:

Bandwidth efficiency, in-band full-duplex, self-interference cancellation,
time-varying channels, underwater acoustic communication.

##  I Introduction

Limited available bandwidth and dynamic transmission media are two serious
issues that challenge underwater acoustic (UWA) communication systems and lead
to a data transmission rate which does not exceed a few tens of kilobits per
second [1]. Full-duplex (FD) communication is an interesting technique to
increase the data rate in UWA channels. By performing transmission and
reception at the same time and on the same frequency band, UWA-FD has the
potential to double the data rate of half-duplex systems. However, it is
challenging to eliminate the self-interference (SI) caused by the local
transmitter. A combination of analog and digital methods can be used for SI
cancellation [2]; it has also been shown that fully-digital SI cancellation
methods are capable of minimizing SI [3]. In fully-digital approaches,
assuming knowledge of the SI channel, the receiver is able to estimate and
eliminate the SI; then, the SI-free signal can be equalized with the purpose
of detecting the data symbols from a remote transmitter. Obviously, this
equalizer is designed to correct for the impulse response of the remote
channel that links the receiver and the remote transmitter. So, channel state
information for both the SI and remote channels is required at the receiver;
however, since the SI signal is much more powerful than the remote signal, an
accurate tracking of the SI channel is critical.

Because of fluctuations in the UWA transmission media, the SI and remote
channels are time-varying [4, 5], and estimates are valid for just a short
time, roughly equivalent to the coherence time of the channel. So, real-time
tracking of both time-varying channels is necessary. Addressing this concern,
in this paper, we propose an adaptive UWA-FD receiver with three stages of
processing at each cycle: channel estimation, SI cancellation, and
equalization. At the channel estimation stage, the SI and remote channels are
jointly estimated; at the next two stages, these estimates are used to cancel
the SI and equalize the SI-free signal, respectively.

There are a few recent related works using the adaptive filtering approach to
perform SI cancellation in UWA channels [6, 7]. In these works, only the SI
channel is estimated by the adaptive filters, and the remote channel is
assumed to be known or estimated in a further separate step. Since the SI is
much stronger than the remote signal, the SI channel is estimated by treating
the remote signal as an additional additive noise. Then, the SI signal is
removed from the received signal based on the estimated SI channel. Obviously,
during this procedure, the SI channel estimation and, consequently, the SI
removal is affected by the unknown remote signal, which leads to a performance
degradation [8]. By contrast, in our proposed receiver, the SI and remote
channels are jointly estimated which implies that, during the channel
estimation stage, both the SI and remote signals are taken into consideration.
As a result, one of the advantages of this approach is that the SI channel
estimation is more accurate which, consequently, results in a better SI
cancellation. The other advantage of the proposed receiver is that the remote
channel is tracked along with the SI channel without using any extra training
sequences. This means that variations in the remote channel are tracked by the
receiver without wasting the bandwidth needed for recurrent training sequences
from the remote transmitter.

The rest of this paper is organized as follows. In Section II, the SI in UWA-
FD system is described and formulated. The proposed UWA-FD receiver is
presented in Section III. Simulation results are given in Section IV and we
conclude the paper in Section V.

##  II Problem Formulation

In Fig. 1, a model of the baseband interference at the reviver is shown; the
SI signal interferes with the desired remote signal. In this model, $i[n]$ and
$x[n]$ are the reference sequences for the local and remote transmitters,
respectively. The local transmitter’s reference, $i[n]$, is the downconverted,
matched filtered, and downsampled (to the symbol rate) output of the local
power amplifier (PA). Because of the non-linearity and noise of the PA [9],
$i[n]$ is not necessarily equal to the local transmitted symbols. However, the
remote transmitter’s reference, $x[n]$, represents the transmitted symbols
from the remote transmitter. The first reason for these considerations is
that, since the SI signal is relatively powerful, the non-linearity and noise
effects of the local PA are not negligible at the receiver; by contrast,
because of the lower power of the remote signal at the receiver, the
distortions of the remote PA are considered trivial [10]. Secondly, it must be
noted that the pulse shapes that are used to shape the transmitted symbols at
the local transmitter are not aligned with the matched filter at the receiver
because the receiver’s matched filter is paired with the pulse shapes from the
remote transmitter. The misalignment between the receiver and the local
transmitter causes an additional distortion on the transmitted symbols from
the local transmitter. In order to take this into account, the matched filter
that is used to generate $i[n]$ must be aligned with the receiver, not with
the local transmitter. According to the above discussion, $i[n]$ can take any
soft value, whereas, $x[n]$ is the actual data symbol generated by the remote
modulator.

Figure 1: Model for baseband interference at the receiver.

In the model presented in Fig. 1, in order to consider the multipath fading in
the UWA transmission media,
$\mathbf{c}[n]={{\left[{{c}_{0}}[n],\ldots{{c}_{M-1}}[n]\right]}^{H}}$and
$\mathbf{h}[n]={{\left[{{h}_{\text{0}}}[n],\ldots{{h}_{L-1}}[n]\right]}^{H}}$
are the time-varying SI and remote channel vectors with length $M$ and $L$,
respectively. By defining the signal vectors
$\mathbf{i}[n]={{\left[i[n],\ldots i[n-M+1]\right]}^{T}}$and
$\mathbf{x}[n]={{\left[x[n],\ldots x[n-L+1]\right]}^{T}}$, the received signal
can be represented as

$y[n]=s[n]+r[n]+w[n],$ (1)

where $w[n]$ denotes the additive ambient noise with power $\sigma_{0}^{2}$;
$s[n]$ and $r[n]$ are respectively the SI and remote signals, with power
$P_{s}$ and $P_{r}$, and can be defined as

$\displaystyle s[n]$ $\displaystyle={\mathbf{c}}^{H}[n]\,\mathbf{i}[n],$
$\displaystyle r[n]$ $\displaystyle={\mathbf{h}}^{H}[n]\,\mathbf{x}[n].$ (2)

Obviously, in (1) and (2), it is assumed that the channel coherence time is
long enough that both the SI and remote channels remain invariant over
$K=\max\\{M,L\\}$ time instances.

The ultimate goal at the receiver is to detect the data symbol $x[n]$ from the
received signal. Since $i[n]$ is known, the simplest approach is to determine
$s[n]$ and remove it from the received signal, $y[n]$, and then employ an
equalizer to detect the data symbols $x[n]$ from the SI-free signal. An
estimate of the SI channel is required for SI cancellation, and an estimate of
the remote channel is needed for equalization. On the other hand, because of
fluctuations in the UWA transmission media, both the SI and remote channels
are time-varying and, hence, the receiver should be able to keep track of
these changes without wasting the bandwidth caused by using frequent training
intervals. To respond to this concern, in this paper, we propose a receiver in
which the variations of the SI and remote channels are jointly tracked without
requiring recurrent training sequences from the remote transmitter. We remove
the necessity of transmitting training by taking feedback from the previously
detected data symbols and using it as a reference sequence for the remote
transmitter.

##  III The proposed UWA-FD receiver

The block diagram in Fig. 2 depicts the proposed UWA-FD receiver.
Specifically, at the $n\text{th}$ cycle, the following tasks are performed: i)
Estimations of both the SI and the remote channel vectors (i.e.
$\mathbf{\hat{c}}[n]$ and $\mathbf{\hat{h}}[n]$) are simultaneously updated.
ii) Using the provided SI channel estimate from the first stage, the $n$th
sample of the SI signal, $\hat{s}[n]$, is estimated and eliminated from the
received signal $y[n]$ to obtain the SI-free signal $\hat{r}[n]$. iii) The
estimated remote channel from the first stage, passes through a damper, which
is used to suppress the undesired fluctuations of the estimates. By using a
damped version of the estimated remote channel, the SI-free signal,
$\hat{r}[n]$, is equalized to detect the data symbol
$\hat{x}[n-\Delta^{\prime}]$, where $\Delta^{\prime}$ is the delay of the
equalizer. Finally, the detected data symbol is used as a reference in the
channel estimation stage at the next cycle with $\Delta=\Delta^{\prime}+1$
delays. In the following, we individually explain the stages of the proposed
receiver.

Figure 2: Proposed UWA-FD receiver at the $n$th cycle.

###  III-A Channel Estimation Stage

Channel estimation is one of the main concerns in UWA-FD systems; with a
precise estimate of the channels, the receiver is able to accurately remove
the SI and, then, detect the data symbols from the SI-free signal. However, it
must be noted that, since the SI channel is significantly stronger than the
remote channel, the accuracy of the SI channel estimation is more critical. In
addition, because of the time-variations of the UWA transmission media, the
channel estimator must provide a real-time estimate of both the SI and remote
channels in order to keep track of the variations [4].

As mentioned before, the SI and remote channels are jointly estimated at the
channel estimation stage. At this stage, the known signal $i[n]$ is used as
the reference sequence for estimating the SI channel; while, for estimating
the remote channel, the remotely transmitted data symbol, $x[n]$, is the
reference. However, since $x[n]$ is unknown, we use the previously detected
data symbols as the reference.

Amongst the various adaptive methods, we use the recursive least squares (RLS)
algorithm which has been shown to be a good estimator for time-varying
channels [11]. At the $n\text{th}$ cycle, the input of the RLS algorithm
includes the received signal $y[n-\Delta]$, the SI reference sequence
$\mathbf{i}[n-\Delta]$, and the previously detected data symbols
$\mathbf{\hat{x}}[n-\Delta]=\left[\hat{x}[n-\Delta],\ldots\hat{x}[n-\Delta-L+1]\right]^{T}$.

The reference sequences $\mathbf{\hat{x}}[n-\Delta]$ and
$\mathbf{i}[n-\Delta]$ are nested in a single column vector as
${\bf{v}}[n-\Delta]={[{\bf{\hat{x}}}[n-\Delta]^{T}\,\,{\mkern 1.0mu}{\mkern
1.0mu}{\mkern 1.0mu}\vdots\,\,{\mkern 1.0mu}{\mkern 1.0mu}{\mkern
1.0mu}{\bf{i}}[n-\Delta]^{T}]^{T}}$ and the estimates of the channel vectors
are also concatenated in a column vector as
${\bf{u}}[n-\Delta]={[{{\bf{\hat{h}}}^{H}}[n-\Delta]\,\,\,\,\vdots\,\,\,{{\bf{\hat{c}}}^{H}}[n-\Delta]]^{H}}$.
The aim of the RLS estimator is to update the previous vector
$\mathbf{u}[n-\Delta-1]$ and obtain $\mathbf{u}[n-\Delta]$. To this end, the
estimated received signal is (see (1) and (2))

$\hat{y}[n-\Delta]={{\mathbf{u}}^{H}}[n-\Delta-1]\,\mathbf{v}[n-\Delta].$ (3)

Then, the error $e=y[n-\Delta]-\hat{y}[n-\Delta]$ is used to update
$\mathbf{u}[n-\Delta-1]$. The RLS algorithm for the joint channel estimation
at the $n\text{th}$ cycle is presented in Algorithm 1, in which $\lambda$ is
the forgetting factor and $\delta$ is a small scalar for initiation.

if $n=0$ then

Set $\lambda$ and $\delta$

Initiate ${\bf{P}}\leftarrow\frac{1}{\delta}{{\bf{I}}_{M+L}}$

else

Input : ${\bf{u}}[n-\Delta-1]$ and ${\bf{v}}[n-\Delta]$

$\hat{y}[n-\Delta]\leftarrow{{\bf{u}}^{H}}[n-\Delta-1]\,{\bf{v}}[n-\Delta]$

$e\leftarrow y[n-\Delta]-\hat{y}[n-\Delta]$

${\bf{\gamma}}\leftarrow{(\lambda+{{\bf{v}}^{H}}[n-\Delta]\,{\bf{Pv}}[n-\Delta])^{-1}}{\bf{Pv}}[n-\Delta]$

${\bf{P}}\leftarrow{\lambda^{-1}}\left({{\bf{I}}-{\bf{\gamma}}{{\bf{v}}^{H}}[n-\Delta]}\right)\,{\bf{P}}$

${\bf{u}}[n-\Delta]\leftarrow{\bf{u}}[n-\Delta-1]+{e^{*}}{\bf{\gamma}}$

Return : ${\bf{u}}[n-\Delta]$

end if

Algorithm 1 : RLS algorithm to jointly estimate the SI and remote channels at
the $n$th cycle

According to (3), because of the equalizer delay in data symbol detection, the
channel estimation is also delayed, such that, at the $n\text{th}$ cycle, the
channel vectors at time $n-\Delta$ are estimated (i.e.
$\mathbf{\hat{c}}[n-\Delta]$ and $\mathbf{\hat{h}}[n-\Delta]$). In other
words, the channel estimation stage is $\Delta$ time instances behind the SI
cancellation and equalization stages, where $\mathbf{\hat{c}}[n]$ and
$\mathbf{\hat{h}}[n]$ are required, respectively. However, recall from (1) and
(2) that the time-variations are assumed to be slow enough so that the SI and
remote channels remain constant over, at least, $K=\max\\{M,L\\}$ time
instances (i.e. $\mathbf{c}[n]\approx\mathbf{c}[n-K]$ and
$\mathbf{h}[n]\approx\mathbf{h}[n-K]$). Therefore, by adjusting $\Delta$ to
satisfy $0\leq\Delta\leq K$, we can use the approximations
$\mathbf{\hat{c}}[n]\approx\mathbf{\hat{c}}[n-\Delta]$ and
$\mathbf{\hat{h}}[n]\approx\mathbf{\hat{h}}[n-\Delta]$, so that the estimated
channel at the channel estimation stage can be utilized at the next two
stages.

As mentioned above, in the RLS algorithm, both the SI and remote signals are
used to construct $\hat{y}[n-\Delta]$, as shown in (3). This implies that for
estimating the SI channel, we not only use the strong SI signal, but we also
take the history of the weak remote signal into account. In the simulations,
we will show that the proposed procedure results in a more accurate SI channel
estimation compared to other adaptive receivers in which the remote signal is
treated as an additive noise during the SI channel estimation.

The other point is that, in the proposed RLS algorithm, there is no need to
transmit frequent training sequences from the remote transmitter; instead, we
use the previously detected data symbols as the reference for estimating the
remote channel. Hence, the bandwidth efficiency is preserved in the proposed
receiver. However, to get a faster initial convergence, in practice, we use a
short training sequence to initiate the RLS algorithm. Based on the simulation
results, the length of the initial training should be a little longer than
$M+L$ symbols, which is small and does not significantly reduce the bandwidth
efficiency of the proposed receiver during a long-term utilization of the
system.

###  III-B Self-Interference Cancellation Stage

The next stage, as shown in Fig. 2, is SI cancellation. By using the SI
reference sequence and the estimate of the SI channel, provided from the first
stage, the SI signal is reconstructed and subtracted from the received signal
as

$\displaystyle\hat{s}[n]={{{\mathbf{\hat{c}}}}^{H}}[n]\,\mathbf{i}[n];$
$\displaystyle\hat{r}[n]=y[n]-\hat{s}[n],$ (4)

where $\hat{s}[n]$ and $\hat{r}[n]$ are the estimates of the SI and the SI-
free signals, respectively. Since the power of the SI signal is significantly
higher than that of the remote signal, it is clear that an accurate estimate
of the SI channel is required in (4). To provide some intuition, consider an
example in which the power of $i[n]$ and $x[n]$ are ${{P}_{i}}={{P}_{x}}=0$
dB. The power of the SI and remote channel vectors are
$\text{E}\left\|\mathbf{c}[n]\right\|^{2}=10$ dB and
$\text{E}\left\|\mathbf{h}[n]\right\|^{2}=-20$ dB, respectively. Therefore,
the power of the SI and remote signals at the receiver become $P_{s}=10$ dB
and $P_{r}=-20$ dB, respectively. A small deviation between the actual and
estimated SI channel,
$\text{E}\left\|\mathbf{c}[n]-\mathbf{\hat{c}}[n]\right\|^{2}=-20$ dB, results
in a residual SI with power ${{P}_{\text{res}}}=-20$ dB, which is in the same
range as the remote signal power, $P_{r}$. This significant residual SI,
caused by a small inaccuracy in the estimated channel, will significantly
deteriorate the performance.

###  III-C Equalization Stage

The last stage in the proposed receiver is equalization, in which the estimate
of the SI-free signal, $\hat{r}[n]$, is equalized with the purpose of
eliminating the multipath effects of the remote channel from the data symbols.
In the proposed receiver, we use a decision feedback equalizer (DFE), which is
a non-linear equalizer; by employing a feedforward filter (FFF) and a feedback
filter (FBF), a DFE outperforms the linear zero forcing (ZF) and minimum mean
squared error (MMSE) equalizers in frequency selective channels [12]. In a
DFE, the FFF with length $\mathcal{L}$ acts as linear equalizer for the first
$\mathcal{L}$ paths of the channel and the FBF removes the inter-symbol-
interference caused by the rest of the paths and the FFF. Also a DFE with
these properties, leads to a $\Delta^{\prime}={\cal L}-1$ symbol delay. So, at
the $n$th cycle, the output of the equalization stage is
$\hat{x}[n-\Delta^{\prime}]$, which is the detected version of
$x[n-\Delta^{\prime}]$.

Since the DFE removes the effects of multipath for remote channel
${\bf{h}}[n]$, the FFF and FBF weights are tuned based on the estimate of this
channel [13]. Recall that the estimation of the remote channel,
${\bf{\hat{h}}}[n]$, is already provided from the channel estimation stage.
However, because of possible errors in the estimated remote channel, the
direct use of ${\bf{\hat{h}}}[n]$ to tune the DFE filters is not optimal.
Since the remote channel is much weaker than the SI channel, in the joint
channel estimation by the RLS algorithm, its estimate is susceptible to the
noise and the errors at the previously detected data symbols. These unwanted
effects lead to incorrect and artificial fluctuations in ${\bf{\hat{h}}}[n]$
at each cycle. Hence, by directly using ${\bf{\hat{h}}}[n]$ to tune the DFE,
the noise and symbol errors propagate to the next cycles and deteriorate the
performance.

In order to prevent this damage and suppress the propagation of the artificial
fluctuations in ${\bf{\hat{h}}}[n]$ to the next cycles, we use
${\bf{\bar{h}}}[n]$ to tune the FFF and FBF at the DFE, where
${\bf{\bar{h}}}[n]$ is the damped version of the estimated channel and is
given as

${\bf{\bar{h}}}[n]=(1-{{\mu}}){\bf{\bar{h}}}[n-1]+{{\mu}}{\bf{\hat{h}}}[n],$
(5)

where $0\leq{{\mu}}\leq 1$ is the damping factor, chosen so that the
artificial and rapid fluctuations in ${\bf{\hat{h}}}[n]$ are damped; however,
the slower natural changes of the channel are not suppressed because the
adjustment of the DFE filters is supposed to follow the natural time-
variations of the channel at each cycle. Having said that, one can see from
(5) that if ${{\mu}}=0$, then ${\bf{\bar{h}}}[n]={\bf{\bar{h}}}[n-1]$, which
means that any changes, including the artificial and natural changes, at
${\bf{\hat{h}}}[n]$ are extremely damped and ${\bf{\bar{h}}}[n]$ remains
constant. On the other hand, if ${{\mu}}=1$,
${\bf{\bar{h}}}[n]={\bf{\hat{h}}}[n]$, indicating that no damping is imposed.
So, the value of ${{\mu}}$ is determined based on the speed of the natural
time-variations of the actual remote channel ${\bf{h}}[n]$. If ${\bf{h}}[n]$
is changing slowly with a long coherence time, ${{\mu}}$ is set to a small
value; however, if ${\bf{h}}[n]$ is rapidly time-varying with a shorter
coherence time, ${{\mu}}$ is set to a larger scalar.

##  IV Simulation Results

In this section we numerically evaluate the performance of the proposed
receiver for UWA-FD system. We consider a bandwidth $B=5\,\text{kHz}$ and use
normalized power BPSK modulation at both the local and remote ends.

The simulations are performed in baseband; however, in order to model the non-
linearity of the local PA, we first generate the passband input of the PA and
then, downconvert the output to baseband. To this end, we consider a carrier
frequency ${{f}_{\text{c}}}=12\,\,\text{kHz}$ and a root raised cosine filter,
with the roll-off-factor $\alpha=0.5$, to shape the signal; the length of the
filter is truncated to $12$ symbols duration. The odd harmonics of the Taylor
series expansion are used to relate the input (i.e. $p(t)$) and output (i.e.
$q(t)$) of the PA as $q(t)=\sum\limits_{m=1,3,5}{{{a}_{m}}p{{(t)}^{m}}}$ [10],
where ${{a}_{1}}=100$, ${{a}_{3}}=5$ and ${{a}_{5}}=10$. The local PA noise,
${w_{{\rm{PA}}}}(t)$ has a Gaussian distribution with power
$\sigma_{\text{PA}}^{2}=10\,\text{dB}$. Then, the local transmitter reference
sequence, $i[n]$, is the downconverted, matched filtered, and downsampled
version of $i(t)=q(t)+{w_{{\rm{PA}}}}(t)$, such that the power of $i[n]$ is
normalized to $P_{i}=0$ dB.

The baseband power delay profile (PDP) of the SI channel is
$\text{PDP}_{{{\rm{SI}}}}[k]={\beta_{1}}{\mkern 1.0mu}f[k]$, where $f[k]$,
shown in Fig. 3, is the estimated PDP for the SI channel in a lake experiment
in [4]. The remote channel is assumed to have an exponential PDP as
$\text{PDP}_{{\text{remote}}}[k]={{\beta}_{2}}\,{{e}^{-0.25\,k}}$. The factors
${{\beta}_{1}}$ and ${{\beta}_{2}}$ are chosen so that the power of the SI and
remote signals are ${{P}_{s}}$ and ${{P}_{r}}$, respectively. The lengths of
the SI and remote channels are $M=30$ and $L=70$, respectively. All paths of
both channels are time-varying with a coherence time $70$ ms, except the
direct path of the SI channel which is relatively time-invariant [4].

Figure 3: The estimated PDP for the SI channel in a lake experiment [4]. The
direct path (first peak) is stable with a long coherence time; however, the
reflected path from the surface (second peak) is time-varying with a coherence
time of $70$ ms.

At the receiver side, the forgetting factor for the RLS algorithm is set to
$\lambda=0.98$, and $\delta=0.0001$. At the DFE, we take the FFF length as
${\cal L}=70$ and the FBF length as 50. With this consideration, the equalizer
delay for delivering the data symbols becomes $\Delta^{\prime}=69$. Based on
the coherence time of the remote channel, we choose the damping factor as
${\mu}={10^{-3}}$. In addition, the initial training sequence from the remote
transmitter contains $130$ symbols, or $130/B=26\,\text{ms}$.

We assume the powers of the SI signal, remote signal and ambient noise are
${{P}_{s}}=0\,\text{dB}$, ${{P}_{r}}=-20\,\,\text{dB}$ and
$\sigma_{0}^{2}=-35$ dB, respectively. We compute the normalized mean square
error (MSE) for the SI channel estimation in the proposed receiver (in which
the SI and remote channels are jointly estimated) and the conventional
receiver (where the SI channel is estimated by ignoring the remote signal) as

$\begin{array}[]{l}{\rho_{{\bf{\hat{c}}}}}=\frac{{{\rm{E}}{{\left\|{{\bf{c}}[n]-{\bf{\hat{c}}}[n]}\right\|}^{2}}}}{{{\rm{E}}{{\left\|{{\bf{c}}[n]}\right\|}^{2}}}}=9.75\times{10^{-5}};\\\
\\\
{\rho_{{\bf{\tilde{c}}}}}=\frac{{{\rm{E}}{{\left\|{{\bf{c}}[n]-{\bf{\tilde{c}}}[n]}\right\|}^{2}}}}{{{\rm{E}}{{\left\|{{\bf{c}}[n]}\right\|}^{2}}}}=2.9\times{10^{-3}},\end{array}$
(6)

where ${{\bf{\hat{c}}}}[n]$ and ${{\bf{\tilde{c}}}}[n]$ are the SI channel
vector estimated by the proposed and the conventional receivers, respectively.
Clearly, ${\rho_{{\bf{\hat{c}}}}}$ is much smaller than
${\rho_{{\bf{\tilde{c}}}}}$, demonstrating the advantage of the proposed
receiver in tracking the SI channel and, eventually, in canceling the SI .

In order to illustrate the damper’s role in suppressing the artificial
fluctuations and reducing the remote channel estimation errors, as an example,
we have plotted one realization of the damper’s input (which is the estimated
remote channel by the RLS algorithm) and output (the channel that is used to
tune the DFE) in Fig. 4. This figure shows the $10$th estimated and damped
delay path of the remote channel denoted by ${{\hat{h}}_{10}}[n]$ and
${{\bar{h}}_{10}}[n]$, respectively, and compares them with that of the true
remote channel denoted by $h_{10}[n]$. As seen, the RLS estimated channel has
artificial fluctuations over the subsequent cycles. However, damping this
estimate kills these artificial fluctuations and passes the smooth and natural
time-variations. By looking at Fig. 4, one can intuitively see that the
channel estimation error for the damper output is lower than that of the
input. To evaluate that, we calculate the normalized MSE for the input and
output of the damper,

$\begin{array}[]{l}{\rho_{{\bf{\hat{h}}}}}=\frac{{{\rm{E}}{{\left\|{{\bf{h}}[n]-{\bf{\hat{h}}}[n]}\right\|}^{2}}}}{{{\rm{E}}{{\left\|{{\bf{h}}[n]}\right\|}^{2}}}}=1.6\times{10^{-2}},\\\
\\\
{\rho_{{\bf{\bar{h}}}}}=\frac{{{\rm{E}}{{\left\|{{\bf{h}}[n]-{\bf{\bar{h}}}[n]}\right\|}^{2}}}}{{{\rm{E}}{{\left\|{{\bf{h}}[n]}\right\|}^{2}}}}=1.6\times{10^{-3}}.\end{array}$
(7)

Comparing ${\rho_{{\bf{\bar{h}}}}}$ with ${\rho_{{\bf{\hat{h}}}}}$ indicates
that the remote channel estimation with RLS suffers from a significant error.
As discussed in Section III-C, this large error is caused by the noise and
errors in the previously detected data symbols. However, the relatively lower
error in the output of the damper confirms that using a damper corrects the
remote channel estimation to some extent, before it is used to tune the DFE.

Figure 4: Amplitude and phase of the $10$th path of the true (dashed red
line), estimated (dash-dotted black line), and damped (solid blue line) remote
channels.

In Fig. 5, we present the BER versus the remote SNR (i.e.
$P_{r}/\sigma_{0}^{2}$), when $P_{s}=0$ dB, $P_{r}=-20$ dB and
$\sigma_{0}^{2}$ is in the range from $-50$ to $-10$ dB. In this figure, the
BER of the proposed and conventional receivers are compared with that of the
ideal receiver, in which both the SI and remote channel state information is
assumed to be perfectly known. Since in the ideal receiver the SI channel is
perfectly known, the SI is completely removed and the SI-free signal is
equalized based on the perfect knowledge of the remote channel. By looking at
Fig. 5, one can see that the BER of the proposed receiver is almost identical
to that of the ideal receiver, implying that the proposed receiver can track
the variations of both SI and remote channels very well. However, the
performance of the conventional receiver is significantly deteriorated. This
is because of two reasons. First, in this receiver, the SI channel estimation
is not as accurate as in the proposed receiver (as shown in (6)) and a
considerable residual SI contaminates the output of the SI cancellation stage.
Second, the variation of the remote channel is not tracked so that the
estimate for the remote channel remains fixed to that obtained by using the
initial training sequence.

In Figs. 6 and 7, we show the BER when $P_{s}=5$ and $10$ dB, respectively. As
seen, the BER of the proposed receiver degrades as the SI power increases. As
a matter of fact, for an even stronger SI, the SI cancellation residual is
larger and leads to performance degradation. In order to investigate the
impact of the SI power on the residual, in Fig. 8, we plot ${\rho_{\hat{r}}}$,
the normalized MSE of $\hat{r}[n]$, versus the SI to remote signal power ratio
($P_{s}/P_{r}$). We assume $P_{r}=-20$ dB and $\sigma_{0}^{2}=-35$ dB, and the
$P_{s}$ range is $-20$ to 20 dB. The normalized MSE for $\hat{r}[n]$ is
calculated as

${\rho_{\hat{r}}}=\frac{{{\rm{E|}}r[n]-\hat{r}[n]{|^{2}}}}{{{\rm{E|}}r[n]{|^{2}}}}.$
(8)

This quantity measures the normalized residual power after SI cancellation. As
seen from Fig. 8, for $P_{s}/P_{r}\leq 20$ dB, ${\rho_{\hat{r}}}$ is very low
for the proposed receiver, indicating that the performance of this receiver is
identical to that of the ideal receiver (as seen in Fig. 5, where
$P_{s}/P_{r}=20$ dB); however, for $P_{s}/P_{r}>20$ dB, a stronger residual
results, which leads to a higher BER (as seen in Figs. 6 and 7, where
$P_{s}/P_{r}=25$ dB and $30$ dB, respectively). In addition,
${\rho_{\hat{r}}}$ of the conventional receiver for $P_{s}/P_{r}\leq 20$ dB is
a significant constant value and surges for $P_{s}/P_{r}>20$ dB. For
$P_{s}/P_{r}=40$ dB, ${\rho_{\hat{r}}}$ for the proposed receiver equals to
that of the conventional receiver.

Figure 5: BER comparison for the proposed, conventional and ideal receivers
($P_{s}=0$ dB and $P_{r}=-20$ dB). Figure 6: BER comparison for the proposed,
conventional and ideal receivers ($P_{s}=5$ dB and $P_{r}=-20$ dB). Figure 7:
BER comparison for the proposed, conventional and ideal receivers ($P_{s}=10$
dB and $P_{r}=-20$ dB). Figure 8: Comparison on ${\rho_{\hat{r}}}$ versus
$P_{s}/P_{r}$ for the proposed and conventional receivers ($P_{r}=-20$ dB and
$\sigma_{0}^{2}=-35$ dB).

##  V Conclusions

In this paper, we proposed a novel receiver for UWA-FD communication with
joint tracking of the SI and remote channels. In order to remove the need to
frequently send training sequences, the channel estimator is fed by feedback
from the previously detected data symbols as a reference for the remote
channel. With this strategy, an accurate SI channel estimation is obtained;
also, the remote channel variations are tracked along with the SI channel
without using extra training. The concern about error propagation due to using
feedback is addressed by employing a damper to suppress the harmful effects of
the incorrectly detected data symbols in the estimated remote channel. Then,
the output of the channel damper is used to tune the equalizer. Simulation
results show that the proposed receiver outperforms conventional methods which
only track the SI channel. The performance of the proposed receiver decreases
as the SI power increases beyond a specific value. This value is determined
based on the channel PDPs, the coherence time and the noise level.

## References

  * [1] S. Sendra, J. Lloret, J. M. Jimenez, and L. Parra, “Underwater acoustic modems,” _IEEE Sensors Journal_ , vol. 16, no. 11, pp. 4063–4071, Jun 2016.
  * [2] L. Li, A. Song, L. J. Cimini, X.-G. Xia, and C.-C. Shen, “Interference cancellation in in-band full-duplex underwater acoustic systems,” in _Proc. OCEANS 2015 - MTS/IEEE_ , Oct 2015, pp. 1–6.
  * [3] G. Qiao, S. Gan, S. Liu, L. Ma, and Z. Sun, “Digital self-interference cancellation for asynchronous in-band full-duplex underwater acoustic communication,” _Sensors_ , vol. 18, no. 6, p. 1700, Jun 2018.
  * [4] M. Towliat, Z. Guo, L. J. Cimini, X.-G. Xia, and A. Song, “Self-interference channel characterization in underwater acoustic in-band full-duplex communications using OFDM,” _arXiv preprint arXiv:2005.10933_ , 2020.
  * [5] B. Tomasi, J. Preisig, G. B. Deane, and M. Zorzi, “A study on the wide-sense stationarity of the underwater acoustic channel for non-coherent communication systems,” in _Proc. 17th European Wireless - Sustainable Wireless Technologies_ , Apr 2011, pp. 1–6.
  * [6] L. Shen, B. Henson, Y. Zakharov, and P. Mitchell, “Digital self-interference cancellation for full-duplex underwater acoustic systems,” _IEEE Trans. on Circuits and Syst. II: Express Briefs_ , vol. 67, no. 1, pp. 192–196, Jan 2020.
  * [7] L. Shen, B. Henson, Y. Zakharov, and P. Mitchell, “Two-stage self-interference cancellation for full-duplex underwater acoustic systems,” in _Proc. OCEANS 2019_ , Jun 2019, pp. 1–6.
  * [8] G. Qiao, S. Gan, S. Liu, and Q. Song, “Self-interference channel estimation algorithm based on maximum-likelihood estimator in in-band full-duplex underwater acoustic communication system,” _IEEE Access_ , vol. 6, pp. 62 324–62 334, Oct 2018.
  * [9] L. Shen, B. Henson, Y. Zakharov, and P. Mitchell, “Robust digital selfinterference cancellation for full-duplex UWA systems: Lake experiments,” in _Proc. Underwater Acoustics Conference and Exhibition_ , Jul 2019, pp. 243–250.
  * [10] D. Bharadia, E. McMilin, and S. Katti, “Full duplex radios,” in _Proc. ACM SIGCOMM_ , pp. 375-386, Aug 2013.
  * [11] S. Haykin, A. H. Sayed, J. R. Zeidler, P. Yee, and P. C. Wei, “Adaptive tracking of linear time-variant systems by extended RLS algorithms,” _IEEE Trans. on Signal Process._ , vol. 45, no. 5, pp. 1118–1128, May 1997.
  * [12] A. Agarwal, S. Sur, A. K. Singh, H. Gurung, A. K. Gupta, and R. Bera, “Performance analysis of linear and non-linear equalizer in Rician channel,” _Procedia Technology_ , vol. 4, pp. 687 – 691, Jan 2012.
  * [13] M. Stojanovic, J. G. Proakis, and J. A. Catipovic, “Analysis of the impact of channel estimation errors on the performance of a decision-feedback equalizer in fading multipath channels,” _IEEE Trans. on Commun._ , vol. 43, no. 2/3/4, pp. 877–886, Feb 1995.