# High-Order Large-Eddy Simulations of a Ducted Wind Turbine

Chi Ding 111Graduate Student in the Department of Mechanical and Aeronautical
Engineering<EMAIL_ADDRESS>Clarkson University, Potsdam, NY, 13676, USA
Bin Zhang 222Research Assistant Professor in the Department of Mechanical and
Aeronautical Engineering<EMAIL_ADDRESS>Clarkson University, Potsdam,
NY, 13676, USA  Chunlei Liang 333Professor in the Department of Mechanical and
Aeronautical Engineering, and AIAA Associate Fellow<EMAIL_ADDRESS>Clarkson University, Potsdam, NY, 13676, USA  Kenneth D. Visser 444Associate
Professor in the Department of Mechanical and Aeronautical Engineering, AIAA
Senior Member<EMAIL_ADDRESS>Clarkson University, Potsdam, NY, 13676,
USA  and Guangming Yao 555Associate Professor in the Department of
Mathematics<EMAIL_ADDRESS>Clarkson University, Potsdam, NY, 13676, USA

###### Abstract

A high-order flux reconstruction method coupled with a high-order sliding mesh
method is applied to analyze the performance of a ducted wind turbine at a
Reynolds number of $1.25\times 10^{6}$. To investigate the impacts of the
duct, axial flow simulations are also performed for the corresponding open-
rotor turbine. It is shown that the ducted turbine has higher thrust and power
outputs than the open rotor. Vorticity and velocity fields of two
configurations are also visualized and analyzed to show the differences. To
evaluate the effects of yawed flows, simulations are carried out for both
turbines under five different yaw angles $\gamma=0^{\circ},\pm 15^{\circ},\pm
30^{\circ}$. The results reveal that the ducted wind turbine has higher power
outputs than the open counterpart at all tested yaw angles.

##  1 Introduction

Wind energy has gained growing popularity due to its clean and renewable
property. One promising way to harvest wind energy is the ducted wind turbine
(DWT), also known as the diffuser-augmented wind turbine (DAWT). Compared with
conventional wind turbines, a DWT is equipped with a duct in order to increase
power output, and has the potential to exceed the Beltz limit [1]. To design
an efficient DWT, accurate estimates of its loads and the flow fields are of
crucial importance. In addition, it is quite often that a DWT operates under
yawed inflow conditions when it is surrounded by buildings in urban
environment, or is installed with yawed misalignments. There is thus a natural
need to study the loads under yawed conditions.

Many efforts have been made to numerically study DWTs. Fletcher [2] used the
Blade Element Momentum (BEM) method to investigate DWTs. Georgalas et al. [3]
combined the lifting-line approach for the rotors and the lifting-surface
method for the ducts to study DWTs with large tip clearance. Recently, one of
the most popular ways to study DWT is the CFD actuator-disk (CFD-AD) method,
which models the rotor by an actuator disk and simulate the flow fields by CFD
methods. The first research using the CFD-AD method was done by Hansen et al.
[4]. This method was also employed to investigate a turbine with a flanged
diffuser [5]. Venters et al.[6] utilized the CFD-AD method to study the
effects of a duct’s geometrical parameters on the rotor’s performance. With
increasing computational power, it is now possible to perform more accurate
CFD simulations of a complete DWT without simplifications of the geometries.

High-order CFD methods (third and above orders) are more accurate and
efficient than low-order methods[7]. The most popular high-order methods
include the Discontinuous Galerkin (DG) methods[8, 9], the Spectral Element
(SE) method[10, 11], the Spectral Difference (SD) method[12, 13, 14, 15], and
the Flux Reconstruction (FR) method[16, 17]. Among them, the SD and the FR
methods are based on the differential form governing equations, and are among
the most efficient high-order methods. The FR method has an additional
advantage that this scheme provides a unifying framework that many existing
high-order schemes, e.g. DG and SD schemes, and even some new schemes can be
recovered by choosing different correction functions[16]. The FR method has
seen a lot of progress recently. Jameson et al.[18] performed a nonlinear
stability analysis for the FR method. Many researchers have successfully
extended FR schemes to mixed elements [19, 20, 21, 22, 23]. More recent
developments about the FR scheme can be found in a review paper[24].

To deal with flows over rotating objects, Zhang et al. [25, 26] introduced the
curved dynamic mortar concept and applied it to the development of high-order
sliding-mesh SD and FR methods. These methods were later extended to sliding-
deforming meshes [27] and to deal with 3D geometries [28]. Zhang et al. [29,
30, 31] further introduced the transfinite mortar concept that has no
geometric error and makes a sliding-mesh method arbitrarily high-order
accurate in space and high-order in time. This method has been applied to
simulate flows around rotating cylinders of different cross-sectional shapes
[32], flapping wings for energy harvesting [33], and more recently the first
high-order eddy-resolving simulation of flow over a marine propeller [34]

This study employs the solver developed by Zhang et al. to study the loads and
flow fields of a wind turbine in a ducted and an open configurations. The
geometries, including the duct and the rotors, are represented accurately
using curved high-order meshes. Both axial flow and yawed flow cases are
simulated in this work. The rest of paper is organized as follows: Section II
introduces the numerical methods including the FR method and the sliding-mesh
method based on the concept of dynamic mortar element. Then the simulation
setup is given in Section III. Following that, we present and analyze the
numerical results in Section IV. Finally, Section V concludes this paper.

##  2 Numerical Methods

###  2.1 The Governing Equations

####  2.1.1 The Physical Equations

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

