The Effect of Iterative Procedures on the Robustness and Fidelity of Augmented Lagrangian SPH

Abstract

The Augmented Lagrangian Smoothed Particle Hydrodynamics (ALSPH) method is a novel incompressible Smoothed Particle Hydrodynamics (SPH) approach that solves Navier–Stokes equations by an iterative augmented Lagrangian scheme through enforcing the divergence-free coupling of velocity and pressure fields. This study aims to systematically investigate the time step size and the number of inner iteration parameters to boost the performance of the ALSPH method. Additionally, the effects of computing spatial derivatives with two alternative schemes on the accuracy of numerical results are also scrutinized. Namely, the first scheme computes spatial derivatives on the updated particle positions at each iteration, whereas the second one employs the updated pressure and velocity fields on the initial particle positions to compute the gradients and divergences throughout the iterations. These two schemes are implemented to the solution of a flow over a circular cylinder at Reynolds numbers of 200 in two dimensions. Initially, simulations are performed in order to determine the optimum time step sizes by utilizing a maximum number of five iterations per time step. Subsequently, the optimum number of inner iterations is investigated by employing the predetermined optimum time step size under the same flow conditions. Finally, the schemes are tested on the same flow problem with different Reynolds numbers using the best performing combination of the aforementioned parameters. It is observed that the ALSPH method can enable one to increase the time step size without deteriorating the numerical accuracy as a consequence of imposing larger ALSPH penalty terms in larger time step sizes, which, overall, leads to improved computational efficiency. When considering the hydrodynamic flow characteristics, it can be stated that two spatial derivative schemes perform very similarly. However, the results indicate that the derivative operation with the updated particle positions produces slightly lower velocity divergence magnitudes at larger time step sizes.

1. Introduction

Smoothed particle hydrodynamics (SPH) method is a mesh free, particle based, Lagrangian computational method that was introduced simultaneously by Lucy [1] and Gingold and Monaghan [2] for solving astrophysics problems. Following its introduction to hydrodynamics studies by Monaghan [3], the SPH method has become a fast growing computational fluid dynamics (CFD) approach. The initial SPH formulation employed a weakly compressible velocity–pressure coupling scheme, in which the pressure field is numerically computed from the density variation through an equation of state. Although this explicit Weakly Compressible SPH (WCSPH) scheme has clear advantages in tackling various types of violent free-surface hydrodynamics problems, which have large and non-linear deformations [4,5,6], it always requires additional numerical treatments to prevent the occurrence of noisy and oscillatory pressure fields [4,7,8].

Therefore, for dealing with incompressible fluid flow accurately, a fully incompressible SPH (ISPH) approach [9,10] is proposed through being inspired from the algorithm originally developed for a mesh based method [11]. In ISPH, the velocity-pressure coupling scheme is based on the projection method, wherein a computationally expensive implicit solution of pressure Poisson equation is required [12]. Many challenging incompressible fluid mechanics problems were successfully simulated with the conventional ISPH method, including free surface flows [13,14], multiphase flows [15,16], electrohydrodynamics [17,18], among others. The ISPH methodology has also been implemented with hybrid Lagrangian–Eulerian discretization approaches [19] and Eulerian SPH [20] in order to achieve a more consistent and stable numerical framework for the SPH method in general [21]. Moreover, in order to reduce the computational costs, explicit incompressible SPH formulations are also present in the literature [22,23,24], which introduce simplifying assumptions for solving the pressure Poisson equation in an explicit fashion. Augmented Lagrangian SPH (ALSPH) is the most recent and novel explicit incompressible SPH method that has been recently introduced and further improved by the authors [25,26].

In general framework, Augmented Lagrangian method is an optimization method [27], which utilizes Lagrange multipliers, together with penalty functions, in order to minimize (or maximize) an objective function. With this combination, constrained optimization problems are solved as unconstrained sets of equations in a manner that the penalty terms are modified at each iteration with respect to the Lagrangian multipliers defined for the mathematical model. Inspired from this mathematical optimization perspective, the ALSPH method utilizes suitable augmented Lagrangian penalty terms on the pressure and velocity fields decoupled by the projection principle [11] and it performs iterations to minimize velocity divergence without solving implicit Poisson equation for pressure.

Fortin and Glowinski pioneered the idea of solving Navier–Stokes equations as an optimization problem with the help of the Augmented Lagrangian method [28]. Fortin and Glowinski [28] introduced augmented Lagrangian algorithms for solving Navier–Stokes equations and other engineering problems. Followingly, Desai and Ito [29] applied the method for optimal control problems that are governed by steady state Navier-Stokes equations. The method is employed on solving 2D lid-driven cavity and channel flow with sudden expansion problems utilizing a grid based discretization approach. Again, with a grid based approach, Vincent and Caltagirone [30] solved 2D and 3D unsteady fluid flow problems with an incompressible projective augmented Lagrangian application. Vincent et al. [31] applied an adaptive augmented Lagrangian method on 3D direct numerical simulation (DNS) with a multiphase, volume of fluid (VOF) framework, where the augmented Lagrangian is implemented in the predictor step of the projection approach. The adaptivity is provided by evaluating the augmented Lagrangian penalty term locally over a conditional penalty coefficient that is based on the dominant flow characteristics. There are also meshless applications of the approach in solid mechanics field for the solution of material and crack discontinuity problems [32] and on elastic solids [33] to enforce boundary conditions.

