Boundary Deficiency Correction for Improved Boundary Conditions at Deformable Surfaces

Smoothed particle hydrodynamics (SPH) is a meshless, Lagrangian CFD method. SPH often utilizes static virtual particles to correct for integral deficiencies that occur near boundaries. These virtual particles, while useful in most cases, can be difficult to implement for objects which experience large deformations. As an alternative to virtual particles, a repulsive force algorithm is presented which loosely emulates the presence of virtual particles in the SPH momentum equation.

Smoothed particle hydrodynamics was fi rst developed for applications in astrodynamics [1][2].SPH was later adapted to fl uid dynamics due to its ability to effi ciently simulate complex free surface fl ows.It is based on the ability of a kernel function to approximate a fi eld value at a point by integrating over surrounding fi eld values.A kernel function is a fi nite area approximation of the Dirac delta function.Th e SPH approximation can be seen by fi rst considering the integration result of the product of a fi eld and the Dirac delta function (Eqn 1a).Th is integration exactly yields the fi eld value at the location of the Dirac delta function.By replacing the discontinuous Dirac delta with a kernel function possessing similar properties (but defi ned over a fi nite area) a useful approximation of fi eld functions can be obtained (Eqn 1b).Th is formula can then be discretized to obtain an equation applicable to Lagrangian particle dynamics (Eqn 1c).In the discretization, W ij is the value of the kernel function between a subject particle (i) and an infl uencing particle (j) with A j as the area or volume associated with the infl uencing particle [8].
Th e integration limits for Eqn 1a can be infi nite with identical results.However, because the integration product is zero everywhere except at x' = x, the integration limits can be reduced to encompass only the point x.Similarly, most kernel functions with compact support domain are defi ned as zero beyond a set radius.It follows that the integration/summation region for Eqns 1b and 1c need only encompass the region in which the kernel function is defi ned to be non-zero.
A classic kernel function, described by Monaghan and Lattanzio, is the Cubic-spline function (Fig. 1) [12].Th e Cubic-spline kernel function is defi ned as nonzero for all r < 2h (where h is a scaling parameter typically based on particle mass).It follows that the integration domain Ω of equation 1b is the region within a radius of 2h of the point x.Similarly, for equation 1c the support domain particles (Ω) consist of all particles within 2h of the subject particle i as illustrated in Fig. 2. Th e SPH method can be applied to fl uid dynamics in a number of ways.A common technique to adapt SPH to fl uid dynamics is Weakly Compressible Smoothed Particle Hydrodynamics (WCSPH).In WCSPH, fl uids typically considered nearly incompressible are allowed a small degree of compressibility.Th is yields a fi nite speed of sound The continuity density formulation instead calculates the velocity divergence for each particle.

Overview of Smoothed Particle Hydrodynamics
Because the velocity divergence of a fluid is equivalent to the time rate-of-change of density this value can be integrated over time to obtain density change.The corresponding particle discretized formula for the continuity density formulation is shown in equation 4 [8].
Particle acceleration is determined from the stress tensor for each particle .The stress tensor can be separated into translational (pressure) and viscous shear stresses.Shear stress can be approximated as the product of fluid's dynamic viscosity (µ) and strain rate tensor .For inviscid cases (µ=0), the particle discretized momentum equation reduces to a simple pressureforce summation.Equation 5 shows the particle discretized SPH momentum equation [8].
where for the fluid and changes the characteristic of the governing equations to allow for uncoupled solving of particle properties.Particle pressures are derived from an equation of state.The parameters of the equation of state determine the numerical speed of sound in that medium.The Tait Equation (Eqn 2) is a commonly used equation of state for water simulations [13].
Two common methods exist to calculate SPH particle density, these are the summation and continuity density formulations.The summation density formulation directly determines density by summing the kernel-weighted densities of all support domain particles (Eqn 3a).Substituting the definition of density (Eqn 3b) into Eqn 3a yields Eqn 3c [8] and tranforms the formula into a mass summation over an implied volume (the kernel function has units of inverse volume). (2) ) Assuming uniformly distributed particles in an unbounded region, the SPH particle discretized governing equations approach the Navier Stokes governing equations as smoothing length approaches zero.The unbounded assumption is required to satisfy that fluid field values exist at all points within a particle's support domain.When a particle's support domain has an insufficient number of particles to adequately approximate the fluid field values within its support radius it is said to be integral or boundary deficient.One way in which this can occur is if particle smoothing radius is set to a small value such that very few particles fall within a support radius, in this case the SPH governing equations will yield poor approximations of field functions.Because of this smoothing radius is typically empirically related to particle mass and density such that an acceptable influence radius is maintained.A second more complicated manner in which integral deficiency can occur is observable when a particle is within close proximity to a fluid boundary.This situation is shown in Figure 3  It is therefore desirable to provide an alternative to virtual particles for boundary deficiency correction of highly deformable objects.While various correction methods exist which achieve similar results, they are often computationally expensive.

