# Stochastic kinetic theory applied to coarse-grained polymer model

Shangren Zhu  Patrick T. Underhill<EMAIL_ADDRESS>Rensselaer Polytechnic
Institute, Troy NY

(8 November, 2023)

###### Abstract

A stochastic field theory approach is applied to a coarse-grained polymer
model that will enable studies of polymer behavior under non-equilibrium
conditions. This article is focused on the validation of the new model in
comparison to explicit Langevin equation simulations under conditions with
analytical solutions. The polymers are modeled as Hookean dumbbells in one
dimension, without including hydrodynamic interactions and polymer-polymer
interactions. Stochastic moment equations are derived from the full field
theory. The accuracy of the field theory and moment equations are quantified
using autocorrelation functions. The full field theory is only accurate for
large number of polymers due to keeping track of rare occurrences of polymers
with a large stretch. The moment equations do not have this error because they
do not explicitly track these configurations. The accuracy of both methods
depends on the spatial degree of discretization. The timescale of
decorrelation over length scales bigger than the spatial discretization is
accurate, while there is an error over the scale of single mesh points.

##  I Introduction

Many interesting phenomena have been observed in polymeric fluid flows,
especially in concentrated systems. Flow can introduce inhomogeneities into
polymer concentrations which causes phase separation [1, 2, 3]. This
concentration inhomogeneity phenomenon is observed in different polymer
systems, including aqueous polymer solutions[4], entangled polymer systems[5],
and charged systems[6]. The membrane industry has also taken advantage of
phase separation in membrane production[7, 8].

To explain the behavior of polymer systems, several theories have been
developed. The Langevin equation model is one of the basic models for particle
or polymer systems involving Brownian dynamics. It calculates the spatial
position of each particle or polymer. This makes it effective in dilute or
semi-dilute system modeling[9, 10]. However, the Langevin equation model
becomes less effective and has high computational cost when applied to
concentrated systems. Over the decades, field theories have been utilized to
create different models in concentrated polymer systems. One example is the
equilibrium field theory[11, 12]. The equilibrium field theory represents the
polymer system as polymer chains interacting with a conjugate field by
Hubbard-Stratonovich transformations. It successfully decouples the polymers
and avoids the massive degree of freedom in typical particle-based theory.

A common implementation of equilibrium field theories entails a mean-field
assumption. This may be a good approximation for a system far from a critical
point (such as system temperature is smaller than critical temperature $T\ll
T_{c}$). However, as the system nears its critical point, the fluctuations
become important, and the mean field value is no longer accurate.[13, 14, 15]
The role of fluctuations can be included by sampling the field theory via
complex Langevin dynamics or via an approximation such as random phase
approximation (RPA). These approaches are limited to equilibrium fluctuations.

One approach to going beyond equilibrium is the single chain in mean-field
(SCMF) method [16]. This method tracks the polymer position with the Langevin
equation and integrates over time while interacting with conjugate fields.
Another approach for nonequilibrium polymer systems, especially for
interactions and phase separations, is the two-fluid model[17, 18]. This model
involves the dynamics of a field for each component, which evolves due to
conservation laws. The two-fluid model requires constitutive equations for
these components, which may not be obtainable in some circumstances or may be
the desired output of a model rather than as an input.

This article uses a field theory that applies to non-equilibrium systems
naturally. It takes the form of a stochastic kinetic theory, which was
previously derived by Dean[19] to describe the density of a stochastic
particle system. This stochastic field approach has been applied to multiple
areas, including interactive particles, bacteria, and active matter [20, 21,
22]. We extend this approach to the dynamics of polymers. The dynamics of the
density is derived from the corresponding Langevin equation for polymers that
contain internal degrees of freedom.

In this article, Section II presents the general formalism of the polymer
bead-spring chain model. The Langevin equations for bead positions describe
the dynamics of the model. This model is used to derive the dynamics of the
density distribution. This article is focused on validating the new stochastic
fluctuation theory in polymer systems. In Section III, the equations are
simplified for a Hookean dumbbell sufficient to examine the characteristics of
polymer dynamics and elasticity. To reduce the dimensionality of the variable
space, moment equations are derived. In Section IV, the accuracy of the
numerical implementation of the stochastic field theory is quantified in one
dimension.

##  II General Formalism: Bead-spring Chains

In this article, polymers are represented by a bead-spring chain model. It is
a coarse-grained model that tracks polymer spatial positions at several beads
along the polymer backbone. The connecting springs represent the entropic
elasticity. Viscous drag forces and Brownian forces on the beads result from
an implicit solvent. Consider $N_{p}$ number of polymers with $N_{b}$ beads in
each polymer chain. In the overdamped limit, the dynamics of the beads are
represented by Langevin equations. If ignoring hydrodynamic interactions
between the beads, the Langevin equation for the position of bead $i$ within
polymer chain $j$ is

$d\bm{r}_{i,j}=\left(\frac{\bm{F}_{int,i,j}+\bm{F}_{ext,i,j}}{\zeta}\right)dt+\sqrt{\frac{2k_{B}T}{\zeta}}d\bm{W}_{i,j}(t)$
(1)

where $\bm{r}_{i,j}$ is the position of bead $i$ within polymer chain $j$,
$\zeta$ is the viscous drag coefficient on a bead, $k_{B}$ is Boltzmann’s
constant, and $T$ is the absolute temperature. We have separated the non-
Brownian forces on the bead into the internal contribution $\bm{F}_{int}$ that
can depend only on the positions of beads within the same polymer chain and
the external contribution $\bm{F}_{ext}$ that depends on the positions of
beads on other polymer chains. The internal contribution includes spring
forces between beads, and drag forces from an external fluid flow. The
external contribution includes all effective interactions between polymer
chains. The $d\bm{W}_{i,j}$ are independent Brownian noise which obey

$\langle d\bm{W}_{i,j}\rangle=0$ (2) $\langle
d\bm{W}_{i,j}(t)d\bm{W}_{m,n}(t^{\prime})\rangle=\delta_{im}\delta_{jn}\bm{I}\delta(t-t^{\prime})dtdt^{\prime}$
(3)

where $\bm{I}$ is the identity tensor, $\delta_{im}$ and $\delta_{jn}$ are the
Kronecker delta, and $\delta(t-t^{\prime})$ is the Dirac delta function.

The density field of polymers $\psi$ is defined as

$\psi(\bm{r}_{1},\cdots,\bm{r}_{N_{b}},t)=\sum_{j=1}^{N_{p}}\delta\left(\bm{r}_{1}-\bm{r}_{1,j}(t)\right)\cdots\delta\left(\bm{r}_{N_{b}}-\bm{r}_{N_{b},j}(t)\right)$
(4)

Since the position of each bead in every polymer follows a Langevin equation,
the density field of polymers is also stochastic. Its dynamical equation can
be derived from the dynamics of the beads using Ito’s Lemma. If the number of
beads per polymer is large, the dimensionality of the function $\psi$ will be
large, which can present challenges to numerical methods.

A common minimal model of polymer elasticity is a dumbbell model ($N_{b}=2$),
which is the focus of the rest of the article. It is useful to change
variables from individual beads positions $\bm{r}_{1}$ and $\bm{r}_{2}$ to the
position of the polymer center mass $\bm{r}_{c}=(\bm{r}_{1}+\bm{r}_{2})/2$ and
spring vector $\bm{Q}=\bm{r}_{2}-\bm{r}_{1}$. In these variables, the field
equation is

$\begin{split}d{\psi}=&\left[-\frac{1}{\zeta}\nabla_{r}\cdot(\bm{F}_{int,c}+\bm{F}_{ext,c})\psi-\frac{1}{\zeta}\nabla_{Q}\cdot(\bm{F}_{int,Q}+\bm{F}_{ext,Q})\psi\right.\\\
&+\left.\frac{k_{B}T}{\zeta}\left(\frac{1}{2}\nabla_{r}^{2}{\psi}+2\nabla_{Q}^{2}{\psi}\right)\right]d{t}+\sqrt{\frac{k_{B}T}{\zeta}}\left[\nabla_{r}\cdot\left(\sqrt{{\psi}}d\bm{U}_{c}\right)+\nabla_{Q}\cdot\left(\sqrt{4{\psi}}d\bm{U}_{Q}\right)\right]\end{split}$
(5)

The variables $\bm{F}_{int,c}$ and $\bm{F}_{ext,c}$ are the internal and
external force contribution to the center mass position given by
$\bm{F}_{int,c}=\frac{\bm{F}_{int,1}+\bm{F}_{int,2}}{2}$ and
$\bm{F}_{ext,c}=\frac{\bm{F}_{ext,1}+\bm{F}_{ext,2}}{2}$. The variables
$\bm{F}_{int,Q}$ and $\bm{F}_{ext,Q}$ are the internal and external force
contribution to the spring vector $\bm{Q}$ given by
$\bm{F}_{int,Q}=\bm{F}_{int,2}-\bm{F}_{int,1}$ and
$\bm{F}_{ext,Q}=\bm{F}_{ext,2}-\bm{F}_{ext,1}$. The internal contributions
depend directly on the variables $\bm{r}_{c}$ and $\bm{Q}$. The external
contributions include interactions through a dependence on $\psi$.

The stochastic variables $d\bm{U}_{c}$ and $d\bm{U}_{Q}$ are independent
Gaussian noise terms with respect to $\bm{r}_{c}$ and $\bm{Q}$ which both obey

$\langle d\bm{U}_{i}\rangle=0$ (6) $\langle
d\bm{U}_{i}(\bm{r}_{c},\bm{Q},t)d\bm{U}_{i}(\bm{r}_{c}^{\prime},\bm{Q}^{\prime},t^{\prime})\rangle=\delta(\bm{r}_{c}-\bm{r}_{c}^{\prime})\delta(\bm{Q}-\bm{Q}^{\prime})\delta(t-t^{\prime})\bm{I}dtdt^{\prime}$
(7)

##  III Independent Hookean Dumbbells at Equilibrium

The equations simplify and have analytical solutions for independent Hookean
dumbbells at equilibrium. This enables a comparison between numerical
solutions and analytical solutions. The $N_{p}$ Hookean dumbbells have spring
constant $H$ and are placed in a system with dimension $d$. The domain for the
center of mass has size $L$ in each dimension and has periodic boundary
conditions. Nondimensionalized variables will be denoted using tildes. Lengths
are nondimensionalized by the characteristic length
$l_{c}=\sqrt{\frac{k_{B}T}{H}}$, so the typical stretch of the springs will be
order one. Times are nondimensionalized by the characteristic time
$t_{c}=\frac{2\zeta}{H}$, which is proportional to the time for the center of
mass to diffuse the stretch of the polymer.

The nondimensional Langevin equations for the position and stretch of polymer
$j$ are

$d\tilde{\bm{r}}_{c,j}=\sqrt{2}d\tilde{\bm{W}}_{c,j}(\tilde{t})$ (8)
$d\tilde{\bm{Q}}_{j}=-4\tilde{\bm{Q}}_{j}d\tilde{t}+\sqrt{8}d\tilde{\bm{W}}_{Q,j}(\tilde{t})$
(9)

The equation for the density $\psi$ can also be nondimensionalized. We define
$\tilde{\psi}=\frac{L^{d}}{N_{p}}\psi$ so that the dimensionless concentration
has a mean of one. The resulting dynamical equation is

$d\tilde{\psi}=\left[\nabla_{Q}\cdot\left(4\tilde{\bm{Q}}\tilde{\psi}\right)+\nabla_{r}^{2}\tilde{\psi}+4\nabla_{Q}^{2}\tilde{\psi}\right]d\tilde{t}+\sqrt{\frac{\tilde{L}^{d}}{N_{p}}}\nabla_{r}\cdot\left(\sqrt{2\tilde{\psi}}d\tilde{\bm{U}}_{c}\right)+\sqrt{\frac{\tilde{L}^{d}}{N_{p}}}\nabla_{Q}\cdot\left(\sqrt{8\tilde{\psi}}d\tilde{\bm{U}}_{Q}\right)$
(10)

For independent polymers, there is a contribution from $\bm{F}_{int,Q}$ in
which the spring force pushes the density towards a smaller stretch.

The function $\tilde{\psi}$ depends on both $\bm{r}_{c}$ and $\bm{Q}$. Many
observables like concentration or polymer alignment are related to moments of
$\tilde{\psi}$. These moments are defined as

