Coherent Vorticity and Discontinuous Flow in Particle-Based Sph Modeling

Smoothed particle hydrodynamics (SPH) is a particle-based method used for computational fluid modeling and simulations. The method uses a numerical representation of the Lagrangian form of the Navier-Stokes equations. As reflected by its name, it relies on smoothing operations to update the fluid particles' properties as time progresses through the simulation run. We found that many of the sequences to be smoothed contain numerous discontinuities, something that causes reduced accuracy if the discontinuities are not taken into consideration. To address this, we present a method for, and the quantified effects of, segmentwise smoothing on particle movement, with particular focus on:

- Identification of discontinuities;

- Vorticity analysis is used to identify particles with discontinuities and coherent vorticity structure;

- Smoothing of sequences containing discontinuities by segmentwise smoothing;

- The effect of segmentwise smoothing on particle movement has been quantified.

The SPH method was introduced in the late 1970s and early 1980s by among others [13, 15, 17], as a particle-based model used in astrophysics, and soon developed into applications for a variety of scientific fields as its flexibility and computational efficiency compared to grid-based models demonstrated improvements in certain cases. Both the gridbased and particle-based methods are still highly relevant today, depending on the specific applications and areas of interest. As these methods are widely used and well presented by numerous publications, we refer the reader to the already existing literature on those topics and limit Section 2 to deal with discontinuity, segmentwise smoothing, and wavelets as they apply to our research presented here. In Section 3, we present the results obtained by the described methods applied to the canonical sloshing tank benchmark case and the quantified effects of these methods on the flow of matter. The discussion section presents the main conclusions of our work and related ideas for future work in this area to benefit the SPH model.

As reflected by its name, smoothing operations in SPH are central parts of the method. The smoothing technique is used to update the particles' properties and propel the flow of the matter at each time step of a simulation run. Smoothing takes place in several of the model's equations. We here focus on the acceleration equation, also referred to as the momentum equation. There are several forms of the acceleration equation used for the numerical implementation of SPH. We used the equation implemented in the open software DualSPHysics [6] of the following form:

where $\mathbf{v}$ is velocity, $\mathbf{g}$ is the external gravity, $m$ is mass, $p$ is pressure, $\rho$ is density, and subscripts $i$ and $j$ refer to particle number. The $\nabla$ denotes the gradient, and $W_{i j}$ is the smoothing function, also referred to as kernel function, as a function of the distance between particles $i$ and $j, j$ representing all of $i$ 's neighbor particles located within the kernel's compactly supported domain $h$. The term $\prod_{i j}$ is the artificial viscosity term introduced in ref. [16].

The two coefficients, $0 \leq \alpha \leq 1$ and $0 \leq \beta \leq 1$, scale the bulk viscosity and the Von Neumann Richtmyer viscosity, respectively. ${ }^{-} c_{i j}$ and ${ }^{-} \rho_{i j}$ are the average speed of sound and average density, respectively, for particles $i$ and $j$. DualSPHysics uses $\beta=0$. Furthermore, $\mathbf{v}_{i j}$ and $\mathbf{r}_{i j}$ are the differences in velocity and distance between particles $i$ and $j$, respectively. With this form of the acceleration equation, the sequence to be smoothed for each fluid particle $i, i=0$, $\ldots, n-1$ at each time step is as follows:

Discontinuities in a function or a sequence can be identified by the use of wavelets as presented in refs. [2, 7] and further demonstrated in ref. [4]. That work did find the first-generation wavelet coiflet to be the most accurate in identifying both the number of discontinuities and their location compared to Daubechies, least symmetric, and symmlets. By performing a discrete wavelet transform (DWT) and applying a thresholding mechanism to the obtained DWT coefficients of the one-dimensional sequence $\left\{s_{(p \rho) \mathrm{j}}\right\}_{j \in h}$, the resulting coefficients after thresholding can be used to identify discontinuities. The method was successfully tested on, among others, a shock wave described by the researchers in ref. [5] and tested in ref. [4], as well as on the actual field data. Several methods for thresholding have been presented, depending on the particular application in question. DWT and thresholding of the resulting coefficients have unique and wide applications in signal and image denoising of white noise that cannot be successfully filtered by frequency-specific filters. We found the false discovery rate (FDR) thresholding mechanism [1] using a confidence level of 0.025 to be the most accurate for the identification of discontinuities. The FDR threshold is expressed as follows:

