Radial basis function Hermite collocation approach for the solution of time dependent convection–diffusion problems

Engineering Analysis with Boundary Elements 29 (2005) 359–370 www.elsevier.com/locate/enganabound

Radial basis function Hermite collocation approach for the solution of time dependent convection–diffusion problems A. La Rocca, A. Hernandez Rosales, H. Power School of Mechanical, Material and Manufacturing Engineering, The University of Nottingham, University Park, Nottingham NG7 2RD, UK Received 9 March 2004; accepted 23 June 2004 Available online 7 April 2005

Abstract This work presents a meshless numerical approach for the solution of time dependent convection–diffusion problems in terms of a Hermite radial basis function interpolation numerical scheme. To test the proposed scheme several numerical examples are analysed including problems with variable convective velocity and reaction coefficient. Comparisons are made with available analytical solutions. q 2005 Elsevier Ltd. All rights reserved. Keywords: Hermite radial basis functions; Convection–diffusion problems; Method of fundamental solution

1. Introduction The use of a mesh is a basic characteristic of traditional numerical approaches for the solution of partial differential equations. In those approaches, assumptions are made for the local approximation of the primitive variables, which require mesh to support them. During recent years, considerable effort has been given to the development of so-called mesh free methods (meshless approach). The aim of this type of approach is to eliminate at least the structure of the mesh and approximate the solution entirely using nodes values inside and/or in the boundary quasi-random distributed in the domain. For instance, the meshlees local Petrov-Galerkin and local boundary integral methods were given by Atluri et al. [1] and [2], respectively. These two methods basically transform the original problem into a local weak formulation over each subdomain, and the shape functions were constructed from using the moving leastsquares approximation to interpolate the solution. The method of fundamental solution (MFS) is one of the most intuitive meshless approaches to solve numerically linear elliptic partial differential equations. In this type of approach an approximate solution is represented in the form of a linear superposition of source functions (fundamental solutions) located outside the problem domain, U. As the fundamental solution satisfies the differential equation at any point except at the source point, it follows that this 0955-7997/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2004.06.005

representation exactly satisfies the governing equation and one seeks to satisfy the imposed boundary conditions approximately at a set of boundary points (collocation points). In Ref. [3] a comprehensive review of the MFS is given, in which this approach is related to a regular indirect single layer boundary element approach (indirect-BEM). One of the pioneer works on the MFS was the work by Kupradze and Aleksidze [4], where they prove the completeness of the potentials in H(U), the set of harmonics functions in U. The implementation of the MFS appears to be quite straightforward. It is only necessary to specify the source points outside the problem domain and the data or boundary points on the contour of the domain then solve the resulting linear system of equations. However, in its numerical application there are several issues that require special attention. Of those some of the most relevant are the location and number of sources and boundary points. Usually the point sources are distributed over an auxiliary surface enclosing the original problem domain and located at a distance R from it. From numerical experiences, it has been found that the distance R has to be of the order of five times the characteristic dimension of the problem domain, when dealing with pure Dirichlet condition, or smaller than five when dealing with Neumann cases. This fixed value of R is recommended because the system of equations coming from the MFS can become highly ill-conditioned as R increased [5,6]. However, it has been observed that despite

360

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370

the ill-conditioning, the accuracy of the numerical solution is largely unaffected [7–9]. Golberg and Chen [3] reported that in these cases direct solvers, as Gaussian elimination, can be used effectively, even if the resulting condition number of the matrix is very large, by ignoring the message of singular system of the numerical routine. The requirement of an auxiliary surface to avoid the singularity of the fundamental solution introduces an additional unknown in the MFS, which is not easy to be determined. Chen and Tanaka [10,11] recently proposed an alternative to the MFS, the boundary knot method (BKM), in which the use the interior fundamental solution of the original partial differential equation, i.e. singular at infinity, instead of the original exterior fundamental, i.e. singular at rZ0. In this way, the fundamental solution can be directly distributed over the problem boundaries and at the same place of the boundary evaluation points. The exterior and interior fundamental solutions for most elliptic differential equations are given in term of radial basis function, therefore these two methods, i.e. MFS and BKM, can be considered as a boundary-only radial basis function free mesh approaches. In this work, we used a truly domain meshfree numerical scheme based upon radial basis functions Hermite interpolation techniques. In recent years, the theory of radial basis functions (RBFs) has undergone intensive research and enjoyed considerable success as a technique for interpolating multivariable data and functions. Although most work to date on RBFs relates to scattered data approximation and in general to interpolation theory there has recently been an increased interest in their use for solving PDEs. This approach, which approximates the whole solution of the PDE by the direct use of RBFs, is very attractive due to the fact that this is truly meshfree technique. Kansa [12,13] introduced the concept of solving PDEs using RBFs direct interpolation approach. He focused upon the multiquadric function and argued that PDEs are intrinsically related to the interpolation scheme from which PDE solvers are derived. The above approach has been applied successfully in several cases (see for example [14–16]), including initial and boundary value problems. However, no existence of solution and convergence analysis is available in the literature and for some cases, it has been reported that the resulting matrix was extremely ill-conditioned and even singular for some distribution of the nodal points (see [17]). Several techniques have been proposed to improve the conditioning of the coefficient matrix and the solution accuracy, as are: replacement of global solvers by block partitioning, LU decomposition schemes, matrix preconditioners, multizone methods etc. (see [18]). Besides some variations of the original Kansa’s idea, which improve the performance of the scheme, have been also suggested. Among them it is worthy to mention the works of Shu and co-workers and Tran-Cond and co-workers. Wu and Shu [19] suggested using differential quadrature approximation to represent the differential operators of

