# High-Order Implicit Large-Eddy Simulation of Flow over a Marine Propeller

Bin Zhang<EMAIL_ADDRESS>Chi Ding  Chunlei Liang  Department of
Mechanical and Aerospace Engineering, The George Washington University,
Washington, DC 20052 USA  Department of Mechanical and Aeronautical
Engineering, Clarkson University, Potsdam, NY 13676 USA

(March 24, 2021)

###### Abstract

We report the first high-order eddy-resolving simulation of flow over a marine
propeller using a recently developed high-order sliding-mesh method. This
method employs the flux reconstruction framework and a new dynamic curved
mortar approach to handle the complex rotating geometries. For a wide range of
working conditions, it is validated to predict the loads very accurately
against experiments. The method’s low-dissipation characteristic has allowed
the capturing of a broad spectrum of turbulence structures for very long
distances even on a very coarse grid. Comparison with a previous low-order
simulation is also carried out to show the low-dissipation advantage of the
present simulations. From detailed load analysis, the major loads and their
distributions and time and frequency scales are identified. Visualizations of
the instantaneous, phase-averaged, and time-averaged flow fields have revealed
the processes of tip vortex formation, major vortex evolutions, and flow
instability developments at different working conditions. The effects of
different fairwaters on the propeller’s overall performance are also
quantitatively assessed.

###### keywords:

high-order method , flux reconstruction , sliding mesh , implicit LES , marine
propeller

##  1 Introduction

Numerical techniques for studying marine propellers have seen immense advances
in the past decades. For example, the lifting surface procedures based on
vortex lattice and panel methods have been widely used in propeller design. A
representative program in this category is the PSF code developed at MIT in
the 1980’s [1]. The vortex lattice and panel methods have very quick
turnaround time, and can predict propeller loads rather accurately when
combined with empirical vortex and wake models. They are, however, incapable
of revealing the fine details of a flow field for in-depth analysis of flow
physics such as flow instability and acoustics. Since the 1990’s, the
Reynolds-averaged Navier–Stokes (RANS) methods have gained popularity in
simulating propeller flows. The RANS methods do provide richer flow field
information than the inviscid vortex lattice and panel methods. But the
averaging nature of a RANS (even an unsteady RANS) method still smooths out a
lot of flow details, especially the small and intermediate instantaneous
eddies. In more recent years, the increasing high-performance computing power
has allowed eddy-resolving techniques, such as detached-eddy simulation (DES)
and large-eddy simulation (LES), to have been performed on marine propellers
in the scale of millions or even trillions of grid elements. For instance,
Muscari et al. [2] and Di Mascio et al. [3] employed DES to study the vortex
dynamics of a propeller in different flow conditions. Verma et al. [4] and
Jang and Mahesh [5] used LES to investigate the flow around a reverse rotating
marine propeller. Balaras et al. [6] applied an immersed boundary based LES
technique to explore the flow around a propeller with and without an upstream
appendage. Kumar and Mahesh [7] systematically studied the wake instability of
a propeller using LES. However, all the aforementioned DES and LES approaches
are of low-order accuracies (at most second-order) that could potentially
introduce large numerical dissipation and dispersion to a flow field.

On the other hand, high-order (third and above) methods, especially polynomial
based high-order methods, are seeing rapid development in recent years. High-
order methods show many advantages over the traditional low-order ones. For
example, for the same number of degrees of freedom, a high-order method
produces solutions with much smaller numerical error than a low-order method
does. Furthermore, a high-order method can employ high-order unstructured
curved meshes that approximate curved geometries much more accurately than the
linear meshes for a low-order method. The most popular high-order methods
include the discontinuous Galerkin (DG) method [8, 9], the spectral element
method [10, 11, 12], the spectral volume method [13, 14], the spectral
difference (SD) method [15, 16, 17, 18, 19, 20], etc. Among these methods, the
SD method solves equations in differential form directly, and is one of the
most efficient high-order methods. Recently, the ideas of collocating solution
and flux points of the SD method and correcting fluxes using higher-degree
polynomials have led to an even more efficient high-order method — the flux
reconstruction (FR) method [21, 22], also known as the correction procedure
via reconstruction (CPR) [23]. Besides its better efficiency, by choosing
different correction polynomials, the FR method can recover many existing
high-order schemes such as DG and SD, and can even produce new schemes that
had never been reported before. The stability of the FR method has been proved
in [24]. The most recent developments on the FR method are summarized in [25].

To the authors’ knowledge, there is still no reported high-order eddy-
resolving simulation of marine propellers. One important reason is the severe
challenge on how to treat the complex rotational geometry of a propeller in a
high-order method without deteriorating the method’s accuracies in both space
and time. To tackle this challenge, Zhang and Liang [26, 27] developed a high-
order sliding-mesh method for the SD and FR methods by introducing the concept
of curved dynamic mortar elements. A parallelization approach was proposed for
this method in [28], and the extension to three-dimensions was achieved in
[29]. More recently, this method was further extended to sliding interfaces
with general nonuniform meshes [30]. An updated version that is high-order in
time, arbitrarily high-order in space, provably conservative, and provably
outflow preservative has also been established [31]. In the present work, we
apply this method to implicit LES (without using an explicit sub-grid-scale
model) [32, 33, 34, 35, 36] of the loads and flow fields of a real marine
propeller at various working conditions. This is the first time a high-order
method being applied to simulate a marine propeller.

The rest of this paper is organized as follows. Section 2 gives a brief
description of the flow equations and the numerical methods. Section 3 details
the simulation setup. Simulation results and discussions are reported in
Section 4. Finally, Section 5 concludes this paper.

##  2 Numerical methods

###  2.1 The physical equations

We numerically solve the following three-dimensional unsteady Navier-Stokes
equations in a conservative form,

$\frac{\partial\mathbf{Q}}{\partial t}+\frac{\partial\mathbf{F}}{\partial
x}+\frac{\partial\mathbf{G}}{\partial y}+\frac{\partial\mathbf{H}}{\partial
z}=\boldsymbol{0},$ (1)

where $\mathbf{Q}$ is the vector of conservative variables, $\mathbf{F}$,
$\mathbf{G}$ and $\mathbf{H}$ are the flux vectors in each coordinate
direction. These terms have the following expressions,

$\displaystyle\mathbf{Q}$ $\displaystyle=[\rho~{}\rho u~{}\rho v~{}\rho
w~{}E]^{\mathsf{T}},$ (2) $\displaystyle\mathbf{F}$
$\displaystyle=\mathbf{F}_{\text{inv}}(\mathbf{Q})+\mathbf{F}_{\text{vis}}(\mathbf{Q},\nabla\mathbf{Q}),$
(3) $\displaystyle\mathbf{G}$
$\displaystyle=\mathbf{G}_{\text{inv}}(\mathbf{Q})+\mathbf{G}_{\text{vis}}(\mathbf{Q},\nabla\mathbf{Q}),$
(4) $\displaystyle\mathbf{H}$
$\displaystyle=\mathbf{H}_{\text{inv}}(\mathbf{Q})+\mathbf{H}_{\text{vis}}(\mathbf{Q},\nabla\mathbf{Q}),$
(5)

where $\rho$ is fluid density, $u$, $v$ and $w$ are the velocity components,
$E$ is the total energy per volume defined as
$E=p/(\gamma-1)+\frac{1}{2}\rho(u^{2}+v^{2}+w^{2})$, $p$ is pressure, $\gamma$
is the ratio of specific heats and is set to 1.4.

The fluxes have been split into inviscid and viscous parts. The inviscid
fluxes are only functions of the conservative variables and have the following
expressions,

$\mathbf{F}_{\text{inv}}=\left[\begin{array}[]{c}\rho u\\\ \rho u^{2}+p\\\
\rho uv\\\ \rho uw\\\ u(E+p)\end{array}\right],\
\mathbf{G}_{\text{inv}}=\left[\begin{array}[]{c}\rho v\\\ \rho uv\\\ \rho
v^{2}+p\\\ \rho vw\\\ v(E+p)\end{array}\right],\
\mathbf{H}_{\text{inv}}=\left[\begin{array}[]{c}\rho w\\\ \rho uw\\\ \rho
vw\\\ \rho w^{2}+p\\\ w(E+p)\end{array}\right].$ (6)

The viscous fluxes are functions of the conservative variables and their
gradients. The expressions are

$\mathbf{F}_{\text{vis}}=-\left[\begin{array}[]{c}0\\\ \tau_{xx}\\\
\tau_{yx}\\\ \tau_{zx}\\\ u\tau_{xx}+v\tau_{yx}+w\tau_{zx}+\kappa
T_{x}\end{array}\right],$ (7)
$\mathbf{G}_{\text{vis}}=-\left[\begin{array}[]{c}0\\\ \tau_{xy}\\\
\tau_{yy}\\\ \tau_{zy}\\\ u\tau_{xy}+v\tau_{yy}+w\tau_{zy}+\kappa
T_{y}\end{array}\right],$ (8)
$\mathbf{H}_{\text{vis}}=-\left[\begin{array}[]{c}0\\\ \tau_{xz}\\\
\tau_{yz}\\\ \tau_{zz}\\\ u\tau_{xz}+v\tau_{yz}+w\tau_{zz}+\kappa
T_{y}\end{array}\right],$ (9)

where $\tau_{ij}$ is shear stress tensor which is related to velocity
gradients as $\tau_{ij}=\mu(u_{i,j}+u_{j,i})+\lambda\delta_{ij}u_{k,k}$, $\mu$
is dynamic viscosity, $\lambda=-2/3\mu$ based on Stokes’ hypothesis,
$\delta_{ij}$ is the Kronecker delta, $\kappa$ is thermal conductivity, $T$ is
temperature that is related to density and pressure through the ideal gas law
$p=\rho\mathcal{R}T$, where $\mathcal{R}$ is the gas constant.

###  2.2 The computational equations

As will be discussed later, we map each moving grid element from the physical
space to a stationary standard cubic element in the computational space where
the equations are solved. Assume the mapping is: $t=\tau$,
$x=x(\tau,\xi,\eta,\zeta)$, $y=y(\tau,\xi,\eta,\zeta)$ and
$z=z(\tau,\xi,\eta,\zeta)$, where $(\tau,\xi,\eta,\zeta)$ are the
computational time and space. It can be shown that the flow equations will
take the following conservative form in the computational space,

$\frac{\partial\widetilde{\mathbf{Q}}}{\partial
t}+\frac{\partial\widetilde{\mathbf{F}}}{\partial\xi}+\frac{\partial\widetilde{\mathbf{G}}}{\partial\eta}+\frac{\partial\widetilde{\mathbf{H}}}{\partial\zeta}=\boldsymbol{0},$
(10)

and the computational variable and fluxes are related to the physical ones
through

$\left[\begin{array}[]{c}\widetilde{\mathbf{Q}}\\\ \widetilde{\mathbf{F}}\\\
\widetilde{\mathbf{G}}\\\
\widetilde{\mathbf{H}}\end{array}\right]=|\mathcal{J}|\mathcal{J}^{-1}\left[\begin{array}[]{c}\mathbf{Q}\vphantom{\widetilde{\mathbf{Q}}}\\\
\mathbf{F}\vphantom{\widetilde{\mathbf{F}}}\\\
\mathbf{G}\vphantom{\widetilde{\mathbf{G}}}\\\
\mathbf{H}\vphantom{\widetilde{\mathbf{H}}}\end{array}\right],$ (11)