$\tilde{\mu}_{0}(\bm{\tilde{r}}_{c},\tilde{t})=\int\tilde{\psi}(\bm{\tilde{r}_{c}},\bm{\tilde{Q}},\tilde{t})d\bm{\tilde{Q}}$
(11)
$\tilde{\bm{\mu}}_{1}(\bm{\tilde{r}}_{c},\tilde{t})=\int{\bm{\tilde{Q}}}\tilde{\psi}(\bm{\tilde{r}}_{c},\bm{\tilde{Q}},\tilde{t})d\bm{\tilde{Q}}$
(12)
$\tilde{\bm{\mu}}_{2}(\bm{\tilde{r}}_{c},\tilde{t})=\int{\bm{\tilde{Q}}\bm{\tilde{Q}}}\tilde{\psi}(\bm{\tilde{r}}_{c},\bm{\tilde{Q}},\tilde{t})d\bm{\tilde{Q}}$
(13)

where $\tilde{\mu}_{0}$, $\tilde{\bm{\mu}}_{1}$, and $\tilde{\bm{\mu}}_{2}$
correspond to scalar zeroth, vector first, and tensor second moment
respectively. These moments only depend on the position of the center of mass
and time. The equations for the dynamics of these moments can be derived by
performing the required integrals over Eqn 10.

The resulting equations are

$d\tilde{\mu_{0}}=(\nabla^{2}\tilde{\mu_{0}})d\tilde{t}+\sqrt{\frac{2\tilde{L}^{d}}{N_{p}}}\nabla\cdot\bm{g}_{c}^{(1)}$
(14)
$d\tilde{\bm{\mu}}_{1}=(\nabla^{2}\tilde{\bm{\mu}}_{1}-4\tilde{\bm{\mu}}_{1})d\tilde{t}+\sqrt{\frac{2\tilde{L}^{d}}{N_{p}}}\nabla\cdot\bm{g}_{c}^{(2)}-\sqrt{\frac{8\tilde{L}^{d}}{N_{p}}}\
\bm{g}_{q}^{(1)}$ (15)
$d\tilde{\bm{\mu}}_{2}=(\nabla^{2}\tilde{\bm{\mu}}_{2}-8\tilde{\bm{\mu}}_{2}+8\bm{I}\tilde{\mu_{0}})d\tilde{t}+\sqrt{\frac{2\tilde{L}^{d}}{N_{p}}}\nabla\cdot\bm{g}_{c}^{(3)}-\sqrt{\frac{8\tilde{L}^{d}}{N_{p}}}(\bm{g}_{q}^{(2)}+(\bm{g}_{q}^{(2)})^{T})$
(16)

where the $\bm{g}$ terms contain the noise. Both $\bm{g}_{c}$ and $\bm{g}_{q}$
are Gaussian noise terms defined by

$\bm{g}_{i}^{(n)}=\int\sqrt{\tilde{\psi}}d\tilde{\bm{U}}_{i}\tilde{\bm{Q}}^{n-1}d\tilde{\bm{Q}}$
(17)

The subscript denotes whether the noise came from a stochastic flux in the
center of mass position or spring vector and the superscript $(n)$ relates to
the number of factors of $\tilde{\bm{Q}}$ that are multiplied via the dyadic
product within the integral.

The $\bm{g}$ noise terms with respect to $\bm{r}_{c}$ and $\bm{Q}$ (i.e.
different subscripts) are independent of each other. However, the noise terms
with the same subscript but different superscripts are correlated, which must
be included when computing the dynamics of the moments. This is done by
relating the covariance of the noise to the moments, which requires a moment
closure approximation when computing the nonlinear dynamics. This will be
shown in detail later for the one-dimensional system.

##  IV Independent Hookean Dumbbells in One Dimension

###  IV.1 Three simulation methods

The validation of the stochastic field theory will be shown in a one-
dimensional system. The Langevin equation for Hookean dumbbells in a one-
dimensional system becomes

$d\tilde{r}_{c,j}=\sqrt{2}d\tilde{W}_{c,j}(\tilde{t})$ (18)
$d\tilde{Q}_{j}=-4\tilde{Q}_{j}d\tilde{t}+\sqrt{8}d\tilde{W}_{Q,j}(\tilde{t})$
(19)

The corresponding density function $\tilde{\psi}$ solves

$d\tilde{\psi}=\left[4\frac{\partial}{\partial\tilde{Q}}(\tilde{Q}\tilde{\psi})+\frac{\partial^{2}\tilde{\psi}}{\partial\tilde{r}_{c}^{2}}+4\frac{\partial^{2}\tilde{\psi}}{\partial\tilde{Q}^{2}}\right]d\tilde{t}+\sqrt{\frac{\tilde{L}}{N_{p}}}\frac{\partial}{\partial\tilde{r}_{c}}\sqrt{2\tilde{\psi}}d{\tilde{U}}_{c}+\sqrt{\frac{\tilde{L}}{N_{p}}}\frac{\partial}{\partial\tilde{Q}}\sqrt{8\tilde{\psi}}d\tilde{U}_{Q}$
(20)

The moment equations simplified to a one-dimensional system are

$d\tilde{\mu_{0}}=\frac{\partial^{2}\tilde{\mu}_{0}}{\partial\tilde{r}_{c}^{2}}d\tilde{t}+\sqrt{\frac{2\tilde{L}}{N_{p}}}\frac{\partial}{\partial\tilde{r}_{c}}g_{c}^{(1)}$
(21)
$d\tilde{\mu_{1}}=\left(\frac{\partial^{2}\tilde{\mu}_{1}}{\partial\tilde{r}_{c}^{2}}-4\tilde{\mu}_{1}\right)d\tilde{t}+\sqrt{\frac{2\tilde{L}}{N_{p}}}\frac{\partial}{\partial\tilde{r}_{c}}g_{c}^{(2)}-\sqrt{\frac{8\tilde{L}}{N_{p}}}g_{q}^{(1)}$
(22)
$d\tilde{\mu_{2}}=\left(\frac{\partial^{2}\tilde{\mu}_{2}}{\partial\tilde{r}_{c}^{2}}-8\tilde{\mu}_{2}+8\tilde{\mu}_{0}\right)d\tilde{t}+\sqrt{\frac{2\tilde{L}}{N_{p}}}\frac{\partial}{\partial\tilde{r}_{c}}g_{c}^{(3)}-2\sqrt{\frac{8\tilde{L}}{N_{p}}}g_{q}^{(2)}$
(23)

where all moments are scalars.

To apply these models in numerical simulation, we discretize them with finite
difference and finite volume methods. The discretized stochastic field theory
and moment equations are given in Appendices A and B respectively.

###  IV.2 Autocorrelation analysis

To quantitatively assess the accuracy of the stochastic fluctuations, we use
time autocorrelation to characterize the behavior of different simulation
trials. Before computing the time autocorrelation, moments are Fourier
transformed from $\tilde{r}_{c}$ to $\tilde{k}$ as

$\hat{\mu}_{n}(\tilde{k},\tilde{t})=\frac{1}{\tilde{L}}\int_{0}^{\tilde{L}}\tilde{\mu}_{n}(\tilde{r}_{c},\tilde{t})e^{-i\tilde{k}\tilde{r}_{c}}d\tilde{r}_{c}$
(24)

We define the autocorrelation as
$C_{n}=\langle\hat{\mu}_{n}(\tilde{k},\tilde{t})\hat{\mu}^{*}_{n}(\tilde{k},\tilde{t}^{\prime})\rangle-\langle|\hat{\mu}_{n}|^{2}\rangle$
where $*$ denotes the complex conjugate. The autocorrelation is a function of
the time difference $\Delta\tilde{t}=|\tilde{t}-\tilde{t}^{\prime}|$.

The autocorrelations can been obtained by solving the linearized equations
analytically to give

$C_{0}=\frac{1-\delta_{\tilde{k}0}}{N_{p}}\left[e^{-\tilde{k}^{2}\Delta\tilde{t}}\right]$
(25) $C_{1}=\frac{1}{N_{p}}\left[e^{-(\tilde{k}^{2}+4)\Delta\tilde{t}}\right]$
(26)
$C_{2}=\frac{1}{N_{p}}\left[(1-\delta_{\tilde{k}0})e^{-\tilde{k}^{2}\Delta\tilde{t}}+2e^{-(\tilde{k}^{2}+8)\Delta\tilde{t}}\right]$
(27)

The time autocorrelation of zeroth and first moments have a single exponential
since their linearized equations are independent of the other moments. The
autocorrelation of the second moment is a combination of two exponential
functions because the zeroth moment appears in Eqn 23. However, it reduces to
a single exponential for $\tilde{k}=0$ since the zeroth moment at
$\tilde{k}=0$ does not fluctuate.

###  IV.3 Dependence on number of polymers

For all three simulation methods, fits to the autocorrelation functions at
small wavevector are performed across a range of $N_{p}$. The system size is
$\tilde{L}=10$, and it is initially divided into $N_{bins}=50$. For the
Langevin equation and moment equations, the timestep is $d\tilde{t}=10^{-3}$.
For the stochastic field theory, the timestep is $d\tilde{t}=8\times 10^{-5}$.
We initially analyze the smallest wavevector that gives a nonzero correlation.
For the zeroth moment, that is $\tilde{k}=\frac{2\pi}{\tilde{L}}$. For the
first and second moments, that is $\tilde{k}=0$. To validate our models with
theoretical results, we fit the numerical solution with

$C_{i}=\frac{1}{N_{p}}a_{i}e^{-\frac{\Delta\tilde{t}}{\tau_{i}}}$ (28)

Fig 1 illustrates the fit parameters for the zeroth moment autocorrelation of
all three models with various polymer numbers $N_{p}$. According to Eqn 25,
the theoretical coefficients are $a_{0}=1$ and $\tau_{0}=1/\tilde{k}^{2}$.
From the simulations, the moment equations model matches with the theory at
all applicable $N_{p}$ within statistical error. The Langevin equation model
also agrees with the theoretical value within statistical error, but the
simulations were only performed up to $N_{p}=10^{4}$ due to the computational
cost. The stochastic field theory deviates from the analytical solution at
small system size $(N_{p}<10^{6})$.

The deviation is due to the larger relative fluctuations for small $N_{p}$
versus large $N_{p}$. The full SFT must keep track of the density of polymers
for each stretch, including the small number of polymers with large stretch.
The explicit integration algorithm given in Appendix A can lead to a negative
density at the end of a timestep. A simple algorithm is used to correct for
this and keep the density positive and the polymer number conserved, but it
leads to some error. For large $N_{p}$, no corrections are needed and the
results match with the theoretical values. For small $N_{p}$, corrections
occur after nearly every timestep, which alters the fluctuations. In
principle, the random fluxes could be sampled from a distribution at higher
computational cost such that negative densities cannot occur. In practice, the
moment equations will be used more often because of the smaller dimensionality
and larger allowable timestep.

Fig 2 and Fig 3 show the fitted parameters to the autocorrelations of the
first and the second moment. The results are similar to the results for the
zeroth moment. The simulations using the Langevin equations for each polymer
and the simulations using the moments match with the theory to within
statistical error. The results for the full SFT deviate from the theory for
small $N_{p}$.

To further investigate the impact of the corrections when the density becomes
negative, we simulate the full SFT with the center of mass domain divided into
different number of bins $N_{bins}$. When the total number of polymers are
spread over a smaller number of bins, the typical number of polymers per bin
increases. This leads to smaller relative fluctations and a smaller chance the
density becomes negative during time-integration. Fig 4 shows the fit
parameters to the zeroth moment autocorrelation with different $N_{bins}$.
With the smaller $N_{bins}$, fewer corrections are made and the errors are
smaller.

###  IV.4 Dependence on wavevector

We compare the performance under different wavevectors $\tilde{k}$ with
smaller timesteps to lower the impact of timestep error at larger wavevector.
The Langevin equation model and moment equation are both performed at
$N_{p}=10^{3}$. Both methods are simulated with a timestep of $10^{-4}$. The
stochastic field theory is simulated at $N_{p}=10^{8}$ with a timestep of
$2\times 10^{-6}$. Fig 5 shows the coefficients extracted for the zeroth
moment with Eqn 28 from all three models. The Langevin equation is consistent
with the analytical solution in the full range of $\tilde{k}$. Both the
stochastic field theory and moment equations deviate in the fitted $\tau_{0}$
from the theory value at large $\tilde{k}$.