a given problem, instead of the direct differentiation used in the Kansa’s method. The proposed differential scheme is the same one used in the dual reciprocity method (DRM) of the boundary element technique to approximate the derivatives of the variables in the volume integrals to be transformed to surface integrals (for more details see Partridge et al. [20]). In this type of differentiation, the basic unknowns are the nodal values of the function instead of the expansion coefficients. Wu and Hon [21] recently solved transient heat diffusion problems using a scheme similar to the one proposed by Wu and Shu. Besides showing the implementation of the approach, they also gave error estimate on the convergence of the scheme. In relation with the use of differential quadrature approximation in meshfree techniques, it is important to comment the work of Golberg and Chen [22]. In this article they extended the MFS to solve nonhomogeneous elliptic partial differential equations, in which the homogeneous problem is solved by using the MFS and the corresponding particular solution is found in terms of the DRM approximation. This technique, original proposed to solve Poisson problems, has been used by them in several articles to solve more complex PDEs. More recently, Hon and Chen [23] and Chen and Hon [24] extended the BKM to solve nonhomogeneous partial differential equations using the above decomposition approach, i.e. homogeneous and particular solutions, but instead of solving the homogeneous problem by the MFS they use the BKM. In these last two articles the authors solve several cases including nontransient nonhomogeneous convection–diffusion problems with constant velocity field. On the other hand, Mai-Duy and Trang-Cong [25] proposed to approximate directly the higher derivative of the differential operator by the radial basis function instead of the function itself and find the lower derivatives and the function via symbolic integrations, instead of finding the derivatives by differentiation of the radial basis function. By this idea, they were able to improve the accuracy of the scheme without increasing the computational cost, in particular the approximation of the higher derivatives. This approach has been also used to solve both transient and steady state problems. However, as is the case of the original Kansa’s approach, none of the above modifications are completely mathematically robust and rigorous proof of their existence of solution is not available. More recently Fasshauer [26] suggested an alternative approach based on Hermite interpolation technique using radial basis functions, which allows not only the interpolation of a given function but also its derivatives. The convergence proof for RBF Hermite-Brikhoff interpolation was given by Wu [27] who also recently proved the convergence of this approach when solving PDEs (see Wu [28] and Schaback and Franke [29]). From a series of simple steady state numerical examples Fasshauer [26] concluded that those methods base on Hermite interpolation performs slightly better than those based on direct interpolation.

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370

361

Jumarhon et al. [30] and Leitao [31] observed similar improvement when using Hermite algorithms instead of direct approaches. More recently Power and Barraco [32] found that direct interpolation methods has some difficulties solving convection–diffusion problems at high Pe´clet number, which do not come across using the Hermite approach. Li and Chen [33] pointed out that these inconveniences that the direct approach has to predict high Pe´clet can be improve by using higher-order radial basis functions and overlapping domain decomposition technique. Furthermore, it is important to mention that most of the works reported in the literature based on the Hermite approach are used to solve steady state non-homogeneous problems with constant coefficients. In the present work, we report how this approach (Hermit collocation) can be implemented to find the numerical solution of time dependent convection–diffusion problems with variable coefficients. In those cases in which the fundamental solutions are given in terms of RBFs, the implementation of the MFS and BKM to solve Dirichlet boundary value problems results in a symmetric matrix system. However, when dealing with mixed boundary value problem the resulting matrix system is unsymmetric. Recently Li and Chen [34] proposed a modification of the BKM, which results in a symmetric system, using a Hermite interpolation in terms of the boundary differential operators applied to the interior fundamental solutions. This work is part of the research carried out by a consortium of European Institutions that are investigating the possibility of the sealing of the subsoil to block seawater intrusion by inducing the formation of crystals inside the porous matrix (CRYSTECHSALIN project, EVK1-CT2000-00055). The use of a Hermitian meshless collocation numerical approach was selected in the project due to its flexibility on dealing with moving boundary problems and their high precision on predicting surface fluxes, which is a crucial part of the diffusion controlled crystallization process that needs to be analysed as part of the CRYSTECHSALIN project.

