Non-homogeneous anisotropic bulk viscosity for
acoustic wave attenuation 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. Recently, AbdulGafoor et al. [34] compared the performance of EDAC and bulk viscosity-enabled classical ACM on static and dynamically distorted grids. They used an isotropic bulk viscosity and observed that the choice of the grid influences the flow field more than the artificial/weakly compressible method used.

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 [35] 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 [36] 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. [37] 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. [38] 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.

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. [38] 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 bluff body. The above advantages would not have been realized with the isotropic and uniform bulk viscosity used throughout the literature. An earlier technical report [39] 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.

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.

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

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, 40, 24].

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

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; The dimensionless bulk viscosity is defined as $\boldsymbol{\mathcal{B}}=(\boldsymbol{\mathcal{B}}^{*}/\rho UL)$, where $\boldsymbol{\mathcal{B}}^{*}$ is the respective dimensional quantity. In the present work, we propose the following form of the anisotropic bulk viscosity tensor.

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

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}}=g(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 [37].

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))

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

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

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

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}}{(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}$ $=\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}}_{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}}_{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.

We apply the three-stage Strong Stability Preserving Runge-Kutta [41] 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

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

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, 42]. 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.

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,

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.

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$.

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. Based on the non-dimensional bulk viscosity $\mathcal{B}$ defined earlier, we can interpret the parameter $\mathcal{B}Re$ to provide the ratio of the bulk viscosity to the dynamic viscosity. For this simulation, $\mathcal{B}Re=78.125$. 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.

To ensure accuracy is maintained with the addition of artificial bulk viscosity, we compared velocity profiles from both solvers with those of Ghia et al. [44]. It can be seen in figure 4 that the plots from GPE and GPE+BV solvers overlap, and they match well with the reference results. This indicates the velocity profiles are intact, and $\mathcal{B}$ only modifies the pressure field.

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 5. 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 [45] and lattice Boltzmann methods [45, 46]. Similar pressure plots presented by [46] 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.

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 [47, 2, 48]. The initial conditions are as follows

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 [49, 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. A direct comparison between the vorticity contours at $t=1$ from both the solvers is presented in figure 6. For the GPE+BV solver, we use, $\mathcal{B}=50\Delta x$ (corresponds to $\mathcal{B}Re=976.5625$). We observe that the contours from both simulations overlap with each other except for the contour $\omega_{z}=0$, for which the GPE+BV solver produces a more accurate result. Moreover, the contours in figure 6 qualitatively resemble the vorticity plots, without any spurious vortices, reported in literature [15, 50, 17]. We stress that, the filled contours from both simulations are indistinguishable. We purposefully included the line contours to highlight such negligible differences.

The significant advantage of using GPE+BV over GPE can be seen from figure 7, which illustrates the mass conservation error, quantified by the divergence of the velocity vector, 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.

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,

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$

The Taylor Green vortices (TGV) [54] 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 formula [17], and the initial conditions are obtained from the same by setting $t=0$.

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 9, which shows the variation of $L_{\infty}$-norm of error with the number of grid points ($N$). For the GPE+BV solver, the convergence study ranges from $\mathcal{B}Re=19.53125$ for $256\times 256$ grid to $\mathcal{B}Re=625.0$ for $8\times 8$ grid. 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.

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 10 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 [55]. 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$.

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 11. 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 11. 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.

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 12. Compared to the results shown in figure 11 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.

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. Also, this corresponds to $\mathcal{B}_{\max}Re=937.50$ and $\mathcal{B}_{\min}Re=31.2525$. 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 13.

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 14. The periodic state is achieved at around $t=400$, which is a significant improvement when compared to the GPE.

In figure 15, 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 16. 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.

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.

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 [55]. Thus, with the proposed approach, we achieved sufficient acoustic wave damping without affecting the solver’s accuracy.

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 16. 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 17 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.

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 18. 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 18. 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.

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 19. For all the variants, we use the constant $\lambda=79$, which indicates $\mathcal{B}_{\max}Re=995.40$ and $\mathcal{B}_{\min}Re=24.885$. The following observations can be made from the figure.

• 1.

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

• 2.

  Both the homogeneous isotropic (figure 19(b)) and non-homogeneous isotropic(figure 19(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 19(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 16.

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

From figure 19, 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.