where $|\mathcal{J}|$ is determinant of the Jacobian matrix and
$\mathcal{J}^{-1}$ is the inverse Jacobian matrix, and their expressions are

$\displaystyle|\mathcal{J}|=\left|\frac{\partial(t,x,y,z)}{\partial(\tau,\xi,\eta,\zeta)}\right|=\left|\begin{array}[]{cccc}1&0&0&0\\\
x_{t}&x_{\xi}&x_{\eta}&x_{\zeta}\\\ y_{t}&y_{\xi}&y_{\eta}&y_{\zeta}\\\
z_{t}&z_{\xi}&z_{\eta}&z_{\zeta}\end{array}\right|,$ (16)
$\displaystyle\mathcal{J}^{-1}=\frac{\partial(\tau,\xi,\eta,\zeta)}{\partial(t,x,y,z)}=\left[\begin{array}[]{cccc}1&0&0&0\\\
\xi_{t}&\xi_{x}&\xi_{y}&\xi_{z}\\\ \eta_{t}&\eta_{x}&\eta_{y}&\eta_{z}\\\
\zeta_{t}&\zeta_{x}&\zeta_{y}&\zeta_{z}\\\ \end{array}\right].$ (21)

Besides the flow equations, the Geometric Conservation Law (GCL) [37] also
needs to be numerically satisfied to ensure free-stream preservation on moving
grids. The GCL equations and the steps for solving them are described in our
previous papers [27, 29].

###  2.3 The flux reconstruction method

The meshes in this work consist of hexahedral elements only. The first step of
the FR method is to map each hexahedral element to a standard cubic element of
unit size, i.e., $0\leq\xi,\eta,\zeta\leq 1$. This can be done using the
following iso-parametric mapping,

$\left[\begin{array}[]{c}x\\\ y\\\
z\end{array}\right]=\sum_{i=1}^{K}M_{i}(\xi,\eta,\zeta)\left[\begin{array}[]{c}x_{i}(t)\\\
y_{i}(t)\\\ z_{i}(t)\end{array}\right],$ (22)

where $K$ is the number of nodes of a hexahedral element, $M_{i}$ (detailed
expressions can be found in, for example, [38]) is the shape function of the
$i$-th node, and $(x_{i},y_{i},z_{i})$ are the coordinates of the $i$-th node.

Next, solution points (SPs, denoted by $X_{s}$) and flux points (FPs, denoted
by $X_{f}$) are defined along each coordinate direction in the standard
element. Fig. 1 shows a schematic of the distribution of the SPs and FPs in
the $\xi$-$\eta$ plane for a fourth-order FR method. Generally, for an $N$-th
order FR scheme, there are $N$ SPs and FPs in each direction, where the SPs
are in the interior and the FPs are on the boundaries of the standard element.
The SPs and FPs are chosen as the Legendre points in this work.

Figure 1: Schematic of solution points (circles) and flux points (squares) in
the $\xi\textnormal{-}\eta$ plane for a fourth-order FR method.

Subsequently, Lagrange interpolation bases are defined at each SP. For
example, at the $i$-th SP we have

$h_{i}(X)=\prod_{s=1,s\neq i}^{N}\left(\frac{X-X_{s}}{X_{i}-X_{s}}\right).$
(23)

The resulting bases also form a basis of polynomials of degrees less than or
equal to $N-1$, i.e., $\boldsymbol{\mathsf{P}}_{N\\!-\\!1}$. These
interpolation bases allow the construction of solution and flux polynomials
inside each element through tensor products. For example,

$\displaystyle\widetilde{\mathbf{Q}}(\xi,\eta,\zeta)$
$\displaystyle=\sum_{k=1}^{N}\sum_{j=1}^{N}\sum_{i=1}^{N}\widetilde{\mathbf{Q}}_{ijk}h_{i}(\xi)h_{j}(\eta)h_{k}(\zeta),$
(24) $\displaystyle\widetilde{\mathbf{F}}(\xi,\eta,\zeta)$
$\displaystyle=\sum_{k=1}^{N}\sum_{j=1}^{N}\sum_{i=1}^{N}\widetilde{\mathbf{F}}_{ijk}h_{i}(\xi)h_{j}(\eta)h_{k}(\zeta),$
(25)

where the subscript $ijk$ denotes the discrete value at the $ijk$-th SP. All
these polynomials are in
$\boldsymbol{\mathsf{P}}_{N\\!-\\!1,N\\!-\\!1,N\\!-\\!1}$.

The above solution and flux polynomials are continuous within each element,
but discontinuous across cell boundaries. Therefore, common values need to be
defined on cell boundaries. There are various ways to calculate these common
values. In this work, the common solution is calculated as the average of the
discontinuous values from the two sides of a boundary; the common inviscid
fluxes are computed using the Rusanov Riemann solver[39]; the common viscous
fluxes are computed from the common solutions and common gradients.

There is one more issue: after taking the first-order spatial derivatives in
the governing equations, the three flux terms become elements in
$\boldsymbol{\mathsf{P}}_{N\\!-\\!2,N\\!-\\!1,N\\!-\\!1}$,
$\boldsymbol{\mathsf{P}}_{N\\!-\\!1,N\\!-\\!2,N\\!-\\!1}$, and
$\boldsymbol{\mathsf{P}}_{N\\!-\\!1,N\\!-\\!1,N\\!-\\!2}$, respectively, and
are inconsistent with the solution term. To fix this issue, the original
fluxes need to be reconstructed. This is done by using correction functions
that are polynomials of degree no less than $N$. Taking the flux in the $\xi$
direction as an example, the reconstructed flux polynomials is

$\widehat{\mathbf{F}}(\xi,\eta,\zeta)=\widetilde{\mathbf{F}}(\xi,\eta,\zeta)+\big{[}\widetilde{\mathbf{F}}^{\text{com}}(0,\eta,\zeta)-\widetilde{\mathbf{F}}(0,\eta,\zeta)\big{]}\cdot
g_{\text{\tiny
L}}(\xi)+[\widetilde{\mathbf{F}}^{\text{com}}(1,\eta,\zeta)-\widetilde{\mathbf{F}}(1,\eta,\zeta)]\cdot
g_{\text{\tiny R}}(\xi)$ (26)

where $\widetilde{\mathbf{F}}(\xi,\eta,\zeta)$ is from (25);
$\widetilde{\mathbf{F}}^{\text{com}}$ represents the common flux on a cell
boundary; $g_{\text{\tiny L}}$ and $g_{\text{\tiny R}}$ are the left and right
correction functions, and are required to at least satisfy

$\displaystyle g_{\text{\tiny L}}(0)$ $\displaystyle=1,\quad g_{\text{\tiny
L}}(1)$ $\displaystyle=0,$ (27) $\displaystyle g_{\text{\tiny R}}(0)$
$\displaystyle=0,\quad g_{\text{\tiny R}}(1)$ $\displaystyle=1,$

which ensures that

$\widehat{\mathbf{F}}(0,\eta,\zeta)=\widetilde{\mathbf{F}}^{\text{com}}(0,\eta,\zeta),\quad\widehat{\mathbf{F}}(1,\eta,\zeta)=\widetilde{\mathbf{F}}^{\text{com}}(1,\eta,\zeta),$
(28)

i.e., the reconstructed fluxes still take the common values at cell
interfaces. In this work, we have employed the $g_{\text{\tiny DG}}$
correction function [21]. The other two fluxes are reconstructed in the same
way.

Finally, the governing equations can be written in the following residual
form,

$\left.\frac{\partial\widetilde{\mathbf{Q}}}{\partial
t}\right|_{ijk}=-\left[\frac{\partial\widehat{\mathbf{F}}}{\partial\xi}+\frac{\partial\widehat{\mathbf{G}}}{\partial\eta}+\frac{\partial\widehat{\mathbf{H}}}{\partial\zeta}\right]_{ijk}=\mathfrak{R}_{ijk},\quad
i,j,k=1,2,\cdots,N,$ (29)

where $\mathfrak{R}_{ijk}$ is the residual at the $(i,j,k)$-th SP. This system
can be time marched using either explicit or implicit temporal schemes.

###  2.4 A sliding-mesh method

In three-dimensions, we consider two types of sliding interfaces as depicted
in Fig. 2: one is annular, and the other is cylindrical. To simplify the
explanation, we have required the mesh to only unmatch in the azimuthal
direction but match in the radial (for annular sliding) or axial (for
cylindrical sliding) direction. We also require equal mesh size in the
azimuthal direction. These restrictions are imposed for explanation purposes
only and can be easily lifted in practice. More details can be found in our
previous papers [26, 27, 28, 29, 30, 31].

Figure 2: Two types of sliding meshes: left, annular sliding; right,
cylindrical sliding.

Since the SPs on the two sides of a nonconforming sliding interface do not
match, we aim to find the best-possible common values of a variable on the two
sides of a sliding interface. This can be achieved by least-squares
projections using mortar elements as the intermediate medium. A mortar element
is formed by the overlapping region of two cell faces. Taking the cylindrical
sliding interface as an example, based on the assumptions we have made, a cell
face has two mortar elements as sketched in Fig. 3. The first step is to map a
cell face and the mortars to standard ones as shown in the same figure, using,
for example, iso-parametric mapping or transfinite mapping.

Figure 3: Map curved cell face and mortars to straight ones.

Let the $(\xi^{\prime},\eta^{\prime})$ denote the mortar space, then the
computational and the mortar spaces are related as

$\xi=o+s\xi^{\prime},~{}~{}~{}\eta=\eta^{\prime},$ (30)

where $0\leq\xi,\eta,\xi^{\prime},\eta^{\prime}\leq 1$, and $o$ and $s$ are
the offset and scaling of a mortar with respect to a cell face.

Let $\phi$ represent the variable of interest, and obviously it can be
represented by the following polynomials on a cell face $\Omega$ and on the
left side of a mortar $\Xi$,

$\displaystyle\phi^{\Omega}(\xi,\eta)$
$\displaystyle=\sum_{j=1}^{N}\sum_{i=1}^{N}\phi_{ij}^{\Omega}h_{i}(\xi)h_{j}(\eta),$
(31) $\displaystyle\phi^{\Xi,L}(\xi^{\prime},\eta^{\prime})$
$\displaystyle=\sum_{j=1}^{N}\sum_{i=1}^{N}\phi_{ij}^{\Xi,L}h_{i}(\xi^{\prime})h_{j}(\eta^{\prime}),$
(32)

where $\phi_{ij}^{\Omega}$ and $\phi_{ij}^{\Xi,L}$ are the discrete values at
the $(i,j)$-th SP on $\Omega$ and the left side of $\Xi$, respectively. The
$(\phi_{ij}^{\Xi,L})$’s are unknown, and can be obtained through the following
projection (refer to Fig. 4(a)),