where C(x,t) is the concentration of the specie at the position x at time t, U is a bounded domain in Rd, d is the dimension of the problem, D the diffusion coefficient, k the reaction or adsorption coefficient and ui the i component of the velocity field. For simplicity we have considered that the Diffusion and reaction coefficients are uniform throughout the domain of solution but the velocity field can be a function of the space. The approach presented in this article can be easily generalized for cases with variable Diffusion and reaction coefficients as well as variable velocity field, both in time and space. The parameter that describes the relative influence of the convective and diffusive components is the Peclet number, PeZuL/D. For small values of Pe Eq. (1) behave as a parabolic differential equation in time, while for large values the equation becomes more like hyperbolic. These changes in the structure of the differential equation according to the values of Pe, have significant effect on its numerical solution. In general, the above differential equation is required to satisfy the following boundary and initial conditions:

2. Convection–diffusion equation

Cðx;t CDtÞKCðx;tÞ " Dt # v2 Cðx;t CDtÞ vCðx;t CDtÞ Zq D Kuj ðxÞ KkCðx;t CDtÞ vxj vx2j " # v2 Cðx;tÞ vCðx;tÞ Cð1KqÞ D Kuj ðxÞ KkCðx;tÞ ð3Þ vxj vx2j

In this article, we present an implementation of the radial basis function Hermite collocation method for the solution of time dependent convection–diffusion problems. The equation of mass conservation of a species in a dilute solution of an incompressible flow with linear chemical reaction or absorption can be written as vCðx; tÞ v2 Cðx; tÞ vCðx; tÞ ZD K uj ðxÞ K kCðx; tÞ; vt vxj vx2j x 2U 3Rd ; tO 0

(1)

ACðx; tÞ C B

vCðx; tÞ Z f ðx; tÞ; x 2G; tO 0 vn

Cðx; tÞ Z C0 ðxÞ; t Z o

(2a)

(2b)

where A and B are known constants and f(x,t) and C0(x) are known functions. When BZ0 and As0, we refer to this type of boundary condition as the Dirichlet type, when AZ0 and Bs0 as the Neumann type and when Bs0 and As0 as the mixed (radiation) or Robin type. In Eq. (1), it is possible to approximate the time derivative of the partial differential operator by a simple forward difference using Crank–Nicholson (q weighted) method, i.e.

In this way, at each time step the transient problem can be seen as a steady state non-homogeneous problem, with the non-homogenous term as function of the solution at

362

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370

the previous time step. " !  # v2 v 1 Cðx;tCDtÞ K kDtC q Dt D 2 Kuj ðxÞ vxj q vxj " !  # v2 v 1 ZðqK1Þ Dt D 2 Kuj ðxÞ Cðx;tÞ K kDtC vxj ðqK1Þ vxj (4)

It is interesting to observe that Duchon’s thin plate splines function with mZ2 corresponds to the fundamental solution commonly used in the BEM technique to solve biharmonic problems. Chen and Tanaka [10] proposed a possible scheme to create RBFs related to other partial differential operators based upon the second Green identity. In a typical interpolation problem we have N pairs of data points fðxj ; Fðxj ÞÞNjZ1 g, which are assumed to be samples of some unknown function F that is to be interpolated by the function f, i.e.

3. Radial basis function meshless Hermite approach f ðxÞ ¼ In recent years, the theory of radial basis functions (RBFs) has undergone intensive research and enjoyed considerable success as a technique for interpolating multivariable data and functions. A radial basis function, JðxK xj ÞZ jðjjxK xj jjÞ depends upon the separation distances of a sub set of data centres, X 3
2mK2

log rðgeneralized thin plate spineÞ;

ðgeneralized multiquadricÞ and e

b$r

2

2 m=2

ðr þ c Þ

ðGaussianÞ

where m is an integer number and rZjjxKxjjj. The Gaussian and the inverse multiquadric i.e. m!0 in the generalized multiquadric function, are positive definite functions, while the thin-plate splines and the multiquadric i.e. mO0 in the generalized multiquadric function, are conditionally positive definite functions of order m, which require the addition of a polynomial term P of order mK1 together with a given homogeneous constraint conditions (see Eq. (7)) in order to obtain an invertible interpolation matrix. The multiquadric functions with values of mZ1 and cZ0 are often referred to as conical functions whilst with mZ3 and cZ0 as Duchon cubics. Duchon [35] derived the thin plate splines (TPS) as an optimum solution to the interpolation problem in a certain Hilbert space via the construction of a reproducing kernel. Therefore, they are the natural generalisation of cubic splines in nO1 dimension. Even though TPS have been considered optimal in interpolating multivariate functions they do, however, only converge linearly (see Powell [36]). On the other hand, the Hardy multiquadratics (MQ) functions converge exponentially and always produce a minimal semi-norm error as proved by Madych and Nelson [37]. However, despite MQ’s excellent performance, it contains a free parameter, c2, often referred to as the shape parameter. When c is small the resulting interpolating surface is pulled tightly to the data points, forming a conelike basis functions. As c increases, the peak of the cone gradually flatters.

