Pablo Eleazar Merino-Alonso, Fabricio Macià, Antonio Souto-Iglesias

1. Mathematical Modeling, Analysis, and Simulation Applied to Engineering Research Group (M2ASAI), Universidad politécnica de Madrid (UPM), School of marine engineering (ETSIN), Madrid, Spain

2. Marine and Hydrodynamical Model Basin Research Group (CEHINAV), Universidad Politécnica de Madrid(UPM), School of Marine Engineering (ETSIN), Madrid, Spain

Abstract: The aim of this work is to study the solution of the smoothed particle hydrodynamics (SPH) discrete formulation of the hydrostatic problem with a free surface. This problem, in which no time dependency is considered, takes the form of a system of linear equations. In particular, the problem in one dimension is addressed. The focus is set on the convergence when both the particle spacing and the smoothing length tend to zero by keeping constant their ratio. Values of this ratio of the order of one, corresponding to a limited number of neighbors, are of practical interest. First, the problem in which each particle has one single neighbor at each side is studied.The explicit expressions of the numerical solution and the quadratic error are provided in this case. The expression of the quadratic error demonstrates that the SPH solution does not converge to the exact one in general under the specified conditions. In this case, the error converges to a residue, which is in general large compared to the norm of the exact solution. The cases with two and three neighbors are also studied. An analytical study in the case of two neighbors is performed, showing how the kernel influences the accuracy of the solution through modifying the condition number of the referred system of linear equations. In addition to that, a numerical investigation is carried out using several Wendland kernel formulas. When two and three neighbors are involved it is found that the error tends in most cases to a small limiting value, different from zero, while divergent solutions are also found in the case of two neighbors with the Wendland Kernel C2. Other cases with more neighbors are also considered. In general, the Wendland Kernel C2 turns out to be the worst choice, as the solution is divergent for certain values of the ratio between the particle spacing and the smoothing length, associated with an ill-conditioned matrix.

Key words: Smoothed particle hydrodynamics (SPH), convergence, water-at-rest problem, hydrostatic tank, free surface, kernel truncation

Introduction

The important topic of convergence in smoothed particle hydrodynamics (SPH) has been scarcely treated in the literature, with few theoretical studies[1-4]. The convergence of the SPH semidiscrete scheme for the Euler equation is treated in detail by Franz and Wendland[5]. In that work, an interesting constructive kernel theory, based on the convergence results, is proposed. Other works[6]have focused on the consistency of the integral formulation in the presence of solid boundaries (being consistency,together with stability, sufficient conditions for convergence).

In order to avoid having to explicitly define the position of free surfaces in SPH, kernel truncation is used instead. The fact that free surfaces do not need to be explicitly tracked in SPH is a major advantage,gaining in simplicity over other numerical methods.The benefits of using the SPH method to simulate free-surface flows have been widely treated in the literature[7]. However, kernel truncation leads to consistency and convergence issues. Some of the problems appearing because of kernel truncation were investigated by Colagrossi et al.[8], where several alternative formulations for both the divergence of the velocity and the gradient of the pressure were proposed.

The integral SPH formulation of the truncated SPH problem was studied by Macià et al.[9]using analytical tools. In that work, the focus is set on the convergence to the linear field of the solution for the pressure. The convergence of the integral SPH solution to the exact pressure field was shown under the assumption that it does not present significant variations at scales smaller than h. The numerical results obtained in that work suggested however that this assumption is not satisfied by the solution of the SPH truncated integral formulation. Furthermore, the numerical solutions that were tested suggested that this spurious influence could be circumvented in the numerics by keeping constant and close to one the ratio between the particle distance and the smoothing length.

Adding to that effort, this work aims to delve into that approach, now focusing in the discrete problem when the ratio between the particle spacing and the smoothing length is kept constant. In Macià et al.[9]the numerics were treated as a complement of the analytical studies. A more detailed study of the numerical solution will be pursued herein. The numerical solution will be studied by considering different ranges of the aforementioned ratio,corresponding to the situations where each particle is influenced by one, two or three neighbor particles in the distance of a kernel radius. Also, other cases with more than three neighbors will be addressed.