where $t_k^*$ is the threshold level, and $\left|z_{\frac{\alpha}{2} * \frac{k}{n / 2}}\right|$ is the right tail of a t-distribution with $\left(\frac{k}{n / 2}-1\right)$ degrees of freedom for a sequence of length $n / 2, n$ being the length of the one-dimensional smoothing sequence. In general, a first-generation wavelet requires equally spaced samples and a dyadic number of observations. The smoothing sequences represented by eqn. (3) were not equally spaced as the distance between neighbor particles varied, and the sequence lengths would only be a dyadic number by coincidence. We tested the coiflet's ability to accurately detect discontinuities under such conditions by resampling the benchmark function Segmentwise Cosine [4] and found that the ability to detect discontinuities did not degrade until less than $55 \%$ of the original 512 samples were removed. Most software packages for wavelet transformation offer a selection of techniques to deal with non-dyadic numbers of observations. We used the periodization option from Python PyWavelets that repeats the sequence's last-term value until a dyadic number is achieved. For the smoothing sequences analyzed in this work, any discontinuity at the far end of a sequence would have an insignificant influence on the results due to the low weights of a kernel function close to the end of the compactly supported domain, and hence any effect of this periodization extension could be ignored. The discontinuities identified were confirmed by visual inspection of the sequences as presented in Section 3.1.

For the vorticity analysis, we used a second-generation wavelet for the DWT of the twodimensional vorticity field. Second-generation wavelets differ from the first-generation wavelets in the multi-resolution structure. The second-generation wavelets' basis functions are not dilatation and transformation of the one and same mother wavelet, and thus these functions adjust to the non-equal spacing. The fulfillment of the multiresolution structure requirements [23] is commonly achieved through the biorthogonality requirement of the basis functions. The vorticity field in the sloshing tank case is highly irregular in terms of rugged edges along part of the surface of the flow and multiple smaller sections inside the field containing no fluid particles. As presented by multiple publications [22,23] and demonstrated in ref. [3], the second-generation wavelets are capable of handling such irregularities. Wavelets have demonstrated their unique ability to identify features in time and space, among others, by revealing small scale properties of significant importance that have not been otherwise possible, such as the continuous wavelet methodology described by [25] and successfully demonstrated by [11, 20]. More specifically, first-generation wavelets have been used for vorticity analysis in grid-based methods of fluid modeling where the sample spacing and the number of samples can be ensured to be regular and dyadic by the generation of the grid. Earlier works [8-10,18] analyzed the nature of vorticity in terms of coherent and incoherent components. They describe the coherent structure as being highly correlated, while the incoherent has the nature of white noise. We build our analysis of the vorticity on those conditions using second generation wavelets to differentiate between its coherent and incoherent structure. We performed a DWT of the two-dimensional vorticity field of the flow at different time steps. The DWT coefficients were subjected to thresholding. Specifically, the DWT coefficients with an absolute value higher than the threshold were taken to represent the coherent structure of the vorticity. The threshold applied was different from the FDR thresholding mechanism used for discontinuity identification. It is instead based on the idea that white noise has a constant variance. The thresholding formula is expressed in eqn. (5). In practice, this variance is unknown and was estimated by iterating until it converged to a constant value. In our case, convergence was typically achieved within four to seven iterations, a result that compared well to the analysis of the grid-based methods cited as less than 10 iterations [21]. The iteration was performed as follows: The first iteration calculated the variance of all DWT coefficients, and this variance was used in eqn. (5) for thresholding. The variance for the DWT coefficients of absolute value less than this threshold was then calculated and used for thresholding of those DWT coefficients. The iteration was repeated until the estimated variance showed a negligible difference from the one of the previous iteration. The last variance value was used for thresholding. The threshold level based on the first variance estimate is for vorticity analysis, commonly referred to as a high cut, while the threshold level based on the latest variance estimate of the iteration is referred to as a low cut.

