Time-dependent Fundamental Solutions for Homogeneous Diffusion Problems

(1)

Abstract

This paper describes the applications of the method of fundamental solutions (MFS) for 1-, 2- and 3-D diffusion equations. The time-dependent fundamental solutions for diffusion equations are used directly to obtain the solution as a linear combination of the fundamental solution of the diffusion operator. The proposed scheme is free from the conventionally used Laplace transform or the finite difference scheme to deal with the time derivative of the governing equation. By properly placing the field points and the source points at a given time level, the solution is advanced in time until steady state solutions are reached. Test results obtained for 1-, 2- and 3-D diffusion problems show good comparisons with the analytical solutions and some with the MFS based on the modified Helmholtz fundamental solutions, thus the demonstration present numerical scheme of MFS with the space–time unification has been demonstrated as a promising mesh-free numerical tool to solve homogeneous diffusion problem.

Keywords: Method of fundamental solutions; Diffusion equation; Diffusion fundamental solution; Multi-dimensions

1. Introduction

In the recent years, the meshless or mesh-free methods have received a considerable attention as alternative numerical schemes to the classical mesh-dependent numeri-cal methods. The boundary element method (BEM) is one of the most common mesh-reduction methods. One great advantage of this method is that the dimensionality of the problem is reduced by one order. Zhu [1], Zhu et al.[2], Bulgakov et al. [3], Zerroukat [4], Sutradhar et al. [5], and Bialecki et al.[6] applied the BEM to solve diffusion equations. However, when the governing equation contains any non-homogeneous term then the governing equation can be treated by the dual reciprocity method, which involves domain discretization as proposed by Nardini and Brebbia [7]. Most of the numerical schemes based on the boundary element method treats the time derivative term in the diffusion equation either using the Laplace transform

[1,2,5] or the finite difference scheme [3] to advance the solution in the time domain. However, the main drawback of the boundary element method is the determination of the singular integrals at the boundaries, which requires a great amount of computational effort especially for the case of 3-D problems.

The problem posed by the boundary element method can be alleviated by the use of the method of fundamental solutions. Detailed accounts of this method for various types of partial differential equations were summarized by Tsai

[8]. The method of fundamental solutions makes use of the fundamental solution to the given problem by satisfying the governing differential equation in the interior of the computational domain under consideration. By means of boundary collocation method, the boundary conditions for the problem are satisfied. This method is free from the singular integral evaluation problem as required by the boundary element method. If a fundamental solution exists for the given governing equation, then the method of fundamental solutions can be utilized to obtain numerical solutions of the variable at a definite number of points in the domain. Golberg [9] and Golberg and Chen [10] used

* Corresponding author. Fax: C886 2 2362 6114. E-mail address: [email protected] (D.L. Young).

1

(2)

the method of fundamental solutions to obtain numerical solution of the Poisson equation. Fairweather and Karageor-ghis[11]and Poullikkas et al.[12,13]solved the harmonic and biharmonic boundary value problems by the method of fundamental solutions. Chen et al. [14] and Balakrishnan and Ramachandran[15]applied the method of fundamental solutions for diffusion equations by using the modified Helmholtz fundamental solution. Walker [16] applied the diffusion fundamental solution combined with the inte-gration in time to solve the diffusion equation. An important issue of the method of fundamental solutions is the locations of source points, which was circumvented with consider-ation of the source positions as unknown variables by Fairweather et al.[11,17].

The works available for the solution of diffusion equation using the MFS utilize either the Laplace transform or the finite difference scheme to deal with the time derivative of the governing equation[14,15]. This is due to the fact that the MFS is always treated in the spatial domain with respect to the location of the source points and the field points. In this paper, the fundamental solution to the diffusion equation is directly used with the method of fundamental solutions without the need for the Laplace transform or the finite difference method to take care of the time-derivative term. Since any diffusion process is a time evolution process, which starts from a set of initial conditions, the transient parts of the diffusion process can be obtained in a number of time steps using the MFS. The proposed method is used to find numerical solution for diffusion processes in 1-, 2- and 3-D geometries. As a first attempt the method is applied for a set of problems with Dirichlet boundary conditions. Analytical solutions are chosen for comparison to test the accuracy of the present method. The paper is organized as follows. In Section 2, the governing equations and the boundary conditions considered are given. The numerical discretization of the method of fundamental solutions is discussed in Section 3. The discussion on the comparison of the present results with the analytical results for the test cases is presented in Section 4. The final conclusions drawn based on the present study are given in Section 5.