To illustrate the reason for the deviation, Fig 6 shows the fitted
coefficients for the zeroth moment derived from the simulations with different
number of spatial bins. In contrast to earlier results, these simulations do
not have any instances of corrections for negative density. Instead, the
influence of the number of bins is to influence the maximum value of
$\tilde{k}$. The maximum $\tilde{k}$ is equal to
$\tilde{k}_{max}=\frac{N_{bins}\pi}{\tilde{L}}$. The fitted time constant
matches with the theoretical result for small $\tilde{k}$, but deviates at the
maximum $\tilde{k}$. The time constant is nearly a factor of two larger than
the theory at the maximum $\tilde{k}$. This is consistently true as the number
of bins changes the value of the maximum $\tilde{k}$. This is likely due to
the finite difference derivative approximation applied to the stochastic
concentration. Under the different conditions, the simulation system with
coarser discretization, which has fewer and larger bins, is less accurate.

Fig 7 is the first moment extraction from the same simulations. Such deviation
effect is also shown in the first moment results. The second moment extraction
(Fig 8) is shown for smaller wavevector up to $\tilde{k}\approx 4$. This is
caused by the nature of the second moment autocorrelation function, which
contains two distinct exponential functions in Eqn 27. Accurately extracting
the coefficients from fits is challenging for larger $\tilde{k}$ when the two
timescales are very similar. To reduce the number of variables needed to fit,
the autocorrelation of the second moment has been divided by the
autocorrelation of the zeroth moment as

$\frac{C_{2}}{C_{0}}=a_{2/0}e^{-\frac{\Delta\tilde{t}}{\tau_{2/0}}}+b_{2/0}$
(29)

This reduces the fitting function to a constant plus an exponential function
which contains only three parameters. However, the fitting remains
challenging. For higher $\tilde{k}$, the zeroth moment autocorrelation
function $C_{0}$ decays in lag-time quickly, which leads to dividing by a
small value that is comparable in size to the statistical error.

##  V Discussion and Outlook

In a dilute polymer system, Langevin equations can be used to track the
position of each coarse-grained polymer segment. A semi-dilute or concentrated
system makes the Langevin equations less effective and more computationally
costly. The stochastic field theory (SFT) approach is applicable in these
systems, while being less accurate at low concentrations. The corresponding
moment equations can be accurate at a range of concentrations in polymeric
systems. For all these models, the spatial discretization of the system
determines the accuracy at small length scales. Even though the error at the
smallest length scale is always about $100\%$, the accuracy at larger length
scales improves as the number of grid points in the system increases. A
typical trade-off occurs in which more grid points improves accuracy but with
higher computational cost. The full SFT has an additional trade-off. Large
number of grid points leads to a small number of polymers in each region on
average. This can cause an error from keeping track of rare conformations with
large stretch.

Although the numerical calculations of the field theory were validated in a
simplified situation with an analytical solution, the full stochastic field
theory and moment equations are designed to work in broader conditions. They
can include anything that can be included in the Langevin equations for bead
positions. For example, finite extensibility can be included via the finitely
extensibile nonlinear elastic (FENE) model. In this situation, the closure in
the moment equations could lead to larger errors than the situations shown
here. The use of higher dimensional systems is necessary for systems pushed
out of equilibrium by flow. The moment equations will be critical in these
cases to reduce the computational cost of the higher dimensional systems.
Finally this approach can be used for understanding the influence of
interactions. This could include Flory-Huggins type interactions for a polymer
and solvent field, interactions within block copolymers, or interactions among
polyelectrolytes including fluctuating ion fields.

## Author Declarations

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Shangren Zhu: Formal Analysis (equal); Software (lead); Writing - original
draft (lead); Writing - review & editing (equal). Patrick T. Underhill:
Conceptualization (lead); Formal Analysis (equal); Funding acquisition (lead);
Writing - review & editing (equal).

## Data Availability

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

###### Acknowledgements.

This work was supported by the donors of ACS Petroleum Research Fund under
Grant 62649-ND9.

##  Appendix A Discretization of stochastic field theory

Consider a two dimensional space with spring vector $\tilde{Q}$ discretized
with labels $i$ and $\widetilde{r_{c}}$ discretized with labels $j$. Boundary
conditions need to be set for both $\tilde{Q}$ and $\widetilde{r_{c}}$ to
construct the simulation space. A periodic boundary condition is applied to
$\widetilde{r_{c}}$, and a no-flux boundary condition is applied to
$\tilde{Q}$ on $[-\tilde{Q}_{Max},\tilde{Q}_{Max}]$. $\tilde{Q}_{Max}$ is set
to $4$ for this work.

For the number density computed at point $i,j$ in such space, the noise terms
corresponding to fluxes in $\tilde{Q}$ and $\widetilde{r_{c}}$ are solved at
half points $(i+\frac{1}{2},j)$, $(i-\frac{1}{2},j)$, $(i,j+\frac{1}{2})$, and
$(i,j-\frac{1}{2})$. Therefore, we define the discretized number density as
the average value around a grid point as
$\tilde{\psi}_{i,j}=\frac{1}{\Delta\widetilde{r_{c}}\Delta\tilde{Q}}\int_{\tilde{r}_{c,j-\frac{1}{2}}}^{\tilde{r}_{c,j+\frac{1}{2}}}\int_{\tilde{Q}_{i-\frac{1}{2}}}^{\tilde{Q}_{i+\frac{1}{2}}}\tilde{\psi}d\tilde{Q}d\widetilde{r_{c}}$.

We get the discretized equation by performing this average over the full
equation. Replacing the density at the boundaries between grid points by the
average of the neighboring grid points gives the stochastic theory as

$\begin{split}d\tilde{\psi}_{i,j}=&\frac{d\tilde{t}}{(\Delta\tilde{r}_{c})^{2}}[\tilde{\psi}_{i+1,j}-2\tilde{\psi}_{i,j}+\tilde{\psi}_{i-1,j}]+\frac{2d\tilde{t}}{\Delta\tilde{Q}}[(\tilde{Q}\tilde{\psi})_{i+1,j}-(\tilde{Q}\tilde{\psi})_{i-1,j})]\\\
&+\frac{4d\tilde{t}}{(\Delta\tilde{Q})^{2}}[\tilde{\psi}_{i+1,j}-2\tilde{\psi}_{i,j}+\tilde{\psi}_{i-1,j}]\\\
&+\sqrt{\frac{2\tilde{L}}{N_{p}}}\frac{1}{\Delta\tilde{r}_{c}}[\tilde{U}_{r,i,j+\frac{1}{2}}-\tilde{U}_{r,i,j-\frac{1}{2}}]+\sqrt{\frac{8\tilde{L}}{N_{p}}}\frac{1}{\Delta\tilde{Q}}[\tilde{U}_{Q,i+\frac{1}{2},j}-\tilde{U}_{Q,i-\frac{1}{2},j}]\end{split}$
(30)

where the Gaussian noise term $\tilde{U}_{{r},i,j+\frac{1}{2}}$ is defined as

$\tilde{U}_{r,i,j+\frac{1}{2}}=\frac{1}{\Delta\tilde{r}_{c}\Delta\tilde{Q}}\int_{\tilde{r}_{c,j}}^{\tilde{r}_{c,j+1}}\int_{\tilde{Q}_{i-\frac{1}{2}}}^{\tilde{Q}_{i+\frac{1}{2}}}\sqrt{\tilde{\psi}}d\tilde{U}_{r_{c}}d\tilde{Q}d\tilde{r}_{c}$
(31)

and the other noise terms are defined similarly. They each have a mean of zero
and variance of

$\langle\tilde{U}_{r_{c},i,j+\frac{1}{2}}\tilde{U}_{r,i,j+\frac{1}{2}}\rangle=\frac{d\tilde{t}}{2\Delta\tilde{r}_{c}\Delta\tilde{Q}}(\tilde{\psi}_{i,j+1}+\tilde{\psi}_{i,j})$
(32)

The error in the numerical density field is determined by the timestep
$d\tilde{t}$. At $\tilde{Q}$ near $\pm\tilde{Q}_{Max}$, the number density is
close to zero. At any step in time, the random noise can make this density
become negative. Even with small timesteps $(d\tilde{t}<10^{-5})$, the
probability of producing negative density is not zero. One approach is to
instead use non-Gaussian stochastic fluxes that still have the correct moments
but do not allow more polymers to leave a location than are at that location.
It is more computationally efficient to correct the error by manually setting
the negative density to zero and renormalizing the density field to obey

$\frac{1}{\tilde{L}}\int_{0}^{\tilde{L}}\int_{-\tilde{Q}_{Max}}^{\tilde{Q}_{Max}}\tilde{\psi}d\tilde{Q}d\tilde{r}_{c}=1$
(33)

The error associated with this correction is the reason for the deviation from
the correct answer at small $N_{p}$ in Figs 1 to 4. In practice, the SFT will
be used in situations where the negative densities are unlikely to occur.

##  Appendix B Discretization of moment equations

The moment equations are also discretized with a similar method, but the
moments only depend on $\tilde{r}_{c}$. Averaging over a grid point gives

$d\tilde{\mu}_{0,j}=\frac{d\tilde{t}}{(\Delta\tilde{r}_{c})^{2}}[\tilde{\mu}_{0,j+1}-2\tilde{\mu}_{0,j}+\tilde{\mu}_{0,j-1}]+\sqrt{\frac{2\tilde{L}}{N_{p}}}\frac{1}{\Delta\tilde{r}_{c}}(g_{c,j+\frac{1}{2}}^{(1)}-g_{c,j-\frac{1}{2}}^{(1)})$
(34)
$d\tilde{\mu}_{1,j}=\frac{d\tilde{t}}{(\Delta\tilde{r}_{c})^{2}}[\tilde{\mu}_{1,j+1}-2\tilde{\mu}_{1,j}+\tilde{\mu}_{1,j-1}]-4d\tilde{t}\tilde{\mu}_{1,j}+\sqrt{\frac{2\tilde{L}}{N_{p}}}\frac{1}{\Delta\tilde{r}_{c}}(g_{c,j+\frac{1}{2}}^{(2)}-g_{c,j-\frac{1}{2}}^{(2)})-\sqrt{\frac{8\tilde{L}}{N_{p}}}g_{q,j}^{(1)}$
(35)
$\begin{split}d\tilde{\mu}_{2,j}=&\frac{d\tilde{t}}{(\Delta\tilde{r}_{c})^{2}}[\tilde{\mu}_{2,j+1}-2\tilde{\mu}_{2,j}+\tilde{\mu}_{2,j-1}]-8d\tilde{t}\tilde{\mu}_{2,j}+8d\tilde{t}\tilde{\mu}_{0,j}\\\
&+\sqrt{\frac{2\tilde{L}}{N_{p}}}\frac{1}{\Delta\tilde{r}_{c}}(g_{c,j+\frac{1}{2}}^{(3)}-g_{c,j-\frac{1}{2}}^{(3)})-2\sqrt{\frac{8\tilde{L}}{N_{p}}}g_{q,j}^{(2)}\end{split}$
(36)

Because the noise terms result from integrating stochastic fluxes, they can be
correlated. For the noise terms $g_{q}$ which are associated with the spring
vector, they are represented with independent Gaussian noises as

$\begin{bmatrix}g_{q,j}^{(1)}\\\
g_{q,j}^{(2)}\end{bmatrix}=\sqrt{\frac{1}{\Delta\tilde{r}_{c}}}\tilde{B}_{j}\begin{bmatrix}d{W_{q,j}}^{(1)}\\\
d{W_{q,j}}^{(2)}\end{bmatrix}$ (37)

where $\tilde{B}_{j}$ is the coefficient matrix. The coefficient matrix
$\tilde{B}_{j}$ is related to the covariance of the noise, and is given by

$\tilde{B}_{j}\tilde{B}_{j}^{T}=\begin{bmatrix}\tilde{\mu}_{0,j}&\tilde{\mu}_{1,j}\\\
\tilde{\mu}_{1,j}&\tilde{\mu}_{2,j}\end{bmatrix}$ (38)

From the moments, the coefficient matrix $\tilde{B}_{j}$ can be analytically
calculated using Cholesky decomposition.

The noise terms $g_{c}$ are associated with the center of mass position and
are evaluated with a similar procedure at the halfway point between grid
points. The noise is written as

$\begin{bmatrix}g_{c,{j+\frac{1}{2}}}^{(1)}\\\ g_{c,{j+\frac{1}{2}}}^{(2)}\\\
g_{c,{j+\frac{1}{2}}}^{(3)}\end{bmatrix}=\sqrt{\frac{1}{\Delta\tilde{r}_{c}}}\tilde{G}_{j+\frac{1}{2}}\begin{bmatrix}d{W_{c,{j+\frac{1}{2}}}}^{(1)}\\\
d{W_{c,{j+\frac{1}{2}}}}^{(2)}\\\ d{W_{c,{j+\frac{1}{2}}}}^{(3)}\end{bmatrix}$
(39)

