- Email: [email protected]
A numerical method for coupled filtration equations in porous media in presence of the inverse compaction factor
Applied Mathematics and Computation 155 (2004) 439–449 www.elsevier.com/locate/amc
A numerical method for coupled filtration equations in porous media in presence of the inverse compaction factor Zahir Seyidmamedov (Muradoglu) Department of Mathematics, Applied Mathematical Sciences Research Center, University of Kocaeli, Ataturk Bulvari, Izmit 41300, Kocaeli, Turkey
Abstract In this paper a mathematical model of the transport and accretion of colloidal type particles in porous media is examined. Different from previously considered ones, the model also includes an inverse compaction factor. We present a numerical method for the solution of the coupled nonlinear differential system of equations, which models the above phenomenon. Accuracy of the presented method is tested on concrete examples. Ó 2003 Elsevier Inc. All rights reserved.. Keywords: Coupled filtration equations; Inverse compaction factor; Implicit finite-difference scheme
1. Introduction In this paper we discuss some mathematical and computational problems related to filtration in porous media. The filtration model of colloid transport was formulated in [2]. In the case of a variable porosity that depends upon the volume of colloids retained by the medium, one-dimensional analogue of this model, which is referred to as HLL-model was developed in [3]. In absence of the inverse compaction factor b 2 ½0; 1Þ the scaled form of the HLL-model is as follows: ðxc0 c þ rÞt ¼ acxx cx ; E-mail address: [email protected] (Z. Seyidmamedov (Muradoglu)). 0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(03)00789-6
ð1Þ
440
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
rt ¼ cF ðrÞ;
ð2Þ
0 < x < 1; t > 0;
x ¼ 1 br:
ð3Þ
Here c ¼ cðx; tÞ (scaled on the colloid concentration scale c0 > 0, which is assumed to be small) is the concentration of the colloid particles, measured in volume of colloids in liquid phase per unit volume of the liquid–colloid suspension, and r ¼ rðx; tÞ denotes the volume of retained, or immobile, colloids per unit volume of porous medium. The function x ¼ xðx; tÞ (0 < x < x0 ) denotes the variable porosity of the medium (0 < x < x0 ), a ¼ kV =D, V > 0 the filtration (Darcy) velocity, D > 0 the dispersion constant and k > 0 the filter coefficient; x0 denotes the free porosity, i.e., the porosity of the porous domain when no particulants are present. The function F ¼ F ðrÞ is the accretion rate. We assume that F ð0Þ ¼ 1;
F ðr Þ ¼ 0;
F ðrÞ P 0;
F 0 ðrÞ < 0;
8r 2 ð0; r Þ:
ð4Þ
An existence of wave front solutions of (1) and (2) in the form of c ¼ cðzÞ, r ¼ rðzÞ, where z ¼ x vt, v > 0 is the wave front speed, was shown in [3] and the analysis of wave fronts for parabolic and hyperbolic systems to examine the interactions between dispersion, advection, and absorption is given in [1,5]. In the case of x ¼ 1;
a ¼ 0;
F ðrÞ ¼ 1 br
ð5Þ
an analytical solution of the problem (1)–(3) was obtained in [2]. Later this particular case was analyzed in [3]. Specifically, for the initial boundary-value problem c0 ct þ rt ¼ cx ; rt ¼ cð1 brÞ; x; t > 0; cðx; 0Þ ¼ rðx; 0Þ ¼ 0; x > 0; cð0; tÞ ¼ 1; t > 0; the analytical solution ebs 1 ebs ; rðx; sÞ ¼ ; b ex þ ebs 1 ex þ ebs 1 b 6¼ 0; s ¼ t c0 x; x > 0 cðx; sÞ ¼
was found. However this solution shows that on the characteristics s ¼ t c0 x between the functions cðx; sÞ and rðx; sÞ the relationship 1 cðx; sÞ þ rðx; sÞ ¼ ; b
b 6¼ 0; x > 0
needs to be satisfied, which may not be in practice.
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
441
In Section 2 we consider the mathematical model which includes the inverse compaction factor b. Numerical method for the solution of the coupled system of equations corresponding to this model, is given in Section 3. In the last section we discuss results of computational experiment and provide a comparison of the results with the particular cases obtained in cited works. The computational results show that the parameter b > 0 cannot be neglected, and has high influence to the function r ¼ rðx; tÞ.
2. Small concentration model in presence of the inverse compaction factor We assume that the colloid concentration is small: c0 1 and restrict the problem to be a bounded spatial domain problem. For this aim let the length scale L > 0 be extent of the domain so that the scaled domain will be 0 < x < 1. Then applying the singular perturbation method [3] to the system of Eqs. (1) and (2) we obtain the following nonlinear coupled system: acxx cx cF ðrÞ ¼ 0; rt ¼ cF ðrÞ; x 2 ð0; 1Þ; t > 0: The assumption that the function F ðrÞ is decreasing (see (4)) corresponds to the physical situation that the rate of accretion decreases as more particles are accreted, as might happen if colloid particles can accrete to the matrix, but not to other accreted colloid particles. The case F 0 ðrÞ > 0 corresponds to the situation when F ðrÞ should increase with increasing the volume of retained (or immobile) colloids per unit volume of porous medium. This means as more and more particles are retained, there is greater possibility that additional particles will be filtered. In this situation there does not exist an equilibrium state for the model (1)–(3) at z ¼ 1, where z ¼ x vt and v > 0 is the wave front speed. A more interesting problem occurs when both colmatage (clogging) and decoltomage occur simultaneously. This situation corresponds to the presence of the inverse compaction factor b > 0 and Eq. (2) to take the form rt ¼ cF ðrÞ br. In this case we obtain the following initial boundary-value problem for nonlinear system of coupled equations: acxx cx cF ðrÞ ¼ 0; ð6Þ rt ¼ cF ðrÞ br; x 2 ð0; 1Þ; t > 0; cð0; tÞ ¼ 1; cx ð1; tÞ ¼ 0; t > 0; ð7Þ rðx; 0Þ ¼ 0; x 2 ð0; 1Þ: In the case, when b ¼ 0, the qualitative analysis and an existence of solution to the above problem is given in [3]. We will refer to the nonlinear model (6) and (7) as the HLL(0)-model, if the parameter b is assumed to be zero, and as the HLL(b)-model, if b 6¼ 0. From the point of view of the mathematical
442
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
problems associated with a filtration model of colloid transport, it is important to study the influence of the parameter b to the solution hcðx; tÞ; rðx; tÞi and the behaviour of the solution on the rate of F ðrÞ. The analysis and numerical solution of the nonlinear model is much more difficult then the case b ¼ 0, and we will discuss here only the numerical method. The mathematical analysis of this model will be discussed in a forthcoming work.
3. Discretization and the numerical method for the nonlinear coupled problem (6) and (7) Finite difference discretization of the first equation (6) is performed by using central difference for the first derivative cx . For the discretization of the nonlinear parabolic equation of the coupled system (6) the implicit finite-difference scheme is used, due to its absolute stability [4]. Let us introduce the uniform mesh wh on [0,1]: wh ¼ fxi 2 ½0; 1 : x0 ¼ 0; xi ¼ hi; h ¼ 1=N ; i ¼ 1; N g, for the space variable x 2 ½0; 1. Denote by s > 0 the time step. Then the pair (xi ; tj ), where tj ¼ sj, j ¼ 1; M, will define the current mesh point (i; j) on the xt-plane. On the uniform mesh wh for each time layer tj the first equation (6) can be approximated as follows: acxx;ij cx;ij cij F ðrij Þ ¼ 0;
i ¼ 1; N ;
j ¼ 1; M;
ð8Þ
where cx;ij ¼ ðcðxiþ1 ; tj Þ cðxi1 ; tj ÞÞ=h is the central difference and cxx;ij is the finite difference approximation of the term cxx [4]. Then for each j ¼ 1; M the finite-dimensional analogue of the first equation (6) with the first boundary conditions (7) has the following form: 8 < Ai ci1;j Bi ðrÞci;j þ Di ciþ1;j ¼ /i ; i ¼ 2; N 1; ð9Þ c ¼ j1 c2;j þ m1 ; : 1;j cN ;j ¼ j2 cN 1;j þ m2 ; with the coefficients Ai ¼ 2a þ h; Di ¼ 2a h; Bi ðrÞ ¼ 4a þ 2h2 F ðri;j Þ; j1 ¼ 0; l1 ¼ 0 j2 ¼ 1; l2 ¼ 0:
/i ¼ 0;
ð10Þ
For the given coefficients (10) the linear system of algebraic equations (9) with the tridiagonal matrix can be easily solved by the factorization method. However the coefficient Bi ðrÞ depends on the function r ¼ rðx; tÞ. Finite difference approximation of the second nonlinear equation is performed by using the implicit finite-difference scheme ri;j sci;j F ðri;j Þ þ bsri;j ri;j1 ¼ 0;
j ¼ 2; M;
ð11Þ
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
443
due to its stability for parabolic equations [4]. For the numerical solution of the nonlinear equation (7) with the initial condition ri;1 ¼ 0;
i ¼ 1; N
ð12Þ
we will use the Newton method, due to the above properties of the function F ¼ F ðrÞ. Denoting by f ðri;j Þ the left hand side of the nonlinear equation (7), we obtain: ðkÞ
ðkþ1Þ ri;j
¼
ðkÞ ri;j
f ðri;j Þ ðkÞ
f 0 ðri;j Þ
;
or in the explicit form ðkÞ
ðkþ1Þ ri;j
¼
ðkÞ ri;j
ðmÞ
ðkÞ
ðkÞ
ðkÞ
ri;j ci;j F ðri;j Þ þ bsri;j ri;j1 ðkÞ
1 þ bsF 0 ðri;j Þ
;
j ¼ 2; M:
ð13Þ
The index m shows the coupling of systems (9) and (11) and its sense will be ðmÞ given below. For a given function ci;j we denote by ri;jðmÞ the solution of Eq. (11). ðmÞ Thus the coefficient B ¼ BðrÞ of Eq. (9) depends on the solution ri;j of the nonlinear equation (11), also, the solution ri;j depends on the solution ci;j of Eq. (9). For this reason, in addition to the iteration process (13), one needs to construct the second (outer) iteration process for the system of Eqs. (9) and (11). For the first time layer j ¼ 1 we have the initial condition (12). Taking into account F ð0Þ ¼ 1 we obtain F ðri;1 Þ ¼ 1 and hence for j ¼ 1 Eq. (9) represents the system of linear algebraic equations with Bi ¼ 4a þ 2h2 . Solving this system we can find ci;1 , i ¼ 1; N . As an initial iteration for the second time layer j ¼ 2 ð0Þ ð0Þ we assume ci;2 ¼ ci;1 . Then substituting this in (11) ci;2 ¼ ci;2 and applying the ð0Þ iteration scheme (13) we may find the solution ri;2ð0Þ ¼ ri;2 ðci;2 Þ of Eq. (11), ð0Þ corresponding to the above iteration ci;2 . Substituting the known ri;2ð0Þ in Eq. ð1Þ (9) and solving it we found ci;2 . Then the value is again used in (11) to find ð1Þ ri;2ð1Þ ¼ ri;2 ðci;2 Þ. The above process is stopped when both conditions kri;2ðk0 1Þ ri;2ðk0 Þ kC;h 6 e;
ðk 1Þ
kci;20
ðk Þ
ci;20 kC;h 6 e;
ð14Þ
hold. Here k kC;h is the discrete analogue of C-norm and e > 0 is the given accuracy. ðk Þ If conditions (14) hold, then the pair hci;2 ; ri;2 i ¼ hci;20 ; ri;2ðk0 Þ i is assumed to be the solution of the system (9) and (11), corresponding to the time layer j ¼ 2. Thus, the iteration algorithm for the nonlinear coupled system of Eqs. (9) and (11) can be summarized as follows:
444
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
(i1) For j ¼ 1 solve the system (9) with Bi ¼ 4a þ 2h2 and compute ci;1 ; ð0Þ (i2) For j P 2 apply NewtonÕs scheme (11) with given ci;j ¼ ci;j1 and compute the solution ri;jð0Þ of Eq. (11); ð1Þ (i3) Substitute ri;jð0Þ in (9) and compute ci;j ; ð1Þ (i4) Use ci;j in (11) to find ri;jð1Þ from the iteration scheme (13); (i5) Repeat the process (i3)–(i4) until the fulfilment conditions (14); (i6) For the time layer j þ 1 the above steps (i2)–(i5) need to repeated; ðk Þ
The found hci;j0 ; ri;jðk0 Þ i is assumed to be the solution of the system (9) and (11), corresponding to the time layer j.
4. Results of computational experiments The numerical results presented here can be divided into three groups: accuracy of the discrete model, the influence of the parameter b to the solution, i.e. comparison of the HLL(0)- and HLL(b)-models, and behaviour of the solution on the rate of decrease of the function F ðrÞ. To analyze the accuracy of the method consider the following test. Example 1 2
cðx; tÞ ¼ ext ðx 1Þ ; 1=ð2nþ1Þ
F ðrÞ ¼ 1 r
;
rðx; tÞ ¼ text ;
r 2 ½0; r Þ:
x 2 ð0; 1Þ; t 2 ð0; 1Þ; ð15Þ
The functions cðx; tÞ, rðx; tÞ satisfy the initial and boundary conditions (7), and the function F ðrÞ, given by (15), satisfies conditions (4). Substituting the above functions in Eq. (6) we get: acxx cx F ðrÞc ¼ g1 ðx; tÞ; rt ¼ cF ðrÞ br þ g2 ðx; tÞ; 0 < x < 1; 0 < t < 1; where n h i o 8 < g1 ðx; tÞ ¼ ext ðx 1Þ2 at2 t 1 þ ðtext Þ1=ð2nþ1Þ þ 2ðx 1Þ½2at 1 þ 2a ; n h io : g2 ðx; tÞ ¼ ext 1 þ bt þ xt ðx 1Þ2 1 ðtext Þ1=ð2nþ1Þ : This test problem was solved on the piecewise-uniform mesh with N ¼ 101 and M ¼ 51 (s ¼ 0:02) with the following given data: a ¼ 1, b ¼ 0:5, n ¼ 1. For the accuracy e ¼ 103 given by (14) the relative errors c ¼ kc ch kC;h =kckC;h and r ¼ kr rh kC;h =krkC;h iterations was found c ¼ r 1:7 102 , for the number of iterations k0 ¼ 100 (see (14)). Here r and c are the exact solutions defined on the above introduced mesh and rh and ch are the numerical solutions found by the above iteration scheme. The number of iterations in
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
445
NewtonÕs scheme was about 2–3. For the coarser time mesh with s ¼ 0:05 (i.e. M ¼ 21) the number of iterations to achieve the same accuracy was k0 ¼ 40. This shows that the algorithm is stable and has enough accuracy. To compare the HLL(0)-model, studied in [3], with the HLL(b)-model, given by the system of Eqs. (6) and (7), consider the following. Example 2 (Influence of the inverse compaction factor). Consider the system of Eqs. (6) and (7), with the test functions given in Example 1, for the values b ¼ 0 and b ¼ 0:5 of the inverse compaction factor. For the function F ðrÞ, given by (15), with n ¼ 3, the numerical solutions cðx; tÞ and rðx; tÞ are shown in Figs. 1 and 2, correspondingly. The solid lines correspond to the cases b ¼ 0 and b ¼ 0:5, for t ¼ 1, and the broken lines––to the cases b ¼ 0 and b ¼ 0:5, for t ¼ 0:5. The dots in Fig. 1 corresponds to the case t ¼ 0, where the solutions cðx; tÞ, corresponding to the above models are same. Since the presence of the parameter b in the first equation of system (6) is in the implicit form, the influence of this parameter to the solution cðx; tÞ is negligible. However, by increasing the time variable t this influence increases and at t ¼ 1 the relative deviation dc ¼ kcðx; 1Þb¼0:5 cðx; 1Þb¼0 kC;h =kcðx; tÞb¼0 kC;h is about dc ¼ 0:01, i.e. 1.1 t=0.0 t=0.5 t=1.0 n=3
1.05
1
β =0 β =0.5 β =0 β =0.5
C(x,t)
0.95
0.9
0.85
0.8
0.75
0.7
β =0 β =0.5 0
0.2
0.4
0.6
0.8
1
1.2
x Fig. 1. Increasing profiles of the colloid concentration cðx; tÞ for HLL(0)- and HLL(b)-models, with b ¼ 0:5, at times t ¼ 0, 0.5, 1; a ¼ 1, F ðrÞ ¼ 1 r1=7 .
446
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449 t=0.5 t=1.0 0.3
n=3
0.25 β =0 0.2
σ(x,t)
β =0.5
0.15 β =0 β =0.5 0.1
0.05
0
0
0.2
0.4
0.6
0.8
1
1.2
x Fig. 2. Increasing profiles of the function rðx; tÞ for the data given in Fig. 1.
1%. The significant differences ware observed in the solutions rðx; tÞb¼0:5 and rðx; tÞb¼0 (see Fig. 2). Thus, for t ¼ 0:5 the relative deviation is about 10% and for t ¼ 1:0––about 20%. The numerical results were obtained on the mesh with parameters N ¼ 101 and M ¼ 51 (s ¼ 0:02). The same results were obtained the time mesh with s ¼ 0:05 (M ¼ 21), which shows more time mesh can be used the finite-difference scheme (11). As it was noted above, many of analytic solutions developed in the literature assume the deep bed filtration hypothesis F ðrÞ ¼ 1. One interesting solution which does not make this assumption was obtained in [2] for the case F ðrÞ ¼ 1 br, which means linearization of the dependence F r. In order to study only the influence of the rate of decrease of the function F ðrÞ to the solution hcðx; tÞ; rðx; tÞi let us analyze the HLL(b)-model with different values n ¼ 0 and n ¼ 3 of the parameter n in the function F ðrÞ. Example 3. Consider the initial value problem (6) and (7) for the values n ¼ 0 and n ¼ 3 of the parameter n in the function F ðrÞ. The results of computational experiments for the different time layers t ¼ 0, 0.5, 1.0 are given in Figs. 3 and 4. The figures show that the values of the solution hcðx; tÞ; rðx; tÞi for the linear function F ðrÞ ¼ 1 r (n ¼ 0) (broken lines) and for the nonlinear
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449 1.1
447
β=0 , n=3 β=0 , n=0
1.05
1
C(x,t)
0.95 t=1 t=0.5 0.9
t=1
0.85
0.8
t=0.5
0.75 t=0 0.7
0
0.2
0.4
0.6
0.8
1
1.2
x Fig. 3. The influence of the rate of decrease of the function F ðrÞ to the colloid concentration cðx; tÞ at times t ¼ 0, 0.5, 1.
function (15), with n ¼ 3, are different essentially. For the colloid concentration cðx; tÞ this distinction monotonically increases by increasing the space variable x. It is also interesting to observe that for the nonlinear function F ðrÞ ¼ 1 r1=7 , independently the function cðx; tÞ, the solution rðx; tÞ is a linear function at each time layer (Fig. 4, solid lines). The obtained results show that the influence of the rate of decrease of the function F ðrÞ on the both functions cðx; tÞ and rðx; tÞ is essential. Note that in all computational experiments the value of the parameter a was assumed to be a ¼ 1. Finally, in order to study both the influences of the inverse compaction factor b and the rate of decrease of the function F ðrÞ to the solution hcðx; tÞ; rðx; tÞi we consider the following. Example 4. Consider problem (6) and (7), for the values n ¼ 1 and n ¼ 3 of the parameter n in the function F ðrÞ, and for the value b ¼ 0:5 of the inverse compaction factor. Fig. 5 shows surfaces of the functions cðx; tÞ and rðx; tÞ.
448
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449 β=0 , n=0 β=0 , n=3 0.6 t=1
σ(x,t)
0.5
0.4
t=0.5
0.3
t=1 0.2 t=0.5 0.1
0
0
0.2
0.4
0.6
0.8
1
1.2
x Fig. 4. The influence of the rate of decrease of the function F ðrÞ to the function rðx; tÞ at times t ¼ 0, 0.5, 1.
Fig. 5. Solutions hcðx; tÞ; rðx; tÞi of problem (6) and (7) for b ¼ 0:5. (a1)–(a2) and (b1)–(b2) correspond to the nonlinear functions F ðrÞ ¼ 1 r1=3 and F ðrÞ ¼ 1 r1=7 , respectively.
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
449
5. Conclusions We investigated the mathematical model of the transport and accretion of colloid type particles in porous media. Different from previously considered filtration models, the model involve the case when both colmatage (clogging) and decolmatage occur simultaneously, i.e. rt ¼ cF ðrÞ br, where b is the inverse compaction factor. In addition, the accretion rate, i.e. the function F ðrÞ, is assumed to be nonlinear. We present a numerical method for solving the corresponding nonlinear coupled system of differential equations. The numerical results permit to estimate an influence of the inverse compaction factor and nonlinearity of the accretion rate to the solution hcðx; tÞ; rðx; tÞi.
Acknowledgements The author would like to thank Professor A. Hasanoglu (Hasanov) and Dr. S. Cohn for useful suggestions.
References [1] J.L. van Duijin, P. Knabner, Solute transport in a porous media with equilibrium and nonequilibrium multiple-site adsorption: travelling waves, J. Reine Angew. Math. 415 (1991) 1–49. [2] J.P. Herzig, D.M. Leclerc, P. Le Goff, Flow of suspensions through porous media––application to deep filtration, Ind. Eng. Chem. 62 (5) (1970) 8–35. [3] D. Logan, S. Cohn, G. Ledder, On filtration equations in porous media, Commun. Appl. Nonlinear Anal. Mech. 6 (1999) 34–45. [4] A.A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, 2001. [5] A.I. Volpert, V.A. Volpert, Travelling Wave Solutions of Parabolic Systems, American Mathematical Society, Providence, 1994.
A numerical method for coupled filtration equations in porous media in presence of the inverse compaction factor Zahir Seyidmamedov (Muradoglu) Department of Mathematics, Applied Mathematical Sciences Research Center, University of Kocaeli, Ataturk Bulvari, Izmit 41300, Kocaeli, Turkey
Abstract In this paper a mathematical model of the transport and accretion of colloidal type particles in porous media is examined. Different from previously considered ones, the model also includes an inverse compaction factor. We present a numerical method for the solution of the coupled nonlinear differential system of equations, which models the above phenomenon. Accuracy of the presented method is tested on concrete examples. Ó 2003 Elsevier Inc. All rights reserved.. Keywords: Coupled filtration equations; Inverse compaction factor; Implicit finite-difference scheme
1. Introduction In this paper we discuss some mathematical and computational problems related to filtration in porous media. The filtration model of colloid transport was formulated in [2]. In the case of a variable porosity that depends upon the volume of colloids retained by the medium, one-dimensional analogue of this model, which is referred to as HLL-model was developed in [3]. In absence of the inverse compaction factor b 2 ½0; 1Þ the scaled form of the HLL-model is as follows: ðxc0 c þ rÞt ¼ acxx cx ; E-mail address: [email protected] (Z. Seyidmamedov (Muradoglu)). 0096-3003/$ - see front matter Ó 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0096-3003(03)00789-6
ð1Þ
440
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
rt ¼ cF ðrÞ;
ð2Þ
0 < x < 1; t > 0;
x ¼ 1 br:
ð3Þ
Here c ¼ cðx; tÞ (scaled on the colloid concentration scale c0 > 0, which is assumed to be small) is the concentration of the colloid particles, measured in volume of colloids in liquid phase per unit volume of the liquid–colloid suspension, and r ¼ rðx; tÞ denotes the volume of retained, or immobile, colloids per unit volume of porous medium. The function x ¼ xðx; tÞ (0 < x < x0 ) denotes the variable porosity of the medium (0 < x < x0 ), a ¼ kV =D, V > 0 the filtration (Darcy) velocity, D > 0 the dispersion constant and k > 0 the filter coefficient; x0 denotes the free porosity, i.e., the porosity of the porous domain when no particulants are present. The function F ¼ F ðrÞ is the accretion rate. We assume that F ð0Þ ¼ 1;
F ðr Þ ¼ 0;
F ðrÞ P 0;
F 0 ðrÞ < 0;
8r 2 ð0; r Þ:
ð4Þ
An existence of wave front solutions of (1) and (2) in the form of c ¼ cðzÞ, r ¼ rðzÞ, where z ¼ x vt, v > 0 is the wave front speed, was shown in [3] and the analysis of wave fronts for parabolic and hyperbolic systems to examine the interactions between dispersion, advection, and absorption is given in [1,5]. In the case of x ¼ 1;
a ¼ 0;
F ðrÞ ¼ 1 br
ð5Þ
an analytical solution of the problem (1)–(3) was obtained in [2]. Later this particular case was analyzed in [3]. Specifically, for the initial boundary-value problem c0 ct þ rt ¼ cx ; rt ¼ cð1 brÞ; x; t > 0; cðx; 0Þ ¼ rðx; 0Þ ¼ 0; x > 0; cð0; tÞ ¼ 1; t > 0; the analytical solution ebs 1 ebs ; rðx; sÞ ¼ ; b ex þ ebs 1 ex þ ebs 1 b 6¼ 0; s ¼ t c0 x; x > 0 cðx; sÞ ¼
was found. However this solution shows that on the characteristics s ¼ t c0 x between the functions cðx; sÞ and rðx; sÞ the relationship 1 cðx; sÞ þ rðx; sÞ ¼ ; b
b 6¼ 0; x > 0
needs to be satisfied, which may not be in practice.
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
441
In Section 2 we consider the mathematical model which includes the inverse compaction factor b. Numerical method for the solution of the coupled system of equations corresponding to this model, is given in Section 3. In the last section we discuss results of computational experiment and provide a comparison of the results with the particular cases obtained in cited works. The computational results show that the parameter b > 0 cannot be neglected, and has high influence to the function r ¼ rðx; tÞ.
2. Small concentration model in presence of the inverse compaction factor We assume that the colloid concentration is small: c0 1 and restrict the problem to be a bounded spatial domain problem. For this aim let the length scale L > 0 be extent of the domain so that the scaled domain will be 0 < x < 1. Then applying the singular perturbation method [3] to the system of Eqs. (1) and (2) we obtain the following nonlinear coupled system: acxx cx cF ðrÞ ¼ 0; rt ¼ cF ðrÞ; x 2 ð0; 1Þ; t > 0: The assumption that the function F ðrÞ is decreasing (see (4)) corresponds to the physical situation that the rate of accretion decreases as more particles are accreted, as might happen if colloid particles can accrete to the matrix, but not to other accreted colloid particles. The case F 0 ðrÞ > 0 corresponds to the situation when F ðrÞ should increase with increasing the volume of retained (or immobile) colloids per unit volume of porous medium. This means as more and more particles are retained, there is greater possibility that additional particles will be filtered. In this situation there does not exist an equilibrium state for the model (1)–(3) at z ¼ 1, where z ¼ x vt and v > 0 is the wave front speed. A more interesting problem occurs when both colmatage (clogging) and decoltomage occur simultaneously. This situation corresponds to the presence of the inverse compaction factor b > 0 and Eq. (2) to take the form rt ¼ cF ðrÞ br. In this case we obtain the following initial boundary-value problem for nonlinear system of coupled equations: acxx cx cF ðrÞ ¼ 0; ð6Þ rt ¼ cF ðrÞ br; x 2 ð0; 1Þ; t > 0; cð0; tÞ ¼ 1; cx ð1; tÞ ¼ 0; t > 0; ð7Þ rðx; 0Þ ¼ 0; x 2 ð0; 1Þ: In the case, when b ¼ 0, the qualitative analysis and an existence of solution to the above problem is given in [3]. We will refer to the nonlinear model (6) and (7) as the HLL(0)-model, if the parameter b is assumed to be zero, and as the HLL(b)-model, if b 6¼ 0. From the point of view of the mathematical
442
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
problems associated with a filtration model of colloid transport, it is important to study the influence of the parameter b to the solution hcðx; tÞ; rðx; tÞi and the behaviour of the solution on the rate of F ðrÞ. The analysis and numerical solution of the nonlinear model is much more difficult then the case b ¼ 0, and we will discuss here only the numerical method. The mathematical analysis of this model will be discussed in a forthcoming work.
3. Discretization and the numerical method for the nonlinear coupled problem (6) and (7) Finite difference discretization of the first equation (6) is performed by using central difference for the first derivative cx . For the discretization of the nonlinear parabolic equation of the coupled system (6) the implicit finite-difference scheme is used, due to its absolute stability [4]. Let us introduce the uniform mesh wh on [0,1]: wh ¼ fxi 2 ½0; 1 : x0 ¼ 0; xi ¼ hi; h ¼ 1=N ; i ¼ 1; N g, for the space variable x 2 ½0; 1. Denote by s > 0 the time step. Then the pair (xi ; tj ), where tj ¼ sj, j ¼ 1; M, will define the current mesh point (i; j) on the xt-plane. On the uniform mesh wh for each time layer tj the first equation (6) can be approximated as follows: acxx;ij cx;ij cij F ðrij Þ ¼ 0;
i ¼ 1; N ;
j ¼ 1; M;
ð8Þ
where cx;ij ¼ ðcðxiþ1 ; tj Þ cðxi1 ; tj ÞÞ=h is the central difference and cxx;ij is the finite difference approximation of the term cxx [4]. Then for each j ¼ 1; M the finite-dimensional analogue of the first equation (6) with the first boundary conditions (7) has the following form: 8 < Ai ci1;j Bi ðrÞci;j þ Di ciþ1;j ¼ /i ; i ¼ 2; N 1; ð9Þ c ¼ j1 c2;j þ m1 ; : 1;j cN ;j ¼ j2 cN 1;j þ m2 ; with the coefficients Ai ¼ 2a þ h; Di ¼ 2a h; Bi ðrÞ ¼ 4a þ 2h2 F ðri;j Þ; j1 ¼ 0; l1 ¼ 0 j2 ¼ 1; l2 ¼ 0:
/i ¼ 0;
ð10Þ
For the given coefficients (10) the linear system of algebraic equations (9) with the tridiagonal matrix can be easily solved by the factorization method. However the coefficient Bi ðrÞ depends on the function r ¼ rðx; tÞ. Finite difference approximation of the second nonlinear equation is performed by using the implicit finite-difference scheme ri;j sci;j F ðri;j Þ þ bsri;j ri;j1 ¼ 0;
j ¼ 2; M;
ð11Þ
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
443
due to its stability for parabolic equations [4]. For the numerical solution of the nonlinear equation (7) with the initial condition ri;1 ¼ 0;
i ¼ 1; N
ð12Þ
we will use the Newton method, due to the above properties of the function F ¼ F ðrÞ. Denoting by f ðri;j Þ the left hand side of the nonlinear equation (7), we obtain: ðkÞ
ðkþ1Þ ri;j
¼
ðkÞ ri;j
f ðri;j Þ ðkÞ
f 0 ðri;j Þ
;
or in the explicit form ðkÞ
ðkþ1Þ ri;j
¼
ðkÞ ri;j
ðmÞ
ðkÞ
ðkÞ
ðkÞ
ri;j ci;j F ðri;j Þ þ bsri;j ri;j1 ðkÞ
1 þ bsF 0 ðri;j Þ
;
j ¼ 2; M:
ð13Þ
The index m shows the coupling of systems (9) and (11) and its sense will be ðmÞ given below. For a given function ci;j we denote by ri;jðmÞ the solution of Eq. (11). ðmÞ Thus the coefficient B ¼ BðrÞ of Eq. (9) depends on the solution ri;j of the nonlinear equation (11), also, the solution ri;j depends on the solution ci;j of Eq. (9). For this reason, in addition to the iteration process (13), one needs to construct the second (outer) iteration process for the system of Eqs. (9) and (11). For the first time layer j ¼ 1 we have the initial condition (12). Taking into account F ð0Þ ¼ 1 we obtain F ðri;1 Þ ¼ 1 and hence for j ¼ 1 Eq. (9) represents the system of linear algebraic equations with Bi ¼ 4a þ 2h2 . Solving this system we can find ci;1 , i ¼ 1; N . As an initial iteration for the second time layer j ¼ 2 ð0Þ ð0Þ we assume ci;2 ¼ ci;1 . Then substituting this in (11) ci;2 ¼ ci;2 and applying the ð0Þ iteration scheme (13) we may find the solution ri;2ð0Þ ¼ ri;2 ðci;2 Þ of Eq. (11), ð0Þ corresponding to the above iteration ci;2 . Substituting the known ri;2ð0Þ in Eq. ð1Þ (9) and solving it we found ci;2 . Then the value is again used in (11) to find ð1Þ ri;2ð1Þ ¼ ri;2 ðci;2 Þ. The above process is stopped when both conditions kri;2ðk0 1Þ ri;2ðk0 Þ kC;h 6 e;
ðk 1Þ
kci;20
ðk Þ
ci;20 kC;h 6 e;
ð14Þ
hold. Here k kC;h is the discrete analogue of C-norm and e > 0 is the given accuracy. ðk Þ If conditions (14) hold, then the pair hci;2 ; ri;2 i ¼ hci;20 ; ri;2ðk0 Þ i is assumed to be the solution of the system (9) and (11), corresponding to the time layer j ¼ 2. Thus, the iteration algorithm for the nonlinear coupled system of Eqs. (9) and (11) can be summarized as follows:
444
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
(i1) For j ¼ 1 solve the system (9) with Bi ¼ 4a þ 2h2 and compute ci;1 ; ð0Þ (i2) For j P 2 apply NewtonÕs scheme (11) with given ci;j ¼ ci;j1 and compute the solution ri;jð0Þ of Eq. (11); ð1Þ (i3) Substitute ri;jð0Þ in (9) and compute ci;j ; ð1Þ (i4) Use ci;j in (11) to find ri;jð1Þ from the iteration scheme (13); (i5) Repeat the process (i3)–(i4) until the fulfilment conditions (14); (i6) For the time layer j þ 1 the above steps (i2)–(i5) need to repeated; ðk Þ
The found hci;j0 ; ri;jðk0 Þ i is assumed to be the solution of the system (9) and (11), corresponding to the time layer j.
4. Results of computational experiments The numerical results presented here can be divided into three groups: accuracy of the discrete model, the influence of the parameter b to the solution, i.e. comparison of the HLL(0)- and HLL(b)-models, and behaviour of the solution on the rate of decrease of the function F ðrÞ. To analyze the accuracy of the method consider the following test. Example 1 2
cðx; tÞ ¼ ext ðx 1Þ ; 1=ð2nþ1Þ
F ðrÞ ¼ 1 r
;
rðx; tÞ ¼ text ;
r 2 ½0; r Þ:
x 2 ð0; 1Þ; t 2 ð0; 1Þ; ð15Þ
The functions cðx; tÞ, rðx; tÞ satisfy the initial and boundary conditions (7), and the function F ðrÞ, given by (15), satisfies conditions (4). Substituting the above functions in Eq. (6) we get: acxx cx F ðrÞc ¼ g1 ðx; tÞ; rt ¼ cF ðrÞ br þ g2 ðx; tÞ; 0 < x < 1; 0 < t < 1; where n h i o 8 < g1 ðx; tÞ ¼ ext ðx 1Þ2 at2 t 1 þ ðtext Þ1=ð2nþ1Þ þ 2ðx 1Þ½2at 1 þ 2a ; n h io : g2 ðx; tÞ ¼ ext 1 þ bt þ xt ðx 1Þ2 1 ðtext Þ1=ð2nþ1Þ : This test problem was solved on the piecewise-uniform mesh with N ¼ 101 and M ¼ 51 (s ¼ 0:02) with the following given data: a ¼ 1, b ¼ 0:5, n ¼ 1. For the accuracy e ¼ 103 given by (14) the relative errors c ¼ kc ch kC;h =kckC;h and r ¼ kr rh kC;h =krkC;h iterations was found c ¼ r 1:7 102 , for the number of iterations k0 ¼ 100 (see (14)). Here r and c are the exact solutions defined on the above introduced mesh and rh and ch are the numerical solutions found by the above iteration scheme. The number of iterations in
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
445
NewtonÕs scheme was about 2–3. For the coarser time mesh with s ¼ 0:05 (i.e. M ¼ 21) the number of iterations to achieve the same accuracy was k0 ¼ 40. This shows that the algorithm is stable and has enough accuracy. To compare the HLL(0)-model, studied in [3], with the HLL(b)-model, given by the system of Eqs. (6) and (7), consider the following. Example 2 (Influence of the inverse compaction factor). Consider the system of Eqs. (6) and (7), with the test functions given in Example 1, for the values b ¼ 0 and b ¼ 0:5 of the inverse compaction factor. For the function F ðrÞ, given by (15), with n ¼ 3, the numerical solutions cðx; tÞ and rðx; tÞ are shown in Figs. 1 and 2, correspondingly. The solid lines correspond to the cases b ¼ 0 and b ¼ 0:5, for t ¼ 1, and the broken lines––to the cases b ¼ 0 and b ¼ 0:5, for t ¼ 0:5. The dots in Fig. 1 corresponds to the case t ¼ 0, where the solutions cðx; tÞ, corresponding to the above models are same. Since the presence of the parameter b in the first equation of system (6) is in the implicit form, the influence of this parameter to the solution cðx; tÞ is negligible. However, by increasing the time variable t this influence increases and at t ¼ 1 the relative deviation dc ¼ kcðx; 1Þb¼0:5 cðx; 1Þb¼0 kC;h =kcðx; tÞb¼0 kC;h is about dc ¼ 0:01, i.e. 1.1 t=0.0 t=0.5 t=1.0 n=3
1.05
1
β =0 β =0.5 β =0 β =0.5
C(x,t)
0.95
0.9
0.85
0.8
0.75
0.7
β =0 β =0.5 0
0.2
0.4
0.6
0.8
1
1.2
x Fig. 1. Increasing profiles of the colloid concentration cðx; tÞ for HLL(0)- and HLL(b)-models, with b ¼ 0:5, at times t ¼ 0, 0.5, 1; a ¼ 1, F ðrÞ ¼ 1 r1=7 .
446
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449 t=0.5 t=1.0 0.3
n=3
0.25 β =0 0.2
σ(x,t)
β =0.5
0.15 β =0 β =0.5 0.1
0.05
0
0
0.2
0.4
0.6
0.8
1
1.2
x Fig. 2. Increasing profiles of the function rðx; tÞ for the data given in Fig. 1.
1%. The significant differences ware observed in the solutions rðx; tÞb¼0:5 and rðx; tÞb¼0 (see Fig. 2). Thus, for t ¼ 0:5 the relative deviation is about 10% and for t ¼ 1:0––about 20%. The numerical results were obtained on the mesh with parameters N ¼ 101 and M ¼ 51 (s ¼ 0:02). The same results were obtained the time mesh with s ¼ 0:05 (M ¼ 21), which shows more time mesh can be used the finite-difference scheme (11). As it was noted above, many of analytic solutions developed in the literature assume the deep bed filtration hypothesis F ðrÞ ¼ 1. One interesting solution which does not make this assumption was obtained in [2] for the case F ðrÞ ¼ 1 br, which means linearization of the dependence F r. In order to study only the influence of the rate of decrease of the function F ðrÞ to the solution hcðx; tÞ; rðx; tÞi let us analyze the HLL(b)-model with different values n ¼ 0 and n ¼ 3 of the parameter n in the function F ðrÞ. Example 3. Consider the initial value problem (6) and (7) for the values n ¼ 0 and n ¼ 3 of the parameter n in the function F ðrÞ. The results of computational experiments for the different time layers t ¼ 0, 0.5, 1.0 are given in Figs. 3 and 4. The figures show that the values of the solution hcðx; tÞ; rðx; tÞi for the linear function F ðrÞ ¼ 1 r (n ¼ 0) (broken lines) and for the nonlinear
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449 1.1
447
β=0 , n=3 β=0 , n=0
1.05
1
C(x,t)
0.95 t=1 t=0.5 0.9
t=1
0.85
0.8
t=0.5
0.75 t=0 0.7
0
0.2
0.4
0.6
0.8
1
1.2
x Fig. 3. The influence of the rate of decrease of the function F ðrÞ to the colloid concentration cðx; tÞ at times t ¼ 0, 0.5, 1.
function (15), with n ¼ 3, are different essentially. For the colloid concentration cðx; tÞ this distinction monotonically increases by increasing the space variable x. It is also interesting to observe that for the nonlinear function F ðrÞ ¼ 1 r1=7 , independently the function cðx; tÞ, the solution rðx; tÞ is a linear function at each time layer (Fig. 4, solid lines). The obtained results show that the influence of the rate of decrease of the function F ðrÞ on the both functions cðx; tÞ and rðx; tÞ is essential. Note that in all computational experiments the value of the parameter a was assumed to be a ¼ 1. Finally, in order to study both the influences of the inverse compaction factor b and the rate of decrease of the function F ðrÞ to the solution hcðx; tÞ; rðx; tÞi we consider the following. Example 4. Consider problem (6) and (7), for the values n ¼ 1 and n ¼ 3 of the parameter n in the function F ðrÞ, and for the value b ¼ 0:5 of the inverse compaction factor. Fig. 5 shows surfaces of the functions cðx; tÞ and rðx; tÞ.
448
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449 β=0 , n=0 β=0 , n=3 0.6 t=1
σ(x,t)
0.5
0.4
t=0.5
0.3
t=1 0.2 t=0.5 0.1
0
0
0.2
0.4
0.6
0.8
1
1.2
x Fig. 4. The influence of the rate of decrease of the function F ðrÞ to the function rðx; tÞ at times t ¼ 0, 0.5, 1.
Fig. 5. Solutions hcðx; tÞ; rðx; tÞi of problem (6) and (7) for b ¼ 0:5. (a1)–(a2) and (b1)–(b2) correspond to the nonlinear functions F ðrÞ ¼ 1 r1=3 and F ðrÞ ¼ 1 r1=7 , respectively.
Z. Seyidmamedov (Muradoglu) / Appl. Math. Comput. 155 (2004) 439–449
449
5. Conclusions We investigated the mathematical model of the transport and accretion of colloid type particles in porous media. Different from previously considered filtration models, the model involve the case when both colmatage (clogging) and decolmatage occur simultaneously, i.e. rt ¼ cF ðrÞ br, where b is the inverse compaction factor. In addition, the accretion rate, i.e. the function F ðrÞ, is assumed to be nonlinear. We present a numerical method for solving the corresponding nonlinear coupled system of differential equations. The numerical results permit to estimate an influence of the inverse compaction factor and nonlinearity of the accretion rate to the solution hcðx; tÞ; rðx; tÞi.
Acknowledgements The author would like to thank Professor A. Hasanoglu (Hasanov) and Dr. S. Cohn for useful suggestions.
References [1] J.L. van Duijin, P. Knabner, Solute transport in a porous media with equilibrium and nonequilibrium multiple-site adsorption: travelling waves, J. Reine Angew. Math. 415 (1991) 1–49. [2] J.P. Herzig, D.M. Leclerc, P. Le Goff, Flow of suspensions through porous media––application to deep filtration, Ind. Eng. Chem. 62 (5) (1970) 8–35. [3] D. Logan, S. Cohn, G. Ledder, On filtration equations in porous media, Commun. Appl. Nonlinear Anal. Mech. 6 (1999) 34–45. [4] A.A. Samarskii, The Theory of Difference Schemes, Marcel Dekker, New York, 2001. [5] A.I. Volpert, V.A. Volpert, Travelling Wave Solutions of Parabolic Systems, American Mathematical Society, Providence, 1994.