Parallel adaptive weakly-compressible SPH for complex moving geometries

We use the dummy particle approach of Adami et al. [61] to enforce the boundary conditions. We enforce no-penetration by using the wall velocity for the dummy particles in the EDAC equation. We enforce the no-slip and free-slip boundary condition by, first, extrapolating the fluid velocity onto the dummy particles using Shepard interpolation:

where $j$ is the index of the fluid particles in the neighborhood of the $i^{\text{th}}$ dummy particle; then, the velocity of the dummy particles is computed from,

where the subscript $w$ denotes the dummy wall particles, $\mathbf{u}_{i}$ is the prescribed wall velocity. To impose the pressure gradient accurately near the wall, the pressure on the dummy wall particles is computed by,

where the subscript $f$ denotes the fluid particles, $\mathbf{a}_{w}$ is the acceleration of the wall, $g$ is the acceleration due to gravity (this is zero in the present work), $r_{wf}=|\mathbf{r}_{w}-\mathbf{r}_{f}|$, and $h_{wf}=(h_{w}+h_{f})/2$.

We use the no-reflection boundary conditions on the inviscid wall of all simulations. For this we use the approach of Lastiwka et al. [62], where we compute the characteristic properties, $J_{1},J_{2}$, and $J_{3}$, op. cit., from the fluid and extrapolate these properties onto the dummy particles describing the boundary using Shepard interpolation, finally, the fluid properties on the dummy particles are set from these characteristic properties.

The forces due to the pressure and viscosity of the fluid on the solid are computed by evaluating,

which in the variable-$h$ SPH discretization is written as,

where the summation index $j$ is over all the fluid particles in the neighborhood of a $i^{\text{th}}$ dummy particle. The drag $C_{d}$, lift $C_{l}$, and pressure $C_{p}$ coefficients used for comparison of our results with reference data are expressed by:

where $U_{\infty}$ is free-stream velocity of the fluid, $\mathbf{e}_{x}$ and $\mathbf{e}_{y}$ are the unit vectors in the $x$ and $y$ directions, respectively, $\mathbf{f}_{\text{p}}+\mathbf{f}_{\text{visc}}=\sum_{j}\mathbf{f}_{j,\text{p}}+\mathbf{f}_{j,\text{visc}}$ is the sum of all forces of the dummy particles, and $D$ is the characteristic length.

We use the Predict-Evaluate-Correct (PEC) integrator to integrate the position, velocity, and pressure. First, we predict the properties at an intermediate time $n+\frac{1}{2}$ by,

next, we evaluate the rate of change of properties at $n+\frac{1}{2}$, then, we correct the properties to get the corresponding values at the new time $n+1$ by,

The time-step is determined by the highest resolution used in the domain, and the minimum of the Courant–Friedrichs–Lewy (CFL) criterion and the viscous condition is taken:

PySPH is an open-source, Python-based, multi-CPU, and GPU framework for SPH simulations. We have used the Python based interfaces PySPH provides to implement most of our adaptive algorithm. PySPH exposes Python-friendly interface to the user and using internal mechanism generates a high performance serial or parallel code based on the choice of the user. For an overview of the PySPH design, see Ramachandran et al. [49] and https://pysph.readthedocs.io for a detailed reference.

We use the Tool interface to split and merge the particles. The Tool interface provides a way for a user defined callback to run before or after the integration stages by overloading the pre_stage or post_stage methods, respectively. It also provides a way to run a user callback before the start of the integration step by overloading the pre_step method. Listing LABEL:lst:tool shows a sample Python code used in our adaptive particle refinement algorithm. The class AdaptiveResolution overloads the pre_step method of Tool to run a user callback self.run which in turn calls other specific functions related to splitting, merging, and shifting.

These functionalities enable one to write the code without having to pay much attention to the low-level high-performance complexities of implementing on different platforms like multi-core CPUs, or GPUs. The respective splitting and merging algorithms are all implemented in pure Python. For complete details of our implementation one can see the source code provided in https://gitlab.com/pypr/asph_motion.

In section 5.4 we solve a case with a stationary cylinder and another with a moving square to demonstrate that we do not require the background particles. In section 5.6 we demonstrate the parallel performance of our algorithms by solving a problem and studying the speedup obtained as we increase the number of threads.

