The unsteady anisotropic-diffusion convection equation of variable coefficients

(1.1)∂∂xidijx∂cx,t∂xj−∂∂xivixcx,t = αx∂cx,t∂t

will be considered. We assume that the flow is compressible, that is

∂vix∂xi≠0.

Equation (1.1) is used to model unsteady diffusion convection process in anisotropic and inhomogeneous (functionally graded) materials. Among the physical phenomena of applications include pollutant transport and heat transfer. In (1.1) the coefficients dijx,vix,αx may represent respectively the diffusivity or conductivity, the velocity of flow existing in the system and the change rate of the unknown variable cx,t.

Nowadays functionally graded materials (FGMs) have become an important issue, and numerous studies on this issue for a variety of applications have been reported. Authors commonly define an FGM as an inhomogeneous material having a specific property such as thermal conductivity, hardness, toughness, ductility, corrosion resistance, etc. that changes spatially in a continuous fashion. Therefore equation (1.1) is relevant for FGMs.

In the last decade investigations on the diffusion-convection equation had been done for finding its numerical solutions. The investigations can be classified according to the anisotropy and inhomogeneity of the media under consideration. For example, Wu et al. [38], Hernandez-Martinez et al. [15], Wang et al. [37] and FendoĞlu et al. [9] had been working on problems of isotropic diffusion and homogeneous media, Yoshida and Nagaoka [39], Meenal and Eldho [22], Azis [5] (for Helmholtz type governing equation) studied problems of anisotropic diffusion but homogeneous media. Rap et al. [26], Ravnik and Škerget [28, 29], Li et al. [21] and Pettres and Lacerda [25] considered the case of isotropic diffusion and variable coefficients (inhomogeneous media). Recently Azis and co-workers had been working on steady state problems of anisotropic inhomogeneous media for several types of governing equations, for examples [7, 36] for the modified Helmholtz equation, [6, 33, 30, 13, 20] for the diffusion convection reaction equation, [32, 10, 16, 19] for the Laplace type equation, [12, 4, 23, 24, 18] for the Helmholtz equation.

Equation (1.1) applies for unsteady problems of anisotropic and inhomogeneous cases, therefore provides a wider class of problems. It covers problems of isotropic and homogeneous media as special cases which occur respectively when d11=d22,d12=0 and the coefficients dij,vi and α are constant.

Zoppou and Knight [40] had been working on finding the analytical solution to the unsteady orthotropic diffusion-convection equation with spatially variable coefficients. The equation considered is almost similar to equation (1.1) but with limitation d11≠d22,d12=0. This paper is intended to extend the recently published works on anisotropic diffusion convection equation with variable coefficients [35, 8, 11, 31, 2] from the steady state to unsteady state equation.

Referred to the Cartesian frame Ox1x2 we will consider initial boundary value problems governed by (1.1) where Ox1x2x=x1,x2. The coefficient [dij](i,j=1,2) is a real positive definite symmetrical matrix. Also, in (1.1) the summation convention for repeated indices holds, so that explicitly (1.1) can be respectively written as

∂∂x1d11∂c∂x1+∂∂x1d12∂c∂x2+∂∂x2d12∂c∂x1+∂∂x2d22∂c∂x2−∂v1c∂x1−∂v2c∂x2=α∂c∂t.

Given the coefficients dijx,vix,αx, the solution cx,t and its derivatives for (1.1) are sought which are valid for time interval t≥0 and in a two-dimensional region Ω with boundary ∂Ω which consists of a finite number of piecewise smooth curves. On ∂Ω1 the dependent variable cx,t is specified, and

(2.1)Px,t=dijx∂cx,t∂xinj

is specified on ∂Ω2 where ∂Ω=∂Ω1∪∂Ω2 and n=(n1,n2) denotes the outward pointing normal to ∂Ω. The initial condition is taken to be

(2.2)cx,0=0