with coefficient matrix given by

$\tilde{G}_{j+\frac{1}{2}}\tilde{G}_{j+\frac{1}{2}}^{T}=\begin{bmatrix}\tilde{\mu}_{0,j+\frac{1}{2}}&\tilde{\mu}_{1,j+\frac{1}{2}}&\tilde{\mu}_{2,j+\frac{1}{2}}\\\
\tilde{\mu}_{1,j+\frac{1}{2}}&\tilde{\mu}_{2,j+\frac{1}{2}}&\tilde{\mu}_{3,j+\frac{1}{2}}\\\
\tilde{\mu}_{2,j+\frac{1}{2}}&\tilde{\mu}_{3,j+\frac{1}{2}}&\tilde{\mu}_{4,j+\frac{1}{2}}\end{bmatrix}$
(40)

The coefficient matrix depends on the third and fourth moments. Since we only
simulate the zeroth to the second moment, these two further moments are
obtained with a closure approximation. We first define scaled moments as
$M_{n}=\frac{\tilde{\mu}_{n}}{\tilde{\mu}_{0}}$. The scaled moments can be
written in terms of cumulants $K_{n}$ as

$M_{1}=K_{1}$ (41) $M_{2}=K_{2}+K_{1}^{2}$ (42)
$M_{3}=K_{3}+3K_{2}K_{1}+K_{1}^{3}$ (43)
$M_{4}=K_{4}+4K_{3}K_{1}+3K_{2}^{2}+6K_{2}K_{1}^{2}+K_{1}^{4}$ (44)

As a closure, we set $K_{n}=0$ for $n\geq 3$. From this condition, $M_{3}$ and
$M_{4}$ can be determined in terms of $M_{1}$ and $M_{2}$, which are
associated with the simulated zeroth to second moments. With the condition of
$M_{2}>M_{1}^{2}$, which is satisfied in most cases, the matrix
$\tilde{G}_{j+\frac{1}{2}}$ can be calculated with Cholesky decomposition.

## References

  * Larson [1992] R. Larson, “Flow-induced mixing, demixing, and phase transitions in polymeric fluids,” Rheologica Acta 31, 497–520 (1992).
  * Liberatore  _et al._ [2009] M. W. Liberatore, N. B. Wyatt, M. Henry, P. L. Dubin, and E. Foun, “Shear-induced phase separation in polyelectrolyte/mixed micelle coacervates,” Langmuir 25, 13376–13383 (2009).
  * Korolkovas [2018] A. Korolkovas, “Giant fluctuations in sheared viscoelastic fluids emerging in a mesoscale simulation,” arXiv preprint arXiv:1808.06381 (2018).
  * Stieger and Richtering [2003] M. Stieger and W. Richtering, “Shear-induced phase separation in aqueous polymer solutions: Temperature-sensitive microgels and linear polymer chains,” Macromolecules 36, 8811–8818 (2003).
  * Hashimoto, Murase, and Ohta [2009] T. Hashimoto, H. Murase, and Y. Ohta, “Shear-induced phase separation of entangled polymer solutions: Formation of optically anisotropic strings as precursor structures of shish-kebab,” in  _Macromolecular symposia_ , Vol. 279 (Wiley Online Library, 2009) pp. 88–95.
  * De Kruif, Weinbreck, and de Vries [2004] C. G. De Kruif, F. Weinbreck, and R. de Vries, “Complex coacervation of proteins and anionic polysaccharides,” Current opinion in colloid & interface science 9, 340–349 (2004).
  * Guillen  _et al._ [2011] G. R. Guillen, Y. Pan, M. Li, and E. M. Hoek, “Preparation and characterization of membranes formed by nonsolvent induced phase separation: a review,” Industrial & Engineering Chemistry Research 50, 3798–3817 (2011).
  * Garcia  _et al._ [2020] J. U. Garcia, T. Iwama, E. Y. Chan, D. R. Tree, K. T. Delaney, and G. H. Fredrickson, “Mechanisms of asymmetric membrane formation in nonsolvent-induced phase separation,” ACS Macro Letters 9, 1617–1624 (2020).
  * Kumar and Larson [2001] S. Kumar and R. G. Larson, “Brownian dynamics simulations of flexible polymers with spring–spring repulsions,” The Journal of Chemical Physics 114, 6937–6941 (2001).
  * Young, Marvin, and Sing [2018] C. D. Young, M. Marvin, and C. E. Sing, “Conformationally averaged iterative brownian dynamics simulations of semidilute polymer solutions,” The Journal of chemical physics 149, 174904 (2018).
  * Fredrickson, Ganesan, and Drolet [2002] G. H. Fredrickson, V. Ganesan, and F. Drolet, “Field-theoretic computer simulation methods for polymers and complex fluids,” Macromolecules 35, 16–39 (2002).
  * Fredrickson  _et al._ [2006] G. Fredrickson _et al._ ,  _The equilibrium theory of inhomogeneous polymers_ , 134 (Oxford University Press on Demand, 2006).
  * Singh [2014] N. Singh, “On the failure of mean-field theories near a critical point,” arXiv preprint arXiv:1402.6841 (2014).
  * Chavanis [2008] P.-H. Chavanis, “Hamiltonian and brownian systems with long-range interactions: V. stochastic kinetic equations and theory of fluctuations,” Physica A: Statistical Mechanics and its Applications 387, 5716–5740 (2008).
  * Baeurle, Efimov, and Nogovitsin [2006] S. A. Baeurle, G. Efimov, and E. Nogovitsin, “Calculating field theories beyond the mean-field level,” EPL (Europhysics Letters) 75, 378 (2006).
  * Müller and Smith [2005] M. Müller and G. D. Smith, “Phase separation in binary mixtures containing polymers: A quantitative comparison of single-chain-in-mean-field simulations and computer simulations of the corresponding multichain systems,” Journal of Polymer Science Part B: Polymer Physics 43, 934–958 (2005).
  * Tree  _et al._ [2017] D. R. Tree, K. T. Delaney, H. D. Ceniceros, T. Iwama, and G. H. Fredrickson, “A multi-fluid model for microstructure formation in polymer membranes,” Soft Matter 13, 3013–3030 (2017).
  * Cromer, Fredrickson, and Gary Leal [2017] M. Cromer, G. H. Fredrickson, and L. Gary Leal, “Concentration fluctuations in polymer solutions under mixed flow,” Journal of Rheology 61, 711–730 (2017).
  * Dean [1996] D. S. Dean, “Langevin equation for the density of a system of interacting langevin processes,” Journal of Physics A: Mathematical and General 29, L613 (1996).
  * Lefevre and Biroli [2007] A. Lefevre and G. Biroli, “Dynamics of interacting particle systems: stochastic process and field theory,” Journal of Statistical Mechanics: Theory and Experiment 2007, P07024 (2007).
  * Tailleur and Cates [2008] J. Tailleur and M. Cates, “Statistical mechanics of interacting run-and-tumble bacteria,” Physical review letters 100, 218103 (2008).
  * Qian, Kramer, and Underhill [2017] Y. Qian, P. R. Kramer, and P. T. Underhill, “Stochastic kinetic theory for collective behavior of hydrodynamically interacting active particles,” Physical Review Fluids 2, 043104 (2017).

Figure 1: Coefficients fitting of zeroth moment autocorrelation $C_{0}$ with
various polymer number ($N_{p}$) at $\tilde{k}=\frac{2\pi}{\tilde{L}}$ for the
Langevin equation (blue circles), moment equations (red squares), and
stochastic fluctuation theory (yellow diamonds). (a) prefactor $a_{0}$ and (b)
time constant $\tau_{0}$. The dashed line represents the theoretical value.

Figure 2: Coefficients fitting of first moment autocorrelation $C_{1}$ with
various polymer number ($N_{p}$) at $\tilde{k}=0$ for the Langevin equation
(blue circles), moment equations (red squares), and stochastic fluctuation
theory (yellow diamonds). (a) prefactor $a_{1}$ and (b) time constant
$\tau_{1}$. The dashed line represents the theoretical value.

Figure 3: Coefficients fitting of second moment autocorrelation $C_{2}$ with
various polymer number ($N_{p}$) at $\tilde{k}=0$ for the Langevin equation
(blue circles), moment equations (red squares), and stochastic fluctuation
theory (yellow diamonds). (a) prefactor $a_{2}$ and (b) time constant
$\tau_{2}$. The dashed line represents the theoretical value.

Figure 4: Coefficients fitting of zeroth moment autocorrelation $C_{0}$ with
various polymer numbers ($N_{p}$) at $\tilde{k}=\frac{2\pi}{\tilde{L}}$ for
stochastic fluctuation theory. (a) prefactor $a_{0}$ and (b) time constant
$\tau_{0}$. The symbols denote a different number of bins used in the SFT,
with $N_{bins}=30$ (blue circles), $N_{bins}=50$ (red squares), and
$N_{bins}=100$ (green diamonds). The dashed line represents the theoretical
value.

Figure 5: Coefficients fitting of zeroth moment autocorrelation $C_{0}$ with
various wavevectors $\tilde{k}$ for the Langevin equation (blue circles),
moment equations (red squares), and stochastic fluctuation theory (yellow
diamonds). (a) prefactor $a_{0}$ and (b) time constant $\tau_{0}$. The dashed
line represents the theoretical value.

Figure 6: Coefficients fitting of zeroth moment autocorrelation $C_{0}$ with
various wavevectors $\tilde{k}$ at $N_{p}=10^{8}$. (a) prefactor $a_{0}$ and
(b) time constant $\tau_{0}$. The symbols denote a different number of bins
used in the moment equations model, with $N_{bins}=30$ (blue circles),
$N_{bins}=50$ (red squares), and $N_{bins}=100$ (green diamonds). The dashed
line represents the theoretical value.

Figure 7: Coefficients fitting of first moment autocorrelation $C_{1}$ with
various wavevectors $\tilde{k}$ at $N_{p}=10^{8}$. (a) prefactor $a_{1}$ and
(b) time constant $\tau_{1}$. The symbols denote a different number of bins
used in the moment equations model, with $N_{bins}=30$ (blue circles),
$N_{bins}=50$ (red squares), and $N_{bins}=100$ (green diamonds). The dashed
line represents the theoretical value.

Figure 8: Coefficients fitting of the ratio of second moment autocorrelation
with zeroth moment autocorrelation $\frac{C_{2}}{C_{0}}$ with various
wavevectors $\tilde{k}$ at $N_{p}=10^{8}$. (a) prefactor $a_{2/0}$ , (b) time
constant $\tau_{2/0}$ and (c) constant $b_{2/0}$ . The symbols denote a
different number of bins used in the moment equations model, with
$N_{bins}=30$ (blue circles), $N_{bins}=50$ (red squares), and $N_{bins}=100$
(green diamonds). The dashed line represents the theoretical value.


# Proper Orthogonal Decomposition Analysis and Modelling of the wake deviation
behind a squareback Ahmed body

Bérengère Podvin  LIMSI, CNRS, Université Paris-Saclay  Stéphanie Pellerin
LIMSI, CNRS, Université Paris-Saclay  Yann Fraigneau  LIMSI, CNRS, Université
Paris-Saclay  Antoine Evrard  ENSTA  Olivier Cadot  University of Liverpool

###### Abstract

We investigate numerically the 3-D flow around a squareback Ahmed body at
Reynolds number $Re=10^{4}$. Proper Orthogonal Decomposition (POD) is applied
to  a symmetry-augmented database in order to describe and model the flow
dynamics. Comparison with experiments at a higher Reynolds number in a plane
section of the near-wake at mid-height shows that the simulation captures
several features of the experimental flow, in particular the antisymmetric
quasi-steady deviation mode. 3-D POD analysis allows us to classify the
different physical processes in terms of mode contribution to the kinetic
energy over the entire domain. It is found that the dominant fluctuating mode
on the entire domain corresponds to the 3-D quasi-steady wake deviation, and
that its amplitude is well estimated from 2-D near-wake data. The next most
energetic flow fluctuations consist of vortex shedding and bubble pumping
mechanisms. It is found that the amplitude of the deviation is negatively
correlated with the intensity of the vortex shedding in the spanwise direction
and the suction drag coefficient. Finally, we find that despite the slow
convergence of the decomposition, a POD-based low-dimensional model reproduces
the dynamics of the wake deviation observed experimentally, as well as the
main characteristics of the global modes identified in the simulation.

##  I Introduction