where $\sigma_{c D}^2$ is the DWT coefficients' estimated variance, and $n / 2$ is the number of coefficients. The second-generation wavelets used in this study were generated from the symmlet 3 wavelet using the MATLAB MathWorks [14] for two-dimensional lifting scheme. The results in terms of the number of significant DWT coefficients were compared to running the same analysis with second-generation wavelets generated from a biorthogonal bior 1.3 wavelet. This resulted in the same number of significant DWT coefficients. For comparison reasons, both the high cut and the low cut threshold were tested in ref. [3], and our findings compared well to the grid-based results of ref. [18].

We built on the segmentwise smoothing method presented by the researchers in ref. [2] where the sequence or function to be smoothed gets divided into segments based on the locations of discontinuities in the sequence and these segments get smoothed individually and independently of each other. Experience from earlier work [4] demonstrated this method's ability to accurately preserve the discontinuities and perform recreation without the inaccuracy commonly represented by the Gibbs phenomenon when discontinuities are not taken into consideration. The fact that each segment is smoothed independently leads, in the SPH case, to only the first segment of a fluid particle's sequence being used for smoothing. This is expressed as follows:

where $j$ represents the first segment's neighboring particles, and $W_{i j, 1}$ is the kernel function with a compactly supported domain, also referred to as smoothing length, $h_1$, and scaled to preserve unity. Variable smoothing lengths and scaling are commonly used in SPH for various reasons, including dealing with cut-offs of the compactly supported domain by edges, controlling the number of neighbor particles both related to cut-offs, and also finally in relation to the more accurate modeling of sections of flow of less smooth nature [12,19,24].

When performing function recreation or smoothing of a function containing discontinuities, the accuracy of the result in the vicinity of the discontinuity is diminished. This is a welldescribed phenomenon in science, often referred to as the Gibbs phenomenon, and it is often observed as oscillations around the location of the discontinuity. Researchers in ref. [5] used the pressure function of a shock wave to research alternative SPH kernels and their ability to improve accuracy at the discontinuity. More accurate results were achieved in ref. [2] by the identification of a function's discontinuities and then the creation of the function's smoother segments between the discontinuities. Researchers in ref. [4] used this two-step method on the shock wave function and achieved a recreation of high accuracy [4]. Here, we present the identification of discontinuities. The results of segmentwise smoothing on particle movement are presented in Section 3.3. Our study focused on a two-dimensional sloshing tank model simulated by the SPH method using the open-source software DualSPHysics [6]. The tank dimensions and fill level were determined according to DualSPHysics' sloshing motion case (refer to Fig. 1).

Most numerical implementations of the acceleration equation, including DualSPHysics, use an artificial viscosity term to avoid the kernel's Laplacian, so the smoothing sequence of interest for this research is ref. (3). This sequence was calculated for every fluid particle $i=0,1, \ldots,(n-1)$ at each time step and smoothed to obtain an estimate of the particles' properties based on the particles' neighbor, $j$ located within the area of a circle of radius $h$ around particle $i$. We used the first-generation coiflet wavelet with two vanishing moments, hence six coefficients, and performed a DWT of the sequence. The DWT coefficients were thresholded, and the resulting significant DWT coefficients were taken to represent discontinuities in the sequence. The FDR thresholding criteria [1] with a confidence level of 0.025 was used to identify the DWT coefficients that represented discontinuities, referred to as significant DWT coefficients. Both the selection of wavelet and the threshold mechanism were based on their superior performance as described in ref. [4]. Examples of sequences with discontinuities and threshold levels are given in Figs. 2 and 3. We found that typical between $24 \%$ and $27 \%$ of the fluid particles had sequences containing at least one discontinuity. Contrary to our expectations, these numbers were seemingly independent of time steps and, as such, seemingly independent of the phase of the flow. This can primarily be explained by the fact that the majority of discontinuities and turbulence were located in the lower section of the flow and is discussed in more detail in Section 3.3. When increasing or decreasing the threshold level, the percentage of discontinuities varied somewhat, but with little or no variation from time step to time step. In an attempt to characterize or differentiate between discontinuities of potentially higher importance to the flow modeling, the fluid particles' pressure, velocity, and vorticity values were analyzed without being able to find such differentiation criteria. However, a more detailed analysis of vorticity by differentiating between the coherent and incoherent components of the vorticity did provide such criteria. The results of our vorticity analysis are presented in the following subsections.