$\int_{0}^{1}\int_{0}^{1}(\phi^{\Xi,L}(\xi^{\prime},\eta^{\prime})-\phi^{\Omega}(\xi,\eta))h_{\alpha}(\xi^{\prime})h_{\beta}(\eta^{\prime})\text{d}\xi^{\prime}\text{d}\eta^{\prime}=0,~{}\forall\,\alpha,\beta=1,2,...,N.$
(33)

Considering the relations in (30), it can be shown that the above two-
dimensional projection is equivalent to the following one-dimensional one,

$\int_{0}^{1}(\phi^{\Xi,L}(\xi^{\prime},X_{j})-\phi^{\Omega}(\xi,X_{j}))h_{\alpha}(\xi^{\prime})\text{d}\xi^{\prime}=0,~{}\forall\,\alpha=1,2,...,N,$
(34)

where $X_{j}$ is the coordinate of the $j$-th SP. Evaluating the above
equation for all the $\alpha$’s, we will get a system of equations about
$\phi_{1:N,j}^{\Xi,L}$. Repeating this process for every $j$, we will obtain
every $\phi_{ij}^{\Xi,L}$. The values on the right side of a mortar, i.e., the
$(\phi_{ij}^{\Xi,R})$’s, can be obtained in the same way.

Figure 4: Projection between face and mortar: (a) from face to the left side
of a mortar, (b) from two mortars back to a face.

After that, we are able to compute a common value on the mortar, e.g., through
averaging or Riemann solver. Let us denote this common value as $\Phi^{\Xi}$.
We then project this common variable back to a cell face from mortars as
demonstrated in Fig. 4(b). And the projection is

$\begin{split}\int_{\xi=0}^{\xi=o_{2}}&\int_{\eta=0}^{\eta=1}(\Phi^{\Omega}(\xi,\eta)-\Phi^{\Xi_{1}}(\xi^{\prime},\eta^{\prime}))h_{\alpha}(\xi)h_{\beta}(\eta)\text{d}\xi\text{d}\eta\\\
&+\int_{\xi=o_{2}}^{\xi=1}\int_{\eta=0}^{\eta=1}(\Phi^{\Omega}(\xi,\eta)-\Phi^{\Xi_{2}}(\xi^{\prime},\eta^{\prime}))h_{\alpha}(\xi)h_{\beta}(\eta)\text{d}\xi\text{d}\eta=0,~{}\forall\,\alpha,\beta=1,2,...,N,\end{split}$
(35)

where $\Phi^{\Omega}$ represents the polynomial of the unknown common variable
on face $\Omega$. Similarly, this projection is equivalent to the following
one-dimensional projection,

$\int_{0}^{o_{2}}(\Phi^{\Omega}(\xi,X_{j})-\Phi^{\Xi_{1}}(\xi^{\prime},X_{j}))h_{\alpha}(\xi)\text{d}\xi+\int_{o_{2}}^{1}(\Phi^{\Omega}(\xi,X_{j})-\Phi^{\Xi_{2}}(\xi^{\prime},X_{j}))h_{\alpha}(\xi)\text{d}\xi=0,~{}\forall\,\alpha=1,2,...,N,$
(36)

from which $\Phi_{1:N,j}^{\Omega}$ are obtained, and then also every
$\Phi_{ij}^{\Omega}$ by repeating this process for $j$.

##  3 Simulation setup

The propeller studied in this work is the DTMB 4119 model designed at the
David Taylor Model Basin [40, 41, 42]. Its parameters are summarized in Tab.
1, and the geometry is given in Tab. 2, where $r$ represents radial position,
$R=D/2$ is propeller radius, $c$ is blade section chord length, $P$ is section
pitch, $\phi_{P}$ is section pitch angle, $t$ is section thickness, and $f$ is
section camber.

parameter |  value
---|---
Number of blades $Z$ |  3
Diameter $D$ [m] |  0.305
Hub diameter ratio $D_{h}/D$ |  $0.2$
Design advance ratio $J$ |  0.833
Rotation |  Right handed
Section thickness |  NACA66 modified
Section mean line |  NACA, $a=0.08$
Table 1: Design parameters of the DTMB 4119 propeller. $r/R$ |  $c/D$ |  $P/D$ |  $\phi_{P}$ [∘] |  $t/c$ |  $f/c$
---|---|---|---|---|---
0.2 |  0.3200 |  1.105 |  60.38 |  0.20550 |  0.01429
0.3 |  0.3625 |  1.102 |  49.47 |  0.15530 |  0.02318
0.4 |  0.4048 |  1.098 |  41.15 |  0.11800 |  0.02303
0.5 |  0.4392 |  1.093 |  34.84 |  0.09016 |  0.02182
0.6 |  0.4610 |  1.088 |  29.99 |  0.06960 |  0.02072
0.7 |  0.4622 |  1.084 |  26.24 |  0.05418 |  0.02003
0.8 |  0.4347 |  1.081 |  23.28 |  0.04206 |  0.01967
0.9 |  0.3613 |  1.079 |  20.88 |  0.03321 |  0,01817
0.95 |  0.2775 |  1.077 |  19.84 |  0.03228 |  0.01631
1.0 |  0.0 |  1.075 |  18.89 |  0.03160 |  0.01175
Table 2: Geometry of the DTMB 4119 propeller.

When the geometry is given, a propeller flow is governed by two nondimensional
parameters: advance ratio and Reynolds number. The advance ratio is defined as

$J=\frac{U_{\infty}}{nD},$ (37)

where $U_{\infty}$ is the incoming flow speed, $n$ is propeller’s revolution
per second (RPS), and $D$ is propeller diameter. The following Reynolds number
is usually adopted in experimental studies of marine propellers

$Re_{c}=\frac{c_{0.7}\,U_{0.7}}{\nu}=\frac{c_{0.7}\sqrt{U_{\infty}^{2}+(2\pi
0.7Rn)^{2}}}{\nu},$ (38)

where $c_{0.7}$ and $U_{0.7}=\sqrt{U_{\infty}^{2}+(2\pi 0.7Rn)^{2}}$ are the
chord length and the relative speed, respectively, at $r/R=0.7$, and $\nu$ is
fluid kinematic viscosity. For numerical simulations, it is more convenient to
define the Reynolds number as,

$Re_{D}=\frac{D\,U_{\infty}}{\nu}.$ (39)

It can be shown that these two Reynolds numbers are conveniently convertible
through the following relation,

$Re_{D}=\frac{Re_{c}}{\frac{c_{0.7}}{D}\sqrt{1+\left(\frac{0.7\pi}{J}\right)^{2}}}.$
(40)

Depending on $c_{0.7}/D$ and $J$, $Re_{D}$ could be either larger or smaller
than $Re_{c}$, but usually not much different.

The geometry is visualized in Fig. 5. As a typical screw-type propeller, DTMB
4119 consists of four components: the shaft, the blades, the hub, and the
fairwater. All components, except the shaft, rotate at an angular speed
$\omega$. Hub and fairwater are usually not part of propeller design. In this
work, the hub has a length $L_{h}=0.5D$, with the three blades installed
evenly along the circumferential direction of the mid-hub. The fairwater in
the figure is a $1:2$ ellipsoid, but cylindrical and hemispherical fairwaters
are also employed in a later section to study their effects.

Figure 5: Two views of the DTMB 4119 propeller: left, side view; right, back
view (pressure side).

The overall computational domain is cylindrical as shown in Fig. 6. It has a
length of $15D$ in the streamwise (i.e., $x$) direction and a diameter of
$12D$. The resulting blockage ratio of this domain with respect to the
propeller is $0.69\%$, which is small enough to guarantee negligible
confinement effects according to the study in [43]. The propeller locates $3D$
downstream from the inlet, with the shaft extends all the way to the inlet.
The blades and the hub are enclosed in a small sliding disk region whose
radius and thickness are $0.75D$ and $0.5D$, respectively. A global view of
the mesh (with a $1/4$ cutout to expose the propeller) is shown on the right
of Fig. 6, where the propeller and the sliding interfaces are colored in red.

Figure 6: Overall computational domain (left) and global view of the mesh
(right).

The overall mesh consists of about $235,000$ quadratic curved hexahedral
elements, of which about $36,000$ are within the small sliding disk region.
The mesh is refined around the propeller as well as in the wake region. The
first layer of the off wall elements on each blade surface has a height of
approximately $0.015D$, and the first off wall solution point is about
$0.0007D$ (for the fifth-order scheme) away from the walls. Two local views of
the meshes on the sliding interfaces and the blade surfaces are shown in Fig.
7. We see that the high-order curved mesh captures the curvatures of the blade
surfaces very well.

Figure 7: Close views of the mesh on the sliding interfaces (left) and on the
propeller surfaces (right).

We treat the inlet as a Dirichlet boundary, the outer cylindrical surface and
the outlet as characteristic farfields that allow waves and flow to leave
without reflection [44], and all solid surfaces as no-slip adiabatic walls.
The incoming freestream flow has a low Mach number of $Ma_{\infty}=0.05$ so
that compressibility effects are small. The Reynolds number is
$Re_{D}=5.59\times 10^{5}$, which is equivalent to $Re_{c}=7.3\times 10^{5}$.
Various advance coefficients are studied but with a focus on the design value
(i.e., $J=0.833$). The nondimensional angular speed of the propeller is
$\omega^{*}=\omega D/U_{\infty}=2\pi nD/U_{\infty}=2\pi/J$, and the rotation
period is $tU_{\infty}/D=J$. The simulations are performed in two ways: to
collect time-averaged statistics, the outer subdomain is fixed, only the inner
sliding region rotates at $\omega$, and a velocity boundary condition is
applied on the fairwater surface; to collect phase-averaged statistics, the
whole domain rotates at angular speed $\omega$, and velocity boundary
condition is applied to the shaft to make it stationary.

A four-stage third-order SSP-RK scheme [45] with a nondimensional time step
size of $\Delta tU_{\infty}/D=2.5\times 10^{-5}$ is adopted for the time
marching. The propeller therefore rotates about $0.01$ degree per time step at
the design condition. For spatial discretization, it is known that high-order
simulation of turbulent flow may experience instabilities due to aliasing
errors [24]. We observed such instabilities on the fifth- and above orders. To
overcome this issue, we have employed the filter reported in [46] (with
strength $\alpha=0.05$) to stabilize the simulations. Meanwhile, we compared
the design-condition results from the fourth-, the fifth-, and the sixth-order
schemes to ensure sufficient resolution. It was observed that the mean loads
from the fourth- and the fifth-order schemes have small differences (around
$3\%$), but the fifth-order scheme resolves the flow structures with more
details. On the other hand, both the mean loads and the flow fields from the
fifth- and the sixth-order schemes are almost indistinguishable, which
indicates that the fifth-order is the optimum choice considering both accuracy
and cost. For this reason, the fifth-order scheme has been used for the
simulations in what follows. It is also worth mentioning that all simulations
were run for a nondimensional time of $tU_{\infty}/D=100$, and phase- and
time-averaging were performed on the last $85$ time units, which represents
approximately 102 revolutions at the design condition.

##  4 Results and discussion

###  4.1 Propeller loads

The loads on a propeller are measured by the thrust and torque coefficients
defined as below,