A surprisingly generic feature of the flow around  symmetric bodies at high
Reynolds numbers is the presence of permanent symmetry-breaking structures in
the wake. These have been observed for the sphere [9], bullet [24] , the flat
plate [2], as well as academic models for ground vehicles such as the Ahmed
body [10] and the Windsor body [18]. The origin and dynamics of these
structures has not been entirely elucidated, although they have been the
object of several experimental and numerical studies. They appear to be
connected to the first bifurcation observed at much lower Reynolds number (see
[7], [8]). At higher Reynolds numbers, rapid switches between different quasi-
stable states can be observed, as described by [22] and [1].

In the remainder of the paper, we will focus on the Ahmed body. Grandemange et
al. [10] have established that the appearance of bi- or multi-stable states
for the flow around a squareback Ahmed body depends on the value of the ground
clearance and the aspect ratio of the body base. Pasquetti and Peres [16]
carried out the first numerical simulation of the Ahmed body that was able to
reproduce the steady wake deviation. However switches were not observed in
Pasquetti and Peres’s simulation, since the typical time separating two
switches is on the order of 1000 $U/H$ where $U$ is the incoming flow speed
and $H$ the body characteristic height (see [14]). In more recent work Dalla
Longa et al. [4] were able to integrate over sufficiently long times to
capture the switch in the wake deviation. However they could only observe one
switch, due to the still longer time scale separating two switches. They also
applied Proper Orthogonal Decomposition and Dynamic Mode Decomposition (DMD,
[25]) to extract the dynamics of the flow. Interaction of the near-wake
recirculation zone with the surrounding shear layers is expected to play a
part in the occurence of switches. They proposed that the switch is triggered
by large hairpin vortices.

Volpe et al. [28] carried out a spectral analysis of the wake behind the
squareback Ahmed body. Three low frequencies were identified in their
experimental data: vortex-shedding modes with Strouhal numbers of 0.13 and
0.19 in the far wake, and wake pumping motion at a Strouhal number of 0.08 in
the recirculation zone. Pavia et al. [17] have applied POD to both the
experimental pressure and velocity field of a Windsor body to identify
structures corresponding to the bi-stable modes and the dynamics of the
switch. They identified pumping motion in the wake with a bi-stable vortical
structure in the streamwise direction and derived a phase-averaged model that
describes the switch between the quasi-stable states. They confirmed the
observation made by Evrard et al. [5] and Cadot et al. [2] that the symmetric
state corresponds to a lower drag.

Drag reduction through control of the flow asymmetry has been the object of
several passive and active control strategies [14], [29], [1], [13], [27],
[6], [23]. An important question is to determine how the different structures
present in the flow contribute to the drag. A low-dimensional description and
modelling of the large-scale structures could be beneficial for understanding,
predicting and ultimately controlling the flow dynamics. Reduced-order models
were developed by [1], [22] in order to control the deviation of the wake.

In the present work, we aim to provide a large-scale description of the flow
by applying Proper Orthogonal Decomposition to a numerical investigation of
the flow around a squareback Ahmed body at a Reynolds number of $Re=10^{4}$.
The spanwise to vertical aspect ratio of the body is 1.18, so that bistability
corresponds to two asymmetric states which are symmetric through reflection of
the vertical mid-plane. Due to the long time scales separating switches, the
change in the wake deviation could not be captured by the simulation. However,
reflection symmetry was enforced artificially, so that we provide a
description of the structures of the flow as a superposition of reflection-
symmetric and reflection-antisymmetric modes.

The paper is organized as follows: section 2 presents the 3-D numerical
configuration, while section 3 describes the specifics of Proper Orthogonal
Decomposition. Section 4 presents a 2-D comparison of the numerical simulation
and an experimental configuration corresponding to the same geometry but a
higher Reynolds number, section 5 presents a 3-D POD analysis of the
structures in the full configuration; a POD-based low-dimensional model is
constructed in section 6 and compared with experimental and numerical results.
A conclusion is given in section 7.

##  II Numerical configuration

Figure 1 presents the numerical configuration. The dimensions of the
squareback Ahmed body are the same as the experimental configuration of Evrard
et al. [5] i.e. $L=1.124m,H=0.297m,W=0.35m$. The ground clearance  (distance
from the body to the lower boundary of the domain) is  0.3H in the simulation.
The Reynolds number based on the incoming velocity $U$, body height $H$ and
fluid viscosity is $10^{4}$. The foremost and upper part body defines the
reference position $(x=0,z=0)$. The plane $y=0$ corresponds to the mid-plane
of the body. The domain extent in the streamwise (x), spanwise (y) and
vertical (z) directions is $[-1,10]H\times[-2,2]H\times[-1.3,1.2]H$.

The code used is SUNFLUIDH, an in-house code developed at LIMSI based on a
second-order finite volume approach which has been described for instance in
[20]. The temporal discretization is based on a second-order backward Euler
scheme. Diffusion terms are treated implicitly and convective terms are solved
with an Adams-Bashforth scheme. The Poisson equation for the computation of
pressure field is solved iteratively. We use $(512\times 256\times 256)$ grid
points in respectively the longitudinal direction $x$, the spanwise direction
$y$, and the vertical direction $z$. The Cartesian grid is refined close to
surfaces. Periodic boundary conditions are used in the spanwise direction. No-
slip velocity boundary conditions on the Ahmed body  and on the ground are
implemented by adjusting the size of the loops. The simulation was initialized
from a uniform condition. About 100 time units based on the upstream velocity
$U$ and body height $H$ were necessary for the flow to develop and statistical
convergence to be reached. We note that all times will be expressed in those
units in the remainder of the paper. All lengths will be made nondimensional
with the body height $H$.

Figure 1: Numerical configuration. The Ahmed body and the recirculation zone
behind it are idenfied by a surface of zero streamwise velocity. The
streamwise, spanwise and wall-normal positions are referred to as $x$, $y$,
$z$.  The length, width and height of the body are respectively $L$, $W$ and
$H$. The origin of the axes is taken at the top and foremost position of the
Ahmed body in the vertical symmetry plane.

##  III Proper Orthogonal Decomposition

The main tool of analysis used in this paper is Proper Orthogonal
Decomposition (POD) [15]. The field on a domain $\Omega$ is written as a
superposition of spatial modes

$\underline{u}(\underline{x},t)=\sum_{n}\tilde{a}_{n}(t)\underline{\phi}_{n}(\underline{x})$
(1)

where the modes are orthogonal (and can be made orthonormal), i.e

$\int\underline{\phi}_{n}(\underline{x}).\underline{\phi}_{m}(\underline{x})d\underline{x}=\delta_{nm}.$

The modes can be ordered by decreasing energy
$\lambda_{1}\geq\lambda_{2}\geq\ldots\geq\lambda_{n}=<\tilde{a}_{n}\tilde{a}_{n}>$,
where $<.>$ represents a time average. The amplitudes $\tilde{a}_{n}$ can be
obtained by projection of the velocity field onto the spatial modes:

$\tilde{a}_{n}(t)=\int_{\Omega}\underline{u}(\underline{x},t).\underline{\phi}_{n}(\underline{x})d\underline{x}.$
(2)

In the remainder of the paper we will consider normalized amplitudes
$a_{n}=\tilde{a}_{n}/\sqrt{\lambda_{n}}.$

In all that follows POD is implemented following the method of snapshots [26]
which is based on computing the autocorrelation between the different samples
of the field obtained at times $t_{i},i=1,\ldots,N$:

$C_{ij}=\frac{1}{N}\int_{\Omega}\underline{u}(\underline{x},t_{i}).\underline{u}(\underline{x}^{\prime},t_{j})d\underline{x},$

and extracting the eigenvalues $\lambda_{n}$ and temporal eigenvectors
$A_{in}=a_{n}(t_{i})$ such that

$CA=\lambda A.$ (3)

The spatial modes can then be reconstructed using

$\tilde{\underline{\phi}}_{n}(\underline{x})=u(\underline{x},t_{i})A_{in}$ (4)

and renormalizing
${\underline{\phi}}_{n}(\underline{x})=\frac{1}{N_{n}}\tilde{\underline{\phi}}_{n}(\underline{x})$
where

$N_{n}^{2}=\int\tilde{\underline{\phi}}_{n}(\underline{x}).\tilde{\underline{\phi}}_{n}(\underline{x})d\underline{x}.$
(5)

Different spatial domains, as well as different quantities will be considered
in the decomposition. In section 4, we will limit the analysis to 2-D velocity
fields in the near-wake region in order to match the data of Evrard et al.
[5]. In section 5, we will apply POD to the full 3-D velocity field to the
full numerical domain.

A key ingredient of the procedure is the definition of the data set. Cross-
sections of the time-averaged streamwise velocity in the mid-height plane
$z=-0.5$ and the vertical mid-span plane $y=0$ are represented in figure 2. We
observe a steady deviation of the wake in agreement with previous observations
[16]. This deviation is breaking the reflection symmetry with respect to the
vertical mid-plane. However the equations are symmetric \- for each flow
realization, the flow obtained by symmetry with respect to the vertical
midplane is also a possible solution. Following recommendations in [11], an
enlarged dataset enforcing the statistical symmetry could then be created as
follows: for each snapshot of the original dataset, a symmetrized snapshot
corresponding to the image of the snapshot by the reflection symmetry was
created. The size of the new data set was therefore twice that of the original
one. By construction, POD modes are thus either symmetric or antisymmetric
with respect to the vertical midplane $y=0$,  so that the amplitude of a POD
mode for a given snapshot of the original dataset is either identical or
opposite to that on the symmetrized snapshot. This allows us to break down
flow patterns into symmetric and antisymmetric components.

|
---|---
Figure 2:  Streamwise velocity contours of the mean flow; left) horizontal
plane at mid-height $z=-0.5$; right) vertical mid-plane at $z=0$.

We emphasize that individual POD modes are different from coherent structures.
A coherent structure typically corresponds to local patterns identified in a
realization. POD modes are defined over the full domain chosen for
decomposition, and every single realization contains a combination of modes. A
coherent structure will therefore correspond to a combination of a few POD
modes, which may be restricted to a portion of a spatial domain. As a
consequence of the data enlargement procedure we have adopted, the mean flow
of the simulation does not correspond to one single mode but to the sum of the
first two modes representing respectively the symmetric and the antisymmetric
part of the mean flow.

##  IV 2-D POD: Comparison with the experiment

In this section we apply POD analysis in the near-wake to the same variables
defined over the same domain in both the experiment and the numerical
simulation, i.e. the streamwise and the spanwise velocity components $u$ et
$v$ defined over the domain $L\leq x\leq L+1.2H$, $-0.6\leq y\leq 0.6H$.
Details are indicated in table 1. The main differences are
(i) the Reynolds numbers considered, as $Re=4$ $10^{5}$ in Evrard’s experiment
and $Re=10^{4}$ in our numerical simulation
(ii) the time resolution - which is $25H/U$ in the experiment versus 0.5H/U in
the simulation.

Type | Re | $\Delta T$ | snapshot number (original data set)
---|---|---|---
Simulation | $10^{3}$ | 0.5 | 150
Experiment |  $4$ $10^{5}$ | 25.25 | 400
Table 1:  Comparison between model and experiment

Results for the POD spectrum are shown in figure 3 and show a good agreement
between the experiment and the simulation. The second mode represents about
35% of the total fluctuating horizontal kinetic energy in the wake
$\sum_{p\geq 2}\lambda_{p}$ and the first eight modes capture more than 50%.
Figure 4 compares the first four POD modes for the experiment and the
simulation. The first mode represents a cross-section of a toroidal-like
structure constituting the recirculation bubble. The second mode corresponds
to a single vortical structure located in the recirculation bubble close to
the rear of the body and sweeping fluid from one side of the wake to the
other. It represents the wake deviation. Overall a good agreeement is observed
between all the modes found in the experiment and those found in the
simulation.

Figure 5 shows for both the experiment and the simulation the corresponding
amplitudes (normalized by the square root of the eigenvalue) of the modes
represented in figure 4. The amplitudes of the first two modes in the
simulation (figure 5 right) display small oscillations around a constant
positive value. In the experiment (figure 5 left), the amplitude of the first
mode oscillates near a constant positive value, while the amplitude of the
second mode changes sign several times, which corresponds to a switch in the
wake asymmetry.

We can also see from figure 4 that mode 3 is symmetric and mode 4 is
antisymmetric. Mode 3 consists of a longitudinal converging (or diverging,
depending on the sign of the amplitude) motion at the extremity of the
recirculation bubble, which could be associated with wake pumping i.e
successive enlargment and shrinking of the recirculation region. Mode 4
consists of three vortical structures - one larger structure extending across
the wake and located close to the rear of the body, and two smaller ones, both
rotating in the opposite direction, on each side of the recirculating bubble.
Its action is therefore to distort the recirculation bubble in the spanwise
direction. Since the PIV results are not resolved in time and the simulation
total time is relatively short, it is difficult to compare directly time
evolutions. However histograms of the amplitudes can be compared in figure 6
(we note that only the amplitudes corresponding to the original set of
snapshots are shown). We can see that there is a relatively good agreement
between the experiment and the simulation.