N X

lj Jðjjx K xj jjÞ þ PmK1 ðxÞ x 2<2

(5)

j¼1

in the sense that Fðxi Þ ¼

N X

lj Jðjjxi K xj jjÞ þ PmK1 ðxÞ

(6)

j¼1

along with the constrains condition N X

lj Pk ðxj Þ ¼ 0 1% k% m K 1

(7)

j¼1

Here the numbers lj, jZ1,2,.,N, are real coefficients and J is a radial basis function. The matrix formulation of the above interpolation problem can be written as AxZb with ! J PmK1 A¼ (8) PTmK1 0 xTZ(l, b) and bTZ(F, 0), where b are the coefficients of the polynomial. Micchelli [38] proved that for a case when the nodal points are all distinct, the matrix resulting from the above radial basis function interpolation is always non singular. Although theoretically the resulting matrix from the above interpolation technique is always invertible, numerical experiments show that the condition number of the matrix obtained with the use of smooth RBFs like Gaussian or multiquadrics are extremely large when compared with those resulting from non-smooth RBFs like the thin-plate splines for low values of m (see Schaback [39]). Similar difficulties to those encountered with smooth functions are found when using non-smooth functions with large values of m. Let us now consider a boundary value problem L½uðxÞ Z f ðxÞ

(9a)

B½uðxÞ Z gðxÞ

(9b)

where the operators L and B are linear partial differential operators on the domain U and at the contour G, respectively. For the above problem, a Hermite RBF collocation method represents the solution u by an

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370

interpolation function of the following type n X

uðxÞ ¼

lk B
x Jðjjx K xk jjÞ þ

k¼1

N X

363

Furthermore to avoid singularity at rZ0 on the resulting differential operators of the matrix A, we use in the representation formula (10) the generalized TPS

lk L
x Jðjjx

j Z r 6 log r

k¼nþ1

K xk jjÞ þ PmK1 ðxÞ

(10)

(13)

together with the corresponding cubic polynomial. PðxÞ Z lNC1 x31 C lNC2 x32 C lNC3 x21 x2 C lNC4 x1 x22

with n as the number of nodes on the boundary of U, NKn the number of internal nodes, L
x and B
x are the adjoint differential operators of those used in (9a,b), but acting on J viewed as a function of the second argument x. This expansion for u leads to a collocation matrix A, which is of the form 0

Bx B
x ½J Bx L
x ½J Bx PmK1

B
A¼B @ Lx Bx ½J Bx PTmK1

Lx L
x ½J Lx PTmK1

1

C Lx PmK1 C A

(11)

0

The matrix (11) is of the same type as the scattered Hermite interpolation matrices and thus non-singular as long as J is chosen appropriately. The convergence of the above approach has been proven by Schaback and Franke [29] in terms of a generalized Fourier transform analysis (see also Wu [28]). Due to the uncertainty regarding which RBF is the best to use in a collocation method for the solution of boundary value problems for partial differential equations in this work we will use the TPS. A possible alternative is to use the kernel RBF for the non-transient convection–diffusion equation with constant velocity suggested by Chen and Tanaka [10] scheme. In the matrix representation (11) of the Hermite collocation numerical solution of Eq. (4), in the cases when boundary conditions of the first, second and mixed kind (Dirichlet, Neumann and Robin) are prescribed, the partial differential operators are given by the following expressions: " !  # v2 v 1 Lx Z q Dt D 2 K uj ðxÞ ; K k Dt C vxj q vxj "

v2 v L
x Z q Dt D 2 C uj ðxÞ vxj vxj N BD x Z 1; Bx Z

!



1 K k Dt C q

#

;

v v nj ðxÞ; BRx Z AðxÞ C DðxÞ n ðxÞ; vxj vxj j

v v
N n ðxÞ; B
Rx Z AðxÞ K DðxÞ n ðxÞ B
D x Z 1; Bx Z K vxj j vxj j (12) In the above relations, the super index D, N, R in the operator B represent the type of boundary conditions implemented, i.e. Dirichlet, Neumann and Robin.

C lNC5 x21 C lNC6 x22 C lNC7 x1 x2 C lNC8 x1 C lNC9 x2 C lNC10

(14)

The non-homogeneous term is obtained by the multiplication ~ ½L~ x B
x ½J; L~ x L
x ½J; L~ x Pm , of the rectangular matrix AZ defined at the internal points, by the l coefficient found at the previous time step, where " !  # 2 v v 1 L~ x Z ðq K1Þ Dt D 2 Kuj ðxÞ K k Dt C vxj ðq K1Þ vxj