$\partialderivative{\mathbf{Q}}{t}+\partialderivative{\mathbf{F}}{x}+\partialderivative{\mathbf{G}}{y}+\partialderivative{\mathbf{H}}{z}=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]^{\top},$ (2) $\displaystyle\mathbf{F}$
$\displaystyle=\mathbf{F}_{\text{inv}}(\mathbf{Q})+\mathbf{F}_{\text{vis}}(\mathbf{Q},\gradient{\mathbf{Q}}),$
(3) $\displaystyle\mathbf{G}$
$\displaystyle=\mathbf{G}_{\text{inv}}(\mathbf{Q})+\mathbf{G}_{\text{vis}}(\mathbf{Q},\gradient{\mathbf{Q}}),$
(4) $\displaystyle\mathbf{H}$
$\displaystyle=\mathbf{H}_{\text{inv}}(\mathbf{Q})+\mathbf{H}_{\text{vis}}(\mathbf{Q},\gradient{\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}}=\begin{bmatrix}\rho u\\\\[2.84526pt] \rho
u^{2}+p\\\\[2.84526pt] \rho uv\\\\[2.84526pt] \rho uw\\\\[2.84526pt]
u(E+p)\end{bmatrix},\quad\mathbf{G}_{\text{inv}}=\begin{bmatrix}\rho
v\\\\[2.84526pt] \rho uv\\\\[2.84526pt] \rho v^{2}+p\\\\[2.84526pt] \rho
vw\\\\[2.84526pt]
v(E+p)\end{bmatrix},\quad\mathbf{H}_{\text{inv}}=\begin{bmatrix}\rho
w\\\\[2.84526pt] \rho uw\\\\[2.84526pt] \rho vw\\\\[2.84526pt] \rho
w^{2}+p\\\\[2.84526pt] w(E+p)\end{bmatrix}.$ (6)

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

$\displaystyle\mathbf{F}_{\text{vis}}=-\begin{bmatrix}0&\tau_{xx}&\tau_{yx}&\tau_{zx}&u\tau_{xx}+v\tau_{yx}+w\tau_{zx}+\kappa
T_{x}\end{bmatrix}^{\intercal},$ (7)
$\displaystyle\mathbf{G}_{\text{vis}}=-\begin{bmatrix}0&\tau_{xy}&\tau_{yy}&\tau_{zy}&u\tau_{xy}+v\tau_{yy}+w\tau_{zy}+\kappa
T_{y}\end{bmatrix}^{\intercal},$ (8)
$\displaystyle\mathbf{H}_{\text{inv}}=-\begin{bmatrix}0&\tau_{xz}&\tau_{yz}&\tau_{zz}&u\tau_{xz}+v\tau_{yz}+w\tau_{zz}+\kappa
T_{z}\end{bmatrix}^{\intercal},$ (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 RT$, where $R$ is the gas constant.

####  2.1.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,

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

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

$\begin{bmatrix}\widetilde{\mathbf{Q}}\\\\[2.84526pt]
\widetilde{\mathbf{F}}\\\\[2.84526pt] \widetilde{\mathbf{G}}\\\\[2.84526pt]
\widetilde{\mathbf{H}}\end{bmatrix}=\absolutevalue{\mathcal{J}}\mathcal{J}^{-1}\begin{bmatrix}{\mathbf{Q}}\vphantom{\widetilde{\mathbf{Q}}}\\\\[2.84526pt]
{\mathbf{F}}\vphantom{\widetilde{\mathbf{F}}}\\\\[2.84526pt]
{\mathbf{G}}\vphantom{\widetilde{\mathbf{G}}}\\\\[2.84526pt]
{\mathbf{H}}\vphantom{\widetilde{\mathbf{H}}}\end{bmatrix},$ (11)

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

$\absolutevalue{\mathcal{J}}=\absolutevalue{\partialderivative{(t,x,y,z)}{(\tau,\xi,\eta,\zeta)}}=\begin{vmatrix}1&0&0&0\\\
x_{\tau}&x_{\xi}&x_{\eta}&x_{\zeta}\\\ y_{\tau}&y_{\xi}&y_{\eta}&y_{\zeta}\\\
z_{\tau}&z_{\xi}&z_{\eta}&z_{\zeta}\end{vmatrix},$ (12)
$\mathcal{J}^{-1}=\partialderivative{(\tau,\xi,\eta,\zeta)}{(t,x,y,z)}=\begin{vmatrix}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{vmatrix}.$ (13)

Besides the flow equations, the Geometric Conservation Law(GCL) [35] 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 [30,
31, 34]

###  2.2 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,

$\begin{bmatrix}x\\\\[2.84526pt] y\\\\[2.84526pt]
z\end{bmatrix}=\sum^{K}_{i=1}M_{i}(\xi,\eta,\zeta)\begin{bmatrix}x_{i}(t)\\\\[2.84526pt]
y_{i}(t)\\\\[2.84526pt] z_{i}(t)\end{bmatrix},$ (14)

where $K$ is the number of nodes of a hexahedral element, $M_{i}$ (detailed
expressions can be found in [36]) 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-\eta$ plane for a fourth-order FR method.

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

$h_{i}(X)=\prod^{N}_{s=1,s\neq i}(\frac{X-X_{s}}{X_{i}-X_{s}}).$ (15)

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)=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{N}\widetilde{\mathbf{Q}}_{ijk}h_{i}(\xi)h_{j}(\eta)h_{k}(\zeta),$
(16)
$\displaystyle\widetilde{\mathbf{F}}(\xi,\eta,\zeta)=\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{k=1}^{N}\widetilde{\mathbf{F}}_{ijk}h_{i}(\xi)h_{j}(\eta)h_{k}(\zeta),$
(17)

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 solver[37]; 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 polynomial is

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

where $\widetilde{\mathbf{F}}$ is form (17);
$\widetilde{\mathbf{F}}^{\text{com}}$ represents the common flux on a cell
boundary; $g_{L}$ and $g_{R}$ are the left and right correction functions, and
are required to satisfy

$\displaystyle g_{\text{\tiny L}}(0)=1,\quad g_{\text{\tiny L}}(1)=0,$ (19)
$\displaystyle g_{\text{\tiny R}}(0)=0,\quad g_{\text{\tiny R}}(1)=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)$
(20)

i.e., the reconstructed fluxes still take the common values at cell
interfaces. In this work, we have employed the $g_{DG}$ correction
function[16]. The other two fluxes are reconstructed in the same way. Finally,
the governing equations can be written in the following residual form,