We study the flow around two stationary cylinders in a side-by-side configuration at $Re=100$. The flow produces distinct wake patterns at different gap spacing between the cylinders. This problem is studied first by Kang [64] and further used as a validation test case by Bao et al. [65] but do not have SPH related studies. The diameter $D$ of each of the cylinders is $1m$. The dimensions of the computational domain are based on $D$, and the schematic diagram of the computational domain is shown in fig. 2. For each inlet, outlet, and wall, we use $10$ layers of particles to get full kernel support. Three different cases of center-to-center cylinder spacing $S=1.2D,1.5D$ and $4D$ are considered for simulations. A constant velocity is prescribed at the inlet, $U_{\infty}=1m/s$, and the kinematic viscosity $\nu$ is given by $U_{\infty}D/Re$. Other geometric and numerical parameters are provided in Table 1.

For all the cases, we deploy particles in the range of $75k$ to $88k$ on a $75$ m $\times\ 60$ m domain. The maximum spacing of fluid particles is $D/\Delta x_{\max}=2$ and the minimum resolution of the fluid particles, which also matches the resolution of the boundary, is $D/\Delta x_{\min}=100$; the ratio of which is $50$ times smaller. Without an adaptive resolution, the problem would need $45$ million particles to simulate at such accuracy using the minimum resolution.

We are able to observe the three kinds of vortex shedding wake patterns that occur for each configuration. The wake patterns are (i) single bluff body at $S=1.2D$, (ii) a random flip-flopping gap at $S=1.5D$, and (iii) anti-phase synchronized wake at $S=4.0D$. Time evolution of the simulated results of total drag and lift force coefficients are presented in figs. 3, 4 and 5. The simulations took a fairly long wall-clock time for $250$ seconds of simulation time, each case requiring an average wall-clock time of $40$ hours on a server with $16$ cores and $32$ threads for the adopted resolution.

Figure 3 shows the comparison of time history of the coefficients of total drag and lift for spacing $S=1.2D$ up to $T=U_{\infty}t/D=250$. The figure illustrates the force characteristics of a single bluff body type with entirely suppressed gap flow observed at $(S=1.2D)$.

Our results show a phase difference which may be caused by the differences in the nature of the solvers, time-stepping used, and boundary conditions. The results are in good agreement with the reference data of [65]. The trend, as well as the magnitude of the curves, is matching well.

Figure 4 shows the force characteristics for the spacing $S=1.5D$. The flow is complicated and belong to a randomly flip-flopping flow for which the wake pattern looks irregular [65]. Our results compare very well at the beginning up to around $75$ seconds of the flow. The results show some discrepancy in the $C_{d}$ afterward but stay in the trend, matching the amplitude and frequency. The $C_{l}$ results match better than the $C_{d}$ in all the cases for both the top as well as bottom cylinders.

Figure 5 displays the comparison of results for the coefficients drag and lift force at spacing $S=4.0D$. The results are matching well with the $C_{d}$ of the reference data in fig. 5 (a) and fig. 5 (b) of the upper and lower cylinders with differences less than $5$ percent. The results of $C_{l}$ are in a very good agreement as shown in fig. 5 (c) and fig. 5 (d) for both the cylinders. The force characteristics of the flow at this gap between the cylinders is anti-phase synchronized where the vortex is shed behind each cylinder independently and symmetrical to the centerline.

Figure 6 exhibits the vorticity distribution corresponding to the three spacing ratios. The results indicate that the flow characteristics and wake structure explained in Figs. 3 to 5 are justified. The flow structure strongly depends on the spacing gap for the given Reynolds number. A suppressed flow, seemingly a single bluff-body wake pattern is shown by the smallest spacing of $S=1.2D$ in fig. 6(a). In fig. 6(b) the gap flow flip-flops for $S=1.5D$ and seems irregular. The wake pattern changes to a synchronized flow as the spacing increases. As shown in fig. 6(c), the instantaneous wake is anti-phase synchronized; the vortex shedding occurs completely behind each cylinder which is symmetric with the wake centerline.

The behaviour of the flows around two stationary cylinders and a single oscillating cylinder, or both is helpful to understand the mechanism of flow around two or more oscillating cylinders. The present SPH model is able to simulate and accurately identify the wake characteristics of the flow. The force quantities are in general in good agreement with the results available in the literature. This accuracy is achieved with a proper particle resolution using a reduced total number of particles.