The method of solution will be to transform the variable coefficient equation (1.1) to a constant coefficient equation, and then taking a Laplace transform of the constant coefficient equation, and to obtain a boundary integral equation in the Laplace transform variable s. The boundary integral equation is then solved using a standard boundary element method (BEM). An inverse Laplace transform is taken to obtain the solution c and its derivatives for all x,t in the domain. The inverse Laplace transform is implemented numerically using the Stehfest formula. The analysis is specially relevant to an anisotropic medium but it equally applies to isotropic media. For isotropy, the coefficients in (1.1) take the form d11=d22 and d₁₂=0 and use of these equations in the following analysis immediately yields the corresponding results for an isotropic medium.

We restrict the coefficients dij,vi,α to be of the form

(3.1)dijx =  d^ijgx
(3.2)vix=  v^igx
(3.3)αx= α^gx

where gx is a differentiable function and d^ij,v^i,α^ are constants. Further we assume that the coefficients dijx, vix and αx are exponentially graded by taking gx as an exponential function

(3.4)gx=expβ0+βixi2

where β₀ and β_i are constants. Therefore if

(3.5)d^ijβiβj+v^iβi−λ=0

where λ is a constant then (3.4) satisfies

(3.6)d^ij∂2g1/2∂xi∂xj+v^i∂g1/2∂xi−λg1/2=0.

Substitution of (3.4), (3.2) and (3.3) into (1.1) gives

(3.7)d^ij∂∂xig∂c∂xj −v^i∂gc∂xi=α^g∂c∂t

Assume

(3.8)cx,t=g−1/2xψx,t

therefore substitution of (3.4) and (3.8) into (2.1) gives

(3.9)Px,t=−Pgxψx,t+g1/2xPψx,t

where

Pgx,t=d^ij∂g1/2x∂xjni Pψx,t=d^ij∂ψx,t∂xjni

And equation (3.7) can be written as

d^ij∂∂xig∂g−1/2ψ∂xj −v^i∂g1/2ψ∂xi=α^g∂g−1/2ψ∂t
d^ij∂∂xigg−1/2∂ψ∂xj +ψ∂g−1/2∂xj   −v^ig1/2∂ψ∂xi+ψ∂g1/2∂xi=α^gg−1/2∂ψ∂t
d^ij∂∂xig1/2∂ψ∂xj+gψ∂g−1/2∂xj  −v^ig1/2∂ψ∂xi+ψ∂g1/2∂xi=α^g1/2∂ψ∂t

Use of the identity

∂g−1/2∂xi=−g−1∂g1/2∂xi

implies

d^ij∂∂xig1/2∂ψ∂xj −ψ∂g1/2∂xj −v^ig1/2∂ψ∂xi +ψ∂g1/2∂xi =α^g1/2∂ψ∂t
g1/2d^ij∂2ψ∂xi∂xj−v^i∂ψ∂xj−ψd^ij∂2g1/2∂xi∂xj+v^i∂g1/2∂xi+d^ij∂ψ∂xj∂g1/2∂xi−d^ij∂ψ∂xj∂g1/2∂xi=α^g1/2∂ψ∂t

Equation (3.6) then implies

(3.10)d^ij∂2ψ∂xi∂xj−v^i∂ψ∂xi−λψ=α^∂ψ∂t

which is a constant coefficients equation. Taking a Laplace transform of (3.8), (3.10), (3.10) with respect to time t, and applying the initial condition (3.10) we obtain

(3.11)ψ*x,s=g1/2xc*x,s
(3.12)Pψ*x,s=P*x,s+Pgxψ*x,sg−1/2x
(3.13)d^ij∂2ψ*∂xi∂xj−v^i∂ψ*∂xi−λ+sα ^ ψ*=0

where s is the variable of the Laplace-transformed domain and the notation * represents the quantity in the Laplace transform framework.

By using Gauss divergence theorem, equation (3.13) can be transformed into a boundary integral equation

(3.14)ηξψ*ξ,s=∫∂Ω P ψ * x,sΦ x,ξ− Pv xΦ x,ξ      +Γx,ξ ψ*x,sdSx

where

Pvx= v^inix

For 2-D problems the fundamental solutions Φ(x,ξ) and Γ(x,ξ) are given as (see for example [3] for the derivation of the fundamental solutions)

Φx,ξ=ρi2πDexp−v.R.2DK0μ˙R˙Γx,ξ=d^ij∂Φx,ξ∂xjni