In order to pursue these objectives, the paper is organized as follows: First, the truncated discrete SPH hydrostatic problem is presented, then, the numerical solution is studied in the three aforementioned cases(one, two and three neighbors), finally, some ideas regarding the general case of any number of neighbors are introduced.

1. The problem and the setting

The object of interest in this article is the following discrete equation

in which the Dirichlet boundary condition corresponds to the dynamic boundary condition for the pressure at the free surface, set at x=0, and whose exact solution is p(x)=-x. The slope imposed to the pressure field in Eq. (2) is set to -1 without lose of generality. Note that the function p(x) can be seen as dimension-less pressure normalized using any given slope so that leads it to satisfy Eq. (2).

The continuous limit of problem (1) is the integral equation

Note that the Dirichlet boundary condition p(0)=0 is imposed in both the numerical SPH problem and in its integral version by extending the pressure field by zero for x>0. This approach was referred to as U0M in Macià et al.[6], and its consistency studied when used to enforce no-slip boundary condition for the velocity in solid boundaries. The extension of the solution by zero leads to a truncation in the integral at the boundary. How this approach affects the convergence of the solution was discussed by Macià et al.[9]Through all this article, Whstands for the SPH kernel. A detailed description of this function is given below.

In the discrete version, Eq. (1), a staggered grid,G:={xa,a=1,…,N} with fixed particles distributed uniformly starting in -Δx/2 and spaced by a distance Δx=H/N is used

This is equivalent to assigning a volume Δx to each point, thus creating a partition of the domain with particles, as in SPH. Also, a finite domain of depth H=1 (in consistent units) is considered. Note the fact that the domain is infinite in the original problem,Eq. (2), and needs to be truncated in order to numerically solve the discrete problem. The following mirroring technique for thepressure field, presented by Bouscasseet al.[10], isused atthe boundary x=-H:

The use of this mirroring technique assures the consistency in the boundary x=-H, as the exact linear solution will be linearly extended for x<-H.This is the usual way to deal with this kind of boundary truncation problems in SPH when a hydrostatic pressure is involved. Ghost particles are used together with this mirroring. The mirrored neighbor of a given particle, xm, is noted xmG. The mirrored discrete values of the pressure field at these ghost particles, calculated using Eq. (5), are noted pmG.

The SPH kernel in one dimension, is of the form

The function W is assumed to be differentiable,even and of compact support attaining its maximum at x = 0. The parameter h is called the smoothing length, and is the main length scale associated with the integral approximation. The radius R of the support of the kernel is proportional to h . The proportionality constant depends on the choice of the function W . In the present paper we choose h such that R= 2h . The requirement of having a positive definite kernel is important in Lagrangian SPH to guarantee the fulfillment of second principle of thermodynamics in dissipative interactions[11-12].However, this requirement can be relaxed in Eulerian SPH to achieve higher order of convergence[13].

In addition, it is assumed that the Fourier transform of W

Examples of such functions are the exponential[14]or the Wendland kernels[15]. Under the hypotheses above,the kernel function always takes values between 0 and W (0)/h .

The derivative of the kernel function appears frequently through the paper. The function

assumed to be bounded and non-negative, will also be used. Note that, in one dimension

2. Numerical analysis

In this section, the discrete solution of Eq. (1) is presented. An analysis of the discrete scheme is performed for different significant ranges of the Δx/ h ratio. These ranges are considered such that the number of neighbors in which the kernel takes values different from zero is fixed. The Δx/ h ratio will be denoted by η in the following. The ranges of η under study are:

(1) 1.0 ≤η < 2.0, in which one has a single neighbor at each side of the considered particle.

(2) 2/3≤η<1.0 in which one has two neighbors at each side of the considered particle.

(3) 0.5 ≤ η< 2/3, in which one has three neighbors at each side of the considered particle.

Some remarks will also be given for cases with more neighbors, i.e., η< 0.5. Figure 1 shows the different positions of the particles within a kernel radius corresponding to the limiting cases of the three aforementioned ranges.

Fig. 1 (Color online) Sketch showing the positions of the neighbors within a kernel radius for limiting cases of the three ranges considered. In purple, η = 2.0 and η =1.0 , upper and lower limits of the range with one neighbor, in orange, the case η = 2/3 , lower limit of the range with two neighbors, and in green the case η = 0.5 , lower limit of the range with three neighbors