Flow around airfoil present difficulties to conduct high accuracy SPH simulation. It can be more difficult when it is thinner, particle distribution is irregular, or if it is moving. It is further challenging in complex motion applications like in a rotational motion, elliptic motion, and motions which combine heaving, oscillation, and pitching. One of the first applications of SPH to model flow around an airfoil was conducted by Shadloo et al. [66] and [67] in which they have shown a comparative study between the weakly compressible and incompressible SPH methods. However, their results have shown poor accuracy even when employing large number of particles at moderate Reynolds number. With the help of multi-resolution and particle shifting technique, Sun et al. [68] indicated an improved success of airfoil SPH simulation for Reynolds number up to $10,000$. Huang et al. [20] illustrated an improved simulation of viscous fluid flow around airfoils using SPH scheme possessing iterative particle shifting technology (IPST) and which does not use kernel gradients. The numerical study of airfoil motion is important for application like in unmanned aerial vehicles (UAVs) for thrust generation and in vertical axis wind turbines (VAWTs) for power optimization.

We discuss three stationary airfoil cases: (i) NACA5515 airfoil at $\alpha=5^{\circ}$ and $Re=420$, (ii) NACA0008 at $\alpha=4^{\circ}$ and $Re=2000$, (iii) NACA0008 at $\alpha=4^{\circ}$ and $Re=6000$. The computational domain is $[-3L,7L]\times[-3L,3L]$, where $L=c=1$ m is the characteristic length (chord) with center at the airfoil midpoint. The algorithm of Negi and Ramachandran [63] is used for packing the airfoils.

The flow around a moving square represents a simple geometry but a flow situation that can be very complex due to the sharp edges. The problem consists of a square geometry moving through a stationary fluid kept inside a rectangular box. The test can generate intense unsteady vorticity and has been a challenging test case for SPH solvers. This simulation is part of a validation test suite compiled by the SPH rEsearch and engineeRing International Community (SPHERIC) [70]. The benchmark results of this test case are simulated using Finite Difference Navier-Stoked (FDNS) solver by [71]. In SPH [43], and [72] studied this problem. We study the time histories for the drag force coefficients (pressure, viscous or total) and the generation and diffusion of vorticity.

The geometric configuration is shown in fig. 11(a) at time $t=0$. The body motion in time is shown in fig. 11(b). The domain size is $L_{y}\times L_{x}=5$ m $\times\ 10$ m, with the square obstacle having a length of 1m. The square motion is prescribed with a smooth acceleration in the $x$-axis starting from rest until it reaches a steady maximum velocity $U_{\text{square}}=1.0\ \text{m/s}$. The boundary conditions are no-slip and no-reflection on all walls of the rectangular domain and slip velocity on the solid square. The simulation is conducted for a Reynolds number of $100$ with no gravity force and zero initial pressure field for the viscous incompressible Newtonian fluid of density $1\,\text{kg m}^{-3}$.

In fig. 12, we show the time evolution of the drag coefficient of a moving square at (a) $Re=100$ using $D/\Delta{x}_{\min}=80$ and (b) at a slightly higher $Re$ of $600$ using $D/\Delta{x}_{\min}=200$. Simulated results obtained from the present method are compared with the results of the FDNS solver [71] and Marrone et al. [72]. The comparison shows, especially with FDNS, very good agreement with all the necessary features.

The simulation requires $44k$ fluid particles using adaptive resolution for the converged solution. The initial spacing of particles on the domain is $\Delta x_{\max}=0.04$; the smallest spacing is $\Delta x_{\min}=0.0125$ near the obstacle. The solid particles are packed using the minimum resolution to create a mapping in the fluid-solid interface. If the $\Delta x_{\min}$ was used in the entire domain, $D/\Delta x_{\min}$ would be $80$ which produces $320k$ as the total number of particles. However, with an adaptive particle resolution, we have used only $44k$ particles to produce accurate results. This substantially reduces the number of particles by $7.273$ times when compared with Marrone et al., [72] and by four times when compared with the results by [71].

Figure 13 shows the variation of drag coefficient for flow around a moving square at different Reynolds numbers. The comparison is made for three Reynolds numbers of $50$, $100$, and $150$. The results of the adaptive EDAC-SPH are in good agreement with [71] results. However, due to numerical reasons [71] show the small amplitude and high-frequency oscillations which are avoided in our case.

Figure 14 show the particle distribution with color indicating pressure, vorticity and velocity magnitude taken at $T=5$ and $T=8$. The qualitative observations illustrate the generation, diffusion, and convection of the unsteady vorticity proving the capability of the adaptive EDAC-SPH scheme.