where

μ˙=v˙/2D2+λ+sα^/D
D=d^11+2d^12ρr+d^22ρr2+ρi2/2
R˙=x˙−ξ˙
x˙=x1+ρrx2,ρix2
ξ.=ξ1+ρrξ2,ρiξ2
v˙=v^1+ρrv^2,ρiv^2
R˙=x1+ρrx2−ξ1−ρrξ22+ρix2−ρiξ22
v˙= v^ 1+ρr v^ 22+ρi v^ 22

where ρ_r and ρ_i are respectively the real and the positive imaginary parts of the complex root ρ of the quadratic equation

d^11+2d^12ρ+d^22ρ2=0

and K₀ is the second kind of the modified Bessel function. Use of (3.11) and (3.12) in (3.14) yields

(3.15)ηg1/2c*=∫∂Ω g −1/2ΦP *+ P g −Pv g 1/2Φ− g 1/2Γ c *dS

Equation (3.15) provides a boundary integral equation for determining the numerical solutions of c^* and its derivatives ∂c*/∂x1 and ∂c*/∂x2 at all points of Ω.

Knowing the solutions c*x,s and its derivatives ∂c*/∂x1 and ∂c*/∂x2 which are obtained from (3.15), the numerical Laplace transform inversion technique using the Stehfest formula is then employed to find the values of cx,t and its derivatives ∂c/∂x1 and ∂c/∂x2. The Stehfest formula (see [34]) is

(3.16)cx,t≃ln2t∑ m=1NVmc*x,sm∂cx,t∂x1≃ln2t∑ m=1NVm∂c*x,sm∂x1∂cx,t∂x2≃ln2t∑ m=1NVm∂c*x,sm∂x2

where

sm=ln2tmVm= −1N2+m∑ k= m+1 2 min m, N 2 kN/2 2k! N2−k!k!k−1!m−k!2k−m!

A simple script has been developed to calculate the values of the coefficients Vm,m=1,2,…,N for any number N. Table 1 shows the values of V_m for N=4,6,8,10.

Table 1 . Values of V_m of the Stehfest formula for N=4,6,8,10.

Vm	N=4	N=6	N=8	N=10
V1	-2	1	-1/3	1/12
V2	26	-49	145/3	-385/12
V3	-48	366	-906	1279
V4	24	-858	16394/3	-46871/3
V5		810	-43130/3	505465/6
V6		-270	18730	-236957.5
V7			-35840/3	1127735/3
V8			8960/3	-1020215/3
V9				164062.5
V10				-32812.5

In order to justify the analysis derived in the previous sections, we will consider several problems either as test examples of analytical solutions or problems without simple analytical solutions.

We assume each problem belongs to a system which is valid in spatial and time domains and governed by equation (1.1) and satisfying the initial condition (3.10) and some boundary conditions as defined in Section. The characteristics of the system which are represented by the coefficients dijx,vix,αx in equation (1.1) are assumed to be of the form (3.4), (3.2) and (3.3) in which gx is an exponential function of the form (3.4).

Standard BEM with constant elements is employed to obtain numerical results. And the value of N in (3.16) for the Stehfest formula is chosen to be N=10. We try to increase the value of N from N=6 to N=12 and obtain the optimized solution when using N=10. Increasing N to N=12 gives less accurate solutions. According to Hassanzadeh and Pooladi-Darvish [14] these worse results are induced by round-off errors. For a simplicity, a unit square (depicted in Figure 1) will be taken as the geometrical domain for all problems. A number of 320 boundary elements of equal length, namely 80 elements on each side of the unit square, are used. For the numerical integration, we use the 10-nodal-point Bode quadrature of error of order O(h11) where 9h is equal to the length of the boundary element (see [1]). This will guarantee the convergence of the solution as the number of boundary elements is, to some extent, increased. A FORTRAN script is developed to compute the solutions and a specific FORTRAN command is imposed to calculate the elapsed CPU time for obtaining the results.

Figure 1. The domain Ω

4.1. Test problems

Other aspects that will be justified are the accuracy and consistency (between the scattering and flow) of the numerical solutions. The analytical solutions are assumed to take a separable variables form