$\evaluated{\partialderivative{\widetilde{Q}}{\tau}}_{ijk}=-\left[\partialderivative{\widehat{\mathbf{F}}}{\xi}+\partialderivative{\widehat{\mathbf{G}}}{\eta}+\partialderivative{\widehat{\mathbf{H}}}{\zeta}\right]_{ijk}=\mathbf{R}_{ijk},\quad
i,j,k=1,2,\cdots,N,$ (21)

where $\mathbf{R}_{ijk}$ is the residual at the $(i,j,k)$-th SP. This system
can be time marched explicitly or implicitly. In the present work, we choose a
four-stage three-order explicit Runge-Kutta method. More details can be found
in [38, 39] .

###  2.3 A Sliding Interface 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 [25, 26, 27, 28, 29, 30, 31, 34].

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 interfaces 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},\quad\eta=\eta^{\prime}$ (22)

where $0\leq\xi,\eta,\xi,\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)=\sum_{i=1}^{N}\sum_{j=1}^{N}\Phi^{\Omega}_{ij}h_{i}(\xi)h_{j}(\eta),$
(23)
$\displaystyle\Phi^{\Xi,L}(\xi^{\prime},\eta^{\prime})=\sum_{i=1}^{N}\sum_{j=1}^{N}\Phi^{\Xi,L}_{ij}h_{i}(\xi^{\prime})h_{j}(\eta^{\prime}),$
(24)

where $\phi^{\Omega}_{ij}$ and $\phi^{\Xi,L}_{ij}$ are the discrete values at
the $(i,j)$-th SP on $\Omega$ and the left side of $\Xi$, respectively. The
$(\phi^{\Xi,L}_{ij})^{,}$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})\mathrm{d}\xi^{\prime}\mathrm{d}\eta^{\prime}=0,\quad\forall\alpha,\beta=1,2,...,N.$
(25)

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

$\int_{0}^{1}(\Phi^{\Xi,L}(\xi^{\prime},X_{j})-\Phi^{\Omega}(\xi,X_{j}))h_{\alpha}(\xi^{\prime})\mathrm{d}\xi^{\prime}=0,\quad\forall\alpha=1,2,...,N.$
(26)

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^{\Xi,L}_{1:N,j}$. Repeating this process for every $j$, we will obtain
every $\phi^{\Xi,L}_{ij}$. The values on the right side of a mortar, i.e.,the
$(\phi^{\Xi,R}_{ij})^{,}$, can be obtained in the same way.

(a)

(b)

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 can compute a common value on the mortar, e.g., through
averaging or Riemann solver. Let us denote this common values as $\Phi^{\Xi}$.
We then project this common variable back to a cell face from mortars as shown
in Fig. 4(b). And the projection is

$\displaystyle\int_{0}^{o^{2}}$
$\displaystyle\int_{0}^{1}(\Phi^{\Omega}(\xi,\eta)-\Phi^{\Xi_{1}}(\xi^{\prime},\eta^{\prime}))h_{\alpha}(\xi)h_{\beta}(\eta)\mathrm{d}\xi\mathrm{d}\eta$
(27)
$\displaystyle+\int_{o^{2}}^{1}\int_{0}^{1}(\Phi^{\Omega}(\xi,\eta)-\Phi^{\Xi_{2}}(\xi^{\prime},\eta^{\prime}))h_{\alpha}(\xi)h_{\beta}(\eta)\mathrm{d}\xi\mathrm{d}\eta=0,\quad\forall\alpha,\beta=1,2,...,N.$

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)\mathrm{d}\xi+\int_{o^{2}}^{1}(\Phi^{\Omega}(\xi,X_{j})-\Phi^{\Xi_{2}}(\xi^{\prime},X_{j}))h_{\alpha}(\xi)\mathrm{d}\xi=0,~{}\forall\alpha=1,2,...,N.$
(28)

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

##  3 Simulation Setup

###  3.1 Geometry

The ducted wind turbine was designed by Prof. Visser at Clarkson University.
Fig. 5 illustrates the duct and rotor in both the front and lateral views. As
shown in Fig. 5, the inlet radius of the duct is 0.456m, and the exit radius
is 1.832m. For the duct section, the width is 0.612m, and the angle of attack
is $25^{\circ}$. The rotor has three 1.5m long blades, and the cylindrical hub
is closed by two hemispherical ends of which the diameter is 0.456m. The rotor
in the open configuration has the same size as that of the ducted
configuration.

Figure 5: The ducted wind turbine geometry: left, front view; right, lateral
view.

###  3.2 Computational Mesh

In this study, unstructured hexahedral meshes are used for the simulation.
Each computation mesh consists of a rotating part and a stationary part. The
rotating part includes the rotor. Fig. 6 shows the surface mesh of the duct,
the rotor, and the front sliding interface, for the ducted case. For this
case, the rotating subdomain has two sliding interfaces. One interface is
visible in Fig. 6(b), and the other interface is behind the duct. The height
of the first mesh layer on the blades is set to $6\times 10^{-3}R_{out}$, and
the first off-wall solution point is located $2.8\times 10^{-4}R_{out}$ away
from the blades (for the fifth-order scheme).

(a)

(b)

Figure 6: Computational mesh of the ducted configuration: (a) rotor and duct,
(b) sliding interface.

Fig. 7 illustrates the surface meshes on the rotor and sliding interfaces for
the open-rotor configuration. It can be seen that the mesh distribution of the
rotor is almost the same as that of the ducted one, and this rotating
subdomain has three sliding interfaces: two in the axial direction, and one in
the azimuthal direction. The radius of this rotating part is $1.5R_{out}$, and
is larger than that of the ducted case because we need more space to achieve a
smooth transition from the blade tip to the azimuthal sliding interface on
which an orthogonal mesh distribution is needed.

(a)

(b)

Figure 7: Computational mesh of the open-rotor configuration: (a) rotor, (b)
sliding interfaces.

