Improvements to the convection scheme in the vortex particle method using a semi-Lagrangian method

Abstract

This study focuses on the convection scheme in the vortex particle method. The usual two-step convection procedure for vortex elements, which comprises the convection of vortex elements on each grid point and redistribution of vorticity into each grid point, is extremely difficult to parallelize because of unavoidable memory-write conflictions in the redistribution step, thereby limiting the computational speed. We apply a semi-Lagrangian method called the cubic semi-Lagrangian method, as the convection scheme in the vortex particle method to improve space–time accuracy and accelerate computations. Two test problems, namely the Taylor–Green vortex and double shear flow, are simulated using the semi-Lagrangian vortex particle method proposed in this study to confirm its suitability. The results show that the accuracy of capturing standing vortices and the total amount of kinetic energy conservation are significantly improved. The performance measurements show that the average execution time for the convection of vortex elements is reduced to one-half, thus verifying the semi-Lagrangian method as a suitable convection scheme for the vortex particle method.

Introduction

The vortex methods^1,2 are Lagrangian numerical method used to solve unsteady flow problems. These methods discretize the vorticity fields into vortex elements and trace the convection of each vortex element to simulate the time evolution of the flows. There are several ways to compute the convection of vortex elements. The Lagrangian vortex method traces the convection of vortex elements according to the Biot–Savart law. A method, known as the vortex in cell (VIC) method, solves the Poisson equation for streamfunction or vector potential on a spatial grid to determine the convection of vortex elements. When vortex elements are redistributed onto each grid point at certain intervals, the VIC method is referred to as the vortex particle method (VPM).

The VPM has some characteristics that are superior to those of the other methods, such as a higher numerical stability and higher affinity for the Eulerian method. When using VPM to simulate flows, Eulerian methods, including the finite difference method (FDM) and the finite element method, can be used to compute non-convection procedures such as viscous diffusion and vortex stretching. Several studies have examined the suitability of the VPM for basic flows such as homogeneous isotropic turbulence, ³ time evolution of three-dimensional (3D) Taylor-Green vortices, ⁴ and flow around an obstacle. ⁵ The VPM has been improved by introducing multi-hierarchy ⁶ and by constructing the redistribution functions to achieve total variation diminishing.^7,8 Areas of application for the VPM have been extended to direct numerical simulation (DNS) of a turbulent channel flow, ⁹ large-scale parallel computation, ¹⁰ multiphase flows, ¹¹ and aeroacoustics. ¹² In recent years, the VPM is successfully used in the field of computer graphics (CG). ¹³

However, the VPM exhibits some problematic characteristics. Numerical oscillations appear near steep vorticity gradients and spread throughout the whole flow field. Parallelization of the redistribution in VPM using shared memory parallel programing features, e.g. OpenMP, on shared memory systems is also quite difficult. These characteristics limit the computational speed of VPM.

The convection of vortex elements in VPM is constructed from the following two steps.

The convection of each vortex element on each grid point

The redistribution of vorticity of vortex elements onto grid points

To prevent numerical oscillations, linear multistep integration methods like the modified Euler method is imposed ¹⁴ as step 1. For step 2, higher order redistribution functions ⁴ and redistribution functions achieving the total variation diminishing^7,8 have been constructed. Rossinelli and Koumoutsakos ¹⁵ implemented the VPM on Graphics Processing Units (GPUs) for fast incompressible flow simulations. However, it is impossible to use their implementation for programs running on Central Processing Units (CPUs) because their implementation has to use a dedicated hardware called “graphics pipeline” on GPUs. Sbalzarini et al. ¹⁶ presented a library, named parallel particle–mesh (PPM), to provide a set of parallelized functions for particle methods including the VPM. PPM is applied to large-scale parallel computations of compressible vortex rings, ¹⁰ achieving the high parallel efficiency. On the other hand, since only few redistribution functions are provided by PPM, the above-mentioned improved redistribution functions^4,7,8 are not available. This means that the parallelization of the VPM and prevention of numerical oscillation cannot be simultaneously achieved when the PPM is utilized.

This study aims to solve the above-mentioned defects of VPM by applying the semi-Lagrangian method as the convection scheme in VPM. Semi-Lagrangian methods are numerical methods used to solve linear convection equations. The semi-Lagrangian method has advantages, like the VPM, of both the Eulerian and the Lagrangian methods. The semi-Lagrangian method is a Lagrangian method in terms of tracing particles that have physical quantities. However, it is also an Eulerian method in terms of handling physical quantities on each grid point, since the semi-Lagrangian method focuses only the particles reaching at grid points. The computational procedure of the semi-Lagrangian method consists of following three steps:

Find upstream points believed to be particle positions for each particle on the grid points.

Estimate the physical quantities of each particle using nearby grid points.

Move each particle to each grid point.

In step 2, an interpolation polynomial is constructed. This polynomial determines the temporal and spatial accuracy of the semi-Lagrangian method. It is known that the interpolation polynomial achieves both a higher order of accuracy and numerical oscillation suppression. In addition to those advantages, the method can be easily parallelized on shared and distributed memory systems.

In the field of CG, semi-Lagrangian methods are widely used for simulations of fluid flows.^17–19 These methods solve Navier–Stokes’ equation. Some simulation methods using “vortex particles”^20,21 or simulation methods called “vortex particle method”^22–24 are used for real-time simulations of smoke, water, and explosions. In fact, these methods are closer to the VIC method than the VPM, according to the classification mentioned above. Since these methods need to simultaneously solve the vorticity and Navier–Stokes’ equations, these methods are not classified into simulation methods based on the vorticity equation.

Only a few studies applying semi-Lagrangian methods to the vorticity equation are found in the field of meteorology. Sawyer ²⁵ applied a first-order semi-Lagrangian method to the two-dimensional vorticity equation and demonstrated that the semi-Lagrangian method is suitable for numerical predictions with long time steps. This study showed qualitative comparisons of isobaric lines obtained from FDM and the proposed method. It also listed correlations between observed and predicted atmospheric pressure distributions. Satoh ²⁶ developed a higher-order semi-Lagrangian method combined with the constrained interpolation profile (CIP) method ²⁷ to trace vortices in the atmosphere. The method simultaneously solves the vorticity and Navier–Stokes’ equations as similar to the methods used in the field of CG.^22–24 Therefore, there are no comprehensive researches on the properties including the accuracy, conservativity, and execution speed of the semi-Lagrangian methods for solving the vorticity equation.

This paper proposes a semi-Lagrangian Vortex Particle method (SLVPM), applying the semi-Lagrangian method as the convection scheme in VPM. SLVPM is applied to simulate two-dimensional test problems to investigate the differences of convection schemes in flow patterns, kinetic energy conservation properties, and execution times.

Numerical methods

Governing equations

The mass and momentum conservation equations for two-dimensional, incompressible, and inviscid flows are expressed as

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$

(1)

$$\{\begin{array}{ll}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}\hfill & \hfill \\ \frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial y}\hfill & \hfill \end{array}$$

(2)

where $(x,y)$ and t represent space and time, u and v are velocity components in x- and y-directions, p is pressure, ρ is density. Although viscous diffusion is not considered in evaluating the differences between the convection schemes, it does not limit the application of SLVPM to inviscid flows only.

Taking the curl of equation (2) and substituting equation (1), the vorticity equation is derived as

$$\frac{\partial \omega }{\partial t}+u\frac{\partial \omega }{\partial x}+v\frac{\partial \omega }{\partial y}=0$$

(3)

where ω is the vorticity around the axis perpendicular to the x–y plane and defined as follows

$$\omega =\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$$

(4)

The streamfunction ψ satisfying the following equations is introduced in order to calculate the velocity from the vorticity

$$\{\begin{array}{ll}u=\frac{\partial \psi }{\partial y}\hfill & \hfill \\ v=-\frac{\partial \psi }{\partial x}\hfill & \hfill \end{array}$$

(5)

Substituting equation (5) into equation (4), the Poisson equation for the streamfunction ψ can be derived. The result is as follows

$$\frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial y^2}=-\omega$$

(6)

Vortex particle method

In VPM, the non-convection procedure is computed by Eulerian methods. The finite difference method (FDM) is widely used as the Eulerian method. The regular grid, which defines all of the physical variables on the same grid points, is usually used to divide the computational domain into the streamfunction-vorticity method and the VIC method. Uchiyama and coworkers^9,14 have successfully clarified that dividing the computational domain by a staggered grid, which defines the physical variables at different points, can ensure the consistency among the discretized equations as well as prevent the numerical oscillation of the flow field. In this study, the staggered grid is used for spatial discretization of the computational domain. As shown in Figure 1, the vorticity and the streamfunction are defined at each grid point. The velocity components are defined at the middle of the grid side.

