- Email: [email protected]
A shock-capturing meshless scheme using RBF blended interpolation and moving least squares
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A shock-capturing meshless scheme using RBF blended interpolation and moving least squares Michael F. Harris a,b,∗, Alain J. Kassab a, Eduardo Divo b a b
Department of Mechanical and Aerospace Engineering, University of Central Florida, Orlando, FL, USA Department of Mechanical Engineering, Embry-Riddle Aeronautical University, Daytona Beach, FL, USA
a r t i c l e Keywords: Meshless methods Radial basis functions Moving least squares Shock-capturing
i n f o
a b s t r a c t Solving the full Navier-Stokes equations requires the use of numerical methods, especially for engineering application involving complex geometry. Numerical methods commonly used today all require a mesh before a solution attempt. Developing the mesh can require expensive preprocessing time, and the quality of the mesh can have major effects on the solution. Meshless methods have become a research interest due to the simplicity of using scattered data points. Radial basis functions (RBF) interpolation is a class of meshless method, and many of these methods have been developed to solve partial differential equations. Schemes using RBF interpolation are capable of multivariate interpolation over scattered data even for data with discontinuities. A meshless method scheme is presented capable of solving the compressible Navier-Stokes equations. A blended RBF interpolation method switching between low and high shape parameter interpolation is used to approximate the inviscid fluxes, while the viscous fluxes are approximated using moving least squares (MLS).
1. Introduction Solutions to nonlinear partial differential equations over complex geometry in most engineering applications require the use of numerical methods. The most common numerical methods used today include finite difference method (FDM), finite volume method, and finite element method (FEM), but these methods require a preprocessing step of developing a mesh. Research in numerical methods has focused on reducing the computational time for the solution while mesh generation remained a bottleneck. Mesh development can be time-consuming for general problems and even more problematic for compressible flow since shock locations are unknown prior to the solution. Mesh adapting methods have been developed to overcome this issue requiring an iterative process of solving and re-meshing capturing the sharp gradients. Generating a high-quality mesh while balancing computation resources while capturing the physics can be difficult. Meshless methods have been a research interest for many years and shown to have advantages over other numerical methods. Many types of meshless methods developed over the years can be used to solve the governing equations found in fluid dynamics, heat transfer, and structural mechanics [1]. Methods such as the Trefftz method [2], the boundary knot method [3], the method of fundamental solutions [4], the method of fundamental solutions in combination with radial basis functions [5,6], and the singular boundary method [7]. Other
∗
methods have been investigated, such as element free Galerkin methods [8–11], Hp cloud method [12] and partition of unity [13]. Among these meshless methods are an important group related to applications of radial basis functions, such as the Kansa Method [14–18]. These numerical methods use a unifying concept of interpolation, but the type of interpolation is different for each method. Meshless methods began from spectral or pseudo-spectral methods based on Legendre and Chebyshev polynomials requiring the point or nodal distribution to be uniform. The use of radial basis functions (RBF) allows the use of irregular spatial distribution of points or nodes [19–21]. Radial basis functions have been a subject of research for decades as a global and local collocation method in solving partial differential equations. The Hardy Multiquadrics RBF was developed for topographical mapping and has proven to be an impressive function for multivariate interpolation [22,23]. Franke determined that the Hardy Multiquadric RBF produces stable results and yields the most accurate results of all other methods studied [24]. Kansa developed the RBF collocation meshless method using Hardy Multiquadric to solve partial differential equations such as the linear advection-diffusion equation and the Poisson equation [14,15]. These numerical experiments demonstrated the efficiency of using RBF multiquadric in contrast to the finite difference method. In the Kansa method, the RBFs are applied globally in strong form to solve partial differential equations, but the method is difficult to manage when using large datasets due to
Corresponding author. E-mail addresses: [email protected] (M.F. Harris), [email protected] (A.J. Kassab), [email protected] (E. Divo).
https://doi.org/10.1016/j.enganabound.2019.08.019 Received 31 March 2019; Received in revised form 21 July 2019; Accepted 6 August 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Fig. 1. Hardy Multiquadrics for varying shape parameter, c.
Table 1 Radial basis functions. Name
Expression
Hardy Multiquadric (MQ) Inverse Hardy Multiquadric (IMQ)
𝜒(𝑟) = (𝑟 + 𝑐 ) 1 𝜒(𝑟) = 1
Gaussian
𝜒(𝑟) =
2
Carlson and Foley observed that the shape parameter value is dependent on the function and is independent of the number and location of data points [48,49]. Their work demonstrates Hardy Multiquadric interpolation on topography that lie along the same plane or track and the distance between the data points can vary by orders of magnitude. As the shape parameter approached zero, they observed the Hardy Multiquadrics interpolation yielding better accuracy in those regions. Carlson and Foley recommended the shape parameters to be reduced enough to remove any oscillations caused by the interpolation due to high gradients in measurements. The suggestion of using a low shape parameter and therefore piecewise linear interpolant may not be favorable strategy since accuracy is lost, but using steep RBF interpolation can provide stability in specific regions of high gradient. Jung and Durante also recommended the notion of reducing the shape parameter to zero for the Multiquadrics RBF near local discontinuities [50]. Their algorithm lowers the shape parameter to zero based on a developed smoothness criterion. The algorithm switches the shape parameter to zero for each collocation point if a jump discontinuity is detected. This iterative procedure proposed by the Jung and Durante works very well, eliminating oscillations in the data and capturing the discontinuities. The Hardy multiquadrics and inverse multiquadrics are the most impressive functions as these radial basis functions produce the best accuracy. These functions depend on a shape parameter that must be chosen properly for accurate approximations. The research in RBF interpolation has focused heavily on overcoming the ill-conditioning of RBF interpolation always seeking flat RBF interpolation for all types of problems. This approach works well for smooth functions, but tends to cause oscillations for data where steep gradients or shocks exist.
Shape Parameter 2
(𝑟2 +𝑐 2 ) 2 2 𝑒−𝑐 𝑟
1 2
c>0 c>0 c>0
ill-conditioning of the interpolation matrix. Domain decomposition [19], matrix preconditioning [25,26] and local RBF collocation methods [20,27–34] have been used to overcome the difficulties of large data sets and ill-conditioned matrices. RBF interpolation using Hardy Multiquadric RBF has shown to provide high accuracy and exponential convergence [35–38]. However, the Hardy multiquadrics and inverse multiquadrics are dependent on the shape parameter, c. These radial basis functions are given in Table 1. The optimal shape parameter for a given problem can be difficult to determine, and many researchers have focused their attention on finding a solution for the shape parameter problem [17,18,37,39–43]. For best accuracy, the shape parameter should be chosen to be a large value causing the condition number of the interpolation matrix to be close to ill-conditioned, rendering the RBF flat [35–37,44,45]. Spectral convergence is reached as the shape parameter c → ∞ [35,46,47]. As c → 0, the Multiquadric RBF becomes steep, and the accuracy is piecewise linear. The multiquadric RBF for different shape parameter values is shown in Fig. 1. Although, low shape parameter can be used to add stability, especially in areas of steep gradients and shocks, high shape parameter interpolation must be used for accuracy. 82
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Reducing the shape parameter to render the RBF steep and the interpolation matrix condition number low eliminates oscillations better, capturing the steep gradients. A new approach for solving hyperbolic PDEs is presented using a localized blended RBF collocation method [51–54]. The blended RBF strategy is to combine the benefits of flat RBF and steep RBF interpolation to form a high-resolution scheme for shock capturing. 2. Radial basis function interpolation Radial basis function interpolation approximates the solution of a function by the expansion 𝑓 (𝑥 ) =
𝑁 ∑ 𝑗=1
𝛼𝑗 𝜒𝑗 (𝑥)
(1)
where 𝛼 j are the expansion coefficients, 𝜒 j (x) is the expansion function, and N is the number of collocation points [14,22,55]. The coefficients 𝛼 j are determined through the collocation of the function f(x) at the points f(xi ) for i = 1…N. 𝑁 ∑ ( ) ( ) 𝑓 𝑥𝑖 = 𝛼𝑗 𝜒𝑗 𝑥𝑖
(2)
{𝑓 } = [C]{α}
(3)
Fig. 2. Local RBF collocation meshless method.
𝑗=1
The collocation matrix [C] is of the form ( ) ( ) … 𝜒𝑁 𝑥 1 ⎤ ⎡ 𝜒1 𝑥 1 ⎥. [C ] = ⎢ ⋮ ⋱ ⋮ ( )⎥ ⎢ ( ) … 𝜒𝑁 𝑥 𝑁 ⎦ ⎣ 𝜒1 𝑥 𝑁
decomposition has been used to solve the difficulties of global methods. Another strategy is to use a local RBF collocation meshless method [19–21]. The localized RBF collocation meshless method (LRCMM) will be described in this section [55]. The LRCMM uses a localized RBF expansion to interpolate the field variables and approximate the derivatives. A set of data points sampled around the data center (xi ,yi ) as shown in Fig. 2. The function f(x, y) can be approximated by the expansion as shown in Eq (10). The expansion can be written in matrix-vector form given by Eq. (11) where 𝜒 j (x,y) is chosen to be the Hardy Multiquadric RBF.
(4)
The expansion function, 𝜒 j (x), is chosen to be the Hardy Multiquadric radial basis function given by ( )1 𝜒𝑗 (𝑥 ) = 𝑟 2 + 𝑐 2 2
(5)
The coefficients, 𝛼 j , are then found by solving the linear system given by Eq. (6) and the function can be reconstructed by recalling Eq. (2). {α} = [C] {𝑓 } −1
𝑓 (𝑥, 𝑦) =
(6)
𝑗=1
The RBF can be augmented by a polynomial expansion given by 𝑓 (𝑥 ) =
𝑁 ∑ 𝑗=1
𝛼𝑗 𝜒𝑗 (𝑥) +
𝑁𝑃 ∑ 𝑗=1
𝛼𝑗 𝑃𝑗 (𝑥)
𝛼𝑗 χ𝑗 (𝑥, 𝑦) +
𝑁𝑃 ∑ 𝑗=1
𝛼𝑗+𝑁 𝑃𝑗 (𝑥)
𝑓 (𝑥, 𝑦) = {𝐶 (𝑥, 𝑦)}𝑇 {α} (7) 𝜒𝑗 (𝑥, 𝑦) =
where N is the number of collocation points and NP is the number of polynomials [55]. The polynomials are Pj (x) can be for example P = {1, x}. For two-dimensional problems, P = {1, x, y}. This addition of the polynomial expansion provides stability in constant and linear fields, but the added polynomial expansion can reduce accuracy in some cases [48,49]. The collocation matrix and the vector of supporting node field values are of the form in Eq. (8) and Eq. (9). ( ) ( ) ( ) ( ) … 𝜒𝑁 𝑥 1 𝑃 1 𝑥 1 … 𝑃𝑁𝑃 𝑥1 ⎤ ⎡ 𝜒1 𝑥 1 ⎢ ⎥ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ( ) ( ) ( )⎥ ⎢ ( ) 𝜒1 𝑥 𝑁 … 𝜒𝑁 𝑥 𝑁 𝑃 1 𝑥 𝑁 … 𝑃𝑁𝑃 𝑥𝑁 ⎥ ⎢ ( ) ( ) [C ] = ⎢ 𝑃1 𝑥1 ⎥ 0 … 0 … 𝑃1 𝑥𝑁 ⎢ … ⎥ ⋱ … ⋮ ⋱ ⋮ ⎢ ⎥ ( ) ( ) ⎣𝑃𝑁𝑃 𝑥1 ⎦ 0 … 0 … 𝑃𝑁𝑃 𝑥𝑁 (𝑁 +𝑁 𝑃 )×(𝑁 +𝑁 𝑃 )
√ (
𝑥 − 𝑥𝑗
)2
(10)
(11)
( )2 + 𝑦 − 𝑦𝑗 + 𝑐 2
(12)
The shape parameter, c, must be chosen wisely to obtain the best accuracy. A linear system of equations are formed from the local support nodes near the data center (xi ,yi ). {𝑓 } = [C]{α}
(13)
{α} = [C]−1 {𝑓 }
(14)
Substituting Eq. (14) into Eq. (11), we can interpolate the field variable value anywhere within the domain by 𝑓 (𝑥, 𝑦) = {𝐶 (𝑥, 𝑦)}𝑇 [𝐶 ]−1 {𝑓 }.
(15)
The derivatives are approximated by applying the derivative operator over the local expansion as { }𝑇 𝜕𝑓 (𝑥, 𝑦) 𝜕𝐶 (𝑥, 𝑦) = (16) [𝐶 ]−1 {𝑓 } 𝜕𝑥 𝜕𝑥
(8) ( ) ⎧ 𝑓 𝑥1 ⎫ ⎪ ⎪ ⎪ ( ⋮ )⎪ ⎪𝑓 𝑥 𝑁 ⎪ {𝑓 } = ⎨ ⎬ ⎪ 0 ⎪ ⎪ ⋮ ⎪ ⎪ 0 ⎪ ⎩ ⎭(𝑁 +𝑁 𝑃 )× 1
𝑁 ∑
and 𝜕𝑓 (𝑥, 𝑦) = 𝜕𝑦
(9)
{
𝜕𝐶 (𝑥, 𝑦) 𝜕𝑦
}𝑇
[𝐶 ]−1 {𝑓 }.
(17)
The vectors {C(x, y)}T [C]−1 {f}, { 𝜕C(𝜕x𝑥,𝑦) }𝑇 [C]−1 and { 𝜕C(𝜕y𝑥,𝑦) }𝑇 [C]−1
can be pre-computed and saved for each data center prior to the computation. Therefore, the derivatives of field variables can be approximated by the product of { 𝜕C(𝜕x𝑥,𝑦) }𝑇 [C]−1 and the field variable values within
Early RBF collocation methods such as the Kansa method are global methods and can be difficult to use for large-scale problems. Domain 83
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
the subdomain. 𝑓 (𝑥, 𝑦) = {𝐶 (𝑥, 𝑦)}𝑇 [𝐶 ]−1 {𝑓 } 𝜕𝑓 (𝑥, 𝑦) = 𝜕𝑥 𝜕𝑓 (𝑥, 𝑦) = 𝜕𝑦
{
{
𝜕𝐶 (𝑥, 𝑦) 𝜕𝑥 𝜕𝐶 (𝑥, 𝑦) 𝜕𝑦
}𝑇 }𝑇
Similarly, the derivatives are approximated by differentiation of the RBF vector { { } } 𝜕 𝐶𝐿 (𝑥, 𝑦) 𝑇 [ ]−1 𝜕 𝐶𝐻 (𝑥, 𝑦) 𝑇 [ ]−1 𝜕𝑓 (𝑥, 𝑦) 𝐶𝐿 {𝑓 } +𝜙 𝐶𝐻 {𝑓 } = (1 − 𝜙) 𝜕𝑥 𝜕𝑥 𝜕𝑥 (33)
(18)
[𝐶 ]−1 {𝑓 }
(19)
[𝐶 ]−1 {𝑓 }
(20)
and 𝜕𝑓 (𝑥, 𝑦) = (1 − 𝜙) 𝜕𝑦
3. Blended RBF collocation method
𝑓𝐿 (𝑥, 𝑦) =
𝑗=1 𝑁 ∑ 𝑗=1
𝛼𝐻 𝑗 𝜒𝐻 𝑗 (𝑥, 𝑦) +
𝛼𝐿 𝑗 𝜒𝐿 𝑗 (𝑥, 𝑦) +
𝑁𝑃 ∑ 𝑗=1
𝑁𝑃 ∑ 𝑗=1
𝛼𝑗+𝑁 𝑃𝑗 (𝑥)
𝛼𝑗+𝑁 𝑃𝑗 (𝑥)
(22)
(23)
2
(
)
𝑓𝐻 𝑥𝑖 , 𝑦𝑖 = ( ) 𝑓𝐿 𝑥𝑖, 𝑦𝑖 =
𝑗=1 𝑁 ∑ 𝑗=1
(
𝛼𝐻 𝑗 𝐶𝐻 𝑗 𝑥𝑖 , 𝑦𝑖 (
𝛼𝐿 𝑗 𝐶𝐿 𝑗 𝑥𝑖 , 𝑦𝑖
)
)
2
𝑟𝑖−1 =
(29)
(36)
2
𝑢𝑖−1 − 𝑢𝑖−2 . 𝑢𝑖 − 𝑢𝑖−1
(38)
The blending parameter in regions where high gradients, shocks or discontinuities are present the minmod limiter reduces a value less between [0,1]. If the value of the successive gradient is less than zero the minmod limiter is equal to zero providing low shape parameter interpolation. Another approach taken that involves using the smoothness indicator developed for the weighted essentially non-oscillatory or WENO scheme [59]. The WENO scheme uses polynomials to interpolate the fluxes at the midpoints of the stencil. Many stencils are used in the formulation to provide high order accuracy, and the smoothness indicator determines the stencil that provides the smoothest flux approximation by comparing the L2 norm values of the derivatives of the polynomial used. This flux is then used to approximate the derivatives in the
Solving for the expansion coefficients by inverting the collocation matrices [CL ] and [CH ], we have the expressions { } [ ]−1 (30) α𝐻 = 𝐶𝐻 {𝑓 } and { } [ ]−1 α𝐿 = 𝐶𝐿 {𝑓 }.
2
The local RBF collocation method gives the versatility of interpolating the field variable for the successive gradients. The minmod limiter is used for the blending parameter and is calculated by ( ) ( ( )) 𝜙 𝑟𝑖 = max 0, min 𝑟𝑖 , 1 (39) ( ) ( ( )) 𝜙 𝑟𝑖−1 = max 0, min 𝑟𝑖−1 , 1 (40)
(27)
[ ]{ } 𝐶𝐿 α𝐿 = {𝑓 }
2
and
(26)
(28)
2
2
( ) ( ) + 𝜙 𝑟𝑖−1 𝑓𝐻 𝑖− 1 − 𝑓𝐿 𝑖− 1
where the blending parameter is a function of the successive gradient r. The successive gradient is evaluated by the expressions 𝑢 − 𝑢𝑖−1 𝑟𝑖 = 𝑖 (37) 𝑢𝑖+1 − 𝑢𝑖
(25)
[ ]{ } C𝐻 α𝐻 = {𝑓 }
2
𝑓𝑖 − 1 = 𝑓𝐿 𝑖 − 1
Next, the weights 𝛼 Hj and 𝛼 Lj must be determined. A system of equations can be formed using all the field variable values within the local topology. 𝑁 ∑
} { [ ]−1 𝜕 𝐶𝐻 (𝑥, 𝑦) 𝑇 [ ]−1 𝐶𝐿 {𝑓 } +𝜙 𝐶𝐻 {𝑓 }. 𝜕𝑦 (34)
A method of determining how to blend between high and low shape parameter is needed. Fortunately, high-resolution schemes exist in computational fluid dynamics providing insight since these schemes switch between low and high order discretization for shock capturing. In CFD, flux or slope limiters are used to limit the gradients in the solver to realistic values. These flux/slope limiters are functions of the successive gradients found within the solution domain. Many flux limiters have been developed over the years [56–58]. For this research, a simple minmod flux limiter will be used. The minmod flux limiter can vary between [0,1] depending on the successive gradient r. For example, the flux is calculated by ) ( ( ) 𝑓𝑖 + 1 = 𝑓𝐿 𝑖 + 1 + 𝜙 𝑟𝑖 𝑓𝐻 𝑖 + 1 − 𝑓𝐿 𝑖 + 1 (35)
Defining two RBF with high and low valued shape parameters cL and cH . The Hardy Multiquadrics RBF is used as an example √ ( )2 ( )2 χ𝐿 𝑗 (𝑥, 𝑦) = (24) 𝑥 − 𝑥𝑗 + 𝑦 − 𝑦𝑗 + 𝑐𝐿2 √ ( )2 ( )2 2 χ𝐻 𝑗 (𝑥, 𝑦) = 𝑥 − 𝑥𝑗 + 𝑦 − 𝑦𝑗 + 𝑐𝐻
}𝑇
4. Blending parameter for shock capturing
where fH may be the smooth function interpolation and fL the steep gradient interpolation. 𝑁 ∑
𝜕 𝐶𝐿 (𝑥, 𝑦) 𝜕𝑦
The term 𝜑 is the blending term where 𝜑 = 1 will provide high shape parameter RBF interpolation and 𝜑 = 0 will give low shape parameter interpolation. A method is needed to determine the appropriate way of determining the blending parameter 𝜑. Methods using flux limiters or smoothness indicator, such a total variation diminishing (TVD) and weighted essentially non-oscillatory (WENO) schemes are good candidates for numerical experiments.
The local RBF collocation method has shown potential to solve partial differential equations over irregularly spaced data points accurately [20,21,27], but the field is usually smooth allowing for the shape parameter to be chosen to give high conditioning providing best accuracy. However, low conditioning allows for steep gradients and discontinuities to be captured using these RBF meshless techniques. Therefore, the RBF interpolation for high and low conditioning must be blended not only to capture the gradients but to keep the solution stable. The local blended RBF scheme is executed by forming two localized expansions using different shape parameter values, one for high conditioning and the other for low conditioning. We begin with the blending expression written as ) ( (21) 𝑓 (𝑥, 𝑦) = 𝑓𝐿 + 𝜙 𝑓𝐻 − 𝑓𝐿
𝑓𝐻 (𝑥, 𝑦) =
{
(31)
Substituting and rearranging, we have a blended RBF interpolation expression as ) ( [ ]−1 [ ]−1 𝑓 (𝑥, 𝑦) = (1 − 𝜙)𝐶𝐿 (𝑥, 𝑦)𝑇 𝐶𝐿 {𝑓 } + 𝜙 𝐶𝐻 (𝑥, 𝑦)𝑇 𝐶𝐻 {𝑓 } . (32) 84
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
governing equations. We can use this approach for the case of RBF interpolation by replacing the method of using different stencils by using different shape parameter RBF and determine the smoothest interpolation. For example, the smoothness indicator is determined by computing the first, second, and third derivative for the RBF collocation for each shape parameter used. The derivatives are found by expressions { } 𝜕 𝑓𝐻 𝜕 C𝐻 𝑇 [ ]−1 = (41) 𝐶𝐻 {𝑓 }, 𝜕𝑥 𝜕𝑥 𝜕 2 𝑓𝐻 𝜕 𝑥2 𝜕 3 𝑓𝐻 𝜕 𝑥3
{ =
=
𝜕 𝑓𝐿 = 𝜕𝑥 𝜕 2 𝑓𝐿 𝜕 𝑥2 𝜕 3 𝑓𝐿 𝜕 𝑥3
𝜕 𝑥2 {
and
{
𝜕 𝐶𝐻 𝜕𝑥
𝜕 𝐶𝐿 𝜕𝑥
{ =
}𝑇
}𝑇
𝜕 2 𝐶𝐿
𝜕 3 𝐶𝐿
}𝑇
}𝑇
𝜕 𝑥3
[ ]−1 𝐶𝐻 {𝑓 },
𝜌𝑢 ⎛ ⎞ ⎛0⎞ ⎟ 𝐹 𝑝 = ⎜𝑃 ⎟ 𝐹𝑐 = ⎜ 𝜌𝑢2 ) ⎟ ⎜( ⎜ ⎟ ⎝ 𝜌𝑒𝑡 + 𝑃 𝑢⎠ ⎝0⎠ ⎛ ⎞ ⎛0⎞ 𝜌𝑢 ⎜ ⎟ ⎜ ⎟ 2 𝜌𝑢 ⎟ 𝐹 = ⎜𝑃 ⎟ 𝐹𝑐 = ⎜ ⎜ 𝜌𝑢𝑣 ⎟ 𝑝 ⎜ 0 ⎟ ⎜(𝜌𝑒𝑡 + 𝑃 )𝑢⎟ ⎜0⎟ ⎝ ⎠ ⎝ ⎠
𝛽𝐿 =
(
[ ]−1 𝐶𝐻 {𝑓 }
⎛ ⎞ ⎛0⎞ 𝜌𝑣 ⎜ ⎟ ⎜ ⎟ 𝜌𝑢𝑣 0 ⎜ ⎟ 𝐺𝑐 = 𝐺 = ⎜ ⎟. ⎜ ⎟ 𝑝 ⎜𝑃 ⎟ 𝜌𝑣2 ( ) ⎜ 𝜌𝑒 + 𝑃 𝑣⎟ ⎜0⎟ ⎝ 𝑡 ⎠ ⎝ ⎠
(43)
[ ]−1 𝐶𝐿 {𝑓 }, [ ]−1 𝐶𝐿 {𝑓 }.
Δ𝑥2
𝜕 𝑓𝐿 𝜕𝑥
)2
( + Δ𝑥4
𝜕 2 𝑓𝐿 𝜕 𝑥2
)2
( + Δ𝑥6
𝜕 3 𝑓𝐿 𝜕 𝑥3
𝛾𝐿 𝛽𝐿 + 𝜀
where 𝜌 is density, u is the x-component of velocity, v is the y component of velocity, P is pressure and et is total energy. The convective flux is expressed in terms of Mach number
(45)
𝜌 ⎞ ⎛ ⎛0⎞ ⎟ 𝐹 = ⎜𝑃 ⎟ 𝐹𝑐 = 𝑀𝑎⎜ 𝜌𝑢 𝑝 )⎟ ⎜( ⎜ ⎟ ⎝ 𝜌𝑒𝑡 + 𝑃 ⎠ ⎝0⎠
)2 (48)
⎛ ⎞ ⎛0⎞ 𝜌 ⎜ ⎟ ⎜ ⎟ 𝜌𝑢 ⎟ 𝐹 = ⎜𝑃 ⎟ 𝐹𝑐 = 𝑀𝑎⎜ 𝑝 ⎜ 𝜌𝑣 ⎟ ⎜0⎟ ⎜(𝜌𝑒𝑡 + 𝑃 )⎟ ⎜0⎟ ⎝ ⎠ ⎝ ⎠
(59)
⎛ ⎞ ⎛0⎞ 𝜌 ⎜ ⎟ ⎜ ⎟ 𝜌𝑢 ⎟ 𝐺 = ⎜0⎟ 𝐺𝑐 = 𝑀𝑎⎜ 𝑝 ⎜ 𝜌𝑣 ⎟ ⎜𝑃 ⎟ ⎜(𝜌𝑒𝑡 + 𝑃 )⎟ ⎜0⎟ ⎝ ⎠ ⎝ ⎠
(60)
We can then write the flux at the face 𝑖 + 𝑓𝑖 + 1 = 𝑓𝑐 𝑖 + 1 + 𝑷 𝑖 + 1 2
(49)
2
2
𝑷 𝑖+ 1 2
2
(62)
2
⎛ 0 ⎞ ⎜ ⎟ = ⎜𝑝𝑖+ 1 .⎟ 2 ⎜ 0 ⎟ ⎝ ⎠
(63)
The polynomials ± and ℘± are used to split the fluxes into left and right running waves for upwinding and are calculated by ⎧ 1 (𝑀 ± |𝑀 |), ⎪2 =⎨ ⎪± 1 (𝑀 ± 1)2 , ⎩ 4
± 𝑖
𝑖𝑓 |𝑀 | > 1 (64) 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
and
(53)
⎧ 1 (𝑀 ± |𝑀 |)∕𝑀, ⎪2 =⎨ ⎪ 1 (𝑀 ± 1)2 (2 ∓ 𝑀 ), ⎩4
2
℘± 𝑖
5. Inviscid flux discretization The advection upstream splitting method or AUSM scheme is used to determine the inviscid fluxes that appear in the inviscid Euler equations and the Navier-Stokes equations. The AUSM scheme begins by noticing the convection, and acoustic waves are physically separate processes, and therefore the flux can be written by splitting the pressure from the inviscid flux [60,61]. 𝐹 = 𝐹𝑐 + 𝐹𝑝
(61)
and
and the flux is then calculated by 2
as
2
𝑓 𝐶 𝑖 + 1 = 𝑚 𝑖 + 1 Φ𝑖 + 1
(52)
𝑓𝑖+ 1 = 𝜙𝐿 𝑓𝐿 𝑖+ 1 + 𝜙𝐻 𝑓𝐻 𝑖+ 1 .
2
1 2
where
(50)
𝜔𝐿 𝜔𝐻 + 𝜔𝐿
(58)
or
where 𝛾 H and 𝛾 L are linear weights. The linear weights are chosen so that the high shape parameter interpolation is weighted higher than the low shape parameter approximation, i.e. 𝛾 H = 0.9 and 𝛾 L = 0.1. The parameter 𝜀 is chosen to be a small number to avoid zero in the denominator, such as 𝜀 = 1e−8 . Finally, the weights can be determined by 𝜔𝐻 𝜙𝐻 = (51) 𝜔𝐻 + 𝜔𝐿 𝜙𝐿 =
(57)
(44)
(46)
The weighting factors can be calculated by 𝛾𝐻 𝜔𝐻 = 𝛽𝐻 + 𝜀 𝜔𝐿 =
(56)
and
Many derivatives can be approximated and used in the smoothness formulation because of the infinitely smooth property of the multiquadric RBF. With the derivatives known, we can attempt to approximate the smoothness indicator for each RBF approximation. √ ( ) ( 2 )2 ( 3 )2 𝜕 𝑓𝐻 2 𝜕 𝑓𝐻 𝜕 𝑓𝐻 𝛽𝐻 = Δ𝑥2 + Δ𝑥4 + Δ𝑥6 (47) 𝜕𝑥 𝜕 𝑥2 𝜕 𝑥3 √
(55)
In two-dimensions, the inviscid flux are
(42)
[ ]−1 𝐶𝐿 {𝑓 },
𝜕 𝑥2 {
=
𝜕 2 𝐶𝐻
}𝑇
where
𝑖𝑓 |𝑀 | > 1
.
(65)
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
The Mach number and pressure at the half nodes are then determined by ( ) ± − ± 𝑚𝑖+ 1 = + (66) 𝑖 + 𝑖+1 , 𝑖 = 𝑀𝑖 2
and ( ) ± − ± 𝑝𝑖 + 1 = ℘+ 𝑗 𝑝𝑗 + ℘𝑗+1 𝑝𝑗+1 , ℘𝑖 = ℘ 𝑀𝑖 .
(54)
2
85
(67)
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Finally, applying upwinding to the fluxes Φ𝑖+ 1 by the criteria {
Φ𝑖 + 1 = 2
used in the topology to give a least squares approximation [55]. To determine the expansion coefficients, 𝛼 j , a least squares minimization is performed on the data set of points N described the expression ] [𝑁 𝑁𝑃 𝑁 ∑ ∑ ( ) ( ) ∑ ( ) ( ) 𝛼𝑗 𝑃𝑖 𝑥𝑘 𝑃𝑗 𝑥𝑘 = 𝑃𝑖 𝑥𝑘 𝑓 𝑥𝑘 . (79)
2
Φ𝑖 ,
𝑖𝑓 𝑚𝑗+ 1 > 1
Φ𝑖+1 ,
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
(68)
2
where
𝑗=1
𝜌𝑎 ⎛ ⎞ Φ = ⎜ 𝜌𝑢𝑎 ⎟ ⎜ ⎟ ⎝𝜌(𝑒𝑡 + 𝑝)𝑎⎠
(69)
2
(70)
𝐶𝑖,𝑗 =
2
The AUSM+ scheme is an update to the AUSM scheme to fix issues near stagnation points. The AUSM scheme discussed in the previous section is a special case of the AUSM+ scheme by allowing the speed of sound to be upwinded with the fluxes as presented before [61]. The AUSM+ scheme works by defining a common speed of sound denoted by 𝑎𝑖+ 1 .
2
𝑓 (𝑥 ) =
(71)
(72)
𝑎∗ . 𝑢𝑗
and 𝜕 𝑓𝑐 = 𝜕𝑥
2
where 𝑎∗ 2 𝑎̃ = . max(𝑎∗ , |𝑢|)
(75)
𝜕 𝑃𝑐 𝜕𝑥
}𝑇 1 ×𝑁𝑃
[𝐶 ]−1 𝑁 𝑃 ×𝑁 𝑃 [𝐶 ]𝑁𝑃 ×𝑁 {𝑓 }𝑁× 1 .
(87)
𝑖𝑓 |𝑀 | > 1 (76)
7. Compressible Euler equations
(77)
Many numerical experiments have been performed using the blended RBF interpolation scheme and are presented next. The code for all the numerical experiments presented are developed in Matlab R2018a. First, the solution of the 1-D inviscid Burgers equation is computed using the blended RBF approach and RBF enhanced finitedifference. This numerical experiment demonstrates the concept of blended RBF interpolation applied to the solution of a PDE with a shock or discontinuity. The governing equation and the initial condition are
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
⎧ 1 (𝑀 ± |𝑀 |)∕𝑀, ⎪2 =⎨ ( ) ⎪ 1 (𝑀 ± 1)2 (2 ∓ 𝑀 ) ± 𝛼𝑀 𝑀 2 − 1 2 , ⎩4
𝑖𝑓 |𝑀 | > 1
.
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
where the recommended values for 𝛼 and 𝛽 are
3 16
and 18 , respectively
𝜕𝑢 𝜕𝑢 +𝑢 =0 𝜕𝑡 𝜕𝑥 𝑢(𝑥, 0) = 1 𝑓 𝑜𝑟 𝑥 < 0.5 𝑢(𝑥, 0) = 0 𝑓 𝑜𝑟 𝑥 ≥ 0.5
6. Viscous flux discretization The localized RBF blended scheme can still produce oscillations in some scenarios. If instabilities occur, the scheme can switch to a moving least square (MLS) scheme to approximate the derivatives [55,62]. The MLS formulation is described in this section. The field variable f(x) can be approximated by the expansion 𝛼𝑗 𝑃𝑗
{
(86)
where the operator {𝐿𝑙𝑠 }𝑇1,𝑁 = {𝐿𝑃𝑐 }𝑇1 ×𝑁𝑃 [𝐶 ]−1 [𝐶 ]𝑁𝑃 ×𝑁 can 𝑁 𝑃 ×𝑁 𝑃 be precomputed and stored in a preprocessing stage similar to RBF interpolation.
and
𝑗=1
(85)
We can rewrite the expression for any linear differential operator { }𝑇 𝐿𝑓 (𝑥) = 𝐿𝑙𝑠 1 ×𝑁 {𝑓 }𝑁× 1 (88)
The equations used to blend the fluxes and pressure based on Mach number are
𝑁𝑃 ∑
(84)
𝑁𝑃
∑ 𝜕 𝑃𝑗 𝜕𝑓 = 𝛼 𝜕𝑥 𝑗=1 𝑗 𝜕𝑥
(73)
To account for the flow is subsonic, Eq. (74) calculates the speed of sound at the half nodes by ( ) (74) 𝑎𝑖+ 1 = min 𝑎̃𝐿 , 𝑎̃𝑅
𝑓 (𝑥 ) =
(83)
The derivatives are approximated by differentiation of the polynomials as
2
℘± 𝑖
𝛼𝑗 𝑃𝑗
{ }𝑇 𝑓𝑐 = 𝑃𝑐 1 ×𝑁𝑃 [𝐶 ]−1 𝑁 𝑃 ×𝑁 𝑃 [𝐶 ]𝑁𝑃 ×𝑁 {𝑓 }𝑁× 1 .
and
⎧ 1 (𝑀 ± |𝑀 |), ⎪2 ± = ⎨ 𝑖 ( ) ⎪± 1 (𝑀 ± 1)2 ± 𝛽 𝑀 2 − 1 2 , ⎩ 2
(82)
Finally, the field variable f(x) can be approximated by the expression
2
2
𝑁𝑃 ∑
and writing in matrix-vector notation, we have the expression { }𝑇 𝑓 (𝑥) = 𝑃𝑐 {𝛼}.
2
(𝛾 + 1)𝑎∗ 𝑎2 1 ℎ𝑡 = + 𝑢2 = 𝛾 −1 2 2(𝛾 − 1)
(81)
Recalling the approximation for the field variable f(x) is
The critical speed of sound then calculates the common speed of sound a∗ which can be expressed in terms on the total enthalpy ht , specific heat ratio 𝛾, and u. For an ideal gas, we have
𝑎𝑖+ 1 =
( ) ( ) 𝑃𝑖 𝑥𝑘 𝑃𝑗 𝑥𝑘 .
𝑗=1
= 𝑚 𝑖 + 1 𝑎 𝑖 + 1 Φ𝑖 + 1 2
𝑘=1
{𝛼}𝑁𝑃 × 1 = [𝐶 ]−1 𝑁 𝑃 ×𝑁 𝑃 [𝑃 ]𝑁𝑃 ×𝑁 {𝑓 }𝑁× 1 .
2
𝑖+ 1 2
𝑁 ∑
The expansion coefficients can be solved by
The flux is determined by the expression 𝑓𝑐
(80)
where
𝑓 𝑖 + 1 = 𝑚 𝑖 + 1 Φ𝑖 + 1 + 𝑷 𝑖 + 1 . 2
𝑘=1
[𝐶 ]𝑁 𝑃 ×𝑁 𝑃 {𝛼}𝑁𝑃 × 1 = [𝑃 ]𝑁𝑃 ×𝑁 {𝑓 }𝑁× 1
Therefore, the flux at the interface is 2
𝑘=1
Writing in matrix-vector form, we have
(89)
The collocation points are randomly distributed across the domain of [0,1] over x as shown in Fig. 3. The local collocation points are specified to be the points i + 1, i, and i − 1 where i is the data center, but the number of collocation points can be increased for better approximations. These local collocation points are used for the RBF interpolation along the subdomain [xi − 1 , xi + 1 ]. A distance Δx is used to interpolate for the value of un at xi − Δx to create an upwind stencil. With the upstream values of un known, a first order explicit scheme
(78)
where 𝛼 j are the expansion coefficients, NP are the number of polynomials, Pj . NP must be less than the number of influence points NF 86
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Table 2 Initial conditions for top hat profile. u(x, 0) = 0 u(x, 0) = 1 u(x, 0) = 0
x < 0.2 0.2 < x < 0.5 x < 0.5
blending scheme, the second derivative was monitored and if the curvature of the function was found to be too high based some criterion, the scheme can switch to low RBF interpolation. This strategy works well in some cases, but a robust blending criterion based on the second derivative of the field proved to be difficult to determine. A better approach is adopted to use the successive gradients and flux limiter. This strategy gave a robust criterion for blending high and low shape parameter RBF. The numerical experiments showed if low shape parameter interpolation is used excessively, the solution becomes dissipative. In the next example, we investigate reducing the amount of dissipation introduced into the solution using blended RBF interpolation. We seek to benefit from the high accuracy of RBF interpolation and form a high-resolution scheme for shock capturing. We observe that the low shape parameter becomes dissipative as the solution progresses. For high shape parameter, the solution does well at keeping the gradients steep with no dissipation, but this strategy tends to cause overshoots and undershoots near steep gradients. This is similar behavior found in mesh based schemes where low order schemes are dissipative and high order schemes are dispersive. Methods like total variation diminishing (TVD) and weighted essentially non-oscillatory schemes seek to minimize or eliminate oscillations by switching schemes near steep gradient and shocks [63]. The blended RBF scheme uses the attributes of the low and high shape parameter. We want to capture the steep gradients as accurately as possible without oscillations. The numerical experiment uses the 1-D advection equation given as
Fig. 3. Initial Condition for the non-linear Burgers’ equation solution.
in time, t and first order finite difference in x is used at some constant distance Δx to calculate the next value of u at the next time step, un + 1 . A high shape parameter value, rendering the RBF flat, is used to demonstrate the instabilities produced by high condition number interpolation near shocks and discontinuities. Condition numbers range from K = 2e7 − 2e10 as the scheme proceeds through the algorithm. A RBF reproduction is evaluated for each data center and within the collocation point of the local domain. As the solution marches in along space, the high shape parameter RBF reproduction causes overshoots near steep gradients. This behavior will eventually lead to instabilities like other high order methods as the schemes marches in time. This is shown in Fig. 4. Next, the local blended RBF collocation method is used, allowing the scheme to switch between high or low shape parameter RBF interpolation vectors. Fig. 5 shows that the condition number of the collocation matrix is high near smooth regions. As the scheme marches towards the discontinuity, the scheme switches to a low shape parameter RBF and low condition number. For Hardy multiquadric RBF, a low shape parameter causes the interpolation to be piecewise linear, providing stability near steep gradients and shocks. In Fig. 5, the condition number varies between K = 6.4 − 6e6 depending on the smoothness of the function. Early in the development of the
𝜕𝑢 𝜕𝑢 +𝑎 = 0. (90) 𝜕𝑡 𝜕𝑥 For the first example, a top-hat profile is used for the initial condition shown in Table 2. The solution is evaluated for three different cases of low, high, and a blended shape parameter as shown in Fig. 6. The minmod flux limiter is used for the blended solution. The shape parameter values used are c = 0.5 and c = 0.001 for the high and low interpolation. The solution of the one-dimensional advection equation using the blended RBF collocation method works well capturing the gradient as close to the high shape parameter solution while using the low shape parameter solution when oscillations are detected. Also, the results
Fig. 4. Solution sweep with a constant shape parameter interpolation approach. Condition numbers are a) At top of discontinuity: 2.5 × 107 b) At bottom of discontinuity: 2.2 × 1010 . 87
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Fig. 5. Solution sweep with blended interpolation approach. Condition numbers a) At the top of discontinuity: 6.4 b) At the bottom of discontinuity: 11.
Fig. 6. 1-D Advection using Example 1 initial conditions.
88
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Fig. 7. 1-D Advection equation using Example 2 initial conditions.
Table 3 Example 2 initial conditions. u(x, u(x, u(x, u(x,
0) = 0 0) = 2 0) = 1 0) = 0
x < 0.1 0.1 ≤ x ≤ 0.3 0.3 < x ≤ 0.5 x > 0.5
provide proof that the scheme is a high-order scheme for hyperbolic partial differential equations. To demonstrate the robustness of the blending approach a different set of initial conditions are used giving the u(x) profile given in Table 3 and the results are shown in Fig. 7. Applying this methodology to solve the inviscid compressible flow governing equations yields stable solutions comparing well to theoretical values. The governing equations for the inviscid Euler equations are 𝜕𝑄 𝜕𝐹 𝜕𝐺 + + =0 𝜕𝑡 𝜕𝑥 𝜕𝑦
Fig. 8. Compressible Euler equations geometry and boundary conditions.
The inlet boundary conditions are Mach number of 2, pressure is 100,000 Pa, density is 1.16 kg/m3 and the temperature is 300 K. The walls of the domain is set to a slip boundary condition forcing the velocity to be tangential to the walls, and the outflow boundary for u-velocity, v-velocity, density, pressure, and temperature are extrapolated from the interior of the domain. The problem setup is shown in Fig. 8. The solution using the blended RBF interpolation with AUSM+ upwinding is shown in Fig. 9 and the numerical comparisons to theoretical oblique shock relations is shown in Table 4. The
(91)
where ⎧ ⎫ ⎧ ⎫ ⎧ ⎫ 𝜌𝑢 𝜌𝑣 ⎪𝜌⎪ ⎪ ⎪ ⎪ ⎪ 2 ⎪ 𝜌𝑢 ⎪ ⎪ ⎪ ⎪ ⎪ 𝜌𝑢𝑣 𝑄 = ⎨ ⎬ 𝐹 = ⎨ 𝜌𝑢 + 𝑃 ⎬ 𝐺 = ⎨ ⎬ 2 ⎪ 𝜌𝑣 ⎪ ⎪( 𝜌𝑢𝑣 ) ⎪ ⎪( 𝜌𝑣 + 𝑃) ⎪ ⎪𝜌𝑒𝑡 ⎪ ⎪ 𝜌𝑒𝑡 + 𝑃 𝑢⎪ ⎪ 𝜌𝑒𝑡 + 𝑃 𝑣⎪ ⎩ ⎭ ⎩ ⎭ ⎩ ⎭ 89
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Fig. 10. Locations for analytical and numerical comparisons.
Fig. 9. Inviscid compressible flow using RBF blended interpolation.
Table 4 Comparison of results for localized blended RBF collocation. Oblique shock relations Location
Mach
Pressure (Pa)
Density (kg/m3 )
Temperature (K)
1 2 3
2.000 1.641 1.285
1.00E+05 1.71E+05 2.80E+05
1.161 1.693 2.405
300.000 351.045 405.958
Fig. 11. Inviscid compressible flow over supersonic airfoil.
Localized blended RBF collocation Location
Mach
Pressure (Pa)
Density (kg/m3 )
Temperature (K)
1 2 3
2 1.641 1.286
1.00E+05 1.71E+05 2.80E+05
1.161 1.694 2.406
300 351.05 405.9
oblique shock and expansion along the surface of the airfoil giving good agreement with theoretical shock relations.
8. Viscous compressible flow results
Finite volume method (FVM) Location
Mach
Pressure (Pa)
Density (kg/m3 )
Temperature (K)
1 2 3
2 1.63953 1.2887
1.00E+05 1.70E+05 2.79E+05
1.16198 1.68475 2.40177
300 350.544 405.132
The blended RBF interpolation scheme gave satisfactory results for the Euler equations, and we would like to add more physics to test the scheme further. The next step is to add viscosity to the fluid flow; therefore, we would like to solve the Navier-Stokes equations. Due to the nature of high Reynolds number flows, the boundary is thin, and it is stated that solving the Navier-Stokes for high-speed flow is one of the most difficult problems in computational fluid dynamics. The governing equations can be written in vector notation as
locations within the domain where comparisons are taken are shown in Fig. 10. Other geometries were used to investigate the robustness of the blended RBF interpolation scheme. Fig. 11 shows the solution for a symmetrical supersonic airfoil. The scheme is capable of capturing the
𝜕 𝐹𝑣 𝜕 𝐺𝑣 𝜕𝑄 𝜕𝐹 𝜕𝐺 + + = + 𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 90
(92)
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
where ⎧ ⎫ ⎧ ⎫ ⎧ ⎫ 𝜌𝑢 𝜌𝑣 ⎪𝜌⎪ ⎪ ⎪ ⎪ ⎪ 2 ⎪ 𝜌𝑢 ⎪ ⎪ ⎪ ⎪ ⎪ 𝜌𝑢𝑣 𝑄 = ⎨ ⎬ 𝐹 = ⎨ 𝜌𝑢 + 𝑃 ⎬ 𝐺 = ⎨ ⎬ 2 ⎪ 𝜌𝑣 ⎪ ⎪( 𝜌𝑢𝑣 ) ⎪ ⎪( 𝜌𝑣 + 𝑃) ⎪ ⎪𝜌𝑒𝑡 ⎪ ⎪ 𝜌𝑒𝑡 + 𝑃 𝑢⎪ ⎪ 𝜌𝑒𝑡 + 𝑃 𝑣⎪ ⎩ ⎭ ⎩ ⎭ ⎩ ⎭ and 𝑒𝑡 = 𝑒 +
) 1( 2 𝑢 + 𝑣2 2
𝑃 = 𝜌(𝛾 − 1)𝑒.
The viscous fluxes can be written in vector form as ⎧ ⎫ ⎧ ⎫ 0 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 𝜏𝑥𝑥 𝜏𝑥𝑥 𝐹𝑣 = ⎨ 𝐺𝑣 = ⎨ ⎬ ⎬. 𝜏𝑥𝑦 𝜏𝑥𝑦 ⎪ ⎪ ⎪ ⎪ ⎪𝑢𝜏𝑥𝑥 + 𝑣𝜏𝑥𝑦 − 𝑞𝑥 ⎪ ⎪𝑢𝜏𝑥𝑦 + 𝑣𝜏𝑦𝑦 − 𝑞𝑦 ⎪ ⎩ ⎭ ⎩ ⎭ The shear stresses are given by Pletcher et al. [64] ( ) 2 𝜕𝑢 𝜕𝑣 τ𝑥𝑥 = 𝜇 2 − 3 𝜕𝑥 𝜕𝑦 τ𝑦𝑦 =
( ) 2 𝜕𝑣 𝜕𝑢 𝜇 2 − 3 𝜕𝑦 𝜕𝑥 (
τ𝑥𝑦 = 𝜇
𝜕𝑢 𝜕𝑣 + 𝜕𝑦 𝜕𝑥
Fig. 13. Domain for Navier-Stokes solution with CHT.
(93)
(94)
) (95)
𝑞𝑥 = −𝑘
𝜕𝑇 𝜕𝑥
(96)
𝑞𝑦 = −𝑘
𝜕𝑇 𝜕𝑦
(97)
where u and v are the velocity components, T is the temperature, q is the heat flux, 𝜇 is the dynamic viscosity, T is temperature, and k is the thermal conductivity. The inlet and outlet boundary conditions remain the same as the Euler equation examples and are provided in Fig. 12. The wall boundary condition is using the no-slip condition. Therefore, the velocity at the wall is zero. A conjugate heat transfer boundary condition is prescibed to the lower fluid domain wall, and the temperature is solved within the solid domain using RBF direct differentiation. The solid material chosen for this case is titanium. The boundary conditions for solid walls are prescribed to be adiabatic, and the lower wall of the solid region is at a constant temperature of 300 K. The solution domain is shown in Fig 13. The Navier-Stokes solution for pressure and temperature is shown in Fig. 14. Points are added near the boundary to capture the boundary layer and the approximate the shear stresses. To overcome the non-uniform
Fig. 14. Pressure (Pa) and temperature(K) results using blended RBF interpolation.
nodal arrangement, the local support of the RBF should be scaled based on the maximum distance of nodes within the local support [65–67]. This approach reduces the shape parameter dependency when using non-uniform node distribution. 9. Conclusions The localized blended radial basis function (RBF) collocation method with Moving Least Squares has shown to solve the Euler and Navier-Stokes equations by solving several numerical experiments. For smooth functions, the RBF interpolation using a high condition
Fig. 12. Boundary conditions for Navier-Stokes with conjugate heat transfer. 91
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
collocation matrix does well in reproducing a test function and determining the derivatives. As the shape parameter is increased the error for the function and derivatives decrease exponentially. For most problems, the high shape parameter RBF should be used frequently. However, increasing the slope of the function causes the high shape parameter RBF interpolation to produce oscillations. We seek to use the RBF interpolation for the solution of PDEs, therefore, the oscillations must be reduced or eliminated for a stable solution. Fortunately, the use of low shape parameter or low condition number RBF interpolation dampens oscillations providing a piecewise linear reproduction of the function. A robust high-resolution scheme is developed benefiting from the spectral accuracy capabilities of the RBF collocation method by allowing the numerical scheme to switch between high and low RBF interpolation,. The blended RBF scheme was able to solve many different PDE example problems, such as the inviscid Burgers equation, two-dimensional advection equation to demonstrate the performance of the blended RBF scheme. Studying the blended RBF scheme using these simple numerical test cases allowed for the experimentation of different blending parameter strategies and ultimately finding a criterion that eliminates oscillation while minimizing the numerical dissipation introduced into the solution. The 2-D inviscid compressible Euler equations and Navier-Stokes equations are solved using the blended RBF scheme. The results for all schemes showed good agreement with analytical shock relations and commercial CFD solutions. The localized blended RBF collocation method provides a stable method that can be used to solve hyperbolic PDEs and introduces a new tool for future research. The many advantages of using meshless methods can be used for hyperbolic PDEs. For example, the capability of refining the solution by adding nodes during the computation can now be used effectively.
[22] Hardy RL. Multiquadric equations of topography and other irregular surfaces. J Geophys Res 1971;76(8):1905–15. [23] Hardy RL. Research results in the application of multiquadric equations to surverying and mapping problems. Surv Mapp 1975;35:321–32. [24] Franke R. Scattered data interpolation: tests of some methods. Math Comput 1982;38(157):181–200. [25] Dyn N, Levin D, Rippa S. Numerical procedures for surface fitting of scattered data by radial basis funtions. SIAM J Sci Stat Computat 1986;7(2):639–59. [26] Karageorghis A, Chen CS, Smyrlos Y-S. A matrix decomposition RBF algorithm: approximation of functions and their derivatives. Appl Numer Math 2007;57:304–19. [27] Sarler B, Vertnik R. Local explicit radial basis funtion collocation method for diffusion problems. Computat Math 2005:1269–82. [28] Yao G, Sarler B, Chen CS. A comparison of three explicit local meshless methods using radial basis funtions. Eng Anal Bound Elem 2011;35:600–9. [29] Lee CK, Liu X, Fan SC. Local multiquadric approximation for solving boundary value problems. Comput Mech 2003;30(5–6):396–409. [30] Gerace S, Erhart K, Kassab A, Divo E. A model-integrated localized collocation meshless method for large scal three dimensional heat transfer problems. Eng Anal 2014;45:2–19. [31] Hon Y-C, Sarler B, Yun D-f. Local radial basis function collocation method for solving thermo-driven fluid-flow problems with free surface. Eng Anal Bound Elem 2015;57:2–8. [32] Kelly J, Divo E, Kassab AJ. Numerical solution of the two-phase incompressible Navier-Stokes equations using a GPU-accelerated meshless method engineering analysis with boundary elements. Eng Anal 2014;40C:36–49. [33] Sarler B, Tran-Cong T, Chen CS. Meshfree direct and indirect local radial basis function collocation formulations for transport phenomena. Bound Elem XVII 2005:417–28. [34] Erhart K, Kassab AJ, Divo E. An inverse localized meshless technique for the determination of non-linear heat generation rates in living tissues. Int J Heat Fluid Flow 2008;18(3):401–14. [35] Cheng AHD, Golberg MA, Kansa EJ, Zammito G. Exponential convergence and H-c multiquadric collocation method for partial differenctial equations. Numer Methods Partial Differ Equ 2003:571–94. [36] Madych WR. Miscellaneous error bounds for multiquadric and related interpolations. Comput Math Applic 1992;24(12):121–38. [37] Buhmann MD, Dinew S, Larsson E. A note on radial basis funtion interpolant limits. IMA J Numer Anal 2010;30:543–54. [38] Madych WR. Bounds on multivariate polynomials and exponential error estimates for multiquadric interpolation. J Approx Theory 1992;70:94–114. [39] Roque CMC, Ferreira AJM. Numerical experiments on optimal shape parameters for radial basis functions. Numer Methods Partial Differ Equ 2009:1–14. [40] Wright GB, Fornberg B. Scattered node compact finite difference-type formulas generated from radial basis funtions. J Comput Phys 2006;212:99–123. [41] Larsson E, Fornberg B. A numerical study of some radial basis function based solution methods for elliptic PDEs. Comput Math Appl 2003;46:891–902. [42] Wang JG, Liu GR. On the optimal shape parameters of radial basis functions used for 2-D meshless methods. Comput Methods Appl Mech Engrg 2002;191:2611–30. [43] Fasshauer GE, Zhang JG. On choosing "optimal" shape parametes for RBF approximations. Numer Algor 2007;45:345–68. [44] Schaback R. Error estimates and condition numbers for radial basis function interpolation. Adv Comput Math 1995;3(3):251–64. [45] Driscoll TA, Fornberg B. Interpolation in the limit of increasingly flat radial basis function. Comput Math Appl 2002;43(3–5):413–22. [46] Buhmann M, Dyn N. Spectral convergence of multiquadric interpolation. In: Proceeding of the Edinburgh mathematical society, Edinburgh; 1993. [47] Larsson E, Fornberg B. Theoretical and compurtational aspects of multivariate interpolation with increasingly flat radial basis funtions. Comput Math Appl 2005;49:103–30. [48] Carlson RE, Foley TA. The parameter R2 in multiquadric interpolation. Comput Math Applic 1991;21(9):29–42. [49] Carlson RE, Foley TA. Interpolation of track data with radial basis methods. Comput Math Applic 1992;24(12):27–34. [50] Jung J-H, Durante VR. An iterative adaptive multiquadric radial basis function method for the detection of local jump discontinuities. Appl Numer Methods 2008;59:1449–66. [51] Harris M, Kassab A, Divo E. Application of an RBF blending interpolation method to problems with shocks. Comput Assisted Methods Eng Sci 2015;22:229–41. [52] Harris M, Kassab E, Divo E. An RBF interpolation blending scheme for effective shock-capturing. Int J Comput Methods Exp Meas 2017;5(3):281–92. [53] Harris MF, Kassab AJ, Divo EA. Meshless method using a RBF interpolation scheme for inviscid compressible flows. In: Second thermal and fluids engineering conference, Las Vegas, NV; 2017. [54] Harris MF, Kassab AJ, Divo EA. Viscous compressible flow soluiton using RBF blended interpolation. In: 3rd thermal and fluids engineering conference (TFEC), Fort Lauderdale, FL; 2018. [55] Pepper DW, Kassab AJ, Divo EA. Introduction to finite element, boundary element, and meshless methods. ASME Press; 2014. [56] Laney CB. Computational gasdynamics. Cambridge University Press; 1998. [57] Chung TJ. Computational fluid dynamics. New York, NY: Cambridge University Press; 2002. [58] Moukalled F, Mangani L, Darwish M. The finite volume method in computational fluid dynamics: an advanced introduction with OpenFOAM and Matlab. Switzerland: Springer International Publishing; 2016.
References [1] Liu GR. Mesh free methods. Boca Raton, FL: CRC Press; 2003. [2] Kolodziej JA, Grabski JK. Many names of the Trefftz method. Eng Anal Bound Elem 2018;96:169–78. [3] Chen W, Hon YC. Numerical convergence of boundary knot method in the analysis of Helmholtz, modified Helmholtz, and convection-diffusion problems. Comput Methods Appl Mech Eng 2003;192:1859–75. [4] Fairweather G, Karageorghis A. The method of fundamental solutions for elliptic boundary value problems. Adv Comput Math 1998;9:69–95. [5] Golberg MA. The method of fundamental solutions for Poisson’s equation. Eng Anal Bound Elem 1995;16:205–13. [6] Grabski JK, Kolodziej JA. Laminar flow of a power-law fluid between corrugated plates. J Mech Mater Struct 2016;11(1):23–40. [7] Chen W, Wang FZ. A method of fundamental solutions without fictitious boundary. Eng Anal Bound Elem 2010;34(5):530–2. [8] Belytschko T, Lu YY, Gu L. Element -free Galerkin methods. Int J Numer Methods Eng 1994:229–56. [9] Atluri T, Zhu T. A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998:117–27. [10] Melenk JM, Babuska I. The partion of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng 1996:289–314. [11] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992:307–18. [12] Duarte CA, Oden JT. Hp clouds- an hp meshless method. Numer Method Partial Diff Eq 1996:673–705. [13] Liu WK, Jun S, Li S, Adee J, Belytschko T. Reproducing kernel particle methods for structural dynamics. Int J Numer Methods Eng 1995:1655–79. [14] Kansa EJ. Multiquadrics - a scattered data approximation scheme with applications to computational fluid dynamics - I. Surface approximations and partial derivative estimates. Comput Math Appl 1990;19(8/9):127–45. [15] Kansa EJ. Multiquadrics - A Scattered data approximation scheme with applications to compuational fluid dynamics-II. Comput Math Appic 1990;19(8/9):147–61. [16] Kansa EJ, Carlson RE. Improved accuracy of multiquadric interpolation using variable shape parameters. Comput Math Applic 1992;24(12):99–120. [17] Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multiquadric radial basis funtions: applications to elliptic partial differential equations. Comput Math Appl 2000;39:123–37. [18] Kansa EJ, Holoborodko P. On the ill-conditioned nature of C∞ RBF strong collocation. Eng Anal Bound Elem 2017;78:26–30. [19] Divo EA, Kassab AJ. Iterative domain decomposition meshless method modeling of incompressible flows and conjugate heat transfer. Eng Anal 2006;30(6):465–78. [20] Divo EA, Kassab AJ. An efficient localized RBF meshless method for fluid flow and conjugate heat transfer. ASME J Heat Transfer 2007;129:124–36. [21] Divo E, Kassab AJ. Localized meshless modeling of natural convective viscous flows. Numer Heat Transfer 2008;53:487–509. 92
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
[59] Liu X, Osher S, Chan T. Weighted essentially nonoscillatory schemes. J Comput Phys 1994;115:200–12. [60] Liou M-S, Steffen CJ Jr. A new flux splitting scheme. J Comput Phys 1993;107:23–39. [61] Liou M-S. A sequel to AUSM: AUSM+. J Comput Phys 1996;129(0256):346–82. [62] Lancaster P, Salkauskas K. Surfaces generated by moving least squares methods. Math Comput 1981;37(155):141–58. [63] Harten A. A high resolution scheme for the computation of weak solutions of hyperbolic conservation laws. J Compuat Phys 1983;49(357).
[64] Pletcher RH, Tannehill JC, Anderson DA. Computational fluid mechanics and heat transfer. CRC Press Taylor & Francis Group, LLC; 2013. [65] Vertnik R, Sarler B. Solution of incompressible turbulent flow by a mesh-free method. Comput Model Eng Sci 2009;44(1):65–95. [66] u Islam S, Vertnik R, Sarler B. Local radial basis funtion collocation method along with explicit time stepping for hyperbolic partial differential equations. Appl Numer Math 2011;67:136–51. [67] Vertnik R, Sarler B. Simulation of continuous casting of steel by a meshless technique. Int J Cast Met Res 2009;22:311–13.
93
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
A shock-capturing meshless scheme using RBF blended interpolation and moving least squares Michael F. Harris a,b,∗, Alain J. Kassab a, Eduardo Divo b a b
Department of Mechanical and Aerospace Engineering, University of Central Florida, Orlando, FL, USA Department of Mechanical Engineering, Embry-Riddle Aeronautical University, Daytona Beach, FL, USA
a r t i c l e Keywords: Meshless methods Radial basis functions Moving least squares Shock-capturing
i n f o
a b s t r a c t Solving the full Navier-Stokes equations requires the use of numerical methods, especially for engineering application involving complex geometry. Numerical methods commonly used today all require a mesh before a solution attempt. Developing the mesh can require expensive preprocessing time, and the quality of the mesh can have major effects on the solution. Meshless methods have become a research interest due to the simplicity of using scattered data points. Radial basis functions (RBF) interpolation is a class of meshless method, and many of these methods have been developed to solve partial differential equations. Schemes using RBF interpolation are capable of multivariate interpolation over scattered data even for data with discontinuities. A meshless method scheme is presented capable of solving the compressible Navier-Stokes equations. A blended RBF interpolation method switching between low and high shape parameter interpolation is used to approximate the inviscid fluxes, while the viscous fluxes are approximated using moving least squares (MLS).
1. Introduction Solutions to nonlinear partial differential equations over complex geometry in most engineering applications require the use of numerical methods. The most common numerical methods used today include finite difference method (FDM), finite volume method, and finite element method (FEM), but these methods require a preprocessing step of developing a mesh. Research in numerical methods has focused on reducing the computational time for the solution while mesh generation remained a bottleneck. Mesh development can be time-consuming for general problems and even more problematic for compressible flow since shock locations are unknown prior to the solution. Mesh adapting methods have been developed to overcome this issue requiring an iterative process of solving and re-meshing capturing the sharp gradients. Generating a high-quality mesh while balancing computation resources while capturing the physics can be difficult. Meshless methods have been a research interest for many years and shown to have advantages over other numerical methods. Many types of meshless methods developed over the years can be used to solve the governing equations found in fluid dynamics, heat transfer, and structural mechanics [1]. Methods such as the Trefftz method [2], the boundary knot method [3], the method of fundamental solutions [4], the method of fundamental solutions in combination with radial basis functions [5,6], and the singular boundary method [7]. Other
∗
methods have been investigated, such as element free Galerkin methods [8–11], Hp cloud method [12] and partition of unity [13]. Among these meshless methods are an important group related to applications of radial basis functions, such as the Kansa Method [14–18]. These numerical methods use a unifying concept of interpolation, but the type of interpolation is different for each method. Meshless methods began from spectral or pseudo-spectral methods based on Legendre and Chebyshev polynomials requiring the point or nodal distribution to be uniform. The use of radial basis functions (RBF) allows the use of irregular spatial distribution of points or nodes [19–21]. Radial basis functions have been a subject of research for decades as a global and local collocation method in solving partial differential equations. The Hardy Multiquadrics RBF was developed for topographical mapping and has proven to be an impressive function for multivariate interpolation [22,23]. Franke determined that the Hardy Multiquadric RBF produces stable results and yields the most accurate results of all other methods studied [24]. Kansa developed the RBF collocation meshless method using Hardy Multiquadric to solve partial differential equations such as the linear advection-diffusion equation and the Poisson equation [14,15]. These numerical experiments demonstrated the efficiency of using RBF multiquadric in contrast to the finite difference method. In the Kansa method, the RBFs are applied globally in strong form to solve partial differential equations, but the method is difficult to manage when using large datasets due to
Corresponding author. E-mail addresses: [email protected] (M.F. Harris), [email protected] (A.J. Kassab), [email protected] (E. Divo).
https://doi.org/10.1016/j.enganabound.2019.08.019 Received 31 March 2019; Received in revised form 21 July 2019; Accepted 6 August 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Fig. 1. Hardy Multiquadrics for varying shape parameter, c.
Table 1 Radial basis functions. Name
Expression
Hardy Multiquadric (MQ) Inverse Hardy Multiquadric (IMQ)
𝜒(𝑟) = (𝑟 + 𝑐 ) 1 𝜒(𝑟) = 1
Gaussian
𝜒(𝑟) =
2
Carlson and Foley observed that the shape parameter value is dependent on the function and is independent of the number and location of data points [48,49]. Their work demonstrates Hardy Multiquadric interpolation on topography that lie along the same plane or track and the distance between the data points can vary by orders of magnitude. As the shape parameter approached zero, they observed the Hardy Multiquadrics interpolation yielding better accuracy in those regions. Carlson and Foley recommended the shape parameters to be reduced enough to remove any oscillations caused by the interpolation due to high gradients in measurements. The suggestion of using a low shape parameter and therefore piecewise linear interpolant may not be favorable strategy since accuracy is lost, but using steep RBF interpolation can provide stability in specific regions of high gradient. Jung and Durante also recommended the notion of reducing the shape parameter to zero for the Multiquadrics RBF near local discontinuities [50]. Their algorithm lowers the shape parameter to zero based on a developed smoothness criterion. The algorithm switches the shape parameter to zero for each collocation point if a jump discontinuity is detected. This iterative procedure proposed by the Jung and Durante works very well, eliminating oscillations in the data and capturing the discontinuities. The Hardy multiquadrics and inverse multiquadrics are the most impressive functions as these radial basis functions produce the best accuracy. These functions depend on a shape parameter that must be chosen properly for accurate approximations. The research in RBF interpolation has focused heavily on overcoming the ill-conditioning of RBF interpolation always seeking flat RBF interpolation for all types of problems. This approach works well for smooth functions, but tends to cause oscillations for data where steep gradients or shocks exist.
Shape Parameter 2
(𝑟2 +𝑐 2 ) 2 2 𝑒−𝑐 𝑟
1 2
c>0 c>0 c>0
ill-conditioning of the interpolation matrix. Domain decomposition [19], matrix preconditioning [25,26] and local RBF collocation methods [20,27–34] have been used to overcome the difficulties of large data sets and ill-conditioned matrices. RBF interpolation using Hardy Multiquadric RBF has shown to provide high accuracy and exponential convergence [35–38]. However, the Hardy multiquadrics and inverse multiquadrics are dependent on the shape parameter, c. These radial basis functions are given in Table 1. The optimal shape parameter for a given problem can be difficult to determine, and many researchers have focused their attention on finding a solution for the shape parameter problem [17,18,37,39–43]. For best accuracy, the shape parameter should be chosen to be a large value causing the condition number of the interpolation matrix to be close to ill-conditioned, rendering the RBF flat [35–37,44,45]. Spectral convergence is reached as the shape parameter c → ∞ [35,46,47]. As c → 0, the Multiquadric RBF becomes steep, and the accuracy is piecewise linear. The multiquadric RBF for different shape parameter values is shown in Fig. 1. Although, low shape parameter can be used to add stability, especially in areas of steep gradients and shocks, high shape parameter interpolation must be used for accuracy. 82
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Reducing the shape parameter to render the RBF steep and the interpolation matrix condition number low eliminates oscillations better, capturing the steep gradients. A new approach for solving hyperbolic PDEs is presented using a localized blended RBF collocation method [51–54]. The blended RBF strategy is to combine the benefits of flat RBF and steep RBF interpolation to form a high-resolution scheme for shock capturing. 2. Radial basis function interpolation Radial basis function interpolation approximates the solution of a function by the expansion 𝑓 (𝑥 ) =
𝑁 ∑ 𝑗=1
𝛼𝑗 𝜒𝑗 (𝑥)
(1)
where 𝛼 j are the expansion coefficients, 𝜒 j (x) is the expansion function, and N is the number of collocation points [14,22,55]. The coefficients 𝛼 j are determined through the collocation of the function f(x) at the points f(xi ) for i = 1…N. 𝑁 ∑ ( ) ( ) 𝑓 𝑥𝑖 = 𝛼𝑗 𝜒𝑗 𝑥𝑖
(2)
{𝑓 } = [C]{α}
(3)
Fig. 2. Local RBF collocation meshless method.
𝑗=1
The collocation matrix [C] is of the form ( ) ( ) … 𝜒𝑁 𝑥 1 ⎤ ⎡ 𝜒1 𝑥 1 ⎥. [C ] = ⎢ ⋮ ⋱ ⋮ ( )⎥ ⎢ ( ) … 𝜒𝑁 𝑥 𝑁 ⎦ ⎣ 𝜒1 𝑥 𝑁
decomposition has been used to solve the difficulties of global methods. Another strategy is to use a local RBF collocation meshless method [19–21]. The localized RBF collocation meshless method (LRCMM) will be described in this section [55]. The LRCMM uses a localized RBF expansion to interpolate the field variables and approximate the derivatives. A set of data points sampled around the data center (xi ,yi ) as shown in Fig. 2. The function f(x, y) can be approximated by the expansion as shown in Eq (10). The expansion can be written in matrix-vector form given by Eq. (11) where 𝜒 j (x,y) is chosen to be the Hardy Multiquadric RBF.
(4)
The expansion function, 𝜒 j (x), is chosen to be the Hardy Multiquadric radial basis function given by ( )1 𝜒𝑗 (𝑥 ) = 𝑟 2 + 𝑐 2 2
(5)
The coefficients, 𝛼 j , are then found by solving the linear system given by Eq. (6) and the function can be reconstructed by recalling Eq. (2). {α} = [C] {𝑓 } −1
𝑓 (𝑥, 𝑦) =
(6)
𝑗=1
The RBF can be augmented by a polynomial expansion given by 𝑓 (𝑥 ) =
𝑁 ∑ 𝑗=1
𝛼𝑗 𝜒𝑗 (𝑥) +
𝑁𝑃 ∑ 𝑗=1
𝛼𝑗 𝑃𝑗 (𝑥)
𝛼𝑗 χ𝑗 (𝑥, 𝑦) +
𝑁𝑃 ∑ 𝑗=1
𝛼𝑗+𝑁 𝑃𝑗 (𝑥)
𝑓 (𝑥, 𝑦) = {𝐶 (𝑥, 𝑦)}𝑇 {α} (7) 𝜒𝑗 (𝑥, 𝑦) =
where N is the number of collocation points and NP is the number of polynomials [55]. The polynomials are Pj (x) can be for example P = {1, x}. For two-dimensional problems, P = {1, x, y}. This addition of the polynomial expansion provides stability in constant and linear fields, but the added polynomial expansion can reduce accuracy in some cases [48,49]. The collocation matrix and the vector of supporting node field values are of the form in Eq. (8) and Eq. (9). ( ) ( ) ( ) ( ) … 𝜒𝑁 𝑥 1 𝑃 1 𝑥 1 … 𝑃𝑁𝑃 𝑥1 ⎤ ⎡ 𝜒1 𝑥 1 ⎢ ⎥ ⋮ ⋮ ⋱ ⋮ ⋱ ⋮ ( ) ( ) ( )⎥ ⎢ ( ) 𝜒1 𝑥 𝑁 … 𝜒𝑁 𝑥 𝑁 𝑃 1 𝑥 𝑁 … 𝑃𝑁𝑃 𝑥𝑁 ⎥ ⎢ ( ) ( ) [C ] = ⎢ 𝑃1 𝑥1 ⎥ 0 … 0 … 𝑃1 𝑥𝑁 ⎢ … ⎥ ⋱ … ⋮ ⋱ ⋮ ⎢ ⎥ ( ) ( ) ⎣𝑃𝑁𝑃 𝑥1 ⎦ 0 … 0 … 𝑃𝑁𝑃 𝑥𝑁 (𝑁 +𝑁 𝑃 )×(𝑁 +𝑁 𝑃 )
√ (
𝑥 − 𝑥𝑗
)2
(10)
(11)
( )2 + 𝑦 − 𝑦𝑗 + 𝑐 2
(12)
The shape parameter, c, must be chosen wisely to obtain the best accuracy. A linear system of equations are formed from the local support nodes near the data center (xi ,yi ). {𝑓 } = [C]{α}
(13)
{α} = [C]−1 {𝑓 }
(14)
Substituting Eq. (14) into Eq. (11), we can interpolate the field variable value anywhere within the domain by 𝑓 (𝑥, 𝑦) = {𝐶 (𝑥, 𝑦)}𝑇 [𝐶 ]−1 {𝑓 }.
(15)
The derivatives are approximated by applying the derivative operator over the local expansion as { }𝑇 𝜕𝑓 (𝑥, 𝑦) 𝜕𝐶 (𝑥, 𝑦) = (16) [𝐶 ]−1 {𝑓 } 𝜕𝑥 𝜕𝑥
(8) ( ) ⎧ 𝑓 𝑥1 ⎫ ⎪ ⎪ ⎪ ( ⋮ )⎪ ⎪𝑓 𝑥 𝑁 ⎪ {𝑓 } = ⎨ ⎬ ⎪ 0 ⎪ ⎪ ⋮ ⎪ ⎪ 0 ⎪ ⎩ ⎭(𝑁 +𝑁 𝑃 )× 1
𝑁 ∑
and 𝜕𝑓 (𝑥, 𝑦) = 𝜕𝑦
(9)
{
𝜕𝐶 (𝑥, 𝑦) 𝜕𝑦
}𝑇
[𝐶 ]−1 {𝑓 }.
(17)
The vectors {C(x, y)}T [C]−1 {f}, { 𝜕C(𝜕x𝑥,𝑦) }𝑇 [C]−1 and { 𝜕C(𝜕y𝑥,𝑦) }𝑇 [C]−1
can be pre-computed and saved for each data center prior to the computation. Therefore, the derivatives of field variables can be approximated by the product of { 𝜕C(𝜕x𝑥,𝑦) }𝑇 [C]−1 and the field variable values within
Early RBF collocation methods such as the Kansa method are global methods and can be difficult to use for large-scale problems. Domain 83
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
the subdomain. 𝑓 (𝑥, 𝑦) = {𝐶 (𝑥, 𝑦)}𝑇 [𝐶 ]−1 {𝑓 } 𝜕𝑓 (𝑥, 𝑦) = 𝜕𝑥 𝜕𝑓 (𝑥, 𝑦) = 𝜕𝑦
{
{
𝜕𝐶 (𝑥, 𝑦) 𝜕𝑥 𝜕𝐶 (𝑥, 𝑦) 𝜕𝑦
}𝑇 }𝑇
Similarly, the derivatives are approximated by differentiation of the RBF vector { { } } 𝜕 𝐶𝐿 (𝑥, 𝑦) 𝑇 [ ]−1 𝜕 𝐶𝐻 (𝑥, 𝑦) 𝑇 [ ]−1 𝜕𝑓 (𝑥, 𝑦) 𝐶𝐿 {𝑓 } +𝜙 𝐶𝐻 {𝑓 } = (1 − 𝜙) 𝜕𝑥 𝜕𝑥 𝜕𝑥 (33)
(18)
[𝐶 ]−1 {𝑓 }
(19)
[𝐶 ]−1 {𝑓 }
(20)
and 𝜕𝑓 (𝑥, 𝑦) = (1 − 𝜙) 𝜕𝑦
3. Blended RBF collocation method
𝑓𝐿 (𝑥, 𝑦) =
𝑗=1 𝑁 ∑ 𝑗=1
𝛼𝐻 𝑗 𝜒𝐻 𝑗 (𝑥, 𝑦) +
𝛼𝐿 𝑗 𝜒𝐿 𝑗 (𝑥, 𝑦) +
𝑁𝑃 ∑ 𝑗=1
𝑁𝑃 ∑ 𝑗=1
𝛼𝑗+𝑁 𝑃𝑗 (𝑥)
𝛼𝑗+𝑁 𝑃𝑗 (𝑥)
(22)
(23)
2
(
)
𝑓𝐻 𝑥𝑖 , 𝑦𝑖 = ( ) 𝑓𝐿 𝑥𝑖, 𝑦𝑖 =
𝑗=1 𝑁 ∑ 𝑗=1
(
𝛼𝐻 𝑗 𝐶𝐻 𝑗 𝑥𝑖 , 𝑦𝑖 (
𝛼𝐿 𝑗 𝐶𝐿 𝑗 𝑥𝑖 , 𝑦𝑖
)
)
2
𝑟𝑖−1 =
(29)
(36)
2
𝑢𝑖−1 − 𝑢𝑖−2 . 𝑢𝑖 − 𝑢𝑖−1
(38)
The blending parameter in regions where high gradients, shocks or discontinuities are present the minmod limiter reduces a value less between [0,1]. If the value of the successive gradient is less than zero the minmod limiter is equal to zero providing low shape parameter interpolation. Another approach taken that involves using the smoothness indicator developed for the weighted essentially non-oscillatory or WENO scheme [59]. The WENO scheme uses polynomials to interpolate the fluxes at the midpoints of the stencil. Many stencils are used in the formulation to provide high order accuracy, and the smoothness indicator determines the stencil that provides the smoothest flux approximation by comparing the L2 norm values of the derivatives of the polynomial used. This flux is then used to approximate the derivatives in the
Solving for the expansion coefficients by inverting the collocation matrices [CL ] and [CH ], we have the expressions { } [ ]−1 (30) α𝐻 = 𝐶𝐻 {𝑓 } and { } [ ]−1 α𝐿 = 𝐶𝐿 {𝑓 }.
2
The local RBF collocation method gives the versatility of interpolating the field variable for the successive gradients. The minmod limiter is used for the blending parameter and is calculated by ( ) ( ( )) 𝜙 𝑟𝑖 = max 0, min 𝑟𝑖 , 1 (39) ( ) ( ( )) 𝜙 𝑟𝑖−1 = max 0, min 𝑟𝑖−1 , 1 (40)
(27)
[ ]{ } 𝐶𝐿 α𝐿 = {𝑓 }
2
and
(26)
(28)
2
2
( ) ( ) + 𝜙 𝑟𝑖−1 𝑓𝐻 𝑖− 1 − 𝑓𝐿 𝑖− 1
where the blending parameter is a function of the successive gradient r. The successive gradient is evaluated by the expressions 𝑢 − 𝑢𝑖−1 𝑟𝑖 = 𝑖 (37) 𝑢𝑖+1 − 𝑢𝑖
(25)
[ ]{ } C𝐻 α𝐻 = {𝑓 }
2
𝑓𝑖 − 1 = 𝑓𝐿 𝑖 − 1
Next, the weights 𝛼 Hj and 𝛼 Lj must be determined. A system of equations can be formed using all the field variable values within the local topology. 𝑁 ∑
} { [ ]−1 𝜕 𝐶𝐻 (𝑥, 𝑦) 𝑇 [ ]−1 𝐶𝐿 {𝑓 } +𝜙 𝐶𝐻 {𝑓 }. 𝜕𝑦 (34)
A method of determining how to blend between high and low shape parameter is needed. Fortunately, high-resolution schemes exist in computational fluid dynamics providing insight since these schemes switch between low and high order discretization for shock capturing. In CFD, flux or slope limiters are used to limit the gradients in the solver to realistic values. These flux/slope limiters are functions of the successive gradients found within the solution domain. Many flux limiters have been developed over the years [56–58]. For this research, a simple minmod flux limiter will be used. The minmod flux limiter can vary between [0,1] depending on the successive gradient r. For example, the flux is calculated by ) ( ( ) 𝑓𝑖 + 1 = 𝑓𝐿 𝑖 + 1 + 𝜙 𝑟𝑖 𝑓𝐻 𝑖 + 1 − 𝑓𝐿 𝑖 + 1 (35)
Defining two RBF with high and low valued shape parameters cL and cH . The Hardy Multiquadrics RBF is used as an example √ ( )2 ( )2 χ𝐿 𝑗 (𝑥, 𝑦) = (24) 𝑥 − 𝑥𝑗 + 𝑦 − 𝑦𝑗 + 𝑐𝐿2 √ ( )2 ( )2 2 χ𝐻 𝑗 (𝑥, 𝑦) = 𝑥 − 𝑥𝑗 + 𝑦 − 𝑦𝑗 + 𝑐𝐻
}𝑇
4. Blending parameter for shock capturing
where fH may be the smooth function interpolation and fL the steep gradient interpolation. 𝑁 ∑
𝜕 𝐶𝐿 (𝑥, 𝑦) 𝜕𝑦
The term 𝜑 is the blending term where 𝜑 = 1 will provide high shape parameter RBF interpolation and 𝜑 = 0 will give low shape parameter interpolation. A method is needed to determine the appropriate way of determining the blending parameter 𝜑. Methods using flux limiters or smoothness indicator, such a total variation diminishing (TVD) and weighted essentially non-oscillatory (WENO) schemes are good candidates for numerical experiments.
The local RBF collocation method has shown potential to solve partial differential equations over irregularly spaced data points accurately [20,21,27], but the field is usually smooth allowing for the shape parameter to be chosen to give high conditioning providing best accuracy. However, low conditioning allows for steep gradients and discontinuities to be captured using these RBF meshless techniques. Therefore, the RBF interpolation for high and low conditioning must be blended not only to capture the gradients but to keep the solution stable. The local blended RBF scheme is executed by forming two localized expansions using different shape parameter values, one for high conditioning and the other for low conditioning. We begin with the blending expression written as ) ( (21) 𝑓 (𝑥, 𝑦) = 𝑓𝐿 + 𝜙 𝑓𝐻 − 𝑓𝐿
𝑓𝐻 (𝑥, 𝑦) =
{
(31)
Substituting and rearranging, we have a blended RBF interpolation expression as ) ( [ ]−1 [ ]−1 𝑓 (𝑥, 𝑦) = (1 − 𝜙)𝐶𝐿 (𝑥, 𝑦)𝑇 𝐶𝐿 {𝑓 } + 𝜙 𝐶𝐻 (𝑥, 𝑦)𝑇 𝐶𝐻 {𝑓 } . (32) 84
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
governing equations. We can use this approach for the case of RBF interpolation by replacing the method of using different stencils by using different shape parameter RBF and determine the smoothest interpolation. For example, the smoothness indicator is determined by computing the first, second, and third derivative for the RBF collocation for each shape parameter used. The derivatives are found by expressions { } 𝜕 𝑓𝐻 𝜕 C𝐻 𝑇 [ ]−1 = (41) 𝐶𝐻 {𝑓 }, 𝜕𝑥 𝜕𝑥 𝜕 2 𝑓𝐻 𝜕 𝑥2 𝜕 3 𝑓𝐻 𝜕 𝑥3
{ =
=
𝜕 𝑓𝐿 = 𝜕𝑥 𝜕 2 𝑓𝐿 𝜕 𝑥2 𝜕 3 𝑓𝐿 𝜕 𝑥3
𝜕 𝑥2 {
and
{
𝜕 𝐶𝐻 𝜕𝑥
𝜕 𝐶𝐿 𝜕𝑥
{ =
}𝑇
}𝑇
𝜕 2 𝐶𝐿
𝜕 3 𝐶𝐿
}𝑇
}𝑇
𝜕 𝑥3
[ ]−1 𝐶𝐻 {𝑓 },
𝜌𝑢 ⎛ ⎞ ⎛0⎞ ⎟ 𝐹 𝑝 = ⎜𝑃 ⎟ 𝐹𝑐 = ⎜ 𝜌𝑢2 ) ⎟ ⎜( ⎜ ⎟ ⎝ 𝜌𝑒𝑡 + 𝑃 𝑢⎠ ⎝0⎠ ⎛ ⎞ ⎛0⎞ 𝜌𝑢 ⎜ ⎟ ⎜ ⎟ 2 𝜌𝑢 ⎟ 𝐹 = ⎜𝑃 ⎟ 𝐹𝑐 = ⎜ ⎜ 𝜌𝑢𝑣 ⎟ 𝑝 ⎜ 0 ⎟ ⎜(𝜌𝑒𝑡 + 𝑃 )𝑢⎟ ⎜0⎟ ⎝ ⎠ ⎝ ⎠
𝛽𝐿 =
(
[ ]−1 𝐶𝐻 {𝑓 }
⎛ ⎞ ⎛0⎞ 𝜌𝑣 ⎜ ⎟ ⎜ ⎟ 𝜌𝑢𝑣 0 ⎜ ⎟ 𝐺𝑐 = 𝐺 = ⎜ ⎟. ⎜ ⎟ 𝑝 ⎜𝑃 ⎟ 𝜌𝑣2 ( ) ⎜ 𝜌𝑒 + 𝑃 𝑣⎟ ⎜0⎟ ⎝ 𝑡 ⎠ ⎝ ⎠
(43)
[ ]−1 𝐶𝐿 {𝑓 }, [ ]−1 𝐶𝐿 {𝑓 }.
Δ𝑥2
𝜕 𝑓𝐿 𝜕𝑥
)2
( + Δ𝑥4
𝜕 2 𝑓𝐿 𝜕 𝑥2
)2
( + Δ𝑥6
𝜕 3 𝑓𝐿 𝜕 𝑥3
𝛾𝐿 𝛽𝐿 + 𝜀
where 𝜌 is density, u is the x-component of velocity, v is the y component of velocity, P is pressure and et is total energy. The convective flux is expressed in terms of Mach number
(45)
𝜌 ⎞ ⎛ ⎛0⎞ ⎟ 𝐹 = ⎜𝑃 ⎟ 𝐹𝑐 = 𝑀𝑎⎜ 𝜌𝑢 𝑝 )⎟ ⎜( ⎜ ⎟ ⎝ 𝜌𝑒𝑡 + 𝑃 ⎠ ⎝0⎠
)2 (48)
⎛ ⎞ ⎛0⎞ 𝜌 ⎜ ⎟ ⎜ ⎟ 𝜌𝑢 ⎟ 𝐹 = ⎜𝑃 ⎟ 𝐹𝑐 = 𝑀𝑎⎜ 𝑝 ⎜ 𝜌𝑣 ⎟ ⎜0⎟ ⎜(𝜌𝑒𝑡 + 𝑃 )⎟ ⎜0⎟ ⎝ ⎠ ⎝ ⎠
(59)
⎛ ⎞ ⎛0⎞ 𝜌 ⎜ ⎟ ⎜ ⎟ 𝜌𝑢 ⎟ 𝐺 = ⎜0⎟ 𝐺𝑐 = 𝑀𝑎⎜ 𝑝 ⎜ 𝜌𝑣 ⎟ ⎜𝑃 ⎟ ⎜(𝜌𝑒𝑡 + 𝑃 )⎟ ⎜0⎟ ⎝ ⎠ ⎝ ⎠
(60)
We can then write the flux at the face 𝑖 + 𝑓𝑖 + 1 = 𝑓𝑐 𝑖 + 1 + 𝑷 𝑖 + 1 2
(49)
2
2
𝑷 𝑖+ 1 2
2
(62)
2
⎛ 0 ⎞ ⎜ ⎟ = ⎜𝑝𝑖+ 1 .⎟ 2 ⎜ 0 ⎟ ⎝ ⎠
(63)
The polynomials ± and ℘± are used to split the fluxes into left and right running waves for upwinding and are calculated by ⎧ 1 (𝑀 ± |𝑀 |), ⎪2 =⎨ ⎪± 1 (𝑀 ± 1)2 , ⎩ 4
± 𝑖
𝑖𝑓 |𝑀 | > 1 (64) 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
and
(53)
⎧ 1 (𝑀 ± |𝑀 |)∕𝑀, ⎪2 =⎨ ⎪ 1 (𝑀 ± 1)2 (2 ∓ 𝑀 ), ⎩4
2
℘± 𝑖
5. Inviscid flux discretization The advection upstream splitting method or AUSM scheme is used to determine the inviscid fluxes that appear in the inviscid Euler equations and the Navier-Stokes equations. The AUSM scheme begins by noticing the convection, and acoustic waves are physically separate processes, and therefore the flux can be written by splitting the pressure from the inviscid flux [60,61]. 𝐹 = 𝐹𝑐 + 𝐹𝑝
(61)
and
and the flux is then calculated by 2
as
2
𝑓 𝐶 𝑖 + 1 = 𝑚 𝑖 + 1 Φ𝑖 + 1
(52)
𝑓𝑖+ 1 = 𝜙𝐿 𝑓𝐿 𝑖+ 1 + 𝜙𝐻 𝑓𝐻 𝑖+ 1 .
2
1 2
where
(50)
𝜔𝐿 𝜔𝐻 + 𝜔𝐿
(58)
or
where 𝛾 H and 𝛾 L are linear weights. The linear weights are chosen so that the high shape parameter interpolation is weighted higher than the low shape parameter approximation, i.e. 𝛾 H = 0.9 and 𝛾 L = 0.1. The parameter 𝜀 is chosen to be a small number to avoid zero in the denominator, such as 𝜀 = 1e−8 . Finally, the weights can be determined by 𝜔𝐻 𝜙𝐻 = (51) 𝜔𝐻 + 𝜔𝐿 𝜙𝐿 =
(57)
(44)
(46)
The weighting factors can be calculated by 𝛾𝐻 𝜔𝐻 = 𝛽𝐻 + 𝜀 𝜔𝐿 =
(56)
and
Many derivatives can be approximated and used in the smoothness formulation because of the infinitely smooth property of the multiquadric RBF. With the derivatives known, we can attempt to approximate the smoothness indicator for each RBF approximation. √ ( ) ( 2 )2 ( 3 )2 𝜕 𝑓𝐻 2 𝜕 𝑓𝐻 𝜕 𝑓𝐻 𝛽𝐻 = Δ𝑥2 + Δ𝑥4 + Δ𝑥6 (47) 𝜕𝑥 𝜕 𝑥2 𝜕 𝑥3 √
(55)
In two-dimensions, the inviscid flux are
(42)
[ ]−1 𝐶𝐿 {𝑓 },
𝜕 𝑥2 {
=
𝜕 2 𝐶𝐻
}𝑇
where
𝑖𝑓 |𝑀 | > 1
.
(65)
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
The Mach number and pressure at the half nodes are then determined by ( ) ± − ± 𝑚𝑖+ 1 = + (66) 𝑖 + 𝑖+1 , 𝑖 = 𝑀𝑖 2
and ( ) ± − ± 𝑝𝑖 + 1 = ℘+ 𝑗 𝑝𝑗 + ℘𝑗+1 𝑝𝑗+1 , ℘𝑖 = ℘ 𝑀𝑖 .
(54)
2
85
(67)
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Finally, applying upwinding to the fluxes Φ𝑖+ 1 by the criteria {
Φ𝑖 + 1 = 2
used in the topology to give a least squares approximation [55]. To determine the expansion coefficients, 𝛼 j , a least squares minimization is performed on the data set of points N described the expression ] [𝑁 𝑁𝑃 𝑁 ∑ ∑ ( ) ( ) ∑ ( ) ( ) 𝛼𝑗 𝑃𝑖 𝑥𝑘 𝑃𝑗 𝑥𝑘 = 𝑃𝑖 𝑥𝑘 𝑓 𝑥𝑘 . (79)
2
Φ𝑖 ,
𝑖𝑓 𝑚𝑗+ 1 > 1
Φ𝑖+1 ,
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
(68)
2
where
𝑗=1
𝜌𝑎 ⎛ ⎞ Φ = ⎜ 𝜌𝑢𝑎 ⎟ ⎜ ⎟ ⎝𝜌(𝑒𝑡 + 𝑝)𝑎⎠
(69)
2
(70)
𝐶𝑖,𝑗 =
2
The AUSM+ scheme is an update to the AUSM scheme to fix issues near stagnation points. The AUSM scheme discussed in the previous section is a special case of the AUSM+ scheme by allowing the speed of sound to be upwinded with the fluxes as presented before [61]. The AUSM+ scheme works by defining a common speed of sound denoted by 𝑎𝑖+ 1 .
2
𝑓 (𝑥 ) =
(71)
(72)
𝑎∗ . 𝑢𝑗
and 𝜕 𝑓𝑐 = 𝜕𝑥
2
where 𝑎∗ 2 𝑎̃ = . max(𝑎∗ , |𝑢|)
(75)
𝜕 𝑃𝑐 𝜕𝑥
}𝑇 1 ×𝑁𝑃
[𝐶 ]−1 𝑁 𝑃 ×𝑁 𝑃 [𝐶 ]𝑁𝑃 ×𝑁 {𝑓 }𝑁× 1 .
(87)
𝑖𝑓 |𝑀 | > 1 (76)
7. Compressible Euler equations
(77)
Many numerical experiments have been performed using the blended RBF interpolation scheme and are presented next. The code for all the numerical experiments presented are developed in Matlab R2018a. First, the solution of the 1-D inviscid Burgers equation is computed using the blended RBF approach and RBF enhanced finitedifference. This numerical experiment demonstrates the concept of blended RBF interpolation applied to the solution of a PDE with a shock or discontinuity. The governing equation and the initial condition are
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
⎧ 1 (𝑀 ± |𝑀 |)∕𝑀, ⎪2 =⎨ ( ) ⎪ 1 (𝑀 ± 1)2 (2 ∓ 𝑀 ) ± 𝛼𝑀 𝑀 2 − 1 2 , ⎩4
𝑖𝑓 |𝑀 | > 1
.
𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
where the recommended values for 𝛼 and 𝛽 are
3 16
and 18 , respectively
𝜕𝑢 𝜕𝑢 +𝑢 =0 𝜕𝑡 𝜕𝑥 𝑢(𝑥, 0) = 1 𝑓 𝑜𝑟 𝑥 < 0.5 𝑢(𝑥, 0) = 0 𝑓 𝑜𝑟 𝑥 ≥ 0.5
6. Viscous flux discretization The localized RBF blended scheme can still produce oscillations in some scenarios. If instabilities occur, the scheme can switch to a moving least square (MLS) scheme to approximate the derivatives [55,62]. The MLS formulation is described in this section. The field variable f(x) can be approximated by the expansion 𝛼𝑗 𝑃𝑗
{
(86)
where the operator {𝐿𝑙𝑠 }𝑇1,𝑁 = {𝐿𝑃𝑐 }𝑇1 ×𝑁𝑃 [𝐶 ]−1 [𝐶 ]𝑁𝑃 ×𝑁 can 𝑁 𝑃 ×𝑁 𝑃 be precomputed and stored in a preprocessing stage similar to RBF interpolation.
and
𝑗=1
(85)
We can rewrite the expression for any linear differential operator { }𝑇 𝐿𝑓 (𝑥) = 𝐿𝑙𝑠 1 ×𝑁 {𝑓 }𝑁× 1 (88)
The equations used to blend the fluxes and pressure based on Mach number are
𝑁𝑃 ∑
(84)
𝑁𝑃
∑ 𝜕 𝑃𝑗 𝜕𝑓 = 𝛼 𝜕𝑥 𝑗=1 𝑗 𝜕𝑥
(73)
To account for the flow is subsonic, Eq. (74) calculates the speed of sound at the half nodes by ( ) (74) 𝑎𝑖+ 1 = min 𝑎̃𝐿 , 𝑎̃𝑅
𝑓 (𝑥 ) =
(83)
The derivatives are approximated by differentiation of the polynomials as
2
℘± 𝑖
𝛼𝑗 𝑃𝑗
{ }𝑇 𝑓𝑐 = 𝑃𝑐 1 ×𝑁𝑃 [𝐶 ]−1 𝑁 𝑃 ×𝑁 𝑃 [𝐶 ]𝑁𝑃 ×𝑁 {𝑓 }𝑁× 1 .
and
⎧ 1 (𝑀 ± |𝑀 |), ⎪2 ± = ⎨ 𝑖 ( ) ⎪± 1 (𝑀 ± 1)2 ± 𝛽 𝑀 2 − 1 2 , ⎩ 2
(82)
Finally, the field variable f(x) can be approximated by the expression
2
2
𝑁𝑃 ∑
and writing in matrix-vector notation, we have the expression { }𝑇 𝑓 (𝑥) = 𝑃𝑐 {𝛼}.
2
(𝛾 + 1)𝑎∗ 𝑎2 1 ℎ𝑡 = + 𝑢2 = 𝛾 −1 2 2(𝛾 − 1)
(81)
Recalling the approximation for the field variable f(x) is
The critical speed of sound then calculates the common speed of sound a∗ which can be expressed in terms on the total enthalpy ht , specific heat ratio 𝛾, and u. For an ideal gas, we have
𝑎𝑖+ 1 =
( ) ( ) 𝑃𝑖 𝑥𝑘 𝑃𝑗 𝑥𝑘 .
𝑗=1
= 𝑚 𝑖 + 1 𝑎 𝑖 + 1 Φ𝑖 + 1 2
𝑘=1
{𝛼}𝑁𝑃 × 1 = [𝐶 ]−1 𝑁 𝑃 ×𝑁 𝑃 [𝑃 ]𝑁𝑃 ×𝑁 {𝑓 }𝑁× 1 .
2
𝑖+ 1 2
𝑁 ∑
The expansion coefficients can be solved by
The flux is determined by the expression 𝑓𝑐
(80)
where
𝑓 𝑖 + 1 = 𝑚 𝑖 + 1 Φ𝑖 + 1 + 𝑷 𝑖 + 1 . 2
𝑘=1
[𝐶 ]𝑁 𝑃 ×𝑁 𝑃 {𝛼}𝑁𝑃 × 1 = [𝑃 ]𝑁𝑃 ×𝑁 {𝑓 }𝑁× 1
Therefore, the flux at the interface is 2
𝑘=1
Writing in matrix-vector form, we have
(89)
The collocation points are randomly distributed across the domain of [0,1] over x as shown in Fig. 3. The local collocation points are specified to be the points i + 1, i, and i − 1 where i is the data center, but the number of collocation points can be increased for better approximations. These local collocation points are used for the RBF interpolation along the subdomain [xi − 1 , xi + 1 ]. A distance Δx is used to interpolate for the value of un at xi − Δx to create an upwind stencil. With the upstream values of un known, a first order explicit scheme
(78)
where 𝛼 j are the expansion coefficients, NP are the number of polynomials, Pj . NP must be less than the number of influence points NF 86
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Table 2 Initial conditions for top hat profile. u(x, 0) = 0 u(x, 0) = 1 u(x, 0) = 0
x < 0.2 0.2 < x < 0.5 x < 0.5
blending scheme, the second derivative was monitored and if the curvature of the function was found to be too high based some criterion, the scheme can switch to low RBF interpolation. This strategy works well in some cases, but a robust blending criterion based on the second derivative of the field proved to be difficult to determine. A better approach is adopted to use the successive gradients and flux limiter. This strategy gave a robust criterion for blending high and low shape parameter RBF. The numerical experiments showed if low shape parameter interpolation is used excessively, the solution becomes dissipative. In the next example, we investigate reducing the amount of dissipation introduced into the solution using blended RBF interpolation. We seek to benefit from the high accuracy of RBF interpolation and form a high-resolution scheme for shock capturing. We observe that the low shape parameter becomes dissipative as the solution progresses. For high shape parameter, the solution does well at keeping the gradients steep with no dissipation, but this strategy tends to cause overshoots and undershoots near steep gradients. This is similar behavior found in mesh based schemes where low order schemes are dissipative and high order schemes are dispersive. Methods like total variation diminishing (TVD) and weighted essentially non-oscillatory schemes seek to minimize or eliminate oscillations by switching schemes near steep gradient and shocks [63]. The blended RBF scheme uses the attributes of the low and high shape parameter. We want to capture the steep gradients as accurately as possible without oscillations. The numerical experiment uses the 1-D advection equation given as
Fig. 3. Initial Condition for the non-linear Burgers’ equation solution.
in time, t and first order finite difference in x is used at some constant distance Δx to calculate the next value of u at the next time step, un + 1 . A high shape parameter value, rendering the RBF flat, is used to demonstrate the instabilities produced by high condition number interpolation near shocks and discontinuities. Condition numbers range from K = 2e7 − 2e10 as the scheme proceeds through the algorithm. A RBF reproduction is evaluated for each data center and within the collocation point of the local domain. As the solution marches in along space, the high shape parameter RBF reproduction causes overshoots near steep gradients. This behavior will eventually lead to instabilities like other high order methods as the schemes marches in time. This is shown in Fig. 4. Next, the local blended RBF collocation method is used, allowing the scheme to switch between high or low shape parameter RBF interpolation vectors. Fig. 5 shows that the condition number of the collocation matrix is high near smooth regions. As the scheme marches towards the discontinuity, the scheme switches to a low shape parameter RBF and low condition number. For Hardy multiquadric RBF, a low shape parameter causes the interpolation to be piecewise linear, providing stability near steep gradients and shocks. In Fig. 5, the condition number varies between K = 6.4 − 6e6 depending on the smoothness of the function. Early in the development of the
𝜕𝑢 𝜕𝑢 +𝑎 = 0. (90) 𝜕𝑡 𝜕𝑥 For the first example, a top-hat profile is used for the initial condition shown in Table 2. The solution is evaluated for three different cases of low, high, and a blended shape parameter as shown in Fig. 6. The minmod flux limiter is used for the blended solution. The shape parameter values used are c = 0.5 and c = 0.001 for the high and low interpolation. The solution of the one-dimensional advection equation using the blended RBF collocation method works well capturing the gradient as close to the high shape parameter solution while using the low shape parameter solution when oscillations are detected. Also, the results
Fig. 4. Solution sweep with a constant shape parameter interpolation approach. Condition numbers are a) At top of discontinuity: 2.5 × 107 b) At bottom of discontinuity: 2.2 × 1010 . 87
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Fig. 5. Solution sweep with blended interpolation approach. Condition numbers a) At the top of discontinuity: 6.4 b) At the bottom of discontinuity: 11.
Fig. 6. 1-D Advection using Example 1 initial conditions.
88
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Fig. 7. 1-D Advection equation using Example 2 initial conditions.
Table 3 Example 2 initial conditions. u(x, u(x, u(x, u(x,
0) = 0 0) = 2 0) = 1 0) = 0
x < 0.1 0.1 ≤ x ≤ 0.3 0.3 < x ≤ 0.5 x > 0.5
provide proof that the scheme is a high-order scheme for hyperbolic partial differential equations. To demonstrate the robustness of the blending approach a different set of initial conditions are used giving the u(x) profile given in Table 3 and the results are shown in Fig. 7. Applying this methodology to solve the inviscid compressible flow governing equations yields stable solutions comparing well to theoretical values. The governing equations for the inviscid Euler equations are 𝜕𝑄 𝜕𝐹 𝜕𝐺 + + =0 𝜕𝑡 𝜕𝑥 𝜕𝑦
Fig. 8. Compressible Euler equations geometry and boundary conditions.
The inlet boundary conditions are Mach number of 2, pressure is 100,000 Pa, density is 1.16 kg/m3 and the temperature is 300 K. The walls of the domain is set to a slip boundary condition forcing the velocity to be tangential to the walls, and the outflow boundary for u-velocity, v-velocity, density, pressure, and temperature are extrapolated from the interior of the domain. The problem setup is shown in Fig. 8. The solution using the blended RBF interpolation with AUSM+ upwinding is shown in Fig. 9 and the numerical comparisons to theoretical oblique shock relations is shown in Table 4. The
(91)
where ⎧ ⎫ ⎧ ⎫ ⎧ ⎫ 𝜌𝑢 𝜌𝑣 ⎪𝜌⎪ ⎪ ⎪ ⎪ ⎪ 2 ⎪ 𝜌𝑢 ⎪ ⎪ ⎪ ⎪ ⎪ 𝜌𝑢𝑣 𝑄 = ⎨ ⎬ 𝐹 = ⎨ 𝜌𝑢 + 𝑃 ⎬ 𝐺 = ⎨ ⎬ 2 ⎪ 𝜌𝑣 ⎪ ⎪( 𝜌𝑢𝑣 ) ⎪ ⎪( 𝜌𝑣 + 𝑃) ⎪ ⎪𝜌𝑒𝑡 ⎪ ⎪ 𝜌𝑒𝑡 + 𝑃 𝑢⎪ ⎪ 𝜌𝑒𝑡 + 𝑃 𝑣⎪ ⎩ ⎭ ⎩ ⎭ ⎩ ⎭ 89
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
Fig. 10. Locations for analytical and numerical comparisons.
Fig. 9. Inviscid compressible flow using RBF blended interpolation.
Table 4 Comparison of results for localized blended RBF collocation. Oblique shock relations Location
Mach
Pressure (Pa)
Density (kg/m3 )
Temperature (K)
1 2 3
2.000 1.641 1.285
1.00E+05 1.71E+05 2.80E+05
1.161 1.693 2.405
300.000 351.045 405.958
Fig. 11. Inviscid compressible flow over supersonic airfoil.
Localized blended RBF collocation Location
Mach
Pressure (Pa)
Density (kg/m3 )
Temperature (K)
1 2 3
2 1.641 1.286
1.00E+05 1.71E+05 2.80E+05
1.161 1.694 2.406
300 351.05 405.9
oblique shock and expansion along the surface of the airfoil giving good agreement with theoretical shock relations.
8. Viscous compressible flow results
Finite volume method (FVM) Location
Mach
Pressure (Pa)
Density (kg/m3 )
Temperature (K)
1 2 3
2 1.63953 1.2887
1.00E+05 1.70E+05 2.79E+05
1.16198 1.68475 2.40177
300 350.544 405.132
The blended RBF interpolation scheme gave satisfactory results for the Euler equations, and we would like to add more physics to test the scheme further. The next step is to add viscosity to the fluid flow; therefore, we would like to solve the Navier-Stokes equations. Due to the nature of high Reynolds number flows, the boundary is thin, and it is stated that solving the Navier-Stokes for high-speed flow is one of the most difficult problems in computational fluid dynamics. The governing equations can be written in vector notation as
locations within the domain where comparisons are taken are shown in Fig. 10. Other geometries were used to investigate the robustness of the blended RBF interpolation scheme. Fig. 11 shows the solution for a symmetrical supersonic airfoil. The scheme is capable of capturing the
𝜕 𝐹𝑣 𝜕 𝐺𝑣 𝜕𝑄 𝜕𝐹 𝜕𝐺 + + = + 𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 90
(92)
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
where ⎧ ⎫ ⎧ ⎫ ⎧ ⎫ 𝜌𝑢 𝜌𝑣 ⎪𝜌⎪ ⎪ ⎪ ⎪ ⎪ 2 ⎪ 𝜌𝑢 ⎪ ⎪ ⎪ ⎪ ⎪ 𝜌𝑢𝑣 𝑄 = ⎨ ⎬ 𝐹 = ⎨ 𝜌𝑢 + 𝑃 ⎬ 𝐺 = ⎨ ⎬ 2 ⎪ 𝜌𝑣 ⎪ ⎪( 𝜌𝑢𝑣 ) ⎪ ⎪( 𝜌𝑣 + 𝑃) ⎪ ⎪𝜌𝑒𝑡 ⎪ ⎪ 𝜌𝑒𝑡 + 𝑃 𝑢⎪ ⎪ 𝜌𝑒𝑡 + 𝑃 𝑣⎪ ⎩ ⎭ ⎩ ⎭ ⎩ ⎭ and 𝑒𝑡 = 𝑒 +
) 1( 2 𝑢 + 𝑣2 2
𝑃 = 𝜌(𝛾 − 1)𝑒.
The viscous fluxes can be written in vector form as ⎧ ⎫ ⎧ ⎫ 0 0 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ 𝜏𝑥𝑥 𝜏𝑥𝑥 𝐹𝑣 = ⎨ 𝐺𝑣 = ⎨ ⎬ ⎬. 𝜏𝑥𝑦 𝜏𝑥𝑦 ⎪ ⎪ ⎪ ⎪ ⎪𝑢𝜏𝑥𝑥 + 𝑣𝜏𝑥𝑦 − 𝑞𝑥 ⎪ ⎪𝑢𝜏𝑥𝑦 + 𝑣𝜏𝑦𝑦 − 𝑞𝑦 ⎪ ⎩ ⎭ ⎩ ⎭ The shear stresses are given by Pletcher et al. [64] ( ) 2 𝜕𝑢 𝜕𝑣 τ𝑥𝑥 = 𝜇 2 − 3 𝜕𝑥 𝜕𝑦 τ𝑦𝑦 =
( ) 2 𝜕𝑣 𝜕𝑢 𝜇 2 − 3 𝜕𝑦 𝜕𝑥 (
τ𝑥𝑦 = 𝜇
𝜕𝑢 𝜕𝑣 + 𝜕𝑦 𝜕𝑥
Fig. 13. Domain for Navier-Stokes solution with CHT.
(93)
(94)
) (95)
𝑞𝑥 = −𝑘
𝜕𝑇 𝜕𝑥
(96)
𝑞𝑦 = −𝑘
𝜕𝑇 𝜕𝑦
(97)
where u and v are the velocity components, T is the temperature, q is the heat flux, 𝜇 is the dynamic viscosity, T is temperature, and k is the thermal conductivity. The inlet and outlet boundary conditions remain the same as the Euler equation examples and are provided in Fig. 12. The wall boundary condition is using the no-slip condition. Therefore, the velocity at the wall is zero. A conjugate heat transfer boundary condition is prescibed to the lower fluid domain wall, and the temperature is solved within the solid domain using RBF direct differentiation. The solid material chosen for this case is titanium. The boundary conditions for solid walls are prescribed to be adiabatic, and the lower wall of the solid region is at a constant temperature of 300 K. The solution domain is shown in Fig 13. The Navier-Stokes solution for pressure and temperature is shown in Fig. 14. Points are added near the boundary to capture the boundary layer and the approximate the shear stresses. To overcome the non-uniform
Fig. 14. Pressure (Pa) and temperature(K) results using blended RBF interpolation.
nodal arrangement, the local support of the RBF should be scaled based on the maximum distance of nodes within the local support [65–67]. This approach reduces the shape parameter dependency when using non-uniform node distribution. 9. Conclusions The localized blended radial basis function (RBF) collocation method with Moving Least Squares has shown to solve the Euler and Navier-Stokes equations by solving several numerical experiments. For smooth functions, the RBF interpolation using a high condition
Fig. 12. Boundary conditions for Navier-Stokes with conjugate heat transfer. 91
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
collocation matrix does well in reproducing a test function and determining the derivatives. As the shape parameter is increased the error for the function and derivatives decrease exponentially. For most problems, the high shape parameter RBF should be used frequently. However, increasing the slope of the function causes the high shape parameter RBF interpolation to produce oscillations. We seek to use the RBF interpolation for the solution of PDEs, therefore, the oscillations must be reduced or eliminated for a stable solution. Fortunately, the use of low shape parameter or low condition number RBF interpolation dampens oscillations providing a piecewise linear reproduction of the function. A robust high-resolution scheme is developed benefiting from the spectral accuracy capabilities of the RBF collocation method by allowing the numerical scheme to switch between high and low RBF interpolation,. The blended RBF scheme was able to solve many different PDE example problems, such as the inviscid Burgers equation, two-dimensional advection equation to demonstrate the performance of the blended RBF scheme. Studying the blended RBF scheme using these simple numerical test cases allowed for the experimentation of different blending parameter strategies and ultimately finding a criterion that eliminates oscillation while minimizing the numerical dissipation introduced into the solution. The 2-D inviscid compressible Euler equations and Navier-Stokes equations are solved using the blended RBF scheme. The results for all schemes showed good agreement with analytical shock relations and commercial CFD solutions. The localized blended RBF collocation method provides a stable method that can be used to solve hyperbolic PDEs and introduces a new tool for future research. The many advantages of using meshless methods can be used for hyperbolic PDEs. For example, the capability of refining the solution by adding nodes during the computation can now be used effectively.
[22] Hardy RL. Multiquadric equations of topography and other irregular surfaces. J Geophys Res 1971;76(8):1905–15. [23] Hardy RL. Research results in the application of multiquadric equations to surverying and mapping problems. Surv Mapp 1975;35:321–32. [24] Franke R. Scattered data interpolation: tests of some methods. Math Comput 1982;38(157):181–200. [25] Dyn N, Levin D, Rippa S. Numerical procedures for surface fitting of scattered data by radial basis funtions. SIAM J Sci Stat Computat 1986;7(2):639–59. [26] Karageorghis A, Chen CS, Smyrlos Y-S. A matrix decomposition RBF algorithm: approximation of functions and their derivatives. Appl Numer Math 2007;57:304–19. [27] Sarler B, Vertnik R. Local explicit radial basis funtion collocation method for diffusion problems. Computat Math 2005:1269–82. [28] Yao G, Sarler B, Chen CS. A comparison of three explicit local meshless methods using radial basis funtions. Eng Anal Bound Elem 2011;35:600–9. [29] Lee CK, Liu X, Fan SC. Local multiquadric approximation for solving boundary value problems. Comput Mech 2003;30(5–6):396–409. [30] Gerace S, Erhart K, Kassab A, Divo E. A model-integrated localized collocation meshless method for large scal three dimensional heat transfer problems. Eng Anal 2014;45:2–19. [31] Hon Y-C, Sarler B, Yun D-f. Local radial basis function collocation method for solving thermo-driven fluid-flow problems with free surface. Eng Anal Bound Elem 2015;57:2–8. [32] Kelly J, Divo E, Kassab AJ. Numerical solution of the two-phase incompressible Navier-Stokes equations using a GPU-accelerated meshless method engineering analysis with boundary elements. Eng Anal 2014;40C:36–49. [33] Sarler B, Tran-Cong T, Chen CS. Meshfree direct and indirect local radial basis function collocation formulations for transport phenomena. Bound Elem XVII 2005:417–28. [34] Erhart K, Kassab AJ, Divo E. An inverse localized meshless technique for the determination of non-linear heat generation rates in living tissues. Int J Heat Fluid Flow 2008;18(3):401–14. [35] Cheng AHD, Golberg MA, Kansa EJ, Zammito G. Exponential convergence and H-c multiquadric collocation method for partial differenctial equations. Numer Methods Partial Differ Equ 2003:571–94. [36] Madych WR. Miscellaneous error bounds for multiquadric and related interpolations. Comput Math Applic 1992;24(12):121–38. [37] Buhmann MD, Dinew S, Larsson E. A note on radial basis funtion interpolant limits. IMA J Numer Anal 2010;30:543–54. [38] Madych WR. Bounds on multivariate polynomials and exponential error estimates for multiquadric interpolation. J Approx Theory 1992;70:94–114. [39] Roque CMC, Ferreira AJM. Numerical experiments on optimal shape parameters for radial basis functions. Numer Methods Partial Differ Equ 2009:1–14. [40] Wright GB, Fornberg B. Scattered node compact finite difference-type formulas generated from radial basis funtions. J Comput Phys 2006;212:99–123. [41] Larsson E, Fornberg B. A numerical study of some radial basis function based solution methods for elliptic PDEs. Comput Math Appl 2003;46:891–902. [42] Wang JG, Liu GR. On the optimal shape parameters of radial basis functions used for 2-D meshless methods. Comput Methods Appl Mech Engrg 2002;191:2611–30. [43] Fasshauer GE, Zhang JG. On choosing "optimal" shape parametes for RBF approximations. Numer Algor 2007;45:345–68. [44] Schaback R. Error estimates and condition numbers for radial basis function interpolation. Adv Comput Math 1995;3(3):251–64. [45] Driscoll TA, Fornberg B. Interpolation in the limit of increasingly flat radial basis function. Comput Math Appl 2002;43(3–5):413–22. [46] Buhmann M, Dyn N. Spectral convergence of multiquadric interpolation. In: Proceeding of the Edinburgh mathematical society, Edinburgh; 1993. [47] Larsson E, Fornberg B. Theoretical and compurtational aspects of multivariate interpolation with increasingly flat radial basis funtions. Comput Math Appl 2005;49:103–30. [48] Carlson RE, Foley TA. The parameter R2 in multiquadric interpolation. Comput Math Applic 1991;21(9):29–42. [49] Carlson RE, Foley TA. Interpolation of track data with radial basis methods. Comput Math Applic 1992;24(12):27–34. [50] Jung J-H, Durante VR. An iterative adaptive multiquadric radial basis function method for the detection of local jump discontinuities. Appl Numer Methods 2008;59:1449–66. [51] Harris M, Kassab A, Divo E. Application of an RBF blending interpolation method to problems with shocks. Comput Assisted Methods Eng Sci 2015;22:229–41. [52] Harris M, Kassab E, Divo E. An RBF interpolation blending scheme for effective shock-capturing. Int J Comput Methods Exp Meas 2017;5(3):281–92. [53] Harris MF, Kassab AJ, Divo EA. Meshless method using a RBF interpolation scheme for inviscid compressible flows. In: Second thermal and fluids engineering conference, Las Vegas, NV; 2017. [54] Harris MF, Kassab AJ, Divo EA. Viscous compressible flow soluiton using RBF blended interpolation. In: 3rd thermal and fluids engineering conference (TFEC), Fort Lauderdale, FL; 2018. [55] Pepper DW, Kassab AJ, Divo EA. Introduction to finite element, boundary element, and meshless methods. ASME Press; 2014. [56] Laney CB. Computational gasdynamics. Cambridge University Press; 1998. [57] Chung TJ. Computational fluid dynamics. New York, NY: Cambridge University Press; 2002. [58] Moukalled F, Mangani L, Darwish M. The finite volume method in computational fluid dynamics: an advanced introduction with OpenFOAM and Matlab. Switzerland: Springer International Publishing; 2016.
References [1] Liu GR. Mesh free methods. Boca Raton, FL: CRC Press; 2003. [2] Kolodziej JA, Grabski JK. Many names of the Trefftz method. Eng Anal Bound Elem 2018;96:169–78. [3] Chen W, Hon YC. Numerical convergence of boundary knot method in the analysis of Helmholtz, modified Helmholtz, and convection-diffusion problems. Comput Methods Appl Mech Eng 2003;192:1859–75. [4] Fairweather G, Karageorghis A. The method of fundamental solutions for elliptic boundary value problems. Adv Comput Math 1998;9:69–95. [5] Golberg MA. The method of fundamental solutions for Poisson’s equation. Eng Anal Bound Elem 1995;16:205–13. [6] Grabski JK, Kolodziej JA. Laminar flow of a power-law fluid between corrugated plates. J Mech Mater Struct 2016;11(1):23–40. [7] Chen W, Wang FZ. A method of fundamental solutions without fictitious boundary. Eng Anal Bound Elem 2010;34(5):530–2. [8] Belytschko T, Lu YY, Gu L. Element -free Galerkin methods. Int J Numer Methods Eng 1994:229–56. [9] Atluri T, Zhu T. A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998:117–27. [10] Melenk JM, Babuska I. The partion of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng 1996:289–314. [11] Nayroles B, Touzot G, Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. Comput Mech 1992:307–18. [12] Duarte CA, Oden JT. Hp clouds- an hp meshless method. Numer Method Partial Diff Eq 1996:673–705. [13] Liu WK, Jun S, Li S, Adee J, Belytschko T. Reproducing kernel particle methods for structural dynamics. Int J Numer Methods Eng 1995:1655–79. [14] Kansa EJ. Multiquadrics - a scattered data approximation scheme with applications to computational fluid dynamics - I. Surface approximations and partial derivative estimates. Comput Math Appl 1990;19(8/9):127–45. [15] Kansa EJ. Multiquadrics - A Scattered data approximation scheme with applications to compuational fluid dynamics-II. Comput Math Appic 1990;19(8/9):147–61. [16] Kansa EJ, Carlson RE. Improved accuracy of multiquadric interpolation using variable shape parameters. Comput Math Applic 1992;24(12):99–120. [17] Kansa EJ, Hon YC. Circumventing the ill-conditioning problem with multiquadric radial basis funtions: applications to elliptic partial differential equations. Comput Math Appl 2000;39:123–37. [18] Kansa EJ, Holoborodko P. On the ill-conditioned nature of C∞ RBF strong collocation. Eng Anal Bound Elem 2017;78:26–30. [19] Divo EA, Kassab AJ. Iterative domain decomposition meshless method modeling of incompressible flows and conjugate heat transfer. Eng Anal 2006;30(6):465–78. [20] Divo EA, Kassab AJ. An efficient localized RBF meshless method for fluid flow and conjugate heat transfer. ASME J Heat Transfer 2007;129:124–36. [21] Divo E, Kassab AJ. Localized meshless modeling of natural convective viscous flows. Numer Heat Transfer 2008;53:487–509. 92
M.F. Harris, A.J. Kassab and E. Divo
Engineering Analysis with Boundary Elements 109 (2019) 81–93
[59] Liu X, Osher S, Chan T. Weighted essentially nonoscillatory schemes. J Comput Phys 1994;115:200–12. [60] Liou M-S, Steffen CJ Jr. A new flux splitting scheme. J Comput Phys 1993;107:23–39. [61] Liou M-S. A sequel to AUSM: AUSM+. J Comput Phys 1996;129(0256):346–82. [62] Lancaster P, Salkauskas K. Surfaces generated by moving least squares methods. Math Comput 1981;37(155):141–58. [63] Harten A. A high resolution scheme for the computation of weak solutions of hyperbolic conservation laws. J Compuat Phys 1983;49(357).
[64] Pletcher RH, Tannehill JC, Anderson DA. Computational fluid mechanics and heat transfer. CRC Press Taylor & Francis Group, LLC; 2013. [65] Vertnik R, Sarler B. Solution of incompressible turbulent flow by a mesh-free method. Comput Model Eng Sci 2009;44(1):65–95. [66] u Islam S, Vertnik R, Sarler B. Local radial basis funtion collocation method along with explicit time stepping for hyperbolic partial differential equations. Appl Numer Math 2011;67:136–51. [67] Vertnik R, Sarler B. Simulation of continuous casting of steel by a meshless technique. Int J Cast Met Res 2009;22:311–13.
93