A Consistent Multi-Resolution Smoothed Particle Hydrodynamics Method
We seek to accelerate and increase the size of simulations for fluid-structure interactions (FSI) by using multiple resolutions in the spatial discretization of the equations governing the time evolution of systems displaying two-way fluid-solid coupling. To this end, we propose a multi-resolution smoothed particle hydrodynamics (SPH) approach in which subdomains of different resolutions are directly coupled without any overlap region. The second-order consistent discretization of spatial differential operators is employed to ensure the accuracy of the proposed method. As SPH particles advect with the flow, a dynamic SPH particle refinement/coarsening is employed via splitting/merging to maintain a predefined multi-resolution configuration. Particle regularity is enforced via a particle-shifting technique to ensure accuracy and stability of the Lagrangian particle-based method embraced. The convergence, accuracy, and efficiency attributes of the new method are assessed by simulating four different flows. In this process, the numerical results are compared to the analytical, finite element, and consistent SPH single-resolution solutions. We anticipate that the proposed multi-resolution method will enlarge the class of SPH-tractable FSI applications.
keywords:smoothed particle hydrodynamics, multi-resolution, refinement and coarsening, fluid-structure interactions
SPH is a Lagrangian method that has been used to solve partial differential equations (PDEs) associated with mass, momentum, and energy conservation laws Monaghan_SPH_2005 . In SPH, these PDEs are spatially discretized employing a set of particles that possess material properties and interact with each other in a fashion controlled by a smoothing, or kernel, function with compact support. These particles move according to inter-particle interactions and external forces. Due to its Lagrangian nature, SPH enables efficient modeling of multiphase flows with moving interfaces Adami2010surface ; Tofighi2013numerical ; Tart2016JCP , complex fluids PanJG2012 ; PanJCP2013 ; SPHPolymer2016 ; SPHNoncolloid2016 , materials subject to large deformation Gray2001elastic ; Pan2013material ; Hu2016dynamic , fluid-structure interactions (FSI) problems Hu2014dynamic ; Schorgenhumer2013interaction ; Pazouki2014FSI , and various transport phenomena Pan_LLNS_2014 ; Tart2015CompGeo ; PanBMC2015 ; SPHProton2015 .
In most SPH simulations, the computational domain is populated with particles of uniform spacing, which leads to a single–resolution scenario. This is suitable when there are relatively small variations in the velocity/pressure gradient over the entire domain of analysis. Otherwise, it is desirable to discretize the domain with SPH particles of multiple resolutions. For instance, for a colloids-in-suspension FSI problem, the near field, i.e., next to each colloid, would call for high resolution; the far field would rely on a coarser resolution. Due to the Lagrangian nature of SPH, a locally refined region is expected to move along with the immersed body leading to an adaptive particle refinement (APR) process. In this regard, APR for SPH serves the same purpose as the adaptive mesh refinement does in the finite element and finite volume methods for FSI. The expectation is that APR will maintain the accuracy of the SPH solution while improving efficiency and increasing the size of the problem solved.
In this context, several multi-resolution SPH methods have been proposed, e.g., Lastiwka2009 ; Oger2006wedge ; Feldman2007refinement ; Lopez2013refinement ; Vacondio2013variable ; Barcarolo2014refinement ; Bian2015DDSPH . In these efforts, SPH particles were distributed according to a pre-defined multi-resolution configuration. What set these approaches apart was the coupling between regions of different resolutions and thereby the accuracy and continuity of the numerical solution at and near the interface. In Omidvar2012variable2D , the same large smoothing length was used for both the low- and high-resolution regions. This enhanced the accuracy attribute of the solution in some degree but led to a large number of neighbors for the particles in the high-resolution region. Thus, performance was hurt, which defeated APR’s purpose. In Vacondio2013variable ; Vacondio2016variable ; Barcarolo2014refinement , each region had its own smoothing length, which reduced the computational cost substantially. However, using an inconsistent SPH discretization, the numerical solution near the interface lacked accuracy, an artifact that became worse as the ratio of the two resolutions increased. This point will be demonstrated in §3.1. To tackle this issue, an approach promoting different-resolution subdomains connected via an overlap region was described in Bian2015DDSPH . Therein, the overlap region was employed to exchange state information and fluxes. At each time step, the solutions in the two subdomains were brought in sync via the overlap region by an iterative process. The authors investigated the sensitivity of numerical accuracy with respect to the thickness of the overlap region; a minimum thickness of the overlap region was determined numerically. On the upside, the approach can employ a different time step for each of the subdomains, hence reducing the computational effort in the low-resolution subdomain. On the downside, bringing the solutions of the two subdomains in agreement via an iterative algorithm remains costly and ad-hoc. Finally, if the two subdomains are moving dynamically, continually establishing an adaptive overlap region poses further difficulties.
Therefore, we propose a different multi-resolution SPH method in this work. In that, each region has its own smoothing length, and particles of different resolutions directly interact without any overlap region. The second-order consistent discretization is employed to discretize the governing equations via local correction matrices computed at each particle Fatehi2011error ; Traska2015second-order ; Pan2017Modeling . It thereby ensures the accuracy and continuity of the solution at the interface. Due to its higher-order accuracy and convergence, our multi-resolution SPH needs fewer number of particles to achieve the same accuracy and thereby further saves computational cost. In addition, we employ a dynamic refinement/coarsening of SPH particles to maintain the predefined multi-resolution configuration as SPH particles advect with the flow. In this respect, we follow in the steps of several other dynamic refinement/coarsening approaches Meglicki19933d ; Kitsionas2002splitting ; Kitsionas2007clump ; Feldman2007refinement ; Vacondio2012split ; Lopez2013refinement ; Vacondio2013variable ; Barcarolo2014refinement ; Vacondio2016variable . The basic underlying idea is that of splitting a “big” particle into several “small” ones according to a certain geometric arrangement centered at the location of the big particle and, conversely, merging nearest small particles into a big one at the center-of-mass of those small particles. To a large extent, what differentiates different approaches is how the splitting and merging occur, for example, upon splitting where the small particles are located and how their masses and kernel lengths are assigned.
Against this backdrop, there are several aspects in which the proposed method is different from the current dynamic refinement/coarsening techniques. First, when a big particle moves out of a low-resolution region, we split it into the same count of small particles according to the ratio of two resolutions. After splitting, the big particle vanishes. In Barcarolo2014refinement the big particle is persistent albeit virtual, and its motion through the high-resolution region is calculated from interpolated information provided by the small particles. When the big particle crosses back into the low-resolution region, it sheds its virtual status to become a SPH particle whose motion is dictated by the governing equations. By avoiding virtual particles, our approach requires much simpler code implementation, uses less memory, and eliminates the overhead of keeping tabs on virtual particles. Second, in Feldman2007refinement , a big particle was split into small particles with non-uniform masses but the same inter-particle spacing and smoothing length. Therein, the mass of each small particle was determined by minimizing the splitting-induced error in the density field, which only partially alleviated the discontinuity of the solution at the interface between regions of different resolutions. Alternatively, an approach for splitting a big particle into small particles with the same mass was proposed in Lopez2013refinement . Therein, the inter-particle spacing and smoothing length were determined by minimizing the splitting-induced error in the density field. In this vein, our approach falls in-between: all small particles have equal mass, spacing , and smoothing length in the splitting process. Moreover, we maintain the accuracy and continuity of the solution at the interface by means of: () the consistent SPH discretization and () interpolation-based splitting/merging process. In (), we use a second-order interpolation scheme to endow each small particle with its own velocity, with a similar interpolation-based approach used for merging. In this regard, our approach is more accurate than the averaging-based technique presented in Vacondio2013variable ; Vacondio2016variable . Finally, compared to these latter techniques, we always split or merge the same number of particles, which leads to the presence of only two families of SPH particles at the low-resolution/high-resolution interface.
Furthermore, in Lagrangian particle-based methods, particle regularity is critical to ensure numerical stability and accuracy of the simulation. However, highly distorted particle distribution may appear as SPH particles advect with flow. Thus, in this method, we employ a particle shifting technique to enforce a uniform particle distribution without sacrificing efficiency and Lagrangian nature. In that, a particle is slightly shifted away from its streamline according to a computed shifting vector, and after shifting hydrodynamic variables are corrected according to its new positions using the second-order interpolation. The multi-resolution configuration of particles is particularly considered in determining the shifting vector. Ultimately, the no-slip boundary condition for the fluid velocity is imposed through ghost solid particles with velocities linearly extrapolated from the fluid velocity TakedaSPHBC1994 ; Morris1997modeling ; Holmes2011Smooth , which ensures the second-order accuracy at the boundaries Macia2011SPHBCs .
This contribution is organized as follows. The next section provides a detailed account of the proposed new multi-resolution SPH method. To that end, it touches on several topics: standard SPH discretization, consistent SPH discretization, dynamic splitting and merging of SPH particles, particle-shifting technique, enforcement of the no-slip boundary condition, and time integration. In Section 3, we demonstrate the accuracy, convergence, and efficiency of the proposed method through numerical tests. We model four different flows – the transient Poiseuille flow, flow around a periodic array of cylinders, flow around a rotating cylinder near a moving wall, and a translating and rotating ellipsoid moving in fluid. The numerical results obtained are compared with analytical, finite element, or consistent SPH single-resolution solutions. We provide concluding remarks and directions of future work in Section 4.
2 Multi-resolution SPH method
2.1 Lagrangian hydrodynamic equations
We consider the fluid flow taking place in the space and the solid (wall or immersed body/bodies) occupying the space . The boundary separates and , i.e., . The time evolution of the fluid is governed by the continuity and Navier-Stokes (NS) equations
which hold for . Here, is the pressure; is the body force per unit mass acting on the fluid; is the kinematic viscosity of the fluid; and are fluid density and viscosity, respectively. In the weakly compressible SPH, the above equations are closed by an equation of state (EOS) for the pressure as:
where is the speed of sound. This EOS guarantees nonzero background pressure and hence no negative pressure in the fluid field Balsara1995StaSPH ; Swegle1995StaSPH ; Morris1996StaSPH . At the fluid-solid boundary, the no-slip boundary condition for the NS equation is imposed as for , where is the velocity of the boundary .
2.2 Spatial discretization
Herein, we rely on SPH for the spatial discretization of the equations of motion. The value of a function at , the position of particle , is then approximated as Monaghan_SPH_2005 :
Here, is the volume of a SPH particle and defined as:
where is the kernel function. We adopt a quintic spline function that has continuous and smooth first and second derivatives, and which is expressed as
Here, ; is the kernel length; ; ; and, for future reference, . The constant assumes the value , and for 1-, 2- and 3D problems, respectively. Only neighbor particles of particle in the set contribute to the summation in the interpolation scheme of Eq. (3), i.e., . According to this equation, the density at can be expressed as
with denoting the mass of particle .
In the proposed multi-resolution SPH method, the computational domain is divided into two regions: a high resolution one, which hosts SPH particles with a small , and a low resolution one in which the SPH particles have a large . As shown in Figure 1, all SPH particles are initially placed on a regular lattice with the spacing of in the high-resolution region and the spacing of in the low-resolution region. Different ratios of are allowed in this method. Here, , and .
In the standard SPH discretization, the gradient and Laplacian of the function are discretized as
two formulas used in the NS equation for the pressure and viscous terms, respectively. Previous studies have demonstrated that this discretization cannot guarantee the exact gradient for linear functions and Laplacian for parabolic functions and hence introduces inconsistent convergence in the numerical accuracy, i.e., the second-order accuracy of SPH can only be recovered by taking sufficiently large compact support of neighbors Quinlan2006truncation ; Basa2009robustness ; Fatehi2011error . We also noticed that for multi-resolution SPH, this discretization leads to inaccurate numerical solution near the interface. In fact, as shown in Section 3.1, the quality of the solution at the interface will deteriorate as the ratio of two resolutions increases. We alleviate this numerical artifact by using a consistent discretization in which the gradient and Laplacian are evaluated as
where “” represents the dyadic product of the two vectors, and “” represents the double dot product of two matrices. The correction matrices and guarantee the exact gradient for linear functions and Laplacian for parabolic functions regardless of the ratio of Fatehi2011error ; Traska2015second-order ; Pan2017Modeling . Hence, the consistent discretization of the NS equation is given as:
Here, both and are symmetric matrices in . Specifically, the component of the inverse matrix of is:
To define , we adopt the Einstein summation convention and seek a solution of the following Traska2015second-order :
where the third-order tensor is defined as:
Here, , , , and denote the th component of vectors , , , and , respectively. The details about how the components of and are computed in 2D are presented in Appendix.
Finally, we employ a renormalization step for calculating the density of SPH particles near the interface of two regions of different resolutions:
2.3 Splitting and merging of SPH particles
Due to the Lagrangian nature of SPH, particles advect with the flow. Thus, as the system evolves, particles from the high-resolution region might leave it to enter a low-resolution region, and vice-versa. This process is anchored by rules facilitating the dynamic splitting and merging of SPH particles. Specifically, when a big particle moves out of its original region, it is split into small particles following the procedure described in Algorithm 1.
Note that in the third step above, based on its new position, each new particle receives a small velocity and density correction via a second-order interpolation that draws on information associated with the original “big” particle:
where the field variable represents or , and is calculated according to Eq. (8). The second-order accuracy of this interpolation scheme matches that of the consistent SPH discretization (Eqs. (8)–(14)) described above. Figure 2 further illustrates this refinement procedure in 2D with ; i.e., splitting a big particle into small particles. The salient attribute of this approach is linear momentum preservation.
Conversely, when small particles move out of a high-resolution region, the nearest ones are grouped and merged into a big particle following the procedure described in Algorithm 2.
If the high/low resolution regions are predefined, the proposed splitting and merging techniques can maintain their geometry over time. In this regard, the techniques proposed are different from the approach in Vacondio2013variable in one key aspect. After a big particle is split into small particles, in Vacondio2013variable the velocity of each small particle is simply set to that of the big particle without any correction step. Merging occurs gradually with two nearest small particles merged at each time step. Consequently, particles of various sizes co-exist in the computational domain, making it difficult to control the numerical accuracy and to implement the periodic boundary condition.
2.4 Particle regularity
SPH particles’ advection can lead to scenarios characterized by high particle disorder and/or regions with high particle depletion/plenitude. We provision against such scenarios, which can undermine the accuracy and stability attributes of the numerical solution, by slightly shifting particles away from streamlines to enforce a uniform particle distribution. The particle shifting technique in incompressible SPH was promoted in Xu2009Accuracy and subsequently found to be effective in sustaining particle regularity and numerical stability in Traska2015second-order ; Xu2009Accuracy ; Pan2017Modeling . This shifting technique was also applied in a weakly compressible SPH (WCSPH) formulation Shadloo2012Robust , and subsequently in a multi-resolution WCSPH method Vacondio2013variable . We follow the latter approach, in which the shifting vector is determined as
where ; ; ; is an adjustable dimensionless parameter and set as in this work; is the maximum velocity of the system; is the time step of simulation. We introduce the particle’s mass in determining the shifting vector in Eq. (16) to maintain the anisotropic configuration of particles near the interface of multi-resolution regions. With this, at the end of each time step, the position of particle is shifted by
2.5 Imposition of no-slip boundary conditions
We enforce the no-slip boundary condition for the fluid velocity via several layers of ghost (fixed) particles (black) in the solid, as illustrated in Figure 4. As such, the SPH approximations of velocity and its spatial derivatives for the fluid particles (red) near the fluid-solid boundary have full support of the kernel contained in the domain (). As in previous studies TakedaSPHBC1994 ; Morris1997modeling , we perform a linear extrapolation of velocities to the ghost particles. If is a fluid particle and a solid ghost particle, the velocity of the latter is computed as
where and denote the closest perpendicular distances to the boundary for the fluid and solid particles, respectively; and, is the velocity at point of the solid boundary. As discussed in Macia2011SPHBCs , the linear extrapolation introduces a local error, which matches the second-order accuracy of the consistent SPH.
There are two caveats. First, to ensure numerical stability, following a recommendation made in Morris1997modeling , we cap the relative velocity as
with . Second, computing and for arbitrary complex geometries can be challenging. Thus, we approximate these distances as in Holmes2011Smooth . To that end, an indicator function is used to differentiate fluid and solid particles as:
The distances to the boundary are then approximated as
where for the quintic spline function used in this work the kernel compact support is . This approach allows the specification of an arbitrarily complex geometry by simply placing a lattice of particles over a domain and labeling particles as either fluid or solid.
2.6 Time integration scheme
Time integration is performed using a second-order, explicit predictor-corrector scheme Monaghan1989problem ; Monaghan1994simulating . If and denote the current and next time steps, respectively, a particle’s velocity and position at the intermediate time step () are first predicted as:
where is the total force per unit mass exerted on particle . The corresponding density and pressure are obtained via Eqs. (6) and (2), respectively. Using this intermediate-step information, is evaluated and subsequently used for correction purposes:
Finally, the particle’s velocity and position at are computed as
while the density and pressure are obtained via Eqs. (6) and (2), respectively. Insofar the step size is concerned, is constrained by the CFL condition Monaghan1992smoothed , magnitude of acceleration , and viscous dispersion and chosen such that
Despite working in a multi-resolution setup, no additional constraints were necessary on selecting . Moreover, given that or for the numerical experiments considered herein, the most stringent constraint was the viscous dispersion, i.e., .
3 Simulation results
Four different flows were used to assess the accuracy, convergence, and efficiency of the proposed multi-resolution SPH method. In this exercise, the numerical solutions were compared to the analytical solution, or numerical solution obtained using the finite element method (FEM), or numerical solution obtained using the consistent, single, and high-resolution SPH approach. The fluid was assumed to be water with and kinetic viscosity . In Sections 3.1 and 3.2, a different kinetic viscosity was used, i.e., , for comparing with results reported in literature. At , the fluid was at rest; the SPH particles were distributed on a regular lattice throughout the computational domain.
3.1 Transient Poiseuille flow
The 2D transient Poiseuille flow took place in a straight channel with two fixed walls oriented along the axis and positioned at and . The flow was subject to the no-slip boundary conditions at the walls and periodic boundary conditions at the remaining boundaries. The flow was driven by a body force acting along the axis.
The simulation domain was partitioned into three regions across the channel width: two high-resolution regions near the walls and one low-resolution region in the middle, see Figure 5. The width of the high-resolution and low-resolution regions is and , respectively. The low resolution was fixed as . In a series of three simulations, the high resolution was varied from to to produce ratios of 1, 2, and 4.
Figure 6 shows the transient velocity profiles across the channel width computed at different times using the standard SPH discretization of Eq. (7). The SPH results deteriorate as time goes on; and, the discrepancy becomes worse as increases. This suggests that two different resolutions cannot be directly coupled using the standard SPH discretization.
Figure 7 shows the transient velocity profiles obtained via the consistent SPH discretization of Eqs. (8)–(14). The proposed approach accurately predicted the velocity profiles for ratios of 1, 2, and 4. No numerical artifact was noticed. Thus, the proposed SPH discretization allows for coupling of different resolutions without using an overlap region.
We next probed the convergence attribute of the multi-resolution SPH method. To that end, we varied the number of particles at each and monitored the relative error of the steady-state velocity. Figure 8 summarizes our findings. The convergence of the multi-resolution SPH is affected by . For small , the second-order convergence of SPH can be sustained. As it increases, the anisotropy of particle configuration increases. Thus, the second-order convergence degrades to the first-order as assumes values close to 2 or higher. The larger is, the earlier the degradation of convergence occurs. This finding about the effect of anisotropic particle configuration on the convergence of SPH is consistent with results reported in the previous studies Fatehi2011error ; Traska2015second-order . Thus, on one hand, it is desirable to use larger as it leads to more computational savings. On the other hand, these savings may sacrifice accuracy. As an optimal cost–accuracy trade-off, we chose for all remaining simulations.
3.2 Flow around a periodic array of fixed cylinders
In this test, a cylinder of radius was positioned at the center of a square domain of length . The flow was driven by a body force along the coordinate, and subject to the no-slip boundary condition at the cylinder and periodic boundary conditions at the remaining boundaries. To examine the long-time particle regularity and numerical stability, the simulation was run for more than . This flow reached its steady-state at about .
We first examined the numerical accuracy of the SPH solution when and ; i.e., in the single resolution case. To that end, we sliced the domain vertically with imaginary lines at and at . Figure 9 compares the SPH results along lines 1 and 2 to those predicted by FEM Morris1997modeling . The agreement between the two numerical solutions is good, with a relative error of 3.5%. Figure 10 pertains to the long-term response. Subfigure (a) shows the particles distribution at , which demonstrates the particle regularity is well maintained in a long-time simulation. Subfigure (b) shows the equi-magnitude contours of steady-state velocity averaged after .
Next, we employed the multi-resolution SPH to simulate the same flow. As illustrated in Figure 11, the computational domain was divided into two regions. The inner, high-resolution, and circular region was located next to the cylinder and had a thickness of . The remaining fluid domain was represented in low resolution. The cylinder and fluid inside the circular region were discretized using ; the low-resolution representation used . For comparison, we also performed a single high-resolution simulation with .
Figure 12 presents the snapshots at for the single- and dual-resolution simulations, respectively. Shown are the particle distributions with color correlated to velocity magnitude. Figure 13 further compares the velocity fields in the single- and dual-resolution simulations via equi-magnitude contours of steady-state velocity averaged after . These figures confirm that the multi-resolution predictions agree very well with those of the single-resolution SPH in the entire fluid domain. Moreover, when comparing the total drag force exerted on the cylinder, the difference is less than 6%. Note that 2943 particles were used in the multi-resolution SPH versus 6400 particles in the single-resolution SPH, a 54% reduction in degrees of freedom. For each 10- (physical time) simulation, the computational time of the single-resolution solution was 1104 , but only 695 for the multi-resolution solution.
Finally, we investigated the effect of the thickness of high-resolution region. For that, we performed simulations with different thicknesses of high-resolution regions. The total drag on the cylinder was then computed and compared to that obtained from the single high-resolution simulation. We found that enlarging the high-resolution region with 17% more particles would lead to about 80% improvement on the accuracy and less than 20% increase in computational time.
3.3 Flow around a rotating cylinder near a moving wall
The computational domain in this test was a square box of length . The center of the cylinder was positioned at and ; its radius was . The cylinder was rotating counter-clockwise at a constant angular velocity . The bottom wall was moving horizontally along at a constant velocity of from left to right; the top wall was fixed. No-slip boundary conditions were imposed at the cylinder and walls; periodic boundary conditions were applied at the remaining boundaries. We divided the computational domain at into two regions: the bottom high-resolution region with ; the remaining fluid domain discretized at a lower resolution . Figure 14 shows the two regions in the initial configuration.
Figure 15 shows a snapshot of the particle distribution: subfigure (a) pertains to the single-resolution SPH; subfigure (b) pertains to the dual-resolution case. The particle color is correlated to the fluid velocity magnitude at . Figure 16 shows the equi-magnitude contour lines of steady-state velocity averaged after . Figure 17 further depicts the distribution of the averaged steady-state velocities at three locations: , , and . All data suggest very good agreements between single and dual-resolution results. Noteworthy, 21,000 particles were used in the dual-resolution simulation, less than a half of 42,400 particles used for the single-resolution solution. For each 10- (physical time) simulation, the computational time of the single-resolution solution was 4168 , but only 2758 for the multi-resolution solution.
Finally, Bian et al. have simulated this flow in Bian2015DDSPH using the domain-decomposition multi-resolution SPH (DD-SPH) method. The results obtained herein confirm their results. Unlike DD-SPH though, our approach directly coupled the two regions of different resolutions without the need for an additional overlap region. Hence, it allowed for a narrower high-resolution region that translates into a reduction in the overall degrees of freedom.
3.4 A moving colloid in the fluid
One salient attribute of SPH is its ability to model bodies moving freely in the fluid. This brings along the prospect of high-to-low-resolution interfaces that advect with the fluid, i.e, cases in which the high and low resolution regions change dynamically. In the problems studied thus far, the two regions of different resolutions were stationary. Against this backdrop, this section focuses on how the proposed multi-resolution SPH method handles a moving refined region. The vehicle for this study was a 2D ellipsoidal colloid moving in fluid. The high-resolution region moved and its geometry was determined based on the location of the moving colloid. This allowed for an accurately resolved flow in the vicinity of the colloid while saving computational cost via a coarser resolution far from the colloid.
The computational domain was a square box of length . The length of the ellipsoid’s minor axis was ; the major axis was . Initially, the center of the ellipsoid was located at and . The high-resolution region, for which , was determined as a moving circular region of radius with its center coincident with that of the moving colloid. The low-resolution region with was the remaining fluid domain. The initial setup of the two regions is shown in Figure 18. At , the ellipsoid was kinematically constrained to move downward as a rigid body with a constant acceleration of along the coordinate. As the ellipsoid reached its maximal translational velocity of , it moved at this constant velocity. Also at , an additional motion was imposed on the ellipsoid – it was rotated as a rigid body counter-clockwise at a constant angular acceleration of . Once the angular velocity reached , the ellipsoid proceeded to rotate at this constant angular velocity. The no-slip boundary condition was imposed at the boundary of the ellipsoid; periodic boundary conditions were applied at the remaining boundaries.
The dual-resolution solution is compared in terms of accuracy with the single high-resolution counterpart. Figures 19–21 present at three different times the distribution of SPH particles; their colors are correlated with their velocity magnitude. The quantitative agreement with the single high-resolution SPH insofar the velocity field is concerned indicates that the solution of multi-resolution SPH with a moving refined region is accurate. In addition, there were 16,727 particles used in the multi-resolution SPH, only about 2/5 of 41,600 particles used in the single resolution. For each 10- (physical time) simulation, the computational time of single-resolution solution is 605 , but only 295 for the multi-resolution solution.
We have presented a new consistent multi-resolution SPH method. It was applied in the context of four different tests: transient Poiseuille flow, flow around a periodic array of cylinders, flow around a rotating cylinder near a moving wall, and an ellipsoidal colloid that translates and rotates in the fluid. The numerical results were compared to analytical, FEM, and consistent single-resolution SPH solutions. We reported good accuracy, second-order convergence, and notable reductions in SPH particle counts. The study was conducted in 2D. Yet, the approach is directly applicable to the 3D case, which insofar efficiency gains are concerned, stands to benefit even more from a multi-resolution approach.
The cornerstone of the proposed method is the second-order consistent SPH discretization employed to ensure the accuracy of the numerical solution when regions of different resolutions are directly coupled. We found that the accuracy and order of convergence were affected by the ratio of the two subdomains’ resolutions (). The second-order convergence can be maintained for . For , our convergence study suggests a gradual decay in the convergence order. As an optimal cost–accuracy trade-off, we chose for all remaining simulations.
A main component of the proposed technique is its splitting and merging strategy for SPH particles leaving or entering the high-resolution region. Both for splitting and merging, the process calls for the state variables to be corrected via the second-order interpolation thus maintaining the second-order accuracy of the consistent SPH discretization. In this setup, the predefined multi-resolution configuration was maintained throughout the simulation despite the SPH particles advecting in and out of these high/low-resolution regions.
Another noteworthy step of the method is a particle-shifting technique employed to enforce particle regularity. The latter is essential in Lagrangian particle-based methods to ensure numerical stability and accuracy. The shifting vector was computed by accounting for the anisotropic particle distribution at the high-to-low-resolution interface. We slightly shifted particles away from streamlines with hydrodynamic variables corrected to account for their new positions using the second-order interpolation. This strategy avoided a particle-spacing distortion, which in turn kept error propagation in check and enhanced stability without compromising the computational efficiency and Lagrangian nature of the solution.
Finally, two additional points, both pertaining to boundary conditions, were important in the overall accuracy of the approach. First, we imposed no-slip boundary conditions via ghost solid particles. Their velocities were linearly extrapolated from the fluid velocity. The second-order accuracy of this approach matches that of the consistent SPH. Second, a “smoothed approximation” was used for calculating SPH particle distance to the boundary, an ingredient that comes into play when imposing no-slip conditions on complex boundaries.
Moving on to aspects related to the computational cost of the proposed method, we first noted that the time step in the explicit integration was determined by and the handling of the NS viscous term. The proposed multi-resolution method does not introduce any additional constraints for the time step. Second, as different resolutions lead to different kernel lengths, big particles have a larger compact support and hence interact with more neighbor particles near the low-to-high-resolution interface. It leads to the computation of more pair interactions. However, this computation accounts for less than of the computational time in all of our simulations. Third, splitting and merging of SPH particles introduce an unavoidable computational overhead. Upon profiling all our tests, we found that this overhead is less than of the total computational time. This suggests that additional tasks mandatory in a multi-resolution solution lead to an insignificant burden. By comparison, the reduction in number of SPH particles was significant in the 2D setup considered here, with the 3D scenario standing to witness even more dramatic reductions in particle counts.
Looking ahead, a 3D implementation represents a future incremental step that will likely pose limited challenges. The less obvious aspects that remain to be addressed is how to develop consistent second-order methods that vary the resolution continuously, and how to account for the presence of features such as a wall, cables, shells, etc., where a non-spherical kernel might yield an improvement in accuracy and/or efficiency gains.
This work was supported in part by the 111 China Project (B16003) of and the National Science Foundation of China grants 11290151 and 11221202. The research was made possible in part by the US National Science Foundation grant GOALI-CMMI 1362583.
- (1) J. J. Monaghan, Smoothed particle hydrodynamics, Reports on Progress in Physics 68 (2005) 1703–1759.
- (2) S. Adami, X. Hu, N. Adams, A new surface-tension formulation for multi-phase SPH using a reproducing divergence approximation, Journal of Computational Physics 229 (13) (2010) 5011–5021.
- (3) N. Tofighi, M. Yildiz, Numerical simulation of single droplet dynamics in three-phase flows using ISPH, Computers & Mathematics with Applications 66 (4) (2013) 525–536.
- (4) A. M. Tartakovsky, A. Panchenko, Pairwise force smoothed particle hydrodynamics model for multiphase flow: Surface tension and contact line dynamics, Journal of Computational Physics 305 (2016) 1119–1146.
- (5) W. Pan, A. M. Tartakovsky, J. J. Monaghan, Smoothed particle hydrodynamics model for ice sheet and ice shelf dynamics, Journal of Glaciology 58 (2012) 216–222.
- (6) W. Pan, A. M. Tartakovsky, J. J. Monaghan, Smoothed particle hydrodynamics non-Newtonian model for ice sheet and ice shelf dynamics, Journal of Computational Physics 242 (2013) 828–842.
- (7) X. Xua, P. Yu, A multiscale SPH method for simulating transient viscoelastic flows using bead-spring chain model, Journal of Non-Newtonian Fluid Mechanics 229 (2016) 27–42.
- (8) A. Vázquez-Quesada, M. Ellero, Rheology and microstructure of non-colloidal suspensions under shear studied with smoothed particle hydrodynamics, Journal of Non-Newtonian Fluid Mechanics 233 (2016) 37–47.
- (9) J. Gray, J. Monaghan, R. Swift, SPH elastic dynamics, Computer Methods in Applied Mechanics and Engineering 190 (49) (2001) 6641–6662.
- (10) W. Pan, D. Li, A. M. Tartakovsky, S. Ahzi, M. Khraisheh, M. Khaleel, A new smoothed particle hydrodynamics non-newtonian model for friction stir welding: Process modeling and simulation of microstructure evolution in a magnesium alloy, International Journal of Plasticity 48 (2013) 189–204.
- (11) W. Hu, Q. Tian, H. Hu, Dynamic fracture simulation of flexible multibody systems via coupled finite elements of ANCF and particles of SPH, Nonlinear Dynamics 84 (4) (2016) 2447–2465.
- (12) W. Hu, Q. Tian, H. Hu, Dynamic simulation of liquid-filled flexible multibody systems via absolute nodal coordinate formulation and SPH method, Nonlinear Dynamics 75 (4) (2014) 653–671.
- (13) M. Schörgenhumer, P. G. Gruber, J. Gerstmayr, Interaction of flexible multibody systems with fluids analyzed by means of smoothed particle hydrodynamics, Multibody System Dynamics 30 (1) (2013) 53–76.
- (14) A. Pazouki, R. Serban, D. Negrut, A high performance computing approach to the simulation of fluid-solid interaction problems with rigid and flexible components, Archive of Mechanical Engineering 61 (2) (2014) 227–251.
- (15) J. Kordilla, W. Pan, A. Tartakovsky, Smoothed particle hydrodynamics model for Landau-Lifshitz-Navier-Stokes and advection-diffusion equations, The Journal of Chemical Physics 141 (22) (2014) 224112.
- (16) A. M. Tartakovsky, N. Trask, K. Pan, B. Jones, W. Pan, J. R. Williams, Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media, Computational Geosciences 20 (4) (2016) 807–834.
- (17) W. Pan, M. Daily, N. A. Baker, Numerical calculation of protein-ligand binding rates through solution of the Smoluchowski equation using smoothed particle hydrodynamics, BMC Biophysics 8 (2015) 7–18.
- (18) S. Liu, J. Savage, G. A. Voth, Mesoscale study of proton transport in proton exchange membranes: Role of morphology, The Journal of Physical Chemistry C 119 (4) (2015) 1753–1762.
- (19) M. Lastiwka, M. Basa, N. J. Quinlan, Permeable and non-reflecting boundary conditions in SPH, International Journal for Numerical Methods in Fluids 61 (7) (2009) 709–724.
- (20) G. Oger, M. Doring, B. Alessandrini, P. Ferrant, Two-dimensional SPH simulations of wedge water entries, Journal of Computational Physics 213 (2) (2006) 803–822.
- (21) J. Feldman, J. Bonet, Dynamic refinement and boundary contact forces in SPH with applications in fluid flow problems, International Journal for Numerical Methods in Engineering 72 (3) (2007) 295–324.
- (22) Y. Reyes López, D. Roose, C. Recarey Morfa, Dynamic particle refinement in SPH: Application to free surface flow and non-cohesive soil simulations, Computational Mechanics 51 (5) (2013) 731–741.
- (23) R. Vacondio, B. D. Rogers, P. K. Stansby, P. Mignosa, J. Feldman, Variable resolution for SPH: A dynamic particle coalescing and splitting scheme, Computer Methods in Applied Mechanics and Engineering 256 (2013) 132–148.
- (24) D. A. Barcarolo, D. Le Touzé, G. Oger, F. De Vuyst, Adaptive particle refinement and derefinement applied to the smoothed particle hydrodynamics method, Journal of Computational Physics 273 (2014) 640–657.
- (25) X. Bian, Z. Li, G. E. Karniadakis, Multi-resolution flow simulations by smoothed particle hydrodynamics via domain decomposition, Journal of Computational Physics 297 (2015) 132–155.
- (26) P. Omidvar, P. K. Stansby, B. D. Rogers, Wave body interaction in 2D using smoothed particle hydrodynamics (SPH) with variable particle mass, International Journal for Numerical Methods in Fluids 68 (6) (2012) 686–705.
- (27) R. Vacondio, B. Rogers, P. Stansby, P. Mignosa, Variable resolution for SPH in three dimensions: Towards optimal splitting and coalescing for dynamic adaptivity, Computer Methods in Applied Mechanics and Engineering 300 (2016) 442–460.
- (28) R. Fatehi, M. Manzari, Error estimation in smoothed particle hydrodynamics and a new scheme for second derivatives, Computers & Mathematics with Applications 61 (2011) 482–498.
- (29) N. Traska, M. Maxeya, K. Kimb, M. Perego, M. L. Parks, K. Yang, J. Xu, A scalable consistent second-order SPH solver for unsteady low Reynolds number flows, Computer Methods in Applied Mechanics and Engineering 289 (2015) 155–178.
- (30) W. Pan, K. Kim, M. Perego, A. M. Tartakovsky, M. L. Parks, Modeling electrokinetic flows by consistent implicit incompressible smoothed particle hydrodynamics, Journal of Computational Physics 334 (2017) 125–144.
- (31) Z. Meglicki, D. Wickramasinghe, G. V. Bicknell, 3D structure of truncated accretion discs in close binaries, Monthly Notices of the Royal Astronomical Society 264 (3) (1993) 691–704.
- (32) S. Kitsionas, A. P. Whitworth, Smoothed particle hydrodynamics with particle splitting, applied to self-gravitating collapse, Monthly Notices of the Royal Astronomical Society 330 (1) (2002) 129–136.
- (33) S. Kitsionas, A. P. Whitworth, High-resolution simulations of clump–clump collisions using SPH with particle splitting, Monthly Notices of the Royal Astronomical Society 378 (2) (2007) 507–524.
- (34) R. Vacondio, B. Rogers, P. Stansby, Accurate particle splitting for smoothed particle hydrodynamics in shallow water with shock capturing, International Journal for Numerical Methods in Fluids 69 (8) (2012) 1377–1410.
- (35) H. Takeda, S. M. Miyama, M. Sekiya, Numerical simulation of viscous flow by smoothed particle hydrodynamics, Progress of Theoretical Physics 92 (5) (1994) 939–960.
- (36) J. P. Morris, P. J. Fox, Y. Zhu, Modeling low Reynolds number incompressible flows using SPH, Journal of Computational Physics 136 (1) (1997) 214–226.
- (37) D. W. Holmes, J. R. Williams, P. Tilke, Smooth particle hydrodynamics simulations of low Reynolds number flows through porous media, International Journal for Numerical and Analytical Methods in Geomechanics 35 (4) (2011) 419–437.
- (38) F. Maciá, M. Antuono, L. Gonzalez, A. Colagrossi, Theoretical analysis of the no-slip boundary condition enforcement in SPH methods, Progress of Theoretical Physics 125 (6) (2011) 1091–1121.
- (39) D. S. Balsara, Von neumann stability analysis of smoothed particle hydrodynamics-Suggestions for optimal algorithms, Journal of Computational Physics 121 (1995) 357–372.
- (40) J. W. Swegle, D. L. Hicks, S. W. Attaway, Smoothed particle hydrodynamics stability analysis, Journal of Computational Physics 116 (1995) 123–134.
- (41) J. P. Morris, A study of the stability properties of smooth particle hydrodynamics, Publications of the Astronomical Society of Australia 13 (1996) 97–102.
- (42) N. J. Quinlan, M. Basa, M. Lastiwka, Truncation error in mesh-free particle methods, International Journal for Numerical Methods in Engineering 66 (13) (2006) 2064–2085.
- (43) M. Basa, N. J. Quinlan, M. Lastiwka, Robustness and accuracy of SPH formulations for viscous flow, International Journal for Numerical Methods in Fluids 60 (10) (2009) 1127–1148.
- (44) R. Xu, P. Stansby, D. Laurence, Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach, Journal of Computational Physics 228 (18) (2009) 6703–6725.
- (45) M. S. Shadloo, A. Zainali, M. Yildiz, A. Suleman, A robust weakly compressible SPH method and its comparison with an incompressible SPH, International Journal for Numerical Methods in Engineering 89 (8) (2012) 939–956.
- (46) J. J. Monaghan, On the problem of penetration in particle methods, Journal of Computational Physics 82 (1) (1989) 1–15.
- (47) J. J. Monaghan, Simulating free surface flows with SPH, Journal of Computational Physics 110 (2) (1994) 399–406.
- (48) J. J. Monaghan, Smoothed particle hydrodynamics, Annual Review of Astronomy and Astrophysics 30 (1) (1992) 543–574.
Appendix: Correction Matrices G and L
The symmetric gradient correction matrix of particle can be expressed in 2D as:
According to Eq. (11), the inverse matrix of can be determined by:
Thus, each component of is computed by inverting the above matrix.
The symmetric Laplacian correction matrix of particle can be expressed in 2D as:
Rewrite Eq. (12) as:
where the matrix is expressed as:
The components of this matrix are defined as:
Here, the components of the symmetric third-order tensor are determined by:
Given and , , , and can be obtained by solving Eq. (21).
Note that computing and for each particle requires no additional neighbor information than the standard SPH operators.