2.1 The case of a single neighbor: 1.0 ≤η< 2.0

The exact solution can be explicitly calculated in this case. Consider the linear system described in Eq.(1). The role of the parameter η is brought out after expressing the problem in terms of F instead of hW′ so, for the sake of clarity, this function will be used in the following.

After applying Eq. (8), the linear system in the case of one neighbor takes the form:

Solving for p2in the first equation of Eq. (9) gives

which is, in general, different from the exact value of the pressure at the second particle

Notice now that the equations corresponding to odd and even particles are decoupled, being the last equation, involving the ghost particle, the only link between the two groups of equations.

Solving for p4in the third equation of Eq. (9) gives:

and after Eq. (10) one gets

Similar considerations lead to p6=3p2and, in general:

The values of every even particle are therefore explicitly determined by Eq. (14). Now, the values of the odd particles need to be calculated after linking both groups of equations through the ghost particle.Note that the mirroring extension plays a major role here. With the considered mirroring, Eq. (5), the extended pressure field in the ghost particle is

and the SPH discrete approximation of the pressure in the last particle is

This expression can be simplified as

The two groups of equations are now linked.Descending through the odd equations, the SPH pressure in the rest of particles follows the rule:

When the number of particles is even, the value of pNcorresponds to an even particle and is thus calculated as pN=p2N/2. The last equation reads

Note that the unknown is now pN-1instead of pN.Following similar considerations to those made above,the SPH pressure at the particle xN-1is obtained,leading to:

From Eqs. (10), (14), (17) and (18), the expression of the general SPH discrete solution for 1.0≤η<2.0 can be written, when the total number of particles is odd, as:

A measure of the numerical solution’s slope can be calculated considering only the even particles. The slope that could be depicted from these particles is given by

Note that the ordering of the particles that are used to compute the slope is p2-p4instead of p4-p2.This is simply a consequence of the choice of the ordering we made for the particles of the grid, which is defined in Eq. (4). This slope is in general different from the exact one, -1. Furthermore, depending on the actual value of η, the evaluation of the kernel takes smaller or bigger values, delivering different numerical slopes. The ratio between the numerical and the exact slopes is

This parameter will be useful in what follows in order to clarify the role of the kernel function in the convergence of the numerical solution.

The expression inside the square root can be developed as

The sum of the second term in the right hand side is

In order to compute the sum of the remaining two terms we need to use the numerical solution. For simplicity, only the case with even number of particles will be considered. The case with odd number of particles can be dealt analogously leading to the same conclusions.

In order to sum the first term in the right hand side of Eq. (23), the odd particles will be considered together with their following even particle. According to the definition given in Eq. (20), the sum of the values of the numerical solution up to theth2n (even)and 2n-1th(odd) particles is

Introducing this in Eq. (23), and reordering the indexes one gets

the first term corresponds to the even particles and the second to the odd ones. After a straightforward manipulation, this can be simplified as

Following a similar strategy, the last term in Eq. (23)reads

Summing up the three terms (24)-(26), and taking into account that N=H/Δx, the quadratic error is found out to be

Equation (27) expresses the behavior of the error depending on the approximation parameters, Δx and h and their ratio η.

Using Eq. (27) we can compute the quadratic error for a fixed η when both Δx and h tend to zero. Note that this limit is not the lowest error we could expect for a given η, as some of the terms in Eq. (27) depending on h and Δx are negative, so non-zero values of the approximation parameters might lead to lower values of E. What this limit tells us is how big the error is when the approximation parameters tend to zero while keeping a constant η ratio. The best approximation might not be achieved in general when Δx and h are closer to 0, but for those values minimizing Eq. (27). This shows that the numerical scheme although non-convergent, might be useful for practical purposes in some cases.

Fig. 2 (Color online) Relative quadratic error in terms of η,given by Eq. (27), with Δx=0 and h=0 Wendland Kernel C2is used

Let us now focus on the Wendland Kernel C2.Consider the relative quadratic error, computed as the quotient between the quadratic error, E, given by Eq.(22) and the quadratic norm of the exact solution