The vorticity analysis was carried out by the use of second generation wavelets as described in ref. [3]. Coherent vorticity was observed both in the earlier phases where flow was seemingly smooth and throughout the simulation run. The initial time steps demonstrated relatively low coherent structure at the ends of the tank, in particular at the outbound end, as can be observed from the graph Significant DWT Coefficients in Fig. 4 that shows the location and magnitude of significant DWT coefficients. As time progressed through the simulation, both vorticity values and significant DWT coefficients displayed the increased absolute values at both ends of the tank, simultaneously, refer the graph Significant DWT Coefficients in Fig. 5. The second-generation wavelet transfer captured approximately $60 \%$ of the fluid particles at each time step. The majority of the remaining $40 \%$ were not captured due to the sparse particle representation in the flow field along the surface, especially inside the flow field observed as white areas in Fig. 6. Such white areas were observed even at the earlier stages that represent relatively smooth flow and correspond to the fact that particles are highly unevenly distributed in the flow field during the flow. The uneven particle distribution has been a topic of discussion from the introduction of the SPH method, as this method represents the matter solely through particles. From a physics point of view, it is due to the difference in the external force influencing each individual particle at a given time, a force that for water in particular is significantly larger than the inertia.

While ideally, we would like to capture a higher percentage of the fluid particles, it may not necessarily have a significant impact where those empty areas border particles with mostly incoherent vorticity. A visual inspection of the thresholded coefficients from Fig. 6b showed both coherent and incoherent vorticity located next to white areas as Fig. 7 shows an example of without further quantification in this work. The effort to, and effect of, capturing a larger percentage of the fluid particles is something we aim at as part of our future work. The percentage of significant DWT coefficients identified by our second-generation wavelet method compares well to the work on grid-based methods cited earlier.

Fluid particles with coherent vorticity and discontinuous smoothing sequences were subjected to segmentwise smoothing, and the results were compared to smoothing without taking the discontinuity into consideration. This selection was based on the fact that coherent vorticity represents correlated structures highly representative of the energy in the flow according to ref. [9]. We found that segmentwise smoothing changed the particle's movement by more than $100 \%$ in the most extreme cases and more typically by about $50 \%$ or less, as related to the change in position the particle would make from one time frame to the next. A relative higher proportion of particles with discontinuities and coherent vorticity were located toward the bottom of the tank, as may be observed from the Significant DWT Coefficients plots in Figs 4 and 5. This was also where the larger change in particle movement as a result of segmentwise smoothing was observed. Figure 8a illustrates the density distributions for the change in particle movement of particles located toward the bottom of the tank, and Fig. 9a illustrates the change in particle movement for the upper part of the flow field. In particular, when comparing the tails of the distributions for the two sections, refer Figs. 8b and 9b, the relative number of particles with a larger change in movement is higher in the lower section of the flow. Our work here has demonstrated the impact on particle movement from one time step to the next when smoothing is performed by taking discontinuities into consideration. This does not reflect the total effect on particles' movement with discontinuous smoothing sequences and coherent vorticity over multiple consecutive time steps. With the methods and results presented here, the evaluation of such effects is part of our future work.

We have here presented the methods for the identification and handling of discontinuities for smoothing operations in SPH modeling and simulations. A sloshing tank case that represents the uneven flow of matter was used. First-generation wavelets proved suitable for the discontinuities identification. Second-generation wavelets served as an efficient tool for vorticity analysis as presented in ref. [3]. The results of this work confirmed quantifiable differences in terms of particle movement when discontinuities were taken into consideration and smoothed according to our segmentwise smoothing method. The observed effects were for individual time steps, and as such, were partially representative of the accumulative effect from multiple time steps. This is particularly true for fluid particles with discontinuities and coherent vorticity structures over consecutive time steps. The plans for our future work further include the following:

- Implementation of discontinuity identification and segmentwise smoothing into a numerical model for particle-based flow simulations.

- Evaluation of alternative threshold levels for vorticity analysis to decide the most optimal tradeoff between computational efficiency and discontinuities with an effect on smoothing, to see if a threshold higher than the low cut used here may be good enough.

Additional future research presented in ref. [3] will be complementary to the work presented here for the development of our numerical model for particle-based flow simulations.