2. Governing equation

Consider a linear diffusion equation that has to be solved over a computational domain U enclosing a boundary G vuð~x; tÞ

vt Z kV

2_uð_{~x; tÞ; ~x 2U} ₍₁₎

in which~x is the general spatial coordinate, t is the time, k is the diffusion coefficient and u is the scalar variable to be determined. The boundary conditions for real diffusion situations may be more complex, varying from simple Dirichlet, Neumann to mixed boundary conditions.

The initial conditions of the diffusion problem is

uð~x; t0Þ Z u0; ~x 2U (2)

with the Dirichlet and Neumann boundary conditions uð~x; tÞ Z ur; ~x 2GD; t0%t% t0CDt (3) vuð~x; tÞ vnG Z qG; ~x 2G N_{; t} 0%t% t0CDt (4)

in which GDand GNstand for the parts of the boundary G where the Dirichlet and Neumann boundary conditions are prescribed, t0is initial time, Dt is positive time increment,

and~nGis the normal vector on the boundary GN. u0, uGand qG

are known functions. Since a new numerical procedure based on MFS is being tested for diffusion problems, a simple boundary condition such as the Dirichlet boundary condition is considered as a first attempt.

3. Numerical method

The fundamental solution of the linear diffusion equation is governed by:

vGð~x; t; ~x; tÞ vt Z kV

2_Gð_{~x; t; x; tÞ Cdð~x K~xÞdðt KtÞ} ₍₅₎

By taking Fourier transform and inverse Fourier transform of Eq. (5), the free space Green’s function or the fundamental solution of the diffusion equation can be obtained as Gð~x; t; ~x; tÞ ¼ exp Kj~xK~xj2 4k½tKt ð4pk½t K tÞn₌₂Hðt K tÞ (6)

where n is the spatial dimension number and H(t) is the Heaviside step function.

Since the diffusion fundamental solution satisfies the homogenous diffusion equation, the solution can be assumed as a linear combination of the fundamental solution of the diffusion operator. According to the method of fundamental solutions, the numerical solutions of the diffusion equation will be of the following form

uð~x; tÞ Z X

N_iCN_b

jZ1

ajGð~x; t; ~xj; tjÞ (7)

where ~x represents the location of the field points and ~x gives the location of the source points. t and t are the time of the field and source points, respectively, Niand Nbare the

number of initial and boundary source points and ajare the

undetermined coefficients. The distributions of the field points and the source points are shown inFig. 1(a) and (b) for 1- and 2-D, respectively. The consideration of the issue of unknown source positions can be found in references of Fairweather et al. [11,17]. In Fig. 1(a), the field points are placed at the boundary portion at tZ(nC1)Dt and at

(3)

the interior domain at tZ(n)Dt. The source points are placed in the same position, but at different time levels. By collocating these field points and using Eqs. (2–4,7) a linear matrix system can be formed as following

½Aijfajg Z fbig (8a)

where Aij¼ exp j~xiK ~xjj 2 4kðtiKtjÞ ! ð4pkðtiKtjÞÞn=2 ; if tiOtj 0; if ti!tj 8 > > > < > > > : (8b)

The matrix {bi} is the combination of initial and

boundary conditions. After inverting the matrix system, the coefficients {aj} can be solved; the function values

inside the time–space box can thus be obtained by Eq. (7). This procedure can be repeated until either the terminal time or steady state solution is achieved.

4. Results and discussions

The proposed numerical scheme is first validated by analyzing the errors for different time increments and different number of node points. Moreover, the error comparisons with the MFS based on the modified Helmholtz fundamental solution are also carried out in this validation example. Since, the present work is the first attempt to apply the MFS directly to treat the time derivative term without resorting either to the Laplace transform or finite difference scheme, simple linear diffusion processes are considered as validation cases. The effectiveness of the method is verified by solving 1-, 2- and 3-D diffusion problems and the results obtained are discussed in the following sections. Since the proposed numerical scheme considers time also as one of the coordinates, the source points and the field points are considered as depicted in

Fig. 1(a) for the case of 1-D diffusion case of x–t domain.