Figure 3:  Near-wake 2-D POD spectrum in the simulation and experiment. n=1 |  n=2 |  n=3 |  n=4
---|---|---|---
|  |  |
|  |  |
Figure 4:  2-D wake POD modes 1 to 4; top row: simulation; bottom row: experiment.  |
---|---
|
Figure 5: Amplitudes of 2-D POD modes 1 to 4 (from top to bottom); left) simulation; right) experiment. n=1 |  n=2 |  n=3 |  n=4
---|---|---|---
|  |  |
|  |  |
Figure 6:  Histogram comparison of 2-D POD modes 1 to 4 (from left to right);
top row: simulation; bottom row: experiment.

The next four order modes in the simulation, shown in figure 7, also present
similarities with the modes in the experiment. Modes 5 and 8 are antisymmetric
and consist of vortical motions respectively dominant on the inner part and
the outer part of the recirculation. Modes 6 and 7 are symmetric and consist
of two counter-rotating vortical structures, extending over the full
recirculation length and located on each side of the wake. Given the
differences in Reynolds number and time resolution, the agreement is
remarkable and suggests that the most energetic structures have common
features over a wide range of Reynolds numbers.

The spectral content of the temporal coefficients is presented in figure 8.
All modes are characterized by low frequencies. The two red and black lines
respectively correspond to the frequencies 0.08 and 0.2. Modes 3 and 7, which
are both symmetric, are characterized by a dominant frequency around 0.08,
while mode 4, 5 and 6, which are antisymmetric, are characterized by a
frequency of 0.2. This frequency, which is associated with vortex shedding,
can be identified most clearly in modes 5 and 6. Mode 8 is characterized by a
mixture of frequencies.

n=5 |  n=6 |  n=7 |  n=8
---|---|---|---
|  |  |
|  |  |
Figure 7:  2-D wake POD modes 5 to 8; top row: simulation; bottom row:
experiment.

Figure 8:  2-D POD amplitude power spectral density $|\hat{a}_{n}|^{2}$ for
modes 3 to 8 - the red and black dashed lines respectively correspond to the
frequencies 0.08 and 0.2.

##  V 3-D POD

To investigate the flow structure and dynamics in more detail, POD is applied
to the full 3-D velocity field over the entire computational domain. The
eigenvalue spectrum is shown in figure 9. The first 3-D mode, which coincides
with the mean flow, is much more energetic than the other modes compared with
2-D analysis (figure 3) since the entire domain contains a large steady
contribution of the kinetic energy outside the wake. Due to the extent of the
spatial domain, the convergence of the fluctuations is much slower than in the
2-D wake measurements. The second mode represents only 8% of the total
fluctuating kinetic energy $\sum_{p\geq 2}\lambda_{p}$,  which is equivalent
to the combined energy of the next eight most energetic modes.  We note that
the modes 3 to 6 have nearly similar eigenvalues. This reflects the fact that
the physical structures are global, three-component modes defined over the
full domain, linking the upstream flow, the four boundary layers along the
body, the near and the far wake. Figure 10 shows the projection of an
instantaneous velocity field on the first ten modes and the difference between
the full field and its projection. It is clear that only the large scales are
captured by the projection. One can see that most of the unresolved modes
(which represent the major portion of the total kinetic energy) are located in
the boundary layers, the shear layers and the far wake.

Figure 9: 3-D POD eigenvalue spectrum.

Figure 10: Streamwise velocity contours on mid-height plane $z=-0.5$; top :
projection of an instantaneous field $\underline{u}$ onto the first ten modes
$u_{recons}=\sum_{n=1}^{10}a_{n}(t)\underline{\phi}(\underline{x})$ where
$a_{n}=\int\underline{u}(\underline{x},t).\underline{\phi}(\underline{x})d\underline{x}.$;
bottom: difference field $|u-u_{recons}|$.

We now investigate the properties of the POD modes. For spatial
characteristics, we represent in figures 11 and 14 the streamwise velocity
contour of the 3-D POD modes in both a horizontal cross-section at mid-height
$z=-0.5$, and in a vertical plane located in an off-center spanwise location
$y=-0.4$ (since asymmetric motions will cancel on the symmetry plane). The
first two spatial modes are shown in figure 11. Although the representation of
the first two POD modes in 2-D and 3-D is different (figures 4 and 11), it can
be shown that the restriction of the first two 3-D POD modes to the cross-
section of the wake is similar to the 2-D POD results, which is not entirely
surprising since they are expected to represent the symmetric and
antisymmetric part of the time-averaged velocity field.

One can see in figure 12 that the temporal evolution of the first two 3-D POD
modes is slightly different from that of the 2-D POD coefficients. Unlike its
2-D counter part (figure 5 right), the amplitude of the first 3-D POD mode is
essentially constant, which shows the statistical convergence of the database.
Moreover, the amplitudes of the second 2-D and 3-D mode present some
discrepancies, but the correlation between the mean deviation coefficient
$a_{2}^{2d}$ and $a_{2}^{3d}$ is 0.6, which indicates that the 2-D
measurements are able to describe reasonably well the evolution of the 3-D
deviation mode. However, although the amplitude of 2-D POD mode 2 changes
sign, the corresponding 3-D POD amplitude $a_{2}^{3d}$ always remains
positive. This means that there are no global switches in the total wake
deviation although a local planar measure may indicate otherwise. In both
cases, their corresponding spectrum indicates fluctuations at low frequencies
thus confirming their very long-time dynamics evolution. In addition to its
permanent asymmetry, the POD mode 2 corresponds to the very low frequency
global mode as reported in [22, 23] for an axisymmetric body and responsible
for the bistable dynamics found in [9] for the Ahmed body. To get a better
insight into the action of the global deviation mode, figure 13 compares
streamlines in the recirculation zone for the symmetrized mean field (mode 1)
and the mean field (modes 1 and 2). We can see that the effect of the global
deviation is to gather the streamlines around one of the base diagonals.

n=1 |  |
---|---|---
n=2 |  |
Figure 11:  Streamwise velocity contours of 3-D POD modes 1 (top row) and 2 (bottom row); left) horizontal section at mid-height $z=-0.5$; right) vertical section at $y=-0.4$. |
---|---
Figure 12:  Left) amplitudes of the first two POD modes; right) power spectral density of the quasi-steady deviation mode $|\hat{a}_{2}|^{2}$ \- the red and black lines respectively correspond to the two frequencies 0.08 and 0.2. |
---|---
Figure 13: Streamlines of the field viewed from downstream left) mode 1
$\sqrt{\lambda_{1}}\underline{\phi}_{1}$; right) combination of mode 1 and
mode 2
$\sqrt{\lambda_{1}}\underline{\phi}_{1}+\sqrt{\lambda_{2}}\underline{\phi}_{2}$
.

Modes 3 and 4 are antisymmetric modes which display a periodicity in the
streamwise direction (figure 14). The time evolution of these modes is quasi-
periodic as shown by the time series and the spectra in figure 15. The
characteristics frequencies lie in the range 0.19-0.24, with a maximum around
0.19. This frequency matches satisfactorily the  Kármán global mode observed
for oscillations in the spanwise direction in (Grandemange et al. [9], Volpe
et al. [28]). Similar observations are made for the symmetric modes 5 and 6 in
figure 14 and the corresponding temporal properties in figure 15. Their
spectra display a maximum around 0.23, that is ascribed to the vortex shedding
(i.e.  Kármán global mode) with oscillations in the vertical direction. The
highest frequency corresponding to the small height and large frequency to the
wider width of the body are in agreement with Grandemange et al. [9], Volpe et
al. [28]. This is a general result for three-dimensional geometries (Kiya and
Abe [12]).

The modes 7 and 8 in figure 14 are symmetric and clearly do not display
streamwise periodicity. The action of mode 7 in the horizontal mid-plane is to
modulate the bubble zone which is either inflated or shrunk, depending on the
sign of $a_{7}$ (figure 15 left). When $a_{7}>0$, figure 14 shows that strong
negative fluctuations are present in the bubble, which delays reattachment,
and vice-versa. The spectra shown in figure 15 (right) show that both modes 7
and 8 are characterized by a strong energetic content at a frequency of 0.08,
which suggests that these modes contribute significantly to the wake pumping,
as characterized by Rigas et al. [24], Volpe et al. [28] and Pavia et al.
[17].

Modes 9 and 10, which are antisymmetric (figure 14), are more difficult to
interpret. As noted earlier, there is no guarantee that an individual POD mode
corresponds to a well-defined physical mechanism. The spectra in figure 15
(right) show that modes 9 and 10 are characterized by a mixture of frequencies
in an intermediate range 0.08-0.2, with a peak for both modes around 0.15. It
is therefore likely that these modes correspond to a superposition of
different physical processes.

As a summary, we find a clear correspondance between the POD modes and the
main global modes that contribute to the wake dynamics reported in the
literature: POD mode 2 is related to the very low frequency asymmetric or
deviation global mode, modes 3 and 4 to vortex shedding with oscillation in
the horizontal direction, modes 5 and 6 to vortex shedding with oscillations
in the vertical direction, and modes 7 and 8 to (symmetric) wake pumping.

|  |
---|---|---
mode |  u contour at $z=-0.5$ |  u contour at $y=-0.4$
n=3 |  |
n=4 |  |
n=5 |  |
n=6 |  |
n=7 |  |
n=8 |  |
n=9 |  |
n=10 |  |
Figure 14:  Streamwise velocity contours of 3-D POD spatial modes 3 to 10 (from top to bottom); left) horizontal section at mid-height $y=-0.5$; right) vertical section at $z=-0.4$. |
---|---
Figure 15:  Left) amplitudes of the modes $a_{n}$, n=3, 5, 7, 9; right) power
spectral density of the 3-D POD mode amplitudes in the simulation
$|\hat{a}_{n}|^{2}$ \- the red and black lines respectively correspond to the
two frequencies 0.08 and 0.2.

The POD decomposition allows to investigate quantitatively the correlations
between the global modes during the wake dynamics. We first compare the
intensity of the deviation global mode given by ($a_{2}^{2}$) with the
magnitude of vortex shedding in the horizontal (spanwise) and the vertical
direction which can respectively be obtained with the sums of the amplitudes
$r_{Kh}^{2}=a_{3}^{2}+a_{4}^{2}$ and $r_{Kv}^{2}=a_{5}^{2}+a_{6}^{2}$. We can
see in figure 16 that minima of $a_{2}$ are associated with high energy in
spanwise vortex shedding. The correlation coefficient between $|a_{2}|$ and
${|a_{3}|^{2}+|a_{4}|^{2}}$ is strongly negative (-0.6). Due to the three-
dimensional nature of vortex shedding, spanwise and vertical shedding motions
are correlated (0.48), so that a negative correlation is also obtained between
$|a_{2}|$ and $a_{5}^{2}+a_{6}^{2}$ (-0.33). The negative correlation is
consistent with the idea that a reduced asymmetry is accompanied by an
increase in vortex shedding, which was observed in control experiments of [1]
and [13].

The next step is to compare the evolution of the POD amplitudes with the base
suction coefficient $C_{B}=-C_{pb}$ where the pressure coefficient $C_{pb}$
corresponds to the integral of the pressure over the base of the body, which
is shown in figure 17 (left). Figure 17 (right) shows that the base suction
coefficient is positively correlated with the amplitude of the wake deviation
$a_{2}^{3d}$ with a positive delay of $\Delta tU/H\sim 2.5$ and with a maximum
correlation coefficient of 0.55 ( we note that the correlation coefficient
with the 2-D deviation mode $a_{2}^{2d}$ is 0.41. with a time delay of 4 time
units). The time delay appears to have some significance as the correlation
without time delay drops to 0.34. This suggests that the variations of the
drag follow those of the mean deviation amplitude. This is in agreement with
the idea that a drag reduction corresponds to a decrease in the quasi-steady
wake asymmetry, as was observed by [2] and [17]. Figure 17 also shows that a
correlation coefficient of -0.55 with a negative delay of $\Delta
tU/H\sim-1.5$ is observed between the pressure and the amplitude of $a_{7}$.
Unlike the previous case, the correlation remains about the same without no
time delay (-0.51), so it is not possible to assign a physical relevance to
the small time delay observed. As seen above, $a_{7}$ is associated with wake
pumping: from figure 14 one can see in the near-wake that when $a_{7}>0$ there
are more negative fluctuations within the zone, so that the size of the
recirculation actually increases. This is consistent with the observation that
the base suction coefficient decreases as the recirculation length increases.