Both the ducted and open-rotor cases have the same cylindrical farfield as
shown in Fig. 8. The inlet is located $6R_{out}$ ahead of the rotor, and its
radius is $30R_{out}$. The outlet is placed $24R_{out}$ behind the rotor. The
blockage ratio is $0.44\%$. For the ducted case, the rotating part has 29,469
hexahedral elements, and the stationary part has 262,941 elements. In total,
the mesh has 292,410 elements. A fifth-order FR scheme in space is employed,
and the resulting number of degrees of freedom is $36.6$ million. The open
case has 38,454 hexahedral elements on the rotating subdomain, and 258,621
elements on the stationary subdomain. The resulting number of degrees of
freedom is $37.1$ million.

For the boundary conditions, the inlet is treated as a Dirichlet boundary, the
outlet, and the outer boundary as characteristic farfields. The duct and blade
surfaces are set as no-slip adiabatic walls. The freestream Mach number
$Ma_{\infty}$ is set to 0.08 to ensure that the compressibility effect is
negligible.

(a)

(b)

(c)

Figure 8: Computational mesh of the ducted wind turbine: (a) overall view, (b)
front view (c) lateral view.

###  3.3 Nondimensional Parameters

Table 1: Operating conditions. Parameter | Value
---|---
Freestream velocity | 10 m/s
Kinematic vicosity | $1.4616\times 10^{-5}$
Rotation axis | x-direction
Rotation speed | -26.18 rad/s

The operating conditions for the wind turbine are listed in Table 1. Based on
this table, two governing nondimensional parameter can be obtained. One is the
Reynolds’ number $Re=U_{\infty}R/v=1.25\times 10^{6}$, where $U_{\infty}$ is
the freestream velocity, $R$ is the duct’s exit radius, and $v$ is kinematic
viscosity. Since the Reynolds number is very high, the flow around the turbine
will become quite turbulent, and the simulation may become unstable due to the
aliasing error[18], To reduce the instability, a filter reported in [40] (with
strength $\alpha=0.5$) is employed during the simulation.

Another governing parameter is the nondimensional angular speed
$\omega^{*}=\omega R/U_{\infty}=-4.79$. For the present work, the
nondimensional time step size is $\Delta tR/U_{\infty}=2.5\times 10^{-5}$,
which means that the turbine rotates 0.0069 degrees per time step and thus the
turbine’s rotation can be simulated accurately.

All the simulations are run over $37.5$ time units (about 27 revolutions).
From $tR/U_{\infty}=25$, we start to collect data to compute the averaged flow
field. The averaging calculation is performed over around $9.5$ revolutions.
For the rotating part, we can obtain the phase-averaged statistics. For the
stationary part, the time-averaged statistics can be computed.

##  4 Results and Discussions

###  4.1 Axial Flow

####  4.1.1 Load Analysis

To study the loads on the turbine quantitatively, two performance
characteristics are calculated: thrust coefficient $C_{T}$, and power
coefficient $C_{P}$. They are defined as below.

$\displaystyle C_{T}=\frac{T}{\frac{1}{2}\rho U_{\infty}^{2}A_{rot}}$ (29)
$\displaystyle C_{P}=\frac{Q\omega}{\frac{1}{2}\rho U_{\infty}^{3}A_{rot}}$
(30)

where $T$ is the thrust acted on the rotor, $Q$ is the torque in the
x-direction, and $A_{rot}$ is the area of the rotor rotation plane. It should
be noted that in this study we only consider the thrust acted on the rotor
rather than on all the components.

(a) Ducted

(b) Open-rotor

Figure 9: Time histories of loads in the two configurations

Fig. 9 demonstrates the time histories of $C_{T}$ and $C_{P}$ for the two
configurations. In both cases, thrust and power outputs have reached a steady
state with high-frequency fluctuations caused by the turbulent flows.

Table 2 lists the mean values and the r.m.s (root-mean-square) deviations of
$C_{T}$ and $C_{P}$, their pressure and viscous components. For both
configurations, the pressure contributions for $C_{T}$ and $C_{P}$ are
dominant, and the r.m.s values are about two magnitudes smaller than
corresponding mean values. It can be seen that the rotor in the ducted case
has larger mean loads than its open-rotor counterpart. Both the mean $C_{T}$
and the mean $C_{P}$ of the ducted case are about $50\%$ higher than that of
the open-rotor case. The ducted case also has a higher root-mean-square
deviations for $C_{T}$ and $C_{P}$ than the open-rotor case, i.e., the load
fluctuations for the ducted one are stronger.

Table 2: Loads for the ducted and open-rotor cases. |  | $C_{T}$ | $C_{T,p}$ | $C_{T,v}$ | $C_{P}$ | $C_{P,p}$ | $C_{P,v}$
---|---|---|---|---|---|---|---
Ducted | mean | 9.47E-01 | 9.47E-01 | 1.97E-04 | 5.60E-01 | 5.61E-01 | -1.19E-03
r.m.s | 1.43E-02 | 1.43E-02 | 1.46E-06 | 8.30E-03 | 8.29E-03 | 1.25E-05
Open-rotor | mean | 6.30E-01 | 6.30E-01 | 1.97E-04 | 3.61E-01 | 3.62E-01 | -1.31E-03
r.m.s | 6.13E-03 | 6.13E-03 | 8.28E-07 | 2.03E-03 | 2.03E-03 | 4.81E-06

To better understand the load distributions on the blades, we plot the phase-
averaged pressure contours on the blade surfaces in Figs. 10 and 11. For the
ducted configuration, the upstream face experiences a higher pressure than the
downstream face. Thus the rotor experience a drag force pointed towards
x-direction. The outboard part of blades, especially the area around the
leading edge, has a higher pressure difference, which means the loads
concentrate in this region.

(a) upstream

(b) downstream

Figure 10: Phase-averaged pressure contours of the blades in the ducted
configuration.

For comparison, the open rotor’s phase-averaged pressure contours are given in
Fig. 11. A similar pressure distribution as the ducted one is observed. One
major difference is that the downstream face has smaller low pressure area
than that of the ducted configuration, which leads to a lower $C_{T}$ than the
ducted one.

(a) upstream

(b) downstream

Figure 11: Phase-averaged pressure contours of the blades in the open-rotor
configuration.