cx,t=g−1/2xhxft

where

hx = exp−0.5+0.1x1+0.4x2

The function g1/2x is

g1/2x = exp−0.25x1−0.15x2

and depicted in Figure 2. We will consider three forms of time variation functions

Figure 2. Function
g x

ft of time domain t=0:5 which are

ft=1−exp−1.75tft=0.15tft=0.12t5−t

We take mutual coefficients d^ij and v^i

for the problems

d^ij=10.250.250.75 v^i=0.1,0.2

so that from (3.5) we have

λ=0.043125

For hx to be a solution to (3.13), the value of α^ has to be

α^=0.016875/s

We also take a set of boundary conditions (see Figure 1)

P is given on side ABc is given on side BCP is given on side CDP is given on sideA D

Problem 1: First, we suppose that the time variation function is

ft=1−exp−1.75t

Function f(t) is depicted in Figure 3. Figure 4 shows the accuracy of the BEM solutions. The errors occur in the fourth decimal place for the c and the derivatives ∂c/∂x1 and ∂c/∂x2 solutions. Figure 5 shows the consistency between the scattering and the flow solutions which verifies that the solutions for the derivatives had also been computed correctly. Figure 6 shows that the solution c changes with time t in a similar way the function ft=1−exp−1.75t does (see Figure 3) and tends to approach a steady state solution as the time goes to infinity, as expected. The elapsed CPU time for the computation of the numerical solutions at 19×19 points inside the space domain which are

Figure 3. Function
f t for Problem 1

Figure 4. The errors of solutions c (top),
∂ c / ∂ x 1 (center),
∂ c / ∂ x 2 (bottom) at t=2.5 for Problem 1

Figure 5. Solutions c and
∂ c / ∂ x 1 , ∂ c / ∂ x 2 at t=2.5 for Problem 1

Figure 6. Solutions c for Problem 1

x1,x2=0.05,0.1,0.15,…,0.9,0.95×0.05,0.1,0.15,…,0.9,0.95

and 11 time-steps which are

t=0.05,0.5,1,1.5,…,3.5,4,4.5,5

is 7883.234375 seconds.

Problem 2: Next, we suppose that the time variation function is (see Figure 7)

Figure 7. Function
f t for Problem 2

ft=0.15t

Figure 8 shows the accuracy of the BEM solutions. The errors occur in the fourth decimal place for the c and the derivatives ∂c/∂x1 and ∂c/∂x2 solutions. Figure 9 shows the consistency between the scattering and the flow solutions. Figure 10 shows that the solution c changes with time t in a manner which is almost similar to as the function ft=0.15t does (see Figure 7), as expected. The elapsed CPU time for the computation of the numerical solutions at 19×19 spatial positions and 11 time steps from t=0.0005 to t=5 is 7912.4375 seconds.

Figure 8. The errors of solutions c (top),
∂ c / ∂ x 1 (center),
∂ c / ∂ x 2 (bottom) at t=2.5 for Problem 2

Figure 9. Solutions c and
∂ c / ∂ x 1 , ∂ c / ∂ x 2 at t=2.5 for Problem 2

Figure 10. Solutions c for Problem 2

Now, we suppose that the time variation function is (see Figure 11)

Figure 11. Function
f t for Problem 3

ft=0.12t5−t

Figure 12 shows the accuracy of the BEM solutions. The errors occur in the fourth decimal place for the c and the derivatives ∂c/∂x1 and ∂c/∂x2 solutions. Figure 13 shows the consistency between the scattering and the flow solutions which again verifies that the solutions for the derivatives had also been computed correctly. Figure 14 shows that the solution c changes with time t in a similar way the function ft=0.12t5−t does. The elapsed CPU time for the computation of the numerical solutions at 19×19 spatial positions and 11 time steps from t=0.0005 to t=5 is 6800 seconds.

Figure 12. The errors of solutions c (top),
∂ c / ∂ x 1 (center),
∂ c / ∂ x 2 (bottom) at t=2.5 for Problem 3

Figure 13. Solutions c and
∂ c / ∂ x 1 , ∂ c / ∂ x 2 at t=2.5 for Problem 3