4. Numerical examples For the validation of the proposed numerical approach, four standard test problems with known analytical solution are considered. In each case, comparison between the corresponding analytical and numerical results is presented. 4.1. Heat equation (pure diffusion) Let us first consider case of the linear flow of heat in a solid bounded by a pair of parallel planes. Basically such solid, characterized by its own thermal and physics proprieties, is forced to respect initial and boundary conditions and the time evolutions of the temperature and heat flow is analysed. The differential equation to be solved is: vT v2 T Zk 2 vt vxj

(15)

where T is the temperature across the solid and k is a constant conductivity coefficient. Let us consider a rectangular region (0!x!2l) and (0!y!h), with zero initial temperature and with the ends surfaces at xZ0 and xZ2l kept at constant temperature T0, for tO0, and zero flux at the lateral walls. The analytical solution of the above problem is provided by Carslaw and Jaeger [40] in terms of the dimensionless variables T/T0, x/l and kt/l2. Due to the symmetricity of the problem, it is possible to analyse only half of the solid with the original temperature T0 given at the surface xZ0 and zero flux at the centre line xZl. To show the convergence of the numerical scheme, two different uniform distributions of surface and internal nodes are analysed in a rectangular domain [10!6] with

364

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370

(a) 6

(b) 6

5

5

4

4

3

3

2

2

1

1

0 0

1

2

3

4

5

6

7

8

9

0

10

0

1

2

3

4

5

6

7

8

9

10

Fig. 1. (a) Mesh with 81 and (b) mesh with 277 collocation points.

Only a longitudinal convection term is considered, i.e.

81 and 277 collocation points, respectively (see Fig. 1). A fully implicit scheme, i.e. qZ1, is implemented with a time step DtZ0.1 s. In Fig. 2, we present the comparison between the solutions obtained with each of the two different uniform distributions of nodes, given in terms of the dimensionless variables defined before. As can be observed as the evolution time increases substantial difference between the results is found. A comparison between the result obtained with the 277 points and the analytical solution is given in Fig. 3. In this case, excellent agreement is found even for larger evolution time. No major improvement was observed by increasing the number of collocation points more than 277 nodes.

vC v2 C vC Z D 2 K u1 vt vx1 vxj

(16)

The analytical solution of the above transport problem in a one-dimensional semi-infinity domain with initial and boundary conditions Cðx; 0Þ Z Ci Z 0

Cð0; tÞ Z C0

vC ðN; tÞ Z 0 vx

(17)

is given by the following expression (see van Genuchten and Alves [41]): Cðx; tÞ Z Ci C ðC0 K Ci ÞAðx; tÞ

(18)

4.2. Convective-diffusion transport equation

where

In this section, the solution of the convection–diffusion equation with constant diffusion coefficient and convection velocity is presented. The numerical solution of this equation presents the challenging of proper representation of the valance between the convective and diffusion processes. As before in this case we compare the numerical results with the analytical solution of the problem.

The semi-infinity one-dimensional domain of the above problem will be treated numerically as a bounded one consisting of a rectangular [20!1] with zero lateral flux (see Fig. 4). Two different types of boundary conditions

    1 x K ut 1 x C ut Aðx; tÞ Z erfc C expðux=DÞerfc 2 2 2ðDtÞ1=2 2ðDtÞ1=2 (19)

1

0.02 277 points 0.04 277 points 0.1 277 points 0.02 81 points 0.04 81 points 0.1 81 points

0.9 0.8 0.7

T/T0

0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x/l Fig. 2. Comparison between the results obtained with a mesh of 81 (solid line) and 277 (dash line) collocation points for different dimensionless time t 0 Zkt/l2.

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370 1 0.9

0.01 0.02 0.03 0.04 0.06 0.08 0.1 0.2 0.3 0.4 0.6 0.8 1 Analytical

0.8 0.7

T/T0

365

0.6 0.5 0.4 0.3 0.2 0.1 0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x/l Fig. 3. Comparison between the analytical (solid line) and numerical (points) solutions for the case of 277 collocation points for different dimensionless time t 0 Zkt/l2.

(Dirichlet and Neumann) at the fictitious far field crosssection (xZ20) will be considered. As a first example, we solve the above problem imposing a zero the flux at xZ20 instead of the condition vCðN; tÞ=vxZ 0. In this way, we are able to analyse the process until the concentration front start to reach the end cross-section at xZ20. In the numerical simulation with a fully implicit scheme, a diffusion coefficient DZ1 m2/s, a velocity uZ0.6 m/s and a time step DtZ0.025 s were used. Comparison between the above analytical solution and our numerical result is given in Fig. 5, showing excellent agreement. As a second case, instead of defining zero flux at the end cross-section, we impose zero concentration, i.e. C(20,t)Z 0. In this case, we use the same diffusion coefficient as well as the same mesh than before, but considered a convective velocity of uZ0.1 m/s. Comparison between numerical and analytical solution, for this case, is shown in Fig. 6. As can be observed excellent results were obtained in both cases, i.e. given flux or concentration at the fictitious crosssection. As before, in this case the solution was stopped before the concentration front reaches the auxiliary contour.