For the ducted case, the phase-averaged pressure contours at different
streamwise locations are shown in Fig. 12. It can be seen that due to the
pressure difference between two sides of the blade, all three blade sections
experience a torque pointed towards the negative x-axis. This torque direction
is the same as the rotation direction, which means a positive work is done on
the rotor and the wind energy is transferred to the rotor’s mechanical energy.
High pressure difference appears on the outboard region, which means most of
torque comes from this region.

(a) x=0.03R

(b) x=0.04R

(c) x=0.05R

Figure 12: Phase-averaged pressure contours of the ducted configuration at
different streamwise locations.

We also plot the phase-averaged pressure field around the open rotor in Fig.
13. A similar trend for pressure distribution is observed. The rotor in this
configuration also experiences a torque in the same direction as its rotation,
and has the main load on the blades’ outboard part.

(a) x=0.03R

(b) x=0.04R

(c) x=0.05R

Figure 13: Phase-averaged pressure contours of the open-rotor configuration at
different streamwise locations.

####  4.1.2 Wake Dynamics

To investigate the wake dynamics of the two configurations, isosurfaces of
instantaneous Q-criterion colored by the streamwise velocity are plotted in
Figs. 14 and 15. For the ducted case, as shown in Fig.14, a lot of vortices
are shed from the duct’s inner surface. The wake area is wrapped by these
small vortex structures from the duct exit to $x=10R$. To obtain more details,
in Fig. 14, we also plot isosurfaces of Q-criterion with the area $r>=0.5R$
blanked. It is seen that a lot of trailing vortices are generated, and they
are well resolved even up to $x=10R$, which means the wake area is filled with
small flow structures.

Figure 14: Isosurface of instantaneous Q-criterion $Q_{cr}R^{2}/U_{\infty}=15$
of the ducted case.

Figure 15: Isosurface of instantaneous Q-criterion $Q_{cr}R^{2}/U_{\infty}=15$
of the open-rotor case.

Fig.15 presents Q-criterion isosurfaces of the open-rotor case. As can be
seen, this case has much stronger and more coherent tip vortices compared with
the ducted case. When the area $r>=0.5R$ is blanked, we can clearly see that a
hub vortex is formed just behind the hub end and remains stable for a very
long distance. Compared with the ducted case, the flow around the open rotor
does not have strong vortices shedding off from blade trailing edges, and thus
has a load with smaller r.m.s deviation values.

The velocity vector can be decomposed into three components: streamwise
velocity $u$, radial velocity $u_{r}$, and azimuthal velocity $u_{\theta}$.
Time-averaged values of these components are visualized in the central $x-y$
plane in Figs. 16, 17, and 18.

(a) Dcuted

(b) Open-rotor

Figure 16: Time-averaged streamwise velocity contours in the central x-y plane

As shown in Fig.16, both the ducted and open-rotor cases have an area with low
streamwise velocity behind the hub end. For the ducted configuration, the low-
velocity area expands in the radial direction as the flow moves downstream.
After $x=4R$, this area begins to shrink. On the contrary, the low-speed area
of the open-rotor case almost keeps its radial size until $x=9R$, and then
gradually vanishes, of which the shape is similar as the hub vortex shown in
Fig. 15.

(a) Dcuted

(b) Open-rotor

Figure 17: Time-averaged radial velocity contours in the central x-y plane

Fig. 17(a) presents the time-averaged radial velocity contour for the ducted
case. The area around the duct has a large positive radial velocity because
the duct’s diverging shape make the flow travels away from the central axis.
Another high radial velocity region appears near the hub end, which means that
the flow expands after it moves past the end. For the open-rotor case, there
is no high radial velocity region like the ducted one.

(a) Dcuted

(b) Open-rotor

Figure 18: Time-averaged azimuthal velocity contours in the central x-y plane

Fig. 18 demonstrates the time-averaged azimuthal velocity fields in the x-y
plane. In the open-rotor case, two high azimuthal velocity zones near the
central line appear, which means that at this place the flow rotates quickly
around the central axis, and is consistent with the hub vortex seen in Fig.
15. For the ducted case, only the flow around the hub end rotates quickly.

(a) x=R

(b) x=3R

(c) x=5R

(d) x=7R

(e) x=9R

Figure 19: Time-averaged streamwise velocity contours for the ducted case at
different streamwise locations

To obtain more details of the streamwise velocity field, we plot the time-
averaged streamwise velocity contours at five different streamwise locations
for the ducted and open-rotor configurations in Figs. 19 and 20, respectively.
From Fig. 19(a), it can be observed that a circular low-velocity region at the
central part, which are circumscribed by two ring regions with higher
velocity. The central circular area corresponds to the wake just behind the
hub end, and the inner ring region corresponds to the wake behind the area
between the duct and the hub. As the flow moves downstream, boundaries between
these regions become unclear due to the momentum exchange between different
regions. At $x=9R$, these areas are totally mixed together.

(a) x=R

(b) x=3R

(c) x=5R

(d) x=7R

(e) x=9R

Figure 20: Time-averaged streamwise velocity contours for the open-rotor case
at different streamwise locations.

In Fig. 20, we can also see three regions for the open-rotor configuration.
However, the area of the central low-velocity circular region is much less
than that in Fig.19, The inner ring region is formed by the flow past the
rotor’s blades. The central area keeps its shape until x=9R. The boundary
between two ring areas become obscure after $x=3R$. However, these region are
not blended together even when the flow moves downstream to $x=9R$. One
explanation for this difference is that the ducted turbine’s wake is filled
with small flow structures, and thus has a higher rate of momentum exchange,
which leads to the difference of velocity between regions disappear more
quickly.

To quantitatively compare the two configurations’ wakes, time-averaged
streamwise velocity profiles are plotted in Fig. 21. Overall these profiles
are symmetrical about $y=0$. The wake behind the ducted turbine has more
velocity loss. When $x=R$, the ducted configuration has an almost stationary
flow in the area $-0.5\leq y/R\leq 0.5$, which corresponds to the hub wake. As
$\absolutevalue{y}$ increases, the velocity of the ducted turbine increases
quickly, and then has two local maximum $\overline{u}/U_{\infty}=0.3$ at
$y/R\approx 1$, i.e. the duct exit radius. The flow at $x=R$ continues to
recover as $\absolutevalue{y}$ increase, and reach the freestream velocity
$\overline{u}/U_{\infty}=1$ when $\absolutevalue{y}/R\geq 1.2$. By contrast,
the open-rotor case has a spike around $x=0$ with a local minimum
$\overline{u}/U_{\infty}=0$, and recovers quickly to
$\overline{u}/U_{\infty}=1$ when $\absolutevalue{y}/R\geq 1$.