Fig. 1(b) shows the source points and field points on the x– y–t domain for 2-D problem.

4.1. Validation example

In order to validate the present numerical method, the following 2-D numerical example in [0,a]![0,b], aZbZ1 is carried out. The initial and boundary conditions are shown as below: uðx;y;tZ0ÞZcos px 2 Ccos py 2 Csin px 2 Csin py 2 (9a)

uð0;y;tÞ¼ 1þcos py 2 þsin py 2 eKkp2t=4 uð1;y;tÞ¼ cos py

2 þ1þsin py 2 eKkp2t=4 uðx;0;tÞ¼ cos px 2 þ1þsin px 2 eKkp2t=4 uðx;1;tÞ¼ cos px 2 þsin px 2 þ1 eKkp2t=4 8 > > > > > > > > < > > > > > > > > : (9b)

And, the analytical solution is: u ¼ cos px 2 þ cos py 2 þsin px 2 þsin py 2 eKðp2=4Þkt (10) Moreover, k is equal to one in all the above equations. The distributions of u along the horizontal x-direction at yZ0.5 and at different time levels are shown in

Fig. 2(a)–(f). An overall look at these figures indicates the expected diffusion evolution of the scalar quantity within the computational domain. Moreover,Fig. 3shows the time evolution history of the solution at (0.5, 0.5) and the maximum absolute error. Wherein, the numerical method demonstrates that better performance is obtained by the present MFS than based on the modified Helmholtz fundamental solution. This is due to the fact that the diffusion fundamental solution is a time-dependent func-tion, which is more capable to capture the transient process (two order accuracy). On the other hand,Fig. 4depicts the error histograms for different time increments and different number of points, in which more points and smaller time increments generally give better results as expected. More-over, Fig. 5 addresses that source locations are not very sensitive (Fig. 5).

(4)

4.2. 1-D diffusion problem

1-D numerical results are obtained for diffusion problem of a scalar quantity u in x-direction with certain imposed initial and boundary conditions.

The equation governing the diffusion process can be written as

vu vtZ k

v2u

vx2; 0%x%1; tR0 (11)

Fig. 3. Comparison for the validation example (left). Time evolution of u at xZyZ0.5 (right). Time evolution history of maximum absolute error (nodes, 225; DtZ0.01).

Fig. 2. Comparison of u-distribution at yZ0.5 at different time levels for the validation example (a) tZ0.04; (b) tZ0.08; (c) tZ0.12; (d) tZ0.16; (e) tZ0.2 and (f) tZ0.24 (nodes, 225; DtZ0.01).

(5)

with the initial condition and boundary conditions given as the following

uðx; 0Þ Z sinðpxÞ; 0%x%1 (12a) uð0; tÞ Z 0; uð1; tÞ Z 0; tR0 (12b) The exact solution of the diffusion problem can be obtained using the following expression

uðx; tÞ Z sinðpxÞeKkp2t (13) where kZ0.1.

It is required to compute the solution of the scalar variable u with time along the unit length of a slab subjected to Dirichlet boundary conditions. The spatial variations of the quantity u in the horizontal x-direction at different time levels are shown in Fig. 6(a)–(f). Since diffusion equation is parabolic, the solution of the variable also will exhibit the parabolic variation as seen from the above figures. The initial value of u is assumed to be sin(px) and hence the value of u decreases continuously satisfying

the enforced Dirichlet boundary conditions as depicted in the above figures. The predicted results for u are in excellent agreement with the analytical results for all the time levels as observed in the figures. The time evolution of the variable u at xZ0.5 and the maximum absolute error are also determined and shown inFig. 7. An excellent agreement of the present results with the analytical solutions demon-strates the capability of the numerical procedure.

The 1-D case was repeated with different boundary conditions. The initial and the boundary conditions were assumed as follows:

uðx; 0Þ Z cosðpxÞ; 0%x%1 (14a) uð0; tÞ Z eKkp2t; uð1; tÞ ZKeKkp2t; t%0 (14b) The exact solution of the diffusion problem is given by uðx; tÞ Z cosðpxÞeKkp2t (15) where kZ0.5.

Fig. 4. Time evolution history of maximum absolute error (left) for different time increments (nodes, 225) (right) for different numbers of node points (DtZ0.01).

(6)

