# Non-homogeneous anisotropic bulk viscosity for
acoustic wave attenuation in weakly compressible methods

Dheeraj Raghunathan  School of Mechanical Sciences
Indian Institute of Technology Goa
Farmagudi-403401, India  Y. Sudhakar

###### Abstract

A major limitation of the weakly compressible approaches to simulate
incompressible flows is the appearance of artificial acoustic waves that
introduce a large mass conservation error and lead to spurious oscillations in
the force coefficients. In this work, we propose a non-homogeneous anisotropic
bulk viscosity term to effectively damp the acoustic waves. By implementing
this term in a computational framework based on the recently proposed general
pressure equation, we demonstrate that the non-homogeneous and anisotropic
nature of the term makes it significantly more effective than the isotropic
homogeneous version widely used in the literature. Moreover, it is
computationally more efficient than the pressure (or mass) diffusion term,
which is an alternative mechanism used to suppress acoustic waves. We simulate
a range of benchmark problems to comprehensively investigate the performance
of the bulk viscosity on the effective suppression of acoustic waves, mass
conservation error, order of convergence of the solver, and computational
efficiency. The proposed form of the bulk viscosity enables fairly accurate
modelling of the initial transients of unsteady simulations, which is highly
challenging for weakly compressible approaches, and to the best of our
knowledge, existing approaches can’t provide an accurate prediction of such
transients.

###### keywords:

Bulk viscosity , Artificial compressibility methods , Weakly compressible
methods , General pressure equation , Incompressible flow , Acoustic waves

††journal: *************††volume: ≈††volume: ≈††volume: missing††volume:
≈††volume: missing††volume: ≈††volume: missing††volume: missing††volume:
missing††volume: missing††volume: missing††volume: ,

##  1 Introduction

Incompressible fluid flows are mathematically described by the Navier-Stokes
equations given by the following non-dimensional form,

$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=0,$ (1a)
$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}$ $\displaystyle=-\nabla
p+\frac{1}{Re}\nabla^{2}\mathbf{u},$ (1b)

where $\mathbf{u}$ and $p$ denote velocity vector and pressure, respectively;
$Re$ is the Reynolds number. The major challenge to solve them numerically is
the lack of an evolution equation for pressure. The most widely used family of
methods [1, 2, 3], by manipulating the above equations, obtain an equation for
pressure or pressure correction ($\phi$) of the following form

$\nabla^{2}\phi=\frac{1}{\Delta t}\nabla\cdot\mathbf{u}^{*},$ (2)

where $\mathbf{u}^{*}$ represents an intermediate velocity vector. It is
computed from the momentum equations but does not satisfy the continuity
equation. The numerical solution of such an elliptic equation is highly
expensive, requires extensive memory, and offers challenges to achieve
scalability on parallel computing architectures [4].

To avoid solving the expensive Poisson equation, Chorin [5] introduced the
concept of pseudo-compressibility and proposed the following pressure
evolution equation that replaces the continuity equation

$\frac{\partial p}{\partial\tau}+\frac{1}{Ma^{2}}\nabla\cdot\mathbf{u}=0.$ (3)

Here, $Ma$ is the Mach number, which is a numerical parameter that quantifies
the pseudo-compressibility, and $\tau$ represents the pseudo-time. The above
equation, together with the momentum equations, can be used to simulate
steady-state incompressible flows [6, 7]. This technique is called the
artificial compressibility method (ACM). The utility of ACM is further
extended to unsteady flows by adopting a dual time-stepping [8] approach, in
which the momentum equations are also appended with a pseudo-time derivative
term, as shown below

$\frac{\partial\mathbf{u}}{\partial\tau}+\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}=-\nabla
p+\frac{1}{Re}\nabla^{2}\mathbf{u}.$ (4)

Such a modification insists that for every physical time step of the solution
process, an extra set of inner iterations needs to be performed until the
pseudo-time derivatives diminish to a predefined tolerance. The dual time-
stepping approach can be easily ported to a parallel computing framework and
is extensively employed to simulate a wide variety of flow problems [9, 10,
11, 12, 13, 14]. However, the presence of the inner iterations makes them
computationally expensive.

The computational cost of the aforementioned approaches to simulate unsteady
incompressible flows is dictated by the Poisson equation or pseudo-time
iterations. Both these operations are highly time-consuming, and their
necessity stems from the unavailability of a separate equation for pressure.
In order to circumvent these expensive steps, there has been a growing
interest in recent years to derive, from conservation laws, a time-evolution
equation for pressure to model incompressible flows. Two notable works that
received considerable attention are the entropically damped artificial
compressiblity method (EDAC) [15] and the general pressure equation (GPE) [16,
17]. While EDAC uses the entropy balance equation to formulate a thermodynamic
constraint equation, GPE is obtained by modifying the energy conservation
equation. We underline that the computational frameworks developed based on
these equations neither require solving a pressure Poisson equation nor
involve dual time schemes. Since they have the potential to be efficiently
parallelised on the modern CPU and GPU supercomputing architectures, these
methods are in a state of rapid development. Several algorithms, utilizing GPE
and EDAC, have been proposed in recent years to handle laminar [15, 17],
turbulent [18, 19, 20, 21], two-phase [22, 23, 24], buoyancy-driven [25] flows
and moving boundary simulations [26]. Another relevant formulation is the
kinetically reduced local Navier-Stokes [27] equation, which does not directly
solve for pressure but uses grand potential from which the pressure field is
obtained.

Similar to the ACM, numerical techniques developed based on these formulations
belong to _weakly compressible_ frameworks. A major limitation of such methods
is the appearance of artificial acoustic (pressure) waves in the flow field
that would have been absent in incompressible flows. Lu and Adams [28]
highlight that there are two mechanisms used in weakly compressible methods to
handle acoustic waves: diffusion term in the pressure (or density) equation
and bulk viscosity. Toutant [17] conducted a preliminary study on the effect
of increasing the amount of pressure diffusion in the GPE to control the
acoustic waves. However, we will show that this provides only a limited
benefit, which comes at the expense of an additional time-step restriction for
stable computations.

The presence of acoustic waves leads to mass conservation errors, extended
time to reach a periodic vortex-shedding phase of bluff body flows, highly
oscillatory flow-induced forces in truly unsteady flow problems, etc. These
limitations need to be addressed in order to make such methods useful in
practical applications. We will demonstrate in this paper that the non-
homogeneous anisotropic form of bulk viscosity is very effective in addressing
all these limitations.

###  1.1 Bulk viscosity in weakly compressible methods

The idea of using bulk viscosity in weakly compressible formulations can be
traced back to the work of Ramshaw and Mousseau [29]. They added the bulk
viscosity to improve the convergence rate of the ACM simulations towards the
steady state. McHugh and Ramshaw [30] employed bulk viscosity in a dual-time
scheme solver with fully implicit discretization. They reported that this term
is effective at large time steps but produced no effect on simulations with
small time steps. Although the solver is unsteady, they reported results only
for steady-state problems. Bulk viscosity has also been used to accelerate the
convergence of low Mach number flows to steady state [31, 32]. A recent study
[33] investigated the influence of bulk viscosity on stability and
computational efforts in an unsteady simulation using the ACM. For the doubly
periodic shear layer problem chosen, they reported that the bulk viscosity
suppressed spurious vortices and enabled stable simulations even on coarse
grids. Huang [23], in his work on the development of a method to handle two-
phase flows using GPE, reported that bulk viscosity is mandatory for stable
computation of three-dimensional flows.

The inclusion of bulk viscosity term is not limited to methods that solve
Navier-Stokes equations. Its use has been reported in other weakly
compressible approaches like Lattice Boltzmann methods (LBM) and Smoother
Particle Hydrodynamics (SPH) techniques.

Dellar [34] reported that the bulk viscosity, when its value is set to be
larger than that of the shear viscosity, damped artifacts in unresolved
simulations of incompressible flows using LBM. The term acts only in regions
of velocity divergence errors, and leaves the other parts affected. Asinari &
Karlin [35] proposed an LBM which contains a free parameter to control the
second viscosity to enhance the numerical stability in the incompressible
limit. They underlined that even though it is not explicitly expressed in the
literature, the stability of LBMs in the incompressible can be attributed to
the bulk viscosity term.

In the weakly compressible SPH framework, Avalos et al. [36] introduced a bulk
viscosity term that can be modelled independently of the shear viscosity. The
proposed term conserves both linear and angular momentum, which are desirable
for simulating flows with free surfaces. In a most recent work, Sun et al.
[37] proposed bulk viscosity as an acoustic damper term in SPH and
demonstrated a drastic reduction in pressure oscillations. They rigorously
tested the effectiveness of this term on test cases involving violent fluid-
fluid and fluid-solid interaction problems.

###  1.2 Summary of existing studies and novelty of the present work

This section provides a brief summary of existing works on weakly compressible
approaches and bulk viscosity. We underline the limitations of these studies
and elaborate on the contributions of the present work.

  1. 1.

By analyzing several weakly compressible approaches, Lu and Adams [28] note
that there are two general mechanisms used for stable computations: pressure
diffusion (also referred to as mass diffusion) and bulk viscosity. Despite
their widespread usage, the relative effectiveness of these two mechanisms
remains unclear. Lu and Adams [28] remarked that both of these are effective.
In this paper, we shed light on, for the first time, the effectiveness of
these two mechanisms. Our results reveal that the bulk viscosity has a
stronger impact in suppressing acoustic waves.

  2. 2.

The existing works on introducing the bulk viscosity either addressed its
effect on the convergence rate of steady state solvers [29], or on the
stability [33]. Sun et al. [37] is the only study that extensively
investigated the artificial damping effect but in the SPH framework. In this
paper, we present a _comprehensive_ study of bulk viscosity on various aspects
of an unsteady solver: damping of acoustic waves, order of convergence, and
mass conservation error.

  3. 3.

The key contribution of the present work is the use of non-homogeneous
anisotropic bulk viscosity. These features enable the rapid annihilation of
artificial acoustic waves. We will show that this leads to accurate
computation of force coefficients on an impulsively started plate, and a
drastic reduction in computation cost to model von-Karman vortex shedding
behind a bluffbody. The above advantages would not have been realized with the
isotropic and uniform bulk viscosity used throughout the literature. An
earlier technical report [38] briefly discusses the use of anisotropic but
uniform bulk viscosity tensor; however, no numerical tests were presented. To
the best of our knowledge, there is no other existing weakly compressible
approach that utilises this idea.

The weakly compressible framework we use in this work is built on GPE. This is
because Toutant [17] showed that the extra convective term in the EDAC
equation offers a negligible influence on the evolution of the flow field.

###  1.3 Structure of the paper

The rest of the paper is arranged as follows. The governing equations and
their spatial as well as temporal discretisation details are presented in § 2,
with an emphasis on the anisotropic nature and the spatial variation of the
bulk viscosity term. A series of numerical experiments conducted to
demonstrate the influence of the bulk viscosity term on different aspects of
the numerical simulation are summarised in § 3. Conclusions drawn from these
studies are described in § 4. Some additional discussions are also included in
A and B to supplement our findings.

##  2 Governing equations and numerical methods

In this section, the governing equations are presented, along with the
numerical discretisation details. We use the following two approaches to
simulate unsteady incompressible flow problems numerically.

###  2.1 GPE without artificial bulk viscosity (denoted as ‘GPE solver’)

GPE is derived from the compressible flow and energy conservation equations in
the low Mach number limit [16], and for an isothermal incompressible flow, GPE
is written in non-dimensional form as

$\frac{\partial p}{\partial
t}+\frac{1}{Ma^{2}}\nabla\cdot\mathbf{u}=\frac{1}{RePr}\nabla^{2}p$ (5)

where, $\mathbf{u}=[\,u\quad v\,]^{\top}$ is the velocity vector ($u$ and $v$
being its components in a 2D Cartesian coordinate system) and $p$ denotes
pressure. The parameters $Ma$ and $Pr$ are the Mach number and Prandtl number,
respectively; $Re$ represents the Reynolds number of the flow. Here, $Ma$ and
$Pr$ are the numerical parameters, and $Re$ is the only physical parameter.