4.3. Convection–diffusion and reaction As third example, a convection–diffusion equation with a linear reaction term is solved. As in the previous example, the convective term is considered only the longitudinal direction and all the coefficients in the governing equation are constants, i.e. vC v2 C vC Z D 2 K u1 K kC vt vx1 vxj

(20)

In this case, we consider the initial boundary value problem in a rectangular domain [0.7!6.0 m], defined by the initial condition CZ0 and boundary condition CZ300 kg/m3 at xZ0 and vC/vnZ0 at xZ6 m. The physical parameters in the governing equation were assumed constant, with diffusivity DZ1 m2/s, a velocity uZ6 m/s and a reaction coefficient kZ0.278 sK1. The computational domain was discretized with a mesh of 385 collocation points and a fully implicit scheme was used with a time step value of 10K3 s. The obtained result of the time evolution of the concentration profile along the centre line of the domain is compared in Fig. 7 with the analytical solution given in

6.0 5.0 4.0 3.0 2.0 1.0 0.0 0

2

4

6

8

10

12

14

Fig. 4. Schematic representation of the computational domain.

16

18

20

366

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370 1 0.9 0.8 1s 2s 3s 4s 5s 6s 7s 8s 9s 10s Analytical

0.7

C/C 0

0.6 0.5 0.4 0.3 0.2 0.1 0 0

2

4

6

8

10

12

14

16

18

20

x (m) Fig. 5. Comparison between the analytical solution (solid line) and the numerical result (points), for the case of constant convective velocity and diffusion coefficient with Neumann boundary condition at xZL.

1 0.9 0.8

1s 2s

0.7

3s 4s

0.6

C/C 0

5s 6s

0.5

7s 8s

0.4

9s 10s

0.3

Analytical 0.2 0.1 0 0

2

4

6

8

10

12

14

16

18

20

x (m) Fig. 6. Comparison between the analytical solution (solid line) and the numerical result (points), for the case of constant convective velocity and diffusion coefficient with Dirichlet boundary condition at xZL.

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370

367

300

250

Concentration (kg/m3)

200

0.1s 0.25s 0.5s 1s 2s Analytical

150

100

50

0 0

1

2

3

4

5

6

x (m) Fig. 7. Comparison between the analytical (solid line) and numerical (points) solutions, for the case of constant convective velocity, diffusion coefficient and reaction term.

Partridge et al. [20], from which we observe excellent agreement between them. To validate our scheme in the case of a two-dimensional convective velocity field, the above problem is solved in the same rectangular domain but rotated by the inclination angle qZp/8, as showed in Fig. 8. The problem is solved in the original (x1, x2) coordinate system using the same initial and boundary conditions but now defined in the inclined geometry. In this way along the problem domain, the velocity filed is constituted by two components u1Zu$cos q and u2Zu$sin q. When presenting the obtained new results along the centre line of the domain, i.e. along an inclined line at an angle qZp/8, we obtain exactly the same results as those reported in Fig. 7.

at all time at infinity. In cylindrical coordinates system, this transient problem is described by the following equation:   v 1 vT vT (21) Z k vr r vr vt The analytical solution of the above problem is provided by Carslaw and Jaeger [40] in terms of the dimensionless variables T/T0, log(r/a) and kt/a2. By expanding the cylindrical Laplacian operator in the above equation, we obtain the following expression k

v2 T k vT vT Z C 2 r vr vt vr

(22)

x2 x'2

4.4. Variable velocity Although only problems with constant velocity have been analysed in the previous examples, our original formulation of the proposed Hermite RBF numerical scheme was defined to deal with variable velocity fields and coefficients. This section discusses the application of the proposed approach to cases with variables convective velocity field. Let us consider the problem of transient heat transfer in an infinite plane plate with a circular hollow at the origin, which is kept at a constant temperature. Initially the plate is at zero temperature and it is required a bounded temperature

x2

x'1 x'1 θ

x'2

x1 x1 Fig. 8. Schematic representation of the rotated rectangular domain.

368

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370 1

0.8

0.1s 0.3s 1s 3s 10s 30s Analytical

T/T 0

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Log (r /a) Fig. 9. Comparison between the analytical (solid line) and numerical (points) solutions, with qZ1 and DtZ0.025 s, for the case of. transient heat transfer in an infinity plane plate with a circular hollow.

1

0.8

12s 20s 60s 200s 400s 800s Analytical

T/T 0

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Log (r /a) Fig. 10. Comparison between the analytical (solid line) and numerical (points) solutions, with qZ1 and DtZ1 s, for the case of. transient heat transfer in an infinity plane plate with a circular hollow.

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370

369

1

0.8

12s 20s 60s 200s 400s 800s Analytical

T/T 0

0.6

0.4