$K_{T}=\frac{T}{\rho n^{2}D^{4}}\qquad\text{and}\qquad K_{Q}=\frac{Q}{\rho
n^{2}D^{5}},$ (41)

where $T$ and $Q$, respectively, represent the force and the torque exerted by
the fluid on a propeller in the axial direction. The efficiency of a propeller
is defined as

$\eta=\frac{TU_{\infty}}{2\pi nQ}=\frac{J}{2\pi}\cdot\frac{K_{T}}{K_{Q}},$
(42)

where the variables have the same meanings as previously explained.

We first compare the loads predicted by the simulation with those measured in
a previous experiment [42] at various advance ratios. As can be seen from Fig.
8, the present numerical results agree very well with the experimental values
under all working conditions: overloading condition (when $J$ is small), near-
design condition, and underloading condition (when $J$ is large). The maximum
difference is observed on the efficiency curve, which is around $5\%$. The
present simulation predicts the highest efficiency around $J=0.9$, which is
close to the design condition that is $J=0.833$ (see Tab. 1). The experiment,
however, shows an optimum performance slightly above $J=1.0$, which is further
away from the design condition. Meanwhile, in the same figure, we also compare
the present results with the latest (also the best available) results on the
same propeller from a low-order simulation [47] using the commercial software
STAR-CCM+. It is evident that the low-order solver predicts the efficiency
well only in a narrow range of working conditions where $J$ is small. When $J$
increases, the low-order prediction becomes much worse. As will be shown in
the next section that when $J$ increases, the flow vortices become weaker,
which are more vulnerable to numerical dissipations. Since the high-order
method introduces much smaller numerical dissipations, it predicts the loads
very accurately under all conditions. In contrast, the large numerical
dissipations of the low-order method have completely demolished the
predictions when $J$ is large.

Figure 8: The mean blade loads of DTMB 4119 at different working conditions.

Since a propeller works around the design condition most of the time,
therefore a detailed look into this condition is presented in what follows.
The instantaneous $K_{T}$ and $K_{Q}$ of the blades at this condition are
plotted in Fig. 9 from $tU_{\infty}/D=15$ to $45$ for about $36$ revolutions.
The two coefficients are seen to fluctuate at small amplitudes about their
means chaotically due to the turbulent nature of the flow. The mean values
(averaged for about 102 revolutions) are compared with the design [40] and the
experimental [41, 42] values in Tab. 3, where the difference is defined as
(simulation/experiment$-1$)$\times 100\%$. The design was based on potential
flow theories. The two experiments were performed in open water at slightly
different Reynolds numbers. The present Reynolds number is chosen to exactly
match that of [41], but it is obvious that the Reynolds number effects are
small at this level. Overall, we see very good agreements between the
simulation and the experiments as well as the design.

Figure 9: Instantaneous loads on the blades of DTMB 4119 at design condition. |  $Re_{D}~{}(\times 10^{5})$ |  $Re_{c}~{}(\times 10^{5})$ |  $\overline{K}_{T}$ (diff.) |  $\overline{K}_{Q}$ (diff.) |  $\overline{\eta}$
---|---|---|---|---|---
present |  5.59 |  7.3 |  0.1514 (-$2.3\%$) |  0.0274 (-$2.3\%$) |  0.7326
design [40] |  - |  - |  0.1540 (-$1.7\%$) |  0.0290 (-$5.5\%$) |  0.7040
 exp. [41] |  5.59 |  7.3 |  0.1500 ( $0.9\%$) |  0.0285 (-$3.9\%$) |  0.6978
 exp. [42] |  7.66 |  10.0 |  0.1460 ( $3.7\%$) |  0.0280 (-$2.1\%$) |  0.6913
Table 3: Mean loads on the blades of DTMB 4119 at design condition.

The simulation allows us to study the loads on different parts of the
propeller. We have summarized the results in Tab. 4, where “mean” denotes the
time-averaged values, “r.m.s” denotes the root-mean-square of the unsteady
components, and the subscripts “$p$” and “$v$” denote pressure and viscous
contribution, respectively. Before proceeding further, it is worth noting that
the thrust on the fairwater needs to be treated carefully. Unlike the blades
which are two-sided and closed, the fairwater is “one-sided”, i.e., it has no
direct upstream counterpart to balance the pressure force on its outer
surface. In real applications, this force will be partially compensated by the
force on an upstream surface with the same cross-sectional area (for example,
part of the hull), resulting in a much smaller net force. The flow in the
upstream region is usually at reduced speed (or even at stagnation) with
higher pressure than that of the freestream. Thus, an underestimated pressure
force on this upstream surface is $p_{\infty}(\pi D_{h}^{2}/4)$, where $D_{h}$
is the hub diameter. We then subtract this force from that on the fairwater to
approximate the net thrust on the fairwater.

From Tab. 4, we see that only the blades experience a thrust (positive
$K_{T}$), while the hub and the fairwater experience drags (negative $K_{T}$).
On the other hand, all three parts experience positive torques. For the
blades, pressure contribution dominates the loads; the r.m.s. values are at
least three magnitudes smaller than the means, suggesting quasi-steady loads.
The $K_{T}$ and $K_{Q}$ of the hub are four magnitudes smaller than those of
the blades, which indicates that the hub has negligible contribution/effect to
the overall performance of the propeller. Furthermore, viscous contribution
dominates the hub loads, which is consistent with the fact that the hub has no
projection in the axial direction and thus pressure has no way to contribute.
In fact, pressure should ideally have zero contribution to the hub loads, and
the present very small pressure contribution implies very small geometric
imperfection of the hub, which obviously has benefited from the high-order
curved representation of the geometry. For the fairwater, the drag is about
$2.1\%$ of the thrust on the blades, and pressure dominates; the torque is
negligibly small, and viscosity dominates due to the geometric symmetry.
Overall, we conclude that the thrust and the torque of the blades, and the
drag of the fairwater are the main factors that affect the propeller’s
performance.

|  |  $K_{T}$ |  $K_{T,p}$ |  $K_{T,v}$ |  $K_{Q}$ |  $K_{Q,p}$ |  $K_{Q,v}$
---|---|---|---|---|---|---|---
blades |  mean |  0.1514 |  0.1519 |  -4.9E-4 |  0.0274 |  0.02722 |  2.4E-4
|  r.m.s. |  3.6E-4 |  3.6E-4 |  5.6E-7 |  7.0E-5 |  7.0E-52 |  2.6E-8
hub |  mean |  -8.5E-5 |  -6.7E-7 |  -8.5E-5 |  3.9E-6 |  3.1E-82 |  3.9E-6
|  r.m.s. |  2.1E-7 |  4.5E-8 |  2.0E-7 |  2.2E-8 |  5.1E-92 |  2.2E-8
fairwater |  mean |  -3.2E-3 |  -3.2E-3 |  -2.4E-5 |  8.8E-7 |  6.9E-10 |  8.8E-7
|  r.m.s. |  1.2E-4 |  1.2E-4 |  4.8E-7 |  9.9E-8 |  9.3E-82 |  3.2E-8
Table 4: Loads on different parts of DTMB 4119 with an ellipsoidal fairwater
at design condition.

The time and frequency scales of the major loads can be identified through the
autocorrelation and the power spectral density (PSD) curves in Fig. 10. It is
worth mentioning that the torque and the thrust of a blade have almost
identical characteristics, the curves for the torque are therefore not
repeated here. The autocorrelation is defined as $\rho(\tau)=R(\tau)/R(0)$,
where $R(\tau)=\left<T(t)T(t+\tau)\right>$ is the autocovariance with a time
lag $\tau$, and “$\left<~{}\right>$” denotes ensemble average. We calculate
the PSD using Welch’s method with a $50\%$ overlapping Hanning window and then
average the results over 102 revolutions.

The narrow main lobe of the autocorrelation of the blade in Fig. 10 implies
small integral time scales of the unsteadiness around the blade. It also
agrees with the corresponding high-frequency peak around $fD/U_{\infty}\approx
7.9$ on the PSD curve. In contrast, the autocorrelation of the fairwater has a
much wider main lobe, which indicates much larger time scales of the dominant
unsteadiness around the fairwater. The corresponding PSD curve is rather
broadband and is dominated by very low-frequency components. A comparison of
the two PSD curves reveals that the fairwater experiences more unsteadiness at
very high-frequencies (e.g., $fD/U_{\infty}>20$) than the blade does. These
time and frequency scales are directly related to flow structures which will
be discussed in the next section.

Figure 10: Autocorrelation and PSD of the blade thrust and fairwater drag.

###  4.2 Flow fields

The flow fields of a propeller have very distinct flow structures that are of
crucial importance to the propeller’s performance. This has made studying the
formation, mutual interaction, and stability of these flow structures a
constant research topic for decades. In this section, we report the details of
the flow fields of DTMB 4119, including the vortices, the velocity field, and
the pressure field.

####  4.2.1 Vortices

The vortical structures in a flow field can be well visualized by isosurfaces
of Q-criterion [48]. The Q-criterion (denoted by $Q_{cr}$) is defined as the
second invariant of the velocity gradient tensor, i.e.,
$Q_{cr}=(\Omega_{ij}\Omega_{ij}-S_{ij}S_{ij})/2$, where
$\Omega_{ij}=(u_{i,j}-u_{j,i})/2$ and $S_{ij}=(u_{i,j}+u_{j,i})/2$ are the
antisymmetric and the symmetric component, respectively, of the velocity
gradient tensor.

When the Reynolds number is given, the only parameter that determines a
propeller’s flow field is the advance ratio $J$. Figure 11 shows the
instantaneous vortical structures as $J$ decreases from $1.1$ to $0.4$. Note
that the nondimensional rotational speed is related to the advance ratio as
$\omega^{*}\\!=2\pi/J$. Thus, a decreasing $J$ is equivalent to an increasing
$\omega^{*}$.

We notice two dramatic changes in the flow field as $J$ decreases: the
increase of vortex strength and the occurrence of flow instabilities. At
$J=1.1$ and $1.0$, the vortices are so weak that they are quickly dissipated
by the wake flow. At $J=0.9$ and $0.8$, the vortices become strong enough to
sustain for a long distance in the wake, and a hub vortex is also well
established. In addition, up to this point the flow remains stable. Obvious
instability occurs when $J$ decreases to $0.7$, and the instability is caused
by mutual interactions between two tip vortices around $x/D=4.8$. At $J=0.6$,
the instability is still caused by mutual tip vortex interactions, but the
occurrence moves upstream to $x/D=3.4$. The occurrence further moves upstream
to $x/D=3.2$ and $2.7$, for $J=0.5$ and $0.4$, respectively. However, the
cause of the instability becomes more complicated. At $J=0.5$, it seems the
instability not only comes from the mutual interaction between the tip
vortices, but also the interaction between tip and hub vortices. Finally, at
$J=0.4$, it looks like the trailing edge vortices have become strong enough to
be the leading cause of the instability. It was conjectured in [7] that blade
trailing edge vortices are an important source of flow instabilities. Based on
our observations here, this is only possible when $J$ is small enough (i.e.,
propeller is at very high relative rotational speed) and when blade trailing
edge vortices are strong enough.