Fig. 8depicts the distributions of u along the length of the slab at different time levels. Since the ends of the slab were subjected to exponential variations of u, the shapes of variations of u for this case are different from the previous case. Any diffusion process will reach a final steady state

solution satisfying the imposed boundary condition. This trend is exactly shown in Fig. 8(f) wherein the diffusion process has come to an end. Similarly,Fig. 9describes the excellent performances of the proposed numerical method by observing the time evolution history.

Fig. 6. Comparison of u-distribution at different time levels for 1-D diffusion example I (a) tZ0.2; (b) tZ0.4; (c) tZ0.6; (d) tZ0.8; (e) tZ1.0 and (f) tZ1.2 (nodes, 31; DtZ0.01).

(7)

The method is extended to study a 2-D diffusion problem in a square slab of size [0,a]![0,b], aZbZ1. The initial

and boundary conditions are shown as below:

uðx; y; t Z 0Þ Z sinðpxÞ CsinðpyÞ (16a)

Fig. 8. Comparison of u-distribution at different time levels for 1-D diffusion example II (a) tZ0.1; (b) tZ0.2; (c) tZ0.3; (d) tZ0.4; (e) tZ0.5 and (f) tZ0.6 (nodes, 31; DtZ0.01).

Fig. 9. Comparison for 1-D diffusion example II (left). Time evolution of u at xZ0.25 (right). Time evolution history of maximum absolute error (nodes, 31; DtZ0.01).

(8)

uð0; y; tÞ Z ðsinðpyÞÞeKkp2t uð1; y; tÞ Z ðsinðpyÞÞeKkp2t uðx; 0; tÞ Z ðsinðpxÞÞeKkp2t uðx; 1; tÞ Z ðsinðpxÞÞeKkp2t 8 > > > > > > < > > > > > > : (16b)

The analytical solution of the problem is given by u Z ðsinðpxÞ C sinðpyÞÞeKkp2t (17) where kZ0.1.

The distribution of u along the horizontal x-direction at yZ0.5 at different time levels are shown in

Figs. 10(a)–(f). An overall look at these figures indicates the expected diffusion evolution of the scalar quantity within the computational domain. As time increases, the gradient between the end surfaces decreases, thus approaching a steady state condition. This trend can be observed in the above figures. Especially, a comparison between Fig. 10(a) and (f) shows that the diffusion of u has taken place almost in the entire domain. A comparison of the results for the distribution of u along x-direction at yZ0.5 with the analytical results indicates

that the present method can correctly and exactly predict the time evolution characteristics of the scalar quantity u. Fig. 11 shows the time evolution of u at a point xZyZ0.5 and the maximum absolute error in the whole domain. The present results show excellent agreement with the analytical solutions as observed in the above figure. The constant u value contours obtained by the present numerical scheme and the analytical solutions at tZ1 are shown in Fig. 12 as solid and dashed lines, respectively. Both the solid and dashed lines coincide with each other, indicating the exact agreement of the present results with analytical solutions.

The proposed numerical method based on the method of fundamental solutions is extended to study the diffusion of u inside a 3-D solid with dimensions [0,a]![0,b]![0,c], aZbZcZ1. It is assumed that the solid is subjected to the following initial and boundary conditions:

uðx; y; z; t Z 0Þ Z sinðpxÞ CsinðpyÞ CsinðpzÞ (18a)

Fig. 10. Comparison of u-distribution at yZ0.5 at different time levels for 2-D diffusion (a) tZ0.2; (b) tZ0.4; (c) tZ0.6; (d) tZ0.8; (e) tZ1.0 and (f) tZ1.2 (nodes, 441; DtZ0.05).

(9)

uð0; y; z; tÞ Z ðsinðpyÞ CsinðpzÞÞeKkp2t uð1; y; z; tÞ Z ðsinðpyÞ CsinðpzÞÞeKkp2t uðx; 0; z; tÞ Z ðsinðpxÞ CsinðpzÞÞeKkp2t

uðx; 1; z; tÞ Z ðsinðpxÞ CsinðpzÞÞeKkp2t

uðx; y; 0; tÞ Z ðsinðpxÞ CsinðpyÞÞeKkp2t uðx; y; 1; tÞ Z ðsinðpxÞ CsinðpyÞÞeKkp2t 8 > > > > > > > > > > > > > < > > > > > > > > > > > > > : (18b)