0.2

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Log (r /a) Fig. 11. Comparison between the analytical (solid line) and numerical (points) solutions, with qZ1/2 and DtZ1 s, for the case of. transient heat transfer in an infinity plane plate with a circular hollow.

which corresponds to a one-dimensional convection–diffusion equation with variable negative convective field, urZKk/r, and a constant diffusion coefficient k. In our numerical simulation, we consider the above onedimensional problem as a two-dimensional one, in a rectangular [1!40 m] with prescribed zero lateral flux. A constant temperature T0 is kept at the circular hollow, represented in our 2D solution by the cross-section at x1Za, while a zero flux condition is prescribed at the auxiliary far field surface defined at x1Zb. The initial temperature is equal to zero all over the domain. In our 2D representation of problem, the convective velocity field is defined as u1ZKk/x1, u2Z0. The analysis is performed using the following geometrical and physical parameters: aZ0.5 m, bZ20.5 m and kZ1 m2/s. Besides, a time step equal to 0.025 s and a value of qZ1 is used. The results obtained for the time evolution of the temperature profile is shown in Fig. 9, given in terms of the dimensionless variables defined before. As it can be seen they are in excellent agreement with the analytical solution. While the fully implicit, i.e. qZ1, unconditionally stable scheme used in the previous examples works perfectly with smaller time steps, a more refined procedure needs to be used for larger time evolution. As the value of Dt increases the fully implicit scheme loses its precision, since its accuracy is only of first order in time. In Fig. 10, we present the results obtained for qZ1 and a time step of DtZ1 s, where we can observe accuracy loss for long run at the region close to the auxiliary far field section. However,

when qZ1/2 (Crank–Nicolsom scheme) our numerical approach is able to reproduce the analytical solution even for time step of DtZ1 s and large time evolution of the process (see Fig. 11). This is due to the second order accuracy in time of the scheme in this case.

5. Conclusions The results present in Section 4 show the versatility of the radial basis function Hermite collocation approach to solve time dependent Convective Diffusion problems involving variable velocity and reaction coefficient. Using the standard Crank-Nicholson weighed method to represent the time derivative in the partial differential equation the original problem reduces to the solution of a steady state non-homogeneous convective diffusion problem at each time step, with the non-homogeneous term proportional to the solution of the previous time step. One of the major contributions of this article is to show the accuracy of the scheme even in cases of strong variable velocity field and considering large evolution times.

Acknowledgements This research was sponsored by the CRYSTECHSALIN project (Contract number EVK1-CT-2000-00055)—part of the FP5, Energy, Environment and Sustainable

370

A. La Rocca et al. / Engineering Analysis with Boundary Elements 29 (2005) 359–370

Development European Commission Programme. The authors also will like to thank the reviewers for their informative comments, which really improve this manuscript.

References [1] Atluri SN, Zhu T. New meshless local Petrov-Galerkin (MLPG) approach in computational mechanics. Comput Mech 1998;22(2): 117–27. [2] Atluri SN, Zhu T. New concepts in meshless methods. Int J Numer Meth Eng 2000;47:537–56. [3] Golberg MA, Chen CS. Discrete projection methods for integral equations. Southampton, UK: Computation Mechanics Publications; 1997. [4] Kupradze VD, Aleksidze MA. The method of functional equations for the approximate solution of certain boundary value problems, USSR; 1964. [5] Kitagawa T. On the numerical stability of the method of fundamental solutions applied to Dirichlet problem. Jpn J Appl Math 1988;35: 507–18. [6] Kitagawa T. Asymptotic stability of the fundamental solution method. J Comput Appl Math 1991;38:263–9. [7] Golberg MA, Chen CS. The method of fundamental solutions for potential, Helmholtz and diffusion problems. In: Colberg M, editor. Boundary integral methods numerical and mathematical aspects. Southampton: WITpress Computational Mechanics Publications; 1999. [8] Bogomolny A. Fundamental solutions method foe elliptic boundary value problems. SIAM J Numer Anal 1985;22:644–69. [9] Golberg MA, Chen CS, Karur SR. Improved multiquadric approximation for partial differential equations. Eng Anal Bound Elements 1996;18:9–17. [10] Chen W, Tanaka M. New insights in boundary-only and domain-type RBF methods. Int J Nonlinear Sci Numer Simul 2000;1(3):145–52. [11] Chen W, Tanaka M. A meshfree, integration-free, and boundary-only RBF technique. Comput Math Appl 2000;43:379–91. [12] Kansa EJ. Multiquadrics—A scattared data approximation scheme with applications to computation fluid-dynamics-I: surface approximations and partial derivatives estimates. Comput Math Appl 1990; 19:127–45. [13] Kansa EJ. Multiquadrics—A scattared data approximation scheme with applications to computation fluid-dynamics-II: solution to parabolic, hyperbolic and elliptic partial differential equations. Comput Math Appl 1990;19:147–61. [14] Dubal MR. Domain decomposition and local refinemet for multiquadric approximations. I: second-order equations in one-dimension. J Appl Sci 1994;1(1):146–71. [15] Hon YC, Mao XZ. An efficient numerical scheme for Burgers’ equations. Appl Math Comput 1998;95:37–50. [16] Zerroukat M, Power H, Chen CS. A numerical method for heat transfer problems using collocation and radial basis functions. Int J Numer Meth Eng 1998;42:1263–79. [17] Dubal MR, Olivera SR, Matzner RA. In: d Inverno R, editor. Approaches to numerical relativity. Cambridge UK: Cambridge University Press; 1993. [18] Kansa EJ, Hon YC. Circumventing the ill conditioning problem with multiquadric radial basis functions: applications to elliptic partial differential equations. Comput Math Appl 2000;39:123–37. [19] Wu YL, Shu C. Development of RBF-DQ method for derivative approximation and its application to simulate natural convection in concentric annuli. Comput Mech 2002;29:477–85.