Figure 2 depicts the limit of the relative quadratic error, Eq. (28), depending on η, when both Δx and h tend to 0 using this particular kernel. It can be seen that the curve presents two minima where the error approaches zero. An explanation to this can be found in Eq. (21). This equation gives the ratio between the slope of the numerical and the exact solutions. From expression (21) it is clear that this ratio only depends on η and on the shape of the function F. This factor can be seen as the ratio between the functions 1/(2s3)and F(s). In Fig. 3 these two functions and their quotient are plotted together. It is clear from the figure that, for two particular values of η, the factor Q (defined in Eq. (21)) is exactly one, and therefore the numerical slope coincides with the derivative of the exact solution. These values of η are precisely those at which the minima in the quadratic error expression are attained, (see Fig. 2).

Fig. 3 (Color online) Derivative of the kernel, F(s), and function 1/(2s3), appearing in the quotient of the numerical and the exact slopes. The ratio between them,Q, is also plotted. The smoothing length is set to h=1.0. Wendland Kernel C2 is used

The other two kernels considered in this work do not present the same behavior. In that cases, the ratio between the numerical and exact slopes stays always greater than one, so convergence to the exact solution is never achieved.With both kernels, the smallest errors are achieved when η=1.0, but they are different from zero anyway. A residual error different from zero is obtained with those kernels in the whole range of η.

2.2 The case of two neighbors

In this range of η, each particle interacts with two neighbors in a kernel radius. The limit cases η=1.0 and η=2/3 correspond to the situations in which the second and the third neighbor, respectively,are placed in the limit of the kernel’s compact support.

Let us first focus on the Wendland Kernel C2,the most used of the three kernels considered. The numerical solution presents two different behaviors:

(1) In the range 2/3<η<0.175 the solution diverges with h/H.

(2) In the range 0.175<η<1.0 the quadratic error converges to a residual error that depends on η.A more detailed classification can be made,establishing subranges of η where different behaviors can be described with precision:

(a) For 0.8<η<1.0 the behavior is similar to that of the case with one neighbor, i.e., a wrong slope is present in the numerical solution, even with h/H small. This prevents the numerical solution from converging to the exact one. Instead of that, the quadratic error converges to a finite quantity. This residue, and the error in the numerical slope, are bigger for values of η close to 1.0. Several examples in this range are shown in Fig. 4, where the dependence of the numerical slope on η is plotted.

(b) In the range 0.78<η<0.793 the numerical solution also presents an incorrect slope. In contrast with the first case addressed, the error in the numerical solution, pk-p(xk), exhibits a change of sign at a certain value x/H that depends on η.Besides, the behavior of the quadratic error is not uniform with h/H: it decreases with h/H above a certain limit, which depends on η, and grows below it. Figure 5 shows the quadratic errors for different values in this range.

(3) For 0.715<η<0.77 the error grows with h/H and converges to a certain residue, which value depends on η. A change in the sign of the error is also found in this case. Figure 6 shows the quadratic errors for several values in this range.

(4) The last range is 2/3<η<0.715. The solution in the limit case η=2/3 is shown in Fig. 6.Note the strong spurious behavior appearing in the numerical solution, with oscillations of the same order of magnitude of the function values. The numerical solutions within the whole range present similar oscillations. The amplitude of this oscillations decays as η grows.

Fig. 4 (Color online) SPH numerical solutions in the range 0.8<η<1.0. 100 particles are used. Wendland Kernel C2 is used

Fig. 5 (Color online) Quadratic error depending on h/H for η constant in the ranges 0.715<η<0.718 and 0.78<η<0.793. Wendland Kernel C2 is used

Fig. 6 Numerical solution for η=2/3. 200 particles are used

Summing up, the numerical solution does not converge to the exact one in the whole range considered when the Wendland Kernel2C is used.Despite of that, the error converges in general to a finite residue except for the lowest values of η. The solution in the cases with η close to 2/3 is found to be particularly spurious, exhibiting divergent quadratic errors. The exact mathematical mechanisms underlying this behavior are not fully understood yet.Nevertheless, some insight into this matter is given next.