Figure 14. Solutions c for Problem 3

4.2. Examples without analytical solutions

Furthermore, we will justify the numerical solutions and show the impact of the anisotropy and the inhomogeneity of the material under consideration on the solutions. We choose

v^i=0.1,0.2 α^=1

Problem 4: For this problem the medium is supposed to be inhomogeneous or homogeneous, anisotropic or isotropic with grading function gx, constant coefficients d^ij and corresponding λ satisfying (3.5) and (3.6) as respectively follows:

• inhomogeneous and anisotropic case
  g1/2x=exp−0.25x1−0.15x2d^ij= 1 0.25 0.25 0.75 λ=0.043125

• inhomogeneous and isotropic case
  g1/2x=exp−0.25x1−0.15x2d^ij= 1 0 0 1 λ=0.085

• homogeneous and isotropic case
  g1/2x=1d^ij= 1 0 0 1 λ=0

• homogeneous and anisotropic case
  g1/2x=1d^ij= 1 0.25 0.25 0.75 λ=0

The boundary conditions are that (see Figure 1)

P=0 on side ABc=0 on side BCP=0 on side CDP=1 on side AD

There is no simple analytical solution for the problem. In fact the system is geometrically symmetric about the axis x₂=0.5. And this had been justified by the results in Figure 15 in which it is observed that anisotropy and inhomogeneity give impact to the values of solution c for being asymmetric about x₂=0.5. Solutions are symmetric only for homogeneous isotropic case, as expected. As also expected, the results (see Figure 16) show that inhomogeneity and anisotropy effects on the values of solution c. Moreover, for all cases the results in Figure 16 indicate that the system has a steady state solution.

Figure 15. Symmetry of solutions c about
x 2 = 0.5 for Problem 4

Figure 16. Solutions c at
x 1 , x 2 = 0.5 , 0.5 for Problem 4

After all, the results show that the anisotropy and inhomogeneity of material effect the values of solution c. This suggests to take both aspects into account in experimental studies.

Problem 5: We consider the inhomogeneous and anisotropic case of Problem 4 again. But we change slightly the set of the boundary conditions of Problem 4 especially on the side AD. Now we use three cases of the boundary condition on the side AD, namely

P=1−exp−1.75t on side ADP=0.15t on side ADP=0.12t5−t on side AD

The results in Figure 17 are expected. The trends of the solutions c mimics the trends of the exponential function 1−exp−1.75t, the linear function 0.15t and the quadratic function 0.12t5−t of the boundary condition on side AD. Specifically, for the exponential function 1−exp−1.75t, as time t goes to infinity, values of this function go to 1. So for big value of t, Problem 5 is similar to Problem 4 of the anisotropic inhomogeneous case. And the two plots of solutions c for Problem 4 and Problem 5 in Figure 17 verifies this, they approach a same steady state solution as t gets bigger.

Figure 17. Solutions c at
x 1 , x 2 = 0.5 , 0.5 for Problem 5

A mixed Laplace transform and standard BEM has been used to find numerical solutions to initial boundary value problems for anisotropic exponentially graded materials which are governed by the diffusion-convection equation (1.1) of compressible flow. The method is easy to be implemented and involves a time variable free fundamental solution. It gives accurate solutions as it does not involve round error propagation. It solves the boundary integral equation (3.15) independently for each specific value of t at which the solution is computed. Unlikely, the methods with time variable fundamental solution may produce less accurate solutions as the fundamental solution sometimes contain time singular points and also solution for the next time step is based on the solution of the previous time step so that the round error may propagate.

As the coefficients dijx,vix,αx do depend on the spatial variable x only and on the same inhomogeneity or grading function gx, it will be of interest to extend the study in the future to the case when the coefficients depend on different grading functions varying also with the time variable t.

In order to use the boundary integral equation (3.15), the values cx,t or Px,t of the boundary conditions as stated in Section of the original system in time variable t$ have to be Laplace transformed first. This means that from the beginning when we set up a problem, we actually put a set of approximating boundary conditions. Therefore it is really important to find a very accurate technique of numerical Laplace transform inversion.