[20] Partridge PW, Brebbia CA, Wrobel LC. The dual reciprocity boundary element method. London: Computational Mechanics Publications, Elsevier Applied Science; 1992. [21] Wu Z, Hon YC. Convergence error estimate in solving free boundary diffusion problem by radial basis functions method. Eng Anal Bound Elements 2003;27:73–9. [22] Golberg MA, Chen CS. The theory of radial basis functions applied to the BEM for inhomogeneous partial differential equations. Bound Elements Commun 1994;5:57–61. [23] Hon YC, Chen W. Boundary knot method for 2D and 3D Helmholtz and convection–diffusion problems under complicated geometry. Int J Numer Methods Eng 2003;56:1931–48. [24] Chen W, Hon YC. Numerical investigation on convergence of boundary knot method in the analysis of homogeneous Helmholtz, modified Helmholtz and convection–diffusion problems. Comput Meth Appl Mech Eng 2003;192:1859–75. [25] Mai-Duy N, tran-Cong T. Numerical solution of Navier-Stokes equations using multiquadric radial basis function networks. Int J Numer Meth Fluids 2001;37:65–86. [26] Fasshauer GE. In: Le Me´chaute´ A, Rabut C, Schumaker LL, editors. Solving partial differential equations by collocation with radial basis functions, proceedings of Chamonix, vol. 1–8. Nashville, TN: Vanderbilt University Press; 1996. [27] Wu Z. Hermite-Birkhoff interpolation of scatterd data by radial basis functions. Approx Theory 1992;8(2):1–11. [28] Wu Z. Solving PDE with radial basis function and the error estimation. In: Chen Z, Li Y, Micchelli CA, Xu Y, Dekker M, GuangZhou, editors. Advances in computational mathematics, lecture notes on pure and applied mathematics, vol. 202. [29] Schaback R, Franke C. Covergence order estimates of meshless collocation methods using radial basis functions. Adv Comput Math 1998;8(4):381–99. [30] Jumarhon B, Amini S, Chen K. The Hermite collocation method using radial basis functions. Eng Anal Bound Elements 2000;24:607–11. [31] Leitao VMA. A meshless method for Kirchoff plate bending problems. Int J Numer Methods Eng 2001;52:1107–30. [32] Power H, Barraco V. A comparison analysis between unsymmetric and symmetric radial basis function collocation methods for the numerical solution of partial differential equations. Comput Math 2002;43:551–83. [33] Li J, Chen CS. Some observations on unsymmetric radial basis function collocation methods for convection–diffusion problems. Int J Numer Methods Eng 2003;57:1085–94. [34] Li J, Chen CS. Some observations on unsymmetric radial basis function collocation methods for convection–diffusion problems. Int J Numer Meth Eng 2003;57:1085–94. [35] Duchon J. Spline minimizing rotation—invariant seminorms in Sobolev spaces. In: Schempp W, Zeller K, editors. Constractive theory of functions of several variables, lecture notes in mathematics, vol. 571. Berlin: Springer; 1977. p. 85–100. [36] Powell MJD. The uniform convergence of thin plate spline interpolation in two dimensions. Numerische Mathematik 1994; 68/1:107–28. [37] Madych WR, Nelson SA. Multivariable interpolation and conditionally positive definite functions, II. Math Comput 1990;54: 211–30. [38] Micchelli CA. Interpolation of scattered data: distance matrices and conditionally positive definite functions. Constr Approx 1986;2: 11–22. [39] Schaback R. Multivariate interpolation and approximation by translates of a basis function. In: Chui CK, Schumaker LL, editors. Approximation theory VIII, 1995. p. 1–8. [40] Carslaw HS, Jaeger JC. Conduction of heat in solids. Oxford: Oxford at the Clarendon press; 1959. [41] van Genuchten M Th, Alves WJ. Analytical solutions of the onedimensional convective–diffusion transport equation, US Department of Agriculture, Technical bulletin No. 1661, USA; 1982.