Muta and Ramachandran [1] used a collection of particles called “background” particles to define the particle spacing $\Delta s$ and reference mass. But in this work, we show that the background particles are not necessary to set this reference mass. Updating the reference mass on the particles themselves produces results closely matching those produced by the original algorithm using background particles. This reduces memory usage by removing background particles’s storage, and a small reduction in the computational time. We perform two simulations to establish the no background particle approach. We first consider the flow past a circular cylinder at $Re=1000$, the simulation parameters and domain sizes are the same as defined in [1]. Furthermore, we compare the results with the established vortex method results [73, 74]. We then consider the flow past a moving square in a box at $Re=150$. Refer to section 5.3 for the simulation parameters. We compare the results to the incompressible FDNS simulation results [71]. In this case, the geometry is moving. In addition, we use solution adaptivity to refine the particles to the lowest resolution where the vorticity exceeds 5% of the maximum vorticity in the simulation.

Figure 15(a) shows the time history of the drag coefficient of the flow past cylinder problem without solution adaptivity. The results of adaptive EDAC-SPH are computed with and without use of background particles, and both results match very well with the reference data from vortex method simulations. Figure 15(b) shows the time history of the drag coefficient of the moving square problem using solution adaptivity. This further confirms a good match when not using the background particles. Figure 16(a) we show the distribution of the spacing $\Delta s$ without using background particles and in fig. 16(b) we show the vorticity distribution. It can be seen that the regions satisfying the solution adaptivity criteria are refined to the highest resolution, and the particles closest to the solid square particles are refined to the highest level.

The numerical examination of flow around an oscillating single cylinder is helpful for the simulation and deeper analysis of two or more oscillating cylinders. The study of flow over a cylinder is important for the study and design of flow over a car, flow over aircraft, flow over airfoil blades, or flow over a pillar for the design of bridges [64, 65, 75, 76, 77].

Bao et al., [65] classify the flow around the oscillatory motion of cylinders into lock-on and unlock-on states which are also characteristic of two cylinder flows. The classification depends on the shedding and oscillating frequencies for which the transition occurs around a frequency ratio of one.

We analyze the oscillatory motion of a single cylinder subject to cross-flow at $Re=100$. The sinusoidal motion of the cylinder is defined by $Y=A\sin(\omega t)$, where $A$ is the amplitude, $\omega=2\pi f$ the angular velocity, and $f$ the frequency of oscillation. Two amplitudes of oscillation (low and high) and frequency ratio $(f_{R})$ ranging from low to high are chosen for the simulation. Here $f_{R}=f/f_{0}$, where $f_{0}$ is a natural vortex shedding which is $0.165$ for a single cylinder oscillating at $Re=100$. The values are $A=(0.25,1.25)D$ and $f=(0.5,1.5)f_{0}$.

Figure 17 shows the time history of the fluctuating lift coefficient for four combinations of amplitude and oscillation frequency. The range includes lock-on and unlock-on states. In fig. 17(a) $(A/D,f_{R})=(0.25,0.5)$, the simulated results initially show faster transition and match with [65] capturing the peak-to-peak amplitude and all trends. In fig. 17(b) $(A/D,f_{R})=(0.25,1.1)$, our simulation depicts pure sinusoidal oscillation from the beginning, matching very well for $t\geq 80$ seconds. The fluctuating lift coefficient beats with a stronger magnitude for higher oscillating frequencies. Figure 17(c) $(A/D,f_{R})=(0.25,1.5)$, takes a combination of the smallest amplitude ratio with the highest value of the oscillating frequency, the results match very well from the beginning. In fig. 17(d) $(A/D,f_{R})=(1.25,1.5)$, both parameters of the amplitude and oscillating frequency are on the higher side. The results of this simulation capture all the trends and peaks of the lift coefficients fluctuation and are in good comparison with the reference data.

In this section, we discuss the parallel performance of the adaptive EDAC-SPH implemented using the framework of PySPH. To show that our adaptive algorithm is fully parallel, we run the moving square test case with the highest resolution of $\Delta x_{\min}=0.0075$ and the lowest resolution of $\Delta x_{\max}=0.2$, giving us 57k particles. As opposed to the gradual start in the moving square problem described in section 5.3 we move the solid with a uniform velocity of 1 ms^-1. We proceed to test this problem on a 6-Core Intel Core i7 CPU. We compare the time taken to run 1000 time-steps on a serial single-core single-thread execution with a parallel (OpenMP) execution on the multi-core CPU with 4-threads and 6-threads.

