# Development and Preliminary Assessment of the Open-Source CFD toolkit SU2
for Rotorcraft Flows

Myles Morelli  Tommaso Bellosta  Alberto Guardone  Department of Aerospace
Science and Technology, Politecnico di Milano, Italy  CFD Laboratory, School
of Engineering, University of Glasgow, United Kingdom

###### Abstract

Computational aerodynamic analyses of rotorcraft main rotor blades are
performed in both hover and forward flight. The open-source SU2 code is used
for rotor performance prediction. The core of the code is the set of RANS
equations, which are solved for determining the flow. In hover, both steady-
state and time-accurate modelling techniques of varying complexity are used
and assessed. Simulation specific parameters which have a significant
influence on the solution are also addressed. In forward flight, the code is
developed to include the main rotor blade kinematics which is a prerequisite
for modelling a trimmed rotor. Two databases are used for the validation of
the rotor performance prediction. The renowned Caradonna-Tung experimental
tests of a model rotor are used to evaluate the pressure distribution along
the blade during hover. The extensive aerodynamic and aeroacoustic data survey
of the AH-1G Cobra helicopter is used to assess the pressure distribution at
different advancing and retreating azimuth angles during forward flight. The
prediction capabilities of the solver in terms of rotor performance are
demonstrated and are overall in good agreement with the measured data.

###### keywords:

Rotorcraft , CFD , Open-Source

††journal: Journal of Computational and Applied Mathematics,

##  1 Introduction

Rotorcraft have the unique ability to be able to vertically take-off and land.
This allows for their operation in highly demanding and challenging flight
scenarios where conventional fixed-wing aircraft cannot pursue. Their frequent
use for demanding operations however can lead to dangerous in-flight
situations which can be seen from the relatively high number of accident
reports from the National Transportation Safety Board (NTSB) [1]. Tools that
allow the analysis of rotorcraft are therefore essential for supporting design
and for reducing accident rates. Sustained computational development over the
past 30 years [2] means that many codes are now capable of modelling
rotorcraft behaviour. Computational advancements have enabled CFD-based
methods to simulate and further understand complex rotorcraft aerodynamics
which are rich in flow physics.

NASA have a long history of rotorcraft CFD code development with OVERFLOW [3]
and FUN3D [4]. Of the two codes, OVERFLOW has been used more extensively for
rotorcraft simulations. OVERFLOW uses body-fitted structured grids near solid
surfaces and automatically generated Cartesian grids in the background and is
a finite-difference node-based solver. FUN3D uses unstructured grids
throughout the domain and is a finite volume node-based solver. Academia also
plays a leading role in rotorcraft CFD code development. In-house codes HMB
from the University of Glasgow [5] and ROSITA from Politecnico di Milano [6]
are both finite-volume solvers and utilize structured Chimera multi-block
grids to account for the blade motion. Each of the codes have been part of
large-scale collaborations for the assessment of their predictive capabilities
such as the GOAHEAD project [7]. Another code developed within academia is the
TURNS research code from the University of Maryland [8, 9]. The TURNS code has
Chimera overset grid capabilities and uses a finite difference numerical
algorithm that evaluates the inviscid fluxes using an upwind-biases flux
scheme. European research centers also have their own codes which have been
heavily developed. French and German research institutes ONERA and DLR have
their own multi-block Chimera based rotorcraft codes. ONERA have the elsA
solver [10] and DLR have the FLOWer solver [11].

Despite all of the codes showing excellent predictive capabilities, none of
the mentioned codes are freely available in the open-source domain which
hinders technological developments. Furthermore, Validation and Verification
(VnV) is shown to be increasingly important so to identify a wide variety of
physical modelling, discretization, and solution errors [12]. To that end, an
open-source code with an active and growing community of users provides a
platform for extensive VnV and innovative new solutions.

This work looks to provide the first open-source and validated rotorcraft CFD
code by developing the well established SU2 code [13]. The main contribution
of this work is the introduction of the blade motion and an approach for
modelling forward flight. The contemporary nature of the code means that it
can benefit from new and effective numerical techniques. The recent
implementation by Gori et al. [14] of a method for dealing with non-conformal
boundary interfaces using the supermesh technique is one example of this.
Another example of this is radial basis function mesh deformation which can
now utilize computationally efficient algorithms for dealing with large-scale
problems. The outline of this paper is as follows; the physical modelling of
the RANS equations is discussed in Section 2, the numerical implementation is
described in Section 3, the validation work and results of an isolated rotor
in hover and forward flight are discussed in Section 4, and finally the main
talking points of the work are concluded in Section 5.

##  2 Physical Modelling

The following section on the physical modelling of the Reynolds-averaged
Navier-Stokes equations inside of the SU2 code is primarily a summary of the
work from Economon et al. [13] and is here to give context to the numerical
implementation in Section 3.

###  2.1 Reynolds-Averaged Navier-Stokes Equations

Within this framework, we are interested in time-accurate turbulent flow
around rotor blades with arbitrary motion. Therefore we are concerned with
compressible flow governed by the Reynolds-averaged Navier-Stokes (RANS)
equations. These mass, momentum, and energy conservation equations can be
expressed in arbitrary Lagrangian-Eulerian differential form as,