Figure 11: Isosurfaces of $Q_{cr}D^{2}/U_{\infty}^{2}=40$ at different working
conditions.

We already saw that the flow fields can be very different at different working
conditions. In the rest of this paper, we focus on the design condition only.
Figure 12 shows an instantaneous view of the flow structures at this
condition. It is seen that the tip vortices are very much equally spaced along
the axial direction, with the distance between two successive vortices being
approximately $0.36D$, which is about one-third the tip pitch (see Tab. 2).
Meanwhile, the surface velocity contours reveal that a tip vortex has lower
streamwise speed on the outer surface, and higher speed on the inner surface.
This indicates that a tip vortex not only revolves helically about the
propeller’s axes, but also about its own core at the same time. The topology
of the root vortices are not very obvious from this instantaneous flow field
due to the many small turbulent structures. These high-frequency small flow
structures are closer to the fairwater than to the blades, and thus have
contributed more to the load unsteadiness of the fairwater than to the blades,
which agrees with our previous observation on the PSD curves in Fig. 10. These
small structures, however, do not dominate the load unsteadiness, which is
likely because of their isotropy that leads to mutual cancellation of the
effects. While the tip vortices break up about $5D$ downstream from the
propeller, the hub vortex stays strong and does not break up even at the
outlet (i.e., $12D$ downstream from the propeller; complete picture not shown
here due to limited space).

Figure 12: Isosurface of instantaneous Q-criterion
$Q_{cr}D^{2}/U_{\infty}^{2}=40$ at design condition.

The phase-averaged isosurfaces of Q-criterion are shown in Fig. 13, and they
reveal the major flow structures, especially the root vortices, more
evidently. The phase-averaged hub vortex is still seen to vary along the axial
direction. In fact, it is the variations of these big structures that dominate
the load unsteadiness of the fairwater. Similarly, the load unsteadiness of
the blades is likely dominated by the unsteadiness of the tip and the trailing
edge vortices.

Figure 13: Isosurface of phase-averaged Q-criterion
$\left<Q_{cr}\right>D^{2}/U_{\infty}^{2}=40$ at design condition.

The FR method is a discontinuous-type of method, and reconstruction must be
performed to make the solutions and fluxes globally continuous. This rule also
applies to the statistics, which are only element-wise continuous unless
reconstructed. However, reconstructing the statistics will impose extra
computational and memory cost to the simulation. For this reason, this process
is not performed in this work, which results in non-smoothness across cell
boundaries as can be seen from Fig. 13, especially in the vicinity of the
propeller where flow changes rapidly. Nevertheless, this simplification should
not alter any of the conclusions here.

The time-averaged flow field is shown in Fig. 14. It is worth noting that
time-averaging is impossible for the sliding region due to the movement of the
propeller, and the flow in this region is thus not shown in the figure. Over
time, the tip vortices form a slightly converging-diverging “duct” in space.
The root vortices, because of the very small instantaneous turbulent
structures, are very difficult to converge in time. Nevertheless, they still
have a tube-like shape over time in space. Unlike the instantaneous and the
phase averaged ones, the time-averaged the hub vortex is very symmetric about
the axis, and almost sees no deviation from the axial direction. This clearly
demonstrates that time-averaging not only helps remove most of the small
unsteadiness, but also the large ones, from the flow field.

Figure 14: Isosurface of time-averaged Q-criterion
$\overline{Q}_{cr}D^{2}/U_{\infty}^{2}=8$ at design condition.

The Q-criterion isosurfaces are able to reveal the most coherent vortical flow
structures. They are, however, inefficient to expose the weak ones like the
trailing edge vortices at the design condition. Additionally, the Q-criterion
cannot reveal the sign of a vortex. For this reason, we have plotted the
streamwise vorticity contours in Figs. 15 and 16 to fill these gaps.

We can clearly see the footprints of the trailing edge vortices (TEVs) in Fig.
15. One end of each TEV connects to a tip vortex, and the other end connects
to the hub or root vortex. As the flow goes downstream, the TEVs tilt more and
more towards downstream, which obviously complies with the wake velocity
distribution. The signs of the vortices reveal that the tip vortices rotate in
the opposite direction to all other vortices as well as the propeller.

Figure 15: Contours of phase-averaged streamwise vorticity in the central
$x$-$y$ plane.

The vorticity contours in Fig. 16 uncover how the vortices develop in the
azimuthal and radial directions as the flow travels downstream. The evolution
of the TEVs is the most prominent: they are elongated and bent in the
azimuthal direction and finally impinge onto the tip vortices. However,
because of flow dissipation and the weak strength of the TEVs, these
interactions do not destabilize the flow. Again, the change fo the TEVs is
closely related to the velocity distribution that will be discussed in the
next section.

Figure 16: Contours of phase-averaged streamwise vorticity in $y$-$z$ planes
at different streamwise locations.

####  4.2.2 Velocity field

We already noticed that the outer and inner surfaces of a tip vortex have
different speeds. The tip vortices thus should be well bounded by velocity
isosurfaces. To confirm this, we plot two isosurfaces of the phase-averaged
streamwise velocity in Fig. 17. It is seen that the isosurfaces of
$\left<u\right>/U_{\infty}=0.9$ and $1.4$ follow the trajectories of the tip
vortices very well. Meanwhile, the increasing gap between the two isosurfaces
also agrees with the decreasing strength of the tip vortices as the flow moves
downstream.

Figure 17: Isosurfaces of phase-averaged streamwise velocity
$\left<u\right>/U_{\infty}=0.9$ (gray) and $1.4$ (yellow).

A closer look of the velocity isosurfaces in the very vicinity of the blades
also reveals the formation of the tip vortices. In Fig. 20, we are looking
towards downstream at the suction side in (a1) and (a2), and towards upstream
at the pressure side in (b1) and (b2). From (a1), it is obvious that each
leading edge (LE) decelerates the incoming flow, resulting in a strand of low-
speed flow along the LE and finally sheds off around the tip. From (a2), each
trailing edge (TE) accelerate the flow and sheds off a strand of high-speed
flow slightly below the tip. When these two strands of flow meet and be
convected downstream, a helical tip-vortex system is generated.

Figure 18: Isosurfaces of $\left<u\right>/U_{\infty}=0.9$ (gray) and $1.4$
(yellow): (a1,a2), suction side; (b1,b2), pressure side.

Figure 19 shows the phase-averaged streamwise velocity contours in the central
$x$-$y$ plane (i.e., $z=0$). It is seen that the blade wake is overall
accelerated, while the fairwater wake is mostly at reduced speed. The flow
immediately downstream of the fairwater has very low speed, suggesting that a
bubble is likely formed in this region. The tip vortices show up as local
velocity min-max pairs along the outer edge of the slipstream. The footprints
of other vortices, such as the TEVs and the hub vortex, are also visible on
the velocity contours.

Figure 19: Contours of phase-averaged streamwise velocity in the central
$x$-$y$ plane.

The Cartesian velocity can be decomposed into three components: streamwise,
radial, and azimuthal, denoted by $u$, $v_{r}$, and $v_{\theta}$,
respectively. Figure 20 shows the time-averaged contours of these components.
We see that the overall slipstream has a converging-diverging shape and is
well contained in the propeller’s swept area, i.e., $r\leq 0.5D$ (note that
$r=|y|$ in the central plane). The radial velocity is small almost everywhere,
except in the very vicinity of the fairwater. The azimuthal speed is large
only in the fairwater wake and is induced by the strong hub vortex. The two
low-speed strips on the azimuthal speed contours around $r/D=0.1$ and $0.4\leq
x/D\leq 1.2$ are footprints of the root vortices.

Figure 20: Contours of time-averaged streamwise, radial and azimuthal velocity
components in the central $x$-$y$ plane.

More detailed velocity profiles are shown in Fig. 21. We see almost no induced
velocity outside the slipstream (i.e., $r/D>0.5$) on all profiles, which
suggests that the propeller introduces very little disturbance to the flow
outside its swept area. The streamwise velocity $\overline{u}$ reaches its
maximum values in the region $0.2<r/D<0.4$, where the flow is accelerated by
more than $30\%$ over $U_{\infty}$. Outside this region, $\overline{u}$
quickly decreases to $U_{\infty}$ around $r/D=0.5$, and also decreases to the
hub-vortex-core speed at $r/D=0$. This velocity distribution leads to the TEV
deformation that is observed in Fig. 15. The streamwise hub-vortex-core speed
is close to zero in the near wake region (e.g., $x/D=0.46$), and then
consistently increases towards downstream. At $x/D=5.0$, it is already
slightly over $U_{\infty}$. The radial velocity $\overline{v}_{r}$ is only
noticeable in the near wake, for example, at $x/D=0.46$ and $0.5$, and then
quickly drops to very small values as the flow travels downstream. The
azimuthal velocity $\overline{v}_{\theta}$ is mostly induced by the hub
vortex. It has almost the same profile at different locations, except in two
small regions: one around $r/D=0.1$ that is affected by the root vortices, and
the other around $r/D=0.42$ that is affected by the tip vortices. The profiles
of $\overline{v}_{\theta}$ are very close to that of a typical Rankine vortex,
and is responsible for the TEV deformation in Fig. 16. The maximum value of
$\overline{v}_{\theta}$ is about $0.62U_{\infty}$ in the near field, and
$0.36U_{\infty}$ in the intermediate wake. These are very large values, and
clearly indicate how strong the hub vortex is.

Figure 21: Time-averaged velocity profiles at different streamwise locations
in the central plane.

As a further validation of the simulation, in Fig. 22 we compare the velocity
profiles at $x/R=0.951$ with a water tunnel measurement from [42]. Overall,
very good agreements are seen between the simulation and the experiment. For
example, in the region $0.3<r/R<1.2$, the maximum difference on $\overline{u}$
is only around $2\%$. However, we do see large discrepancies in the region
that is close to $r/R=0.2$. This is because of the setup difference between
the experiment and the simulation. In the experiment, the propeller shaft is
actually at downstream, making it a stationary wall surface at
$(x/R,r/R)=(0.951,0.2)$. In contrast, for the simulation, the shaft is at
upstream, and it is a flow region at the same location.

Figure 22: Comparison of mean velocity profiles at $x/R=0.951$ with
experiment.

Based on the Kutta-Joukowski theorem, for an inviscid flow the “lift” on a
unit span of a body (such as a propeller blade) is proportional to the
circulation. Following the derivation in [49, 50], the circulation around a
blade section at $r$, denoted by $\Gamma(r)$, is related to the
circumferential speed of the slipstream as

$\Gamma(r)\approx-2\pi rv_{\theta}(r)/Z,$ (43)

where $Z$ is the number of blades of a propeller. Applying the above relation
to the mean flow, we can define the following nondimensional circulation for a
blade,

$G(r)=-\frac{1}{Z}\frac{r}{R}\frac{\overline{v}_{\theta}(r)}{U_{\infty}}\approx\frac{\Gamma(r)}{2\pi
RU_{\infty}}.$ (44)