To this end, the ALSPH method is the first meshless implementation of the augmented Lagrangian approach to the solution of the Navier–Stokes equations. The augmented Lagrangian penalty term in the ALSPH method is designed to mimic the bulk viscosity term of the Navier–Stokes equations, which should normally be equal to zero for an incompressible flow problem. In this regard, the ALSPH method tries to ensure that this term be zero through its iterative scheme. The penalty term is modified (or augmented) at each iteration while the bulk viscosity term diminishes. Fatehi et al. [25] performed simulations with the ALSPH method on 2D flow over the backward facing step and 2D pressure jump problems, when comparing the results with the weakly compressible SPH (WCSPH) method. In our recent study [26], along with algorithm enhancements and a simple adaptive scheme for the penalty term estimation, we have further investigated the performance of the ALSPH method through solving a challenging incompressible flow problem, namely, 2D flow past a circular cylinder under low to moderate Reynolds numbers. The authors [25,26] demonstrated the advantage of ALSPH over WCSPH in terms of rendering smoother pressure fields and smaller velocity divergences at the expense of a larger computational cost. The findings of our study [26] have further indicated the accuracy of the ALSPH method, especially in higher Reynolds numbers, as compared to the WCSPH method at same particle resolutions and time step sizes. In our previous study [26], the effects of the time step size or the maximum number of iterations per time step were not investigated. When considering higher particle resolution requirements of the WCSPH method at higher Reynolds numbers [34,35,36], our results have also pointed out the advantage of the ALSPH method over WCSPH in terms of computational costs, even without optimizing these temporal and iterative discretization parameters [26].

Given the fact that the ALSPH is a relatively new approach to enforcing the incompressibility condition in particle based simulations, it is important to scrutinize its all aspects to be able to make it a versatile pressure-velocity coupling algorithm. To this end, it is essential to understand the effect of various parameters (i.e., time step, number of iterations and locations of derivative operations) on the accuracy of the ALSPH results. Hence, the main focus of the current study is to determine the optimum time step size, and a practical iteration procedure for the ALSPH method, which involves the number of inner iteration for each time step and the identification of the most effective particle location for computing spatial derivatives during each inner iteration. With n representing temporal discretization in the scheme, the first approach computes the spatial derivatives at $n+1$ by utilizing the displaced particle position results of the previous iteration [25]. In the second approach, the derivatives are always computed with respect to the initial particle positions at n [26]. The velocity and pressure values are approximated with the same SPH procedures in both of the schemes.

Simulations are performed on the benchmark problem of 2D incompressible flow past a circular cylinder using the in-house C++ code in order to examine the numerical performances of both ALSPH schemes and to obtain the optimum time step and number of inner iteration values. The results are compared in terms of the dynamic properties (force and pressure coefficients), the wake kinematics (Strouhal number analysis), and the velocity divergence as a measure of error in enforcing incompressibility condition in the flow domain. Initially, a series of simulations employing the two derivative estimation schemes are performed to determine an optimum dimensionless time step size by limiting the maximum number of iterations per time step to five. Subsequently, utilizing the best applicable dimensionless time step size, the schemes are tested for the iteration number of two, five, and twenty. Finally, the two schemes are examined under Reynolds numbers varying from 100 to 300 utilizing the optimum combination of the aforementioned parameters. Ultimately, the novelty of this paper lies in the detailed and systematic investigations of the above stated parameters, so that the ALSPH method can be reliably used to model incompressible flow problems. The key findings of this research are elaborated within the numerical results section and concisely stated in the conclusion.

2. Governing Equations and ALSPH Methodology

Continuity and linear momentum equations which govern the fluid motion are written as;

$$\frac{D\rho }{Dt}=-\rho \nabla ·\mathbf{u},$$

(1)

$$\rho \frac{D\mathbf{u}}{Dt}=-\nabla p+\nabla ·\mathbf{\tau }+\rho \mathbf{g},$$

(2)

where $D/Dt$ is the material time derivative, t is the time, $\rho$ is the density, $\mathbf{u}$ is the velocity vector, p is the thermodynamic pressure, $\mathbf{\tau }$ is the viscous stress tensor, and $\mathbf{g}$ is the gravitational acceleration vector. The viscous force expressed as the divergence of the viscous stress tensor $\mathbf{\tau }$ in Equation (2) can be written for a Newtonian fluid, as

$$\nabla ·\mathbf{\tau }=\mu {\nabla }^2\mathbf{u}+\kappa \nabla \left(\nabla ·\mathbf{u}\right),$$

(3)