Figure 1.

Two-dimensional staggered grid system.

The convection of the vortex elements are estimated by Lagrangian computation of the following equation

$$\{\begin{array}{ll}\frac{d{x}_v}{dt}=u(x_v,y_v)\hfill & \hfill \\ \frac{d{y}_v}{dt}=v(x_v,y_v)\hfill & \hfill \end{array}$$

(7)

where (x_v, y_v) is the position of each vortex element. The convection velocity is obtained by interpolating the velocity defined at the middle of the grids. After determining the convection of the vortex elements by equation (7), the vorticity of the vortex elements is redistributed onto the grid points by

$$\omega ={\displaystyle \sum _{\text{g}=1}^{{\text{N}}_{\text{g}}}{\displaystyle \sum _{\text{e}=1}^{{\text{N}}_{\text{v}}}{\omega }^{\text{e}}}}M\left(\frac{x^{\text{g}}-x_v^{\text{e}}}{\Delta x}\right)M\left(\frac{y^{\text{g}}-y_v^{\text{e}}}{\Delta y}\right)$$

(8)

where ${\text{N}}_{\text{v}},{\text{N}}_{\text{g}}$ are the number of vortex elements and grid points, respectively. $(x_v^{\text{e}},y_v^{\text{e}})$ is the position vector for vortex element $\text{e}$ , $(x^{\text{g}},y^{\text{g}})$ is the position for the grid point $\text{g}$ , and $\Delta x$ and $\Delta y$ are the grid width of the x– and y– directions, respectively. For the redistribution function M, the following B-spline function ³ is applied

$$M(\sigma )=\{\begin{array}{ll}1-\frac{5}{2}{\sigma }^2+\frac{3}{2}|\sigma {|}^3\hfill & \left|\sigma \right|<1\hfill \\ \frac{1}{2}{\left(2-\left|\sigma \right|\right)}^2\left(1-\left|\sigma \right|\right)\hfill & 1\le \left|\sigma \right|\le 2\hfill \\ 0\hfill & \left|\sigma \right|>2\hfill \end{array}$$

(9)

When the flow at t is known, the time evolution is simulated by the following procedure

Calculate the convection of the vortex elements from equation (7).

Redistribute ${\omega }^{\text{e}}$ onto grid points from equation (8).

Calculate ψ by solving equation (6).

Calculate velocity from equation (5).

According to equations (8) and (9), the vorticity ${\omega }^{\text{e}}$ of a vortex element is redistributed onto 16 grid points around the position of the vortex element. When redistributing the vorticity ${\omega }^{\text{e}}$ onto $(x^{\text{g}},y^{\text{g}})$ , the memory-write conflicts known as data races occur because it is necessary to accumulate the redistributed vorticities of all vortex elements near $(x^{\text{g}},y^{\text{g}})$ . This “scatter”-type algorithm results in difficulty for parallelization on shared memory systems. To avoid this difficulty, this study applies the semi-Lagrangian method to VPM.

Semi-Lagrangian method

Semi-Lagrangian methods are numerical methods used to compute linear convection equations. Semi-Lagrangian methods use characteristic of linear convection equations for computation. The linear convection equation with constant convection velocity c is expressed as

$$\frac{\partial f}{\partial t}+c\frac{\partial f}{\partial x}=0$$

(10)

where f is a physical quantity. The exact solution of equation (10) is $f=f(x-ct,t)$ . This solution shows that the quantity f convects at a distance ct with velocity c and maintains its spatial distribution. Semi-Lagrangian methods postulate that the distribution of f at $t+\Delta t$ is the same as the distribution of f at t with convection at a distance ct; $f(x_{\text{i}},t+\Delta t)=f(x_{\text{i}}-ct,t)$ , as shown in Figure 2, where $x_{\text{i}}$ is a position of the grid point i. Since the upstream point of the grid point i, namely $x_{\text{i}}-ct$ , does not fall on a grid point in almost all cases, the value of f at $x_{\text{i}}-ct$ is estimated by nearby values from the interpolation polynomial F. This study adopts the third-order interpolation polynomial with upstream weighting. Assuming that $c\ge 0$ and $c\Delta t<\Delta x$ , the interpolation polynomial is defined by equation (11)

$$F_i^n=c_3{(x-x_{\text{i}})}^3+c_2{(x-x_{\text{i}})}^2+c_1(x-x_{\text{i}})+f_i$$

(11)

where superscript n is the time step. The coefficients $c_1,c_2,c_3$ are obtained under the conditions $F_i^n(x_{\text{i}-1})=f_{i+1}^n,F_i^n(x_{\text{i}})=F_i^n$ , and $F_i^n(x_{\text{i}+1})=f_{i+1}^n$

$$\begin{array}{l}{c}_3=\frac{f_{i+1}^n-3F_i^n+3f_{i+1}^n-f_{i-2}^n}{6\Delta x^3}\hfill \\ c_2=\frac{f_{i+1}^n-2F_i^n+f_{i+1}^n}{2\Delta x^2}\hfill \\ c_1=\frac{2f_{i+1}^n+3F_i^n-6f_{i+1}^n+f_{i-2}^n}{6\Delta x}\hfill \end{array}$$

(12)

Figure 2.

A schematic of the semi-Lagrangian method.

A semi-Lagrangian method which estimates the upstream f from equations (11) and (12) is known as the cubic semi-Lagrangian method. Rewriting equation (12) to clearly show the similarity with equations (8) and (9) has been shown as follows

$$f^{n+1}(x_{\text{i}}-ct)={\displaystyle \sum _{\text{m}=\text{i}-2}^{\text{i}+1}{f}_m^n}W\left(-\frac{c\Delta t}{\Delta x}\right)$$

(13)

$$f_{\text{i}}^{n+1}=f^{n+1}(x_{\text{i}})=f^{n+1}(x_{\text{i}}-ct)$$

(14)

$$W(\sigma )=\{\begin{array}{ll}-{\sigma }^3+\sigma \hfill & \text{m}=\text{i}-2\hfill \\ 3{\sigma }^3+3{\sigma }^2-6\sigma \hfill & \text{m}=\text{i}-1\hfill \\ -3{\sigma }^3-6{\sigma }^2+3\sigma +1\hfill & \text{m}=\text{i}\hfill \\ {\sigma }^3+3{\sigma }^2+2\sigma \hfill & \text{m}=\text{i}+1\hfill \\ 0\hfill & \text{elsewhere}\hfill \end{array}$$

(15)

When the convection velocity is negative, the interpolation polynomial is constructed from points $\text{i}-1,\text{i},\text{i}+1,\text{i}+2$ . In a two-dimensional simulation, the value of f at upstream of grid point $(\text{i},\text{j})$ is estimated from 16 points around $(\text{i},\text{j})$ and updated by

$$\begin{array}{c}{f}^{n+1}(x_{\text{i}}-ct,y_{\text{j}}-ct)\\ ={\displaystyle \sum _{\text{n}=\text{j}-2}^{\text{j}+1}{\displaystyle \sum _{\text{m}=\text{i}-2}^{\text{i}+1}{f}_{m,n}^n}}W\left(-\frac{c_x\Delta t}{\Delta x}\right)W\left(-\frac{c_y\Delta t}{\Delta y}\right)\end{array}$$

(16)

$$f_{\text{i},\text{j}}^{n+1}=f^{n+1}(x_{\text{i}},y_{\text{j}})=f^{n+1}(x_{\text{i}}-c_{x}t,y_{\text{j}}-c_{y}t)$$

(17)

where c_x and c_y are the velocity components of the x and y directions, respectively. It is noted that equation (16) is formulated for $c_x>0,c_y>0$ .

In the cubic semi-Lagrangian method, the convection of a particle at the grid point $(\text{i},\text{j})$ can be parallelized easily and simply on the shared memory systems because the value of f at the 16 points around the particle is referred to only. This “gather”-type algorithm avoids the memory-write conflicts. In addition, numerical oscillations can be suppressed through cubic interpolation. Although small overshoots and undershoots appear near steep vorticity gradients, the numerical oscillations do not spread throughout the whole flow field.

Numerical tests

To evaluate the suitability of the semi-Lagrangian vortex particle method for incompressible and inviscid flows, two test problems, Taylor-Green vortex and double shear flow, ²⁸ are simulated.

Taylor–Green vortex

The Taylor–Green vortex is one of the exact solutions of the Euler and the Navier–Stokes equation. In the Taylor–Green vortex, vortices with different rotating direction are arranged periodically in the x and y directions. It is assumed that the flow field has a periodicity of length $2\pi$ in both directions, and the streamfunction ψ and the vorticity ω are expressed as

$$\psi =\sin x\sin y$$

(18)

$$\omega =2\sin x\sin y$$

(19)

In the inviscid flow, the vortices maintain the above initial conditions. The computational domain with side lengths $L_x\times L_y=[0,2\pi ]\times [0,2\pi ]$ is divided by staggered grids into ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=32\times 32,64\times 64,128\times 128$ , and 256 × 256. Time evolutions of the Taylor–Green vortex are simulated in a duration from t = 0 to 10 s. The Courant numbers based on the maximum velocity $U_{\max }$ and grid width $\Delta x$ are set to 0.01 and 0.1, to consider the effect of the time interval $\Delta t$ . The periodic boundary condition is imposed on all boundaries. The Poisson equation for the streamfunction is solved using a Fast Fourier Transformation (FFT). A FFT subroutine named DFTI included in the Intel Math Kernel Library is used. The simulations applying FDM and VPM as the convection schemes are also performed in order to compare the results. For the FDM, the convection term is discretized so as to conserve the total amount of kinetic energy and the first-order Euler method is used for time integrations. In this case, the total amount of kinetic energy is not affected by the convection term, but kinetic energy varies throughout time integrations due to the truncation errors resulting from the Euler method. For VPM, the first-order Euler method is also applied to equation (7).

All simulations have been performed on a personal computer (Windows 7, CPU Intel Xeon E5–1620 3.5GHz, 32GB memory) and the simulation programs are compiled by Intel Parallel Studio XE2017 Composer for Fortran with options/O3/QaxCORE-AVX512/QxHost/Qmkl:parallel. The programs run on all CPU cores (four physical and four logical cores) with eight threads to achieve the maximum performance.

Figure 3 shows vorticity distributions at t = 10. In Figure 3, the vorticity distribution at t = 0 is depicted for comparison. Comparing Figure 3(a) and (b) shows that FDM completely conserves the initial vorticity distribution. In the result from VPM, vortices are diffused and the decrement of those peak values is also observed in Figure 3(d). The SLVPM, shown in Figure 3(c), captures the peak value of vortices, although the vorticity gradient becomes smaller.

Figure 3.

Vorticity distributions of Taylor–Green vortices at t = 10 s ( ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=32\times 32$ ).

Figure 4 shows the time variation of the total amount of kinetic energy $K_{\text{total}}$ . $K_{\text{total}}$ is calculated and normalized as

$$K_{total}^n=\frac{{\displaystyle {\sum }_{\text{j}=1}^{{\text{N}}_{\text{y}}}{\displaystyle {\sum }_{\text{i}=1}^{{\text{N}}_{\text{x}}}{\omega }_{i,j}^n}}{\omega }_{i,j}^n\Delta x\Delta y}{{\displaystyle {\sum }_{\text{j}=1}^{{\text{N}}_{\text{y}}}{\displaystyle {\sum }_{\text{i}=1}^{{\text{N}}_{\text{x}}}{\psi }_{\text{i},\text{j}}^0}}{\omega }_{\text{i},\text{j}}^0\Delta x\Delta y}$$

(20)

where ${\psi }^0$ and ${\omega }^0$ are the initial values of ψ and ω, respectively. FDM conserves $K_{\text{total}}$ within the range of round-off error of the floating-point number. In contrast to this, VPM and SLVPM do not conserve and spurious dissipation occurs. The dissipation rate of $K_{\text{total}}$ for both VPM and SLVPM is almost the same.

Figure 4.

Time variation of the total value of kinetic energy, $K_{\text{total}}$ ( ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=32\times 32$ ).

Figure 5 shows the energy conservation errors versus the grid width $\Delta x$ . Hereafter, a positive value of the error is called “production” and a negative one is called “dissipation.” It should be noted that the results of FDM do not appear in Figure 5 because the FDM conserves $K_{\text{total}}$ as mentioned above. The decrement of the dissipation of $K_{\text{total}}$ with decreasing $\Delta x$ shows VPM and SLVPM both have the property of grid convergence. The gradient of the dissipation in VPM is proportional to $\Delta x$ at both Courant numbers. In the result from SLVPM, the gradient is proportional to $\Delta x$ at the Courant number 0.1. However, the gradient is proportional to $\Delta x$ to the power of 1.5 at the Courant number 0.01. This means that the substantial accuracy in space of SLVPM for the conservation of the total amount of kinetic energy is higher than that of VPM.

Figure 5.

Energy conservation errors in the Taylor–Green vortex at t = 10 s.

Double shear flow

In this section, the double shear flow is simulated to evaluate the suitability of SLVPM for flows with strong nonlinearity. In the two-dimensional domain with side lengths $L_x\times L_y=[0,2\pi ]\times [0,2\pi ]$ , the velocity distribution is given by the following equation

$$\begin{array}{l}u=\{\begin{array}{ll}\tanh \left[\left(y-\frac{L_y}{4}\right)/\epsilon \right]& y\le \frac{L_y}{2}\\ \tanh \left[\left(\frac{3L_y}{4}-y\right)/\epsilon \right]& y>\frac{L_y}{2}\end{array}\\ v=\delta \sin \left(\frac{2\pi x}{L_x}\right)\end{array}$$

(21)

where ε is the thickness of the shear layers and δ is the amplitude of perturbation in the shear layers. Taking the curl of equation (21), the initial vorticity is obtained below

$$\omega =\{\begin{array}{ll}\frac{2\pi }{L_x}\delta \cos \left(2\pi x/L_x\right)\hfill & \hfill \\ -\frac{1}{\epsilon {\cosh }^2\left[\left(y-\frac{L_y}{4}\right)/\epsilon \right]}\hfill & y\le \frac{L_y}{2}\hfill \\ \frac{2\pi }{L_x}\delta \cos \left(2\pi x/L_x\right)\hfill & \hfill \\ +\frac{1}{\epsilon {\cosh }^2\left[\left(\frac{3L_y}{4}-y\right)/\epsilon \right]}\hfill & y>\frac{L_y}{2}\hfill \end{array}$$

(22)

where $\epsilon =L_y/30$ and $\delta =0.05$ are used according to the existing study. ²⁸ The domain is divided by the staggered grid into ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=128\times 128,256\times 256$ , and 512 × 512. The simulations are performed in a duration from the non-dimensional time, defined as $U_{\max }t/L_x,t^{*}=0$ to 1.8 ²⁸ with the Courant numbers based on the maximum velocity in the x direction $U_{\max }$ and grid width $\Delta x$ defined as 0.01, 0.1. The periodic boundary condition is imposed on all boundaries. FFT is used to solve the Poisson equation of the streamfunction.

Figures 6 and 7 are the vorticity distributions at $t^{*}=0.6$ and 1.8 with ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=256\times 256$ and Courant number of 0.01. At $t^{*}=0.6$ , the distributions simulated by FDM and SLVPM are almost parallel, as shown in Figure 6(a) and (b). The result of VPM shown in Figure 6(c) is slightly different from both FDM and SLVPM. The shear layers roll up and strong thin shear layers appear at $t^{*}=1.8$ . Once shear layers start to roll up, numerical oscillations appear and spread throughout the whole domain by FDM, causing the simulation to eventually blow up. On the vorticity distribution in VPM, numerical oscillations appear near the stretched-thin shear layers. As with the result of the Taylor–Green vortex, the vorticity distribution is diffused. SLVPM favorably captures the shear layers and successfully suppresses numerical oscillations. This vorticity distribution is comparable to those of previous existing results. ²⁸

Figure 6.

Vorticity distributions of double-shear flow at $t^{*}=0.6$ ( ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=256\times 256$ ). (a) Finite Difference, (b) Present and (c) Vortex Particl.

Figure 7.

Vorticity distributions of double-shear flow at $t^{*}=1.8$ ( ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=256\times 256$ ). (a) FDM, (b) Present and (c) Vortex Particle.

The time variations of the total amount of kinetic energy $K_{\text{total}}$ are shown in Figure 8. The results for ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=256\times 256$ and Courant number of 0.1 are shown. No schemes conserve $K_{\text{total}}$ . $K_{\text{total}}$ increases with time marching when FDM and VPM are used. On the other hand, $K_{\text{total}}$ decreases in SLVPM. Since the production of kinetic energy leads to the divergence of simulations, this characteristic of SLVPM is considered to be an advantage in terms of stabilization of simulations. FDM initially conserves total kinetic energy; however, $K_{\text{total}}$ increases suddenly just after the shear layers roll up. For VPM, the positive gradient of $K_{\text{total}}$ becomes larger after the shear layers roll up. Although the negative gradient of $K_{\text{total}}$ for SLVPM also changes after the shear layers roll up, the magnitude of the gradient is smaller than those resulting from VPM.

Figure 8.

Time variations of the total value of kinetic energy, $K_{\text{total}}$ ( ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=256\times 256$ ).

Figure 9 shows the energy conservation errors at $t^{*}=1.8$ versus the grid width $\Delta x$ . When the grid is coarse, the error from VPM is smaller than that from SLVPM. The error increases with decreasing width $\Delta x$ . This phenomenon means that VPM does not have grid convergence for kinetic energy conservation. Therefore, VPM’s suitability for flows with strong nonlinearity becomes questionable. In contrast to this, the error of SLVPM decreases with decreasing width $\Delta x$ , demonstrating that the grid convergence for kinetic energy conservation is sustained.

Figure 9.

Energy conservation errors of double-shear flow at $t^{*}=1.8$ .

Comparison of execution time

The execution times for the simulations applying the three different convection schemes are measured using Windows system time, and the results are listed in Table 1. The measurements are performed for double-shear flow with ${\text{N}}_{\text{x}}\times {\text{N}}_{\text{y}}=512\times 512$ and a Courant number of 0.01. The execution time during $t^{*}\le 1.8$ was measured and time per computational step is obtained. The execution times for non-optimized program are also listed in Table 1. The execution times are measured when the programs are executed on all CPU cores (four physical and four logical cores) with eight threads. The execution times using SLVPM are shorter than those using VPM, regardless of program optimization. The differences between execution times for VPM and SLVPM can be attributed to the convection algorithm. In VPM, after determining the convection of a vortex element, the vorticity of the vortex element is redistributed onto 16 grid points around the element. When multiple vortex elements are redistributed with their vorticity onto the same grid point, memory-write conflictions occur. Therefore, the parallelization of this “scatter”-type algorithm on shared memory system is difficult. SLVPM avoids such conflictions and can easily be parallelized according to its “gather”-type algorithm. It should be noted that the difference between execution times per step is shorter than 0.1 ms when assignment operations to the array handling discretized vorticity in the subroutine of VPM and SLVPM are commented out. Consequently, the difference in execution times is caused by algorithm differences rather than by optimization level of each subroutine.

Table 1.

Execution times of numerical convection schemes.

	Exection time (ms)		Speedup ratio (−)
	Non-optimized	optimized
Finite difference	5.12	757 $\times {10}^{-3}$	6.8
Vortex particle	41.1	10.3	4.0
Present	17.5	4.65	3.8

The details of the execution time are shown in Figure 10. In the legends in Figure 10, “Convection”, “Poisson eq.”, and “Velocity” represent the convection of vortex elements, solving the Poisson equation of streamfunction by FFT, and the computation of convection velocities, respectively. The execution times of the convection of vortex elements in SLVPM is less than half that of those for VPM, and the total execution time is reduced to 60%, demonstrating the effectiveness of SLVPM.

Figure 10.

Computational time per step and details.

Conclusion

In this study, a semi-Lagrangian vortex particle method is proposed to solve the defects in the vortex particle method, such as numerical oscillations and difficulty in parallelization on shared memory systems. The method uses the cubic semi-Lagrangian method for determining the convection of vortex elements. The cubic semi-Lagrangian method has advantages of both Eulerian and Lagrangian methods, like the vortex particle method. The proposed method is applied to simulate basic test problems to comprehensively investigate the properties of the proposed method. The results are as follows:

In flow field with standing vortices like the Taylor–Green vortex, SLVPM favorably captures the peak of vortices. The method does not conserve the total amount of kinetic energy, but the errors of kinetic energy conservation converge with decreasing grid width. The magnitude of the error is comparable with that from VPM.

In a flow with strong nonlinearity, the kinetic energy conservation errors of the SLVPM converge with decreasing grid width. By contrast, the error from VPM does not converge. This characteristic of SLVPM ensures that the error of the kinetic energy conservation is improved as the grid becomes finer, demonstrating the suitability of SLVPM for flows with a high Reynolds number.

SLVPM is two times faster than VPM. The advanced speed is due to the memory-write conflicts, which appear in the redistribution step in VPM but do not occur in SLVPM.