$\left\\{\begin{matrix}\begin{array}[]{cclllc}\mathcal{R}(\boldsymbol{U})&=&\frac{\partial
U}{\partial
t}+\triangledown\cdot\boldsymbol{F}_{ale}^{c}-\triangledown\cdot(\mu_{tot}^{1}\boldsymbol{F}^{\upsilon
1}+\mu_{tot}^{2}\boldsymbol{F}^{\upsilon
2})-\boldsymbol{Q}=0&\textrm{in}&\Omega,&t>0\\\
\boldsymbol{\upsilon}&=&\boldsymbol{u}_{\Omega}&\textrm{on}&S&\\\
\partial_{n}T&=&0&\textrm{on}&S&\\\
(W)_{+}&=&W_{\infty}&\textrm{on}&\Gamma_{\infty}&\end{array}\end{matrix}\right.$
(1)

with the vector of conservative variables being represented by,
$\boldsymbol{U}=\left\\{\rho,\,\rho\boldsymbol{\upsilon},\,\rho
E\right\\}^{T}$, inside the flow domain, $\Omega$. The term,
$\boldsymbol{\upsilon}=\boldsymbol{u}_{\Omega}$, denotes the no-slip condition
on the surface, $S$. The expression, $\partial_{n}T=0$, represents the
adiabatic condition on the surface, $S$. The final condition,
$(W)_{+}=W_{\infty}$, is the characteristic-based boundary condition at the
far-field, $\Gamma_{\infty}$. The convective fluxes, viscous fluxes and source
terms within Eq. 1 can then be respectively expressed as,

$\boldsymbol{F}_{ale}^{c}=\begin{Bmatrix}\rho(\boldsymbol{\upsilon}-\boldsymbol{u}_{\Omega})\\\
\rho\boldsymbol{\upsilon}\otimes(\boldsymbol{\upsilon}-\boldsymbol{u}_{\Omega})+\bar{\bar{I}}_{p}\\\
\rho
E(\boldsymbol{\upsilon}-\boldsymbol{u}_{\Omega})+p\boldsymbol{\upsilon}\end{Bmatrix},\,\boldsymbol{F}^{\upsilon
1}=\begin{Bmatrix}\cdot\\\ \bar{\bar{\tau}}\\\
\bar{\bar{\tau}}\cdot\boldsymbol{\upsilon}\end{Bmatrix},\,\boldsymbol{F}^{\upsilon
2}=\begin{Bmatrix}\cdot\\\ \cdot\\\ c_{p}\triangledown
T\end{Bmatrix},\,\boldsymbol{Q}=\begin{Bmatrix}q_{\rho}\\\
\boldsymbol{q}_{\rho\upsilon}\\\ q_{\rho E}\end{Bmatrix}$ (2)

where the fluid density, flow velocity vector, static pressure, temperature,
and specific heat are respectively given by the orthodox notation
$\rho,\,\boldsymbol{\upsilon},\,p,\,T,\,c_{p}$. Of the remaining variables;
$E$, is the total energy per unit mass; $\bar{\bar{\tau}}$ is the viscous
stress tensor; $\boldsymbol{u}_{\Omega}$ is the grid velocity.

###  2.2 Turbulence Models

Here we are concerned with rotorcraft which during forward flight, climb,
descent, and manoeuvre can be characterised by unsteady turbulent flows. With
this in mind, the solution of the unsteady Reynolds average Navier-Stokes
equations need to be solved which requires the inclusion of a turbulence
model. Using the Boussinesq hypothesis [15], the effect of turbulence is
represented as an increase in the viscosity. The total viscosity is then
separated into laminar and turbulent viscosity and can be respectively denoted
as $\mu_{dyn}$ and $\mu_{tur}$. The laminar viscosity is determined based upon
Sutherland’s law [16]. The total viscosity as part of the momentum and energy
equations in Eq. 1 is then substituted by

$\mu_{tot}^{1}=\mu_{dyn}+\mu_{tur},\quad\mu_{tot}^{2}=\frac{\mu_{dyn}}{Pr_{l}}+\frac{\mu_{tur}}{Pr_{t}}$
(3)

where the laminar and turbulent Prandtl numbers are given by $Pr_{l}$ and
$Pr_{t}$.

The turbulent viscosity is computed using a suitable turbulence model which
itself is dependent upon the flow state and a new set of variables,
$\hat{\nu}$, to represent the turbulence, such that,
$\mu_{tur}=\mu_{tur}(\boldsymbol{U},\,\hat{\nu})$. One of the most widely used
models for aeronautical attached flows is the one-equation Spallart-Allmaras
(SA) turbulence model [17].

###  2.3 Rotating Frame of Reference

The unique ability of rotorcraft to be able to hover as well as to climb and
descend vertically in axial flight helps to distinguish them from other
aircraft. During these specific operational flight conditions, the flow around
the main rotor can be considered as a steady rotation. Under this assumption,
it is then possible to transform the unsteady problem into a steady problem to
improve the efficiency of the simulation. This is possible by transforming the
system of governing equations in Eq. 1 into a rotating frame of reference
which rotates at the constant rotational velocity of the main rotor blades.
This modification to the system was implemented by Economon et al. [18] inside
the SU2 code and can be written as,

$\frac{\partial\boldsymbol{U}}{\partial
t}=0,\quad\boldsymbol{u}_{\Omega}=\boldsymbol{\omega}\times\boldsymbol{r},\quad\boldsymbol{Q}=\begin{Bmatrix}\cdot\\\
-\rho(\boldsymbol{\omega}\times\boldsymbol{r})\\\ \cdot\end{Bmatrix}$ (4)

where the rotational velocity vector of the rotating frame of reference is
specified by
$\boldsymbol{\omega}=\\{\omega_{x},\,\omega_{y},\,\omega_{z}\\}^{T}$ and the
radial distance from the center of rotation is specified by $\boldsymbol{r}$.

##  3 Numerical Implementation

Alongside the physical modelling in the paper from Economon et al. [13] there
is also a full description of the numerical implementation. In detail, it
discusses the spatial and temporal integration which are the foundation for
more application-specific methods such as the harmonic balance and non-
conformal boundary interface treatment. The concurrent implementation for
modelling rotor blade kinematics will finally be introduced.

###  3.1 Spatial Integration

The flow equations are solved numerically via a finite volume method [19]
which is applied on unstructured grids with an edge-based structure.
Integrating the governing equations over a control volume and using the
divergence theorem to obtain the semi-discretized form gives,

$\int_{\Omega_{i}}\frac{\partial\boldsymbol{U}}{\partial
t}d{\Omega}+\sum_{j\in\mathcal{N}(i)}(\tilde{\boldsymbol{F}_{ij}^{c}}+\tilde{\boldsymbol{F}_{ij}^{v}})\Delta
S_{ij}-\boldsymbol{Q}\left|\Omega_{i}\right|=\int_{\Omega_{i}}\frac{\partial\boldsymbol{U}}{\partial
t}d{\Omega}+\mathcal{R}_{i}(\boldsymbol{U})=0$ (5)

where $\mathcal{R}_{i}(\boldsymbol{U})$ includes the convective and viscous
fluxes integrated over the surface area of a control volume and any source
terms. The numerical approximations of the convective and viscous fluxes are
expressed respectively by, $\tilde{\boldsymbol{F}_{ij}^{c}}$ and
$\tilde{\boldsymbol{F}_{ij}^{v}}$. The area of the face belonging to the edge
$ij$ is represented by $\Delta S_{ij}$ and the set of neighbouring vertices to
vertex $i$ is expressed by $\mathcal{N}(i)$. The volume of control volume $i$
in the domain is denoted by $\left|\Omega_{i}\right|$.

The convective and viscous fluxes are evaluated at the mid-point of an edge.
Centered or upwind schemes such as the Jameson-Schmidt-Turkel scheme [20] or
the approximate Riemann solver of Roe [21] are then used to discretize the
convective fluxes. Second-order of accuracy of upwind schemes is achieved by
using Monotone UpStream-Centered schemes for conservation Laws (MUSCL) [22] to
reconstruct the variables on the cell interfaces. Green-Gauss or weighted
least-squares methods are used approximate the spatial gradients of the flow
at the cell faces to determine the viscous fluxes.

###  3.2 Time Integration

In conjunction with the spatial discretization, there is also the requirement
for temporal discretization of the governing equations and so Equation 5 is
further discretized over a control volume, $\left|\Omega_{i}\right|$, in time
such that it becomes,

$\frac{d}{dt}(\left|\Omega_{i}\right|\boldsymbol{U})+\mathcal{R}_{i}(\boldsymbol{U})=0$
(6)

and a backward Euler scheme can be used to evaluate the solution state at the
updated time for steady problems.

The temporal discretization for unsteady simulations is achieved using a dual-
time stepping approach which allows for second order accuracy in time [23].
Using this approach, an additional fictitious time is introduced so that the
unsteady problem becomes a series of pseudo-steady problems. The fictitious
time is introduced in front of the governing equations such that,

$\left|\Omega_{i}\right|\frac{\partial\boldsymbol{U}_{n}}{\partial\tau}+\mathcal{R}^{*}(\boldsymbol{U}_{n})=0$
(7)

with the subscript $n$ denoting the physical time level and to achieve second-
order backward difference in time,

$\mathcal{R}^{*}(\boldsymbol{U}_{n})=\frac{2}{3\Delta
t}\boldsymbol{U}_{n}+\frac{1}{\left|\Omega_{i}\right|_{n+1}}\left(\mathcal{R}(\boldsymbol{U}_{n})-\frac{2}{\Delta
t}\left|\Omega_{i}\right|_{n}\boldsymbol{U}_{n}+\frac{1}{2\Delta
t}\left|\Omega_{i}\right|_{n-1}\boldsymbol{U}_{n-1}\right)$ (8)

where the physical and fictitious time are expressed by $\Delta t$ and $\tau$
respectively. The convergence of every physical time in pseudo time results in
the modified residual equates to $\mathcal{R}^{*}(\boldsymbol{U}_{n})=0$
meaning it is equivalent to finding the state
$\boldsymbol{U}_{n}=\boldsymbol{U}_{n+1}$.

###  3.3 Harmonic Balance

The high computational cost of time-accurate methods means that it is often
desirable to use reduced-order models for simulating unsteady problems.
Adopting the harmonic balance method is beneficial as it can be used for
quasi-periodic flows dominated by a specific set of frequencies which need not
be integral multiples of each other [24]. The subsequent harmonic balance
method was implemented by Ref. [25] in SU2 and will briefly be summarised.

In the harmonic balance method the time operator, $\mathcal{D}_{t}$, is
introduced and is approximated using spectral interpolation. Applying the
spectral operator to the vector of conservative variables,
$\tilde{\boldsymbol{U}}$, which are now to be evaluated at $\mathcal{N}$ time
instances, one obtains,

$\mathcal{D}_{t}(\boldsymbol{U})\approx\mathcal{D}_{t}(\tilde{\boldsymbol{U}})$
(9)

and the harmonic balance operator $\mathcal{D}_{t}$ can eventually take the
form,

$\mathcal{D}_{t}(\tilde{\boldsymbol{U}})=\boldsymbol{H}\tilde{\boldsymbol{U}}$
(10)

where $\boldsymbol{H}$ is the spectral operator matrix. While now considering
$\tilde{\boldsymbol{U}}$ as the vector of conservative variables evaluated at
$\mathcal{N}$ time instances, Equation 7 can be rewritten for a single time
instance as,

$\left|\Omega_{i}\right|\frac{\Delta\boldsymbol{U}_{n}^{q+1}}{\Delta\tau_{n}}+\boldsymbol{H}\tilde{\boldsymbol{U}}+\mathcal{R}(\boldsymbol{U}_{n}^{q+1})=0$
(11)

where $q+1$ is the physical time step index and where $n$ is now the time
instance. The expression can then be linearised and a semi-implicit approach
can be used to solve for each time instance. An unsteady problem can then be
characterized by $K$ frequencies.

###  3.4 Non-Conformal Boundary Interface Treatment

The numerical discretization of complex rotorcraft geometries is a non-trivial
task. It requires the main and tail rotor blades to move in relative motion to
the fuselage. To permit this, the computational domain should be split into
separate sub-regions, thus introducing non-conformal boundary interface.
Suitable treatment of each artificial boundary interface is achieved through
the supermesh technique [26] implemented by developers Gori et al. [14] of
SU2. The supermesh acts an additional auxiliary grid between two non-conformal
interfaces and is of size $n-1$ dimensional elements where $n$ is the
dimension of the computational domain. The algorithm used to construct the
supermesh of two generic parent mesh is shown in Fig. 1 and subsequently
outlined:

  1. 1.

Identify the nodes of the control volume element $\mathcal{E}_{B}$ contained
within $\mathcal{E}_{A}$.

  2. 2.

Determine the points of intersection between the edges of $\mathcal{E}_{A}$
and $\mathcal{E}_{B}$.

  3. 3.

Decompose the overlap region into triangles.

  4. 4.

Calculate the overlapping area of $\mathcal{E}_{B}$ inside of
$\mathcal{E}_{A}$ and the weight $W_{i}$.

this procedure is repeated for each of the overlapping neighbouring elements
and its contribution to the numerical flux is assembled by
$\sum^{N_{sf}}_{q=1}{W_{q}}$. Where the number of supermesh faces that are
mapped to each control volume is denoted by $\mathcal{N}_{sf}$. The flux
balance across the supermesh interface is obtained through incorporating the
neighbouring cells either side of the boundary into Eq. 5, which can then be
re-written in the semi-discrete form in the cell center of each control
volume, $k$ as,

$\frac{\partial\boldsymbol{U}_{k}}{\partial
t}=-\frac{1}{\left|\Omega_{k}\right|}\left[\sum_{p=1}^{\mathcal{N}_{f}}\boldsymbol{F}(\boldsymbol{U}_{k},\,\boldsymbol{U}_{p})\Delta
S_{p}+\sum_{p=1}^{\mathcal{N}_{bf}}\Delta
S_{p}\left(\sum_{q=1}^{\mathcal{N}_{sf}}W_{q}\boldsymbol{F}(\boldsymbol{U}_{k},\,\boldsymbol{U}_{q})\right)\right]$
(12)

where the number of internal faces and boundary faces of each control volume
is represented by $\mathcal{N}_{f}$ and $\mathcal{N}_{bf}$. The area of the
face is represented by $\Delta S$. The vector of fluxes are represented as
$\boldsymbol{F}$. The conservative variables are once again represented by
$\boldsymbol{U}$ and the subscripts $p$ and $q$ denote the neighbouring and
boundary cells.

Figure 1: Schematic showing the steps for the construction of the supermesh.

###  3.5 Main Rotor Blade Kinematics

The method in which conventional main-rotor/tail-rotor helicopters generate
their vertical and propulsive forces, as well as moments to control the
attitude and position of the helicopter in three-dimensional space, will now
be introduced and its implementation within SU2 discussed. A schematic
depicting the blade motion is shown in Fig. 2. Unlike fixed-wing aircraft, the
rotor blades alone must provide both the lifting forces and the control
meaning the aerodynamics and dynamics of the blades are coupled so that it is
hard to understand one without the other. During hover and axial flight, the
velocity variation along the blade is azimuthally axisymmetric and depends
solely upon the relationship of $\omega r$. In forward flight, however, this
is not the case and a component of the free-stream velocity alters the blade
velocity such that it now depends upon its azimuthal position, $\psi$ as given
by,

$M_{\textrm{n}}(\psi)=M_{\text{tip}}\frac{r}{R}+M_{\infty}\sin\psi=M_{\text{tip}}\left(\frac{r}{R}+\mu\sin\psi\right)$
(13)

where $M_{\textrm{n}}$ represents the blade normal Mach number and
$M_{\infty}$ is the freestream Mach number. The advance ratio, $\mu$, can be
considered as the relationship of $M_{\infty}/M_{\text{tip}}$ . A repercussion
of local variations in velocity is that each blade is now required to be
individually controlled to eliminate rolling and pitching moments. With this,
each blade is typically connected to a rotor head by a set of hinges which
allows the blades to move independently. In reality, the rotor blades are also
highly elastic in nature, however, for simplification this work considers the
blades to be rigid with the ramifications of this being understood. Rigid
blades are then able to move with respect to the hinge positions. Each of the
hinges will now be introduced: a flap hinge allows the blade to move in the
plane containing the blade and the shaft; a pitch hinge allows the blade to
move around the quarter chord and spanwise axis; finally a lead-lag hinge
allows the blade to move in-plane forward or backwards. The flapping $\beta$,
lead-lag $\delta$, and pitching $\theta$ motion can then be described as a
function of the azimuthal position $\psi$ of the blades as they rotate,

$\displaystyle\begin{aligned} \psi&=\omega t\\\
\beta(\psi)&=\beta_{0}-\beta_{1s}\sin(\psi)-\beta_{1c}\cos(\psi)-\beta_{2s}\sin(2\psi)-\beta_{2c}\cos(2\psi)-...\\\
\delta(\psi)&=\delta_{0}-\delta_{1s}\sin(\psi)-\delta_{1c}\cos(\psi)-\delta_{2s}\sin(2\psi)-\delta_{2c}\cos(2\psi)-...\\\
\theta(\psi)&=\theta_{0}-\theta_{1s}\sin(\psi)-\theta_{1c}\cos(\psi)-\theta_{2s}\sin(2\psi)-\theta_{2c}\cos(2\psi)-...\end{aligned}$
(14)

where $\beta_{0}$ denotes the coning angle and $\theta_{0}$ denotes the
collective pitch applied equally to all the blades. The lateral and
longitudinal components of the motion are represented by subscripts $c$ and
$s$ respectively.

A rotation matrix then transforms the blade Cartesian coordinates,
$\boldsymbol{x}$, from the fixed frame of reference into the hub frame of
reference. The hub reference system is then transformed into each of the blade
reference systems and the blade motion laws described in Equation 14 are
applied. If we now consider the definition of the reference system where; the
rotation occurs around the $z$-axes, the flapping occurs around the $y$-axes,
the lead-lag motion around the $z$-axes, and the pitching occurs around the
$x$-axes as depicted in Fig. 2 the following transformation matrices can be
introduced to prescribe the blade motion,

$\displaystyle\begin{aligned}
C_{\textrm{rot}}&=\begin{pmatrix}\cos\psi&-\sin\psi&0\\\
\sin\psi&\cos\psi&0\\\ 0&0&1\end{pmatrix},\\\ C_{\textrm{lead-
lag}}&=\begin{pmatrix}\cos\delta&-\sin\delta&0\\\ \sin\delta&\cos\delta&0\\\
0&0&1\end{pmatrix},\end{aligned}$ $\displaystyle\begin{aligned}
C_{\textrm{flap}}&=\begin{pmatrix}\cos\beta&0&-\sin\beta\\\ 0&1&0\\\
\sin\beta&0&\cos\beta\end{pmatrix}\\\
C_{\textrm{pitch}}&=\begin{pmatrix}1&0&0\\\ 0&\cos\theta&-\sin\theta\\\
0&\sin\theta&\cos\theta\end{pmatrix}\end{aligned}$ (15)

The transformation from the fixed $\rightarrow$ hub reference system using
$C_{\textrm{rot}}$ is implemented by a rigid rotation of the entire grid as
the azimuthal rotation of the blades far exceeds the blade deflections. The
rate of rotation is currently constrained and is maintained constant such
that,

$C_{\textrm{rot}}=\begin{pmatrix}\cos\psi&-\sin\psi&0\\\
\sin\psi&\cos\psi&0\\\ 0&0&1\end{pmatrix}$ (16)

The transformation from the hub $\rightarrow$ blade reference system can then
be combined and applied to the blade surface. This is achieved through
deforming the mesh at each physical time step and is required as each blade
has its own independent motion. The subsequent transformation matrix,
$C_{\textrm{flap-leadlag-pitch}}$ can then be introduced,

$C_{\textrm{flap-leadlag-
pitch}}=\begin{pmatrix}\cos\beta\cos\delta&\sin\beta\sin\theta-\cos\beta\cos\theta\sin\delta&\cos\theta\sin\beta+\cos\beta\sin\delta\sin\theta\\\
\sin\delta&\cos\delta\cos\theta&-\cos\delta\sin\theta\\\
-\cos\delta\sin\beta&\cos\beta\sin\theta+\cos\theta\sin\beta\sin\delta&\cos\beta\cos\theta-\sin\beta\sin\delta\sin\theta\end{pmatrix}$
(17)

The motion can then be applied to the coordinates,
$\boldsymbol{x}=\\{{x,y,z}\\}^{T}$ at an arbitrary point, $i$ on the blade as,

$\boldsymbol{x}_{i}=C_{\textrm{rot}}C_{\textrm{flap-leadlag-
pitch}}(\boldsymbol{x}_{i}-\boldsymbol{x}_{\textrm{hinge}})$ (18)

where $\boldsymbol{x}_{\textrm{hinge}}$ is the position of the hinge about
which the flapping, lead-lag and pitching occurs. Furthermore the grid
velocity is updated using a second order finite difference scheme,

$\boldsymbol{u}_{\Omega_{i}}=\frac{3\boldsymbol{x}_{i}^{n+1}-4\boldsymbol{x}_{i}^{n}+\boldsymbol{x}_{i}^{n-1}}{2\Delta
t}$ (19)

Figure 2: Schematic of the main rotor blade motion and respective rotations.

###  3.6 Radial Basis Function Mesh Deformation

To allow the rotor blades to move independently requires the mesh to deform at
each physical time step. Rendall and Allen, to the authors’ knowledge, were
the first and only group to demonstrated how Radial Basis Functions (RBF)
could be used to account for the blade motion [27]. The outlook of using RBF
mesh deformation for rotor blade motion was promising however one major
drawback was the high computational cost for large meshes. Despite this, RBF
mesh deformation methods are robust and preserve high-quality mesh even during
large deformations. They also have the unique property that they do not
require the grid connectivity meaning that even for three-dimensional problems
they are relatively simple to implement.

The term radial basis function refers to a series of functions whose values
depends on their distance to a supporting position. In the most general of
forms, radial basis functions can be written as,
$\phi(\textbf{r},\textbf{r}_{i})=\phi\left(\|\textbf{r}-\textbf{r}_{i}\|\right)$,
where the distance corresponds to the radial basis centre, $\textbf{r}_{i}$.
This distance is frequently taken as the Euclidean distance, meaning it
becomes the spatial distance between two nodes.

An interpolation function, $f(\textbf{r})$ can be introduced as a method for
describing the displacement of a set of nodes in space and can be approximated
by a weighted sum of basis functions. However, the interpolation relies on the
weight coefficients of the basis points, $\alpha$. The interpolation takes the
form

$f(\textbf{r})=\sum_{i=1}^{N}\alpha_{i}\phi\left(\left\|\textbf{r}-\textbf{r}_{i}\right\|\right)$
(20)

The primary expense of RBF mesh deformation is associated to the solving the
linear system to obtain the weight coefficients. To address the concerns of
the computation cost of RBF for rotor blades, multi-level greedy surface point
selection algorithms [28] and volume point reduction methods [29] for large
scale problems are introduced which greatly improve the efficiency. These
numerically enhanced methods were implemented by Ref. [30] in SU2 and are now
used within this work.

##  4 Results

###  4.1 Hover

The benchmark Caradonna and Tung experimental tests of a model rotor in hover
[31] are used for hover performance validation. The experiment was carried out
in the U.S. Army Aeromechanics Laboratory’s hover test facility which included
special ducting designed to eliminate room recirculation. The rotor employed
two manually adjustable cantilever-mounted blades attached to a drive shaft at
the top of a column test stand. The blades were untwisted and untapered with
an aspect ratio of 6. The profile of the blades used a symmetric NACA0012
airfoil. The rotor radius was $1.143$ m and the blade chord was $0.191$ m. The
data recorded during the tests consisted of pressure measurements along the
chord of the blade at various radial positions. The tests were conducted at a
variety of collective pitch angles and tip Mach numbers. The conditions used
for the validation of the blade loading predictions are detailed in Table 1.

Table 1: Hover test conditions taken from the Caradonna-Tung experiment [31]. | Collective Pitch
---
(∘)
| Rotational Speed
---
(RPM)
| Tip Mach No.
---
( )
8 | 1250 | 0.439

The problem is transformed into a rotating frame of reference for
computational efficiency and the RANS equations are solved to first assess the
influence of the radial position on the surface pressure distribution.
Computed and experimental surface pressure distributions are compared at
select radial positions $r/R=0.68\%$, $r/R=0.80\%$, $r/R=0.89\%$ and
$r/R=0.96\%$ as represented by Fig. 3.

The standard one-equation SA turbulence model is used with the RANS equations.
The SA turbulence variable is convected using a first-order scalar upwind
method. The convective fluxes are computed using the Roe scheme and second-
order accuracy is achieved using the MUSCL scheme. Spatial gradients for the
viscous fluxes and second order reconstruction are approximated using the
Green-Guass method. The flow solution was considered converged when there was
a reduction of 6 orders of magnitude on the density residual.

The predicted results in Fig. 4 are in close agreement with the measured data
at all radial positions. Both the upper and lower surface pressure profiles
follow the measured data. The suction peak on the upper surface is captured
well. There are however discrepancies at the trailing edge and this is
expected to be caused by modelling the trailing edge as flat which was done to
aid with the mesh generation process and convergence of the solver. The
results reaffirm the calculations from Ref. [32].

An assessment of the most influential parameters during hover simulations are
displayed in Fig. 5. Each of the parameters are compared against the measured
data at $r/R=0.80\%$. The spatial mesh resolution is assessed in Fig. 5a. Four
levels of mesh refinement are used ranging from coarse $\rightarrow$ very
fine. The ‘coarse’ mesh has 1.56 million volume elements and 22 thousand
surface elements. The ‘medium’ mesh has 2.62 million volume elements and 41
thousand surface elements. The ‘fine’ mesh has 5.19 million volume elements
and 80 thousand surface elements. The ‘very fine’ mesh has 11.09 million
volume elements and 157 thousand surface elements. The coarse mesh exhibits
instabilities in the pressure profile at the leading-edge. The medium mesh
reduces this and shows a smoother profile however the suction peak remains
slightly below the measured data. The fine mesh improves the suction peak
compared to the medium mesh. There appear to be no discernible differences
between the fine and very fine mesh suggesting the mesh has reached
convergence.

It is frequently being desirable to model the rotor as a time-accurate problem
where the blades are physically moving. In this scenario, the number of
revolutions the blades have completed when the flow is initialized from the
freestream is important and is shown in Fig. 5c. It shows that the time for
convergence of the solution is at least 3 full rotor revolutions.

The influence of the choice of steady or time-accurate modelling on the
pressure profiles is shown in Fig. 5e. All methods show good agreement with
the measured data. The rotating frame method shows the closest agreement.
There is a very slight difference on the upper surface for the time-accurate
sliding mesh, rigid motion and harmonic balance methods. Furthermore, all of
the time-accurate methods appear overlapping.

Finally the flow regime is evaluated in Fig. 5g. It depicts the influence of
inviscid, laminar and turbulent flow on the pressure distribution. As excepted
there appears to be significant difference between inviscid flow and laminar
and turbulent flow regimes.

A contour map of the pressure coefficient on both the suction and pressure
sides of the blades for time-accurate problems is shown in Fig. 6. The
approach using unsteady rigid motion is displayed at $10^{\circ}$ azimuth
increments in Fig. 6a. The approach using the harmonic balance is displayed at
3 time instances with the input frequencies of
$\boldsymbol{\omega}=(0,\,\pm\omega_{1})$ in Fig. 6b. Both sets of results
show that the pressure coefficient during hover is axisymmetric about the out-
of-plane axis.

The iso-surface of the Q criterion visualizing the near-field wake and blade
tip vortices during hover are displayed in Fig. 7. The wake exhibits no
interaction effects due to the preceding blade tip vortices passing above
previous blade tip vortices and due to the isolated blades being modelled
without a hub or test stand.

###  4.2 Forward Flight

The extensive report and data survey from Cross and Watts on tip aerodynamics
and acoustics [33] is used for forward flight performance validation. The
report describes the Tip Aerodynamic and Acoustics Tests (TAAT) carried out at
the NASA Ames Research Center. The TAAT used a highly instrumented AH-1G cobra
helicopter and measured the rotor airloads at multiple radial locations. The
AH-1G is a two-person single-engine helicopter. It has a two-bladed teetering
rotor. The blades are untapered and have a linear twist of $-10^{\circ}$ from
root to tip. The rotor radius is 6.71 m and the blade chord is 0.686 m
resulting in an aspect ratio of 9.8. The profile of the blades use a symmetric
highly modified 540 airfoil section stemming from the NACA0012 family. In the
aerodynamics phase of testing several specific thrust coefficients, tip Mach
numbers, and advance ratios were flown to study the correlation between these
parameters and the pressure distributions.

The flight tests chosen for the present validation work are detailed in Table
2. Two subsets of the flight tests are used. The first being a low-speed test
which correlates to counter number 2157 from the database and the second being
a high-speed test which correlates to counter number 2152 from the database.
During the 150 km/h low-speed condition the rotor is operated at a tip Mach
number of 0.65 and an advance ratio of 0.19. The collective pitch of
$11.7^{\circ}$ is set to be congenial with the measured rotor thrust of
$\textrm{C}_{\textrm{T}}=0.00464$ and the cyclic pitch is trimmed to eliminate
pitching and rolling moments. The sine and cosine cyclic pitch angles are
respectively set to $1.7^{\circ}$ and $-5.5^{\circ}$. The sine and cosine
cyclic flap angles are respectively set to $-0.15^{\circ}$ and $2.13^{\circ}$.
During the 290 km/h high-speed condition the rotor is operated at a tip Mach
number of 0.64 and an advance ratio of 0.24. The collective pitch of
$18.0^{\circ}$ is set to be congenial with the measured rotor thrust of
$\textrm{C}_{\textrm{T}}=0.00464$ and the cyclic pitch is trimmed to eliminate
pitching and rolling moments. The sine and cosine cyclic pitch angles are
respectively set to $3.6^{\circ}$ and $-11.8^{\circ}$. The sine and cosine
cyclic flap angles are respectively set to $1.11^{\circ}$ and $2.13^{\circ}$.

Table 2: Forward flight conditions taken from the TAAT data survey [33]. Variable |  | Low Speed
---
Counter 2157
| High Speed
---
Counter 2152
$\textrm{V}_{\infty}$ | 150 km/h | 290 km/h
$\textrm{M}_{\infty}$ | 0.12 | 0.24
$\textrm{M}_{\textrm{T}}$ | 0.65 | 0.64
$\textrm{Re}_{\textrm{T}}$ | $9.7\times 10^{6}$ | $10.2\times 10^{6}$
$\mu$ | 0.19 | 0.38
$\theta_{0}$ | $11.7^{\circ}$ | $18.0^{\circ}$
$\theta_{0.75}$ | $5.7^{\circ}$ | $12.0^{\circ}$
$\theta_{s}$ | $1.7^{\circ}$ | $3.6^{\circ}$
$\theta_{c}$ | $-5.5^{\circ}$ | $-11.8^{\circ}$
$\beta_{0}$ | $2.75^{\circ}$ | $2.75^{\circ}$
$\beta_{c}$ | $2.13^{\circ}$ | $2.13^{\circ}$
$\beta_{s}$ | $-0.15^{\circ}$ | $1.11^{\circ}$
$\textrm{C}_{\textrm{T}}$ | 0.00464 | 0.00474

The spatial discretization is achieved using a single-zone mesh of the two
main rotor blades. The outer far-field is placed 5 radii away from the blades.
A mixed-element grid composed of 30.25 million elements and 760 thousand
vertices is used. Each blade has 105 thousand surface elements. A density
region surrounding the rotor blades is used to sufficiently resolve the near-
field wake. Elements within the density region have a maximum element size of
$0.0075\,x/c$.

The standard one-equation SA turbulence model is used with the RANS equations.
The SA turbulence variable is convected using a first-order scalar upwind
method. The convective fluxes are computed using the Roe scheme and second-
order accuracy is achieved using the MUSCL scheme. Spatial gradients are
approximated using the Green-Guass numerical method. The flow solution was
considered converged when there was a reduction of 6 orders of magnitude on
the density residual. Per rotor revolution, there are 360 physical time-steps
equating to $1^{\circ}$ azimuth increments. At each physical time-step, there
are a maximum of 30 internal pseudo time-steps. Each simulation was run for a
total of 5 rotor revolutions corresponding to 1800 physical time steps.

The blade motion is prescribed using the implementation introduced in Section
3.5 with the values outlined in Table 2. The mesh is updated at each physical
time-step using the RBF mesh deformation strategy. Multi-level surface point
reduction and volume point reduction algorithms are used to improve the
efficiency. Four levels are employed with the surface point selection error
reduction rates and efficiency being shown in Fig. 8. It shows that the multi-
level method for updating the blade position is highly effective in minimizing
the computational cost while simultaneously reducing the surface error. The
control points and associated normalized surface error of the computed
displacement at each level are shown in Fig. 9. The minimum orthogonality
angle of the grid prior to deformation was $40.979^{\circ}$. The minimum
orthogonality angle of the grid post 1800 deformations is $40.997^{\circ}$.
This indicates there is no detrimental effect on the grid quality from the
continuous deformation.

The pressure coefficient distributions at selected advancing and retreating
azimuth positions as well as radial positions are obtained from the data
survey and are used for validation. The low-speed condition uses data at
$r/R=0.60$ for the advancing side of the rotor at azimuth angles
$\psi=30^{\circ}$, $\psi=90^{\circ}$, and $\psi=180^{\circ}$ with the results
being shown in Fig. 10. Data at $r/R=0.91$ on the retreating side of the rotor
at azimuth angles $\psi=270^{\circ}$, $\psi=285^{\circ}$, and
$\psi=300^{\circ}$ are also used and the results are shown in Fig. 11. In
general, the trends are relatively well captured when compared against the
measured data. On the advancing side of the rotor the on the upper surface the
suction peaks are marginally over-exaggerated. On the retreating side of the
rotor this is improved slightly.

The high-speed conditions uses data at $r/R=0.86$ for the advancing side of
the rotor at azimuth angles $\psi=70^{\circ}$, $\psi=90^{\circ}$, and
$\psi=110^{\circ}$ with the results being shown in Fig. 12. Data at $r/R=0.96$
on the retreating side of the rotor at azimuth angles $\psi=250^{\circ}$,
$\psi=270^{\circ}$, and $\psi=290^{\circ}$ are also used and the results are
shown in Fig 15. The higher flight speed initiates the presence of shocks on
the advancing side of the blade and dynamic stall on the retreating side of
the blade which makes replicating the measured data more challenging than the
low-speed condition. Consequently there are significantly different pressure
coefficient distributions. The suction peaks on the upper surface of the blade
on the advancing side of the rotor are represented well. The shock positions
however show discrepancies when compared to the measured data. The pressure
coefficient distributions on the retreating side of the blade close to the tip
appear to be better represented.

A contour map of the pressure coefficient on both the suction and pressure
sides of the blades for both the low- and high-speed conditions are shown in
Fig. 16. Both sets of results show that the pressure coefficient during
forward flight is no longer axisymmetric about the out-of-plane axis. To
achieve the higher flight speed an increase in the collective pitch is
required. This however results in a low speed stalled region on the retreating
side of the rotor near the blade root where there is a high angle-of-attack.
This effect can be seen in Fig. 16b which is not present in Fig. 16a.

The iso-surface of the Q-criterion visualizing the near-field wake and blade
tip vortices is shown in Fig. 17. It displays the different flow field
behaviour in low-speed and high-speed forward flight. In low-speed forward
flight, the vortex upstream of the preceding blade passes below the advancing
blade close to the tip. Due to the low flight speed, the vortex downstream
remains close to the rotor blades. In high-speed forward flight the vortex
upstream of the preceding blade passes below the advancing blade closer to the
root of the blade. With the flight speed being much higher the vortex
downstream is no longer in the vicinity of the rotor blades. Due to both cases
being straight and level flight there are no blade-vortex interactions
present. It is of note that even with a density region and such a large grid
that it is challenging to reduce artificial numerical dissipation of the blade
tip vortices.

##  5 Conclusion

This work provides a unique approach for the modelling of rotorcraft
aerodynamics using the open-source toolkit SU2. The performance prediction
capabilities are demonstrated on two test cases for preliminary validation.
The first test case of the Caradonna-Tung model rotor in hover allowed for the
assessment of the pressure coefficient distribution along four different
radial positions. Models with varying fidelity ranging from steady-state to
fully time-accurate were assessed alongside significant influential
parameters. The second test case of the AH-1G rotor in forward flight required
the introduction of the blade motion. Low and high speed forward flight tests
were modelled. The low-speed test condition was found to be easier to simulate
and the performance predictions are in good agreement. The high-speed test
condition was more challenging due to the presence of a shock near the blade
tip on the advancing side and dynamic stall on the retreating side of the
blade near the root. Despite the results of this work being complementary to
the measured data, it is not yet sufficient to suggest this is satisfactory
for complete validation. It is import that more tests are carried out with
different configurations and test conditions. Moving forward, likely areas of
research will include; rotor-fuselage interaction effects, blade design
optimization, rotor acoustics, mesh adaptation, and the extension of rigid
blades to elastic blades.

## Figures

Figure 3: Blade stations where the static pressure measurements were taken
during Ref. [31]. Select radial positions $r/R=0.68\%$, $r/R=0.80\%$,
$r/R=0.89\%$ and $r/R=0.96\%$ are used for the assessment of numerical
predictions.

(a) $r/R=0.68\%$.  (b) $r/R=0.80\%$.

(c) $r/R=0.89\%$.  (d) $r/R=0.96\%$.

Figure 4: Pressure distributions at different radial positions computed using
the rotating frame method during hover. Predictions compared against measured
data taken from Caradonna-Tung [31].

(a) Grid sensitivity.  (b) Region around the leading-edge.

(c) Rotor revolutions.  (d) Region around the leading-edge.

(e) Steady-state vs. time accurate approaches.  (f) Region around the leading-
edge.

(g) Inviscid, laminar and turbulent-SA.  (h) Region around the leading-edge.

Figure 5: Comparison of different influential parameters on the performance
prediction during hover. Pressure distributions shown at $r/R=0.80\%$ and
compared against measured data taken from Caradonna-Tung [31].

(a) Unsteady rigid motion at $10^{\circ}$ azimuth increments.

(b) Harmonic balance at 3 time instances with the input frequencies of
$\boldsymbol{\omega}=(0,\,\pm\omega_{1})$.

Figure 6: Contour maps of the pressure coefficient on the upper and lower
blade surface during hover. Comparing the reduced order harmonic balance
method against a fully time-accurate method. Figure 7: Iso-surface of the
Q-criterion visualizing the near-field wake and blade tip vortices as they are
convected downstream in the wake of the rotor during hover. Figure 8: Multi-
level greedy surface point selection error reduction rates for a single
deformation on the AH-1G rotor during forward flight.

(a) Level 1 - avg. 6 control points.  (b) Level 2 - avg. 21 control points.

(c) Level 3 - avg. 42 control points.  (d) Level 4 - avg. 136 control points.

Figure 9: Evolution of the multi-level greedy surface point selection
algorithm used for the RBF mesh deformation. Displaying the control points as
black dots and the contours of the normalized surface error displacement.
Where the normalized error is computed as the difference between the computed
displacement and the actual displacement.

(a) $\psi=30^{\circ}$.  (b) $\psi=90^{\circ}$.  (c) $\psi=180^{\circ}$.

Figure 10: Pressure distributions at $r/R=0.60$ on the advancing side of the
rotor during low-speed forward flight. Predictions compared against the
measured data taken from the TAAT [33].

(a) $\psi=270^{\circ}$.  (b) $\psi=285^{\circ}$.  (c) $\psi=300^{\circ}$.

Figure 11: Pressure distributions at $r/R=0.91$ on the retreating side of the
rotor during low-speed forward flight. Predictions compared against the
measured data taken from the TAAT [33].

(a) $\psi=70^{\circ}$.  (b) $\psi=90^{\circ}$.  (c) $\psi=110^{\circ}$.

Figure 12: Pressure distributions at $r/R=0.86$ on the advancing side of the
rotor during high-speed forward flight. Predictions compared against the
measured data taken from the TAAT [33].

(a) $\psi=250^{\circ}$.  (b) $\psi=270^{\circ}$.  (c) $\psi=290^{\circ}$.

Figure 13: Pressure distributions at $r/R=0.96$ on the advancing side of the
rotor during high-speed forward flight. Predictions compared against the
measured data taken from the TAAT [33].

(a) $r/R=0.75$.  (b) $r/R=0.86$.

Figure 14: Predictions compared against the measured data taken from the TAAT
[33].

(a) $r/R=0.75$.  (b) $r/R=0.86$.

Figure 15: Predictions compared against the measured data taken from the TAAT
[33].

(a) Low speed forward flight.  (b) High speed forward flight.

Figure 16: Contour maps of the pressure coefficient on the upper blade surface
during forward flight. Displaying the both low-speed and high-speed test
conditions. Figure 17: Iso-surface of the Q-criterion visualizing the near-
field wake and blade tip vortices. Displaying the different flow field
behaviour in low speed and high speed forward flight.

###### List of Figures

  1. 1 Supermesh schematic
  2. 2 Main rotor schematic
  3. 3 Blade static pressure measurements
  4. 4 $C_{p}$ distributions at various radial positions during hover
  5. 5 Influential parameters on the performance prediction during hover
  6. 6 $C_{p}$ distribution on the upper and lower blade surface during hover
  7. 7 Q-criterion iso-surface during hover
  8. 8 Multi-level convergence history on the AH-1G rotor
  9. 9 Selected control points on the AH-1G rotor
  10. 10 Low-speed performance prediction on the advancing side
  11. 11 Low-speed performance prediction on the retreating side
  12. 12 High-speed performance prediction on the advancing side
  13. 13 High-speed performance prediction on the retreating side
  14. 14 Predictions compared against the measured data taken from the TAAT [33].
  15. 15 Predictions compared against the measured data taken from the TAAT [33].
  16. 16 $C_{p}$ distribution on the blade surface in forward flight
  17. 17 Q-criterion iso-surface during forward flight

###### List of Tables

  1. 1 Hover test conditions taken from the Caradonna-Tung experiment
  2. 2 Forward flight conditions taken from the TAAT data survey

## References

## References

  * [1] Aviation Accident Reports - National Transportation Safety Board, https://www.ntsb.gov/investigations/AccidentReports/Pages/aviation.aspx (Accessed: May, 2020).
  * [2] R. C. Strawn, F. X. Caradonna, E. P. Duque, 30 years of rotorcraft computational fluid dynamics research and development, Journal of the American Helicopter Society 51 (1) (2006) 5–21, https://doi.org/10.4050/1.3092875.
  * [3] D. Jespersen, T. Pulliam, P. Buning, Recent enhancements to overflow, in: 35th Aerospace Sciences Meeting and Exhibit, Reno, Navada, U.S.A, 06-09 January, 1997, p. 644, https://doi.org/10.2514/6.1997-644.
  * [4] W. K. Anderson, D. L. Bonhaus, An implicit upwind algorithm for computing turbulent flows on unstructured grids, Computers & Fluids 23 (1) (1994) 1–21, https://doi.org/10.1016/0045-7930(94)90023-X.
  * [5] R. Steijl, G. Barakos, K. Badcock, A framework for cfd analysis of helicopter rotors in hover and forward flight, International journal for numerical methods in fluids 51 (8) (2006) 819–847, https://doi.org/10.1002/fld.1086.
  * [6] M. Biava, L. Vigevano, Simulation of a complete helicopter: A cfd approach to the study of interference effects, Aerospace Science and Technology 19 (1) (2012) 37–49, https://doi.org/10.1016/j.ast.2011.08.006.
  * [7] A. Antoniadis, D. Drikakis, B. Zhong, G. Barakos, R. Steijl, M. Biava, L. Vigevano, A. Brocklehurst, O. Boelens, M. Dietz, et al., Assessment of cfd methods against experimental flow measurements for helicopter flows, Aerospace Science and Technology 19 (1) (2012) 86–100, https://doi.org/10.1016/j.ast.2011.09.003.
  * [8] G. R. Srinivasan, J. Baeder, S. Obayashi, W. McCroskey, Flowfield of a lifting rotor in hover-a navier-stokes simulation, AIAA journal 30 (10) (1992) 2371–2378, https://doi.org/10.2514/3.11236.
  * [9] G. Srinivasan, J. Baeder, Turns: A free-wake eule/navier-stokes numerical method for helicopter rotors, AIAA journal 31 (5) (1993) 959–962, https://doi.org/10.2514/3.49036.
  * [10] M. Gazaix, A. Jolles, M. Lazareff, The elsa object-oriented computational tool for industrial applications, in: 23rd Congress of ICAS, Toronto, Canada, 08-13 September, 2002, p. 220, Corpus ID: 56582723.
  * [11] J. Raddatz, J. K. Fassbender, Block structured navier-stokes solver FLOWer, MEGAFLOW - Numerical Flow Simulation for Aircraft Design. Notes on Numerical Fluid Mechanics and Multidisciplinary Design (NNFM) Edition, Vol. 89, Springer, Berlin, Heidelberg, 2005, https://doi.org/10.1007/3-540-32382-1.
  * [12] W. L. Oberkampf, F. G. Blottner, Issues in computational fluid dynamics code verification and validation, AIAA journal 36 (5) (1998) 687–695, https://doi.org/10.2514/2.456.
  * [13] T. D. Economon, F. Palacios, S. R. Copeland, T. W. Lukaczyk, J. J. Alonso, Su2: An open-source suite for multiphysics simulation and design, Aiaa Journal 54 (3) (2016) 828–846, https://doi.org/10.2514/1.J053813.
  * [14] G. Gori, E. van der Weide, A. Guardone, On conservation in compressible flow simulations using sliding mesh coupling, in: VII International Conference on Computational Methods for Coupled Problems in Science and Engineering conference, Rhodes Island, Greece, June 12-14, 2020.
  * [15] D. C. Wilcox, et al., Turbulence modeling for CFD, 2nd Edition, DCW industries La Canada, CA, 1998, ISBN 13: 9780963605153.
  * [16] W. Sutherland, The viscosity of gases and molecular force, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 36 (223) (1893) 507–531, https://doi.org/10.1080/14786449308620508.
  * [17] P. Spalart, S. Allmaras, A one-equation turbulence model for aerodynamic flows, in: 30th aerospace sciences meeting and exhibit, Reno, Nevada, U.S.A, 06-09 January, 1992, p. 439, https://doi.org/10.2514/6.1992-439.
  * [18] T. D. Economon, F. Palacios, J. J. Alonso, A viscous continuous adjoint approach for the design of rotating engineering applications, in: 21st AIAA computational fluid dynamics conference, San Diego, California, U.S.A, 24-27 June, 2013, p. 2580, https://doi.org/10.2514/6.2013-2580.
  * [19] H. K. Versteeg, W. Malalasekera, An introduction to computational fluid dynamics: the finite volume method, 2nd Edition, Pearson education, 2007, ISBN: 978-0-13-127498-3.
  * [20] A. Jameson, W. Schmidt, E. Turkel, Numerical solution of the euler equations by finite volume methods using runge kutta time stepping schemes, in: 14th fluid and plasma dynamics conference, Palo Alto, California, U.S.A, 23-25 June, 1981, p. 1259, https://doi.org/10.2514/6.1981-1259.
  * [21] P. L. Roe, Approximate riemann solvers, parameter vectors, and difference schemes, Journal of computational physics 43 (2) (1981) 357–372, https://doi.org/10.1016/0021-9991(81)90128-5.
  * [22] B. Van Leer, Towards the ultimate conservative difference scheme. v. a second-order sequel to godunov’s method, Journal of computational Physics 32 (1) (1979) 101–136, https://doi.org/10.1016/0021-9991(79)90145-1.
  * [23] A. Jameson, Time dependent calculations using multigrid, with applications to unsteady flows past airfoils and wings, in: 10th Computational Fluid Dynamics conference, Honolulu, Hawaii, U.S.A., p. 1596.
  * [24] K. C. Hall, J. P. Thomas, W. S. Clark, Computation of unsteady nonlinear flows in cascades using a harmonic balance technique, AIAA journal 40 (5) (2002) 879–886, https://doi.org/10.2514/2.1754.
  * [25] A. Rubino, M. Pini, P. Colonna, T. Albring, S. Nimmagadda, T. Economon, J. Alonso, Adjoint-based fluid dynamic design optimization in quasi-periodic unsteady flow problems using a harmonic balance method, Journal of Computational Physics 372 (2018) 220–235, https://doi.org/10.1016/j.jcp.2018.06.023.
  * [26] E. Rinaldi, P. Colonna, R. Pecnik, Flux-conserving treatment of non-conformal interfaces for finite-volume discretization of conservation laws, Computers & Fluids 120 (2015) 126–139, https://doi.org/10.1016/j.compfluid.2015.07.017.
  * [27] T. Rendall, C. Allen, Parallel efficient mesh motion using radial basis functions with application to multi-bladed rotors, International journal for numerical methods in engineering 81 (1) (2010) 89–105, https://doi.org/10.1002/nme.2678.
  * [28] G. Wang, H. H. Mian, Z.-Y. Ye, J.-D. Lee, Improved point selection method for hybrid-unstructured mesh deformation using radial basis functions, AIAA Journal 53 (4) (2015) 1016–1025, https://doi.org/10.2514/1.J053304.
  * [29] L. Xie, H. Liu, Efficient mesh motion using radial basis functions with volume grid points reduction algorithm, Journal of Computational Physics 348 (2017) 401–415, https://doi.org/10.1016/j.jcp.2017.07.042.
  * [30] M. Morelli, T. Bellosta, A. Guardone, Efficient radial basis function mesh deformation methods for aircraft icing, in: In the 7th European Seminar on Computing, Pilsen, Czech Republic, June 8-12, 2020.
  * [31] F. X. Caradonna, C. Tung, Experimental and analytical studies of a model helicopter rotor in hover, in: Presented at the 6th European Rotorcraft and Powered Lift Aircraft Forum, Bristol, England, September 16-19, 1980.
  * [32] F. Palacios, T. D. Economon, A. Aranake, S. R. Copeland, A. K. Lonkar, T. W. Lukaczyk, D. E. Manosalvas, K. R. Naik, S. Padron, B. Tracey, et al., Stanford university unstructured (SU2): Analysis and design technology for turbulent flows, in: In 52nd Aerospace Sciences Meeting, National Harbor, Maryland, January 13-17, 2014, p. 0243, https://doi.org/10.2514/6.2014-0243.
  * [33] J. L. Cross, M. E. Watts, Tip aerodynamics and acoustics test: a report and data survey, National Aeronautics and Space Administration, Scientific and Technical Information Division, Reference Publication No. 1179, Ames Research Center, Moffett Field, California, December, 1988.


# Why (and When and How) Contrastive Divergence Works

Ian E. Fellows
Fellows Statistics
http://www.fellstat.com
<EMAIL_ADDRESS>
###### Abstract

Contrastive divergence (CD) is a promising method of inference in high
dimensional distributions with intractable normalizing constants, however, the
theoretical foundations justifying its use are somewhat weak. This document
proposes a framework for understanding CD inference, including how and when it
works. It provides multiple justifications for the CD moment conditions,
including framing them as a variational approximation. Algorithms for
performing inference are discussed and are applied to social network data
using an exponential-family random graph models (ERGM). The framework also
provides guidance about how to construct MCMC kernels providing good CD
inference, which turn out to be quite different from those used typically to
provide fast global mixing.

##  1 Introduction to Contrastive Divergence

Performing inference in the high dimensional data setting presents unique
computational challenges as traditional inference methods such as maximum
likelihood quickly become intractable. An alternative inference procedure is
contrastive divergence (CD) Hinton (2002), which has the advantage of being
relatively easy to compute stochastically resulting in it becoming a widely
used technique in the context of Boltzmann machines (Fischer and Igel, 2014).
This paper will focus on providing a theoretical foundation for CD, with
particular attention devoted to the case of exponential-family models.

Let $Y$ be a multi-dimensional random variable with realization $y$, $p(Y)$ be
the true distribution of the data generating process, and $q(Y)$ be a family
of of distributions that we are attempting to fit to the data. For example,
$q$ may be defined as being from an exponential-family

$q(Y)=\frac{1}{z(\eta)}e^{\eta g(Y)+o(Y)}~{}~{}~{}y\in{\cal Y},$

where $\eta$ is a vector of parameters, i$g$ is a vector of sufficient
statistics, $o$ is an offset statistic, and $z$ is the normalizing constant
defined as $z(\eta)=\sum_{y^{\prime}in{\cal Y}}e^{\eta
g(y^{\prime})+o(y^{\prime})}$Barndorff-Nielsen (1978). As the dimensionality
of $y$ increases, the sum in $z$ quickly becomes intractable to compute
exactly (except in special cases), thus the actual numeric value of $q$ is
difficult to evaluate. For every possible $\eta$, the distribution may also be
parameterized by a unique and equivalent set of mean value parameters
$\mu_{\eta}=E_{q}(g(Y))$Barndorff-Nielsen (1978).

Suppose we have one observation $y$ drawn from the true distribution $p$. It
is well known that maximum likelihood inference is equivalent to minimizing
the Kullbeck-Liebler (KL) divergence between $p$ and $q$ with the KL
divergence is defined as $\rm{KL}(p||q)=\sum_{y}\log(\frac{p(y)}{q(y)})p(y)$.
Minimizing the divergence yields $\rm{KL}(p||q)\propto-
E_{p}(log(q(Y)))\approx-\log(q(y)),$ which is the sample negative log-
likelihood.

Unfortunately, in many cases it is impossible to evaluate the log-likelihood
directly because it contains an intractable normalizing constant. As a
solution to this problem CD introduces a related distribution
$T_{q}(Y|y^{\prime}),~{}y^{\prime}\in{\cal Y})$, which is defined to be a
valid Gibbs transition kernel for $q$ (i.e.
$q(y)=\sum_{y^{(0)}}T_{q}(y|y^{(0)})q(y^{(0)})=E_{q}(T_{q}(y|Y))$. The
contrastive divergence objective function is defined as

$\rm{KL}(p||q)-\rm{KL}(T_{q}(Y|y)||q).$ (1)

The first term in the objective function is identical to that of maximum
likelihood inference. The second is the deviation between the Markov chain
distribution starting at the observed data vector $y$ and $q$. Typically
augmented objective functions of this form are justified by observing that
when $p=q$ the objective function is minimized, or nearly minimized (Lyu,
2011). Clearly this is true with the first term of the equation, but it is not
necessarily true for the second. When $p=q$ the objective function becomes
$-\rm{KL}(T_{p}(Y|y^{(0)}=y)||p)$, implying that the function reaches its
minimum at $p$ if $T$ and $q$ are maximally divergent at that point in the
neighborhood of $p$. The second term can be made to be arbitrarily close to 0
by choosing $T_{q}$ to be the result of many steps of an MCMC process.
Specifically, suppose $T_{q}^{k}$ represents the transition probability
resulting from applying k steps of an MCMC process, then Cover and Thomas
(2006) show that $\rm{KL}(T_{q}^{k}||q)\leq\rm{KL}(T_{q}^{k-1}||q),$ with
convergence assured under regularity conditions. This assures that provided
the chain length of CD is sufficiently large, the method becomes equivalent to
maximum likelihood. Unfortunately, even though the convergence of MCMC is
exponential, it may take an extraordinarily long time to reach the equilibrium
distribution. In many common implementations a chain length of 1 is used,
resulting in a distribution that is starkly different from $q$.

If $q$ belongs to an exponential family, then the gradient of the objective
function is

$\frac{\delta}{\delta\eta_{j}}(\rm{KL}(p||q)-\rm{KL}(T_{q}||q))=-E_{p}(g_{j}(Y))+E_{T_{q}}(g_{j}(Y))-\frac{\delta
T_{q}}{\delta\eta_{j}}\frac{\delta\rm{KL}(T_{q}||q)}{\delta T_{q}}.$ (2)

The first term on the RHS of equation 2 can be approximated using the observed
values from the data ($g_{j}(y)$), and the second expectation can be
approximated by sampling from $T_{q}$ which is relatively inexpensive
computationally. The third term is problematic in that it can not be estimated
without evaluating $q$. Hinton (2002) suggested dropping it based on
simulations that suggested that it “is small and seldom opposes the resultant
of the other two terms.” This then results in the contrastive divergence
moment conditions

$g(y)=E_{T_{q}}(g(Y)).$ (3)

There is, however, no theoretical reason to believe that the third term is
ignorable in general. In the context of Restricted-Boltzmann machines, the
bias of the gradient approximation in equation 3 relative to the maximum
likelihood gradient goes to 0 as the number of MCMC steps in $T_{q}$ goes to
infinity Bengio and Delalleau (2009), Fischer and Igel (2011). However, this
result relies upon the chain $T_{q}$ approaching $q$ as the number of steps
$k$ increases, which is unrealistic as the typical number of steps is very
small.

##  2 Reframing the Problem

There are two main theoretical problems with the current development of CD.
Firstly, it is questionable whether the objective function in equation 1 finds
a $q$ close to $p$ when it is minimized,. Secondly, it is unknown whether the
minimization is adequately achieved when ignoring the third term of equation
2. In order to address these issues, we develop an alternative objective
function which provides a more principled foundation for the CD gradient
approximation. First we will introduce some notation. Let $y$ be indexed as
$y_{i}$ for $i\in(1,...,m)$, $b$ be a vector of sets of indices of $y$, with
union $a=\bigcup_{i}b_{i}$. Let $r(b)$ represent the probability of selecting
indices $b$. We can then define the joint probability of an observed Gibbs
chain as
$T_{q}(Y^{(1:k)},B|y^{(0)})=\prod_{i=1}^{k}q(Y^{(i)}_{B_{i}}|Y_{\backslash
B_{i}}^{(i-1)})r_{i}(B_{i}),$ where $B_{i}$ are the indices of $Y$ which have
non zero probability of changing at step $i$.

The conditional probability of the last step of the MCMC process, given $A$ is

$q^{*}(Y^{(k)}|y^{(0)},a)=\frac{1}{\pi(a)}\sum_{a=\bigcup_{i}b_{i}}\sum_{Y^{(1:k-1)}}T_{q}(Y^{(1:k)},b|y^{(0)})$

where $\pi(A)=\sum_{A=\bigcup_{i}b_{i}}\prod_{i}^{k}r(b_{i})$. Because the
Gibbs update is a valid transition kernel for the unconditional distribution,
$q^{*}$ is a valid kernel for the conditional distribution given $a$ in that
it satisfies detailed balance.

$\displaystyle\sum_{y^{(0)}_{a}}q^{*}(Y^{(k)}|y^{(0)},a)q(y^{(0)}_{a}|y^{(0)}_{\backslash
a})$ $\displaystyle=$
$\displaystyle\frac{1}{\pi(a)}\sum_{a=\bigcup_{i}b_{i}}(\prod_{i}^{k}r(b_{i}))\sum_{y^{(0)}_{a}}q(y^{(0)}_{a}|y^{(0)}_{\backslash
a})\sum_{Y^{(1:k-1)}}\prod_{i=1}^{k}q(Y^{(i)}_{B_{i}}|Y_{\backslash
B_{i}}^{(i-1)})$ $\displaystyle=$
$\displaystyle\frac{1}{\pi(a)}\sum_{a=\bigcup_{i}b_{i}}(\prod_{i}^{k}r(b_{i}))q(Y^{(k)}_{a}|Y^{(k)}_{\backslash
a}=y^{(0)}_{\backslash a})$ $\displaystyle=$ $\displaystyle
q(Y^{(k)}_{a}|Y^{(k)}_{\backslash a}=y^{(0)}_{\backslash a})$

Thus, while the MCMC distribution $T_{q}$ is unlikely to be close to the
equilibrium distribution $q$ after a small number of steps, if the MCMC
sampler is chosen appropriately $q^{*}$ may be very close to
$q(Y_{a}|y_{\backslash a})$ because of the drastically reduced dimensionally
(i.e. $|a|<<m$).

###  2.1 Augmented Divergence

Now let us consider an augmented objective function for a single subset $a$.
We define this to be

$\displaystyle d_{a}(p,q)$ $\displaystyle=$
$\displaystyle\rm{KL}(p||q)-\rm{KL}(p_{m}||q_{m})+\sum_{y}\log(\frac{q(Y_{a}|y_{\backslash
a})}{q^{*}(Y_{a}|y_{\backslash a},y^{(0)})})p(y)$ $\displaystyle=$
$\displaystyle\rm{KL}(p||q^{*}p_{m})\geq 0$

where $p_{m}(y_{\backslash a})=\sum_{y_{a}}p(y)$ and $q_{m}(y_{\backslash
a})=\sum_{y_{a}}q(y)$ are the marginal distributions of $p$ and $q$.
Similarly, let us define the conditional distributions as
$p_{c}=p(Y_{a}|y_{\backslash a})$ and $q_{c}=q(Y_{a}|y_{\backslash a})$.

To see why $d_{a}$ is a good augmentation let us consider the first two terms.
From Huber (1985) and Lyu (2011) we know that
$\rm{KL}(p||q)-\rm{KL}(p_{m}||q_{m})\geq 0$, with equality when $p=q$. The
third term represents the discrepancy between the conditional distribution of
$q$ and the approximation $q^{*}$ averaged over the true distribution, and
when $p=q$ it reduces to $d_{A}(p,p)=\rm{KL}(p||p^{*}p_{m})\approx 0.$ Thus,
$d$ reduced to near its lower bound when $q$ fits the data generating
distribution perfectly. Further, if $p^{*}=p_{c}$ then the equality is
achieved. This equation can also be used to justify setting $y^{(0)}$ to be
from the sample data, as starting the MCMC chain from the equilibrium
distribution $p$ will cause (at least approximately) the resulting chain
$q^{*}$ to be closest to $q_{c}$ when $q=p$. This implies that even if the
lower bound is not achieved at $p=q$, the minimum will be nearly achieved, and
the divergence may be approximated with

$\displaystyle\operatorname*{\arg\\!\min}_{q}d_{a}(p,q)$ $\displaystyle=$
$\displaystyle\operatorname*{\arg\\!\min}_{q^{*}}\rm{KL}(p||q^{*}p_{m})$
$\displaystyle=$
$\displaystyle\operatorname*{\arg\\!\max}_{q^{*}}E_{p}(\log(q^{*}(Y_{a_{i}}|y_{\backslash
a},y^{0})))$ $\displaystyle\approx$
$\displaystyle\operatorname*{\arg\\!\max}_{q^{*}}\log(q^{*}(y_{a_{i}}|y^{0}))$
$\displaystyle=$
$\displaystyle\operatorname*{\arg\\!\max}_{q^{*}}\ell_{a}(q^{*}|y).$

So the augmented divergence is approximately minimized at by maximizing the
observed log likelihood of the transitional distribution $q^{*}$.

###  2.2 Combined Augmented Divergence

$d_{a}$ is a good measure of the local information about $p$ within $a$, but
we wish to have a more global criterion. Now we will combine all possible
subsets to generate a global augmented divergence by simply taking a weighted
average.

$cd(p,q)=\sum_{A}d_{A}(p,q)\pi(A)=E_{\pi}(\ell_{A}(q^{*}|y)).$ (6)

Since each one of the individual $d_{A}$ obtains a minimum at or near 0 when
$p\approx q$, the combined augmented divergence $cd$ is also a reasonable
measure of the divergence between $p$ and $q$ aggregating the local
information in each subset $A$. Small $cd$ values close to 0 indicate good
agreement between $p$ and $q$ over the subsets, and large values indicate poor
model fit. We may then minimize $cd$ subject to $q$ in order to find a good
model fit.

If $q$ is exponential family, and $q^{*}\approx q_{c}$ (which is not
unreasonable even for modest numbers of MCMC steps), then taking the gradient
and setting it to zero yields

$\displaystyle 0$ $\displaystyle=$ $\displaystyle E_{\pi}(\frac{\delta}{\delta
q}\log(q^{*}(y_{a}|y))$ $\displaystyle\approx$ $\displaystyle
E_{\pi}(g(y)-E_{q^{*}(Y_{A}|y_{\backslash A})}(g(Y)))$ $\displaystyle=$
$\displaystyle g(y)-E_{T}(g(Y)).$

Thus, provided that our MCMC approximation is reasonable within the restricted
sample space of $a$, the minimum of $cd$ occurs approximately when the
contrastive divergence moment conditions $g(y)=E_{T}(g(Y))$ are satisfied.

If $q$ is not exponential family or if $q^{*}$ is appreciably different from
$q$, then we may still have some hope that the moment conditions will provide
good inferences. The quantity in equation 6 can be recognized as the empirical
analog of a set of generalized method of moments conditions
$E_{T_{q}}(E_{\pi}(\frac{\delta}{\delta q}\log(q^{*}(Y_{A}|y^{(0)}))))=0$.
This implies that minimizing $cd$ is equivalent to fitting $T$ via moment
equations. In general it is not true that these moment equations reduce to the
CD moment conditions, but using the CD moment conditions will provide a fit of
$T$ which matches the data and thus will likely drive $cd$ to near its
minimum.

###  2.3 A Variational Approximation to $\ell_{A}$

Even if $q^{*}$ is not close to $q_{c}$ in distribution, the contrastive
divergence moment conditions may be derived as a variational approximation to
$cd$. Let $q$ belong to an exponential family. We may approximate a single
augmented divergence with:

$\displaystyle\log(q^{*}(y^{(k)}|a))$ $\displaystyle=$
$\displaystyle\sum_{b:a=\bigcup_{i}b_{i}}\sum_{y^{(1:(k-1))}}(\log(\frac{q^{*}(y^{(1:k)},b|a)}{f(y^{(1:(k-1))},b|a)})f(y^{(1:(k-1))},b|a)$
$\displaystyle+\log(\frac{f(y^{(1:(k-1))},b|a)}{q^{*}(y^{(1:(k-1))},b|y^{(k)},a)})f(y^{(1:(k-1))},b|a))$
$\displaystyle\gtrsim$
$\displaystyle\sum_{b:a=\bigcup_{i}b_{i}}\sum_{y^{(1:(k-1))}}\log(\frac{q^{*}(y^{(1:k)},b|a)}{f(y^{(1:(k-1))},b|a)})f(y^{(1:(k-1))},b|a)$
$\displaystyle=$ $\displaystyle
E_{f}(log(q_{c}(y^{(k)}_{B_{k}}|Y^{(k-1)}_{\backslash B_{k}})))$
$\displaystyle-\rm{KL}(f(y^{(1:k-1)},b|a)||q^{*}(y^{(1:k-1)},b|a))$
$\displaystyle=$ $\displaystyle U(y^{(k)},q,f)$

where $f$ is an arbitrary distribution, and with equality in 2.3 obtained when
$q^{*}(y^{(1:(k-1))},b|Y^{(k)},a)=f(y^{(1:(k-1))},b|a)$. The derivative of
equation 2.3 is then

$\frac{\delta
U(y^{(k)},q,f)}{\delta\eta}=g(y^{(k)})-E_{f}(E_{q_{c}}(g(Y^{(k)})|Y^{(k-1)}_{\backslash
B_{k}}))+\sum_{i=1}^{k-1}E_{f}(g(Y^{(i)}))-E_{f}(E_{q_{c}}(g(Y^{(i)})|Y^{(i-1)}_{\backslash
B_{i}}))$

Substituting $f=q^{*}(y^{(1:k-1)},b|a)$ into the derivative causes the terms
in the summation to cancel out yielding

$\frac{\delta
U(y^{(k)},q,f=q^{*})}{\delta\eta}=g(y^{(k)})-E_{q^{*}}(g(Y^{(k)})|a).$

Note that the inequality in 2.3 is tight as
$\rm{KL}(q^{*}(y^{(1:(k-1))},b|Y^{(k)},a)||q^{*}(y^{(1:(k-1))},b|a))$ is
small, especially when $k$ is large. We may combine the individual results by
minimizing equation 6, resulting in $E_{\pi}(\log(q^{*}(y^{(k)}|a)))\gtrsim
E_{\pi}(U(y^{(k)},q,f=q^{*})).$ Taking taking the derivative, and setting it
to 0 gives us the fixed point $g(y^{(k)})=E_{T_{q}}(g(Y^{(k)})),$ which we
recognize as the contrastive divergence moment conditions.

So we have three principled ways of justifying our use of the contrastive
divergence moment conditions within the combined augmented divergence
framework. First, if $q^{*}$ is close enough to $q_{c}$ to be approximately an
exponential family, then the moment conditions fall rather naturally out of
$cd$. Otherwise, by observing that minimizing $cd$ is equivalent to finding a
GMM fit of $T$, the CD moment conditions can be justified as an alternate GMM
fit. Further, the conditions may be derived from a variational approximation
to $cd$.

##  3 Special Cases

For most cases, the solution to equation 6 is not tractable analytically,
requiring us to use the CD moment conditions to find an approximate solution.
There are however some cases where we can solve the equation directly, and
these examples will provide us some intuition about when CD is likely to work
well, and when it has limited utility.

Composite Likelihood: Suppose that $T$ is a single step of a blocked Gibbsian
kernel such that $T(Y,A|y^{(0)})=q_{c}(Y_{A}|y^{(0)}_{\backslash A})r(A).$ We
can see that in this case, our $q^{*}$ simplifies to
$q^{*}(Y^{(k)}|y^{(0)},a)=q(Y_{A}|y^{(0)}_{\backslash A})$, and if we let
$y^{(0)}=y$, then the objective function becomes

$\operatorname*{\arg\\!\max}_{q}\sum_{A}\log(q_{c}(Y_{A}=y_{A}|y_{\backslash
A}))r(A),$

which we recognize as the objective function for composite likelihood. If $q$
is exponential family, the first derivative is
$\frac{\delta}{\delta\eta}=g(y)-E_{T}(g(Y))$ with a second derivative of
$\frac{\delta^{2}}{\delta\eta_{i}\delta\eta_{j}}=cov_{T}(g_{i}(Y),g_{j}(Y)).$
So in this case, the CD moment conditions are identical to the moment
conditions for composite likelihood.

If $r$ is chosen to be the uniform distribution over $A$ such that $A$ only
contains a single element, then the objective function becomes
$\operatorname*{\arg\\!\max}_{q}\sum_{i=1}^{m}\log(q(Y_{i}=y_{i}|y_{\backslash
i})),$ which is known as the pseudo-likelihood. Pseudo-Likelihood is a simple
form of composite likelihood inference and will serve as a baseline to
evaluate CD inference in section 5.

The more elements present in $A$, the closer the composite likelihood estimate
will be to the maximum likelihood estimate. On the other hand, the
computational complexity of sampling from $q_{c}$ increases exponentially as
more indices are added. Alternatively, but relatedly, suppose that $T$ selects
a subset $A$ with probability $\pi(A)$, and then runs an MCMC sampler within
$A$ on the conditional distribution $q_{c}$ until the chain reaches
equilibrium. Then the transition kernel can be written as

$T(Y^{(k)},A|y^{(0)})=\pi(A)\sum_{i=1}^{k-1}\prod_{i=1}^{k}t_{i}(Y^{(i)}_{B_{i}}|Y^{(i-1)},B_{i}.y^{(0)})r_{i}(B_{i}|Y^{(i-1)})=q_{c}(Y_{A}=y_{A}|y^{(0)}_{\backslash
A})\pi(A)$

which by identical argument results in a composite likelihood solution with
weights defined by $\pi$. So, if $T$ results in the exact equality of $q^{*}$
and $q_{c}$, either though direct sampling via a single blocked Gibbs step, or
though running a sampler to convergence within each $A$, then the CD objective
function is equivalent to the composite likelihood objective function.

Single Scan Gibbs: Suppose that $T$ is a Gibbsian kernel sequential updating
$y_{s}:y_{s+k}$ with $k\leq m$, $s$ chosen at random, and indices wrapping
when needed. Then the log likelihood for a given $A$ is

$\ell_{A}(q^{*}|y)=\log(q^{*}(y^{i}_{A}|y^{i}_{\backslash
A},y^{(0)}_{A}))=\sum_{i=1}^{k}\log(q(Y_{i}=y_{i}|y_{\backslash
A},y_{s:s+i-1},y^{0}_{(i+1):(s+k)}))$

and setting our initial MCMC state equal to our observation (i.e. $y^{(0)}=y$)
we arrive at
$\ell_{A}(q^{*}|y)=\sum_{i=1}^{k}\log(q(Y_{i}=y_{i}|y^{(0)}_{\backslash i}).$
Applying equation 6 results in
$\operatorname*{\arg\\!\max}_{q}\sum_{i=1}^{m}\log(q(Y_{i}=y_{i}|y_{\backslash
i})),$ which is identical to the pseudo-likelihood function, indicating that
our longer chain gave no additional improvement over a single step of random
scan Gibbs. Indeed, any k-step Gibb update which makes only one pass though
$A$ will result in a pseudo-likelihood objective function, so if we wish to
make non-pseudo-likelihood inferences, the kernel must revisit indices that it
has already sampled.

Conditional Independence: Suppose that $T$ is a single step blocked Gibbs
kernel. Let $r$ be such that it chooses $A$ to be the set of two indices at
random $i,j$ such that they are conditionally independent (i.e.
$q(Y_{i},Y_{j}|y_{\backslash\\{i,j\\}})=q(Y_{i}|y_{\backslash\\{i,j\\}})q(Y_{j}|y_{\backslash\\{i,j\\}})$).
With $y^{0}=y$ we obtain $\ell_{A}(q^{*}|y)=\log(q(Y_{i}=y_{i}|y_{\backslash
i})q(Y_{j}=y_{j}|y_{\backslash j})),$ and minimizing equation 6 results in the
objective function
$\operatorname*{\arg\\!\max}_{q}\sum_{i=1}^{m}\log(q(Y_{i}=y_{i}|y_{\backslash
i}))\sum_{j=1}^{m}r(\\{i,j\\}),$ which is again the objective function for
pseudo-likelihood (weighted by $\sum_{j=1}^{m}r(\\{i,j\\})$), so if non-
pseudo-likelihood inferences are the goal, then $A$ must be chosen such that
there is dependence among its members.

One and a Half Pass Gibbs: Suppose that $q$ represents an exponential family,
and that $m=2$. Let $T$ be a Gibbs sampler sampling $y_{1}$ then $y_{2}$, then
back to $y_{1}$. In this case, with $y^{0}=y$ we obtain

$\ell_{A}(q^{*}|y)=\log(q(Y_{1}=y_{1}|y_{\backslash
1})\sum_{y_{1}^{\prime}}q(Y_{1}=y_{1}^{\prime}|y_{\backslash
1})q(Y_{2}=y_{2}|y_{\backslash\\{1,2\\}},y_{1}=y_{1}^{\prime})).$

Taking the derivative yields

$\frac{\delta\ell_{A}(q^{*}|y)}{\delta\eta_{i}}=g_{i}(y)-E_{Y_{1}|y_{\backslash
1}}(g_{i}(Y))+2E_{Y^{\prime}|y_{\backslash
1}}(g(Y^{\prime})f-g(Y^{\prime}))-E_{Y^{\prime}|y_{\backslash
1}}(fE_{Y_{2}|Y_{\backslash
2}^{\prime}}(g(\\{Y_{1}^{\prime},Y_{2},y_{3:m}\\})))$

where
$f=\frac{q(Y_{2}=y_{2}|Y_{1}=y_{1}^{\prime})}{E(q(Y_{2}=y_{2}|Y_{1}=Y_{1}^{\prime}))}$,
and $Y^{\prime}=\\{Y_{1}^{\prime},y_{2}\\}$. Each of the expectations in the
derivative may be approximated by sampling from $T$, however it is apparent
that even with this simple Gibbsian updating, the likelihood becomes difficult
to express. More complex $q^{*}$ representing longer MCMC chains will only
increase this complexity. Thus it is useful (even for the above simple case)
to use the CD moment conditions to approximately optimize the objective.

##  4 Learning Algorithms

In the context of maximum likelihood inference in exponential families, one
can perform maximum likelihood inference using the stochastic approximation
algorithm (Robbins and Monro, 1951), where the update equation

$\eta^{i+1}=\eta^{i}+\gamma_{i}(g(y)-E_{q}(g(Y)|\eta^{i}))$

converges to the maximum likelihood solution when the gradient coefficient is
allowed to decrease geometrically (i.e. $\gamma_{i}=\frac{a}{i}$ for some
$a>0$). The challenge here from a computation perspective is that though
$E_{q}$ can be theoretically estimated by MCMC, the chain length required for
convergence is too long to make this practical. This is precisely the reason
that the CD moment conditions are attractive, as they do not require full
convergence of the MCMC chain, rather only a limited number of steps ($k$).

The CD update step takes the simplistic approach of replacing $q$ in the
stochastic approximation with $T$ yielding

$\eta^{i+1}=\eta^{i}+\gamma_{i}(g(y)-E_{T}(g(Y)|\eta^{i})).$

$T$ is very easy and quick to sample from, since each sample takes just a few
steps from an MCMC process to compute. The expectation can then be
approximated by simply taking the sample expectation
($\frac{1}{n}\sum_{i}^{n}g(y^{i})$ where $y^{i}$ are independent samples from
$T$). Yuille (2004) showed that this update equation converges under some
regularity conditions, and provided some conditions under which the converged
upon value is the maximum likelihood solution. Unfortunately, the conditions
required for convergence to the maximum likelihood are unlikely to be met
unless the MCMC chain is run long enough to reach equilibrium.

Stochastic approximation approaches are similar to gradient descent, in that
convergence is generally linear in the number of iterations. There may be some
advantage to be gained by using a Newton-like update. If $q^{*}$ is
sufficiently close to $q_{c}$ to be approximately exponential family, then the
second derivative within equation 6 is $H(\eta)=cov_{T}(g(Y),g(Y)),$ which is
identical to the composite likelihood hessian in section 3. We can then use a
Newton-like update

$\eta^{i+1}=\eta^{i}-[H(\eta_{i})]^{-1}(g(y)-E_{T}(g(Y)|\eta^{i})).$

Both the hessian and expectation can be approximated with the same set of
samples from $T$, so this algorithm takes similar computational effort to the
stochastic gradient algorithm at each step. If the hessian is a good
approximation, this algorithm will see quadratic convergence to the CD
solution. The down side is that convergence is not guaranteed, especially when
the hessian approximation is poor.

##  5 Experiments

One popular class of exponential families are Exponential-family Random Graph
Models (ERGM) (Frank and Strauss, 1986). These models have been used widely to
investigate the structure of social interactions (Hunter et al., 2008), and
are typically fit either by finding the maximum pseudo-likelihood estimate
(MPLE) (van Duijn et al., 2009) or using MCMC methods to find a maximum
likelihood estimate (MCMC-MLE) (Hunter and Handcock, 2006). As such, they
provide a useful case study for approximate inferential methods which have
less computational complexity than MCMC-MLE and more accurate inference than
MPLE. ERGMs model the presence or absence of a relationship between a set of
nodes which typically represent individual people. Let the dyad $y_{i,j}$ be 1
if there is a connection between node $i$ and $j$, and 0 if there is none.
Connections are considered to be undirected (i.e. $y_{i,j}=y_{j,i}$).

The dataset we use comes from the National Longitudinal Study of Adolescent
Health (Add Health), which in addition to collecting data on health related
behaviors, also recorded information about the social networks of subjects
(Harris et al., 2003). One particular high-school surveyed by the study had
1270 students, with 361, 309, 346, and 254 students in 9th, 10th, 11th and
12th grade respectively. Students were asked to select up to five close male
and five close female students as friends. A relationship is considered to be
present if and only if both students nominate each other. We will consider a
simple four parameter ERGM. The first term in the model is the number of edges
in the network (${\rm edges}=\sum_{i>j}y_{i,j}$). The second term is the
number of students with no connections (${\rm
isolates}=\sum_{i}I(\sum_{j}y_{i,j}>0)$, where $I$ is the indicator function).
We expect that more connections are present between members of the same grade,
so the third term will model this with the count of connections between
students of the same grade (${\rm
nodeMatch}=\sum_{i>j}I(grade_{i}=grade_{j})y_{i,j}$). Another important social
relationship is transitivity, which is the tendency of the friend of my friend
to also be my friend. Simplistic statistics used to model this have the
tendency to display poor statistical properties such as phase
transitionsHandcock (2003). For this reason, a more robust measure known as
the geometrically weighted edgewise shared partner (GWESP) statistic is used.
See Hunter and Handcock (2006) for the mathematical specification of this
statistic as well as a more thorough justification. GWESP statistics also
require the specification of a curve parameter $\alpha$, which we set to value
of $\frac{2}{3}$ based on the MLE goodness of fit. Finally, since students can
only select up to a maximum of 10 friends, we should incorporate this
constraint into the model. This is done by including an offset which is set to
$-\infty$ if any node has more than 10 connections, and 0 otherwise.

The special cases outlined in section 3 lead us to believe that the choice of
MCMC algorithm will be important to improving our estimate over the MPLE. The
dimensionality of $y$ is 805,815, so a sequential scan Gibbian update would
require a huge number of steps before it revisits any elements of $y$.
Similarly we might expect a random scan Gibbsian update to perform similarly
due to the low probability of revisiting dyads. For this reason, we consider
three alternate MCMC schemes based on the intuition that dyads connected to
the same node are unlikely to be conditionally independent. The Node-$s$ MCMC
kernel proceeds by selecting (for each chain) a node at random, and then
selecting $s$ dyads incident to that node. A simple random scan Gibbs is then
performed for $k$ steps within the selected dayds.

Figure 1 shows the mean value parameters for CD inference using Node-$s$
kernels with $s=$ 200, 500 and 1269 (Full) at different chain lengths up to
$2^{14}$. Also shown are the results for random scan Gibbs, and the MLE and
MPLE estimates. Each CD fit was done using the Newton-like updates in section
4, leading to much faster convergence than the gradient update. Mean values
$\mu$ are shown rather than natural parameters $\eta$ because $\mu$ is more
useful in determining how well a parameter set matches the MLE. The pseudo-
likelihood estimate fits well for the isolates, edges and nodeMatch terms as
the mean values are all within one standard deviation of the MLE’s mean value
among networks simulated from the MLE model. The MPLE grossly underestimates
the amount of transitivity in the network as represented by the GWESP term,
with a mean value more than 3 standard deviations lower than the MLE. The
distribution of the GWESP term under the MLE is displayed in the density band
along the left side of the plot, which shows that the average GWESP term for
networks generated under the MPLE model are virtually unheard of in networks
generated by the MLE. This indicates that there is significant room for
improvement on this fit.

Figure 1:  Mean value parameter estimates for different MCMC methods and step
lengths. The distributions on the left edges of the plots represent the
distributions of the statistics $g(Y)$ at the MLE.

The Random Scan Gibbs solutions are indistinguishable from pseudo-likelihood,
even with a $k$ of 16,384. This confirms our theoretical result that the
asymptotics of $T$ approaching $q$ as $k$ increases provides a poor foundation
for CD-k inference, as a hugely massive $k$ would be required in order to
reach equilibrium in the high dimensional space of $y$. However, by
restricting our Gibbs sampling to be within subsets of $y$, significant
improvement can be made over pseudo-likelihood. For the GWESP term, we see
that as $s$ and $k$ increase, the solutions become closer to the MLE, though
there still exists some downward bias. We also see some improvement in the
edges and nodeMatch terms relative the the MPLE, and all estimates are close
to the MLE for the isolates parameter. Node-Full reaches near its best
estimate when $k\geq 2048$, which is less than 2 times larger than the
dimensionality of the subset (1269); Far too few to reach full equilibrium,
but enough so that elements of $y$ have a high probability of being revisited.

##  6 Conclusion: Why, When and How

This paper explored why contrastive divergence works using a novel combined
augmented divergence framework. We illustrated how this new $cd$ objective
function alleviates some of the theoretical concerns with the original
objective function and provided three justifications of the CD moment
conditions. From the examples in section 3 we may surmise that MCMC kernels
work for CD when they visit and revisit highly related sections of $y$ in a
modest number of steps, otherwise CD inference will have no advantage over the
Pseudo-Likelihood solution. Finally, we showed how to implement CD using a
Newton-like update equation. The simulation study conducted used the Newton-
like update with success, and validated our heuristics about what types of
kernels are best.

## References

  * Hinton [2002] Geoffrey E. Hinton.  Training products of experts by minimizing contrastive divergence.  _Neural Computation_ , pages 1771–1800, 2002.
  * Fischer and Igel [2014] Asja Fischer and Christian Igel.  Training restricted boltzmann machines: An introduction.  _Pattern Recognition_ , pages 25–39, 2014.
  * Barndorff-Nielsen [1978] Ole E. Barndorff-Nielsen.  _Information and Exponential Families in Statistical Theory_.  Wiley, New York, 1978.
  * Lyu [2011] Siwei Lyu.  Unifying non-maximum likelihood learning objectives with minimum kl contraction.  In _NIPS_ , pages 64–72, 2011.
  * Cover and Thomas [2006] Thomas Cover and Joy Thomas.  _Elements of Information Theory_.  Wiley-Interscience, 2006.
  * Bengio and Delalleau [2009] Y. Bengio and O. Delalleau.  Justifying and generalizing contrastive divergence.  _Neural Computation_ , 21:1601–1621, 2009.
  * Fischer and Igel [2011] Asja Fischer and Christian Igel.  Bounding the bias of contrastive divergence learning.  _Neural Computation_ , 23(3):664–673, 2011.
  * Huber [1985] Peter J Huber.  Projection pursuit.  _The annals of Statistics_ , pages 435–475, 1985.
  * Robbins and Monro [1951] H. Robbins and S. Monro.  A stochastic approximation method.  _The Annals of Mathematical Statistics_ , 22(3):400, 1951.
  * Yuille [2004] Alan L Yuille.  The convergence of contrastive divergences.  In _Advances in Neural Information Processing Systems_ , pages 1593–1600, 2004.
  * Frank and Strauss [1986] Ove Frank and David Strauss.  Markov graphs.  _Journal of the American Statistical Association_ , 81(395):832–842, 1986.
  * Hunter et al. [2008] David R Hunter, Mark S Handcock, Carter T Butts, Steven M Goodreau, and Martina Morris.  ergm: A package to fit, simulate and diagnose exponential-family models for networks.  _Journal of Statistical Software_ , 24(3):nihpa54860, 2008.
  * van Duijn et al. [2009] MA van Duijn, KJ Gile, MS Handcock, et al.  A framework for the comparison of maximum pseudo likelihood and maximum likelihood estimation of exponential family random graph models.  _Social networks_ , 31(1):52–62, 2009.
  * Hunter and Handcock [2006] David R. Hunter and Mark S. Handcock.  Inference in curved exponential family models for networks.  _Journal of Computational and Graphical Statistics_ , 2006.
  * Harris et al. [2003] Kathleen M. Harris, Francesca Florey, Peter S. Bearman, Jo Jones, and J. Richard Udry.  The national longitudinal of adolescent health: Research design, 2003\.  URL http://www.cpc.unc.edu/projects/addhealth/design.
  * Handcock [2003] Mark S. Handcock.  Assessing degeneracy in statistical models of social networks.  Working paper #39, Center for Statistics and the Social Sciences, University of Washington, 2003.  URL {http://www.csss.washington.edu/Papers}.