where $\mu$ and $\kappa$ represent the dynamic and bulk viscosity coefficients, respectively. With an incompressible flow assumption, the divergence of velocity becomes zero and the bulk viscosity (volume viscosity) vanishes. Accordingly, the conservation of linear momentum for incompressible flow is written as

$$\rho \frac{D\mathbf{u}}{Dt}=-\nabla p+\mu {\nabla }^2\mathbf{u}+\rho \mathbf{g}.$$

(4)

2.1. Derivation of the ALSPH Formulation

In the ALSPH method, the velocity and pressure fields are decoupled via the projection approach [11]. The Helmholtz-Hodge decomposition theorem states that any arbitrary sufficiently smooth vector field ${\mathbf{u}}^{*}$ can be written as the sum of a divergence-free vector field and the gradient of a scalar field, such that

$${\mathbf{u}}^{*}=\mathbf{u}+\frac{\Delta t}{\rho }\nabla p.$$

(5)

Differing from the conventional projective ISPH schemes, ALSPH utilizes the orthogonal projection operator on the momentum balance equation in Equation (2), rather than the one in Equation (4), such that

$$\frac{{\mathbf{u}}^{*}-\mathbf{u}}{\Delta t}=\nu {\nabla }^2\mathbf{u}+\kappa \nabla (\nabla ·\mathbf{u})+\mathbf{g},$$

(6)

where $\nu$ is the kinematic viscosity that is defined as $\nu =\mu /\rho$. As a result of the projection operation, the pressure gradient term presented in Equation (2) vanishes since it is perpendicular to the divergence-free subspace. Hence, the predicted velocity field ${\mathbf{u}}^{*}$ is an intermediate velocity field without the contribution of pressure forces. The aim of this decomposition approach is to obtain the velocity field $\mathbf{u}$ that ensures the divergence-free condition. Therefore, upon taking divergence of Equation (5), with the assumption of $\nabla ·\mathbf{u}=0$, the pressure Poisson equation is obtained as

$${\nabla }^{2}p=\frac{\rho }{\Delta t}\nabla ·{\mathbf{u}}^{*}.$$

(7)

The augmented Lagrangian method is employed onto Equation (6) by replacing the bulk viscosity coefficient $\kappa$ with an augmented Lagrangian penalty coefficient $r_{AL}$, transforming the volume viscosity force into the augmented Lagrangian penalty term that diminishes through iterations as the gradient of velocity divergence declines. With the superscripts n and m representing the time step and iterations, respectively, initializing the iterations with ${\mathbf{u}}^{n+1,m}={\mathbf{u}}^n$ and $p^{n+1,m}=p^n$, the ALSPH algorithm predicts the intermediate velocity ${\mathbf{u}}^{*}$ as

$${\mathbf{u}}^{*}-{\mathbf{u}}^n=r_{AL}\Delta t{〈\nabla \left(\nabla ·{\mathbf{u}}^{n+1,m}\right)〉}^{n+1,m}+\Delta t\nu {〈{\nabla }^2{\mathbf{u}}^{n+1,m}〉}^{n+1,m}+\Delta t\mathbf{g}.$$

(8)

Here, $r_{AL}$ is the penalty coefficient that is defined as $r_{AL}=C_{BV}{u}_{max}^2\Delta t$, while $u_{max}$ being the maximum velocity magnitude in the domain and $C_{BV}$ is a constant with a value between 1 and 100 [25]. In this study, $C_{BV}$ is taken as 100. The value of $r_{AL}$ is only updated at the beginning of time step and it stays constant through iterations. The angle bracket “〈〉” in Equation (8) instructs that the derivative operation is to be performed based on the particle positions at the time n and iteration m frames that are designated by its superscripts.

ALSPH utilizes an equation of state approach to fulfill the divergence-free condition in order to avoid an implicit solution of Equation (7). Approximating $\rho =c^{-2}p$ in Equation (1), a velocity divergence is obtained, as

$$\nabla ·\mathbf{u}=-\frac{1}{c^2\rho }\frac{Dp}{Dt},$$

(9)

where c is the speed of sound parameter. Taking divergence of Equation (5) and replacing $\nabla ·\mathbf{u}$ with Equation (9), prospective pressure field is obtained in the iterative form as;

$$p^{n+1,m+1}-p^{n+1,m}=-r_{AL}\rho \left({〈\nabla ·{\mathbf{u}}^{*}〉}^{n+1,m}-\frac{\Delta t}{\rho }{〈\nabla ·\nabla p^{n+1,m}〉}^{n+1,m}\right).$$

(10)

In obtaining Equation (10), one can clearly see that the $r_{AL}$ value replaces $c^2\Delta t$. Upon setting $c^2\Delta t$ equal to the above definition of $r_{AL}=C_{BV}{u}_{max}^2\Delta t$, it can be shown that $C_{BV}={(c/u_{max})}^2$. For any flow to be treated incompressible, the ratio of the speed of the flow to the speed of sound should be smaller than 0.3. To be on the safe side, the ratio of 0.1 is utilized, which leads to $C_{BV}$ of 100.

After obtaining new pressures, the velocity field is updated utilizing Equation (5), as follows:

$${\mathbf{u}}^{n+1,m+1}={\mathbf{u}}^{*}-\frac{\Delta t}{\rho }{〈\nabla p^{n+1,m+1}〉}^{n+1,m}.$$

(11)

Respectively, new positions of particles are obtained by

$${\mathbf{r}}_{\mathbf{i}}^{n+1,m+1}={\mathbf{r}}_{\mathbf{i}}^n+\frac{1}{2}\left({\mathbf{u}}_{\mathbf{i}}^n+{\mathbf{u}}_{\mathbf{i}}^{n+1,m+1}\right)\Delta t,$$

(12)

where ${\mathbf{r}}_{\mathbf{i}}$ represents the position vector of the particle $\mathbf{i}$.

The convergence criteria is imposed either ensuring a satisfactory velocity divergence value that is defined by the parameter $ϵ$;

$${\left|{〈\nabla ·{\mathbf{u}}^{n+1,m+1}〉}^{n+1,m+1}\right|}_{max}\le ϵ,$$

(13)

or reaching a predefined maximum number of iterations per time step. In this study, $ϵ$ is taken as ${10}^{-3}\left[s^{-1}\right]$. Section 2.4 provides the iterative algorithm for a given time step for the ALSPH method.

2.2. Spatial Derivatives of Field Functions

The value of any field function in the computational domain can be approximated by utilizing SPH interpolation, as follows [37];

$$f_{\mathbf{i}}^s=\sum _{\mathbf{j}=1}^N{V}_{\mathbf{j}}{f}_{\mathbf{j}}^s{W}_{\mathbf{ij}}.$$

(14)

Here, $f_{\mathbf{i}}^s$ represents the value of any vector valued or scalar function at the spatial coordinates of particle $\mathbf{i}$. Subscript $\mathbf{j}$ designates the variables associated with the neighbor particles of $\mathbf{i}$, which reside within a support domain that is limited by the radius h. Hereby, the neighbors constitute a total number N that varies for each particle $\mathbf{i}$. $W_{\mathbf{ij}}$, or in full form $W_{\mathbf{ij}}\left({\mathbf{r}}_{\mathbf{ij}},h\right)$ is a smoothing kernel/an interpolation weight factor, which is a function of relative particle positions ${\mathbf{r}}_{\mathbf{ij}}={\mathbf{r}}_{\mathbf{i}}-{\mathbf{r}}_{\mathbf{j}}$ and the smoothing length h. As another interpolation factor, $V_{\mathbf{j}}$ is the volume of each neighbor particle computed by the inverse of the summed kernel function as $V_{\mathbf{i}}=1/{\sum }_{\mathbf{j}=1}^N{W}_{\mathbf{ij}}$. This study employs the quintic kernel function used in the reference [38].

In the ALSPH method, flow domain is discretized with particles. Accordingly, the gradient and the Laplacian terms in the ALSPH formulation are evaluated by a corrected SPH approach [39], as follows;

$$\frac{\partial f_{\mathbf{i}}^s}{\partial x_{\mathbf{i}}^k}{\alpha }_{\mathbf{i}}^{kl}=\sum _{\mathbf{j}=1}^N{V}_{\mathbf{j}}\left(f_{\mathbf{j}}^s-f_{\mathbf{i}}^s\right)\frac{\partial W_{\mathbf{ij}}}{\partial x_{\mathbf{i}}^l},$$

(15)

$$\frac{\partial }{\partial x_{\mathbf{i}}^k}\left(\frac{\partial f_{\mathbf{i}}^s}{\partial x_{\mathbf{i}}^k}\right){\alpha }_{\mathbf{i}}^{sl}=8\sum _{\mathbf{j}=1}^N{V}_{\mathbf{j}}\left(f_{\mathbf{i}}^s-f_{\mathbf{j}}^s\right)\frac{r_{\mathbf{ij}}^s}{r_{\mathbf{ij}}^2}\frac{\partial W_{\mathbf{ij}}}{\partial x_{\mathbf{i}}^l},$$

(16)

$${\alpha }_{\mathbf{i}}^{sl}=\sum _{\mathbf{j}=1}^N{r}_{\mathbf{ji}}^s{V}_{\mathbf{j}}\frac{\partial W_{\mathbf{ij}}}{\partial x_{\mathbf{i}}^l},$$

(17)

where ${\alpha }_{\mathbf{i}}^{sl}$ is a second rank correction tensor.

2.3. Artificial Particle Displacement

Before moving to the next time step, the particle positions are shifted by the artificial particle displacement method [8,40] in order to prevent non-physical particle clustering effect. The corrected particle position ${\widehat{\mathbf{r}}}_{\mathbf{i}}$ is calculated as:

$${\widehat{\mathbf{r}}}_{\mathbf{i}}={\mathbf{r}}_{\mathbf{i}}+\delta {\mathbf{r}}_{\mathbf{i}},\delta {\mathbf{r}}_{\mathbf{i}}=\sum _{\mathbf{j}=1}^N\frac{{\mathbf{r}}_{\mathbf{ij}}}{r_{\mathbf{ij}}^3}{r}_0^2{u}_v\Delta t,u_v=\frac{∥{\sum }_{\mathbf{j}=1}^N\left({\mathbf{u}}_{\mathbf{i}}-{\mathbf{u}}_{\mathbf{j}}\right)W_{\mathbf{ij}}∥}{{\sum }_{\mathbf{j}=1}^N{W}_{\mathbf{ij}}},$$

(18)

where $r_{\mathbf{ij}}=∥{\mathbf{r}}_{\mathbf{ij}}∥$, $r_0={\sum }_{\mathbf{j}=1}^N{r}_{\mathbf{ij}}/N$ and $\delta {\mathbf{r}}_{\mathbf{i}}$ is the position displacement vector. Accordingly, the velocity and pressure of the particle are also corrected, as [25]:

$${\widehat{\mathbf{u}}}_{\mathbf{i}}={\mathbf{u}}_{\mathbf{i}}+{\widehat{\mathbf{r}}}_{\mathbf{i}}·{〈\nabla ·\mathbf{u}〉}_{\mathbf{i}},{\widehat{p}}_{\mathbf{i}}=p_{\mathbf{i}}+{\widehat{\mathbf{r}}}_{\mathbf{i}}·{〈\nabla p〉}_{\mathbf{i}}.$$

(19)

2.4. Derivative Estimation Schemes

The original ALSPH algorithm that was developed by Fatehi et al. [25] is explained at the beginning of this section. In our previous study [26], we proposed dropping the computationally expensive, repeating neighbor searching operation during the iterations. Furthermore, the algorithm on the estimation of the spatial derivatives was modified. The aim of the present study is to compare two derivative estimation schemes, performing a single neighbor searching operation in each time step.

The first scheme can be considered to be a modified version (without a repetitive neighbor search) of the scheme that was introduced by Fatehi et al. [25], utilizing the updated particle positions as ${\langle \rangle }^{n+1,m}$. The second one is the same scheme employed in our previous study [26], which uses the particle positions at the beginning of the time step as ${\langle \rangle }^n$. The aforementioned schemes will respectively be referred to as Algorithms 1 and 2 hereafter. The superscripts belonging to the angle brackets in Equations (8), (10) and (11), corresponding to Algorithm 1, are changed from ${\langle \rangle }^{n+1,m}$ to the form ${\langle \rangle }^n$, hence constructing Algorithm 2, as given below as pseudo codes.

Algorithm 1 The pseudo code for the algorithm that computes spatial derivatives on updated particle positions

1: Time step initialization with ${\mathbf{r}}^n$, ${\mathbf{u}}^n$, $p^n$ 2: Generate ghost particles for wall boundaries (Figure 1 and Figure 2) 3: Neighbor particle search 4: Initialize iterations with: $m=0$, ${\mathbf{r}}^{n+1,m}={\mathbf{r}}^n$, ${\mathbf{u}}^{n+1,m}={\mathbf{u}}^n$, $p^{n+1,m}=p^n$ 5: while$m\le$ max iterations do 6:     compute ${\mathbf{u}}^{*}$ (Equation (8))       ▷ compute derivative terms using ${\mathbf{r}}^{n+1,m}$ 7:     compute $p^{n+1,m+1}$ (Equation (10))  ▷ compute derivative terms using ${\mathbf{r}}^{n+1,m}$ 8:     update pressures of cylinder particles (Equation (20)) 9:     compute ${\mathbf{u}}^{n+1,m+1}$ (Equation (11))  ▷ compute derivative terms using ${\mathbf{r}}^{n+1,m}$ 10:     compute ${\mathbf{r}}^{n+1,m+1}$ (Equation (12)) 11:     check for convergence (Equation (13)) 12:     if converged then 13:         apply artificial particle displacement (Equations (18) and (19)) 14:         finalize time step 15:     else 16:         $m=m+1$

Algorithm 2 The pseudo code for the algorithm that computes spatial derivatives on initial particle positions

1: Time step initialization with ${\mathbf{r}}^n$, ${\mathbf{u}}^n$, $p^n$ 2: Generate ghost particles for wall boundaries (Figure 1 and Figure 2) 3: Neighbor particle search 4: Initialize iterations with: $m=0$, ${\mathbf{r}}^{n+1,m}={\mathbf{r}}^n$, ${\mathbf{u}}^{n+1,m}={\mathbf{u}}^n$, $p^{n+1,m}=p^n$ 5: while$m\le$ max iterations do 6:     compute ${\mathbf{u}}^{*}$ (Equation (8))       ▷ compute derivative terms using ${\mathbf{r}}^n$ 7:     compute $p^{n+1,m+1}$ (Equation (10))  ▷ compute derivative terms using ${\mathbf{r}}^n$ 8:     update pressures of cylinder particles (Equation (20)) 9:     compute ${\mathbf{u}}^{n+1,m+1}$ (Equation (11))  ▷ compute derivative terms using ${\mathbf{r}}^n$ 10:     compute ${\mathbf{r}}^{n+1,m+1}$ (Equation (12)) 11:     check for convergence (Equation (13)) 12:     if converged then 13:         apply artificial particle displacement (Equations (18) and (19)) 14:         finalize time step 15:     else 16:         $m=m+1$

