A compact RBF-FD based meshless method for the
incompressible Navier-Stokes equations
P P Chinchapatnam', K Djidjeli'^*, P B Nair^, and M Tan^*
'Centre for Medical Image Computing, University College London, UK
^Computational Engineering and Design Group, University of Southampton, UK ^Fluid Structure Interactions Group, University of Southampton, UK
The manuscript was received on 2 February 2009 and was accepted after revision for publication on 16 April 2009. DOI: 10.I243/I4750902JEMEI51
Abstract: Meshless methods for solving fluid and fluid-structure problems have become a promising alternative to the finite volume and finite element methods, h i this paper, a mesh-free computational method based on radial basis functions i n a finite difference mode (RBF-FD) has been developed for the incompressible Navier-Stokes (NS) equations i n stream function vorticity f o r m . This compact RBF-FD formulation generates sparse coefficient matrices, and hence advancing solutions will i n time be of comparatively lower cost. The spatial discretization of the incompressible NS equations is done using the RBF-FD method and the temporal discretization is achieved by explicit Euler time-stepping and the Crank-Nicholson method. A novel ghost node strategy is used to incorporate the no-slip boundary condifions. The performance of the RBF-FD scheme w i t h the ghost node strategy is validated against a variety of benchmark problems, including a model fluid-structure interaction problem, and is f o u n d to be i n a good agreement vnth the exisfing results. I n addition, a higher-order RBF-FD scheme (which uses ideas f r o m Hermite interpolafion) is then proposed for solving the NS equations.
Keywords: meshless method, radial basis funcfions, finite difference, incompressible Navier-Stokes equations, fluid-structure interaction, stream f u n c d o n
1 INTRODUCTION
I n recent years, there has been an upsurge i n interest i n the development of so-called mesh-free methods as an alternative to the mesh-based methods, due to their potential i n alleviating the mesh-generation complexities arising i n traditional methods such as finite difference (FD), finite volume (FV), and finite element (FE) methods [1]. One c o m m o n character-istic among all the mesh-free methods is that they can construct the functional approximation or interpolation entirely f r o m the information at a set of scattered/random nodes or points. The earliest example of meshless methods is perhaps the gen-eralized finite difference scheme [2]. Some of the well-known meshless methods are smooth particle
'Corresponding author: School of Engineering Sciences, Uni-versity of Soutiiampton, UniUni-versity Road, S017 IB], UK. email: kkel@soton.ac.uk
hydrodynamics (SPH) method [3], diffuse element method (DEM) [4], element free Galerkin method (EFGM) [51, reproducing Kernel particle method (RKPM) [6], partidon of unity method (PUM) [7], finite point method (FPM) [8], and local Petrov-Galerkin method (MLPG) [9]. One of the m a i n advantages of mesh-free methods is that it is computationally easy to add or remove nodes f r o m a pre-existing set of nodes, which is not the case for mesh-based methods, where the addition or removal of a point/element would lead to heavy remeshing and hence computationally be difficult to implement.
For many years, radial basis functions (RBFs) have been synonymous w i t h scattered data approxima-fion, especially i n higher dimensions. I n recent years, there has also been an increased interest i n their use i n solving partial differential equations (PDEs) on irregular domains by a global collocation approach (see, for example, references [10] to [13]). Despite the excellent results of eariier works related
to the use of RBFs for the numerical solution of PDEs, the traditional RBFs are globally defmed functions, which result i n a dense linear system. This hinders tiie application of the RBFs to solve large-scale fluid dynamics problems (such as the Navier-Stokes (NS) equations), as they w i l l be computationally intensive, owing to tiie resulting matrices being fully populated. Also for numerical problems, which necessitate large numbers of points, these fully populated matrices tend to be ill-condi-tioned. Some attempts have been made to resolve tills problem [14-17] (and references therein). Among the recent proposed method is the local RBF-FD method, w h i c h is f o u n d to give spectral accuracy for a sparse, better-conditioned hnear system and more flexibility for handling non-linearities when used to solve elliptic problems [18]. This idea of using RBFs i n a finite-difference mode (RBF-FD), which can be seen as generalized finite differences w i t h arbitrary/random points instead of a regular grid system, was proposed by Wright and Fomberg [18], and Tolstykh and Shrio-bokov [19] independently i n the literature.
I n tills paper, the local RBF-FD collocation metiiod is extended for solving incompressible flows i n a stream function vorticity f o r m . A novel ghost centre strategy is employed to satisfy the boundary conditions. The performance of the local RBF-FD scheme and its higher-order w i t h the ghost node strategy is validated against a variety of benchmark problems, and compared vdth the existing results i n the literature. I n addition, a numerical study of near-bed submarine pipelines under current is presented to demonstrate the applicability of the RBF-FD scheme to model fluid-structure interactions.
which can be expressed i n terms of the stream function as
dé dé
Sy dx (3)
In the RBF-FD metiiod, tiie complete domain is represented by a set of scattered nodes present i n the interior and on the boundary. For each interior node, a supporting region/stencil is identified by choosing N nearest nodes. Then at each node, a local RBF interpolation problem is set up to determine the RBF-FD weights for each derivative.
The standard RBF interpolation is of the f o r m
u{x)Ks{x)=J2m\x-Xi\\) + [i ( 4 )
(=1
where ^ ( | | . | | ) is the multiquadric RBF and /? is a constant.
In Lagrangian f o r m , equation (4) can be wnritten as
1=1
(5)
where Z ( I | J : - X , | | ) is of tiie f o r m (4) and satisfies tiie usual cardinal conditions, i.e.
(6)
To derive RBF-FD formula at the node, say Xi (see Fig. 1), the differential operator is approximated using the Lagrangian f o r m of the RBF interpolant, i.e.
2 R B F - F D FOR INCOMPRESSIBLE NAVIER-STOKES EQUATIONS
The non-dimensional governing equations f o r un-steady incompressible NS equations expressed i n terms of vorticity (OJ) and stream f u n c t i o n W are given by dco dw doj -dï^^'-d^^'Ty 1 fd^co d^co\ " Re [dx^ ^ dy2) dx^ (1) (2)
where Re is the Reynolds number and u and v denote the components of velocity i n the x and y directions,
o 0 0 0 0 O 0 0 O 0 ° 0 ;• 'é o^- ° 0 0 0 0 O O 0 0 \ ° O .® 0 0 0-" 0-" \ ' Suppo/ting faglo 0 0 0 0 O n t e A C i 0 0 0 0 \ ° O 0 0 0 O 0 0 O ° O 0 0
Fig. I Schematic diagram of a higher-order RBF-FD stencil. The circle indicates the supporting region/stencil for the node Xi
Cu{xx) *Cs{x,) = C-/X\\xi -Xi\\)u(Xi) (7)
1=1
Equation (7) can be rewritten as a FD f o r m u l a of the form
1=1
(8)
where the RBF-FD weights are formally given by the operator £ apphed on the Lagrange f o r m of the basis functions, i.e.
11^(1.1)= ^^yX\\x\-Xi
' 0 e' w' " £ 0 i " .e^ OJ . / ' . 0
In practise, the weights are computed by solving the linear system
(10)
where £<I)i denotes the evaluation of the c o l u m n vectors C<^=[£<j){\\x-Xi\\) C(j>(\\x-X2\\) ••• C(l>(\\x-x„\\)]'^ at the node A T I . Here, / i is a scalar value and enforces the condition
n
(=1
which ensures that the stencil is exact for all constants.
Once the RBF-FD method is applied to discretize the spatial derivatives i n the governing equations, equations (1) to (2), at any interior node x,- the following are obtained
Re
7=1 J = l
(11)
and
(12)
where Nis the total number of interior and boundary nodes which lie i n the supporting region/stencil for the node x,-, and jt^g.,, zt/g,, i t / g j are the
RBF-FD weights obtained f r o m the system of equation (10) w i t h the corresponding differential operator iê/dx, d/dy, d^/dx^, S ^ / d f ) applied to the basis func-tions o n the right-hand side.
The system of ordinary differential equations obtained for vorticity after spatial discretization, equation (11), is advanced i n time using the basic Euler time-stepping scheme. Denoting the value of any physical quantity at ï = f" vwth the superscript n
dt
'Ye
7=1 7=1
U=i
(13)
where öt is the time-step. Similarly equation ( I I ) is temporally discretized using a 0-weighting scheme (0 < Ö < 1), the discretised equation at the node x, reads as öt . y=i + {l-0) R e Z ^ +"'(/.;) J'^y 7 = " 7=1 (14)
Equations (13) and (14) need to be supplemented by the boundary condition for vorticity. The value of vorticity at the boundary is obtained by higher-order finite difference expressions [20]; see Table 1. Here, the subscript b refers to the value of the quantity on the boundary and subscript 1 refers to the interior node, w h i c h is locally orthogonal to the boundary and at a distance Ii f r o m the boundary.
Once the value of vorticity i n the whole domain is obtained, the governing equation for the stream function, equation (2), is solved w i t h Dirichlet boundary conditions to update the stream function. This process is repeated u n t i l convergence
Jnew-ftJoldllz 0)
new||2
(15)
Table 1 0{h^) wall boundary conditions [20]
Left wall: a>i,=
-Right wall: o>i,= l : h + ( ^ ) - ^ ] Bottom wall:
Top wall: t'>b =
-where e is a predetermined convergence limit. The complete procedure is outlined in Table 2.
3 GHOST NODE STRATEGY FOR
INCORPORATING BOUNDARY CONDITIONS I n the previous section, a locally orthogonal grid at the boundary was used to enforce the no-slip boundary conditions. This restriction o n the nodes near the boundary makes the implementation of boundary conditions very straightforward. A con-siderable amount of work would be needed, how-ever, to ensure a locally orthogonal grid near curved surfaces, and hence this approach would be cum-bersome for complex geometries. I n this section, a method is proposed for implementing the no-slip boundary conditions based o n ghost nodes. This ghost node strategy enables randomly placed points near the boundary and is still able to satisfy the boundary conditions accurately. Sample point dis-tributions used i n the locally orthogonal grid and the ghost node strategies are shown i n Fig. 2.
The no-slip boundary conditions at a boundary F are given by
\]/ = Ci xeV
oil, (16)
dn = C2 x e r
where Ci and C2 are constants and n is the outward normal direcdon f r o m the boundary. I n the pro-posed strategy, each boundary node is associated w i t h a support region/stencil, w h i c h also includes a ghost node placed outside the computational
do-Locally Orthogonal Grid Ghost Nodes Grid
Fig. 2 Schematic figure depicting the locally orthogo-nal boundary and the ghost nodes. Note that the ghost nodes are represented as grey shaded circles
main. The RBF-FD discretization is carried out to approximate the normal derivative at the boundary node Xj, i.e.
cn = E ' ^ S i ^ + " S h o . , ' A g h o (17)
where N is the number of supporting points inside and on the boundary. The value of the stream f u n c t i o n at the ghost node is evaluated by substitut-ing the no-slip boundary condition equation (16) and the stream f u n c t i o n values of the interior nodes evaluated at the previous time step i n equation (17). The value of vorticity on the boundaries can then be evaluated by RBF-FD discretization of equation (2) at the boundary node x,-.
4 NUMERICAL STUDIES
This section presents numerical studies conducted o n two test problems using the modified RBF-FD scheme w i t h the ghost node strategy. Further, a numerical study of a near-bed submarine pipeline is presented to illustrate the applicability o f t h e RBF-FD method to solve fluid-structure interaction problems.
Table 2 RBF-FD algorithm for incompressible Navier-Stokes equations Given an initial guess ij/" and w° and a particular node conflguration:
1. For each interior node, determine the support/stencil size.
2. Obtain the RBF-FD weights by solving the RBF interpolaüon problem.
3. Advance the vorticity solution to the next step using a suitable time-stepping algorithm.
4. Calculate the vorticity on the boundary using Table 1.
5. Solve equation (2) with Dirichlet boundary conditions for stream function to obtain the new stream function values. 6. Check for convergence. If converged, stop, else go to step 3.
4.1 Square driven cavity flow
Numerical studies conducted o n the lid-driven cavity flow problem i n a square [0, 1] x [0, 1] domain are first presented. The boundary conditions for this problem are given by
IIj = 0,'^=0 on . ï = O a n d x = l ^!, = 0,
1^=0
on y = 0<P = 0,'^ = l on y = l
(18)
The results obtained using the presented RBF-FD formulation are validated against the benchmark multigrid finite difference results obtained i n refer-ence [21].
The time-dependent f o r m of the governing equa-tions is used i n stream function-vorticity f o r m . The spatial discretization is done using the RBF-FD scheme, while the temporal discretization is carried out using the Crank-Nicholson method (equation (14) w i t i i 0 = 0.5), w i t h a time step öt=om. Both u n i f o r m and random point distributions are con-sidered and the flow problem is solved for three different Reynolds numbers (i?e= |100, 400, 1000}). Nine supporting points are used i n each RBF-FD stencil for discretizations of the f u n c t i o n derivatives and the value of shape parameter is obtained using the leave-one-out optimization strategy [22], for each RBF-FD stencil. To apply the no-slip boundary conditions, the ghost node strategy proposed i n section 3 is employed on the boundary RBF-FD stencils. This facilitates a complete random point distribution i n the interior of the domain.
The accuracy of the proposed ghost node strategy is first examined. Figure 3 shows the waU vorticity distribution obtained on the moving l i d for the
square lid-driven cavity flow problem at two differ-ent Reynolds numbers. A complete random distribu-tion of points without any restricdistribu-tion at the nodes near the boundary was considered for obtaining the results. The results are compared w i t h the wall vorticity values obtained by Ghia et al. [21] for the purpose of validation. From the figures, it can be seen that the obtained vorticity distribution agrees well w i t h the benchmark results.
I n Fig. 4, the stream f u n c t i o n and vorticity con-tours obtained using the RBF-FD method for Re = 100 are shovm. The results displayed are generated using 4 1 x 4 1 randomly spaced points. From the plot of the stream f u n c t i o n contours it can be seen that the secondary and tertiary vortices near the b o t t o m wall are also captured. It is worth noting that the global features of the flow were captured w i t h relatively small 21 x 2 1 distribution of points.
The comparison of velocity components at the horizontal and vertical centres o f the cavity w i t h those obtained by Ghia et al. [21] are displayed i n Fig. 5. The velocity profiles obtained using 3 1 x 3 1 and 4 1 x 4 1 u n i f o r m and random distribution of points are presented i n Fig. 5. It can clearly be seen that the velocity profiles are captured accurately as the number of points i n the domain is increased.
Figure 6 shows the stream f u n c t i o n and vorticity contours obtained for Re = 400 and 5 1 x 5 1 random point distribution.The comparison of velocity pro-files is presented i n Fig. 7. To capture the velocity profiles f o r Re = 400 accurately, a larger number of points (51x51) were needed as compared w i t h those required for Re =100 (41x41). However, the points required were much less than that of the second-order finite difference method, w h i c h required about 129x129 points i n order to capture the velocity profiles [21]. This may be ovdng to RBF basis w i t h an optimized shape parameter employed i n the RBF-FD
0 1 0 2 0 3 0 4 0 5 0 6 0 7 Ofl 0 9
(a)
1 0 1 0 ? 0 3 0 4 0 5 0 6 0 7 0 8
(b)
Fig. 3 Comparison of wall vorticity obtained using ghost nodes on the moving boundary of the square driven cavity flow with reference [21]. (a) Re = 100, (b) Re = 400
-oja -005
0 0 1 0 2 0 3 04 0 5 0 5 0 7 0 0 0 9 1
(a)
Fig. 4 Square driven cavity: Re = 100, contours of stream function and vorticity obtained using 41x41 random point distribution, (a) Stream function, (b) vorticity
method as compared w i t h the polynomial basis used i n reference [21].
Similar results for Re= 1000 are shown i n Figs 8 and 9. The primary, secondary, and tertiary vortices are captured satisfactorily. The velocity profiles obtained using 5 1 x 5 1 and 6 1 x 6 1 u n i f o r m and random point distributions are displayed i n Fig. 9.
I n comparison v^ath the global RBF collocation method developed by the authors [13], the RBF-FD method is able to provide similar accuracy, but at m u c h lower computational cost. This reduction i n computational cost is mainly attributable to the sparse structure of the coefficient matrices. It is also observed that, although the sensitivity of the shape parameter is reduced for the RBF-FD method, the shape parameter still influences the accuracy of the obtained solution, particularly when random node
^k- 31X31 Uii*>mi - - 31X31 RBfkl.ni 41X41 U r ^ f m — 41X41 RsnAOT A GI»aBia1 (a)
stencils are used for spatial discretization. Finally, additional computational cost incurred by this method vdth respect to a second-order finite difference method is the calculation o f the weights for each stencil, which involves an optimization strategy involving typically around 20 inversions of a 9 x 9 matrix i n the present paper.
4.2 Rectangular driven cavity flow
I n this subsection, the RBF-FD approach is applied w i t h the ghost node strategy for solving the driven cavity flow i n a rectangular cavity w i t h aspect ratio 2. The problem is defined and solved i n the rectangle 0 < A : < l , 0 < y < 2 . This problem is solved for three different Reynolds numbers of 100, 400, and 1000.
31X31 Urifotm — - 31X31 Random — 41X41 U n i f o m — 41X41 Random G H a e t a i . 01 02 03 04 0 5 oa 07 oh 09 (b)
Fig.5 Square driven cavity: comparison of velocity profiles obtained on the vertical and horizontal centre-lines using RBF method with reference [21] for Re= 100. (a) Vertical centre-line, (b) horizontal centre-line
s -0,03 -0,03 -0.05 -0C5 ,09' -0,07 9 -tl i ^ ê p ? p » S ' ' P P ^Qog
-^ 1
Fig. 6 Square driven cavity: Re = 400, contours of stream function and vorticity obtained using 51x51 random point distribution, (a) Stream funcdon, (b) vorticity
- * - 4 t X 4 t U n l t a m - - 41X41 Random 51X51 Uniform - - 51X51 Random <} G h i a « a l - * - 41X41 Uniform - - 41X41 Random — 51X51 Uniform 51X51 Random O G h i a e l al.
Fig. 7 Square driven cavity: comparison of velocity profiles obtained on the vertical and horizontal centre-lines using RBF method udth reference [21] for i?<? = 400. (a) Verdcal centre-line, (b) horizontal centre-line
Fig. 8 Square driven cavity: Re = 1000, contours of stream function and vorticity obtained usine 61x61 random point distribution, (a) Stream fiinction, (b) vorticity
JEMElSl
(a) (b)
Fig. 9 Square driven cavity: comparison of velocity profiles obtained on the vertical and horizontal centre-lines using RBF method with reference [21] for Re = 1000. (a) Vertical centre-line, (b) horizontal centre-line
The results are validated against those obtained by Gupta and Kalita [23].
The stream f u n c t i o n contours obtained for die three Reynolds numbers are shown i n Figs 10 and 11. From the figures, it can be observed that there are two rotating primary vorfices as well as second-ary vortices i n the bottom corners of the rectangular cavity. The top primary vortex properties are reported i n Table 3, and are compared w i t h those obtamed by Bruneau and Jouron [24]. It can be seen that the RBF-FD method results are i n close agreement w i t h the benchmark results.
01 0 2 0 3 04 Oil 06 0 7 Ofl 0 9
Fig. 10 Rectangular driven cavity: streamline patterns obtained for Re =100 using 41x81 uniform point distribution
4.3 Near-bed pipeline under current
I n titis section, the RBF-FD scheme is appUed to conduct a stability analysis of a near-bed pipeline under current. When the pipeline is placed near the seabed the fluid current across the pipeline is asymmetrical, thus causing a fluid force to exert on die structure, resulting i n the deformation of the pipeline. This problem can be simplified as a beam w i t h different boundary conditions and a non-linear fluid forcing term (see Fig. 12). Based on the Timoshenko beam theory, the general governing equations of the pipehne can be written as
d r /'dw \ ox) dx' + / = 0 dw dx + 0 =0 (19) where w is the deflection of the pipeline, 0 the rotation, G the shear modulus, A cross-section area, its the shear correction coefficient, E the elasticity modulus, I the moment of inertia. f[x) is the fluid force caused by the current pressure difference. Although i t has an analytical solution, i t converges slowly. Lam et al. [25] proposed an approximate rational expression for this fluid force i n order to speed up computations and is given by
f{x)=-pAU^cid) (20)
where p is the mass density, UQ the current velocity, c
' -"-Vos ^ " „1 . -oji?:" -"0? \ 01 0 ? 0 3 OA 0 5 0 6 0 ? (a) » 01 0 2 0 3 0 4 0 5 0 6 0 7 OH 0 9 1 tb)
I
i
Fig. 11 Rectangular driven cavity: streamline patterns obtained for Re = 400 using 81x161 uniform points and Re = 1000 using 101 x201 uniform points, (a) Re = 400, (b) Re = 1000
Table 3 Rectangular driven cavity: top primary vortex strength and location and comparison with reference [24)
Reynolds number ipmin location
100 [24] -0.1033 (0.6172, 1.7344) 400 Present method -0.1030 (0.625, 1.721) 400 [24] -0.1124 (0.5547, 1.5938) 1000 Present method -0.1120 (0.555, 1.6125) 1000 [241 -0.1169 (0.5273, 1.5625) Present method -0.1185 (0.525, 1.57) relation c(d) = 2.23rf^-H2.54f/-1-0.02 0.77rf' + 0.44rf2-f-0.02rf (21)
where d is defined as rf= (Do - wix) - i?s)/(2i?s). and Do is the distance between the central line of pipehne at inidal status and the seabed, R^ the outer pipe radius. For the present numerical study, fixed boundary condirions are employed. Several physical and material parameters employed are obtained from Li et al. [26]. The spatial derivatives are approximated using RBF-FD scheme and the dis-cretized non-hnear problem is solved using fixed-point iteration. Figure 13 presents the results ob-tained by conducting the stability analysis. The
ft
critical velocity (f/c/, = 12.3m/s) of the current obtained using RBF-FD for instability failure analysis is very near to the one obtained i n reference [26] ( f / c , = 12.2m/s), thus demonstrating the ability of the RBF-FD m e d i o d to be applicable for problems involving fluid-structure interaction.
5 H I G H E R - O R D E R R B F - F D SCHEMES
Previous sections have demonstrated the applicabil-ity of RBF-FD schemes f o r die incompressible NS equations. A higher-order version of the RBF-FD scheme is now explored using ideas from Hermite interpolation. This higher-order discretization method using RBFs can be regarded as a general-ization of the Mehrstellenvarfahren introduced by Collatz [27] and later developed into compact FD formulae by Lele [28]. I n tiie compact FD methodol-ogy, for example, tiie partial derivative of an unknown f u n c t i o n w i t h respect to the x-coordinate at any grid point (/, j) is given by
Pu CX
('•y> t e { ( - l , / , , - + l }
Fig. 12 Schematic diagram for model fluid-structure
interaction problem ks{i-U+l]
ikj]
(22)
JEMElSl
Fig. 13 (a) Deflection and (b) stress distribudon along the pipeline (when Uo = 10 m/s and Do = 0.7 m); and (c) instability failure analysis carried out using RBF-FD method
The accuracy of the FD approximation is increased by adding the second term (derivative information), as shown i n equation (22). Note that this additional term does not change the stencil size at the grid point (/, J) i n the compact finite difference medio-dology. Higher-order RBF-FD methods were earlier used for the solution of linear and non-Unear Poisson problems [29]. I n diis work, this approach is extended to develop higher-order schemes for the incompressible NS equations. The higher-order accuracy of the presented f o r m u l a t i o n is demon-strated by solving the steady-state incompressible Navier-Stokes equations.
5.1 Basic formulation
The RBF-FD method generates a local RBF inter-polant f o r expressing the f u n c d o n derivatives at a node as a linear combination of the f u n c t i o n values on the nodes present i n the support region of the considered node. I n the spirit of compact FD schemes, the accuracy of the RBF-FD discretization can be increased by considering not only the fiinction values but also the derivative values on the nodes present i n the supportmg region. The weights of the higher-order stencil are computed using tiie Hermite interpolation technique.
First a brief introduction to the Hermite mterpolation metiiod is given. Let £ be an arbitirary linear differential operator and let // be a vector containing some combination of m < n distinct numbers from the set [1, 2, n}. The fiinction values u[x^\ are specified at each of tiie
n
distinct data points {Xi}"^•^. I n addition, data corresponding to die differential operator operat-ing on tiie fimction, £ii(jc„,), are specified at m points { x , „ } ) l i . Note tiiat tiie point set { x , „ } " j is a subset of tiie set {jc,};'=i. Then, die interpolant passing tiirough all the data can be written asu{x)^s{x)=Y,Ximx-Xi\\)
1=1
m
+ £ l , £ 2 < ^ ( | | x - . ï , „ | | ) - ^ / i (23) / = i
where £ 2 i ^ ( | | ||) is a basis f u n c t i o n derived by the functional £ acting o n the muhiquadric basis i^(ll.ll) as a f u n c t i o n of the second variable (centre) and ^ is a constant. The u n k n o w n coefficients are obtained by enforcing the conditions s{x,) = «(x,), i=\, n; £ s ( x , „ ) = £ u ( x j , / = l , • • • , m ; and E " = i ^ w = 0. I m -posing tiiese conditions leads to the following block linear system of equations
" CD £2'I> e" ' X ll £ 0 ££2<1) 0 I Cu 0^ 0 J . 0 where = H\\xi-xi\\), i j = h--,n C2<t>ij = £ 2 ( ^ ( | | ^ / - ^ . , ; | | ) . 1 = 1. • • - , « . j = L C<b,j = /:4>{\\x„-xj\\). i = l , y = ! , • • • , « £ £ 2 « > i j = CC2<l>{lx„-x,„\\), i=l, - -,m, j=l, --.rn
and e, = 1,1 = 1, n- Equation (24) is solved using a backwards substitution routine.
I n Lagrange f o r m , the Hermite mterpolant can also be w r i t t e n as 5 ( A : ) = X ^ X ( I I ^ - ^ . - | I ) " ( ^ / ) 1=1 m + J:X{\\^-X4)CU{X,„) (25) /=i
where x(\\x-Xi\\) and 2 ( | | A : - X , „ | | ) are of die f o r m equation (23) and satisfy the cardinal conditions, i.e.
z ( i k . - x , i i ) = { J ; f c = i , - , « (26)
^ ; ^ ( I K - ^ ' l l ) = 0 ' ^ ' = 1 ' (27)
and
H\\^k-x,„\\)=0, k=l, •••,n (28)
^ X ( | K - A : , „ | | ) = | ^ k=l....,m (29)
Equation (25) is the basis for deriving higher-order RBF-FD stencils. Consider the node Xi w i t h its support region containing ?i points, denoted by a dashed circle around Xi, as i n Fig. 1. The goal is to obtain a higher-order RBF-FD discretization of £ii{xi). The nodes i n the support region, which are shaded grey, are those nodes where both the f u n c t i o n values iii(.x)) and the functional values (Cu{x)) are used, i.e. the set ;/ of cardinality, say,
The higher-order RBF-FD discretization f o r Cu{xi) is given by the Lagrange f o r m o f the interpolant, i.e.
Cu{xi)^£s{xi)=Y,CyX\\x,~xi\\)u{xi) 1=1
m
-f ^ £ 2 ( | k i - x „ | | ) £ » ( x , „ ) (30) ; = i
Equation (30) can be rewritten as a compact FD f o r m u l a of the f o r m
" m
Cuixi)^
E
+
E"fi.O-^^K')
(31) '•=1 / = !where tiie weights for higher-order RBF-FD ^""^ { " ' f i - ' ) } / ^ ! " ^ " ^ 8iven by
< , , - ) = ^ z ( l k i - - i ; , i l ) , wfij^=£x{\\xi-xj) (32)
where the superscript £ on tiie weights denote tiiat the higher-order RBF-FD weights are computed f o r that particidar operator.
I n practise, the weights are computed by solving the linear system
£2$
e"w'
£(D
££2a>
0iv
= £*4)i0^ 0 0
where C*<S>i and £ * ö i denote the evaluation o f the c o l u m n vectors £*a) = [ £ 0 ( | | x - j c i | | ) £ ( ^ ( | | x - X 2 | | ) •••
m\x-x„\\f,
and r Ö = [ £ £ 2 < ^ ( | | A ; - J C , „ | | )£ £ 2 ^ ( | | x - . - c „ J | ) ••• £ £ 2 . ^ ( | | x - j c „ J | ) ] ^ at tiie node
Xl. Here, / ( i s a scalar value related to the constant /i
i n equation (23) and enforces the condition
n
E"'(i')=o
1=1
w h i c h ensures that the stencil is exact f o r all constants.
Once the weights are computed by solving equa-tion (33) for each node, they can be stored and used to discretize the partial differential equation i n a similar manner as i n the compact FD schemes.
6 HIGHER ORDER R B F - F D FOR T H E INCOMPRESSIBLE NS EQUATIONS
The previous section outiined the idea of using the Hermite interpolation technique to obtain a higher-order RBF-FD discretisation f o r each RBF-FD stencil. It can be observed that a family of higher-order schemes can be derived by defining w h i c h operator £ one is using for the higher-order RBF-FD stencil. This section presents one such formulation for the steady-state incompressible NS equations. This f o r m u l a t i o n has the advantage of easier i m p l e m e n -tation of the no-slip boundary conditions. First the governing equations of steady incompressible flows i n streamfunction (i/^) - vorticity (OJ) f o r m u l a t i o n are recalled
V^ij,= -w (34)
V o j = fle(^u-+._j (35)
where Re is the Reynolds number and (h , V) are tiie Cartesian velocity components of the flow.
As usual, the entire d o m a i n is discretized into a set of interior and boundary nodes and determine the stencil at each node. Each higher-order RBF-FD stencil contains n nodes and the vector 17 of cardinality m ^ n. The stencil i n f o r m a t i o n consists
JEME151
of the function values (tp{x) or (o{x)} at each of the n nodes and the functional information {C\l/{x) or Co}{x)) on each of the m nodes. Note that the operator £ is arbitrary.
The higher-order RBF-FD discretizations for the Laplacian of stream f u n c t i o n at each interior node Xj is given by
n m
f 2 1 "
where the higher-order RBF-FD weights Y^Ji.j)}.^^ and {"^^w)}"' J are obtained using equation (24) w i t h the operator C = V^.
Similarly the discretization of Laplacian of vorti-city can be obtained as
It m
;=1 /=1
For tiie sake of brevity, the value of any physical quantity at node x, is denoted by the subscript j. I n equation (36), the second term (quantity i n curly brackets) can be replaced by the right-hand side of equation (34). Similarly, i n equation (37), the second term can be replaced by the right-hand side of equation (35).
The m o d i f i e d discretizations now become
where zï,- and i?,- are the current estimates of components of the velocity vector. Substituting the derived higher-order RBF-FD discretizations f o r Laplacian of streamfunction (equation (38)) and vorticity (equation (39)) into equation (40)
and
-Re
7=1 7=1
= 0
(42)
Note that i n equation (42), the vorticity gradients ^ ^ were discretized using the RBF-FD method. The velocity components (z7, v) i n equation (42) are obtained using the higher-order RBF-FD dis-cretizations given by " - dé * E O y + E H - , / ) ö ^ /• ; = 1 / = 1 ^ (43) 'll V=l Z=l (38) and " „2 ^ , / d(o dco\ (39)
The solution of the governing NS equations via a fixed-point iteration scheme is now returned to. Denoting the iteration number fc w i t h a superscript fc o n the physical variable, the governing equations at iteration fc-(-1 for the node x,- are given by
2,/,t+l _
Ui C(0 k+l Ö0)\t+r
(40) "ax n _ dé L/=i '=1 (44)where ij/ is the current estimate of the stream f u n c t i o n .
The iteration procedure is explained for a problem with slip boundary conditions. Recall that the no-slip boundary condition consists of a Dirichlet and a Neumann condition f o r stream f u n c t i o n at each boundary point, see equation (16). Now, given an initial guess for stream f u n c t i o n and vorticity, i t is possible to solve the system o f equations arising from satisfying equation (41) at all interior nodes along w i t h the Dirichlet boundary conditions f o r the stream fiinction to obtain the new estimate for stream func-tion (lA). To obtain the new velocity vector esthnate, the system of equations arising from satisfying equations (43) and (44) at all interior nodes is solved. Note that, whenever the support point x„, for a node is on the boundary, the Neumann condition for stream
f u n c t i o n is used thus facilitating an easier implemen-tation of no-slip boundary conditions. Next, the new estimate of vorticity on the boundary is obtained using the ghost node strategy proposed i n section 3. Once the velocity vector estimate is known, the linear system of equations arising from satislying equation (42) is solved w i t h the Dirichlet vorticity conditions obtained using the ghost node strategy. Once the physical quantities are obtained, advance to the next iteration. This procedure is repeated until convergence.
6.1 Numerical results
Numerical results obtained f o r the higher-order RBF-FD method f o r the steady convection diffrision equation and the incompressible NS equations are now presented.
First the higherorder RBFFD approach is i l l u -strated for the steady-state convection-diffusion equation of the f o r m
(45)
i n the domam [0, l ] x [ 0 , 0.6] w i t h the boundary conditions « = 1 o n x = 0 , u = 2 o n j : = l fu cy = 0 o n y = 0. du f y = 0 o n y = l (46) (47)
The exact solution for this problem is given by
_ l - e x p [ P e ( x - l ) ]
" ' " ' ^ ' - ^ l - e x p ( - P , ) (48)
Figure 14 presents the convergence plots of the RBF-FD and higher-order RBF-RBF-FD f o r the convection diffusion problem, equation (45), f o r two Peclet numbers 1.0 and 10.0. The operator C is taken as
dx^^By^ "^'dx
and the higher-order RBF-FD weights are obtained using equation (33). From Fig. 14, it can be seen that the results obtained using the higher-order method are at least two orders more accurate than that of RBF-FD method. It is also worth m e n t i o n i n g that f o r Pe = 100.0, w i t h a u n i f o r m discretization of 201x201 nodes, the error n o r m observed f o r the computed solution using higher-order RBF-FD was £ = 0(10"^).
Next, t u r n i n g to the incompressible NS, the f o r m u l a t i o n outUned i n section 6 is used to solve the square lid-driven cavity flow. Figure 15 presents the stream f u n c t i o n contours obtained f o r the cavity flow at Re = 100 w i t h a u n i f o r m distribution o f 3 1 x 3 1 nodes. The left subfigure solution is obtained using the RBF-FD method and the right subfigure is obtained using the higher-order RBF-FD method. From Fig. 15, it can be clearly seen that the higher-order method captures the solution more accurately w i t h a small number of 31 x31 points.
Fig. 14 Comparison of convergence behaviours of RBF-FD and higher-order RBF-FD for a model steady state convecdon-diffiision equation, h is the mesh spacing and e denotes the L.^ norm between exact and computed solutions, (a) = 1.0, (b) Pe = 10.0
h -003 -005 -o.os -ooj ^ -0.0' -0 03 -0.05 —oo?_
I
I
01 0 2 0 3 0 4 0 5 0 6 0 / 0 0 0 9 (a)Fig. 15 Comparison of convergence beiiaviours of (a) RBF-FD and (b) higher-order RBF-FD for a model steady-state driven cavity problem
Figure 16 estimates the performance of higher-order RBF-FD and RBF-FD i n terms of accuracy. O n the X-axis, the mesh spacing h is plotted o n a log scale i n the reverse direcdon. O n the y-axis, the m i n i m u m value of stream f u n c t i o n i/'min i n the whole domain (strength of the p r i m a r y vortex) is plotted. The benchmark value obtained by Ghia et al. [21] is shown as a horizontal dotted line i n the figure. F r o m Fig. 16, i t can cleariy be observed that the higher-order method captures die true solution at con-siderably less points (h = 0.033) as compared w i ü i the original RBF-FD m e ü i o d .
7 CONCLUDING REMARKS
The RBF-FD method is presented f o r solving incompressible Navier-Stokes equadons. This method approximates the f u n c t i o n derivatives at a node i n terms of the f u n c t i o n values o n a scattered set of points present i n support region o f the node. The RBF-FD method uses local interpolation pro-blems and hence generates sparse and well-condi-tioned matrices. It also has the property of decreased sensitivity vdth respect to shape parameter value i n comparison w i t h the RBF coUocadon method. A
-0.07 -0.075 -0.08 -0.085 -0.09 -0.095 -0.105 -0.11 0.1 Higher Order RBF-FD RBF-FD Ghia et al. 0.0333 0.025 0.02 0.0167
h
0.01Fig. 16 Convergence of i^-mm for higher-order RBF-FD and RBF-FD at Reynolds number 100
ghost node strategy employed for incorporadng no-slip boundary conditions removes the hmitations of having a locally orthogonal grid near the boundary and thus makes the method more suitable f o r complete random node discretisations. Numerical studies conducted o n the driven cavity flow pro-blems and the flidd-structure interaction of near-bed submarine pipelines using the RBF-FD method show that this method achieves accurate results w h i c h are i n good agreement w i t h the benchmark results.
A higher-order RBF-FD method is explored for solving partial differential equadons. The higher-order method is obtained by using Hermite RBF interpolation method to construct the f u n c t i o n approximation at each node i n the domain. A higher-order formulation for steady incompressible Navier-Stokes equations is presented. The accuracy of the higher-order method is investigated by solving f o r a model convection diffusion equation and square lid-driven cavity flow. Numerical results obtained indicate that this method indeed is a higher-order method w i t h a higher capability of spatial resolution w i t h respect to the RBF-FD method.
© Authors 2009
R E F E R E N C E S
1 Belytschko, T., Krongauz, Y., Organ, D., Fleming, M., and Krysl, P. Meshless methods. An overview and recent developments. Computer Methods Appl. Mech. Engng, 1996, 139, 49-74.
2 Perrone, N. and Kao, R. A general finite difference method for arbitrary meshes. Computers and Structs, 1975, 5, 45-47.
3 Lucy, L. B. A numerical approach to the testing of fission hypothesis. Astronomical J., 1977, 8, 1013-1024.
4 Nayroles, B., Touzot, G., and Villon, P. General-izing the finite element method: diffuse approx-imation and diffuse elements. Comput. Mech., 1992, 10, 307-318.
5 Belytschko, T., L u , Y., and Gu, L. Element free Galerkin methods. Itit. J. Numer. Methods Engng, 1994, 37, 229-256.
6 Liu, W. K., Jun, S., and Belytschko, T. Reproducing Kernel particle metiiods. Int. J. Numer. Methods in Fluids, 1995, 20, 1081-1106.
7 Melenk, J. M. and Babuska, 1. The partition of unity finite element method: basic theory and applications. Technical Report, Texas Institute for Computational and Applied Mathematics, 1996. 8 Onate, E . , Idelsohn, S., Zienkiewicz, O. C , Taylor,
R. L., and Sacco, C. A stabilized finite point method
for analysis of fluid mechanics problems. Compu-ter Methods Appl. Mech. Engng, 1996, 139, 315-346. 9 Atluri, S. N. and Zhu, T. New meshless local Petrov-Galerkin approach in computational me-chanics. Comput. Mech., 1998, 22, 117-127. 10 Kansa, E . J . Multiquadrics-a scattered data
ap-proximation scheme with applications to compu-tation fluid dynamics-I. Suriface approximations and partial derivatives estimates. Computers Math, with Applie, 1990, 19(8/9), 127-145.
11 Wang, J. G. and Liu, G. R. On optimal shape parameters of radial basis functions used for 2D meshless methods. Computer Methods in Appl. Mech. Engng, 2002, 191, 2611-2630.
12 Shokri, M. A numerical method for KdV equation using collocation and radial basis functions. Non-linear Itynamics, 2007, 50, 111-120.
13 Chinchapatnam, P. P., Djidjeli, K., and Nair, P. B. Radial basis function meshless method for the steady incompressible Navier-Stokes equations. Int. I. Computer Math., 2007, 84, 1509-1521. 14 Fasshauer, G. E . Solving partial differential
equa-tions with radial basis funcequa-tions: multilevel meth-ods and smoothing. Adv. in Comput Math., 1999, 11, 139-159.
15 Kansa, E. J. and Hon, Y. C. Circumventing the i l l -conditionmg problem with multiquadric radial basis functions: apphcations to elliptic partial differential equations. Comput Math. Appl, 2002, 39, 123-137.
16 Ling, L. and Kansa, E. J. A least-squares precondi-tioner for radial basis functions collocation meth-ods. Adv. Comput Math., 2005, 23, 31-54. 17 Chinchapatnam, P. P., Djidjeli, K., and Nair, P. B.
Domain decomposition for time-dependent pro-blems using radial based meshless methods. Numer. Methods for Partial Differential Equations, 2006, 23, 38-59.
18 Wright, B. W. and Fomberg, B, Scattered node compact finite difference-type formulas generated from radial basis functions, ƒ. Comput Phys., 2006, 212 99-123.
19 Tolystykh, A. 1. and Shirobokov, D. A. On using radial basis functions in a finite difference mode with applications to elasticity problems. Comput Mech., 2003, 33, 68-79.
20 Spotz, W. F. High-order compact flnite dijference schemes for computational mechanics, PhD Thesis, University of Texas at Austin, 1995.
21 Ghia, U., Ghia, K. N., and Shin, C. T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, ƒ. Comput Phys., 1982, 48, 387-411.
22 Chinchapatnam, P. P. Radial basis function based meshless methods for fluid flow problems, PhD Thesis, Southampton University, 2006.
23 Gupta, M. M. and Kalita, J. C. A new paradigm for solving Navier-Stokes equations: stream function-velocity formulation, ƒ. Comput. Phys., 2005, 207, 52-68.
24 Bruneau, C. H. and Jouron, C. An efficient scfieme for solving steady incompressible Navier-Stokes equations, ƒ. Comput. Phys., 1990, 89, 389^13.
25 Lam, K. Y., Wang, Q. X., and Zong, Z. A nonlinear fluid-structure interaction analysis of a near-bed submarine pipeline in a current. /. Fhiids and Structs, 2002, 16, 1177-1191.
26 Li, H., Cheng, J. Q., Ng, T. Y., Chen, J., and Lam, K. Y. A meshless Hermite-cloud method for nonlinear fluid-structure analysis of near-bed submarine
pipelines under current. Engng Structs, 2004, 26, 531-542.
27 Collatz, L. The numerical treatment of differential equations, 1960 (Springer Verlag, BerUn).
28 Lele, K. Compact finite difference schemes with spectral-like resolution, ƒ. Comput. Phys., 1992, 103, 16-42.
29 Wright, G. B. and Fomberg, B, Scattered node compact finite difference-type formulas generated from radial basis functions. /. Comput. Phys., 2006, 212, 99-123.