This variable can be employed to measure the load distribution on a blade.
Moreover, $G$ is conservative for inviscid flows, and should be roughly
conservative for high Reynolds number flows (where viscous effects are small).
Figure 23 shows the circulation profiles at different streamwise locations.
The curves in (a) and (b) start around $r/R=0.2$ because of the presence of
the fairwater at these two locations. An experimental measurement from [42] is
also shown in (c), which agrees well with the simulation result. Overall, we
see that all profiles have very similar shapes and amplitudes, which confirms
that the circulation is indeed roughly conserved. Nevertheless, viscous
effects are still evident since the local narrow peaks are gradually smoothed
out as the flow travels downstream. These curves also signify that the load is
mostly concentrated around the mid-section (i.e., $r/R=0.5$) of each blade.
The large peak around $r/R=0.9$ in the near filed indicates that the propeller
is also heavily loaded around the tips, which is consistent with the strong
tip vortices that we have observed in the flow field.

Figure 23: Flow circulation at different streamwise locations.

####  4.2.3 Pressure field

Figure 24 shows the phase-averaged pressure field in the central $x$-$y$
plane. Comparing with the isosurfaces of Q-criterion in Fig. 13, we see that
the tip vortices show up as local pressure minima along the edge of the
slipstream. Meanwhile, the fairwater wake, especially the hub vortex, is a
very low pressure region, and is responsible for the large drag on the
fairwater (see Tab. 4).

Figure 24: Contours of phase-averaged pressure in the central $x$-$y$ plane.

Of great importance is the pressure distribution around the blades, which
directly affects the thrust and torque on the propeller. To see the pressure
effects on the thrust more clearly, we have plotted in Fig. 25 several slices
in the $x$-$y$ plane through the top blade at different spanwise (i.e., $z$)
locations. From (a) to (h), we are moving from the leading edge to the
trailing edge of the blade (refer to Fig. 5) along the $z$ direction. The
suction side is on the left, and the pressure side is on the right. As
expected, we see that most part of the suction side is in a low pressure
region, and the pressure side is in a high pressure region. When we go from
the leading edge to the trailing edge, the size of the low pressure region
first increases and then decreases, with the maximum size around the mid-span,
i.e., $z=0$. In contrast, the size of the high pressure region first
decreases, and then increases. These pressure distributions apparently suggest
that the thrust load is more concentrated on the trailing portion ($z>0$) of
the blade.

Figure 25: Contours of phase-averaged pressure in different $x$-$y$ planes
through the top blade.

Similarly, the contribution from different parts of a blade to the torque can
be visualized through the pressure distribution in different $y$-$z$ planes
through the blades as shown in Fig. 26. Again, from (a) to (h) we are moving
from the leading edge to the trailing edge (refer to Fig. 5), but along the
$x$ direction this time. Taking the top blade for example, the suction side is
on the left, and the pressure side is on the right. From this perspective, the
suction side is almost always in a low pressure region. In (a)-(d), a large
portion of the pressure side actually has low surface pressure (although

Figure 26: Contours of phase-averaged pressure in different $y$-$z$ planes
through the blades.

there is a large high-pressure “bubble”, but it is detached from the blade
surface). In contrast, the pressure side has increased pressure in (e)-(h).
This pressure distribution results in higher torque on the trailing portion of
the blade. We also notice that the tips of the cross-sections have the largest
pressure difference most of the time. Considering the fact that the tips also
has the largest arms in the cross-sections, the torque is therefore also very
heavily loaded around the edge of each blade.

The time-averaged pressure field in the central $x$-$y$ plane is shown in Fig.
27, and a series of pressure profiles are given in Fig. 28. From the contours,
it is seen that on average the propeller generates an obvious high-pressure
region in the near field $x/D<1.0$ (inside and outside of the slipstream).
Other than this region, the propeller’s effects on the pressure field are
mostly contained within the slipstream. The hub vortex represents the pressure
minima of the whole flow field. From the profiles, we notice three local
pressure minima around $r/D=0$, $r/D\approx 0.12$ (only in the very near
field), and $r/D\approx 0.42$. They actually correspond to the hub vortex, the
root vortices, and the tip vortices, respectively. The pressure recovery in
the blade wake ($0.1<r/D<0.5$) is evident as the flow goes downstream. In
contrast, we do not see consistent pressure recovery in the fairwater wake
($r/D<0.1$) due to the strong hub vortex.

Figure 27: Contours of time-averaged pressure in the central $x$-$y$ plane.
Figure 28: Profiles of time-averaged pressure at different streamwise
locations in the central $x$-$y$ plane.

###  4.3 Fairwater effects

The fairwater is usually not considered in propeller design. A user has the
freedom to choose a fairwater based on their preference or the availability of
parts. The quantitative effects of fairwater shape have rarely been studied.
In this section, we briefly study two more fairwater shapes: cylindrical and
hemispherical, to compare with the ellipsoidal one from the previous sections.
For a fair comparison, we require the three fairwaters to have the same
surface area so that they contact with the same amount of fluids. This means
that the fairwaters will have different lengths. The 1:2 ellipsoidal fairwater
has a length of $0.2D$ as shown in Fig. 5. The geometries and sizes of the
other two are shown in Fig. 29. As marked in the figure, the cylindrical
fairwater has a length of $0.121D$, and the hemispherical one has a length of
$0.171D$ (a hemisphere of radius $0.1D$ sitting on top of a cylinder whose
height is $0.071D$). Overall, the shape is more elongated (streamlined) as the
shape changes from cylindrical to hemispherical and then to ellipsoidal.

Figure 29: DTMB 4119 with cylindrical fairwater (left) and hemispherical
fairwater (right).

The loads for the above two configurations are summarized in Tabs. 5 and 6,
respectively. The overall loads on the hubs and the torques on the fairwaters
are once again negligibly small. Comparing with the ellipsoidal configuration
(see Tab. 4), we notice that the blades in both configurations here have
smaller thrust and torque. However, the efficiencies of the blades (excluding
fairwater contributions) in all three configurations stay almost unaffected as
shown in the first row of Tab. 7, where the efficiency differences are around
$0.3\%$. The drags on the fairwaters, on the other hand, are dramatically
different for the three configurations. The cylindrical fairwater has the
largest drag, followed by the hemispherical one, and then the ellipsoidal one.
The overall efficiencies (with fairwater contributions included) are also
summarized in Tab. 7. We see that the cylindrical, the hemispherical, and the
ellipsoidal fairwaters reduce the propeller’s overall efficiency by $3.3\%$,
$2.9\%$, and $2.1\%$, respectively. We need to emphasize that these numbers
are not small for a propulsion system, and we also need to repeat that these
numbers are underestimated based on the assumptions that we made in Sec. 4.1.
Since pressure contribution dominates the drag on a fairwater, it is worth
looking into the pressure fields to see how the fairwaters affect the pressure
distributions.

|  |  $K_{T}$ |  $K_{T,p}$ |  $K_{T,v}$ |  $K_{Q}$ |  $K_{Q,p}$ |  $K_{Q,v}$
---|---|---|---|---|---|---|---
blades |  mean |  0.1503 |  0.1507 |  -4.8E-4 |  0.0271 |  0.02682 |  2.4E-4
|  r.m.s. |  1.4E-4 |  1.4E-4 |  1.1E-7 |  2.8E-5 |  2.8E-52 |  5.8E-9
hub |  mean |  -9.8E-5 |  -8.5E-7 |  -9.7E-5 |  4.0E-6 |  2.9E-82 |  3.9E-6
|  r.m.s. |  1.1E-7 |  3.9E-8 |  9.4E-8 |  1.1E-8 |  2.2E-92 |  1.0E-8
fairwater |  mean |  -4.9E-3 |  -4.9E-3 |  -2.9E-5 |  1.5E-6 |  1.3E-11 |  1.5E-6
|  r.m.s. |  6.0E-5 |  6.0E-5 |  1.2E-7 |  1.9E-8 |  7.9E-12 |  1.9E-8
Table 5: Loads on different parts of DTMB 4119 with a cylindrical fairwater. |  |  $K_{T}$ |  $K_{T,p}$ |  $K_{T,v}$ |  $K_{Q}$ |  $K_{Q,p}$ |  $K_{Q,v}$
---|---|---|---|---|---|---|---
blades |  mean |  0.1502 |  0.1507 |  -4.8E-4 |  0.0271 |  0.02692 |  2.4E-4
|  r.m.s. |  3.2E-4 |  3.2E-4 |  5.6E-7 |  6.2E-5 |  6.2E-52 |  3.0E-8
hub |  mean |  -9.1E-5 |  -9.2E-7 |  -9.0E-5 |  3.9E-6 |  2.9E-82 |  3.9E-6
|  r.m.s. |  1.8E-7 |  5.2E-8 |  1.6E-7 |  2.9E-8 |  3.4E-92 |  2.9E-8
fairwater |  mean |  -4.4E-3 |  -4.4E-3 |  -3.9E-5 |  1.3E-6 |  3.5E-11 |  1.3E-6
|  r.m.s. |  1.6E-4 |  1.6E-4 |  2.8E-7 |  4.0E-8 |  5.8E-10 |  4.0E-8
Table 6: Loads on different parts of DTMB 4119 with a hemispherical fairwater. |  cylindrical |  hemispherical |  ellipsoidal
---|---|---|---
blades efficiency (excl. FW) |  0.7353 |  0.7348 |  0.7326
overall efficiency (incl. FW) |  0.7113 |  0.7133 |  0.7171
efficiency loss (from FW) |  -3.3% |  -2.9% |  -2.1%
Table 7: Blade efficiency and overall efficiency for different fairwater (FW)
configurations.

The time-averaged pressure fields for the above two configurations are shown
in Fig. 30. Comparing with the ellipsoidal configuration in Fig. 27, we see
that the three configurations overall have very similar pressure distribution
in the wake, which explains why the blade efficiencies are not affected much.
The differences are mostly limited to the region immediately downstream of the
fairwater. The cylindrical fairwater generates a very large low-pressure
region at its end. The hemispherical fairwater has a low-pressure region of
similar size to that of the ellipsoidal one. However, at the junction of the
hemispherical fairwater and the hub, there is another small but very low-
pressure region due to the geometric change, which leads to an overall larger
drag for this configuration than the ellipsoidal one.

Figure 30: Time-averaged pressure in the central $x$-$y$ plane of DTMB 4119
with cylindrical and hemispherical fairwaters.

##  5 Summary

The first high-order eddy-resolving simulation of a marine propeller has been
successfully performed in this work using a recently developed sliding-mesh
method. This method combines the flux reconstruction framework and a new
dynamic curved mortar approach to deal with the complex rotating geometry of a
propeller without sacrificing the high-order accuracy at all. Even on a very
coarse mesh with less than one-fourth million cells, the method predicts the
propeller loads very accurately under a wide range of working conditions, and
also captures the flow structures with a lot of details. Moreover, this method
allows both phase and time averaging on the same set of grid, and thus can
provide more information about a flow field.

Through visualization of vortical flow structures, it is revealed that when
the advance coefficient $J$ decreases, the strengths of the major vortices
grow and flow instability gradually develops. The instability first comes from
tip-tip vortex interaction, then tip-tip as well as tip-hub vortex
interactions, and finally the trailing-edge vortices become strong enough and
start playing an important rule when $J$ is sufficiently small. At the design
condition, the sources of each tip vortex are identified through velocity
isosurfaces to be a strand of decelerated flow from the leading edge and a
strand of accelerated flow from the trailing edge of each blade.