In contrast, the correlation of the drag with the intensity of vortex shedding
appears to be weaker (and negative).

The observations made above suggest the following picture: the flow is
characterized by a quasi-steady wake deviation, vortex shedding modes and low-
frequency wake pumping modes. The drag coefficient depends on the global
(symmetric) size of the bubble, which is associated with wake pumping, as a
longer bubble corresponds to a lower drag. It also depends on the magnitude of
the quasi-steady wake asymmetry: a lower drag corresponds to a decrease of
asymmetry, and to an increase in the intensity of vortex shedding, a negative
correlation of vortex shedding in spanwise and vertical directions with both
the quasi-steady deviation and the base suction coefficient.

|
---|---
Figure 16:  Left) energy of POD modes corresponding to quasi-steady deviation
$a_{2}^{2}$, spanwise ($r_{Kh}^{2}=a_{3}^{2}+a_{4}^{2}$) and vertical
($r_{Kv}^{2}=a_{5}^{2}+a_{6}^{2}$) vortex shedding intensities - vertical
lines correspond to minima of $a_{2}^{2}$; right) correlation coefficient
between vortex shedding energy and quasi-steady deviation energy.
---
Figure 17:  Left) base suction coefficient $C_{B}=-C_{p}$ obtained by
integrating the pressure over the rear of the body; right) correlation between
the base suction coefficient and the POD amplitudes corresponding to the
steady deviation ($a_{2}$, top), wake pumping mode ($a_{7}$, middle) and
vortex shedding modes ($r_{Kh,v}^{2}$, bottom).

##  VI  Low-dimensional model

We now examine whether it is possible to model the behavior of the largest
scales using a POD-based model, even if the scales considered represent only a
fraction of the total fluctuating kinetic energy. Following the approach
described in [19], [21], we build a low-dimensional model to reproduce the
dynamics. We use a Galerkin approach to project the Navier-Stokes equations
onto the basis of spatial modes for a selected truncation, and obtain a set of
ordinary differential equations for the normalized amplitudes $a_{n}(t)$. The
model is of the form

$\dot{a}_{n}=L_{nm}a_{m}+Q_{nmp}a_{m}a_{p}+T_{n}.$ (6)

where

  * •

the linear terms contain the viscous dissipation

$L_{nm}=\int\nu\Delta\underline{\phi}_{m}.\underline{\phi}_{n}d\underline{x}.$
(7)

For the ten-mode truncation they form a diagonal matrix $L\sim-0.05I$ as
indicated in table 2.

  * •

$T_{n}$ is a closure term representing the contribution of the unresolved
stresses (associated with the modes excluded from the truncation) to the
evolution of the amplitude $a_{n}$.

  * •

the quadratic terms $Q_{nmp}$ can be written in symmetric form as

$Q_{nmp}=\frac{\sqrt{\lambda_{m}\lambda_{p}}}{\lambda_{n}}\frac{1}{2}(2-\delta_{mp})\int(\underline{\phi}_{p}.\nabla\underline{\phi}_{m})+\underline{\phi}_{m}.\nabla\underline{\phi}_{p}).\underline{\phi}_{n}d\underline{x}.$
(8)

For the evolution equation of the amplitude $a_{n}$, $1<n\leq 10$, shown in
equation ( 6), the interaction coefficient of $a_{n}$ with the mean mode,
$Q_{nn1}$, is essentially independent of $n$ and its value is about 0.2. A
physical interpretation of this is that each of the modes interacts directly
and equally with the mean mode, in particular its strong shear layers.

Generally speaking, the magnitudes of the quadratic coefficients $Q_{nmp}$
provide insight into the interactions between the different modes. Table 3
contains the values of the quadratic coefficients $Q_{nm1}$ for $3\leq n,m\leq
10$. The dominant values which will be kept for the model are indicated in
bold.

###  VI.1 Simplified model

We first consider a simplified version of the model by making the following
assumptions:

  * •

we assume that the first two modes $a_{1}$ and $a_{2}$ are constant.

  * •

we neglect quadratic terms of small magnitude.

  * •

we model the energy transfer to the unresolved modes by assuming that their
effect is to compensate for the production term i.e the interaction with the
mean shear (mode 1). This means that

$T_{n}=(-L_{n}-Q_{1nn}a_{1})a_{n}.$ (9)

This leads to the following form for the model, which we will refer to as
$M_{0}$:

$\displaystyle\dot{a}_{1}$ $\displaystyle=$ $\displaystyle 0$ (10)
$\displaystyle\dot{a}_{2}$ $\displaystyle=$ $\displaystyle 0$ (11)
$\displaystyle\dot{a}_{3}$ $\displaystyle=$ $\displaystyle 0.91a_{4}$ (12)
$\displaystyle\dot{a}_{4}$ $\displaystyle=$ $\displaystyle-0.94a_{3}$ (13)
$\displaystyle\dot{a}_{5}$ $\displaystyle=$
$\displaystyle-1.09a_{6}+0.49a_{9}$ (14) $\displaystyle\dot{a}_{6}$
$\displaystyle=$ $\displaystyle 1.22a_{5}$ (15) $\displaystyle\dot{a}_{7}$
$\displaystyle=$ $\displaystyle-0.42a_{8}$ (16) $\displaystyle\dot{a}_{8}$
$\displaystyle=$ $\displaystyle-0.55a_{5}+0.44a_{7}$ (17)
$\displaystyle\dot{a}_{9}$ $\displaystyle=$
$\displaystyle-0.6a_{4}-0.86a_{10}$ (18) $\displaystyle\dot{a}_{10}$
$\displaystyle=$ $\displaystyle-0.88a_{3}+0.98a_{9}$ (19)

We can see that there are only interactions between modes with the same
parity: (3,4,9,10) on the one hand and (5,6,7,8) on the other hand. The
interactions between the normalized amplitudes $(a_{2i-1},a_{2i}),2\leq i\leq
5$ are of the form
$a_{2i-1}=-\omega_{i}(\frac{\lambda_{2i}}{\lambda_{2i-1}})^{1/2}a_{2i}$,
$a_{2i}=\omega_{i}(\frac{\lambda_{2i-1}}{\lambda_{2i}})^{1/2}a_{2i-1}$, which
correspond to  propagative oscillatory solutions for the amplitudes
$\sqrt{\lambda_{n}}a_{n}$, in agreement with the convective dynamics expected
for the corresponding  Kármán modes.

The averaged frequencies $|\omega_{i}|=2\pi f_{i}$ identified for the pairs
(3,4), (5,6), (7,8) and (9,10) are 0.92, 1.15, 0.43, 0.9, which correspond to
frequencies (or Strouhal numbers) of 0.15, 0.18, 0.07 and 0.14. This is in
good agreement with the main frequencies identified in the simulation. We
emphasize that this prediction of the relevant time scales is based
exclusively on the spatial structure of the modes, which are extracted from a
set of samples arbitrarily separated in time. Since the computation is based
on the derivatives of the spatial modes, some uncertainty exists in the
determination of the time scales.

The model was integrated from a random initial condition, and the amplitudes
$a_{i},3\leq i\leq 10$ are represented in spectral space in figure 18. As
expected, the frequencies of the model coefficients agree well with the
dominant frequencies identified in the previous section for the amplitudes of
the modes in the simulation. The amplitudes of the modes are also close to
their expected values, as shown in table 2. These results indicate that the
main temporal dynamics of the flow can be recovered from the quadratic
interactions between spatial POD modes, even if the snapshots are obtained
with large separations, which is evidence of the predictive abilities of the
POD-based model.

Figure 18: Power spectral density of the 3-D POD mode amplitudes 3 to 10
$|\hat{a}_{n}^{3D,M_{0}}|^{2}$ in the simplified model \- the red and black
lines respectively correspond to the two frequencies 0.08 and 0.2.

###  VI.2 Switches

We now examine if and how switches can be reproduced by a more complex version
of the model, which we derive by relaxing the assumptions of the simplified
model as follows:

  * •

the second mode is allowed to vary.

  * •

feedback is provided between the unresolved terms and the modes of the
truncation.

As shown in [21], we assume that the rate of energy transferred to the small
scales depends on the energy available in the largest scales. If more energy
is available in the large scales, then more energy is extracted by the small
scales and conversely. This leads to us to introduce a time-varying viscosity
term, so that the effect of the unresolved terms is modelled as

$T_{n}=A_{n}a_{n}+\epsilon_{n}$ (20)

where

  1. 1.

the time-averaged value $<A_{n}>$ satisfies

$<A_{n}>=-L_{n}-Q_{nn1},$ (21)

  2. 2.

$A_{n}$ contains a linear part and a quadratic part

$A_{n}=<A_{n}>-\alpha_{n}\sum_{p=2}^{10}\lambda_{p}+\alpha_{n}\sum_{p=2}^{10}\lambda_{p}|a_{p}|^{2},$
(22)

As has been shown in [21], we have

$T_{n}=(<A_{n}>+\alpha_{n}\sum_{p\geq 1}^{N}\lambda_{p}-\alpha_{n}\sum_{p\geq
1}^{N}|a_{p}|^{2})a_{n}+\epsilon_{n}$ (23)

where

$\alpha_{n}=-\frac{<A_{n}>}{2\sum_{p\geq 2}^{N}\lambda_{p}}.$

The value of $\alpha_{i}$ was evaluated to be around 0.5. Since the model was
not found to be largely insensitive to the exact values of $\alpha_{i}$, in
what follows a constant value of $\alpha_{i}=\alpha=0.5$ was used for all
modes for the sake of simplicity. Examination of equation (20) and of the POD
eigenvalues shows that about 50% of the turbulent viscosity is dependent on
the amplitude of mode $a_{2}$. In that sense the structure of the model
displays similarities with both [22] and [3]’s model, which includes a cubic
term in $a_{2}$. However it is derived from entirely different physical
arguments.

  3. 3.

$\epsilon_{n}$ represents Gaussian noise. The amplitude of the noise used to
integrate the model was determined using $\epsilon_{n}\approx|A_{n}|$. We used
$\epsilon_{2}=0.09$ and $\epsilon_{i}=0.04$ for $3\leq i\leq 10$. [A Gaussian
variable was generated every $2.5$ convective time units and linearly
interpolated was used for the integration].

The modified model reads:

$\displaystyle\dot{a}_{1}$ $\displaystyle=$ $\displaystyle 0$ (24)
$\displaystyle\dot{a}_{2}$ $\displaystyle=$
$\displaystyle(L_{2}-\alpha\sum_{2}^{10}a_{p}^{2}\lambda_{p})a_{2}+\epsilon_{2}$
(25) $\displaystyle\dot{a}_{3}$ $\displaystyle=$
$\displaystyle(L_{3}-\alpha\sum_{2}^{10}a_{p}^{2}\lambda_{p})a_{2}+0.71a_{4}+\epsilon_{3}$
(26) $\displaystyle\dot{a}_{4}$ $\displaystyle=$
$\displaystyle(L_{4}-\alpha\sum_{2}^{10}a_{p}^{2}\lambda_{p})a_{2}-0.74a_{3}+\epsilon_{4}$
(27) $\displaystyle\dot{a}_{5}$ $\displaystyle=$
$\displaystyle(L_{5}-\alpha\sum_{2}^{10}a_{p}^{2}\lambda_{p})a_{2}-1.09a_{6}+\epsilon_{5}$
(28) $\displaystyle\dot{a}_{6}$ $\displaystyle=$
$\displaystyle(L_{6}-\alpha\sum_{2}^{10}a_{p}^{2}\lambda_{p})a_{2}-1.2a_{5}+\epsilon_{6}$
(29) $\displaystyle\dot{a}_{7}$ $\displaystyle=$
$\displaystyle(L_{7}-\alpha\sum_{2}^{10}a_{p}^{2}\lambda_{p})a_{2}-0.42a_{8}+\epsilon_{7}$
(30) $\displaystyle\dot{a}_{8}$ $\displaystyle=$
$\displaystyle(L_{8}-\alpha\sum_{2}^{10}a_{p}^{2}\lambda_{p})a_{2}-0.55a_{5}+0.44a_{7}+\epsilon_{8}$
(31) $\displaystyle\dot{a}_{9}$ $\displaystyle=$
$\displaystyle(L_{9}-\alpha\sum_{2}^{10}a_{p}^{2}\lambda_{p})a_{2}-0.6a_{4}-0.86a_{10}+\epsilon_{9}$
(32) $\displaystyle\dot{a}_{10}$ $\displaystyle=$
$\displaystyle(L_{10}-\alpha\sum_{2}^{10}a_{p}^{2}\lambda_{p})a_{2}-0.88a_{3}+0.98a_{9}+\epsilon_{10}.$
(33)