Together with GPE, the momentum equations (equation (1b)) govern the evolution
of the flow field in the incompressible limit. This approach has been used to
simulate various laminar and turbulent incompressible flow problems [17, 21,
20, 39, 24].

###  2.2 GPE with artificial bulk viscosity (denoted as ‘GPE+BV solver’)

In this approach, we solve the same GPE equation (equation (5)) to determine
pressure. In contrast to the GPE solver, the momentum equation is appended
with an artificial bulk viscosity tensor ($\boldsymbol{\mathcal{B}}$). The
modified momentum equation is presented below

$\frac{\partial\mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}=-\nabla
p+\frac{1}{Re}\nabla^{2}\mathbf{u}+\nabla\cdot\left(\boldsymbol{\mathcal{B}}\nabla\cdot\mathbf{u}\right)$
(6)

where the artificial bulk viscosity term can be expanded as
$\nabla\cdot\left(\boldsymbol{\mathcal{B}}\nabla\cdot\mathbf{u}\right)=\boldsymbol{\mathcal{B}}\nabla\left(\nabla\cdot\mathbf{u}\right)+(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$.
The second term
$(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$
arises due to the non-homogeneous nature of the bulk viscosity as proposed in
this work. Note that the above equation is in non-dimensional form, and
$\boldsymbol{\mathcal{B}}=(\boldsymbol{\mathcal{B}}^{*}/\rho UL)$ represents
an analogous Reynolds number based on the dimensional
$\boldsymbol{\mathcal{B}}^{*}$. In the present work, we propose the following
form of the anisotropic bulk viscosity tensor.

$\boldsymbol{\mathcal{B}}=\begin{bmatrix}\mathcal{B}^{\mathcal{X}}&0\\\
0&\mathcal{B}^{\mathcal{Y}}\end{bmatrix}$ (7)

where $\mathcal{B}^{\mathcal{X}}$ and $\mathcal{B}^{\mathcal{Y}}$ are the
components of the tensor in their usual notation. We underline that these
components are grid-size dependent and, hence, non-homogeneous in the
computational domain. We propose the following functional form of these
components

$\displaystyle\mathcal{B}^{\mathcal{X}}$ $\displaystyle=\lambda\Delta x,$ (8)
$\displaystyle\mathcal{B}^{\mathcal{Y}}$ $\displaystyle=\lambda\Delta y,$ (9)

where $\lambda>0$ is a constant, $\Delta x$ and $\Delta y$ represent the grid
size in $x$ and $y$ directions respectively. In a non-uniform Cartesian mesh,
the components of bulk viscosity, for the chosen form, are functions of the
respective coordinates only, i.e., $\mathcal{B}^{\mathcal{X}}=f(x)$ and
$\mathcal{B}^{\mathcal{Y}}=f(y)$. In a uniform mesh with $\Delta x=\Delta y$,
the tensor will be isotropic and uniform, represented by
$\boldsymbol{\mathcal{B}}=\mathcal{B}\bf{I}$ where $\mathcal{B}$ is a scalar
and $\bf{I}$ is the identity tensor. We use the boldface symbol
$\boldsymbol{\mathcal{B}}$ to denote bulk viscosity tensor in a non-uniform
mesh and $\mathcal{B}$ to denote its scalar value on a uniform mesh. We
emphasize while the bulk viscosity is a physical parameter for compressible
flows, it is a free parameter in our work, and hence, its values can be set
independent of the shear viscosity [36].

The choice of the functional form for $\boldsymbol{\mathcal{B}}$ is driven by
the following argument. To completely eliminate artificial waves, the time
scale associated with the bulk viscosity should be much smaller than the
acoustic time scale. For a 1D flow, in non-dimensional settings, this implies
(refer to equations (16) and (20))

$\frac{\Delta x^{2}}{\mathcal{B}}\ll\frac{Ma\Delta x}{U}.$ (10)

However, ensuring the above condition requires a reduction in the time step of
the simulation, which leads to a significant increase in computational cost.
So, our objective is to introduce maximum diffusion of acoustic waves without
introducing additional restrictions on the time step. This leads to

$\gamma\frac{\Delta x^{2}}{\mathcal{B}}=\frac{Ma\Delta x}{U},$ (11)

where $\gamma$ is a positive constant, such that $0\leq\gamma\leq 1$,
depending on the problem. From the above, we get

$\mathcal{B}=\frac{\gamma U\Delta x}{Ma},$ (12)

and by defining $\lambda=\gamma U/Ma$, we obtain the expression defined
earlier, $\mathcal{B}=\lambda\Delta x$.

Figure 1: The computational stencil, indicating the arrangement of variables
on a staggered grid.

###  2.3 Spatial Discretization

The governing equations are recast into integral form, and a finite volume
approach is adopted to discretize the governing equations. We use the
Cartesian coordinate system in which the axes are denoted by $x$ and $y$.
Figure 1 shows the computational stencil used for discretisation. In order to
avoid odd-even decoupling, the variables are staggered. Pressure is stored at
the cell centre, and the velocity components are on the cell face centres. All
spatial derivatives appearing in the governing equations are approximated
using the central difference scheme similar to [17, 20]. The convective terms
are treated in their conservative form. In the following, we present only the
discretization details of the bulk viscosity term.

  1. 1.

$x$-component of
$\boldsymbol{\mathcal{B}}\nabla\left(\nabla\cdot\mathbf{u}\right)$

$\displaystyle\int_{B}^{\mathcal{X}}\frac{\partial}{\partial
x}(\nabla\cdot\mathbf{u})d$
$\displaystyle\frac{\mathcal{B}^{\mathcal{X}}_{i,j}}{(0.5\Delta
x_{i}+0.5\Delta x_{i-1})}$
$\displaystyle\biggl{[}\left(\frac{u_{i+1,j}-u_{i,j}}{\Delta
x_{i}}+\frac{v_{i,j+1}-v_{i,j}}{\Delta y_{j}}\right)-$
$\displaystyle\left(\frac{u_{i,j}-u_{i-1,j}}{\Delta
x_{i-1}}+\frac{v_{i-1,j+1}-v_{i-1,j}}{\Delta
y_{j}}\right)\biggr{]}\Delta\mathcal{X}_{i,j}$

where, $\mathcal{B}^{\mathcal{X}}_{i,j}=\lambda(0.5\Delta x_{i}+0.5\Delta
x_{i-1})$ and $\Delta\mathcal{X}_{i,j}=(0.5\Delta x_{i}+0.5\Delta
x_{i-1})\Delta y_{j}$.

  2. 2.

$y$-component of
$\boldsymbol{\mathcal{B}}\nabla\left(\nabla\cdot\mathbf{u}\right)$

$\displaystyle\int_{B}^{\mathcal{Y}}\frac{\partial}{\partial
y}(\nabla\cdot\mathbf{u})d$
$\displaystyle\frac{\mathcal{B}^{\mathcal{Y}}_{i,j}}{(0.5\Delta
y_{j}+0.5\Delta y_{j-1})}$
$\displaystyle\biggl{[}\left(\frac{u_{i+1,j}-u_{i,j}}{\Delta
x_{i}}+\frac{v_{i,j+1}-v_{i,j}}{\Delta y_{j}}\right)-$
$\displaystyle\left(\frac{u_{i+1,j-1}-u_{i,j-1}}{\Delta
x_{i}}+\frac{v_{i,j}-v_{i,j-1}}{\Delta
y_{j-1}}\right)\biggr{]}\Delta\mathcal{Y}_{i,j}$

where, $\mathcal{B}^{\mathcal{Y}}_{i,j}=\lambda(0.5\Delta y_{j}+0.5\Delta
y_{j-1})$ and $\Delta\mathcal{Y}_{i,j}=\Delta x_{i}(0.5\Delta y_{j}+0.5\Delta
y_{j-1})$.

  3. 3.

$x$-component of
$(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$

$\displaystyle\int_{\partial\mathcal{B}^{\mathcal{X}}}{\partial
x}(\nabla\cdot\mathbf{u})d\biggl{[}$ $\displaystyle\left(\frac{\lambda\Delta
x_{i}-\lambda\Delta x_{i-1}}{0.5\Delta x_{i}+0.5\Delta
x_{i-1}}\right)\left(\nabla\cdot\mathbf{u}\right)^{\mathcal{X}}_{i,j}\biggr{]}\Delta\mathcal{X}_{i,j}$

where,

$\displaystyle\left(\nabla\cdot\mathbf{u}\right)^{\mathcal{X}}_{i,j}=$
$\displaystyle\left(\frac{u_{i+1,j}-u_{i-1,j}}{\Delta x_{i}+\Delta
x_{i-1}}\right)+$ $\displaystyle\frac{1}{\Delta y_{j}}$
$\displaystyle\left(\frac{v_{i,j+1}\Delta x_{i-1}+v_{i-1,j+1}\Delta
x_{i}}{\Delta x_{i}+\Delta x_{i-1}}-\frac{v_{i,j}\Delta
x_{i-1}+v_{i-1,j}\Delta x_{i}}{\Delta x_{i}+\Delta x_{i-1}}\right)$

  4. 4.

$y$-component of
$(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$

$\displaystyle\int_{\partial\mathcal{B}^{\mathcal{Y}}}{\partial
y}(\nabla\cdot\mathbf{u})d\biggl{[}$ $\displaystyle\left(\frac{\lambda\Delta
y_{j}-\lambda\Delta y_{j-1}}{0.5\Delta y_{j}+0.5\Delta
y_{j-1}}\right)\left(\nabla\cdot\mathbf{u}\right)^{\mathcal{Y}}_{i,j}\biggr{]}\Delta\mathcal{Y}_{i,j}$

where,

$\displaystyle\left(\nabla\cdot\mathbf{u}\right)^{\mathcal{Y}}_{i,j}=\frac{1}{\Delta
x_{i}}$ $\displaystyle\left(\frac{u_{i+1,j-1}\Delta y_{j}+u_{i+1,j}\Delta
y_{j-1}}{\Delta y_{j}+\Delta y_{j-1}}-\frac{u_{i,j-1}\Delta
y_{j}+u_{i,j}\Delta y_{j-1}}{\Delta y_{j}+\Delta y_{j-1}}\right)+$
$\displaystyle\left(\frac{v_{i,j+1}-v_{i,j-1}}{\Delta y_{j}+\Delta
y_{j-1}}\right)$

Here, $i$ and $j$ indicate the indices in $x$ and $y$ directions,
respectively, whereas $\Delta x_{i}$ and $\Delta y_{j}$ represent the grid
size of a particular control volume as shown in figure 1.

###  2.4 Temporal Discretization

We apply the three-stage Strong Stability Preserving Runge-Kutta [40] method
for the time integration. This procedure is briefly explained below for an
ordinary differential equation with the dependent variable $\Phi$. The
momentum equation and the GPE can be recast into the following form after the
spatial discretization

$\frac{d\Phi}{dt}=\mathit{\mathbb{L}}(\Phi),$ (13)

where $\mathit{\mathbb{L}}$ is the spatial discretisation operator. We follow
the steps below to update $\Phi$ from time level $n$ to $n+1$.

$\displaystyle\begin{split}\Phi^{(1)}=&\>\>\>\Phi^{(k)}+\frac{\Delta
t}{}\mathit{\mathbb{L}}(\Phi^{(n)})\\\
\Phi^{(2)}=&\>\>\>\frac{3}{4}\Phi^{(k)}+\frac{1}{4}\Phi^{(1)}+\frac{1}{4}\frac{\Delta
t}{}\mathit{\mathbb{L}}(\Phi^{(1)})\\\
\Phi^{(n+1)}=&\>\>\>\frac{1}{3}\Phi^{(k)}+\frac{2}{3}\Phi^{(2)}+\frac{2}{3}\frac{\Delta
t}{}\mathit{\mathbb{L}}(\Phi^{(2)})\end{split}$ (14)

The time-step $\Delta t$ for the simulations is determined by considering the
diffusive stability requirements and the Courant– Friedrichs–Lewy (CFL)
condition arising from the convective velocity and the artificial acoustic
wave speed [21, 41]. Moreover, the inclusion of the bulk viscosity tensor in
the GPE+BV solver adds another restriction on the time step. Theoretically,
for two-dimensional simulations, a minimum over the following bounds
determines the time step, as given below.

$\Delta t=f_{s}\cdot\textrm{min}(\Delta t_{acs},\Delta t_{ad},\Delta
t_{diff},\Delta t_{visc},\Delta t_{bv}),$ (15)

where $0<f_{s}\leq 1$ is the safety factor. The criteria appearing in the
above equation are listed here.

  1. 1.

CFL criterion based on the acoustic wave speed

$\Delta t_{acs}=\min\limits_{i,j}\left[\frac{U/Ma}{\Delta
x_{i}}+\frac{U/Ma}{\Delta y_{j}}\right]^{-1}$ (16)

  2. 2.

CFL criterion based on the convective velocity

$\Delta t_{ad}=\min\limits_{i,j}\left[\frac{|u_{i,j}|}{\Delta
x_{i}}+\frac{|v_{i,j}|}{\Delta y_{j}}\right]^{-1}$ (17)

  3. 3.

Stability criterion based on the viscous term in the momentum equation

$\Delta
t_{visc}=\min\limits_{i,j}\left[\frac{1}{2}\frac{\Delta_{i,j}^{2}}{(1/Re)}\right]$
(18)

  4. 4.

Stability criterion based on the diffusion term in the GPE

$\Delta
t_{diff}=\min\limits_{i,j}\left[\frac{1}{2}\frac{\Delta_{i,j}^{2}}{(1/RePr)}\right]$
(19)

  5. 5.

Time-step restriction due to the bulk viscosity term

$\Delta
t_{bv}=\min\limits_{i,j}\left[\frac{1}{2}\frac{\Delta_{i,j}^{2}}{\mathcal{B}}\right]$
(20)

In the above expressions, $\min\limits_{i,j}$ represents the minimum of all
values computed for each cell on the mesh. Note that all the parameters in the
above equations are in non-dimensional form. $U$ denotes the characteristic
velocity of the flow and $\Delta$ represents the grid-size factor for 2D grids
given by,

$\Delta^{2}_{i,j}=\frac{\Delta x_{i}^{2}\Delta y_{j}^{2}}{\Delta
x_{i}^{2}+\Delta y_{j}^{2}}.$

From the above equations, it is evident that the time-step depends on the flow
parameter $Re$ and the three numerical parameters such as $Ma$, $Pr$ and
$\mathcal{B}$. Unlike $Re$, which is problem-dependent, the numerical
parameters can be adjusted. For weakly compressible approaches, $Ma\ll 1$ so
as to closely represent incompressible flows. Hence, $\Delta t_{acs}$ in most
simulations dictates the time step of the simulation. A lower value of $Pr$
helps in damping the acoustic waves [17] but at the cost of a lower time-step
requirement. Similarly, a large value of $\mathcal{B}$ implies better damping,
but requires a reduction in $\Delta t$. Hence, there is a trade-off between
the computational cost and the rate of acoustic damping. In this work, we
chose a large permissible value of $\mathcal{B}$ for the given allowable
$\Delta t_{acs}$. This ensures the maximum benefit of including the bulk
viscosity without increasing the computational cost.

##  3 Numerical experiments

The present paper aims to comprehensively investigate the effect of artificial
bulk viscosity on weakly compressible methods. In this regard, we present a
few benchmark examples covering free-shear flows and external flows over a
solid. For each example, we perform simulations without (denoted as GPE
solver) and with the artificial bulk viscosity (denoted as GPE+BV solver).

Each test case reported in this section investigates a different influence of
the inclusion of bulk viscosity. We start with the simple lid-driven cavity
problem in order to provide evidence for the damping of acoustic waves by the
bulk viscosity. Then, we discuss how the present approach yields reduced mass
conservation error using the doubly periodic shear layer (DPSL). Further, we
present the simulation data from the Taylor-Green vortex (TGV) problem to
showcase how the inclusion of this additional term helps to obtain a uniform
order of convergence. Next, we demonstrate for the flow over a square cylinder
that the proposed form of the bulk viscosity term helps us to quickly reach
the periodic vortex shedding, thereby achieving improved computational
efficiency. Finally, for an impulsively started plate, a challenging unsteady
problem for weakly compressible approaches, we provide evidence for better
force computation using the proposed approach. It is to be noted that the
first three problems are simulated on uniform meshes, on which the bulk
viscosity is a scalar; the last two test cases are on stretched meshes, and
the advantages of using the tensorial form of the bulk viscosity are
underlined.

Toutant [17] reported that by setting $Ma=0.02$ and $Pr=1$, GPE produced
accurate results for incompressible flows. Unless otherwise stated, we use
these values in our simulations.

###  3.1 Lid-driven cavity flow: Evidence of damping by bulk viscosity

Figure 2: Geometric specifications of lid-driven cavity flow problem. P is the
point at which the evolution of pressure with respect to time is studied. It
is located at $x=L/2$ and $y=H/2$.

We illustrate the presence of the acoustic pressure waves and their
suppression using the bulk viscosity within the well-known lid-driven cavity
problem. Although this is a steady-state problem, our objective is to observe
the evolution of the flow field until the steady-state is reached. The
schematic of the problem and the boundary conditions are shown in figure 2.
Throughout the domain, both the velocity and pressure fields are initialised
to zero. A grid of $64\times 64$ and a time-step of $10^{-4}$ is used. We
performed simulations at two different Reynolds numbers, viz. $Re=100$ and
$Re=400$.

Figure 3: Comparison of pressure evolution at the centre point of the domain,
with GPE and GPE+BV solvers for the lid-driven cavity flow problem (a)
$Re=100$ (b) $Re=400$.

Firstly, we simulated the flow without considering the artificial bulk
viscosity term. i.e. with the GPE solver. Figure 3 shows the time variation of
the pressure at the cavity centre (denoted by point P in figure 2). In these
plots, the time is scaled with $t_{c}$, which is the time required to achieve
the steady state. The highly oscillatory nature of the pressure field, seen in
the figure, can be attributed to the inherent acoustic waves associated with
the weakly compressible family of methods. It is clear that oscillations of
high amplitude are visible at the beginning of the simulation, which die out
eventually.

The cavity problem was simulated again with the GPE+BV solver to check the
effect of incorporating the artificial bulk viscosity term. The pressure
plots, thus obtained, are presented in figure 3. We used $\mathcal{B}=50\Delta
x$, with $\Delta x$ being the size of the uniform mesh. We observe a rapid
damping of pressure oscillations with this minor modification of the solver.
When the bulk viscosity term is ignored, the presence of oscillations is
visible even at $t/t_{c}=0.3$ for $Re=100$. In contrast, with the GPE+BV
solver, complete damping is attained at around $t/t_{c}=0.05$. Similarly, for
$Re=400$, complete damping is observed at $t/t_{c}=0.3$ and $t/t_{c}=0.03$ for
GPE and GPE+BV solvers, respectively. A previous study [29] mentioned such
acoustic wave annihilation using bulk viscosity in the lid-driven cavity
problem. However, they did not present any pressure evolution plots to
visualise the damping.

Figure 4:  Effect of artificial bulk viscosity parameter ($\lambda$) for the
lid-driven cavity flow problem (a) $Re=100$ (b) $Re=400$.

It is of interest to investigate the influence of $\mathcal{B}$ on the
effectiveness of damping. Since $\Delta x$ is constant throughout the domain,
we changed $\lambda$ to vary $\mathcal{B}$. The time-evolution of pressure at
the cavity centre for different values of $\lambda$ are presented in figure 4.
It is directly evident that increasing the value of $\lambda$ leads to the
elimination of oscillations at a faster pace. However, as can be seen from
equation (20), the allowable $\Delta t$ is inversely proportional to
$\mathcal{B}$. For the current problem, setting $\lambda=50$ does not show any
stability issues for the chosen time step. However, for $\lambda=100$ and
$\lambda=500$, we need to reduce the time-step to $10^{-5}$ and $10^{-6}$,
respectively. This increases the computational cost but provides a quicker
damping rate. Hence, the value of $\mathcal{B}$ should be based on a pay-off
between the computational expenses and the rate of damping the acoustic
pressure waves. In all other simulations, we choose $\mathcal{B}$ so that the
simulation is stable for $\Delta t$ decided by the $Ma$ as given in equation
(16). This is done so as to get the maximum benefit of acoustic wave
suppression without introducing additional computational overhead by reducing
the time step.

Pressure oscillations in the lid-driven cavity flow have been presented
previously using similar weakly compressible flow solvers such as classical
ACM [42] and lattice Boltzmann methods [42, 43]. Similar pressure plots
presented by [43] indicate that the pressure oscillations are visible even at
$t/t_{c}=0.7$ for $Re=100$, whereas for GPE and GPE+BV solvers, the
oscillations are eliminated at $t/t_{c}=0.3$ and 0.05, respectively. We can
conclude from this that the pressure diffusion term in GPE provides damping of
acoustic waves as envisaged by Toutant [16], but the present modification
enhances the damping significantly.

Figure 5: Vorticity contours for the doubly periodic shear layer test case
with a grid of $512\times 512$ at $t=1$ using the GPE+BV solver.

###  3.2 Doubly periodic shear layer: Effect on velocity divergence

In the previous test case, which was a steady state problem, we found that the
GPE+BV solver can more effectively damp the artificial acoustic waves. We
extend our analysis further to investigate how the addition of bulk viscosity
affects the flow field of an unsteady problem. For this purpose, the doubly
periodic shear layer (DPSL) on a unit square domain [$1\times 1$] is chosen
[44, 2, 45]. The initial conditions are as follows

$\displaystyle u$
$\displaystyle=\begin{cases}\tanh(\rho(y-0.25))&\text{if}\hskip
5.69046pt{y\leq 0.5}\\\ \tanh(\rho(0.75-y))&\text{otherwise}\end{cases}$ (21a)
$\displaystyle v$ $\displaystyle=\delta\sin(2\pi(x+0.25))$ (21b)
$\displaystyle p$ $\displaystyle=0.$ (21c)

In the above equation, the value of $\rho$ determines the thickness of the
shear layer, and $\delta$ represents the amplitude or strength of the initial
perturbation. We set $\rho=80$ and $\delta=0.05$ and use a uniform grid of
$512\times 512$, which is necessary to avoid spurious vortices [46, 17]. We
set $Re=10^{4}$, and the simulation is run until $t=1$, with $\Delta
t=10^{-5}$.

Due to the initial perturbation, the shear layers roll up into two separate
vortical structures, and these are presented in figure 5 obtained using GPE+BV
solver with $\mathcal{B}=50\Delta x$ at $t=1$. This qualitatively resembles
the vorticity plots, without any spurious vortices, reported in literature
[15, 47, 17].

The significant advantage of using GPE+BV over GPE can be seen from figure 6,
which illustrates the mass conservation error in the whole domain at $t=1$.
While GPE produced an error of the order $10^{-2}$, the artificial bulk
viscosity present in GPE+BV reduced this error to $10^{-4}$. Thus, a reduction
in error of two orders of magnitude is achieved by the addition of
$\mathcal{B}$. A major limitation of weakly compressible methods, like GPE, is
the large mass conservation error, as these methods don’t enforce strict
incompressibility. Here, by adding a simple term to the momentum equation and
without any complicated modifications, we are able to downscale this error by
two orders of magnitude. To the best of our knowledge, such a large reduction
in mass conservation error is not reported in the literature.

Figure 6: Contours of velocity divergence for the doubly periodic shear layer
test case with a grid of $512\times 512$ at $t=1$ using (a) GPE solver (b)
GPE+BV solver.  Figure 7: Evolution of average kinetic energy ($E_{k}$) for
the doubly periodic shear layer test case with a grid of $512\times 512$. The
computed values are compared with that reported by Minion and Brown [45] and
Clausen [15].

Having demonstrated the reduced mass conservation error, it is also important
to verify that the addition of $\mathcal{B}$ does not induce additional error
in any other form. In order to check this, we computed the evolution of
average kinetic energy,

$E_{k}=\frac{1}{}\int_{\left(\mathbf{u}\cdot\mathbf{u}\right)}{2}d$ (22)

where
$correspondstothetotalvolumeofthecomputationaldomain.Thevariationof$E_k$withtimeobtainedfromGPE+BVandreferenceresultsarepresentedinfigure~{}\ref{dpsl_keVsTime}.TheexactagreementofresultsoftheGPE+BVsolverandtheotherweaklycompressibleapproachwithoutthebulkviscosity~{}\cite[cite]{[\@@bibref{Number}{Clausen2013}{}{}]},andfromthesolutionofincompressibleNavier-
Stokesequation~{}\cite[cite]{[\@@bibref{Number}{minion1997}{}{}]}underlinesthattheadditionof$B$doesnotintroduceanyothererror.Thisisduetothefactthatbulkviscositydirectlyworkswiththedivergenceofvelocityfieldonly,leavingallotherincompressiblecomponentsintact~{}\cite[cite]{[\@@bibref{Number}{sun2023}{}{}]}.\par
Thevelocitydivergenceoftheweaklycompressibleapproachescannotbeimprovedwithmeshrefinement,whichwasdiscussedbyToutant~{}\cite[cite]{[\@@bibref{Number}{Toutant2018}{}{}]}.Toaddressthis,heratherchosetoincreasethediffusionofartificialacousticwavesbyreducing$Pr$.Fromhisparametricstudyon$Pr$,theminimumvalueofvelocitydivergenceachievedwasoftheorder$10^-3$with$Pr=0.01$.However,asgiveninequation~{}\eqref{time_step_eqn_Pr},thismandatesareductionin$Δt$.Incontrast,withtheGPE+BVsolverwith$Pr=1$,weareabletoattainvelocitydivergenceoftheorder$10^-4$withoutreducing$Δt$.Thisconfirmsthatthebulkviscosityservesasacomputationallymoreefficientproceduretoannihilateartificialacousticwaveswhencomparedtothepressurediffusionterm.\par
TheFEMcommunityhasalsousedasimilarapproachofincludingthebulkviscosityterm~{}\cite[cite]{[\@@bibref{Number}{de_mulder1998,rasthofer2018}{}{}]},whichtheycall`grad-
div^{\prime}stabilisation.Evenfortheincompressibleflowformulations,ithasbeenreportedthatsuchamodificationgivesareducedvelocitydivergenceerror~{}\cite[cite]{[\@@bibref{Number}{olshanskii2009}{}{}]}.\par$

###  3.3 Laminar Taylor-Green vortex problem: Effect on the order of
convergence

The Taylor Green vortices (TGV) [51] is a widely used test case for verifying
the order of accuracy of numerical algorithms. Our objective here is to
investigate how the artificial bulk viscosity term affects the order of
convergence. The analytical solution for this problem, in a doubly-periodic
unit square, is given by the following equation [17], and the initial
conditions are obtained from the same by setting $t=0$.

$\displaystyle\text{u}(x,y,t)$ $\displaystyle=\text{cos}(2\pi
x)\text{sin}(2\pi y)e^{-\frac{8\pi^{2}}{Re}t}$ (23a)
$\displaystyle\text{v}(x,y,t)$ $\displaystyle=-\text{sin}(2\pi
x)\text{cos}(2\pi y)e^{-\frac{8\pi^{2}}{Re}t}$ (23b)
$\displaystyle\text{p}(x,y,t)$ $\displaystyle=-\frac{1}{4}(\text{cos}(4\pi
x)+\text{cos}(4\pi y))e^{-\frac{16\pi^{2}}{Re}t}$ (23c)

Figure 8: Convergence plots for the Taylor-Green vortex problem, with the GPE
and GPE+BV solvers (a) $x-$component of velocity (b) $y-$component of velocity
(c) pressure.

The TGV case is simulated using GPE and GPE+BV solvers. Similar to the
previous test problems, we use a constant value of $\mathcal{B}=50\Delta x$.
We set $Re=100$ and $\Delta t=10^{-5}$.

Convergence plots at $t=1$ for velocity components and pressure are presented
in figure 8, which shows the variation of $L_{\infty}$-norm of error with the
number of grid points ($N$). It can be seen that velocity components exhibit
second-order convergence, as expected, for both solvers. However, the order of
convergence for pressure exhibits some irregularity for coarser grids when the
GPE solver is used. In contrast, we observe uniform second-order convergence
for pressure with the GPE+BV solver for the entire range of grid sizes under
consideration.

Toutant [17] also reported the same behaviour for the GPE solver. He deduced
the convergence rate of velocity and pressure to be 2 and 1.87, respectively.
In comparison, we find that GPE+BV enables a consistent second-order
convergence rate for pressure also.

###  3.4 Flow over a square cylinder: Time to reach periodic vortex shedding

The unsteady problems considered in the earlier sections were free-shear
flows. In this section, we consider the fluid flow over a square cylinder, an
example that involves large-scale boundary layer separation. In this and the
next example, the domain is covered with a stretched mesh; thus, the
artificial bulk viscosity becomes non-homogeneous and anisotropic as given by
equation (7). The objective is to investigate the proposed form of
$\boldsymbol{\mathcal{B}}$ on (i) the effectiveness in suppressing the
acoustic waves, thereby reaching the periodic shedding state quickly, and (ii)
accuracy of computing integral quantities over a solid surface, such as drag
and lift forces.

Figure 9 shows the problem’s schematic and boundary conditions. We use a non-
uniform mesh of $526\times 442$, and the complete mesh details are given in
table 1. The geometry, boundary conditions and mesh size around the cylinder
are the same as that of Sharma and Eswaran [52]. The flow field is initialised
to $u=v=p=0$, and the time step of $10^{-4}$ is used to simulate the flow
until the periodic state is established. The flow is defined by $Re=100$ and
the characteristic velocity, $U=1$.

Figure 9: Geometric specifications and boundary conditions of flow over a square cylinder. Table 1: Details of the grid for the case of flow over a square cylinder. Grid size along $\boldsymbol{x}$ | Grid size along $\boldsymbol{y}$
---|---
$\boldsymbol{x}$ | 0 | 8.25 | 9.75 | 26 | $\boldsymbol{y}$ | 0 | 9.25 | 10.75 | 20
$\boldsymbol{\Delta x}$ | 0.25 | 0.008334 | 0.008334 | 0.25 | $\boldsymbol{\Delta y}$ | 0.25 | 0.008334 | 0.008334 | 0.25

We first discuss the results obtained from the GPE solver. After some initial
transients, we observe the vortex shedding, which can be graphically
visualised from the lift coefficient ($C_{L}$) and the drag coefficient
($C_{D}$) plots, shown in figure 10. The high frequency oscillations present
in the plots indicate the acoustic waves visually, and with time they get
damped out. To examine whether the periodic state is achieved, we plot $C_{L}$
and $C_{D}$ at different time intervals and are presented in figure 10. We
observe that the pressure waves are highly prominent around $t=200$, which
gets almost filtered out by $t=900$. However, it is only around $t=1200$ that
a complete periodic flow is attained. It can be seen that $C_{D}$,
predominantly contributed by the pressure drag, requires a longer time to
stabilise than $C_{L}$. Thus, the delay in achieving periodic flow is
attributed to the slower rate of damping the acoustic pressure oscillations.

Figure 10: The evolution of lift and drag coefficients with respect to time
for the flow over a square cylinder problem using the GPE solver ($Pr=1$). (a)
Plots till $t=500$ (b) Plots at different time intervals to check periodicity.

The GPE solver itself has the following mechanisms to dissipate the acoustic
waves: viscosity, numerical diffusion, and pressure diffusion. Out of these,
only the pressure diffusion term can be specified as an input parameter by
adjusting the value of $Pr$ (equation (5)). Before inspecting the prospects of
using the GPE+BV solver, it is important to examine the effectiveness of
increasing the pressure diffusion in annihilating the acoustic waves.
Accordingly, we simulated the test case with $Pr=0.01$. The corresponding
$C_{L}$ and $C_{D}$ plots are presented in figure 11. Compared to the results
shown in figure 10 with $Pr=1$, we observe a slightly faster damping of the
pressure oscillations. i.e., the time taken to reach the periodic state is
reduced to $t=900$, from $t=1200$ when $Pr=1$. However, this improvement is
achieved with a higher computational cost since a lower time-step ($\Delta
t=10^{-5}$) is necessary for stability when $Pr=0.01$. Evidently, faster
damping can be achieved with a lower $Pr$ but at a higher computational cost,
which is undesirable.

Figure 11: The evolution of lift and drag coefficients with respect to time
from the simulation of the flow over square cylinder problem using GPE solver
($Pr=0.01$).

Next, we examine the capability of the GPE+BV solver. As discussed earlier,
for a non-uniform mesh, $\boldsymbol{\mathcal{B}}$ is non-homogeneous and
anisotropic. In order to isolate the influence of the non-homogeneity and the
anisotropic nature, we performed simulations with the following variants of
$\boldsymbol{\mathcal{B}}$.

  1. 1.

Homogeneous isotropic function,
$\boldsymbol{\mathcal{B}}=\lambda\Delta_{\textrm{min}}\bf{I}$, where
$\Delta_{\textrm{min}}$ is the smallest mesh size of the entire domain. For
the considered mesh, $\Delta_{\textrm{min}}=0.008334$ (refer table 1)

  2. 2.

Non-homogeneous isotropic function,
$\boldsymbol{\mathcal{B}}=\tilde{\lambda}\sqrt{\Delta x^{2}+\Delta
y^{2}}\bf{I}$, where $\Delta x$ and $\Delta y$ denote the size of the control
volume under consideration.

  3. 3.

Non-homogeneous anisotropic function,
$\boldsymbol{\mathcal{B}}=\begin{bmatrix}\mathcal{B}^{\mathcal{X}}&0\\\
0&\mathcal{B}^{\mathcal{Y}}\end{bmatrix}$, where,
$\mathcal{B}^{\mathcal{X}}=\lambda\Delta x$ and
$\mathcal{B}^{\mathcal{Y}}=\lambda\Delta y$.

We use $\lambda=37.5$, the highest value that can be used without reducing the
time step. In order to ensure stability for the variant-2, we set
$\tilde{\lambda}=\lambda/AR_{max}$, where $AR_{max}$ is the maximum aspect
ratio of the grid. For both the non-homogeneous variants, the distribution of
bulk viscosity in the domain is presented in figure 12.

Figure 12: Spatial variation of artificial bulk viscosity coefficient. (a)
Non-homogeneous isotropic function,
$\boldsymbol{\mathcal{B}}=\tilde{\lambda}\sqrt{\Delta x^{2}+\Delta
y^{2}}\bf{I}$ (b) Components of non-homogeneous anisotropic function
$\mathcal{B}^{\mathcal{X}}=\lambda\Delta x$,
$\mathcal{B}^{\mathcal{Y}}=\lambda\Delta y$.

The first variant we considered is the homogeneous isotropic function. We
observe a drastic reduction in pressure oscillations and faster damping
compared to the GPE solver, as evident from figure 13. The periodic state is
achieved at around $t=400$, which is a significant improvement when compared
to the GPE.

Figure 13: Plots of $C_{D}$ and $C_{L}$ for the flow over a square cylinder
problem using the homogeneous isotropic variant of GPE+BV solver. (a) Plots
till $t=500$ (b) Plots at different time intervals to check periodicity.

In figure 14, we compare $C_{D}$ plots of all three variants between $t=200$
and $t=220$. Both the non-homogeneous isotropic as well as the non-homogeneous
anisotropic variants achieve the periodic state in this time frame. This
implies significantly more effective damping of the acoustic waves when
compared to the isotropic homogeneous version, which is employed by all the
existing weakly compressible approaches.

Although variants-2 and 3 achieved periodic states approximately at the same
time, inspection of $C_{D}$ variation during the earlier stages of flow
development underlines the superior damping characteristics of the anisotropic
variant, as evident from figure 15. It can be seen that while the non-
homogeneous isotropic variant still exhibits unwanted oscillations, the
anisotropic variant has almost completely eliminated the pressure acoustics.
The consequence of this enhanced performance will be demonstrated for a truly
unsteady problem in the next section.

Figure 14: Comparison of $C_{D}$ for all three variants of
$\boldsymbol{\mathcal{B}}$, for the flow over square cylinder problem. Figure
15: $C_{D}$ plot for the flow over a square cylinder problem from $t=100$ to
$t=180$, using the GPE+BV solver obtained using non-homogeneous isotropic and
non-homogeneous anisotropic variants of $\boldsymbol{\mathcal{B}}$.

The computational cost is directly proportional to the time taken to achieve
the steady state ($t_{Periodic}$) and the time step. The effect of damping the
acoustic waves directly translates into $t_{Periodic}$, and table 2 summarises
the same for different variants of $\boldsymbol{\mathcal{B}}$. We normalised
the computational costs of all the cases with that incurred for the proposed
anisotropic variant. The non-homogeneous variants, which are proposed for the
first time in this paper, lead to a reduction of computational cost by at
least six times compared to the existing GPE solver.

Table 2: Comparison of time to reach periodic phase using various schemes. Scheme | $\boldsymbol{\Delta t}$ | $\boldsymbol{t_{Periodic}}$ | Cost
---|---|---|---
GPE | $10^{-4}$ | 1200 | 6
GPE+BV ($\boldsymbol{\mathcal{B}}=\lambda\Delta_{\textrm{min}}\bf{I}$) | $10^{-4}$ | 400 | 2
GPE+BV ($\boldsymbol{\mathcal{B}}=\tilde{\lambda}\sqrt{\Delta x^{2}+\Delta y^{2}}\bf{I}$) | $10^{-4}$ | 200 | 1
GPE+BV ($\mathcal{B}^{\mathcal{X}}=\lambda\Delta x$ and $\mathcal{B}^{\mathcal{Y}}=\lambda\Delta y$) | $10^{-4}$ | 200 | 1

Before concluding, we quantify the accuracy of the proposed GPE+BV solver. For
this, we compare the average drag coefficient ($C_{D_{avg}}$), the RMS value
of the lift coefficient ($C_{L_{rms}}$) and the Strouhal number ($St$)
obtained for all the variants discussed in this section. These values
presented in table 3 show that all three variants produce accurate results.
Specifically, the non-homogeneous anisotropic variant exhibits a deviation of
around $2.6\%$, $0.1\%$ and $0.8\%$ respectively for $C_{D_{avg}}$,
$C_{L_{rms}}$ and St, when compared to Sharma and Eswaran [52]. Thus, with the
proposed approach, we achieved sufficient acoustic wave damping without
affecting the solver’s accuracy.

Table 3: Computed values of average drag coefficient ($C_{D_{avg}}$), RMS of lift coefficient ($C_{L_{rms}}$) and Strouhal number ($St$) for the flow over square cylinder problem. $\mathbf{Reference}$ | $\boldsymbol{C_{D_{avg}}}$ | $\boldsymbol{C_{L_{rms}}}$ | $\boldsymbol{St}$
---|---|---|---
Sohankar [53] | 1.4770 | 0.1560 | 0.1460
Sharma and Eswaran [52] | 1.4936 | 0.1922 | 0.1488
Sen et al. [54] | 1.5287 | 0.1928 | 0.1452
Zhao et al. [55] | 1.47 | 0.156 | 0.146
Saha and Shrivastava [56] | 1.524 | - | 0.163
GPE ($Pr=1$) | 1.5270 | 0.1986 | 0.1480
GPE ($Pr=0.01$) | 1.5322 | 0.1997 | 0.1480
GPE+BV ($\boldsymbol{\mathcal{B}}=\lambda\Delta_{\textrm{min}}\bf{I}$) | 1.5311 | 0.1948 | 0.15
GPE+BV ($\boldsymbol{\mathcal{B}}=\tilde{\lambda}\sqrt{\Delta x^{2}+\Delta y^{2}}\bf{I}$) | 1.5341 | 0.1935 | 0.15
GPE+BV ($\mathcal{B}^{\mathcal{X}}=\lambda\Delta x$ and $\mathcal{B}^{\mathcal{Y}}=\lambda\Delta y$) | 1.5336 | 0.1924 | 0.15

###  3.5 Flow over an impulsively started thin plate: Effectiveness in a truly
unsteady problem

The previous case was focused on the periodic vortex shedding behind a
cylinder, which occurs far away from the initial time. Hence, the effect of
artificial acoustic waves in the initial transients did not play a role in the
desired periodic state, as can be interpreted from figure 15. In this section,
we investigate the evolution of twin vortices behind an impulsively started
thin plate. Our focus is to study in detail the initial stages of development
of the flow field and the rapidly varying force coefficients. This is a highly
challenging test case for weakly compressible approaches, and we will show
that the anisotropic nature of bulk viscosity proposed in this paper provides
the best possible results.

Figure 16: Geometric specifications and boundary conditions of flow over an
impulsively started plate.

Figure 16 shows the schematic of the problem, which is adapted from Kiris and
Kwak [11]. The $x$-component of velocity, $u$, is initialized with a constant
value of $U=1$ throughout the flow domain to enforce the impulsive starting,
and the simulation is carried out at $Re=126$. We use a stretched grid of
$679\times 922$, and the complete mesh details are given in table 4. The mesh
is sufficiently finer around the thin plate to capture the vortex structures
properly. The problem is simulated till time $t=10$ with $\Delta t=10^{-5}$.
In order to validate the results from our solver, we simulated this example
using the popular open-source software OpenFOAM (denoted as INS in all the
plots), which solves the incompressible Navier-Stokes equations. This serves
as a reliable reference since the artificial acoustic effects are absent in
the ‘INS’ results.

Table 4: Details of the grid for the case of flow over an impulsively started thin plate. Mesh size along $\boldsymbol{x}$ | Mesh size along $\boldsymbol{y}$
---|---
$\boldsymbol{x}$ | 0 | 3.75 | 4.25 | 4.78 | 5.53 | 20 | $\boldsymbol{y}$ | 0 | 4.5 | 5.25 | 6.75 | 7.5 | 12
$\boldsymbol{\Delta x}$ | 0.1 | 0.025 | 0.0025 | 0.0025 | 0.025 | 0.1 | $\boldsymbol{\Delta y}$ | 0.1 | 0.025 | 0.0025 | 0.0025 | 0.025 | 0.1

Firstly, we simulated the test case with the GPE solver. The plot of bubble
length ($L_{B}$) with respect to time is presented in figure 17. Even though
the GPE solver is able to predict $L_{B}$ well, the plot exhibits minor
fluctuations. This is due to the acoustic waves appearing in the flow field,
which induces drastic oscillations in the drag coefficient as shown in figure
17. As mentioned in § 1, we can see that the artificial sound waves cause
spurious oscillations in $C_{D}$ to the extent that even its overall variation
is untraceable.

Figure 17: Evolution of (a) bubble length ($L_{B}$) compared with [57, 58];
and (b) drag coefficient ($C_{D}$) for the flow over an impulsively started
thin plate case, using the GPE solver.

Similar to the previous test cases, reducing $Pr$ leads to suppression of
acoustic waves, but at the cost of increased computational time and reduced
accuracy. These results are summarised in B.

Next, we simulate this test case using all three variants of bulk viscosity
presented in the previous section and the results are presented in figure 18.
For all the variants, we use the constant $\lambda=79$. The following
observations can be made from the figure.

  * 1.

All three variants accurately capture $L_{B}$ as can be seen from figure
18(a).

  * 2.

Both the homogeneous isotropic (figure 18(b)) and non-homogeneous
isotropic(figure 18(c)) are ineffective.

  * 3.

The non-homogeneous anisotropic variant of $\boldsymbol{\mathcal{B}}$
completely eliminated the artificial acoustic waves after $t=2$ as is evident
from figure 18(d).

We stress that although for the previous test case, the non-homogeneous
isotropic variant produced the same periodic state as that of the non-
homogeneous anisotropic case, to capture the dynamics in the early stages of
flow development for this truly unsteady problem, the anisotropic nature of
bulk viscosity is crucial. This is consistent with the results reported in
figure 15.

Figure 18: Evolution of (a) bubble length; and drag coefficient for the flow
over an impulsively started thin plate case using the GPE+BV solver with
different variants of bulk viscosity coefficient ($\boldsymbol{\mathcal{B}}$)
such as (b) homogeneous isotropic (c) non-homogeneous isotropic (d) non-
homogeneous anisotropic.

The pressure contours at $t=2$ presented in figure 19 provide visual evidence
for the enhanced performance of the anisotropic variant of bulk viscosity.
Consistent with figure 18, while the anisotropic variant fully suppressed the
acoustic waves, the signature of these waves is clearly seen for GPE and the
other isotropic variants.

Figure 19: Pressure field at $t=2$ for the flow over an impulsively started
thin plate test case using (a) GPE solver; and the GPE+BV solver with
different variants of bulk viscosity coefficient ($\boldsymbol{\mathcal{B}}$)
such as (b) homogeneous isotropic (c) non-homogeneous isotropic (d) non-
homogeneous anisotropic.

From figure 18, we observe that as time progresses, the computed $C_{D}$
values obtained from the GPE+BV solver approach those of the reference solver.
To provide a quantitative estimate of the error, we listed $C_{D}$ obtained
from GPE+BV and the reference INS simulation in table 5. For the entire time
range reported, our method provides an accurate drag coefficient.

Table 5: Percentage deviation of drag coefficient computed using GPE+BV compared to that obtained from INS. $\boldsymbol{t}$ |  $\boldsymbol{C_{D}}$ (INS) |  $\boldsymbol{C_{D}}$ (GPE+BV) | % Deviation
---|---|---|---
2 | 2.2833 | 2.2787 | 0.199
3 | 2.0106 | 2.0341 | 1.167
4 | 1.8881 | 1.9089 | 1.102
5 | 1.8206 | 1.8327 | 0.665
6 | 1.7707 | 1.7822 | 0.648
7 | 1.7357 | 1.7443 | 0.497
8 | 1.7096 | 1.7168 | 0.422
9 | 1.6894 | 1.6955 | 0.358
10 | 1.6734 | 1.6786 | 0.314

To qualitatively analyse the unsteady flow field, the evolution of vorticity
contours is illustrated in figure 20. This plot enables a robust one-to-one
comparison between results obtained from GPE+BV and INS solvers. We observe
that two counter-rotating vortices appear in the wake region behind the plate.
Over time, these vortices increase in size as they are constantly fed by the
shear layer appearing at the tip of the plate [58]. An excellent agreement
between these results can be seen after $t=2$. A similar observation is made
by Sun et al. [37] for the viscous flow past an inclined elliptic cylinder
problem where the lift and drag coefficient plots exhibit high fluctuations
within non-dimensional time of $10$, even with the inclusion of an (isotropic)
acoustic damper term in their smoothed-particle hydrodynamics solver.

Figure 20: Comparison of vorticity contours from GPE+BV and INS solvers, at
various time instances for the flow over an impulsively started thin plate
test case. The contours plotted are for $\omega_{z}=\pm 1,\pm 2,\pm 3,\pm
4,\pm 5,\pm 6,\pm 10$.

To provide additional confirmation for the effectiveness of the proposed
approach, we investigate its performance at a higher $Re$. We simulated the
impulsively started plate at $Re=1000$ on the same mesh and with the same
time-step using GPE and GPE+BV solver with non-homogeneous anisotropic bulk
viscosity; the evolution of $L_{B}$ and $C_{D}$ are plotted in figure 21. The
accuracy of our method at this high $Re$ is also evident from the figure. This
proves the robustness of the proposed approach. Similar to $Re=126$, the
isotropic variants of bulk viscosity were ineffective.

It is to be noted that, compared to the damping mechanism using pressure
diffusion, the bulk viscosity approach is better suited for high Re flows.
This is because the magnitude of the pressure diffusion, which is decided by
its coefficient, $1/RePr$, reduces when $Re$ increases, whereas
$\boldsymbol{\mathcal{B}}$ is unaltered with higher $Re$.

Figure 21: Evolution of (a) bubble length and (b) drag coefficient for the
flow over an impulsively started thin plate test case, with $Re=1000$.

For the effective damping of acoustic waves, the time scale associated with
bulk viscosity should be smaller than the acoustic time scale. In contrast,
for the above examples, we chose both these time scales are of the same order,
so that the same $\Delta t$ can be used for both GPE and GPE+BV solvers.
Despite this fact, the non-homogeneous anisotropic GPE+BV variant suppresses
acoustic waves effectively, except for $t<2$. Although this can be viewed as a
drawback of the present scheme, we highlight that our approach outperforms the
isotropic bulk viscosity that is widely used in weakly compressible solvers.

##  4 Conclusions

Weakly compressible approaches to simulate incompressible flows have gained a
lot of popularity in recent years because of their potential to achieve
scalability in HPC implementations. A major limitation of such methods is the
presence of artificial acoustic waves that act as the source of large mass
conservation errors and highly oscillatory pressure fields. Although pressure
diffusion and bulk viscosity are employed in many works to address this issue,
their relative effectiveness remains unclear. In this paper, we demonstrated
using various numerical examples that the bulk viscosity is more effective in
annihilating acoustic waves than the pressure diffusion. We presented a
comprehensive study of bulk viscosity on various aspects of weakly
compressible methods, including damping of pressure oscillations, order of
convergence and mass conservation error. A key contribution of this work is
the formulation of a non-homogeneous anisotropic form of the bulk viscosity.
We proved that the chosen form exhibited superior performance in annihilating
the acoustic waves, even in the initial transients of unsteady flows. Such an
example is highly challenging for weakly compressible methods, and we
underline that the proposed bulk viscosity was significantly better in
predicting the forces acting on a solid body. We carried out our
investigations using a specific weakly compressible method called the general
pressure equation. However, our findings are not restricted to this method
alone, as the results we obtained from the entropically damped artificial
compressibility based solver (not presented in this paper) also lead to
similar conclusions. Moreover, bulk viscosity has been commonly used in the
lattice Boltzmann methods and smoothed particle hydrodynamics based solvers.
We envisage that these widely used weakly compressible methods will also
benefit from the non-homogeneous anisotropic form of the bulk viscosity
proposed in this paper.

## CRediT authorship contribution statement

Dheeraj Raghunathan: Data curation, Investigation, Methodology, Software,
Formal analysis, Validation, Visualization, Writing – original draft. Y.
Sudhakar: Conceptualization, Methodology, Resources, Supervision, Writing –
review & editing.

## Declaration of competing interest

The authors declare that they have no known competing ﬁnancial interests or
personal relationships that could have appeared to inﬂuence the work reported
in this paper.

## Data availability

The data supporting this study’s findings are available from the corresponding
author upon reasonable request.

## Acknowledgments

We acknowledge the financial support provided by IIT Goa in the form of a
startup grant. Simulations are run on computational resources setup from the
DST-SERB Ramanujan fellowship (SB/S2/RJN-037/2018). The first author would
like to acknowledge the contribution endowed by Mr. Smruti Ranjan Jena for
simulating the flow over an impulsively started thin plate using OpenFOAM.

## References

  * [1] S. Patankar, Numerical heat transfer and fluid flow.  Taylor & Francis, 2018.
  * [2] D. L. Brown and M. L. Minion, “Performance of under-resolved two-dimensional incompressible flow simulations,” Journal of Computational Physics, vol. 122, no. 1, pp. 165–183, 1995.
  * [3] J. Perot, “An analysis of the fractional step method,” Journal of Computational Physics, vol. 108, pp. 51–58, Sept. 1993.
  * [4] S. Borok, S. Ansumali, and I. V. Karlin, “Kinetically reduced local Navier-Stokes equations for simulation of incompressible viscous flows,” Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, vol. 76, no. 6, pp. 1–9, 2007.
  * [5] A. J. Chorin, “A numerical method for solving incompressible viscous flow problems,” Journal of Computational Physics, vol. 2, pp. 12–26, Aug. 1967\.
  * [6] E. Turkel, “Preconditioned methods for solving the incompressible and low speed compressible equations,” Journal of Computational Physics, vol. 72, pp. 277–298, Oct. 1987.
  * [7] R. Peyret and T. D. Taylor, Computational methods for fluid flow.  Springer Berlin Heidelberg, 1983.
  * [8] C. Merkle, “Time-accurate unsteady incompressible flow algorithms based on artificial compressibility,” 8th Computational Fluid Dynamics Conference, 1987.
  * [9] W. Soh, “Time-marching solution of incompressible Navier-Stokes equations for internal flow,” Journal of Computational Physics, vol. 70, pp. 232–252, May 1987.
  * [10] S. E. Rogers and D. Kwak, “An upwind differencing scheme for the incompressible Navier–Stokes equations,” Applied Numerical Mathematics, vol. 8, no. 1, pp. 43–64, 1991.
  * [11] C. Kiris and D. Kwak, “Aspects of unsteady incompressible flow simulations,” Computers and Fluids, vol. 31, no. 4-7, pp. 627–638, 2002\.
  * [12] R. Nourgaliev, T. Dinh, and T. Theofanous, “A pseudocompressibility method for the numerical simulation of incompressible multifluid flows,” International Journal of Multiphase Flow, vol. 30, pp. 901–937, July 2004.
  * [13] W. W. Kim and S. Menon, “An unsteady incompressible Navier-Stokes solver for large eddy simulation of turbulent flows,” International Journal for Numerical Methods in Fluids, vol. 31, no. 6, pp. 983–1017, 1999.
  * [14] R. Ranjan, M. K. Venkataswamy, and S. Menon, “Dynamic one-equation-based subgrid model for large-eddy simulation of stratified turbulent flows,” Physical Review Fluids, vol. 5, p. 064601, June 2020.
  * [15] J. R. Clausen, “Entropically damped form of artificial compressibility for explicit simulation of incompressible flow,” Physical Review E, vol. 87, no. 1, pp. 1–12, 2013.
  * [16] A. Toutant, “General and exact pressure evolution equation,” Physics Letters, Section A: General, Atomic and Solid State Physics, vol. 381, no. 44, pp. 3739–3742, 2017.
  * [17] A. Toutant, “Numerical simulations of unsteady viscous incompressible flows using general pressure equation,” Journal of Computational Physics, vol. 374, pp. 822–842, 2018.
  * [18] A. Kajzer and J. Pozorski, “Application of the Entropically Damped Artificial Compressibility model to direct numerical simulation of turbulent channel flow,” Computers & Mathematics with Applications, vol. 76, pp. 997–1013, Sept. 2018.
  * [19] W. Trojak, N. Vadlamani, J. Tyacke, F. Witherden, and A. Jameson, “Artificial compressibility approaches in flux reconstruction for incompressible viscous flow simulations,” Computers & Fluids, vol. 247, p. 105634, Oct. 2022\.
  * [20] D. Dupuy, A. Toutant, and F. Bataille, “Analysis of artificial pressure equations in numerical simulations of a turbulent channel flow,” Journal of Computational Physics, vol. 411, p. 109407, 2020.
  * [21] X. Shi and C. A. Lin, “Simulations of wall bounded turbulent flows using general pressure equation,” Flow, Turbulence and Combustion, vol. 105, no. 1, pp. 67–82, 2020.
  * [22] A. Kajzer and J. Pozorski, “A weakly compressible, diffuse-interface model for two-phase flows,” Flow, Turbulence and Combustion, vol. 105, pp. 299–333, Aug. 2020.
  * [23] J.-J. Huang, “Numerical simulation of two-phase incompressible viscous flows using general pressure equation,” 2020.
  * [24] H. Bodhanwalla, D. Raghunathan, and Y. Sudhakar, “A general pressure equation based method for incompressible two-phase flows,” May 2023.  arXiv:2305.17885 [physics].
  * [25] M. Sharma, K. Srikanth, T. Jayachandran, and A. Sameen, “DNS of buoyancy-driven flows using EDAC formulation solved by high-order method,” Computers & Fluids, vol. 265, p. 105997, Oct. 2023.
  * [26] M.-P. Bolduc, R. Ghoreishi, and B. C. Vermeire, “A high-order entropically-damped artificial compressibility approach on moving and deforming domains,” Computers & Fluids, vol. 257, p. 105839, May 2023\.
  * [27] S. Ansumali, I. V. Karlin, and H. C. Öttinger, “Thermodynamic theory of incompressible hydrodynamics,” Physical Review Letters, vol. 94, no. 8, pp. 1–4, 2005.
  * [28] J. Lu and N. A. Adams, “General mechanisms for stabilizing weakly compressible models,” Physical Review E, vol. 107, p. 055306, 2023.
  * [29] J. Ramshaw and V. Mousseau, “Accelerated artificial compressibility method for steady-state incompressible flow calculations,” Computers & Fluids, vol. 18, pp. 361–367, Jan. 1990.
  * [30] P. McHugh and J. Ramshaw, “Damped artificial compressibility iteration scheme for implicit calculations of unsteady incompressible flow,” International journal for numerical methods in fluids, vol. 21, no. 2, pp. 141–153, 1995.
  * [31] J. Ramshaw and V. Mousseau, “Damped artificial compressibility method for steady-state low-speed flow calculations,” Computers & fluids, vol. 20, no. 2, pp. 177–186, 1991.
  * [32] K. Mazaheri and P. L. Roe, “Bulk viscosity damping for accelerating convergence of low mach number euler solvers,” International journal for numerical methods in fluids, vol. 41, no. 6, pp. 633–652, 2003.
  * [33] T. Yasuda, I. Tanno, T. Hashimoto, K. Morinishi, and N. Satofuka, “Artificial compressibility method using bulk viscosity term for an unsteady incompressible flow simulation,” Computers & Fluids, vol. 258, p. 105885, May 2023.
  * [34] P. J. Dellar, “Bulk and shear viscosities in lattice Boltzmann equations,” Physical Review E, vol. 64, p. 031203, Aug. 2001.
  * [35] P. Asinari and I. V. Karlin, “Quasiequilibrium lattice Boltzmann models with tunable bulk viscosity for enhancing stability,” Physical Review E, vol. 81, p. 016702, Jan 2010.
  * [36] J. B. Avalos, M. Antuono, A. Colagrossi, and A. Souto-Iglesias, “Shear-viscosity-independent bulk-viscosity term in smoothed particle hydrodynamics,” Physical Review E, vol. 101, p. 013302, Jan. 2020.
  * [37] P. Sun, C. Pilloton, M. Antuono, and A. Colagrossi, “Inclusion of an acoustic damper term in weakly-compressible SPH models,” Journal of Computational Physics, vol. 483, p. 112056, June 2023.
  * [38] J. Ramshaw, P. O’Rourke, and A. Amsden, “Acoustic damping for explicit calculations of fluid flow at low Mach number,” Tech. Rep. LA-10641-MS, 6100813, University of North Texas Libraries, UNT Digital Library, January 1986\.
  * [39] D. Pan, “A high-order finite volume method solving viscous incompressible flows using general pressure equation,” Numerical Heat Transfer, Part B: Fundamentals, vol. 82, pp. 146–163, Nov. 2022.
  * [40] S. Parameswaran and J. C. Mandal, “A novel Roe solver for incompressible two-phase flow problems,” Journal of Computational Physics, vol. 390, pp. 405–424, 2019.
  * [41] K. Yang and T. Aoki, “Weakly compressible Navier-Stokes solver based on evolving pressure projection method for two-phase flow simulations,” Journal of Computational Physics, vol. 431, p. 110113, Apr. 2021.
  * [42] X. He, G. D. Doolen, and T. Clark, “Comparison of the lattice Boltzmann method and the artificial compressibility method for Navier–Stokes equations,” Journal of Computational Physics, vol. 179, pp. 439–451, July 2002.
  * [43] K. Nagata, N. Ikegaya, and J. Tanimoto, “Consideration of artificial compressibility for explicit computational fluid dynamics simulation,” Journal of Computational Physics, vol. 443, p. 110524, Oct. 2021.
  * [44] J. B. Bell, P. Colella, and H. M. Glaz, “A second-order projection method for the incompressible Navier-Stokes equations,” Journal of Computational Physics, vol. 85, pp. 257–283, Dec. 1989.
  * [45] M. L. Minion and D. L. Brown, “Performance of under-resolved two-dimensional incompressible flow simulations, II,” Journal of Computational Physics, vol. 138, pp. 734–765, Dec. 1997.
  * [46] A. Shah, L. Yuan, and A. Khan, “Upwind compact finite difference scheme for time-accurate solution of the incompressible Navier-Stokes equations,” Applied Mathematics and Computation, vol. 215, no. 9, pp. 3201–3213, 2010.
  * [47] T. Hashimoto, T. Yasuda, I. Tanno, Y. Tanaka, K. Morinishi, and N. Satofuka, “Multi-GPU parallel computation of unsteady incompressible flows using kinetically reduced local Navier–Stokes equations,” Computers & Fluids, vol. 167, pp. 215–220, May 2018.
  * [48] T. De Mulder, “The role of bulk viscosity in stabilized finite element formulations for incompressible flow: A review,” Computer Methods in Applied Mechanics and Engineering, vol. 163, pp. 1–10, Sept. 1998.
  * [49] U. Rasthofer and V. Gravemeier, “Recent developments in variational multiscale methods for large-eddy simulation of turbulent flow,” Archives of Computational Methods in Engineering, vol. 25, pp. 647–690, July 2018.
  * [50] M. Olshanskii, G. Lube, T. Heister, and J. Löwe, “Grad–div stabilization and subgrid pressure models for the incompressible Navier–Stokes equations,” Computer Methods in Applied Mechanics and Engineering, vol. 198, pp. 3975–3988, Nov. 2009.
  * [51] T. G. Ingram and G. A. Edward, “Mechanism of the production of small eddies from large ones,” in Proceedings of Royal Society A. Lond., vol. A158, pp. 499–521, 1937.
  * [52] A. Sharma and V. Eswaran, “Heat and fluid flow across a square cylinder in the two-dimensional laminar flow regime,” Numerical Heat Transfer, Part A: Applications, vol. 45, pp. 247–269, Feb. 2004.
  * [53] A. Sohankar, C. Norberg, and L. Davidson, “Low-Reynolds-number flow around a square cylinder at incidence: study of blockage, onset of vortex shedding and outlet boundary condition,” International Journal for Numerical Methods in Fluids, vol. 26, pp. 39–56, Jan. 1998.
  * [54] S. Sen, S. Mittal, and G. Biswas, “Flow past a square cylinder at low Reynolds numbers,” International Journal for Numerical Methods in Fluids, vol. 67, pp. 1160–1174, Nov. 2011.
  * [55] M. Zhao, L. Cheng, and T. Zhou, “Numerical simulation of vortex-induced vibration of a square cylinder at a low Reynolds number,” Physics of Fluids, vol. 25, p. 023603, Feb. 2013.
  * [56] A. K. Saha and A. Shrivastava, “Suppression of vortex shedding around a square cylinder using blowing,” Sadhana, vol. 40, pp. 769–785, May 2015.
  * [57] S. Taneda and H. Honji, “Unsteady flow past a flat plate normal to the direction of motion,” Journal of the Physical Society of Japan, vol. 30, no. 1, pp. 262–272, 1971.
  * [58] P. Koumoutsakos and D. Shiels, “Simulations of the viscous flow normal to an impulsively started and uniformly accelerated flat plate,” Journal of Fluid Mechanics, vol. 328, pp. 177–227, Dec. 1996.

##  Appendix A Significance of the additional bulk viscosity term

Due to its non-homogeneous nature, the artificial bulk viscosity term expands
to
$\nabla\cdot\left(\boldsymbol{\mathcal{B}}\nabla\cdot\mathbf{u}\right)=\boldsymbol{\mathcal{B}}\nabla\left(\nabla\cdot\mathbf{u}\right)+(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$.
Here, the second term
$(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$ is
zero for a uniform mesh; on a non-uniform mesh, its magnitude depends on the
stretching ratio, which, for all practical computations, is maintained close
to 1. This implies $\nabla\cdot\boldsymbol{\mathcal{B}}\ll 1$; hence, this
term makes only a non-significant contribution. In this section, we examine
this presumption using the two test cases for which we use stretched grids.
Figure 22 shows the $C_{L}$ and $C_{D}$ plots during the periodic phase
obtained from the simulations of flow over a square cylinder, with and without
the additional term. Similarly, in figure 23, we present the evolution of
$C_{D}$ for the case of flow over an impulsively started thin plate. In both
cases, we observe only negligible contribution from
$(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$.
This is also evident from table 6.

Figure 22: Effect of second bulk viscosity term for the flow over a square
cylinder test case. Figure 23: Effect of second bulk viscosity term for the
flow over an impulsively started thin plate test case. Table 6:  Comparison of
simulation data with and without
$(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$ for
the flow over an impulsively started thin plate test case.  Square Cylinder
---
Parameter |  with $(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$ |  without $(\nabla\cdot\boldsymbol{\mathcal{B}})\left(\nabla\cdot\mathbf{u}\right)$ | % Deviation
$C_{D_{avg}}$ | 1.5336 | 1.5344 | 0.052
$C_{L_{rms}}$ | 0.1924 | 0.1917 | 0.36
Impulsively started plate
$C_{D}$ at $t=10$ | 1.678682 | 1.679478 | 0.047

##  Appendix B Effect of pressure diffusion term on damping the acoustic waves

Figure 24: Effect of Prandtl number on the evolution of bubble length
($L_{B}$) and drag coefficient ($C_{D}$) for the flow over an impulsively
started thin plate test case.

We conducted a detailed parametric study on $Pr$ for the flow over impulsively
started thin plate case. The values used are $Pr$=0.1, 0.01 and 0.001, and the
plots of $L_{B}$ and $C_{D}$ with respect to time are presented in figure 24.
From the plots, we conclude that the acoustic waves reduce as $Pr$ reduces.
The most effective damping is observed for the lowest value of $Pr$
($Pr=0.001$). However, the corresponding $L_{B}$ curve deviates from the
reference values. The maximum error for $L_{B}$ at $t=10$ is estimated to be
4.4%. Moreover, considering stability issues, reducing $Pr$ results in a
compulsory need to reduce the time step. For the simulations with $Pr=0.1$ and
$Pr=0.01$, the time-step required is $10^{-5}$ and $10^{-6}$, respectively,
whereas, for $Pr=0.001$, we need to reduce it further to $10^{-7}$. Evidently,
the rate of damping with $Pr=0.001$ is still not sufficient, and further
reducing $Pr$ results in both increased computational cost and reduced
accuracy. To reiterate, any attempt to damp the acoustic waves by adjusting
the numerical parameters in GPE negates the advantage of GPE over other
conventional schemes in terms of computational efficiency.


#  $N-N$, $P_{T}-N$ and $P_{T}-P_{T}$ fluctuations in nucleus-nucleus
collisions at the NA61/SHINE experiment

Evgeny Andronov for the NA61/SHINE Collaboration  Saint Petersburg State
University, ul. Ulyanovskaya 1, 198504, Petrodvorets, Saint Petersburg, Russia
<EMAIL_ADDRESS>
###### Abstract

The NA61/SHINE experiment aims to discover the critical point of strongly
interacting matter and study the properties of the onset of deconfinement. For
these goals a scan of the two dimensional phase diagram ($T-\mu_{B}$) is being
performed at the SPS by measurements of hadron production in proton-proton,
proton-nucleus and nucleus-nucleus interactions as a function of collision
energy. This paper presents preliminary results from Be+Be collisions on
pseudorapidity dependences of transverse momentum and multiplicity
fluctuations expressed in terms of strongly intensive quantities. It is shown
that non-trivial effects evolve from the Poissonian-like fluctuations for
small pseudorapidity intervals with expansion of the acceptance. These
fluctuations are supposed to be sensitive to the existence of the critical
point. The results will be compared to the predictions from the EPOS model.

##  1 Introduction

The NA61/SHINE experiment [1] is a multi-purpose fixed target experiment at
the Super Proton Synchrotron (SPS) of the European Organization for Nuclear
Research (CERN). The strong interactions programme of NA61/SHINE consists of
studies of the onset of deconfinement (OD) [2] in nucleus-nucleus collisions
and search for the critical point (CP) [3] of strongly interacting matter.

NA61/SHINE performs measurements of hadron production in collisions of protons
and various nuclei (p+p, p+Pb, Be+Be, Ar+Sc, Xe+La, Pb+Pb) in a range of beam
momenta (13A \- 150/158A GeV/c). It is expected that there will be a non-
monotonic dependence of fluctuations of a number of observables on energy and
system size in this scan due to the phase transition of strongly interacting
matter and the possible existence of the CP [4]. Some hints of such behaviour
were already observed by the NA49 experiment [5].

##  2 $P_{T}-N$ fluctuation measures

The volume of the system created in ultrarelativistic ion collisions
fluctuates from event to event. In order to suppress contribution from this
”trivial” fluctuations strongly intensive observables are used [6, 7]:

$\displaystyle\Delta[A,B]=\frac{1}{C_{\Delta}}\biggl{[}\langle
B\rangle\omega[A]-\langle A\rangle\omega[B]\biggr{]}$ (1)
$\displaystyle\Sigma[A,B]=\frac{1}{C_{\Sigma}}\biggl{[}\langle
B\rangle\omega[A]+\langle A\rangle\omega[B]-2\bigl{(}\langle AB\rangle-\langle
A\rangle\langle B\rangle\bigr{)}\biggr{]},$ (2)

where $\langle\cdots\rangle$ stands for averaging over all events and
$\omega[X]$ is the scaled variance of any quantity $X$ defined as $(\langle
X^{2}\rangle-{\langle X\rangle}^{2})/\langle X\rangle$. In case of joint
fluctuations of the total transverse momentum $P_{T}$ and multiplicity $N$ one
defines: $A=P_{T}=\sum_{i=1}^{N}p_{T_{i}}$, $B=N$,
$C_{\Delta}=C_{\Sigma}=\langle N\rangle\omega[p_{T}]$, where $\omega[p_{T}]$
is the scaled variance of the inclusive $p_{T}$ spectrum.

The advantage of the strongly intensive quantities is that in the models of
independent sources [8] and in the models of independent particle production
the contribution from volume fluctuations is eliminated, allowing to probe the
genuine CP signals. Moreover, in the models of independent particle production
$\Delta[P_{T},N]=\Sigma[P_{T},N]=1$. These equalities also hold for the ideal
Boltzmann gas both in the grand canonical and in the canonical ensemble
formulations. Another trivial property of these quantities is that they vanish
in case of absence of fluctuations of $A$ and $B$.

In recent measurements by the NA61/SHINE collaboration [9, 10] of the strongly
intensive quantities $\Delta[P_{T},N]$ and $\Sigma[P_{T},N]$ no anomaly
attributable to the CP was observed neither in p+p nor in forward energy
selected Be+Be and Ar+Sc collisions. This observation led to the extension of
the analysis to studies of the pseudorapidity dependence of strongly intensive
observables. By changing the width of the chosen pseudorapidity interval one
can probe different values of baryochemical potential as evident from the
different shapes of proton and antiproton rapidity spectra [11]. Therefore,
the new analysis allows to probe additional regions of the phase diagram of
strongly interacting matter.

In this paper new results were obtained for Be+Be collisions at 150A GeV/c
selected for the smallest 8$\%$ of forward energies. This selection of the
most central events by the forward energy was done using information from the
NA61/SHINE forward calorimeter PSD. This detector consists of 44 square
modules. A specific set of modules is selected in order to maximize the signal
from projectile nucleon spectators and to suppress background from produced
particles. Analysis was performed for the 120k selected events containing
2000k quality selected tracks in total.

Fluctuations were studied in 9 pseudorapidity intervals defined in the
laboratory reference frame using all charged particles produced in strong and
electromagnetic processes from the primary vertex. The forward edge of the
intervals was fixed at 5.2 units of pseudorapidity, with the backward edge
changing from 3 up to 4.6 units. The choice of lower bounds was motivated by
the small azimuthal angle acceptance at smaller pseudorapidities. The upper
bound was introduced in order to suppress possible quasi-diffractive or
electromagnetic effects which become important at larger pseudorapidities.
Transverse momenta of all charged particles were restricted to $0<p_{T}<1.5$
GeV/c. Moreover, the NA61/SHINE acceptance map [12] was applied.

Figure 1 shows preliminary results on the dependence on the width of the
pseudorapidity interval for $\Delta[P_{T},N]$ and $\Sigma[P_{T},N]$ and
comparisons to the EPOS1.99 model [13] predictions. Both quantities change
monotonously for the data in contrast to the EPOS1.99 results for
$\Delta[P_{T},N]$ which show a minimum for intermediate width and lie
significantly above the measurements. EPOS1.99 also fails to describe
$\Delta[P_{T},N]$ for the full NA61/SHINE pseudorapidity acceptance in
inelastic p+p interactions [9]. Another observation is that for smaller
windows these quantities approach unity (independent particle production
limit) as the number of particles in the interval gets small. Moreover, the
inequalities $\Delta[P_{T},N]<1$ and $\Sigma[P_{T},N]\geq 1$, previously
observed for all studied systems at all collision energies [9, 10], also hold
when modifying the width of the pseudorapidity window of the measurement. In
general, no traces of the possible critical point of strongly interacting
matter are visible.

(a)

(b)

Figure 1: Dependence of $\Delta[P_{T},N]$ (a) and $\Sigma[P_{T},N]$ (b) on the
width of the pseudorapidity window. Red squares are preliminary NA61/SHINE
measurements for 0-8$\%$ Be+Be collisions at 150A GeV/c. Blue band represents
the EPOS1.99 model predictions.

The results were not corrected for detector inefficiencies and trigger biases
as simulations have shown that their effect estimated using the GEANT3 package
does not exceed 5$\%$. In contrast the trigger bias modifies results
significantly for p+p collisions [9].

The statistical uncertainties were determined using the sub-sample method. The
systematic uncertainties were estimated to be smaller than 10$\%$ for
$\Delta[P_{T},N]$ and $\Sigma[P_{T},N]$. Further analysis of the systematic
uncertainties will follow.

##  3 $N-N$, $P_{T}-N$ and $P_{T}-P_{T}$ fluctuations in two rapidity
intervals

By analogy with the $P_{T}-N$ case described in the previous section one can
introduce a number of other strongly intensive observables. Joint fluctuations
of two quantities in two separated pseudorapidity intervals are of special
interest because they are closely connected with studies of forward-backward
correlations that have a long history of measurements [14, 15]. Forward-
backward correlations are usually quantified by the correlation coefficient
which is not strongly intensive and is sensitive to the employed centrality
selection procedures in nucleus-nucleus collisions [16]. Use of strongly
intensive quantities can suppress these trivial fluctuations and allows to
study the intrinsic properties of particle emitting sources. Therefore, for
studies of joint fluctuations of multiplicities ($N_{F}$, $N_{B}$) and the sum
of transverse momenta ($P_{TF}$, $P_{TB}$) in two windows of rapidity the
following quantities are defined [17]:

$\Sigma[N_{F},N_{B}]=\frac{1}{\langle N_{B}\rangle+\langle
N_{F}\rangle}\cdot\biggl{[}\langle N_{B}\rangle\omega[N_{F}]+\langle
N_{F}\rangle\omega[N_{B}]-2\bigl{(}\langle N_{F}N_{B}\rangle-\langle
N_{F}\rangle\langle N_{B}\rangle\bigr{)}\biggr{]}$ (3)
$\Sigma[N_{F},P_{TB}]=\frac{1}{(\langle N_{B}\rangle+\langle
N_{F}\rangle)\cdot\overline{p_{T}}_{B}+\langle
N_{F}\rangle\omega[p_{T}]_{B}}\cdot$ $\biggl{[}\langle
P_{TB}\rangle\omega[N_{F}]+\langle N_{F}\rangle\omega[P_{TB}]-2\bigl{(}\langle
N_{F}P_{TB}\rangle-\langle N_{F}\rangle\langle P_{TB}\rangle\bigr{)}\biggr{]}$
(4) $\Sigma[P_{TF},P_{TB}]=\frac{1}{\langle
P_{TB}\rangle\left(\overline{p_{T}}_{F}+\omega[p_{T}]_{F}\right)+\langle
P_{TF}\rangle\left(\overline{p_{T}}_{B}+\omega[p_{T}]_{B}\right)}\cdot$
$\biggl{[}\langle P_{TB}\rangle\omega[P_{TF}]+\langle
P_{TF}\rangle\omega[P_{TB}]-2\bigl{(}\langle P_{TF}P_{TB}\rangle-\langle
P_{TF}\rangle\langle P_{TB}\rangle\bigr{)}\biggr{]}$ (5)

Here, averaging over all tracks in the forward (backward) window is denoted
as:

$\overline{p_{T}}_{F,B}=\frac{1}{N_{tr}}\sum_{tr_{F,B}}p_{T}\ .$ (6)

Moreover, the scaled variance of inclusive $p_{T}$ spectra
$\omega[p_{T}]_{F,B}$ is introduced for both forward and backward windows.

Typically, forward-backward correlations are studied not for the sum of
transverse momenta but for the event mean transverse momentum [18], but in
order to construct a strongly intensive quantity from two observables $A$ and
$B$ both of these quantities have to be extensive [6].

The analysis was performed for 7 pairs of pseudorapidity intervals, with
forward window fixed at (4.7,5.2) and the backward window moving from (3,3.5)
to (4.2,4.7). Figure 2 shows preliminary results for the dependence of
quantities (3)-(5) on the separation between the two windows. It is peculiar
that all three quantities show qualitatively the same behaviour which is quite
well reproduced by the EPOS1.99 model. Moreover, analogous behaviour was
predicted by the quark-gluon string model [19] for p+p interactions at LHC
energies. The PHOBOS collaboration studied joint fluctuations of
multiplicities in two windows quantified by the variance of
$C=\frac{N_{F}-N_{B}}{\sqrt{N_{F}+N_{B}}}$ [20] which has similar properties
like $\Sigma[N_{F},N_{B}]$, although it is not strongly intensive. The
variance of $C$ showed a similar trend with separation between rapidity
windows in Au+Au collisions at $\sqrt{s_{NN}}=200$ GeV.

(a)

(b)

(c)

Figure 2:  Dependence of $\Sigma[N_{F},N_{B}]$ (a), $\Sigma[N_{F},P_{TB}]$ (b)
and $\Sigma[P_{TF},P_{TB}]$ (c) on the separation between the two
pseudorapidity windows. Red squares are preliminary NA61/SHINE measurements
for 0-8$\%$ Be+Be collisions at 150A GeV/c. Blue band represents the EPOS1.99
model predictions.

Similarly to the results in the previous section, no corrections for detector
inefficiencies and trigger biases were applied as simulations have shown that
their effect does not exceed 5$\%$.

The statistical uncertainties were estimated using the sub-sample method, and
the systematic uncertainties were estimated to be smaller than 10$\%$.

##  4 Conclusions

Dependence of the quantities $\Delta[P_{T},N]$ and $\Sigma[P_{T},N]$ on the
width of the pseudorapidity interval were measured by the NA61/SHINE
experiment in forward energy selected Be+Be collisions at 150A GeV/c. No
indications of the critical point of strongly interacting matter were
observed. However, a significant qualitative discrepancy between the
experimental results on $\Delta[P_{T},N]$ and the EPOS1.99 model predictions
is clearly visible.

NA61/SHINE preliminary results for $\Sigma[N_{F},N_{B}]$,
$\Sigma[N_{F},P_{TB}]$ and $\Sigma[P_{TF},P_{TB}]$ on the separation of the
forward and backward pseudorapidity intervals were also presented. It is the
first measurement of this type at SPS energies. Experimental results are in
good agreement with the EPOS1.99 model predictions. The observed trend to
increase with separation is also similar to the predictions of the quark-gluon
string model for LHC energies.

This work was supported by the Russian Science Foundation under grant
17-72-20045.

## References

## References

  * [1] Abgrall N et al. [NA61 Coll.] 2014 JINST 9 06005
  * [2] Gazdzicki M and Gorenstein M 1999, Acta Phys. Polon. B 30 2705
  * [3] Fodor Z and Katz S 2004 J. High Energy Phys. JHEP0404(2004)050
  * [4] Stephanov M A, Rajagopal K and Shuryak E V 1999 Phys. Rev. D 60 114028
  * [5] Grebieszkow K [for the NA49 Coll.] 2009 Nucl. Phys. A 830 547C
  * [6] Gorenstein M and Gazdzicki M 2011 Phys. Rev. C 84 014904
  * [7] Gazdzicki M, Gorenstein M and Mackowiak-Pawlowska M 2013 Phys. Rev. C 88 p 024907
  * [8] Bialas A, Bleszynski M and Czyz W 1976 Nucl. Phys. B 111 461
  * [9] Aduszkiewicz A et al [NA61 Coll.] 2016 Eur. Phys. J. C 76 no. 11 p 635
  * [10] Andronov E [for the NA61/SHINE Coll.] 2017 Acta Phys. Polon. B Proceedings supplement 10 no. 3 p 449
  * [11] Aduszkiewicz A et al [NA61 Coll.] 2017 Eur. Phys. J. C 77 no. 10 p 671
  * [12] Czopowicz T Acceptance map used in paper https://edms.cern.ch/document/1549298/1
  * [13] Pierog T and Werner K 2009 Nucl. Phys. Proc. Suppl. 196 102-105
  * [14] Derrick M et al. 1986 Phys. Rev. D 34 3304
  * [15] Adam J et al. [ALICE Coll.] 2015 J. High Energy Phys. JHEP1505(2015)097
  * [16] Abelev B et al. [STAR Coll.] 2009 Phys. Rev. Lett. 103 172301
  * [17] Andronov E 2015 Theor. Math. Phys. 185 no.1 p 1383
  * [18] Braun M A, Kolevatov R S, Pajares C and Vechernin V V 2004 Eur. Phys. J. C 32 535
  * [19] Vechernin V and Andronov E 2017 Strongly intensive observable for multiplicities in forward and backward windows in string model (XII Workshop on Particle Correlations and Femtoscopy, Amsterdam, 12-16 June 2017) https://indico.cern.ch/event/539093/contributions/2568069/
  * [20] Back B et al. [PHOBOS Coll.] 2006 Phys. Rev. C 74 011901