When the flow moves downstream, the ducted case’s velocity profiles has only
one minimum for $x\geq 3R$, and is still smaller than that of the open-rotor
one when $x\leq 7R$. When $x=9R$, the two profiles have a similar shape with a
minimum value $\overline{u}/U_{\infty}\approx 0.2$. Among other regions, the
ducted case’s velocity is smaller than the open-rotor one. The open-rotor
configuration keeps its spike at $y=0$ until $x=7R$. The locations where two
profiles reach the steady state remain unchanged.

(a) x=R

(b) x=3R

(c) x=5R

(d) x=7R

(e) x=9R

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

###  4.2 Yawed Flow

This section focuses on the yawed flows. To facilitate the description, two
geometric parameters: yaw angle and azimuthal angle, are defined and shown in
Fig. 22. In this study, yaw angle is denoted as the angle between the x-axis
and the flow direction in the central x-y plane. Azimuthal angle is defined as
angle between the blade axis and the horizontal line shown in Fig.22. Initial
azimuthal angles for three blades are $30^{\circ}$, $150^{\circ}$, and
$270^{\circ}$.

(a) Yaw angle

(b) Azimuthal angle

Figure 22: Two angle definitions0 for yawed flow.

For the yawed flow, we consider five different yaw angles: $0^{\circ}$, $\pm
15^{\circ}$, and $\pm 30^{\circ}$ for two configurations. Fig 23 compares the
time-averaged $C_{T}$ of two configurations for different yawed conditions. It
can be seen that both the ducted turbine and the open-rotor turbine’s
$\overline{C}_{T}$ are symmetrical about $\gamma=0^{\circ}$, and decrease with
the growing yaw angle magnitude. Even though the ducted turbine experiences a
larger drop on $\overline{C}_{T}$(about $0.2$) from $0^{\circ}$ to
$30^{\circ}$ than the open-rotor one (about $0.04$), the ducted rotor
experiences a larger thrust on every yawed angle. To compare the relative
change, in Fig. 23(b) we normalize both configurations’ $\overline{C}_{T}$ by
their own $\overline{C}_{T}$ at $\gamma=0^{\circ}$. It can be found that
relative changes of $\overline{C}_{T}$ for the open-rotor cases are smaller
than that of the ducted ones. Especially, for $\gamma=\pm 30^{\circ}$, the
open rotor only experiences about $6\%$ thrust drop, much less than the 20$\%$
for the ducted case.

(a) $C_{T}$

(b) Normalized $C_{T}$

Figure 23: Time-averaged $C_{T}$ for two configurations under yawed flow.

As shown in Fig. 24, the $\overline{C}_{P}$ curves for yawed conditions have a
similar trend as the $\overline{C}_{T}$ curves, which have a symmetry about
$\gamma=0^{\circ}$ and see a drop as the yaw angle grows. The ducted turbine
has more power output than the open rotor in all yawed cases, and experiences
a drop by $0.12$ on its $\overline{C}_{P}$ value when $\gamma=\pm 30^{\circ}$.
In comparison, the open rotor only decreases by $0.05$. In terms of the
relative change in $\overline{C}_{P}$, the ducted configuration has a much
larger drop (about $30\%$) for $\gamma=\pm 30^{\circ}$, whereas the open-rotor
case only drops by about $18\%$. The $\overline{C}_{P}$ values are symmetric
with respect to $\gamma=0^{\circ}$, and decrease as $\gamma$ increases. The
ducted rotor have more power output than the open rotor on all yawed angles.
For $\gamma=\pm 15^{\circ}$, relative changes for two configurations are
similar (around 5%). However, the ducted configuration has a much larger
relative decrease on power output (about $30\%$) for $\gamma=\pm 30^{\circ}$.
The open-rotor case only experiences around $18\%$ decrease on
$\overline{C}_{P}$.

(a) $C_{P}$

(b) Normalized $C_{P}$

Figure 24: Time-averaged $C_{P}$ for the two configurations under yawed flows.

$\widetilde{C}_{T}$, i.e., the azimuthal-averaged $C_{T}$ can be computed by
averaging $C_{T}$ values at the time instant when the azimuthal angle is equal
to a given value, like $90^{\circ}$. Fig. 25 shows $\widetilde{C}_{T}$ of
blade 1 for $\gamma=0,\pm 15^{\circ}$, which is normalized by the
corresponding $\overline{C}_{T}$. We see that all yawed cases have a sine-like
fluctuation during one period. One major cause is the blade advancing and
retreating effects reported in [41]. When the turbine blade rotates in yawed
flows, it will move upstream during half of the rotation period, and move
downstream during the other half of the period, which results in the variation
of the relative velocity and then the variation of angle of attack. Thus a
fluctuation in thrust is seen during one revolution. When $\gamma=0^{\circ}$,
$C_{T}$ remain stable around $\overline{C}_{T}$ for both configurations. Both
cases have a $180^{\circ}$ phase difference for $C_{T}$ values between
$\gamma=15^{\circ}$ and $\gamma=-15^{\circ}$. For $\gamma=-15^{\circ}$, ducted
cases has a maximum at $\theta=180^{\circ}$, whereas open-rotor cases has its
maximum at around $\theta=120^{\circ}$ and has a larger amplitude (about
$1.1$).

(a) Ducted

(b) Open-rotor

Figure 25: Azimuthal-averaged $C_{T}$ of the blade 1 for yaw angle
$\gamma=0,\pm 15^{\circ}$.