Let us examine the behavior of the solution when the other kernels are considered: the Wendland Kernels C4and C6. In contrast with the previous kernel analyzed, the solution is not divergent in these cases for η close to 2/3. The quadratic error converges instead to a finite residue in the whole range of η. In Fig. 7 the quadratic errors for different values of η when the Wendland Kernel4C is used are shown. Note that the error always converges. Also,it is a surprising fact that the case η=2/3 exhibits the smaller errors, in contrast with the divergence found with the Wendland Kernel2C. It is also noticeable that the values of η delivering the best approximations are not the same observed in the case of the Wendland Kernel2C. The behavior when the Wendland Kernel6C is used is found to be similar.

The kernel has therefore a great influence on the numerical solution. The mechanism under this influence can be explored by analyzing the matrix of problem (1), that will be note by A. The general form of this matrix when two neighbors are considered is the following The factor C depends on the kernel, but it only acts as a scaling parameter, without modifying the structure of the matrix. In view of that, it can be concluded that this structure depends on the kernel only through the parameter λ. The influence of the kernel on the numerical solution can therefore be studied by analyzing how the parameter λ influences the properties of the matrix. We refer mainly to the condition number of the matrix, κ(A),a parameter expressing the numerical solvability of the linear system.

Fig. 7 (Color online) Norm of the error for η in the range 2/3<η<1.0. Wendland Kernel 2C is used

The condition number depends on the choice of the matrix norm. Through all this work, the matrix’s largest singular value (2-norm) is used as matrix norm.Consider now the matrix A, as described by Eq. (29).This matrix is parametrized by two numbers, C and λ. The condition number of the matrix does not depend on C, as this is simply a scaling factor. The dependence of the condition number on the parameter λ is shown in Fig. 8. This picture shows that the matrix is ill-conditioned for values of λ close to two.In addition to that, Fig. 9 shows the dependence of λ,given by Eq. (31), on the parameter η for the three kernels considered. From this picture is clear that,when the Wendland kernel2C is used, lower values of η in the considered range will lead to an illconditioned system, as the value of λ approaches λ=2.0 when the ratio η is made smaller. This is the reason of the spurious behavior reported above in the case of this kernel for lower values of η in the range. The solution shown in Fig. 6 is an example of this phenomenon. Note that the horizontal axis in Fig. 9 is limited at 0.9 instead of 1.0, which is the limit of the considered range. It has been numerically checked that λ grows very quickly in that range of η, and that the variations of the condition number within that range of λ are lower than 10-6.

Fig. 8 (Color online) Condition number of the matrix A depending on λ. Matrices of size 400×400 are considered

Fig. 9 (Color online) Dependence of λ on η for the three kernels tested

2.3 The case of three neighbors: 0.5<η<2/3

In the range 0.5<η<2/3 each particle interacts with three neighbors within a distance of a kernel radius. The norm of the error has been computed numerically for different values of η and using the three kernels under study. When the Wendland kernels C4and C6are used, the solutions are always convergent, while the error tends to a limiting value different from zero. When the Wendland Kernel C2is used the norms of the error also converge with h/H. Despite of that, it is found that, within certain ranges of h/H, the norm of the error grows unexpectedly and decays again when this kernel is used.

When the Wendland Kernel4C is used, the quadratic error converges in general to a small residue when h/H tends to zero, as it is shown in Fig. 10.When the Wendland Kernel6C is used, similar results are obtained. The norm of the errors tends to a finite residue in the whole range of η as it is shown in Fig. 11.

The case of the Wendland Kernel2C is slightly more complicated to analyze, as has been already introduced. The general trend is that the norm of the error converges to some small residue. However,some spurious cases are found, for which the norm of the error explodes, in specific ranges of h/H,decaying again when h/H is made smaller. This behavior seems to be also linked to ill-conditioned matrices, however further inspection is required to understand exactly the mechanisms under this phenomenon. Similarly to the other ranges of η considered, the Wendland Kernel2C exhibits the worst results.

Fig. 10 (Color online) Norm of the error depending on h/H for η in the range 0.5<η<2/3. Wendland Kernel C4 is used

Fig. 11 (Color online) Norm of the error depending on h/H for η in the range 0.5<η<2/3. Wendland Kernel C6 is used

2.4 Other cases with more neighbors: η<0.5