Figure 1. Computational domain.

Figure 1. Computational domain.

Figure 2. Free slip ghost particles for channel boundaries [26].

Figure 2. Free slip ghost particles for channel boundaries [26].

3. Problem Definition

In this study, 2D incompressible channel flow past circular cylinder is simulated to determine the best possible ALSPH parameters in terms of accuracy and computational costs. Channel geometry is given in Figure 1, where D is the cylinder diameter and $L=14D$, $H=13D$, and $L_U=4D$.

Reynolds number is defined as $Re=UD/\nu$ with uniform inlet velocity U. The initial particle discretization is realized by uniform orthogonal fluid particles with equal particle distances of $\Delta x$ in both x and y axes. Channel walls are represented by ghost fluid particles with a free-slip boundary condition, as illustrated in Figure 2 [26]. Two buffer zones with lengths of 6 h are defined along the inlet and outlet regions to ensure mass conservation within the channel (Figure 1). Particles of the inlet buffer zone enters into the flow domain with the constant velocity U. Fluid particles reaching the outlet buffer zone are enforced to preserve the normal component of their velocity to the boundary until they travel beyond 6 h thick outlet buffer zone. Fixed solid cylinder particles are radially distributed with respect to the same $\Delta x$ scale. An additional pressure term is imposed upon the cylinder particles in order to ensure no-slip and Neumann boundary conditions for flow around the solid obstacle;

$$\frac{\partial p}{\partial n}=\rho \left[\frac{D{\mathbf{u}}_s}{Dt}·\mathbf{n}\right],p_s=p+\frac{\partial p}{\partial n}{r}_s,$$

(20)

where ${\mathbf{u}}_s$ is a pseudo-velocity for the solid particles computed using Equation (8), $\mathbf{n}$ is the outwards facing surface normal, $p_s$ is the modified pressure, $\partial p/\partial n=\nabla p·\mathbf{n}$ and $r_s$ is the pseudo-displacement along the normal direction.

The instantaneous values of drag coefficient $C_D\left(t\right)=2F_D\left(t\right)/\left(\rho U^{2}A\right)$ and lift coefficient $C_L\left(t\right)=2F_L\left(t\right)/\left(\rho U^{2}A\right)$ for the cylinder are computed by integrating the forces acting on the solid particles of the cylinder as $F^k\left(t\right)={\sum }_{\mathbf{i}=1}^{N_{cyl}}{V}_{\mathbf{i}}{\rho }_{\mathbf{i}}{a}_{\mathbf{i}}^k\left(t\right)$, where $a_{\mathbf{i}}^k\left(t\right)=D{u}_{\mathbf{i}}^k\left(t\right)/Dt$ and $N_{cyl}$ is the number of cylinder particles. Similarly, the pressure coefficient is computed over the cylinder surface as $C_P(t,\theta )=2\left(p(t,\theta )-p_0\left(t\right)\right)/\left(\rho U^2\right)$, where $\theta$ is the angular coordinate with respect to the center of the cylinder (Figure 1) and $p_0\left(t\right)$ is the instantaneous mean inlet pressure. Because the flow is oscillatory, the mean values of the drag and pressure coefficients as well as the root mean square of lift coefficient, respectively, $C_D^{{}^{\prime }}$, $C_P^{{}^{\prime }}$, and $C_L^{{}^{\prime }}$ are computed over the averaging period T in accordance with the following relations;

$$C_D^{{}^{\prime }}=\frac{1}{T}{\int }_t^{t+T}{C}_D\left(t\right)dt,$$

(21)

$$C_P^{{}^{\prime }}\left(\theta \right)=\frac{1}{T}{\int }_t^{t+T}{C}_P(t,\theta )dt,$$

(22)

$$C_L^{{}^{\prime }}={\left[\frac{1}{T}{\int }_t^{t+T}{C}_L^2\left(t\right)dt\right]}^{\frac{1}{2}}.$$

(23)

The mean velocity divergence in the domain is computed as ${\int }_t^{t+T}\left({\sum }_i^{N_f}\left|〈\nabla ·{\mathbf{u}}_{\mathbf{i}}\left(t\right)〉\right|\right)dt$, where $N_f$ is the total number of fluid particles.

4. Numerical Results