Table 3 shows the time taken by various algorithms described in this manuscript along with the time taken by the adaptive EDAC-SPH scheme. The splitting and merging component constitutes the Update Spacing algorithm 1, splitting algorithm 2, and merging algorithm 3; the removing and adding of particles constitutes for the time taken by the removal of the merged particles’ indices from the data structure, and addition of new particles to the data structure; the EDAC scheme constitutes the overall time taken by the SPH formulation.

Table 3 shows that the multi-thread implementation scales well in parallel with a scale-up of 3.13x on 4-threads and 3.73x on 6-threads. The time taken for removing and adding particles to the data structure constitutes less than 2.5% of the adaptive algorithm and this process is serial in execution. The present scheme again scales well with an increase in number of threads. The scale-up is 3.32x on 4-threads and 4.05x on 6-threads. The parallel execution of the adaptive algorithm constitutes 10–11% of the adaptive EDAC-SPH scheme.

The results of this section validate our formulation with the existing numerical simulations.

The acceleration of the center of mass of the ellipse, $a_{\text{cm}}(t)$, and the instantaneous pitch angle with respect to the horizontal, $\theta(t)$, are given by,

where $b_{0}=2.8209512$, $b_{1}=0.525652151$, $c=0.14142151$, $d=-2.55580905e-08$, and the pitch amplitude $\theta_{0}=10^{\circ}$ are the constants.

We simulate the translating and pitching ellipse in a rectangular domain of size $L_{y}\times L_{x}=5c$ $\times\ 10c$ at a Reynolds number of $550$. The ellipse is started with a smooth acceleration in the $x$-axis starting from rest until it reaches a steady maximum velocity of 1 m s^-1. The schematic of the domain is shown in fig. 18. The wall boundary condition is the same as that used in the moving square problem of section 5.3.

Figure 20 shows the vorticity distribution around the ellipse at different times. The top figure is at $T=1.3$, where the ellipse is $0^{\circ}$ from the horizontal. The middle figure is at $T=4.5$, where is ellipse is at $10^{\circ}$ from the horizontal and the bottom figure is at the final time of $T=8$.

Kinematics of a flapping motion describes a vertical translation (heaving), a horizontal oscillation, and pitching rotation. The amalgamation of these motions constructs an elliptical trajectory (path) of the pitching ellipse. This kind of kinematics is seen in the flying birds, aquatic animals, and has potential applications in high-efficiency Micro Air Vehicles (MAVs), and other energy harvesting vertical axis wind turbine (VAWT) designs [78, 79, 80, 81]. The pitching orientation and angle of attack can also vary according to the type of motion undertaken for high-level performance. The schematic diagram of ellipse motion is shown in 21.

The motion equations for the vertical plunge, $(h(t)$, and the horizontal sinusoidal motion, $w(t)$, of the ellipse are defined by,

where $h_{0}=2$, $\omega=2\pi/5$ is the angular velocity, and $b=0.5$ is a factor relating a horizontal motion with the vertical counterpart.

Figure 23 shows the vorticity plots of a plunging ellipse undergoing the mixed vertical-horizontal oscillatory motions forming an elliptical trajectory. The motion of the solid ellipse is described using the equations in 40, of a stationary fluid with boundary conditions discussed in 5.3. The schematic diagram of the problem is illustrates in fig 21 (a). In fig. 23(a) the ellipse is at the start of upward stroke with backward motion. In fig. 23(b) ellipse starts an upward stroke with a forward motion at period $T/4$. At period $T/2$, the ellipse starts the downward stroke with a forward motion. Similarly, figs. 23(c) and 23(d) show the motion of the ellipse at different time periods of $3T/2$, and $T$ respectively.

We consider a counter-clockwise rotating S-shaped body at $Re=2000$. The outer diameter $D$ is 2 m and fig. 24 shows the domain schematic. The minimum-resolution $D/\Delta x_{\min}$ is 10 and the maximum-resolution $D/\Delta x_{\min}$ is 200. The body is rotated at a frequency of $0.1$ sec^-1. We simulate the problem for $t=10$ sec. Figure 26 shows the vorticity distribution of the particles for the rotating S-shape at $t=2,5,8,10$ secs.