The analytical solution of the problem is given by

u Z ðsinðpxÞ C sinðpyÞ C sinðpzÞÞeKkp2t (19)

where kZ0.1.

The time evolution of u at the center node inside the domain (xZyZzZ0.5) and the maximum absolute error in the whole domain are displayed in Fig. 13. Since the boundary conditions involve exponential and sinusoidal functions, the variations of u with time are as expected. A comparison of the present results with the analytical solution indicates that the present method computes excellently the variation of u with time for 3-D diffusion also. The constant u-contours on the horizontal x–y plane at zZ0.5 obtained by the present method and the analytical solutions are shown in Fig. 14(a) and (b) for tZ3 and 6, respectively. The solid lines and dashed lines indicate the results obtained by the MFS and the analytical solutions, respectively. For both the above time levels, the results obtained by the MFS are in exact agreement with the analytical solutions as indicated by the excellent coincidence of the solid and dashed lines in the above figures.

5. Conclusions

The transient diffusion problems in multi-dimensions are solved using the method of fundamental solutions directly without using any kind of time transformation or discretization. The generally adopted Laplace transform and finite difference scheme to treat the time derivative term while using the method of fundamental solutions to solve diffusion equations are not required in the proposed numerical solution procedure. The independent time variable was treated as one of the solution domains in the fundamental solutions and this space–time unification

Fig. 11. Comparison for 2-D diffusion (left). Time evolution of u at xZyZ0.5 (right). Time evolution history of maximum absolute error (nodes, 441; DtZ0.05).

Fig. 12. Comparison of constant value u-contour for 2-D diffusion at tZ1 (nodes, 441; DtZ0.05).

(10)

method has made possible to make use of the method of fundamental solutions to obtain directly the numerical solutions of diffusion equation without transformation or discretization for the time domain. By properly locating the source points and the field points at every time level, the solution were advanced in time to get results until the system reaches the steady state conditions. The numerical procedure developed in the present work was validated by comparing the time evolutions and spatial distributions of a scalar variable with the results obtained by analytical solutions for 1-, 2- and 3-D diffusion problems under Dirichlet boundary conditions. The exact

agreement of the present results with the analytical results indicates the effectiveness of the present method to solve diffusion equations directly by the method of fundamental solutions and without requiring any time transformation or time discretization methods.

Acknowledgements

The National Science Council of Taiwan is gratefully acknowledged for providing financial support to carry out the present work.

Fig. 13. Comparison for 3-D diffusion (left). Time evolution of u at xZyZzZ0.5 (right). Time evolution history of maximum absolute error (nodes, 343; DtZ0.3).

(11)

Elem 1999;23:201–9.

[5] Sutradhar A, Paulino GH, Gray LJ. Transient heat conduction in homogeneous and non-homogeneous materials by the Laplace transform Galerkin boundary element method. Eng Anal Bound Elem 2002;26:119–32.

[6] Bialecki RA, Jurgas P, Kuhn G. Dual reciprocity BEM without matrix inversion for transient heat conduction. Eng Anal Bound Elem 2002; 26:227–36.

[7] Nardini D, Brebbia CA. A new approach to free vibration analysis using boundary elements. In: Brebbia CA, editor. Boundary element methods in engineering. Berlin: Springer; 1982.

[8] Tsai CC, Meshless numerical methods and their engineering applications. PhD Thesis. Department of Civil Engineering, National Taiwan University, Taipei, Taiwan; 2002.

[13] Poullikkas A, Karageorghis A, Grorgiou G. The method of fundamental solutions for inhomogeneous elliptic problems. Comput Mech 1998;22:100–7.

[14] Chen CS, Golberg MA, Hon YC. The method of fundamental solutions and quasi-Monte-Carlo method for diffusion equations. Int J Numer Meth Eng 1998;43:1421–35.

[15] Balakrishnan K, Ramachandran PA. The method of fundamental solutions for linear diffusion-reaction equations. Math Comput Model 2000;31(2/3):221–37.

[16] Walker SP. Diffusion problems using transient discrete source superposition. Int J Numer Meth Eng 1992;35:165–78.

[17] Fairweather G, Karageorghis A, Martin PA. The method of fundamental solutions for scattering and radiation problems. Eng Anal Bound Elem 2003;27:759–69.