SPH Boundary Behavior
Feldman and Bonet developed one such method to correct boundary deficiency for straight and corner boundaries by generating a curve fit to boundary deficiency accelerations [11].
While the goal is to achieve a boundary deficiency correction for arbitrary boundaries, it is advantageous to first consider the simple case of a 1-dimensional boundary.By observation of the boundary deficiency of a one dimensional fluid as in Fig. 3, it apparent that the boundary deficiency is similar in shape to the smoothing function.
Analyzing equation 5 acting at a boundary and assuming constant pressure, density, and mass, the momentum equation can be reduced as shown in equation 6.
Equation 6 shows that the acceleration boundary deficiency is proportional to the total area of the deficiency as well as fluid pressure.Because cases of interest typically involve non-constant pressures equation 6 is not strictly suitable for use as a boundary correction.However proportionality to the kernel function for boundary deficiency can stil be observed even in cases in which a pressure gradient is present.Fig. 5 illustrates one such case.Shown in Fig. 5 is the vertical acceleration of a unity density fluid influenced by a downward body acceleration of unity magnitude.
The pressure at the top of the water column is zero (gauge pressure), with pressure in the fluid varying as dp/dz=−ρg.The boundary deficiency in acceleration at the bottom of the water column still strongly correlates to the kernel function even in the presence of a pressure gradient.It is interesting to note that a boundary deficiency exists not only at the bottom boundary but also at the free surface.However, the low magnitude boundary deficiency at the free surface does not introduce substantial error as it typically results in only minor particle clumping.
While an exact acceleration boundary deficiency correction could be gained by discerning the pressure gradient normal to a surface and analyzing the geometry of the deficiency, such an approach would increase computational complexity.Instead a simple -if inexact-correction is suggested in which the pressure is assumed to be nearly constant over the scale of the smoothing length.Then by assuming a known boundary acceleration, the relative fluidboundary acceleration in the absence of boundary forces can be calculated.If this relative acceleration is assumed to be the result of a boundary defi ciency then it can be corrected by applying a repulsive acceleration distributed as CW ij away from the boundary.Where C is determined by choosing a sample particle, determining the relative particleboundary normal acceleration and dividing by W ij .To improve robustness it is advisable to average C over a small sample of nearboundary particles.Th is reduces correction error due to spurious pressure fl uctuations.Equation 7shows this one dimensional acceleration boundary defi ciency correction with j representing a boundary particle index.
where Th e assumption that a boundary acceleration is known is suitable for fi xed boundaries but requires an approximation when applied to freely moving boundaries.For objects with much greater density than the surrounding fl uid, forward interpolation of acceleration is acceptable as object acceleration will change only slowly relative to the simulated timescales when acted on by fl uid forces alone.However, as relative object-fl uid density decreases, this approximation worsens.Further work is necessary to assess the impact of this assumption on low density object dynamic behavior.
While the correction presented in equation 7 is suffi cient for one dimensional cases where a single boundary particle governs a boundary defi cient region, extension to higher dimensions requires a blending of corrections from multiple boundary particles.An empirical weighting function to blend correction values for two dimensional cases is shown in equation 8. Fig. 6 shows a visualization of the resultant weights.
(where Ω´ is the set of all boundary particles within r infl uence of particle i) (where 'spacing' is the local boundary particle spacing) Fig. 6 illustrates the weighting scheme applied to simple concave and convex boundaries.Th e weights of the three boundary particles are shown by the blue, green, and red shades respectively.Th e boundary defi ciency correction is applied to two test cases (one with a rectangular boundary, one with an elliptical boundary) each with constant pressure fl uid.Fig. 7 shows the rectangular tank case.( Similar boundary acceleration correction work performed by Feldman and Bonet [11] does not require an assumption of a known boundary acceleration, but instead determines an acceleration correction by calculating the fl uid pressure gradient and assuming a known boundary geometry (straight or straight-corner).Further work is required to compare the relative performance of the two methods when applied to various cases.
Virtual particles are normally used to correct errors due to integral defi ciencies that appear in the governing equations near boundaries.Virtual particle behavior for deformable objects can be diffi cult to implement due to particle clumping after deformation.A simple repulsive correction which loosely emulates the presence of virtual particles in the momentum equation has been derived.An empirical weighting function to extend the theoretical boundary correction to higher dimension cases has been presented.but tends to slightly overcorrect at sharp corners and undercorrect near regions with high curvature.

Fig
Fig. 2. SPH Particle Support Domain