A Stable SPH Discretization of the Elliptic Operator with Heterogeneous Coefficients

Smoothed particle hydrodynamics (SPH) has been extensively used to model high and low Reynolds number flows, free surface flows and collapse of dams, study pore-scale flow and dispersion, elasticity, and thermal problems. In different applications, it is required to have a stable and accurate discretization of the elliptic operator with homogeneous and heterogeneous coefficients. In this paper, the stability and approximation analysis of different SPH discretization schemes (traditional and new) of the diagonal elliptic operator for homogeneous and heterogeneous media are presented. The optimum and new discretization scheme is also proposed. This scheme enhances the Laplace approximation (Brookshaw's scheme (1985) and Schwaiger's scheme (2008)) used in the SPH community for thermal, viscous, and pressure projection problems with an isotropic elliptic operator. The numerical results are illustrated by numerical examples, where the comparison between different versions of the meshless discretization methods are presented.


Introduction
Smoothed particle hydrodynamics (SPH) was developed a few decades ago to model inviscid fluid and gas flow dynamics in astrophysical problems [1][2][3][4]. The SPH is an interpolation-based numerical technique that can be used to solve systems of partial differential equations (PDEs) using either Lagrangian or Eulerian descriptions. The nature of SPH method allows to incorporate different physical and chemical effects into the discretized governing equations with relatively small code-development effort. In addition, geometrically complex and/or dynamic boundaries, and interfaces can be handled without undue difficulty. The SPH numerical procedure of calculating state variables (i.e., density, velocity, and gradient of deformation) is computed as a weighted average of values in a local region.
Despite a few advantages of SPH method, this method is not free from disadvantages. For example, for fluids, gases, or solids with non-trivial boundaries there is incompleteness of the kernel support combined with the lack of consistency of the kernel interpolation in conventional meshless methods which results in fuzzy boundaries. In some cases, this can be fixed by an automatic incorporation of the boundary condition [1]. The completeness of mesh-free particle methods was discussed in [33]. However, care must be taken to ensure that variables whose values do not approach zero at boundaries are accurately represented.
It has been observed in the literature that meshless methods (e.g., SPH) are not free from instabilities, especially in the modeling solid mechanics problems. For example, the tensile instabilities were identified in [34] by performing the Neumann analysis of the one-dimensional governing equations (conservation laws). Therefore, different stabilization techniques have been developed. It is important to note at this point that it is very difficult to perform the general Neumann analysis in two-and three-dimensions. Furthermore, the high-frequency instability results from the low order discretization (rank deficient) of the divergence operator. The tensile instability results from the interaction between the second derivative of the Eulerian kernel (i.e., computed kernel in the Eulerian coordinates) and the tensile stress. The tensile instability only occurs for the Eulerian kernels because the Eulerian kernels depend on both the stress and the second derivative of the kernel. It has been shown [35] that in the case where the kernel is a function of the Lagrangian coordinates (the Lagrangian kernel), tensile instability does not occur. A comprehensive analysis on this subject can be found in [35]. Recently, some regularization and stabilization of SPH schemes were proposed in [36,37]. However, there is a lack of the mathematical analysis of these numerical techniques.
Since its introduction, SPH has been successfully used to model a wide range of fluid flows and the behavior of solids subjected to large deformations. For example, the SPH method was applied to simulate high energy explosions [38] and impact [39,40], most notably free surface flows and collapse of dams [41], elastoplasticity [42][43][44]40], to model low Reynolds number flows [45,46], to study pore-scale flow and dispersion [47,48], and for thermal problems [43,49,50]. A comprehensive overview of SPH applications to incompressible or nearly incompressible flow can be found in [51].
In different applications including, but not limited to, fluid flow related problems, it is required to have the stable and accurate discretization of the following operator (diagonal elliptic operator), which is the research subject of this paper: L (u) = −∇ · (M (r) ∇u (r)) − g (r) , ∀r ∈ Ω ⊂ R n , M (r) ∈ {diag[m(r)] : m(r) ∈ R n + , m (r) ∈ L 2 (Ω)}, diag : R n → R n×n , where u (r) is the unknown scalar or vector variable field, M αβ (r) is the diagonal matrix of the mobility field. For the reminder of the paper, lower-case Greek indices are used to describe the components of any tensor, the upper-case Latin indices are used for particle indices and lower-case Latin indices are used to describe the index of the Cartesian coordinates. One example of the mobility field includes M αβ (r) = m (r) δ αβ , α, β = 1, . . . , n; where m (r) is the mobility (or diffusivity, viscosity) scalar field, δ αβ is the Kronecker symbol, n = 1, 2, 3 is the spatial dimension. It is important to note that we do not make any distinction between contravariant and covariant components of the any tensor in this paper by assuming that these tensors are defined in a Cartesian coordinate system. The sink/source term g (r) is assumed to be zero in this paper. Consider the operator in the expression (1) with piecewise continuous coefficients M (r) in Ω.
It has been noticed [52] that some of the numerical methods for elliptic equations may violate the maximum principle (i.e. lead to spurious oscillations). Therefore, proposed methods must satisfy a discrete maximum principle to avoid any spurious oscillations. This is also applicable to meshless discretizations. Usually, the oscillations are closely related to the poor approximation of the variable gradient ∇u in the flux computation. In this paper, different numerical discretizations of the elliptic isotropic operator are analyzed. The objective of this paper is to develop numerical scheme satisfying the two-point flux approximation nature in the form where J is the neighboring particle of the particle I, V r J is the J-particle volume, Ψ is the special kernel. This structure (2) allows to have a better analysis of the stability and monotonicity. Furthermore, this will allow to apply upwinding strategy during the solution of nonlinear PDEs. The new optimum (i.e., high order scheme with monotonicity property) discretization scheme of the shape (2) is proposed in this paper. This scheme is based both on the Laplace approximation (Brookshaw's scheme [53]) and on a gradient approximation commonly used in the SPH community for thermal, viscous, and pressure projection problems. The proposed discretization scheme is combined with mixed corrections, which ensure linear completeness. The mixed correction utilizes Shepard Functions in combination with a correction to derivative approximations. In corrected meshless methods, the domain boundaries and field variables at the boundaries are approximated with the improved accuracy comparing to the conventional SPH method. The resulting new scheme improves the particle deficiency (kernel support incompleteness) problem. The outline of the paper is as follows. In Section 2, the existing discretizations of the Laplace operator with the building blocks necessary for these methods are discussed. In this section, the SPH kernel and its gradient properties are also discussed. The description of meshless transmissibilities and their connections to the existing mesh-dependent discretization schemes is given in Section 3 including the construction of a new meshless discretization scheme. The approximation, stability, and monotonicity analysis are performed in Section 4. The numerical analysis of different boundary value problems is presented in Section 5. The paper is concluded by Section 6.

SPH discretization of the Laplace operator
The meshless approximations (SPH approximations) to the operator (1) with heterogeneous and homogeneous coefficients are presented in this section. Let us consider a rectangle in R n , n = 1, 2, 3: as the numerical domain. Here, a α are the center coordinates of the rectangular and l α are the side lengths. The following L p (Ω) , L p (Ω) norms are used to quantify the accuracy of different approximations for the entire numerical domain Ω: where f is the approximated physical quantity, V r K is the volume of the particle located at r K ,Ω denotes entire domain including boundary particles,Ω denotes only internal part of the domain Ω. The particle volumes V r K are selected to satisfy during numerical simulations only The proposed discretization schemes should also provide a good accuracy of the operator (1) with a discontinuous m (r) (or piecewise continuous function) coefficients since these coefficient cannot be differentiable in the classical sense. The standard SPH spatial discretization of the Laplace operator (1) arises from the following relations with g (r) = 0, ∀r ∈ Ω r I ,h I and W ∈ H 1 ( Ω r I ,h I ) (i.e., we consider kernel functions from Sobolev space with L 2 -norms of the function itself and its gradient): where W ( r − r I ,h I ) is the kernel that weakly approximates the Dirac delta function δ (r − r I ) but with finite characteristic widthh I =h(r I ) (or physical smoothing length) around the particle I (subscript is omitted ifh refers to any particles). The effective characteristic widthh I is defined by (6) using smoothing particle length h I around the particle I. The control volumes in the meshless discretization are the patches, which are interior to the support of the kernels and r, r I are points in Euclidean space R n , ''supp'' indicates the support of the kernels. Additionally, it can be required that the kernels are radially symmetric and compactly supported ) leading to following definitions: where ∆x IJ is the distance between particles I and J (we assume that particles do not overlap), ∆x I,m is the radius of space around particle I which is not filled by particles (∆x m refers to ∆x I,m at any particle I), ∆x I is the maximum distance between particles I and J inside control volume Ω r I ,h I (∆x refers to ∆x I at any particle I),h I,m is the inner smoothing length (h I,m ≤h I ),h IJ is the effective smoothing length between particles I and J (which will be defined in Section 2.2). However, it is important to note that, in the case of the homogeneous particle distribution, it is common to useh I =h I,m = h = f · ∆x = const, where ∆x I,m = ∆x is the inter-particles distance, f is the scaling factor. Hence, it is important that pairs of length scales (h I,m , ∆x I,m ) for heterogeneous particle distribution and (h, ∆x) for homogeneous particle distribution are distinguished during the analysis.
In relation (5), if ∇u (r) is smooth and kernels are radially symmetric and compactly supported then ∇u (r) can be expanded in a Taylor series about r I and it can be shown using integration by parts that the smoothing error due to (similar to [54]) in the largest disk in R n where the Taylor series converges: ) .
From the definition of the Ω r I ,h I , it follows that there is an infinite cover of the numerical domain Ω: According to the Heine-Borel theorem, there is a finite subcover (since we consider only compact numerical domain Ω), that is where N is the number of particles (i.e., patches) in the numerical discretization.
The SPH spatial discretization of the integral (7) is defined over the control volumes Ω r I ,h I to obtain the final discretization of the Laplace operator. The final step in the particle method is to approximate the integral relation on the right-hand side of the (7) using Monte-Carlo expressions or any cubature rules [55,43,33,15] and which is known as the particle approximation step and resulting in the approximate form of the operator ⟨L (u (r I ))⟩, where integration is substituted by summation. Hence, using SPH discretization technique and assuming M αβ (r) = δ αβ , u (r) = r, the expression (5) can be written for ∀r I ∈ Ω r I ,h I and ∀{r K , where N is the SPH correction vector, which shows the deviation from the zero vector. The strategy of selecting V r K plays a crucial role depending on the particle position r K , kernel W ( r K − r I ,h I ) and the smoothing lengthh I . In the following section, we formulate the limiting condition for the particle volumes V r K ∀K which we call the SPH continuity condition.

Remark 1.
A sufficient condition for (10) to hold (i.e., double limit and the iterated limits exist and are equal) is Moore- converges uniformly forh I ̸ = 0 then the double limit and the iterated limits exist and are equal.
The assumption that the boundary term from the integration by parts is zero in (5) is valid only in regions where the kernel has a full support, or the function, or the gradient of the function itself is zero. For particles near free surfaces or boundary, the neglect of these terms leads to significant errors for boundary value problems. Several techniques have been developed to address these errors through various correction methods, e.g., by calculating the boundary integrals [15]. In addition, this can be corrected as it will be shown below by applying normalized corrected meshless methods in the derivative approximations. In the following section, the commonly used SPH kernels and its gradients are considered.

SPH kernel and its gradient
A central point of the SPH formalism is the concept of the interpolating function (or kernel) through which the continuum properties of the medium are recovered from a discrete sample of N points (8) with prescribed mass (for conventional Lagrangian methods) or volume V I (for fully Eulerian methods). In the Lagrangian formulation, these points move according to the specified governing laws, whereas these points are fixed in space for the Eulerian formulation. A good interpolating kernel must satisfy a few basic requirements: it must weakly tend to the delta function in the continuum limit and has to be a continuous function with piecewise first derivatives at least.
From a more practical point of view it is also advisable to deal with symmetric finite range kernels, the latter to avoid O calculations. In addition, instead of using theoretical finite characteristic smoothing lengthh I in the definition of the kernel, the effective smoothing lengthh IJ between particles I and J is used directly in the definition of the kernel expression in practical calculations. Cubic and quintic splines are the commonly used kernels in SPH formulations [56,15]. Since the quintic spline does not provide the numerical advantages, the cubic spline is used in this paper: Definition 1. The kernel function centered at particle r I is defined as where z J =   r J − r I   2 /h IJ is the dimensionless variable,Ŵ is the kernel in a dimensionless form, and Ξ is the normalization factor equal to 3/2, 10/(7π ), and 1/π in 1D, 2D, and 3D, respectively. In continuous case, the particle index is omitted as z for dimensionless variable.
The different choices of computing effective smoothing length between particles used in this paper will be discussed in the upcoming subsection. Although, the kernel is normalized in a continuous sense, it is important to note that does not satisfy the normalization condition in the discrete space due to the particle distribution and incomplete kernel support near the boundary and, hence, the discretized-normalized kernel function can be considered: where ν (r I ) is the specific concentration of particle r I (i.e., it is approximately the inverse of the compact volume support multiplied by the particle I volume and ν (r I ) ∼ O (1)) which has a larger value in a dense particle region than in a dilute particle region. In regions of the high particle density, the denominator in (13) is high resulting in lower values of the kernel W . Thus, the denominator normalizes the kernel function to ensure that the kernel W forms a local partitioning of unity regardless of the particle distribution within the Ω r I ,h I = supp W ( r J − r I ,h I ) . The discretized-normalized kernel function will also be used in the discretization schemes below. At this point, all possible options of computing ∇W (r J − r I ,h IJ ) are listed: where ν(r I ) is the specific volume of the particle located at the point r I .
Additionally, two options can be written as where ν(r I ) is assumed to be a constant during the differentiation with respect to r J .
The alternative case is when we assume that ν(r I ) can be considered as a function of the location of other particles leading to the following relations Where ∇W (r J − r I ,h IJ ) and ∇ r I in all definitions above represent the nabla operator with respect to r I . Assuming this for the rest of the paper, the r I index is omitted for the remainder of this manuscript.

Remark 2.
In general, we assume in this paper some variations of the smoothing lengthh(r) and, hence, definitions of the kernel gradient should include ∇h. However, we consider the numerical examples where ∥∇h∥ ≪ 1 and terms with ∇h are omitted.
Similar to the Moving Least Squares (MLS) method [57], Eqs. (16) and (18) are two different forms of ''full derivatives''. At the same time, Eqs. (15) and (17) are two different forms of ''diffuse derivatives''. From these options, it follows that ''full derivatives'' are connected with the differentiation of a discrete function and ''diffuse derivatives'' are connected with the differentiation of an exact function. Two types of derivatives are different in treating the normalization factor ν(r I ) where ''full derivatives'' assume ν(r I ) as function and ''diffuse derivatives'' assume ν(r I ) as constant. Both types of derivatives have some advantages and disadvantages and the choice depends on the application. The impact of different options on numerical results will be shown below.
Finally, it is important to note the following relation for the kernel gradients (15), (17), and (18): where function G(z J , r I ) ≤ 0 has a different shape depending on the definitions (15), (17) and (18). Using SPH discretization technique and assuming M αβ (r) = δ αβ , u (r) = ∥r − r I ∥ 2 , r ∈ R n , the expression (5) can be written as The sum (20) contains only positive numbers, hence, the necessary condition to fulfill this relation (20) is to satisfy the SPH continuity condition for selecting particle volumes given their positions and the smoothing length: The condition (21) provides an upper bound for the particle volumes V r J given the dimension of the space n, particle and smoothing lengthh I . It is important at this point to introduce the following quantity which will be used in the error analysis of the proposed SPH schemes: The condition (22) may not converge uniformly forh I ̸ = 0 for some set of particle volumes {V r J }. Hence, the relation between ∆x I andh I (orh I,m ), kernel W and set of particle volumes {V r J } defines the convergence conditions of (10) and (20). The SPH method shows good approximation properties in regions where the kernel has full support. For particles near stress-free surfaces or boundary, the SPH method shows a poor approximation. Several techniques have been developed to address these errors through various correction methods, e.g., by applying normalized-corrected meshless methods in the derivative approximations [33,15,39,40], which requires normalized-corrected definitions of the kernel gradient as follows: where ∇ * α W is also defined with respect to r I and the elements C αβ are defined by the matrix C (r I ): where the summation by repeating Greek indices is only assumed throughout this paper, ∇ * α W is the normalized-corrected gradient of the kernel, and C αβ is the correction symmetric tensor [59]. It was shown [60,61] that the value of the minimum eigenvalue λ C (r I ) of the matrix C −1 based on the discretized-normalized kernel function depends on the particle distribution within the domain Ω r I ,h I = suppW ( r I − r,h I ) . When going away from the Ω r I ,h I domain this eigenvalue tends theoretically to zero, while inside this domain the eigenvalue tends theoretically to one. This important information allows determining regions of the continuum media where free-surfaces are located [60,61]. In this paper, the discretizations of the Laplace operator are based on the gradient of the kernel. Several methods of Laplace discretizations were proposed [62][63][64][65] using second derivatives of the variable u. However, second-order derivatives can often be avoided entirely if the PDE is written in a weak form (5). It is important to note that approximations using second-order derivatives of the kernel are often noisy and sensitive to the particle distribution, particularly for spline kernels of lower orders.

SPH symmetrization of smoothing length
The true particle smoothing length h I may vary both in space and time in general. Therefore, in general case, each particle has its own smoothing length h I . Considering the case where h I ̸ = h J for two different interacting particles I and J and the kernel support based onh I and located in I which covers the particle J but the kernel support located in J does not cover the particle I. In this case, the particle I acting on particle J (produces, e.g., a flux or a force) without particle J acting on the particle I, which leads to a violation of fundamental laws (e.g., mass conservation or Newton's third law) for a closed system of particles. This problem has been resolved by introducing the symmetrization of the smoothing length. In this study, the following symmetrization option is used: It is clear that h =h I =h I,m , h =h IJ ∀I, J and ∆x I,m = ∆x ∀I for homogeneous particle distribution.
This completes the description of basic properties of SPH method allowing to construct all necessary building elements of SPH discretization such as list of neighbors, kernel values, and kernel gradients. The following sections describe the traditional and newly proposed discretization schemes for the Laplace operator. Without loss of generality, the discretization schemes below are build for −⟨L (u (r I ))⟩. Hence, there are two sources of the discretization error: (a) the error (smoothing error) associated with the weak approximation of the Dirac delta function δ (r − r I ), which is related to the smoothing particle lengthh I and (b) the error (discretization error) associated with Monte-Carlo expressions or any cubature rules controlled by ∆x I , which is related to the number of particles and particle distribution within the control volume Ω r I ,h I .

Brookshaw's scheme (1985)
Brookshaw proposed [53] an approximation of the Laplacian for an inhomogeneous scalar field m (r), i.e., M αβ (r) = m (r) δ αβ , α, β = 1, . . . , n; that only includes first order derivatives: where V r J is the volume of the particle J, ∥•∥ is the Euclidean norm throughout this paper, u (r) is the unknown scalar or vector field (e.g., pressure p or velocity the field coefficients. This scheme can be derived by applying the particle approximation step to the right-hand side of (5) with the following assumptions ∇u ( Some special words need to be said about the mobility approximation, which comes in the form The factor 2 is introduced to compensate the factor of 1/2 in the second leading term of the Taylor expansion of the relation (28). Furthermore, the relation (29) allows to capture a heterogeneous mobility field distribution. The theoretical error analysis presented below in this paper requires smoothness properties of the considered functions, e.g., u (r), W and m (r). The cubic kernel W (12) and m (r) have limited smoothness properties and impact of this will be covered below.
Theorem 1. The discretization scheme (27) at a given point r I ∈ Ω ⊂ R n sufficiently far away from the boundary ∂Ω (i.e., where the kernel W has a full support and particle volumes satisfy the SPH continuity condition (21)) has an error (discretization error) bounded by C Proof. Note that the SPH continuity condition (21) is the necessary condition to fulfill the relation (20) and limits the values of V r J ∀J. Using Taylor series expansions about a point r I , the following relations can be written: Since scalar mobility field m (r) ∈ C 1 (Ω) ∩ H 2 (Ω), the following relation can be written: Substituting relations (30)-(31) into the scheme (27), it leads to the following relations: where N α (r I ) = ∑ (this is a generalization of the definition (11)) and tensor Γ αγ is defined by Here, symbols α, γ , ω represent indices such that α, β, γ = 1, . . . , n. While obtaining (32), the following relation has been used (which is not valid for kernel gradient (16) subject to the random particle distribution): The leading term of the error function E in (32) is defined by the following expression where U αβγ (r I ) is defined by Typically, 1D (one dimension) error analysis allows to gain a deep insight into the nature of discretization error due to the smoothing integral in (5) and a discrete version of a given scheme (e.g., widely used scheme (27)) and kernel gradient ∇W . To compute exact error in 1D SPH approximation of the following coefficients N α , Γ αγ , and U αβγ (indices will be omitted) at a particle I in terms of the kernel function and its derivatives at r I , the smoothing length, the particle distributions, the below expression derived in [54] is used for numerical domain of length 2h.
where ∆x, h are the particles volume and smoothing length ∀J, F is the approximated data function with its first F ′ and third F ′′′ derivatives in 1D, W ′ is the kernel derivative in 1D, B p+2 are the Bernoulli numbers, p is the highest integer such that the pth derivative and all lower derivatives are zero at the edges of the compact support of the kernel W (cubic spline (12) has p = 2),Ŵ (m) z=a is the mth derivatives of dimensionless kernelŴ at the point z = a. The conditions for expression (37) to be valid [54]: (a) uniformly spaced particle distribution in 1D, (b) smooth data functions F with its derivatives, and (c) kernel functions which are smooth (except at the support boundary), even and normalized. In addition, it is required that the particle volumes span the compact support without gaps or overlaps.
The following functions are substituted into (37): This leads to the following expressions: From expressions (38), (39), (40), and (7), it follows that the error of discretization only in 1D is the sum of O (smoothing error) and (error due to discretization scheme). Hence, as goes to zero, the accuracy is controlled by smoothing error O . At the same time, as smoothing error h goes to zero, the accuracy and parameter p defines the behavior of the error in 1D subject to the limitations of expression (37).
To examine the effects of non-uniform particle distribution, it is assumed that the numerical domain of length 2h is partitioned into a series of subintervals, where each of which represents the volume of a particle as was done in [54]. Each of these volumes has a value ∆x I and is centered at a locationx I which does not necessarily coincide with the particle location, x I . These assumptions were also made in [54] to derive the following relation in 1D: Dimensionless variables δ J , ∆z J can be replaced by average valuesδ and ∆z and moved outside the summations (introducing a higher order error into the analysis itself) [54]. As a result, we can obtain the following asymptotic for N (x I ), Γ (x I ), and U (x I ) similar to the uniform particle distribution case.
Asymptotic relations (42), (43), (44) show the interaction between parameters h,z andδ. The leading terms are secondorder in , first-order in the nonuniformity parameterδ, and order 1/h I,m in smoothing length. These properties are independent of the smoothness of the kernel. The relation (42) shows that discretization error of the scheme (27)  In general (in R n ) and taking into account limited smoothness properties of W (12), the first term in the expression (32) can be bounded in the following way (by applying triangle inequality): Combining the results of conservatives estimates (45), (46) and 1D analysis, we can conclude that the total error (discretization and smoothing errors) can be bounded by C of the error is consistent with results presented in [58].
As noted in [50], since the Hessian matrix H = ( u ,αγ (r I ) ) is coupled with the tensor Γ αγ via inner product, it is difficult to separate them. However, using the identity [H : Γ ] = tr [HΓ T ], along with the fact that both H and Γ are symmetric, one can write the relation as in [50] The Q αβ is a tensor with zero trace: Note that tensor Γ defined by (33) is identical to C −1 , where C is defined by (24). However, a more general definition of Γ (which reduces in some cases to (24)) will be proposed in the upcoming sections. Hence, the maximum accuracy is achieved when which is difficult to fulfill simultaneously for different kernel gradients. Therefore, the overall discretization error is covering the worse and best cases. Note that in order to maintain the results of the above theorem close to the boundary points, additional ghost particles must be added outside the numerical domain. □
for any particles distribution but not the conditions (34) and (49)(a). One may decide to use ∇ * W (r J − r I ,h IJ ) in the discretization scheme (27) which leads to the error with the leading term This leads to the incorporation of the correction factor into Brookshaw's approximation to reduce the impact of the following matrix Γ αγ in (32) for the Laplacian approximation.
A different correction factor has been introduced and investigated in [50]. Alternatively, as was discussed in [50], the correction multiplier can be introduced in (27) Fig. 1 shows Laplacian values for the function ∇ 2 (x 2 + y 2 ) using original Brookshaw's approximation (27) with (a) conventional kernel ∇ γ W , (b) corrected kernel ∇ * γ W and corrected Brookshaw's approximation (51) with the conventional kernel ∇ γ W . In addition, Fig. 1 shows that the discretization scheme (27) with the ∇ * W (r J − r I ,h IJ ) kernel gradient is less accurate near the boundary due to the magnitude of the (50) which is related to the associated leading term of the error.
The scheme (27) is widely (almost unconditionally) used in the SPH modeling community. For example, it was used for a thermal conduction [53,66,67], for modeling a viscous diffusion [45], for a vortex spin-down [68] and Rayleigh-Taylor instability, for simulating Newtonian and non-Newtonian flows with a free surface [69] for the comparison of weakly compressible and truly incompressible algorithms, for macroscopic and mesoscopic flows [70], for a simulation of a solid-fluid mixture flow [71]. Recently, it has been used to model electrokinetic flows [72], Dam-break problem and Taylor-Green vortex [36].
There are different numerical SPH schemes used in numerical simulations. High order accuracy approximations can also be derived by using the SPH discretization on the higher order Taylor series expansion [55,15,50,73]. However, it is usually required that the discrete numerical schemes can reproduce linear fields [59,74,39,40] or polynomials up to a given order [75,14].

Schwaiger's scheme (2008)
The correction terms to the Brookshaw formulation which improve the accuracy of the Laplacian operator near boundaries were proposed by Schwaiger in [50]: n where n = 1, 2, 3 is the spatial dimension and the tensor Γ αβ should be defined in general as follows (differently to the expression (33)) The gradient ⟨∇ α u (r I )⟩ is the corrected gradient which can reproduce linear fields [59,74]. For multi-dimensional problems, the correction tensor Γ αβ (r I ) is a matrix. If the particle r I has entire stencil support (i.e., the domain support for all kernels W is entire and symmetric) then Γ αβ (r I ) ≈ δ αβ . Lemma 1. The following identities are valid for ∇ γ W ( r J − r I ,h IJ ) described by (15) and (17) and (18): Proof. These identities can be established by applying chain rule differentiation. □ Proposition 1. The correction tensors Γ (described by (55)) and C −1 (where C is defined by (24)) are the symmetrical and identical tensors for ∇ γ W described by (15) and (17) and (18).
Proof. The following relation can be established (using Lemma 1): It follows that the scheme (52)-(54) requires to use the trace of Γ −1 αβ (r I ) or C αβ (r I ). The trace of Γ αβ (r I ) provides to us the sum of all eigenvalues. It does not give us any range for the eigenvalues which is necessary to estimates C αα (r I ).
Even if we have a positive trace we can still theoretically have several small and large eigenvalues with an opposite sign. Before we get to the analysis of the eigenvalues, it is convenient to introduce a condition for matrices known as strictly diagonally dominant (SDD). The diagonal elements of a matrix are noted by indices {αα}, 1 ≤ α ≤ n.
Hence, the following Lemma can be formulated.  (15) and (17) and (18) with G(z J , r I ) ≤ 0 ∀z J and Γ ββ (r I ) ≤ ∑ Proof. It follows from the fact that: where the summation by repeating Greek indices is assumed throughout this paper. In order to figure out what range the eigenvalues of the strictly diagonally dominant matrix Γ αβ (r I ) (i.e., - would be in, we use Gershgorin's theorem: Remark 3. Since strictly diagonally dominant matrices are always nonsingular, the correction matrix C αβ (r I ) can be defined and C αα ( The correction matrix Γ αβ is strictly diagonally dominant for the homogeneous and symmetrical particles distribution subject to the full compact support of the kernel. In addition, it is important to note that in the case of using only the corrected gradient ∇ * β W and, hence, Γ −1 αα n = 1. However, ∇ * γ W does not satisfy relation (56).
For multi-dimensional problems, the correction tensor Γ αβ (r I ) is a matrix. If the particle r I has entire stencil support (i.e., the domain support for all kernels W is entire and symmetric) then Γ αβ (r I ) ≈ δ αβ . Unfortunately, Γ αβ (r I ) deviates from δ αβ for the provided algorithm and, hence, it is important to minimize this deviation from δ αβ in the new methods.
Finding neighbors of a particle in SPH is computationally costly major task in meshless methods. Hence, this task has a significant impact on the overall computational cost and specifically in dynamics simulations involving repositioning of the particles with the time [76]. The typical complexity of the nearest neighbor search algorithm, which can be used in any particle method is of O (N log N) [77]. To calculate coefficients in the scheme (52)-(55) is a trivial task and it is O (N). However, in general, it should be performed at each Newton-Raphson iteration in the nonlinear case (i.e., m = m (u (r I ))). It also requires additional efforts to compute the correction matrices C αβ and Γ αβ (inversion of n × n matrices per each particle, where n = 1, 2, 3 is the spatial dimension) when they are different (please see Proposition 1) and storage cost of ∇ α W ( , and corresponding Γ −1 αα per each particle. Theorem 2. The discretization scheme (52)-(55) at a given point r I ∈ Ω ⊂ R n sufficiently far away from the boundary ∂Ω (i.e., where the kernel W has a full support and particle volumes satisfy the SPH continuity condition (21)) has an error (discretization error) bounded by C Proof. Using Lemma 1 and Proposition 1, additional terms proposed by Schwaiger [50] reduce to in the vicinity of the point r I with the errorẼ. Where The expression (62) has the leading term outlined in (32). Subtraction of these terms from (32) leads to where and performing the normalization lead to the announced in the theorem result. Note that in order to maintain the results of the above theorem close to the boundary points additional ghost particles must be added outside the numerical domain. □ It follows that the usage of Γ * αβ instead of Γ αβ in (52) will improve the accuracy of the scheme (52)- (55). Furthermore, if one uses ∇W (r J − r I ,h IJ ) in place of ∇ * W (r J − r I ,h IJ ) in the first term of the discretization scheme (52) then the definition for N α has to be modified in accordance of (50) to maintain the higher order discretization accuracy, for example, as: which reduces to the conventional N α in the case of ∇W (r J − r I ,h IJ ) ̸ = ∇ * W (r J − r I ,h IJ ) and subject to the conditions of Lemma 1. Fig. 2 shows the Laplacian for the function ∇ 2 (x 2 + y 2 ) using Schwaiger's approximation (52)-(55) with (a) conventional kernel ∇ γ W and (b) corrected kernel ∇ * γ W .

Meshless transmissibilities
The well-known two-point flux approximation (TPFA) is a mesh dependent numerical scheme used in solving elliptic equation (1) whereT JI is the transmissibility between cells J and I, q is the total flux through the boundary of the control volume located at the point r I . The transmissibilityT JI defined at an interior face f between cells J and I is calculated as order of accuracy for any mobility tensor field M. The expression (67) ensures that the flux into the adjoining region is continuous [66]. The TPFA scheme (66) is an unconditionally monotone scheme [52].
It is clear that the expression (52) cannot be written in the form (66) due to terms ⟨∇ α (m (r I ) u (r I ))⟩N α and u (r I ) ⟨∇ α m (r I )⟩N α . Hence, it is only possible in this case to introduce a definition of a partial meshless transmissibility between particles r J and r I as follows: It is important to note that transmissibilities T P JI andT P JI have different physical units. Furthermore, it raises the question wherever the scheme (52)-(55) is monotone.
Hence, there is an additional term that has not been taken into account in (52)- (55). The following section describes an alternative numerical scheme for the heterogeneous Laplace operator. Some initial attempts were also made in [78].

New scheme
The correction terms to the Brookshaw [53] and Schwaiger [50] formulations which improve the accuracy of the Laplacian operator near boundaries can be done as follows: n where n = 1, 2, 3 is the spatial dimension. The tensorΓ αβ can be defined by the relation (55). It is important to stress at this point that it is also applicable to useΓ αβ = Γ αβ without losing the monotonicity and stability quality of the scheme. However, the higher accuracy can be achieved by taking into account higher order corrections in Γ αβ (r I ) as per relation (64). Hence, the following definition can be used and Following (69), we only need to compute the trace of the matrixΓ −1 αβ (r I ). Furthermore, it is clear that one can apply results of Lemma 2 to Γ * ββ (r I ) with following function The function G * (z J , r I ) is negative if The reason for having the correction factor in the form (70)-(71) is that Γ * αβ (r I ) = 0 in some cases, where particles have the incomplete kernel support (e.g., at the corners and boundaries of the numerical domain). For multi-dimensional problems, the correction tensorΓ αβ (r I ) is also a matrix. If the particle r I has entire stencil support (i.e., the domain support for all kernels W ( r J − r I ,h IJ ) is completed and symmetric) thenΓ αβ (r I ) ≈ δ αβ . The proposed correction matrix deviates less from the unit matrix compare to (55).
The computation complexity of the proposed scheme is very similar to the scheme (52)- (55). Although, the overall complexity over considered in this paper SPH schemes can be bounded by O (N log N) due to nearest neighbor search algorithm, in all problems consider in this paper the nearest neighbor search algorithm is performed once at the preprocessing step. The rest of the computations has the complexity of O (N). We do not consider here the complexity of the linear solver required to solve resulting algebraic system of equations.
Theorem 3. The discretization scheme (69)-(71) at a given point r I ∈ Ω ⊂ R n sufficiently far away from the boundary ∂Ω (i.e., where the kernel W has a full support and particle volumes satisfy the SPH continuity condition (21)) has an error (discretization error) bounded by C Proof. Using Lemma 1, additional term proposed by new scheme (69) subject to the diagonal mobility tensor M αβ (r) = m (r) δ αβ reduces to in the vicinity of the point r I with the errorÊ which has the leading term outlined in (32). Subtraction of these terms from (32) and performing the normalization lead to the announced in the theorem result. Note that in order to maintain the results of the above theorem close to the boundary points, additional ghost particles must be added outside the numerical domain. □ The scheme has the two-point flux approximation nature and can be written in the form of (2), which can be proved using the arguments above. The scheme (69)-(71) is in line with an alternative formulation for continuum mechanics called the peridynamic model [79,80], which was proposed several years ago.
All presented schemes in this paper do not require exact expressions for the gradient (i.e., spatial derivatives) of the mobility field ∇ γ m (r) to keep a higher order of accuracy for any mobility field. Hence, this scheme can be used with the discontinuous (or piecewise continuous) mobility field m (r) ∈ L 2 (Ω). It is important to note that Brookshaw [53] and Schwaiger [50] schemes can also be written for the diagonal mobility matrix M αβ (r) by substituting M αβ into (27) and (52) instead of m (r) and performing summation by repeating indices.

Approximation, stability, and monotonicity
The approximation, stability and monotonicity are important properties of numerical schemes which provide and quantify the confidence in the numerical modeling and results from corresponding simulations. Therefore, in order to be confident that the proposed numerical schemes provide the adequate accuracy of the elliptic operator (1), several numerical analyses to identify the order of approximation, stability and monotonicity have been performed. It is important to recall that all numerical schemes are characterized by two length scalesh and ∆x m :h = f · ∆x m . Hence, while investigating the approximation of the meshless discretization scheme, it is important to distinguish two cases: (a) the neighborhood number of particles is fixed f = const with varying the inter-particle distance, (b) the inter-particle distance is fixed ∆x m = const with varying the neighborhood number of particles. The first case will be analyzed by looking at the error defined by where L [·] is the analytical Laplacian at the particle r J and ⟨L⟩ [·] is the approximation of the Laplacian at the particle r J , ∥E∥ L 2 Ω is the averaged error over the entire domain. In the following paragraph, the numerical analysis is performed for various functions, particle distributions, and media properties. Following the work [50], the ability of the discretization to reproduce the Laplacian was tested for several functions in R n , n = 1, 2, 3: for selecting these testing polynomials is as follows. Since the functional space L p (Ω) is separable, the above polynomials form the everywhere dense subset of the L p (Ω). Hence, any function u ∈ L p can be approximated in L p using linear combination of the above polynomials leading to the following relations: The approximation error produced by the discretization schemes and considered in this paper for the polynomials u s (x) and u s m (x) gives the information about the error growth for the arbitrary function u(r). In each test, the homogeneous and heterogeneous particle distribution varying the smoothing lengthh and inter-particle distance ∆x m are used to study the approximation properties of the proposed discretization scheme.

Isotropic homogeneous media
Note that in the limit of the homogeneous isotropic media (i.e., M αβ (r) = m (r) δ αβ , m (r) = 1) and without the source term, the aforementioned operator L (u) in (1) reduces to the conventional Laplacian operator (i.e., L (u) ≡ ∇ 2 u). The domain for the patch test is a unit square similar to [50] for n = 2: In all of the following tests, the results are displayed along the cross-section y = 2.5 in R n , n = 2 (similar to [50] for n = 2). In each test, the following discretizations are compared: corrected Brookshaw's scheme (CB-SPH) (51), Schwaiger's scheme (S-SPH) (52)- (55), and new proposed scheme (M-SPH) (69)- (71). Note that the Schwaiger's scheme was tested against several schemes published in [62,49] and show better accuracy; hence, schemes in [62,49] are not considered in this paper. The comparison of different schemes starts with test functions of the form u s (x) described by (79)(a) in R n , n = 2. Plots of the Laplacian approximations and the relative errors defined by (78) for the case m = 3 are shown in Fig. 3. Note that it is important for the accurate solution of different boundary value problems to have an accurate approximation of the considered Laplace operator in the internal particles only. The additional relation will be supplied to approximate boundary conditions in the following section, e.g. (100).
The new scheme (M-SPH) has the greatest accuracy at the boundary and is accurate to machine precision ε in the interior. The same behavior is observed in 3D where the M-SPH scheme with the ∇ α W is the most accurate scheme. Furthermore, the new scheme should be tested for the functions requiring cross-derivatives as was reported in [50]. To examine the effect of the cross-derivative terms, the same suite of tests was run with the function u s m (x) described by (79)(b) in R n , n = 2, 3. Relative errors for (xy) m on a same array as in Fig. 3 along y = 2.5 were computed and again the proposed scheme (69)-(71) was uniformly more accurate for the various experiments. The behavior of each discretization is similar to that shown in Fig. 3. The proposed new scheme (M-SPH) performed nearly as well as the Schwaiger's scheme at the boundary. The CB-SPH and S-SPH forms also perform with greater accuracy than all other forms in the interior except in the case with the highest exponent.
An additional concern is that although it deviates from the exact solution near boundaries, it acquires no off-diagonal terms due to the alignment of the array of particles and the boundaries with the coordinate axes. To test the accuracy of the new approximations when there are off-diagonal terms, an array with particles rotated 45 • was used with the test function (xy) m (similar to [50]). The new proposed scheme formulation performs consistently well for lower exponents.

Isotropic heterogeneous media
In isotropic heterogeneous media, i.e., M αβ (r) = m (r) δ αβ , the Laplace operator takes the general form written in (1). The numerical domain is the same as in the previous section with the same number of particles and particle length scales:h = f · ∆x m , ∆x m = 0.05, f = 1.2. Plots of the Laplacian approximations and the relative errors defined by (78) for the case of heterogeneous mobility with m = 1 (see (80)) are shown in Fig. 4. The new scheme (M-SPH) has the best accuracy at the boundary and is accurate to machine precision ε in the interior. The error has also been computed for test functions with exponents from m = 2, . . . , 5. For these functions, the proposed scheme (69)-(71) is uniformly more accurate for the different increasing exponents. The same behavior is observed in 3D where M-SPH scheme with the ∇ α W is the most accurate scheme.

von Neumann stability analysis
For linear PDEs, there is the Lax-equivalence theorem which connects the consistency and stability with the convergence. The idea of the von Neumann stability analysis is to study the growth of waves λe ik·r (similar to Fourier methods). After applying one of the discretization methods above to the Laplace operator (1), the following relation can be written: where τ is the iteration parameter (e.g., time step), n is the iteration index (e.g., time index), {u} n = ( leads to the expression for the von Neumann growth factor subject to (81) and linear Laplace operator: , r ∈ R n , n = 1, 2, 3; where { e ik·r } = ( e ik 1 ·r 1 , . . . , e ik N ·r N ) and λ j (τ ) is the von Neumann growth factor. For the discretization to be stable, it is required that In case of pseudo random particle distribution, the von Neumann growth factor has both real and imaginary parts forming complex shape but satisfying the requirement ⏐ ⏐ λ j (τ ) ⏐ ⏐ ≤ 1, ∀j almost everywhere (i.e., it could be some problems at the boundary particles).

Numerical stability analysis
The results of von Neumann stability analysis are validated using the following 2D numerical example considering a laminar flow of an incompressible (∂v x /∂x + ∂v y /∂y = 0) Newtonian fluid through a parallel plate slit. The governing equations in this case are given by: where w = p/ρ 0 is the specific (per unit mass) thermodynamic work, ρ 0 is the density, and t is the time, v x and v y are the velocity in the x-longitudinal and y-transverse directions of a parallel plate slit, p is the pressure. Eq. (84) is written in dimensionless form. The dimensionless Reynolds number is defined as Re = (U · L · ρ) /µ, where U is the characteristic velocity, L is the characteristic length, and µ is the dynamic viscosity of the fluid. We consider a square of 2 by 2 in R 2 . The homogeneous particle distribution is used to discretize the computational domain. The initial and boundary conditions for v = {v x , v y } required to solve (84) are given by: and v(t, y = ±1) = 0; (85) where v * and v * * are prescribed velocities as some functions of time. The specific thermodynamic work is defined up to a constant.
We consider flow through a truncated section of a pipe with prescribed specific (per unit mass) thermodynamic work ∆w/L gradient governed by: The solution to (86) subject to the conditions in (85) and infinite in the longitudinal direction pipe yields to the following analytical velocity profile: where ∆w/L is the imposed gradient of the specific thermodynamic work and w(x = 0) = w 1 , w(x = 2) = w 2 , ∆w = w 2 − w 1 . We will use this v a x profile to set up the boundary conditions for inlet and outlet particles velocity, which are placed at the boundary of the computational domain.
The temporal discretization is as follows v n+1 where v n+1 , v n+1 y } and v n = {v n x , v n y } are the divergence-free velocities at time n + 1 and n, respectively and L () is the meshless discretization of the Laplace operator using M-SPH method. We introduce the maximum error at time step n as follows: where J is the particle index and N is the number of particles. As a result, the M-SPH spatial discretization is able to correctly simulate flow in a parallel plate slit. As an example, Table 1 shows the error log 10 E n decay for the explicit method (88) during the numerical simulation obtained using the M-SPH method (69)-(71) and different characteristic diffusivity ( Re ·h 2

I,m
) /∆t with Re = 1, ∆w/L = 1. The error curve is interpolated by using cubic splines to obtain values at the selected times presented in Table 1. From the previous section, it follows that the stability requirement is met for ν ≤ 0.25. Table 1 shows that the solution is stable under this condition and has smaller error around ν = 0.2 approaching steady-state solution in the optimum way. The solution is strictly speaking unstable for ν > 0.25 and as it shown in the last column of Table 1.

Monotonicity and convergence
In real life applications, the numerical domains are large and, therefore, many degree of freedom (i.e., unknowns) is required. The resulting matrix is usually sparse, but because of fill-in a direct method requires significant amount of memory and time in general case. Hence, iterative solvers are used to overcome these issues. The iterative solvers converge only if the matrix satisfies certain properties related to the property of the diagonal dominance. However, these properties are only sufficient conditions for the method to converge.
The proof of the convergence of linear meshless schemes applied to a linear elliptic boundary value problem can be done in the following steps (see, Bouchon (2007) [81], Bouchon and Peichl (2007) [82], Matsunaga and Yamamoto (2000) [83], Thomée (2001) [84]: if it is shown that the truncation error ε tends to 0 as the maximum smoothing length h max = max I∈NhI goes to 0, then the linear system Aδu = ε that couples the variable error δu with ε proves the convergence of the schemes, given that the matrix of the discretized meshless operator is monotone.
The M-matrix is monotone since its inverse is positive, i.e. A −1 ≥ 0. Furthermore, the monotone schemes do satisfy a discrete maximum principle producing solutions without spurious oscillations. It is clear that the new scheme (69)- (71) can be written in the form (66), where meshless transmissibility between particles r J and r I can be defined as T IJ with: whereÑ α (r I ) is defined by (53). Applying the relation (19) valid for kernel gradients (15), (17), (18) and the definition of the corrected kernel gradient (23), (24) leading tō The sign of the transmissibility coefficient is described by the following Lemma.
where sgn is the sign function.
where ξ β I is the vectorÑ α (r I ) C (r I ) αβ , ∆r JI is the vector

Solution of boundary value problems
Following the work by [50], the numerical tests for inhomogeneous Dirichlet, Neumann and mixed boundary conditions are considered for homogeneous and heterogeneous media with the characteristics M (r). To illustrate the performance of the proposed scheme, a modeling of a single phase steady-state fluid flow in fully anisotropic porous media with different type of boundary conditions is also presented in this section. The square 2D and 3D domains (3) are considered. The relative error used to quantify the accuracy of the proposed schemes in the subsections 2.3, 2.4, and 3.1 during numerical simulation in this section is given by: where u (r) is the analytical or reference solution field and ⟨u (r K )⟩ is the approximated solution field. Several numerical results using considered in this paper schemes for uniform and pseudo random particle distributions are shown in the following sections that confirm the theoretical results from the previous sections.

Inhomogeneous boundary condition
Firstly, the homogeneous properties of the porous media M (r) = I are assumed for simplicity of the derivation. The analytical solution of (1) subject to the assumption that g (r) ≡ 0 for r ∈ Ω ⊂ R 2 and constant boundary conditions is the following: where ψ 1 is the boundary at y = 0, ψ 2 at x = L, ψ 3 at y = H, and ψ 4 at x = 0. Tables 2, 3 show the convergence rate of different schemes for uniform particle distribution for different number of degrees of freedom (DoF). The solutions with MPFA and Mimetic schemes were obtained using MATLAB Reservoir Simulation Toolbox (MRST) [85]. It is clear from the convergence results that schemes (52)- (55) and (69)-(71) are identical and they over perform the scheme (51). This is due to the fact that a regular grid of particles does not activate the Table 2 The error of convergence for different schemes (uniform particle distribution) and f = 0.5005.  Table 3 The error of convergence for different schemes (uniformed particle distribution) and f = 1.001.
The convergence rate as was predicted theoretically is at least O (h ω ) , 1 ≤ ω ≤ 2. Tables 4, 5 shows the convergence rate for different schemes and pseudo random distribution. The reason for starting from f = 0.5005 is that we want to force the same number of non-zero entries in the matrix as for MPFA or Memetic methods. This procedure may lead to a badly conditioned matrix since This distribution of particles was generated by perturbing the regularly distributed particles using a uniform random variable varying between +10% and −10% of the maximum smoothing lengthh max = 1.001. Since random realizations are used to compute errors, the statistical data about 30 realizations are used to compute mean values and standard deviations. This data is presented in Tables 4, Table 5. The positive trends of the scheme remain the same as for the case of uniform particle distributions. The scheme (51) has higher error compare to schemes (52)-(55) and (69)- (71). However, the dispersion of the approximation error is higher in the scheme (69)- (71). Fig. 7 shows the numerical solution of the boundary value problem obtained using different discretization methods. The solution was obtained using 40 particles/cells in each direction.
It is important to stress that the new scheme is not always outperforming Schwaiger scheme (52)-(55) for random particle distribution. In Table 5 the new scheme shows higher accuracy, compared with Schwaiger scheme. In order to gain insight into the error distribution, we should notice from 1D analysis that leading term of the correction U * ωβα · N ω term of tensor Γ αβ has the following asymptotic: and having independent of smoothing length term and, hence, should improve the accuracy if we take this into account while computing the correction factor. On the other hand, we cannot control the error introduced to the tensor Q αγ (essential to the magnitude of the error) defined by (47), which can be sometimes larger for new scheme.

Inhomogeneous mixed boundary condition test
The general steady-state solution for a Dirichlet condition along the base of the plate and Neumann conditions elsewhere is given in [50] by where ψ 1 is the boundary value of u at y = 0, ψ 2 is the flux at x = L, ψ 3 is the flux at y = H, and ψ 4 is the flux at The following parameters (dimensionless) were chosen: ψ 1 = 150, ψ 3 = 150, ψ 4 = 200, ψ 2 = 90. Again, for uniform and regular particles distribution, MPFA method reduces to the TPFA. Similar to the previous section, the regular and pseudo random particles distributions are considered to access the convergence properties of the proposed schemes.   Table 7 The error of convergence for different schemes (uniform particle distribution) and f = 1.001.
The following smoothing multiplication factors f = 0.5005, f = 1.001, and f = 2.002 are considered. Tables 6,   Table 7 show the convergence rate for different schemes using the uniform and pseudo random particle distribution. The convergence rate as was predicted theoretically is at least O (h ω ) , 1 ≤ ω < 2. Interestingly, the scheme (51) shows a very bad convergence rate for f = 0.5005. This explains by the fact that we have one particle from each side in the kernel support almost next to the boundary of the kernel support leading to   ∇W ( r J − r I , h )  ≪ 1 and, hence, to a badly conditioned matrix. This cannot be observed for other schemes due to the normalization coefficient. As a result, the scheme (51) should be used with more than one particle in the compact support in each direction leading to a larger bandwidth in the matrix. Tables 8, 9 show the convergence rate for different schemes with pseudo random particle distribution. The pseudo random distribution of particles was again generated by perturbing the regularly distributed particles using the uniform random variable varying between +10% and −10% of the maximum smoothing lengthh max = 1.001. Similar to the above case, the statistical data about 30 realizations are used to compute mean values and standard deviations. These data are presented in Tables 8, 9. The general trends of the scheme prosperities remain the same as for uniform particle distributions. The scheme (51) has higher error compared to schemes (52)-(55) and (69)- (71). However, the dispersion of the approximation error is higher in the scheme (69)- (71). Results of the numerical solutions are shown in Fig. 8.

Remark 5.
The solution of boundary value problems with inhomogeneous Dirichlet only and mixed boundary conditions shows that new scheme has identical accuracy with Schwaiger (2008) scheme (52)-(55) for uniform particle distribution as expected theoretically due tõ and a regular grid of particles does not activate the additional correction terms in both schemes. However, for random particle distributions subject to the ∇ γ W ( r J − r I ,h IJ ) described by (15), (17) and (18), we havẽ Finally, we must stress that advantages of the new scheme are related to desirable properties (i.e., satisfying the twopoint flux approximation nature in the form (2), stability and proven monotonicity but not the accuracy which has the same order as Schwaiger (2008) scheme (52)- (55) in general and better in some cases. Furthermore, desirable properties of the new scheme will allow to apply upwinding strategy during the solution of nonlinear PDEs.

SPE10
In addition, we investigate the accuracy of the numerical schemes using a well-known SPE10 benchmark [86] with the Laplace equation: L (u) = −∇ (M (r) ∇u (r)) = 0, ∀r ∈ Ω ⊂ R n , M αβ (r) = K αβ (r) , K αβ (r) = 0 ∀α ̸ = β; α, β = 1, . . . , n; (101) where u (r) is the unknown pressure field, K αβ (r) is the diagonal permeability field. The original model contains 85 layers, where layers from 1 to 35 have smooth permeability with lognormal distribution and layers 36 to 85 have channelized formations that are considered to be significantly more challenging for numerical simulations. The subset of this model is defined by the global Cartesian indices I, J, K . The 85 layer was used as a permeability field for 2D simulation with Cartesian indices I = 1 : 60 and J = 1 : 60. In 3D, the subsection Cartesian indices I = 1 : 60, J = 1 : 60, and K = 1 : 60.
The permeability fields K 11 for both cases are shown in Fig. 9. The boundary conditions correspond to the unit pressure drop over the entire domain in J-direction (i.e., y min = 0 and y max = 1). The numerical results using new scheme for the SPE10 cases are presented in Fig. 10. The relative error distribution is also shown for 2D and 3D cases, where the error is computed using (95) and numerical solution based on TPFA. Figs. 11 and 12 compare the convergence rates of the various discretizations subject to different preconditioners. Furthermore, if we notice the convergence comparison for this test in    55)). The reason for such convergence is the similarity between condition numbers of the linear system of equations resulting from these methods.  55)). The linear system of equations is built using 3D case.
Another interesting observation in this numerical test is the small values of the relative error between the TPFA solution and proposed new method for both 2D and 3D cases. This observation exhibits the efficiency of the proposed meshless discretization scheme, which explains the higher accuracy of linear reproduction.

Conclusion
In this paper, the new stable SPH discretization of the elliptic operator for heterogeneous media is proposed. The scheme has the two-point flux approximation nature and can be written in the form of (2). Using this structure, it was possible to make some theoretical monotonicity analysis, which is difficult to perform for other schemes (e.g., Schwaiger's method). The paper also discusses the monotonicity and convergence properties of the new proposed scheme. The sufficient monotonicity condition is also formulated leading to the constraint on the kernel shape and particles distribution (see, Lemma 3 and Theorem 4). Furthermore, it follows from the Taylor's series analysis that proposed scheme is the optimum (i.e., high order with monotonicity property) one in this class (2) for a diagonal matrix of the operator coefficients. In addition, the proposed scheme allows to apply upwinding strategy during the solution of nonlinear PDEs.
The new scheme is based on a gradient approximation commonly used in thermal, viscous, and pressure projection problems and can be extended to include higher-order terms in the appropriate Taylor series. The proposed new scheme is combined with mixed corrections which ensure linear completeness. The mixed correction utilizes Shepard Functions in combination with a correction to derivative approximations. Incompleteness of the kernel support combined with the lack of consistency of the kernel interpolation in conventional meshless method results in fuzzy boundaries. In the presented meshless method, the domain boundary conditions and internal field variables are approximated with the default accuracy of the method. The resulting new scheme offers a higher order of accuracy, but also minimize the impact of the particle deficiency (kernel support incompleteness) problem. Furthermore, different kernel gradients and their impact on the property of the scheme and accuracy are discussed. The model was tested by solving an inhomogeneous Dirichlet and mixed boundary value problems for the Laplacian equation with good accuracy confirming our theoretical results. The accuracy of Schwaiger's scheme and new scheme is the same for homogeneous particle distribution and different for the distorted particles. The stability analysis shows that von Neumann growth factor has both real and imaginary parts forming complex shape for general particle distribution but satisfying the stability requirement almost everywhere.
As was previously mentioned in the introduction, several methods have been proposed to address the difficulties involved in calculating second-order derivatives with SPH. In contrast to the present formulation, many of these methods achieve a high accuracy through fully calculating the Hessian or requiring that the discrete equations exactly reproduce quadratic functions. The primary attraction of the present method is that it provides a weak formulation for Darcy's law which can be of use in further development of meshless methods. The SPH model was previously used to model threedimensional miscible flow and transport in porous media with complex geometry, and we are planning to use this model in future work for large (field) scale simulation of transport in porous media with general permeability distributions. Hybrid correction tensor to the SPH Hessian discretization