The quadratic terms in the expression for $a_{2}$ are not included in the
model as their magnitude was small (we checked that including these terms in
the equations did not change the dynamics reported below).

The effect of the feedback term on the dynamics of the model is that if there
is less energy in mode $a_{2}$, the higher-order modes will extract less
energy, which will allow mode $a_{2}$ to grow. Conversely, if mode $a_{2}$
becomes too large, the energy transfer to the higher-order modes will be
increased, which will in turn affect the energy of mode $a_{2}$.

Figures 19 and 20 shows results of the model integration for a noise amplitude
of about 0.15. Figure 19 (left) shows the 3-D coefficient $a_{2}^{3d}$
predicted with the model M along with the 2-D coefficient $a_{2}^{2d}$ in the
experiment, which appears a relevant comparison since, as shown in the
previous section, there is a reasonably good correlation between $a_{2}$ in
3-D and 2-D (0.6). The model displays time scales of $O(1000)$ between
switches, in agreement with experimental observations. Figure 19 (right) shows
that the histogram of the amplitude is similar to that observed in the
experiment in figure 6. This shows that the model is able to reproduce
deviations in a way that is consistent with experiments. As figure 20
indicates, modes 3 to 10 are characterized by relatively fast oscillating time
scales and slower amplitude variations. The frequencies of the amplitudes are
shown in figure 20 (right) and compare relatively well with those measured in
the simulation, given the crudeness of the truncation. Table 2 (last line)
shows that the magnitude of the normalized POD amplitudes is relatively well
estimated by the model with values of about 0.2 to 2 for the last modes of the
truncation (we emphasize that the model contains only 20% of the total
fluctuating energy). The main dynamics of the largest scales are therefore
captured by the model.

n | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
$L_{i}$ |  | -0.05 | -0.05 | -0.05 | -0.05 | -0.05 | -0.05 | -0.05 | -0.05 | -0.05
$<a_{n}^{2}>_{M_{0}}$ | 1 | 1 | 0.93 | 0.93 | 0.86 | 0.86 | 1.45 | 1.51 | 1.05 | 1.05
$<a_{n}^{2}>_{M}$ | 1 | 1.05 | 0.29 | 0.29 | 0.18 | 0.18 | 1.01 | 1.03 | 2.05 | 2.33
Table 2:  Model linear coefficients and predicted energy $Q_{1jm}$ | m=3 | m=4 | m=5 | m=6 | m=7 | m=8 | m=9 | m=10
---|---|---|---|---|---|---|---|---
j=3 | 0.21 | 0.91 |  |  |  | -0.22 |  | 0.37
j=4 | -0.94 | 0.23 |  |  |  | 0.0 | 0.29 | 0.17
j=5 |  |  | 0.22 | 1.09 | 0. | 0.34 |  |
j=6 |  |  | -1.22 | 0.20 | -0.12 | -0.10 |  |
j=7 |  |  |  | 0.17 | 0.16 | -0.42 |  |
j=8 |  |  | 0.55 |  | 0.4 | 0.17 |  |
j=9 | 0.32 |  | 0.59 |  |  |  | 0.22 | -0.86
j=10 | -0.88 | -0.18 |  |  |  |  | 0.98 | 0.24
Table 3: Quadratic interaction coefficients with the main mode $Q_{j1m}$; Only
coefficients larger than 0.1 are indicated; only coefficients larger than 0.4
(indicated in bold) are included in the model.

Figure 19:  Left) POD amplitude $a_{2}$; top: model M (3D); bottom: experiment
(2D); right) histogram of $a_{2}$

Figure 20:  Left) amplitudes of 3-D POD amplitudes of modes 3, 5 , 7 and 9 in
the model; right) power spectral density of POD amplitudes in the model
$|\hat{a}_{n}^{3d,M}|^{2}$ \- the red and black lines respectively correspond
to the two frequencies 0.08 and 0.2.

##  VII Conclusion

We have applied Proper Orthogonal Decomposition to the 3-D numerical
simulation of the flow behind an Ahmed body at $Re=10^{4}$. Reflection
symmetry was applied to the computed data set in order to compensate for the
relatively short time of the simulation, which precludes the observation of
switches in the wake deviation. As a consequence of the enforced statistical
symmetry, the flow can be decomposed into symmetric and antisymmetric
structures. The mean flow consists of a symmetric recirculation bubble and an
antisymmetric deviation, the effect of which is to gather flow streamlines
around one of the base diagonals. 2-D POD analysis was performed in the near-
wake in the simulation and compared with experimental results obtained for the
same geometry. Despite the discrepancy in Reynolds number between the
simulation and the experiment, an excellent agreement was observed for both
the spatial structure and temporal statistics of the POD modes.

2-D results were then confronted with a 3-D approach. The energetic importance
of the quasi-steady wake deviation was established. The evolution of this
global 3-D deviation mode was relatively well captured by 2-D measurements in
the near-wake. The next most energetic patterns are associated with vortex
shedding and wake pumping. Characteristic frequencies were identified for each
type of structure. Structures associated with wake pumping were characterized
by a low frequency of about 0.08. Both symmetric and antisymmetric structures
associated with vortex shedding in respectively the vertical and spanwise
direction were characterized by dominant frequencies of 0.19 and 0.23 in the
far wake. A strong negative correlation was noted between the intensity of the
vortex shedding modes and the magnitude of the quasi-steady deviation. In
addition, increases in the base drag were found to correspond to an increase
of the quasi-steady deviation, along with a shrinkage of the recirculation
zone associated with wake pumping.

Finally, POD-based low-dimensional models were derived for the largest scales
of the flow. The energy content of the modes was correctly captured by the
model. A simplified model was able to single out the main frequencies of the
POD amplitudes observed in the simulation from the spatial modes, regardless
of the time separation between the snapshots used to compute POD. This
predictive ability of the model is remarkable in view of the slow convergence
of the decomposition, which reflects the complexity of the flow. A more
elaborate version of the model was also considered. The approach is consistent
with Rigas et al.’s model [22] and in particular the structure of the model
would be the same if the POD truncation was limited to two modes. By
accounting for the effect of the unresolved modes with a feedback term, the
POD-based model was able to reproduce the characteristics of the switches in
the wake deviation. The success of the model supports the idea that wake
switching is triggered by higher-order modes.

## References

  * [1] R.D. Brackston, J.M. Garci De La Cruz, A. Wynn, G. Rigas, and J.F. Morrison.  Stochastic modelling and feedback control of bistability in a turbulent bluff bodywake.  J. Fluid Mech., 802:726–749, 2016.
  * [2] O. Cadot.  Stochastic fluid structure interaction of three-dimensional plates facing a uniform flow.  J. Fluid Mech., 794:726–749, 2016.
  * [3] O. Cadot, A. Evrard, and L. Pastur.  Imperfect supercritical bifurcation in a three-dimensional turbulent wake.  Phys. Review E., 91(6), 2015.
  * [4] L. Dalla Longa, O. Evstafyeva, and A. S. Morgans.  Simulations of the bi-modal wake past three-dimensional blunt bluff bodies.  Journal of Fluid Mechanics, 866:791–809, 2019.
  * [5] A. Evrard, O. Cadot, V. Herbert, D. Ricot, R. Vigneron, and J. Delery.  Fluid force and symmetry breaking modes of a 3d bluff body with a base cavity.  Journal of Fluids and Structures, 61:99–114, 2016.
  * [6] O. Evsafyeva, A. Morgans, and L. Dalla Longa.  Simulation and feedback control of the ahmed body flow exhibiting symmetry breaking behaviour.  J. Fluid Mech., 817, 2017.
  * [7] D. Fabre, F. Auguste, and J. Magnaudet.  Bifurcations and symmetry breaking in the wake of axisymmetric bodies.  Phys. Fluids, 20(0517), 2017.
  * [8] M. Grandemange, O. Cadot, and M. Gohlke.  Reflectional symmetry breaking of the separated flow over three-dimensional bodies.  Physical Review E, 86(3):035302, 2012.
  * [9] M. Grandemange, M. Gohlke, and O. Cadot.  Turbulent wake past a three-dimensional blunt body. part 1 - experimental sensitivity analysis.  J. Fluid Mech., 752:439–461, 2013.
  * [10] M. Grandemange, M. Gohlke, and O. Cadot.  Turbulent wake past a three-dimensional blunt body. part 1 - global modes and bi-stability.  J. Fluid Mech., 722:1–26, 2013.
  * [11] P. Holmes, J.L. Lumley, and Gal Berkooz.  Turbulence, Coherent Structures, Dynamical Systems and Symmetry.  Cambridge University Press, 1996.
  * [12] M. Kiya and Y. Abe.  Turbulent elliptic wakes.  J. Fluids Struct., 13:1041–1067, 1999.
  * [13] R. Li, D. Barros, J. Borée, O. Cadot, B.R. Noack, and L. Cordier.  Feddback control of bimodal wake dynamics.  Exp Fluid, 51(4):158, 2016.
  * [14] J.M. Lucas, O. Cadot, V. Herbert, S. Parpais, and J. Délery.  A numerical investigation of the asymmetric wake mode of a squareback ahmed body - effect of a base cavity.  J. Fluid Mech., 831:675–697, 2017.
  * [15] J.L. Lumley.  The structure of inhomogeneous turbulent flows.  In A.M Iaglom and V.I Tatarski, editors, Atmospheric Turbulence and Radio Wave Propagation, pages 221–227. Nauka, Moscow, 1967.
  * [16] R. Pasquetti and N. Peres.  Simulation and feedback control of the ahmed body flow exhibiting symmetry breaking behaviour.  Computers and Fluids, 114:203–217, 2015.
  * [17] G. Pavia, M. Passmore, and C. Sardu.  Evolution of the bi-stable wake of a square-back automotive shape.  Exp. in Fluids, 59:2742, 2018.
  * [18] A.K. Perry, G. Pavia, and M. Passmore.  Influence of short rear end tapers on the wake of a simplified square-back vehicle: wake topology and rear drag.  Experiments in Fluids, 57(11), 2016.
  * [19] B. Podvin.  A pod-based model for the wall layer of a turbulent channel flow.  Phys. Fluids, 21(1):015111, 2009.
  * [20] B. Podvin and Y. Fraigneau.  A few thoughts on proper orthogonal decomposition in turbulence.  Physics of FLuids, 29:020709, 2017.
  * [21] B. Podvin and A. Sergent.  Precursor for wind reversal in a square rayleigh-bénard cell.  Phys. Rev. E, 95:013112, 2017.
  * [22] G. Rigas, A.S Morgans, R.D. Brackston, and J.F. Morrison.  Diffusive dynamics and stochastic models of turbulent axisymmetric wakes.  J. Fluid Mech., 778(R2), 2015.
  * [23] G. Rigas, A.S Morgans, and J.F. Morrison.  Stability and coherent structures in the wake of axisymmetric bluff bodies.  Fluid Mechanics and its Applications, 107:143–148, 2015.
  * [24] G. Rigas, A.R. Oxlade, A.S Morgans, and J.F. Morrison.  Low-dimensional dynamics of a turbulent axisymmetric wake.  J. Fluid Mech., 755:159, 2014.
  * [25] P.J. Schmid.  Dynamic mode decomposition of numerical and experimental data.  J. Fluid Mech., 656:5–28, 2010.
  * [26] L. Sirovich.  Turbulence and the dynamics of coherent structures part i: Coherent structures.  Quart. Appl. Math., 45(3):561–571, 1987.
  * [27] E. Varon, Y. Eulalie, S. Edwige, P. Gilotte, and J.L. Aider.  Chaotic dynamics of large-scale structures in a turbulent wake.  Phys. Review Fluids ., 2(034604), 2017.
  * [28] R. Volpe, P. Devinant, and A. Kourta.  Experimental characterization of the unsteady natural wake of the full-scale square back ahmed body: flow bi-stability and spectral analysis.  Exp. in Fluids, 56(5):1–22, 2015.
  * [29] E. Wassen, S. Eichinger, and F. Thiele.  Simulation of active drag reduction for a square-back vehicle.  Notes on Numerical Fluid Mechanics and Multidisciplinary Design, 108:241–255, 2014.