Following the ideas presented above, other cases with more neighbors have been explored. The most remarkable finding is that when the Wendland Kernel C2is used, the pattern adopted by the numerical solution for η close to 2/3, shown in Fig. 6, appears regularly for different values of η. This leads to divergent solutions, as it is shown in Fig. 12. The values of η presenting this behavior, within the considered cases up to 12 neighbors, correspond to the lower limits of the ranges of η in which and even number of neighbors is involved. These limits are of the form η=2/(2n+1), n=1,2,…. The reason why these particular values of η produce divergent solutions is still to be understood. However, this phenomenon is linked to matrices with high condition number, as it is shown in Fig. 13. In addition to that,Fig. 13 also shows how the condition number grows in other ranges of η, different to those previously addressed. This fact suggests that further inspection into the structure of the matrix could enlighten the reasons of this spurious behavior.

Fig. 12 (Color online) Norm of the error depending on h/H for fixed η. Wendland Kernel 2C is used

Fig. 13 (Color online) Condition number of the matrix depending on η. Wendland Kernel 2C is used and 250 particles are considered in the calculation

3. Conclusions

The convergence of the solution to the discrete SPH approximation of the hydrostatic equation for the pressure has been investigated in this paper. The canonical first order differential operator has been used and no kernel correction has been considered at the free surface, where the Dirichlet boundary condition is imposed. At the discrete level, the problem takes the form of a system of linear equations.The analysis has been performed by making the particle spacing, Δx and the smoothing length, h,tend to zero, keeping constant their ratio η=Δx/h .The cases with one, two and three neighbors have been studied in detail.

In the case of one neighbor, an explicit expression for the numerical solution has been given.Furthermore, an explicit expression for the quadratic error has also been provided, demonstrating the relation between the quadratic error and the approximation parameters and their ratio. It has been shown that the numerical solution does not converge in general to the exact pressure field in this case.

In the case of two neighbors it has been numerically shown that the Wendland Kernels C4and C6behave better than the Wendland Kernel C2. For the three kernels tested, it has been found that the error tends in general to a small residue, different from zero, and that divergent solutions appear only when the Wendland Kernel C2is used. In that case,the numerical solution presents a particularly spurious behavior when ηapproaches its lower limit, η=2/3.In addition to that, a link between this divergence in the solution and an ill-conditioned matrix has been established. A parametric study of the matrix involved in the linear system has been presented, showing how the choice of the kernel affects the structure of the matrix.

When a three-neighbor interaction is considered,it has been found that the error also tends to a limiting residue, different from zero, and smaller in general than the residue when two neighbors are used. In this case, numerical approximations are convergent in general with the three kernels. In the case of the Wendland Kernel2C the convergence is non monotonic.

A particularly relevant consequence of the present analysis is the fact that the kernel that is used has a great influence in the solution and its accuracy.Furthermore, it has been shown that non convergent cases are linked to ill-conditioned matrices. The exact mathematical link between these high condition numbers and the choice of the kernel is to be established.

These results show that the belief that keeping a constant and close to one η ratio leads to convergent solutions is not true in general. On the contrary, it has been shown that for wide ranges of η, the numerical solution do not converge to the exact one.

It is interesting the fact that more regular kernels,as it is the case of Wendland Kernels C4and C6,lead to better solutions when the ratio η is kept constant. However, when the continuous integral form of the problem is considered, less regular kernels, i.e.,the Wendland Kernel2C, provide better results, as is addressed in a previous work by the authors.

Many interesting questions arise from these results regarding the theoretical analysis of the numerical SPH solution and are left for future work.The relationship between the condition number and the convergence of the solution, when the ratio η is kept constant, will be investigated. It will also be assessed how the choice of the kernel influences the structure of the matrix, and thus, the convergence. In addition to that, the behavior of the scheme in the context of an evolution problem, for which the solution presented here is the steady state, is to be studied. This analysis will be also applied to the symmetrized and antisymmetized operators, more popular among SPH practitioners than the one presented here. This study will provide guidance on the choice of operator depending on the particular setting of the problem.

Acknowledgements

This work was supported by the Spanish ministry of Innovation and Universities (MCIU) (Grants Nos.MTM2017-85934-C3-3-P, RTI2018-096791-B-C21,“Hidrodinámica de elementos de amortiguamiento del movimiento de aerogeneradores flotantes”), P.E.Merino-Alonso is supported during the completion of his Ph. D. Thesis by MeyFP (Grant No.FPU17/05433).