A comparison between the present high-order simulation and a previous low-
order one on the same propeller has clearly demonstrated the low-dissipation
advantage of the high-order method, which has allowed accurate prediction of
the loads under all working conditions. In contrast, the high-dissipation of
the low-order method completely failed the mission for large $J$ that
generates weak flow vortices. Detailed load analysis at the design condition
has revealed that the major loads are the blade thrust and torque as well as
the fairwater drag, and pressure contribution dominates these loads. The
pressure field and the circulation distribution show that the blade loads
concentrate more on the trailing portion as well as the radial mid-section of
each blade.

By studying three fairwaters of different shapes, it is found that these
fairwaters do not have obvious effects on the blade performance in the present
setups. They, however, do dramatically change the pressure distribution on
their surfaces, resulting in different induced drags and different performance
degradation to the overall propulsion system. More specifically, we see an
efficiency loss of at least $3.3\%$, $2.9\%$, and $2.1\%$, from the
cylindrical, the hemispherical, and the ellipsoidal fairwater, respectively.
It remains to be investigated whether there is an optimum fairwater shape that
can minimize the efficiency loss.

## Acknowledgment

This work was supported by: the Office of Naval Research (Grant
N00014-18-1-2265), the George Washington University, and Clarkson University.
The computing resources were provided by the Navy DoD Supercomputing Resource
Center through the High Performance Computing Modernization Program. The
authors would like to thank Dr. Thad Michael from the Naval Surface Warfare
Center Carderock Division for providing the propeller geometry and for helpful
discussions. Thanks are also due to Dr. Zihua Qiu for helping on part of the
mesh generation.

## Dedication

This paper is dedicated to Prof. Antony Jameson for celebrating his 85th
birthday and his immense and continuing contributions to the many aspects of
computational fluid dynamics.

## References

  * Greeley and Kerwin [1982] D. S. Greeley, J. E. Kerwin, Numerical methods for propeller design and analysis in steady flow, Transactions-Society of Naval Architects and Marine Engineers 90 (1982) 415–453.
  * Muscari et al. [2013] R. Muscari, A. Di Mascio, R. Verzicco, Modeling of vortex dynamics in the wake of a marine propeller, Computers & Fluids 73 (2013) 65–79.
  * Di Mascio et al. [2014] A. Di Mascio, R. Muscari, G. Dubbioso, On the wake dynamics of a propeller operating in drift, Journal of Fluid Mechanics 754 (2014) 263–307.
  * Verma et al. [2012] A. Verma, H. Jang, K. Mahesh, The effect of an upstream hull on a propeller in reverse rotation, Journal of Fluid Mechanics 704 (2012) 61.
  * Jang and Mahesh [2013] H. Jang, K. Mahesh, Large eddy simulation of flow around a reverse rotating propeller, Journal of Fluid Mechanics 729 (2013) 151–179.
  * Balaras et al. [2015] E. Balaras, S. Schroeder, A. Posa, Large-eddy simulations of submarine propellers, Journal of Ship Research 59 (4) (2015) 227–237.
  * Kumar and Mahesh [2017] P. Kumar, K. Mahesh, Large eddy simulation of propeller wake instabilities, Journal of Fluid Mechanics 814 (2017) 361–396.
  * Reed and Hill [1973] W. H. Reed, T. R. Hill, Triangular mesh methods for the neutron transport equation, Tech. Rep. LA-UR-73-479; CONF-730414-2, Los Alamos Scientific Laboratory, 1973\.
  * Cockburn et al. [2011] B. Cockburn, G. E. Karniadakis, C.-W. Shu (Eds.), Discontinuous Galerkin Methods: Theory, Computation and Applications, vol. 11 of _Lecture Notes in Computational Science and Engineering_ , Springer, New York, 2011\.
  * Patera [1984] A. T. Patera, A spectral element method for fluid dynamics: laminar flow in a channel expansion, Journal of Computational Physics 54 (3) (1984) 468–488.
  * Karniadakis and Sherwin [2005] G. E. Karniadakis, S. J. Sherwin, Spectral/hp Element Methods for Computational Fluid Dynamics, Oxford University Press, Oxford, 2 edn., 2005.
  * Kopriva [2009] D. A. Kopriva, Implementing spectral methods for partial differential equations: Algorithms for scientists and engineers, Springer Science & Business Media, 2009\.
  * Wang [2002] Z. J. Wang, Spectral (finite) volume method for conservation laws on unstructured grids: basic formulation, Journal of Computational Physics 178 (2002) 210–251.
  * Wang and Liu [2002] Z. J. Wang, Y. Liu, Spectral (finite) volume method for conservation laws on unstructured grids: II. Extension to two-dimensional scalar equation, Journal of Computational Physics 179 (2) (2002) 665–697.
  * Kopriva and Kolias [1996] D. A. Kopriva, J. H. Kolias, A conservative staggered-grid Chebyshev multidomain method for compressible flows, Journal of Computational Physics 125 (1996) 244–261.
  * Kopriva [1998] D. A. Kopriva, A staggered-grid multidomain spectral method for the compressible Navier-Stokes equations, Journal of Computational Physics 143 (1998) 125–158.
  * Liu et al. [2006] Y. Liu, M. Vinokur, Z. J. Wang, Spectral difference method for unstructured grids I: Basic formulation, Journal of Computational Physics 216 (2006) 780–801.
  * Wang et al. [2007] Z. J. Wang, Y. Liu, G. May, A. Jameson, Spectral difference method for unstructured grids II: Extension to the Euler equations, Journal of Scientific Computing 32 (2007) 45–71.
  * Jameson [2010] A. Jameson, A proof of the stability of the spectral difference method for all orders of accuracy, Journal of Scientific Computing 45 (2010) 348–358.
  * Balan et al. [2012] A. Balan, G. May, J. Schöberl, A stable high-order spectral difference method for hyperbolic conservation laws on triangular elements, Journal of Computational Physics 231 (5) (2012) 2359–2375.
  * Huynh [2007] H. T. Huynh, A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods, AIAA paper 2007-4079, 2007.
  * Huynh [2009] H. T. Huynh, A reconstruction approach to high-order schemes including discontinuous Galerkin for diffusion, AIAA paper 2009-403, 2009.
  * Wang and Gao [2009] Z. J. Wang, H. Gao, A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids, Journal of Computational Physics 228 (2009) 8161–8186.
  * Jameson et al. [2012] A. Jameson, P. E. Vincent, P. Castonguay, On the non-linear stability of flux reconstruction schemes, Journal of Scientific Computing 50 (2) (2012) 434–445.
  * Wang and Huynh [2016] Z. J. Wang, H. T. Huynh, A review of flux reconstruction or correction procedure via reconstruction method for the Navier-Stokes equations, Mechanical Engineering Reviews 3 (1) (2016) 15–00475.
  * Zhang and Liang [2015a] B. Zhang, C. Liang, A simple, efficient, high-order accurate sliding-mesh interface approach to FR/CPR method on coupled rotating and stationary domains, AIAA paper 2015-1742, 2015a.
  * Zhang and Liang [2015b] B. Zhang, C. Liang, A simple, efficient, and high-order accurate curved sliding-mesh interface approach to spectral difference method on coupled rotating and stationary domains, Journal of Computational Physics 295 (2015b) 147–160.
  * Zhang et al. [2016] B. Zhang, C. Liang, J. Yang, Y. Rong, A 2D parallel high-order sliding and deforming spectral difference method, Computers & Fluids 139 (2016) 184–196.
  * Zhang and Liang [2016] B. Zhang, C. Liang, A high-order sliding-mesh spectral difference solver for simulating unsteady flows around rotating objects, in: 31st Symposium on Naval Hydrodynamics, Monterey, CA, 2016.
  * Zhang et al. [2018] B. Zhang, Z. Qiu, C. Liang, A flux reconstruction method with nonuniform sliding-mesh interfaces for simulating rotating flows, AIAA paper 2018-1094, 2018.
  * Zhang and Liang [2020] B. Zhang, C. Liang, A conservative high-order method utilizing dynamic transfinite mortar elements for flow simulation on curved sliding meshes, arXiv preprint arXiv:2009.04400, 2020.
  * Parsani et al. [2010] M. Parsani, G. Ghorbaniasl, C. Lacor, E. Turkel, An implicit high-order spectral difference approach for large eddy simulation, Journal of Computational Physics 229 (2010) 5373–5393.
  * Uranga et al. [2011] A. Uranga, P.-O. Persson, M. Drela, J. Peraire, Implicit large eddy simulation of transition to turbulence at low Reynolds numbers using a discontinuous Galerkin method, International Journal for Numerical Methods in Engineering 87 (1-5) (2011) 232–261.
  * Beck et al. [2014] A. D. Beck, T. Bolemann, D. Flad, H. Frank, G. J. Gassner, F. Hindenlang, C.-D. Munz, High-order discontinuous Galerkin spectral element methods for transitional and turbulent flow simulations, International Journal for Numerical Methods in Fluids 76 (8) (2014) 522–548.
  * Bull and Jameson [2015] J. R. Bull, A. Jameson, Simulation of the Taylor-Green vortex using high-order flux reconstruction schemes, AIAA Journal 53 (9) (2015) 2750–2761.
  * Vermeire et al. [2016] B. C. Vermeire, S. Nadarajah, P. G. Tucker, Implicit large eddy simulation using the high-order correction procedure via reconstruction scheme, International Journal for Numerical Methods in Fluids 82 (5) (2016) 231–260.
  * Thomas and Lombard [1979] P. D. Thomas, C. K. Lombard, Geometric conservation law and its application to flow computations on moving grids, AIAA Journal 17 (1979) 1030–1037.
  * Bathe [2006] K.-J. Bathe, Finite Element Procedures, Klaus-Jurgen Bathe, 2006.
  * Rusanov [1961] V. V. Rusanov, Calculation of interaction of non-steady shock waves with obstacles, Journal of Computational and Mathematical Physics USSR 1 (1961) 267–279.
  * Denny [1968] S. B. Denny, Cavitation and open-water performance tests of a series of propellers designed by lifting-surface methods, Tech. Rep. No. 2878, David W. Taylor Naval Ship Research and Development Center, 1968.
  * Jessup et al. [1984] S. D. Jessup, C. Schott, M. Jeffers, S. Kobayashi, Local Propeller Blade Flows in Uniform and Sheared Onset Flows Using LDV Techniques, Fifteenth Symposium on Naval Hydrodynamics, Hamburg, Germany (1984) 221–235.
  * Jessup [1989] S. D. Jessup, An experimental investigation of viscous aspects of propeller blade flow, Ph.D. thesis, The Catholic University of America, 1989.
  * Wilson [1994] R. E. Wilson, Aerodynamic behavior of wind turbines, in: D. A. Spera (Ed.), Wind Turbine Technology: Fundamental Concepts of Wind Turbine Engineering, ASME Press, 1994.
  * Jameson [1983] A. Jameson, Solution of the Euler equations for two dimensional transonic flow by a multigrid method, Applied Mathematics and Computation 13 (3-4) (1983) 327–355.
  * Spiteri and Ruuth [2002] R. J. Spiteri, S. J. Ruuth, A new class of optimal high-order strong-stability-preserving time discretization methods, SIAM Journal on Numerical Analysis 40 (2002) 469–491.
  * Fischer and Mullen [2001] P. Fischer, J. Mullen, Filter-based stabilization of spectral element methods, Comptes Rendus de l’Academie des Sciences Series I Mathematics 332 (3) (2001) 265–270.
  * Sezen and Kinaci [2019] S. Sezen, O. K. Kinaci, Incompressible flow assumption in hydroacoustic predictions of marine propellers, Ocean Engineering 186 (2019) 106138.
  * Hunt et al. [1988] J. Hunt, A. Wray, P. Moin, Eddies, streams, and convergence zones in turbulent flows, in: Proceedings of CTR Summer Program, Center for Turbulence Research, 193–208, 1988.
  * Kerwin [1982] J. E. Kerwin, Flow Field Computations for Non-Cavitating and Cavitating Propellers, in: 14th ONR Symposium on Naval Hydrodynamics, Ann Arbor, Michigan, 1982\.
  * Wang [1985] M.-H. Wang, Hub effects in propeller design and analysis, Ph.D. thesis, Massachusetts Institute of Technology, 1985.