There are three main investigations in the present study. The first one is the comparison of the spatial derivative Algorithms 1 and 2. To this end, all other investigations are realized for both Algorithms 1 and 2. Two different particle resolutions, namely, $D/\Delta x=20$ and $D/\Delta x=30$ are utilized for all relevant simulations in order to understand the effect of spatial resolution on the numerical results. The second one is to determine the maximum applicable time step size, written in a dimensionless form as $\Delta tU/\Delta x$. Accordingly, time step size is incrementally tested. Finally, after determining an optimum dimensionless time step size for the problem at hand with five iterations per time step, the effect of the number of iterations on the accuracy of the results are examined by utilizing the iteration number of two and 20 per time step. All three investigations to find the optimum parameters are realized at $Re=200$. Ultimately, having determined the best performing time step size and the number of iterations, both 1 and 2 algorithms are tested under $Re=100-300$ with $D/\Delta x=30$ to examine how these parameters respond to the Reynolds number of the flow. It should be emphasized that the algorithms are equivalent in terms of computational cost. In the figures of the following subsections, the results of Algorithms 1 and 2 will be designated as (1) and (2), respectively.

4.1. Time Step Size

The mean force coefficients and vortex shedding characteristics are investigated and presented in Figure 3 and Figure 4, respectively, in order to capture both dynamics and kinematics of the problem. It must be restated that the simulations in this subsection are always performed with a maximum of five iterations per time step.

Figure 3 shows the values of mean drag and lift coefficients varying with the dimensionless time step sizes ($\Delta tU/\Delta x$). Algorithms 1 and 2 are in tune with each other, as it can be seen from both $C_D^{{}^{\prime }}$ and $C_L^{{}^{\prime }}$ results. As for the kinematics, periodic vortex shedding characteristics are compared over Strouhal number "$Sr$” versus dimensionless time step size in Figure 4, where $Sr=fD/U$, with f being vortex shedding frequency in Hz. Similarly, Algorithms 1 and 2 yield coherent results in this aspect. The noticeable difference between the force coefficient values for $D/\Delta x=20$ and $D/\Delta x=30$ indicates that $D/\Delta x=20$ cannot provide a sufficient particle resolution. Yet, it can be inferred from Figure 5 that the mean velocity divergence, hence the continuity error, is comparable for these two particle resolutions. Figure 5 also indicates that increasing the time step sizes results in a better achievement of overall incompressibility in the domain. This is due to the fact that, since the time step size is a term in the penalty coefficient as $r_{AL}=C_{BV}{u}_{max}^2\Delta t$, the larger the time step size, the larger is the penalty coefficient, which affects the iterative guesses of the intermediate velocity field $u^{*}$ and the pressure field $p^{n+1}$ through Equations (8) and (10), respectively. Consequently, the penalty coefficient naturally influences the convergence characteristics of the iterative process. It should be stated that for simulations with five iterations per time step, the utilization of the time step beyond the maximum value that is represented in Figure 5 yields diverging simulations for both particle resolutions and both algorithms. Algorithms 1 and 2 show generally similar trend with increasing time step sizes. However, Algorithm 1 performs slightly better than Algorithm 2 in larger time step sizes.

It should also be noted that the simulations in our previous study [26] were performed with a lower dimensionless time step size as $\Delta tU/\Delta x=6.35\times {10}^{-4}$, which corresponds to the lower ends on the x axes of Figure 3, Figure 4 and Figure 5. Therefore, it is reasonable to conclude that the algorithm distinction does not noticeably affect the results of the previous study [26] due to the unoptimized time step sizes.

For the sake of completeness, in

Table 1. Comparison with two-dimensional (2D) numerical simulation results of literature for $Re=200$.

Source of the Result	${\mathit{C}}_{\mathit{D}}^{{}^{\prime }}$	${\mathit{C}}_{\mathit{L}}^{{}^{\prime }}$	$\mathit{S}\mathit{r}$
Algorithm 1	1.427	0.478	0.213
Algorithm 2	1.428	0.482	0.211
Zhang et al. [41]	1.423	0.529	0.203
Rajani et al. [42]	1.336	0.424	0.196
Marrone et al. [34]	1.380	0.680 *	0.200

* Mean lift amplitude.

, the results of two algorithms are compared with the 2D numerical findings of the literature in terms of mean drag and root mean square of lift coefficient as well as the Strouhal number for $Re=200$, where the simulations are performed with $D/\Delta x=30$ and $\Delta tU/\Delta x=5.08\times {10}^{-3}$.

At inner iterations, Algorithm 1 updates the particle positions by computing the derivative terms on previously predicted particle positions, whereas Algorithm 2 calculates the derivatives at initial particle positions with the updated field function values of $\mathbf{u}$ and p, as explained in Section 2.4. The effect of derivative locations on the accuracy is negligibly small because both schemes compute the new positions of particle for the inner iteration using Equation (12) based on ${\mathbf{r}}^n$.

4.2. Number of Iterations

In the previous subsection, it was concluded that the best performing dimensionless time step size is the largest applicable one for the simulations with five iterations per time step. Thus, the effect of the maximum iteration number on the results is investigated, utilizing this dimensionless time step value, namely $\Delta tU/\Delta x=5.08\times {10}^{-3}$ for the forthcoming simulations.