As shown in Fig. 26, For $\gamma=\pm 30^{\circ}$, these $C_{T}$ curves have
the same phase difference between the positive yaw angle and the negative yaw
angle. Both configurations’ maximum values appear at $\theta=120^{\circ}$.
However the ducted one has a larger amplitude (about $1.3$), and has another
local maximum at around $\theta=250^{\circ}$ due to the duct’s effects.

(a) Ducted

(b) Open-rotor

Figure 26: Azimuthal-averaged $C_{T}$ of the blade 1 for yaw angle
$\gamma=0,\pm 30^{\circ}$.

##  5 Conclusion

In this study, we perform high-order large eddy simulations for a wind turbine
under two configurations: ducted and open-rotor.

For the axial flow situation, a comparison of the thrust coefficient $C_{T}$
and the power coefficient $C_{P}$ between the two cases is made. It is found
that the ducted one has higher mean values for thrust and power output than
its bare counterpart, and also has higher r.m.s deviation values, which means
stronger load fluctuations. To better understand the load distribution, time-
averaged pressure fields on the blade surfaces and near-blade regions. We see
that the two configurations experience a drag force pointed towards the x
axis, and a torque pointed towards the negative x axis, and the load
concentrates on the outboard part of the blades, especially around the leading
edges. To investigate the wake dynamics, Q-criterion isosurfaces are plotted.
It is found that the wake of the ducted turbine is filled with small flow
structures, whereas the open rotor has a coherent hub vortex and weak tip
vortices that are dissipated quickly by the flow. The time-averaged velocity
contours for three components: $u$, $u_{r}$, and $u_{\theta}$ are also
presented. We found that the ducted case has a larger low-velocity wake area
behind the hub end. In comparison, the open-rotor configuration has a thinner
band region with low streamwise velocity and high azimuthal velocity.
Streamwise velocity contours at different streamwise locations reveal that the
wake behind the ducted rotor has a faster rate of momentum exchange.

For the yawed flow situations, we compare $C_{T}$ and $C_{P}$ of the two
configurations. For all tested yaw angles, the ducted one has a higher thrust
and power output than its bare counterpart. However, the ducted case
experiences a larger drop on $C_{T}$ and $C_{P}$ in both absolute and relative
senses. We also compute $\widetilde{C}_{T}$, i.e., the azimuthal-averaged
$C_{T}$, for one blade. Due to the blade advancing and retreating effect, each
yawed case’s curve has some fluctuations. And there is a $180^{\circ}$ phase
difference between two flows with the same magnitude but opposite signs.

To better understand the performance of the ducted wind turbine, we will
perform simulations for different tip speed ratios in the future.

## Acknowledgments

The authors would like to express our acknowledgments to Clarkson University
for financial support through the IGNITE Fellowship Program.