11institutetext: Skobeltsyn Institute of Nuclear Physics, Lomonosov Moscow
State University, Moscow, Russia

# A Scientific Workflow System for Satellite Data Processing
with Real-Time Monitoring

Minh Duc Nguyen 11<EMAIL_ADDRESS>
###### Abstract

This paper provides a case study on satellite data processing, storage, and
distribution in the space weather domain by introducing the Satellite Data
Downloading System (SDDS). The approach proposed in this paper was evaluated
through real-world scenarios and addresses the challenges related to the
specific field. Although SDDS is used for satellite data processing, it can
potentially be adapted to a wide range of data processing scenarios in other
fields of physics.

##  1 Introduction

The ultimate goal of the research in space weather is to create complex
operational magnetosphere-ionosphere models which allow predicting the
radiation risk for space satellites in different orbits, and also estimating
the occurrence of risks of technological disasters due to magnetic storms and
charged particle precipitation. The models use space-related data sets
acquired by multiple satellites of various purposes and capabilities as well
as ground stations to produce forecasts. Such models require all data sets to
be collected at one place and also to be converted into a unified data format
used by them as the input. Such task is challenging and time-consuming because
data formats, storage methods, datasets and values of registered factors
differ between satellites. The size of data sets can be huge so that the disk
capacity of a personal computer might be insufficient for data storage. There
is a high demand for creating an automatic system for collecting data from
various satellites, processing them, converting them to a unified canonical
format and making them available in a standardized machine-readable format.

Currently, there are several existing solutions to the problem such as SWX
smdc , xv_SAT xv_sat ; pthread_xv_sat , the Herschel Data Processing System
herschel , and NJOY njoy . The Satellite Data Downloading System (SDDS)
presented in this paper has several advantages in comparison with the systems
listed above. Firstly, SDDS does not act as a specific data processing system
that targets a particular data source and has just one workflow. SDDS is a
framework that eases the task of creating such systems that can have different
workflows in data processing. Secondly, SDDS provides a web interface for data
scientists, a RESTful API and a Python library for developers so they can
build their systems upon SDDS. Thirdly, from the first day, SDDS has been
developed as a modular system that can scale easily. Components of SDDS can
run either on the same machine or different machines connected via a local
network or the Internet. Different versions of a component can run
interchangeably. Each component can be replaced with a new one just by editing
the proper configuration file as long as the new component follows the
convention. And finally, SDDS has a built-in real-time monitoring subsystem
that helps operators to identify the scope of an error quickly and prevent it
in the future. To our knowledge, SDDS is the first framework dedicated to
satellite data processing that provides API for developers.

##  2 Overview

SDDS has been designed to implement the whole workflow of a satellite data
processing cycle. A cycle consists of six steps: 1) Connecting to data
sources. A data source is a server where data of a satellite are collected and
can be downloaded; 2) Downloading actual data (so-called raw data). Raw data
of a satellite can be a big one file in which data of all instruments
installed on the satellite are packed or a set of files each of which contains
data of an instrument; 3) Extracting binary (level-0) data of each instrument
from the raw data; 4) Processing level-0 data and producing scientific
(level-1) data that are ready to be loaded into a database; 5) Loading level-1
data to a database; 6) Compressing and moving raw, level-0, and level-1 data
to storage.

SDDS consists of two parts: a backend and a frontend. A detailed view of the
architecture is shown in figure 1. Backend components have been developed
using the Django Framework written in Python. Frontend components have been
developed using the React library created by Facebook, Javascript, HTML, and
CSS.

Figure 1: The architectural overview of SDDS

The backend part consists of the following components:

  * •

A set of satellite controllers, each of which is responsible for a satellite
and implement the whole data processing cycle. Each satellite controller uses
decoders provided by satellite developers to decode level-0 data and produce
level-1 data;

  * •

The Live Monitor subsystem which is responsible for monitoring all activities
of all components in real time. The Live Monitor subsystem records each event
during the data processing cycle and sends it to the RabbitMQ message broker;

  * •

The DBIface interface which is responsible for database interaction. The
DBIface interface supports the following engines: PostgreSQL, MySQL, and
Oracle;

  * •

The RabbitMQ message broker which is responsible for broadcasting all recorded
events (short text messages) to all subscribed clients. A subscribed client is
an instance of the Live Monitor’s web interface that has established a
WebSocket connection to the RabbitMQ message broker. The messages are sent
using the STORM protocol;

  * •

A web server, which could be either Apache or Nginx, to serve RESTful
requests.

On the frontend, there are two base components available currently:

  * •

The Live Monitor’s web interface shows recorded events (or current conditions)
during the data processing cycle;

  * •

The Data Access services allow users to acquire data from the database. Data
can be retrieved either by using the web interface or the RESTful API.

Each component of SDDS is a micro service with a RESTful API. Components can
work together locally on the same computer or can be distributed across
different machines connected via a local network or the Internet. Almost every
component of SDDS has been developed as a library. Hence it gives users of
SDDS the possibility to customize certain components to meet their
requirements. It is also possible to change the behavior of a component by
editing its configuration file which is written in the JSON format.

##  3 Workflow automation

The main component of SDDS is the satellite controller which is purely a
Python class called SDDownloader. Each satellite controller in SDDS is
dedicated for processing data of a certain satellite. Each step of the data
processing cycle is implemented as a function of the SDDownloader class. Each
function, in turn, is a composition of smaller functions (routines). Using the
object oriented programming principles, one can create a customized class
which is inherited from the SDDownloader class and redefine certain functions
to meet their requirements.

Usually, data of a satellite are duplicated in several data sources. In this
case, it is enough to download and process data from the first data source.
But in many cases, data are distributed across data sources, and they need to
be merged after processing. The SDDownloader class can handle all described
scenarios. If a data source is temporarily or permanently unavailable, the
satellite controller will retry and continue downloading data when the data
source becomes available.

Different data sources have different security constraints and thus use
different types of connection. At the moment the following connection types
are supported by SDDS out-of-box: HTTP(S), FTP(S), PPTP, IPsec, HTTP(S)/FTP(S)
over PPTP or IPsec. To support a new connection type one just needs to
implement a function for this purpose and pass it as a parameter to the
function of the SDDownloader class used for establishing connections to data
sources. The SDDownloader class also supports adding additional routes to the
routing table after a connection has been established. This kind of operation
is required mostly when the data source is located in another subnet that
differs from the subnet leased by the VPN gateway.

Currently, the SDDownloader class supports decoders written in Python, PHP,
and Javascript programming languages. C and Fortran programs compiled
statically are also supported. To call a certain decoder, one must describe
how the decoder should be executed in the configuration file. The SDDownloader
class also supports post-processing. The mechanism of executing a post-
processing decoder is the same as the decoder execution. To control how a
decoder runs, the SDDownloader class records the whole output produced by the
decoder and its return code. All the recorded events are then passed to the
RabbitMQ message broker so that the satellite operators could watch the
execution process in real time.

The SDDownloader class uses the meta information from the instrument
configuration files to access and load data into the database using the
DBIface interface. But one can implement a customized way to load data by
passing a program as a parameter to the SDDownloader class. The execution
mechanism is the same as the decoder execution.

The SDDownloader class requires the storage to be mounted to the file system.
By default, the SDDownloader class compresses data before moving them to the
storage. All raw, level-0, and level-1 data are stored. Currently, this
behavior can be changed only by editing the global configuration file of the
satellite.

##  4 Conclusion

A data processing, storage, and distribution system has been created using the
SDDS framework and deployed at the Space Monitoring Data Center of Skobeltsyn
Institue of Nuclear Physics, Lomonosov Moscow State University. The system is
available at https://downloader.sinp.msu.ru. During eighteen months of
operation, the system has been functioning correctly and showing its
effectiveness and stability. Using the API provided by SDDS to add support for
new satellites to the system is easy and requires a minimum of coding. In
particular, adding support for the DSCOVR satellite and three satellites of
the GOES-series has taken 14 hours and 12 hours accordingly. SDDS has been
approved at SINP MSU as the base framework to be used in satellite service
development in the future.

## Acknowledgment

I would like to thank the Space Monitoring Data Center’s team for
collaboration in the development and support of RDB, metadata, and data access
services. I would like to thank Dr. Vladimir Kalegaev for helpful discussions
on satellite data processing and clear problem statements. Also, I would like
to thank Dr. Alexander Kryukov for his support in all aspects. This project is
supported by RSF grant #16-17-00098.

## References

  * (1) V.O. Barinova, O.G. Barinov, S.Y. Bobrovnikov, S.A. Dolenko, V.V. Kalegaev, D. Mottl, E.A. Muraveva, L.R. Mukhametdinova, and I.N. Myagkova, _Proceedings of the 8th Solar System Plasma Conference_ (Space Research Institute of RAS, Moscow, Russia, 2013) p. 132
  * (2) V.A. Egorov, V.O. Ilyin, E.A. Loupian, A.A. Mazurov, A.A. Proshin, and E.V. Flitman, Current problems in remote sensing of the earth from space 1, 1, 431–436 (2004)
  * (3) D.A. Kobets, A.M. Matveev, A.A. Mazurov, and A.A. Proshin, Current problems in remote sensing of the earth from space 12, 1, 145–155 (2015)
  * (4) S. Ott, J. Bakker, J. Brumfitt, N. de Candussio, L. Dwedari, A. M.Heras, S. Leeks, A.P. Marston, J.-J. Mathieu, J. Pizarro, and H. Siddiqui, Astronomical Data Analysis Software and Systems XV ASP Conference Series 351, 516–519 (2006)
  * (5) R. E. MacFarlane, D. W. Muir, R. M Boicourt, A. C. Kahler, and J. L. Conlin, _The NJOY Nuclear Data Processing System_ (https://github.com/njoy/NJOY2016-manual/raw/master/njoy16.pdf, Version 2016)