The same analysis pattern as the previous Section 4.2 is adopted in Figure 6, Figure 7 and Figure 8, investigating, respectively, force coefficients, Strouhal number, and the mean velocity divergence against the number of iterations per time step. Figure 6 and Figure 7 do not indicate any clear distinction between the two algorithms both in terms of force coefficients and Strouhal numbers with an increasing number of iterations. In Figure 8 a decreasing trend of mean continuity error is observed with an increasing number of iterations per time step for both resolutions.

In Figure 9, spatial variations of the velocity divergence and the magnitude of vorticity at similar lift extremities (around $tU/D\approx 88$ for each case) are given for $D/\Delta x=30$. Despite the almost identical results of vorticities for Algorithms 1 and 2, there are distinguishable differences for velocity divergences as the mean values that are depicted in Figure 8 suggest. Notwithstanding the fact that twenty iterations per time step yields the best results for both algorithms in general, five iterations per time step is an expedient choice, given the trade off between admissible accuracy and computational cost.

In Figure 10, pressure coefficients for Algorithms 1 and 2 with two, five, and 20 iterations are compared over the simulations with $D/\Delta x=30$. The results of $C_P^{{}^{\prime }}$ for the Algorithm 1 are unaffected by the number of iterations per time step, yielding overlapping mean pressure coefficients along the cylinder surface. Additionally, the differences between the results of Algorithms 1 and 2 are trivial for a kind of flow characteristic, which is fairly unsteady. It can be noticed that the mean drag and lift coefficient results of resolution $D/\Delta x=30$ in Figure 6 also support the findings of pressure coefficients, since they both demonstrate almost insensitive responses to the number of iterations per time step for Algorithm 1. The finite volume method results of Rajani et al. [42] and experimental results of Thom [43] are also included in this comparison to present supplementary information.

4.3. Reynolds Number

Following a thorough investigation of flow characteristics at $Re=200$, the two algorithms are compared at Reynolds numbers that range from 100 to 300. All of the simulations are performed utilizing $\Delta tU/\Delta x=5.08\times {10}^{-2}$ with a maximum of five iterations per time step, except for the $Re=300$ cases, which produce unstable simulations. Thus, for only at $Re=300$ cases, time step sizes are reduced to $\Delta tU/\Delta x=3.59\times {10}^{-2}$.

In Figure 11, the mean drag and lift coefficients against Reynolds numbers are presented. It can be observed that both $C_D^{{}^{\prime }}$ and $C_L^{{}^{\prime }}$ series are almost identical. Higher particle resolutions are required with increasing Reynolds numbers, which can be inferred through comparing the results of the current study with that of Zhang et al. based on the finite difference approach [41]. Moreover, the Strouhal number results of the two algorithms at different Reynolds numbers are compared in Figure 12. Correspondingly to the results of force coefficients (Figure 11), the Strouhal number results of Algorithms 1 and 2 match with each other. The mean velocity divergences of the flow fields display similar behaviour with the overall outlook of the results that were classified in the previous subsection. Figure 13 shows that Algorithm 1 produces slightly lower mean velocity divergence in the flow domain at all Reynolds numbers covered.

5. Conclusions

The aim of the current study is to investigate in detail the effect of time step size, number of inner iterations, and two different algorithms as to the calculation of spatial derivatives on the robustness, accuracy, and computational cost of the augmented Lagrangian SPH method. To be clear, Algorithm 1 utilizes the particle positions of the previous iteration to compute gradients and divergences of pressure and velocity fields, whereas Algorithm 2 employs the initial particle positions of the time step for these computations. As for the time step, a smaller numerical error is observed at larger time step sizes for both algorithms, which is associated with fact that the time step size is embedded in the penalty coefficient $r_{AL}$, which augments the prediction of the optimum pressure and velocity fields at each iteration. The largest applicable and best performing dimensionless time step size is found to be around $\Delta tU/\Delta x=5.08\times {10}^{-2}$ for $Re=200$. The further increase in the dimensionless time step size leads to instability and divergence in both of the algorithms, hence blowing up the simulation. The findings of the present study indicate that the dimensionless time step size that is utilized in our previous study based on the Algorithm 2 is in a safe range, such that it does not produce any difference from the Algorithm 1. It is observed that increasing the number of inner iterations improves the accuracy of simulations for both algorithms, as expected. Nevertheless, the accuracy that is gained by increasing the maximum number of iterations from five to twenty is not notably high to justify the additional computational cost.

Finally, the most efficient combination of the dimensionless time step size and the maximum number of iterations per time step, which are $\Delta tU/\Delta x=5.08\times {10}^{-2}$ and five iterations respectively, are tested at $Re=100$ to 300. Algorithms 1 and 2 yield overlapping results in terms of both force coefficients acting on the cylinder and periodic vortex shedding characteristics of the downstream. Overall, the results of the two algorithms do not meaningfully differ in terms of the dynamics and kinematics of the flow. However, Algorithm 1 can be deemed to be better, since it produces slightly smaller computational error at larger time step sizes.

Although in this study, the ALSPH method is extensively tested against time step size, the number of iterations, and locations of derivative operations for incompressible laminar flow, it is expedient to test it further for challenging free surface and turbulent flow problems, which will form the future direction of the current study.