## References

  * Oka et al. [2013] Oka, N., Furukawa, M., Yamada, K., Kawamitsu, K., Kido, K., and Oka, A., “Aerodynamic design optimization of wind-lens turbine,” _ASME 2013 Fluids Engineering Division Summer Meeting_ , 2013.
  * Fletcher [1981] Fletcher, C. A., “Computational analysis of diffuser-augmented wind turbines,” _Energy Conversion and Management_ , Vol. 21, No. 3, 1981, pp. 175–183.
  * Georgalas et al. [1991] Georgalas, G., Koras, A., and Raptis, S., “Parameterization of The Power Enhancement Calculated for Ducted Rotors with Large Tip Clearance.” _Wind Engineering_ , Vol. 15, 1991.
  * Hansen et al. [2000] Hansen, M. O. L., Sørensen, N. N., and Flay, R. G. J., “Effect of Placing a Diffuser around a Wind Turbine,” _Wind Energy_ , Vol. 3, No. 4, 2000, pp. 207–213.
  * Ohya et al. [2012] Ohya, Y., Uchida, T., Karasudani, T., Hasegawa, M., and Kume, H., “Numerical Studies of Flow around a Wind Turbine Equipped with a Flanged-Diffuser Shroud Using an Actuator-Disk Model,” _Wind Engineering_ , Vol. 36, No. 4, 2012, pp. 455–472.
  * Venters et al. [2018] Venters, R., Helenbrook, B. T., and Visser, K. D., “Ducted wind turbine optimization,” _Journal of Solar Energy Engineering_ , Vol. 140, No. 1, 2018\.
  * Wang et al. [2013] Wang, Z. J., Fidkowski, K., Abgrall, R., Bassi, F., Caraeni, D., Cary, A., Deconinck, H., Hartmann, R., Hillewaert, K., Huynh, H. T., et al., “High-order CFD methods: current status and perspective,” _International Journal for Numerical Methods in Fluids_ , Vol. 72, No. 8, 2013, pp. 811–845.
  * Reed and Hill [1973] Reed, W. H., and Hill, T. R., “Triangular mesh methods for the neutron transport equation,” Tech. rep., Los Alamos Scientific Lab., N. Mex.(USA), 1973\.
  * Cockburn et al. [2012] Cockburn, B., Karniadakis, G. E., and Shu, C.-W., _Discontinuous Galerkin methods: theory, computation and applications_ , Vol. 11, Springer Science & Business Media, 2012.
  * Patera [1984] Patera, A. T., “A spectral element method for fluid dynamics: laminar flow in a channel expansion,” _Journal of computational Physics_ , Vol. 54, No. 3, 1984, pp. 468–488.
  * Karniadakis and Sherwin [2013] Karniadakis, G., and Sherwin, S., _Spectral/hp element methods for computational fluid dynamics_ , Oxford University Press, 2013.
  * Kopriva and Kolias [1996] Kopriva, D. A., and Kolias, J. H., “A conservative staggered-grid Chebyshev multidomain method for compressible flows,” _Journal of computational physics_ , Vol. 125, No. 1, 1996, pp. 244–261.
  * Kopriva [1996] Kopriva, D. A., “A conservative staggered-grid Chebyshev multidomain method for compressible flows. II. A semi-structured method,” _Journal of computational physics_ , Vol. 128, No. 2, 1996, pp. 475–488.
  * Kopriva [1998] Kopriva, D. A., “A staggered-grid multidomain spectral method for the compressible Navier–Stokes equations,” _Journal of computational physics_ , Vol. 143, No. 1, 1998, pp. 125–158.
  * Liu et al. [2006] Liu, Y., Vinokur, M., and Wang, Z. J., “Spectral difference method for unstructured grids I: Basic formulation,” _Journal of Computational Physics_ , Vol. 216, No. 2, 2006, pp. 780–801.
  * Huynh [2007] Huynh, H. T., “A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods,” AIAA paper 2007-4079, 2007.
  * Huynh [2009] Huynh, H. T., “A reconstruction approach to high-order schemes including discontinuous Galerkin for diffusion,” AIAA paper 2009-403, 2009.
  * Jameson et al. [2012] Jameson, A., Vincent, P. E., and Castonguay, P., “On the non-linear stability of flux reconstruction schemes,” _Journal of Scientific Computing_ , Vol. 50, No. 2, 2012, pp. 434–445.
  * Castonguay et al. [2011] Castonguay, P., Williams, D., Vincent, P., Lopez, M., and Jameson, A., “On the development of a high-order, multi-GPU enabled, compressible viscous flow solver for mixed unstructured grids,” _20th AIAA Computational Fluid Dynamics Conference_ , 2011, p. 3229.
  * Haga et al. [2010] Haga, T., Gao, H., and Wang, Z., “A high-order unifying discontinuous formulation for 3-D mixed grids,” _48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition_ , 2010, p. 540.
  * Haga et al. [2011] Haga, T., Gao, H., and Wang, Z. J., “A high-order unifying discontinuous formulation for the Navier-Stokes equations on 3D mixed grids,” _Mathematical Modelling of Natural Phenomena_ , Vol. 6, No. 3, 2011, pp. 28–56.
  * Lu et al. [2012] Lu, Y., Yuan, X., and Dawes, W., “Investigation of 3D internal flow using new flux-reconstruction high order method,” _Turbo Expo: Power for Land, Sea, and Air_ , Vol. 44748, American Society of Mechanical Engineers, 2012, pp. 2195–2216.
  * Vincent et al. [2015] Vincent, P., Farrington, A., Witherden, F., and Jameson, “An extended range of stable-symmetric-conservative Flux Reconstruction correction functions,” _Computer Methods in Applied Mechanics and Engineering_ , Vol. 296, 2015, pp. 248–272.
  * Wang and Huynh [2016] Wang, Z. J., and Huynh, H. T., “A review of flux reconstruction or correction procedure via reconstruction method for the Navier-Stokes equations,” _Mechanical Engineering Reviews_ , Vol. 3, No. 1, 2016, pp. 15–00475–15–00475.
  * Zhang and Liang [2015a] Zhang, B., and Liang, C., “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] Zhang, B., and Liang, C., “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_ , Vol. 295, 2015b, pp. 147–160.
  * Zhang et al. [2016] Zhang, B., Liang, C., Yang, J., and Rong, Y., “A 2D parallel high-order sliding and deforming spectral difference method,” _Computers & Fluids_, Vol. 139, 2016, pp. 184–196.
  * Zhang and Liang [2016] Zhang, B., and Liang, C., “A high-order sliding-mesh spectral difference solver for simulating unsteady flows around rotating objects,” _31st Symposium on Naval Hydrodynamics_ , Monterey, CA, 2016.
  * Zhang et al. [2018] Zhang, B., Qiu, Z., and Liang, C., “A flux reconstruction method with nonuniform sliding-mesh interfaces for simulating rotating flows,” AIAA paper 2018-1094, 2018.
  * Zhang and Liang [2020] Zhang, B., and Liang, C., “A conservative high-order method utilizing dynamic transfinite mortar elements for flow simulation on curved sliding meshes,” arXiv preprint arXiv:2009.04400, 2020.
  * Zhang and Liang [2021] Zhang, B., and Liang, C., “A high-order sliding-mesh method for flow simulation about complex 3D rotating geometries,” Under review, 2021.
  * Zhang and Liang [2019a] Zhang, B., and Liang, C., “High-order numerical simulation of flows over rotating cylinders of various cross-sectional shapes,” AIAA paper 2019-3430, 2019a.
  * Zhang and Liang [2019b] Zhang, B., and Liang, C., “High-order numerical simulation of flapping wing for energy harvesting,” AIAA paper 2019-3338, 2019b.
  * Zhang et al. [2021] Zhang, B., Ding, C., and Liang, C., “High-Order Implicit Large-Eddy Simulation of Flow over a Marine Propeller,” _Computers & Fluids_, 2021, p. 104967.
  * Thomas and Lombard [1979] Thomas, P. D., and Lombard, C. K., “Geometric conservation law and its application to flow computations on moving grids,” _AIAA Journal_ , Vol. 17, 1979, pp. 1030–1037.
  * Bathe [2006] Bathe, K.-J., _Finite Element Procedures_ , Klaus-Jurgen Bathe, 2006.
  * Rusanov [1961] Rusanov, V. V., “Calculation of interaction of non-steady shock waves with obstacles,” _Journal of Computational and Mathematical Physics USSR_ , Vol. 1, 1961, pp. 267–279.
  * Spiteri and Ruuth [2002] Spiteri, R. J., and Ruuth, S. J., “A new class of optimal high-order strong-stability-preserving time discretization methods,” _SIAM Journal on Numerical Analysis_ , Vol. 40, 2002, pp. 469–491.
  * Ruuth [2006] Ruuth, S., “Global optimization of explicit strong-stability-preserving Runge-Kutta methods,” _Mathematics of Computation_ , Vol. 75, No. 253, 2006, pp. 183–207.
  * Fischer and Mullen [2001] Fischer, P., and Mullen, J., “Filter-based stabilization of spectral element methods,” _Comptes Rendus de l’Academie des Sciences Series I Mathematics_ , Vol. 332, No. 3, 2001, pp. 265–270.
  * Schulz et al. [2017] Schulz, C., Letzgus, P., Lutz, T., and Krämer, E., “CFD study on the impact of yawed inflow on loads, power and near wake of a generic wind turbine,” _Wind Energy_ , Vol. 20, No. 2, 2017, pp. 253–268.