- Email: [email protected]
Approximate analytical solution for seepage flow with fractional derivatives in porous media
Computer methods in applied mechanics and engineering Comput.
ELSEVIER
Approximate
Methods
Appl. Mech. Engg.
167 (1998) 57-68
analytical solution for seepage flow with fractional derivatives in porous media Ji-Huan He
Shanghai
University,
Shanghai Institute of Applird Mathematics Received
17 December
and Mechanics,
1997; accepted
Shanghai 200072, People’s
12 January
Republic of China
1998
Abstract In this paper, a new and more exact model for seepage flow in porous media with fractional derivatives has been proposed, which has modified the well-known Darcy law and overcome the continuity assumption of seepage flow. A new kind of analytical method of nonlinear problem called the variational iteration method is described and used to give approximate solutions of the problem. The results show that the proposed iteration method, requiring no linearization or small perturbation, is very effective and convenient. 0 1998 Elsevier Science S.A. All rights reserved.
1. Introduction It is very difficult to solve a nonlinear problem, either numerically or theoretically, and even more difficult to establish a real model for the nonlinear problem. Much assumption has to be made artificially or unnecessarily to make the practical engineering problems solvable, leading to loss of most important information. The purpose of this paper is manifold. First, we propose a new kind of analytical method of nonlinear problem called the variational iteration method, which, not requiring small parameter in an equation as the perturbation techniques do, has been shown [l-4] to solve effectively, easily, and accurately, a large class of nonlinear problems including the nonlinear fractional differential equations with approximations converging rapidly to accurate solutions. Second, we consider procedures for evaluation of solutions for the differential equations of seepage flow, and its fractional differential partner as well. Fractional derivative [5] of complex order is a generalization of the ordinary concept of derivative of integer order, and the fractional differential model is much more well-suited to physical problems than its differential partner since it makes less unnecessary or over-restricted assumptions which may change the problem being solved, sometimes seriously. In this paper, a fractional differential model for seepage flow in porous will be constructed, attempting to modify the well-known Darcy law and overcome the over-restricted assumption of continuity of seepage flow. Although the subject of the fractional differential equations is very old, it does not appear to have been considered systematically in modern literature, and its applicability range is severely limited, due possibly to the difficulty of solving the problems with fractional derivatives even with numerical simulation. So, most problems are reduced to partial differential equations, leading to loss of most significant information. In this paper, we first begin with the traditional partial differential equations of seepage flow in porous, then consider its fractional differential partner, and their solution procedures via the variational iteration method.
2. The partial differential
equations
of percolation
The partial differential equations for incompressible single phase percolation continuity and Darcy law can be written down in general as follows [6]: 0045-7825/9X/$19.00 0 1998 Elsevier Science S.A. All rights reserved PII: SOO45-7825(98)OOlOS-X
flow under the hypotheses
of
.I. He I Cotnput.
m
Methods
Appl. Mrch.
Etlgrg.
I67 (1998)
57-68
x, y>z)I,=,,= q:,cLy, z)
(4)
where K>, K, and K_ are the percolation and R denotes the percolation domain,
coefficients along the x, y and z direction, respectively, and S 1 + S2 covers all its percolation domain.
3. An introduction
iteration method
to the variational
P is pressure,
The variational iteration method is a generalization of the general Lagrange method proposed by Inokuti et al. [7] in 1978, in the late method, an extremely accurate approximation at some special point can be obtained, but not an analytical one, for details, please see [7] and references cited therein. Now, we will illustrate its basic idea of the variational iteration method by following general partial differential equation: L,u + L>u + L,u + Lp + Nu = g(t, x, y, z)
(5)
where L,, L., , L, and L, are linear operators above t, x, y, and Z, respectively, According to the variational method [l-4], we can write down respectively, in t-, x-, y- and z-directions, respectively, as follows:
and N is a nonlinear operator. following correction functionals
I
u,,+ , 0, x, y, z) = U,,O,x, y>z) + u,,+ , 0, x, y>z) = u,,0, x, y>z) +
I I’A,@,, + I,]’4@p,~+ I 0
/i{L,u,,
+ (L, + L, + L; + N)I?,, - g} dr
(4 + L, + L: + N)& -
s>65
(6b)
CL,+ L, + L: + NG,, -
s>dl
(‘5~)
0
u,,+, 0, .c )‘, z) = u,,(t,x, y, z) +
u,,. ,(L x, .Y>z) = q,(t, x, .Y, z) +
0
(6a)
~,{L ds
(6d)
where A, (i= 1, 2, 3, 4) are general Lagrange multipliers [7], and can be identified optimally by the variational theory [8], ~1” is an initial approximation or trial-function with possible unknowns, i,, are considered as restricted variations [8], i.e. 6U,, =O. For a linear ordinary differential equation, no nonlinear terms exist, so the multiplier can be identified exactly without using the concept of restricted variations, its exact solution, therefore, can be obtained by only one iteration, for example y’+y=sint+t, its correction
Y,, + ,(I)
y(0) = A
functional
Y,,(f) +
(7)
can be written down as follows:
I’
h{)‘,‘,(r) + y,,(r) -
sin7 ~ 4
dr
(8)
(1
Making
the above correction
functional
stationary,
i.e.
J. He
sY,l+,(f) = SY,,(f) + 6
yields following
stationary
I Comput.
Methods
Appl.
Mech.
Engrg.
167 (1998)
59
57-68
h{y,!,(r) + y,,(r) - sin r - r} dr
conditions:
- h’(7) + A(r) = 0
(9a)
1 + A(r)],=, = 0
(9b)
Herein, the prime denotes the derivative with respect to 7. The general solution of (9a) is h(r) = C e’ with a suitable constant C, which can be determined by imposing the ‘boundary condition’ (9b): C = - eP’, so the following iteration formula can be obtained Y,,
+
]
0) = Y,, (t) -
e”-”
{y,:(r) + y,,(r) - sin r - r} dr
We begin with y,,(t) =y(O)=A,
by the above iteration
formula
(10) (lo),
we have
ecrm”{A - sin r - r} dr
y,(t)=A-
1 1 + T (sin t - cos t) + t - 1
(11)
which is the exact solution. For nonlinear problems, due to the approximate identification of the multiplier, its exact solution can be obtained by successive iterations. To illustrate its general process, we will consider the following two examples. EXAMPLE
1
(12) Its t-direction
correction
functional
can be constructed
as follows: (13)
The nonlinear term in the correction stationary conditions read
functional
(13) is considered
as a restricted
variation,
therefore
its
ah(r) ----0 (14) 1 1 :A(r),~=~ The multiplier
= 0
can be identified
as A= - 1, and the following
iteration
formula in t-direction
can be obtained
(15) We start with an initial approximation: successively its approximations:
u,, =u(x, O)=O, by the above iteration
formula
(15), we can obtain
60
J. Hr
I Comput.
Methods
Mech. Engrg.
167 (1998)
57-68
I d7=i2j-fXit3
1
-*‘+,(2xr)’
u2(x, t) = x2t -
Appl.
The results are the same as those obtained by Adomian’s decomposition approximations, we can assume that the solution has the following form:
method
U(X,t) = x’f(t) Substitution
[9]. From the above
(16)
of ( 16) in ( 12) yields
f’(j) +,f’ - I = 0,
f(0) = 0
(17)
It is easy to solve the above nonlinear
ordinary
differential
equation (18)
f(t) = 1 - -&s thus, we obtain
-&J
.(,,t,=xfl
(19)
which is an exact solution. EXAMPLE 2. Duffing equation d’u htL+LL+&=o U(0) = 1, Its correction
(20) U’(0) = 0
functional
can be expressed
u,,+,(t) = u,,(t) +
as follows:
+ u,,(r) + bud:
where U,, is considered as a restricted variation. Making the correction functional (2 1) stationary,
+,
8U,, (f) = 61*,,(t)+ 6
, I [ I,
h
I
(21)
dr
noticing
that &,, = 0,
d’u (7) n + U,,(T) + G;(T) dr dr2 I
= h,,(t) + A(T) &A,‘JT)~,=, - A’(T) ~u,~(T)/_ + ,: (A” + A) Su,, dr = 0 I yields the following
stationary
conditions
A”(7) + A(?-) = 0
(22)
A(r)I = 0 1 A’(r)\,_, = 0 1 The multiplier,
therefore,
A = sin(r - t)
can be determined
as follows: (23)
.I. He I Comput. Methods Appl. Mech. Engrg.
Substitution
61
167 (1998) 57-68
of (23) in (21) results in (24)
Supposing
its initial approximation
has the form
uo(t) = cos cut where cy= a(&) is an unknown function of s with (Y(O)= 1. The initial approximation satisfies all its initial conditions.
By the iteration
formula
The constant cy can be identified by various methods such as method of weighted method, method of collocation, Galerkin method), in this paper we set (Y= \/l + 3&/4 only for simplicity.
(24), we have
residuals
(least square
(27)
The result (26) then leads to the identity
u,(t)=cosLYr+$
I0
sin(T - r) cos 3er dr 6
=
cos at -
(cos 3at - cos t)
4(1-9a’)
with cy defined as (27), which can be approximated cz = Jl Therefore,
+ 3s/4 = 1 + 3&/8 + O(E*) (28) can be approximated
(28) by (29)
as
which is the same as that obtained by perturbation methods. It should be specially pointed out that Eq. (30) is valid only for small 6 while Eq. (28) is uniformly valid for any values of F. The maximal error of the approximate period obtained from Eq. (28) T = 2Am in O<&
4. Variational We consider
iteration formulae the following
dP d2P a*P t +2 +2 = 0, dY
for seepage flow
2-D example (x, y) E f2,
t> 0
(31)
P(t, 0, y) = - 4t + y*
(32a)
-& Pk 0, Y) = 0
(32b)
P(t,x,O)=
(32~)
-4t+x2
-$ P(t, x, 0) = 0
(32d)
P(0, x, y) =x2 + y2
(32e)
62
J. Hr
I Cotnpt.
Methods
Appl.
Med.
Engr,q.
167 (199X) 57-68
the above differential equations are studied by Huang et al. [6] with the Adomian’s Its variational iteration formulae in t-, x- and y-directions can be readily obtained
decomposition
method.
(33)
(35) We start with P,,=P(O,
x, y)=x’
+vl,
by the variational
iteration
formula
in r-direction
(33)
we get (36)
which is an exact solution! We can also assume its initial approximation
has the form
P,#, x, .v) = P(t, 0, J) = - 4t + $ By the variational P,(r,x,~)=
iteration
formula
-4t++2
I
(37) in x-direction
$x)dc=
(34), we obtain
-4t+~‘+x’
which is also an exact solution. It is also very convenient to start with a more general trial-function
with unknowns
like this:
P,,(t, x. v) = C,t + c,x + c,x’ + c&Y + CJ
(38)
where C, (i= 1, 2, 3, 4) are unknowns to be determined. The variational iteration formulae (33) (34) or (35) can be used, or alternately have P,(t, x, y) = c, t + c,x + c,x” + c,y + c,y’ = c,x + c,x’ + c,v + c$ or alternatively
-
used. For example, by (33) we
’ {C, + 2C, + 2C,} dr I0
~ (2C, + 2C,)t
by (35), we get
P,(t,x,~)=C,t+C,x+C,x’+C,~+C,~‘+ = c, t + C&x3
I,;
(s - YHC,+ 2C3 + 2C,l ds
+ c,p - ; (C, + 2C,)$
Imposing the boundary/initial conditions leads also to its exact solution. We can also start with Po(f, x, J) = -4t+_y’+ Ax+&, where A and B are unknowns the iteration formula in x-direction (34), we obtain P, (t, x, y) = - 4t + $ + Ax + (2B - 1)x’
to be determined.
By
(39)
The unknowns can be identified by imposing the boundary condition (32~) and initial condition (32e): A =0 and B = 1, leading also to the exact solution. The approximate solution obtained by Huang et al. in ]6] via Adomian’s decomposition method reads P=
2 2 3+7+.
(
.
2 +7+.
. . (xz++44t) 1
(40)
J. He I Coinput.
Methods
Appl.
Mech.
Engrg.
167 (1998)
57-68
63
though the above series (40) converges to the exact solution, the procedure for evaluation of solutions needs some special techniques and seems to be very cumbersome, while the proposed method is extremely simple and convenient.
5. The fractional
partial differential
equations
of percolation
The above percolation equations (l-4) are deduced under the assumptions of continuity of seepage flow and Darcy law. Generally speaking, these two assumptions are not valid for real seepage flow. The author proposes following modified Darcy law or generalized Darcy law with Riemann-Liouville’s fractional derivatives: 9, = 4
da’P -
axCr’ ’
dW2P qy = Y -
dY”’ ’
O
(414
o
(41b)
O
(4lc)
a*'p q;=&i)z”‘,
in case of (Y1 = cu2= cu3= 1, (4 1a-c) correspond The Riemann-Liouville’s fractional derivative
to Darcy law. is generally defined as [5, lo]
a".44 1 d’ ax"= rtl_ aj dx I o (x -6”4t)dt For example
!!IgxA)= un+1) I‘(h-a+l)X Under the assumption
A-u
,
of continuity
A>-1 of seepage flow, we have following
fractional
differential
equation:
(424 or
(42b) If seepage flow is considered
as rigid body motion, the continuity
equation
can be written down as follows: (43)
with definition ,?‘ulax’=u. Actually the seepage flow is neither continued, be expressed as follows:
nor rigid, so the more general equation
for seepage flow can
(44) where O<,Bl, ,G2, p3<1 It is more exciting that the variational iteration method is also valid for such fractional differential equations, which are very difficult to solve even with numerical simulation. No contribution exists, as far as the author knows, concerning the approximation analytical solution of such problem. Due to the fact that fractional
J. He I Compur. Methods
64
Appl.
Mrch.
Eqrg.
167 (1998)
Z-68
differential equation can excellently describe the natural phenomena, an approximate approach to it has caught much attention by numerous mathematicians and engineers. We will apply the variational iteration method to obtain an analytical solution for a fractional differential equation. Consider first the following system a”u T=f(x+)
(l<(y<2)
(45)
I u(a) = h Though the system is extremely simple, it is remarkably difficult to obtain its approximate analytical solution by traditional nonlinear analytical techniques. According to the variational iteration method, we can construct the following correction functional: u,i + , (xl = q,w
+ I”m)
(46)
where U”(X) is an initial approximation, integrate, defined as follows [5,10]:
F(x) = h[a”u,, li)x” -,f(x, u,,)], and I” is Riemann-Liouville’s
fractional
For example r(A+
I”XA
1)
,i+Ct
zz
r(A+cu+
l)X
,
cu>O,h>
-1
To identify approximately the Lagrange multiplier, we apply restricted variations to nonlinear terms and also to a”ulax” when there exist derivatives with integer orders; But when there exists no derivative with integer order, as far as the author knows, there exists no way to obtain the stationary conditions directly from a functional with Riemann-Liouville’s fractional integrate, failing to determine the multiplier, so, in order to identify approximately the multiplier, some approximation must be made. To do so, we find a minimal integer n = ceil(a) > (Y,or find a maximal integer n = floor(a) < N. For example, in case (Y= 4/3, ceil(u) = 2, floor(a) = 1, so the correction functional (46) can be approximately expressed as follows.
I.4,,+
,(-d = u,,(x) +
1,’{A, [3
-fk
Jj dx
u,,)
(47)
01
The multiplier can be easily determined by the same manipulation as before, substituting Lagrange multipliers, respectively, into (46) results in the following iteration procedures:
U,,+,(X) = U,,(X)+ I”{(B, A, + P*&)[a”%/,xO
-f(x,
where 0, and & are weighted factors with 0, + p, = 1. To illustrate its general process, we consider following EXAMPLE D’l’u
%,I>
the identified
(51)
two examples.
1 II? = x 1/2(u 5-r
+ I),
u(O)=O,O
1
(52)
1. He I Comput. Methods Appl. Mech. Engrg.
its correctional
functional
with Riemand-Liouville’s
fractional
l/2 u,+,(x)
=
u,,(x)
+
D”‘un
Z”2
For ceil( l/2) = 1, its approximate
- & rr
(u, + 1)2
correctional
167 (1998) 57-68
integrate
reads
11
functional
65
(53)
can be expressed
as follows: (54)
and its stationary
conditions
can be readily obtained
d/\(r) 0 a?= 0 1 1 + A(r)I,=, so the multiplier
(55)
can be identified
as: A= - 1, the substitution
of it in (53) results in
l/2
D”‘u-k(u+ Tr
q,+,(x)=u,(x)-Z”2
We start with ~1~(0)=u(O)=O,
1 3 =2x+-x2+Px3 8 while its exact solution u(x)=(l
1)2
by the variational
(56)
> iteration
formula
we have
15 192 [lo] is
-x))“2-
1 =;x++x’+O(x~)
(57)
So, we can see clearly that we can obtain an ideal approximation approximation converges to its exact solution. EXAMPLE
(56)
2. Consider
a more complex
D2u D”2 -+ -+-t++E113=o, Dt2 Dt”’ x(0) = 1, dx(O)ldt = 0,
by its second approximation.
and its nth
system:
O<&<<
1 (58)
The above fractional differential equation can find wide applications in engineering, for example applied to nonlinear oscillator with memory-type damping [lo], or to strong ground vibrations [ll]. Its correction functional can be written down as follows:
it can be
(59) Here, the linear term is also considered evaluation of solutions. Under this condition, as
as a restricted variation in order to simplify the procedures for the natural conditions of above correction functional are expressed
A”( 5) = 0 &$)I,=,
= 0
1 - w315=,
(60) = 0
.I. He I Compur. Mrthd
66
The Lagrange
multiplier
can be identified
Appl.
Mech.
as A= t-r,
I&,( 5) + D”‘u,,W DC’ Supposing
the initial approximation
D5”’
Engrg.
167 (1998)
.57-r%
and the following
iteration
+ u,,( 5) + ~uf( 5)
d5
formula can be obtained:
(61)
has the form
~,~(t) = A e”’ + B e”
(62)
where A, B, a and /3 are unknowns. In view of (61), we have u,(r) = A e”’ + B e” + ,j’([-r)A[(ru2+a”‘+ I + The constant
I
cl’([-t)B[(,B+,f?i”+
l)ePr
l)e”r+,e3’r’]d,$ + E e3”‘] d[ +
’ (5 ~ t){3ABs ec”+B’t (A e”’ + B e”‘)} d[ i0
Q and p can also be identified by the weighted residuals method as mentioned
cyz + (Y”? + I + 3A’Bse’“‘“”
above, herein we set
=0
(63)
and +3,@Ee’“-““=0
@+pr/?+,
(64)
only for simplicity. The above equations (63) and (64) are actually characteristic pseudopolynomial [lo] of linearized (58) when E=O, and can be approximately solved by perturbation methods. We therefore have u,(t) = A emi + Be” =
A
The unknowns A = 1,
e”’
+
B
equation
+ AE ,, ( LJ- t) e3”< dt + BE c: (5 - t) eiBr d< I I
e”’
A and B can be determined
by imposing
the boundary/initial
conditions:
B=++)
Such that we obtain its first approximation without secular terms. It can be found out from (63) and (64) that the natural circular frequency will depend, not only on the amplitude, but also on the time. We now consider 2-D continuous seepage flow which follows the generalized Darcy law (41)
Its correction
functionals
in t-, X- and y-direction
can be written as follows: (67)
(68)
(69) were p,, are considered as restricted variations. The Lagrange multipliers can be identified as follows:
J. He I Comput.
Methods
Appl.
Mech.
Engrg.
167 (1998)
67
57-68
-1
A,=
A2 =
‘-i”? -I
h = 3
ceil(a3)
s-y, 1 -1
floor(a3)
= 1 = 0
Thus, we can obtain its variational
iteration
formulae
in t-, x- and y-direction
P,,+,(t,X,).)=P,~(t,r.y)+jl:{~+~(~)+~(~)}dr
as follows:
(70)
I P,, +*(6 x3 Y) = P,,(L x, Y) +
18,(5-x~-1321{~+~(~)
I0
+$($)}G
(71)
? p,,
+
, 09x> Y> = P,,O> x, Y) +
i 0
(72)
,,,,,-Y,-,,{~+~(~)+~(~)}d~
where p, and /3, are weighted factors with p, + & = 1. As an example we consider the following 1 -D seepage flow !$+A
ax
d’12P ax”2
( )=O
(73)
P(0, x) =x3’?
(74)
3V?i P(t, 0) = - 4 t
(75)
$
(76)
P(t, 0) = 0
its iteration
formula
in t-direction
can be written down as follows: (77)
We start with P,(t, x) = P(0, x)=x~‘~, by the iteration P,(t,x)=x”2-Tt
3J;;;
formula
(77), we get, (78)
which is the exact solution.
6. Conclusion In this paper, the author proposes a very effective and convenient method called variational iteration method to solve (fractional) differential equations, and presents a generalized Darcy law with fractional derivatives and a more general fractional differential equation for seepage flow without the assumption of continuity.
References [I] J.H. He, A new approach to nonlinear partial differential equations, Comm. Nonlinear Sci. Numer. Simul. 2(4) (1997) 230-235, IP address in Internet: http://jiangt.pku.edu.cn/ -nlscomm/, htp://www.pku.edu.cn./docs/xb/xby.html [2] J.H. He, Variational iteration method for delay differential equations, Comm. Nonlinear Sci. Numer. Simul. 2(4) (1997) 235-236. [3] J.H. He, A variational iteration method for nonlinearity and its applications (in Chinese), Mech. Practice 20(I) (1998) 30-32.
68
J. He I Cornput.
Method.\
Appl.
Mrch.
Eqrg.
167 (1998)
57-68
[4] J.H. He, Variational iteration approach to 2.spring system, Mech. Sci. Tech&. 17(2) (1998) 221-223. [51 D. Delbosco, Existence and uniqueness for a nonlinear fractional differential equation, J. Math. Anal. Appl. 204(2) (1996) 609-612. [6] A.X. Huang et al., A new decomposition for solving percolation equations in porous media, 3rd Int. Symp. on Aerothermodynamics of Internal Flows, Beijing, China (1996) 417-420. [7] M. Inokuti et al., General use of the Lagrange multiplier in nonlinear mathematical physics, in: S. Nemat-Nasser, ed., Variational Method in the Mechanics of Solids (Pergamon Press, 1978) 156-162. [8] B.A. Finlayson, The Method of Weighted Residuals and Variational Principles (Academic Press, 1972). [9] G. Adomian, A review of the decomposition method in applied mathematics, J. Math. Anal. Applic. I35 (1988) 510-513. [IO] L.M.B.C. Campos, On the solution of some simple fractional differential equations, Int. J. Math. Math. Sci. l3(3) (1990) 481-496. [I I] R.T. Duarte, A new technique for the analysis of strong ground vibrations and the quantification of earthquake actions, Proc. 7th Europ. Conf. Earth. Eng., Athens, 1983.
ELSEVIER
Approximate
Methods
Appl. Mech. Engg.
167 (1998) 57-68
analytical solution for seepage flow with fractional derivatives in porous media Ji-Huan He
Shanghai
University,
Shanghai Institute of Applird Mathematics Received
17 December
and Mechanics,
1997; accepted
Shanghai 200072, People’s
12 January
Republic of China
1998
Abstract In this paper, a new and more exact model for seepage flow in porous media with fractional derivatives has been proposed, which has modified the well-known Darcy law and overcome the continuity assumption of seepage flow. A new kind of analytical method of nonlinear problem called the variational iteration method is described and used to give approximate solutions of the problem. The results show that the proposed iteration method, requiring no linearization or small perturbation, is very effective and convenient. 0 1998 Elsevier Science S.A. All rights reserved.
1. Introduction It is very difficult to solve a nonlinear problem, either numerically or theoretically, and even more difficult to establish a real model for the nonlinear problem. Much assumption has to be made artificially or unnecessarily to make the practical engineering problems solvable, leading to loss of most important information. The purpose of this paper is manifold. First, we propose a new kind of analytical method of nonlinear problem called the variational iteration method, which, not requiring small parameter in an equation as the perturbation techniques do, has been shown [l-4] to solve effectively, easily, and accurately, a large class of nonlinear problems including the nonlinear fractional differential equations with approximations converging rapidly to accurate solutions. Second, we consider procedures for evaluation of solutions for the differential equations of seepage flow, and its fractional differential partner as well. Fractional derivative [5] of complex order is a generalization of the ordinary concept of derivative of integer order, and the fractional differential model is much more well-suited to physical problems than its differential partner since it makes less unnecessary or over-restricted assumptions which may change the problem being solved, sometimes seriously. In this paper, a fractional differential model for seepage flow in porous will be constructed, attempting to modify the well-known Darcy law and overcome the over-restricted assumption of continuity of seepage flow. Although the subject of the fractional differential equations is very old, it does not appear to have been considered systematically in modern literature, and its applicability range is severely limited, due possibly to the difficulty of solving the problems with fractional derivatives even with numerical simulation. So, most problems are reduced to partial differential equations, leading to loss of most significant information. In this paper, we first begin with the traditional partial differential equations of seepage flow in porous, then consider its fractional differential partner, and their solution procedures via the variational iteration method.
2. The partial differential
equations
of percolation
The partial differential equations for incompressible single phase percolation continuity and Darcy law can be written down in general as follows [6]: 0045-7825/9X/$19.00 0 1998 Elsevier Science S.A. All rights reserved PII: SOO45-7825(98)OOlOS-X
flow under the hypotheses
of
.I. He I Cotnput.
m
Methods
Appl. Mrch.
Etlgrg.
I67 (1998)
57-68
x, y>z)I,=,,= q:,cLy, z)
(4)
where K>, K, and K_ are the percolation and R denotes the percolation domain,
coefficients along the x, y and z direction, respectively, and S 1 + S2 covers all its percolation domain.
3. An introduction
iteration method
to the variational
P is pressure,
The variational iteration method is a generalization of the general Lagrange method proposed by Inokuti et al. [7] in 1978, in the late method, an extremely accurate approximation at some special point can be obtained, but not an analytical one, for details, please see [7] and references cited therein. Now, we will illustrate its basic idea of the variational iteration method by following general partial differential equation: L,u + L>u + L,u + Lp + Nu = g(t, x, y, z)
(5)
where L,, L., , L, and L, are linear operators above t, x, y, and Z, respectively, According to the variational method [l-4], we can write down respectively, in t-, x-, y- and z-directions, respectively, as follows:
and N is a nonlinear operator. following correction functionals
I
u,,+ , 0, x, y, z) = U,,O,x, y>z) + u,,+ , 0, x, y>z) = u,,0, x, y>z) +
I I’A,@,, + I,]’4@p,~+ I 0
/i{L,u,,
+ (L, + L, + L; + N)I?,, - g} dr
(4 + L, + L: + N)& -
s>65
(6b)
CL,+ L, + L: + NG,, -
s>dl
(‘5~)
0
u,,+, 0, .c )‘, z) = u,,(t,x, y, z) +
u,,. ,(L x, .Y>z) = q,(t, x, .Y, z) +
0
(6a)
~,{L
(6d)
where A, (i= 1, 2, 3, 4) are general Lagrange multipliers [7], and can be identified optimally by the variational theory [8], ~1” is an initial approximation or trial-function with possible unknowns, i,, are considered as restricted variations [8], i.e. 6U,, =O. For a linear ordinary differential equation, no nonlinear terms exist, so the multiplier can be identified exactly without using the concept of restricted variations, its exact solution, therefore, can be obtained by only one iteration, for example y’+y=sint+t, its correction
Y,, + ,(I)
y(0) = A
functional
Y,,(f) +
(7)
can be written down as follows:
I’
h{)‘,‘,(r) + y,,(r) -
sin7 ~ 4
dr
(8)
(1
Making
the above correction
functional
stationary,
i.e.
J. He
sY,l+,(f) = SY,,(f) + 6
yields following
stationary
I Comput.
Methods
Appl.
Mech.
Engrg.
167 (1998)
59
57-68
h{y,!,(r) + y,,(r) - sin r - r} dr
conditions:
- h’(7) + A(r) = 0
(9a)
1 + A(r)],=, = 0
(9b)
Herein, the prime denotes the derivative with respect to 7. The general solution of (9a) is h(r) = C e’ with a suitable constant C, which can be determined by imposing the ‘boundary condition’ (9b): C = - eP’, so the following iteration formula can be obtained Y,,
+
]
0) = Y,, (t) -
e”-”
{y,:(r) + y,,(r) - sin r - r} dr
We begin with y,,(t) =y(O)=A,
by the above iteration
formula
(10) (lo),
we have
ecrm”{A - sin r - r} dr
y,(t)=A-
1 1 + T (sin t - cos t) + t - 1
(11)
which is the exact solution. For nonlinear problems, due to the approximate identification of the multiplier, its exact solution can be obtained by successive iterations. To illustrate its general process, we will consider the following two examples. EXAMPLE
1
(12) Its t-direction
correction
functional
can be constructed
as follows: (13)
The nonlinear term in the correction stationary conditions read
functional
(13) is considered
as a restricted
variation,
therefore
its
ah(r) ----0 (14) 1 1 :A(r),~=~ The multiplier
= 0
can be identified
as A= - 1, and the following
iteration
formula in t-direction
can be obtained
(15) We start with an initial approximation: successively its approximations:
u,, =u(x, O)=O, by the above iteration
formula
(15), we can obtain
60
J. Hr
I Comput.
Methods
Mech. Engrg.
167 (1998)
57-68
I d7=i2j-fXit3
1
-*‘+,(2xr)’
u2(x, t) = x2t -
Appl.
The results are the same as those obtained by Adomian’s decomposition approximations, we can assume that the solution has the following form:
method
U(X,t) = x’f(t) Substitution
[9]. From the above
(16)
of ( 16) in ( 12) yields
f’(j) +,f’ - I = 0,
f(0) = 0
(17)
It is easy to solve the above nonlinear
ordinary
differential
equation (18)
f(t) = 1 - -&s thus, we obtain
-&J
.(,,t,=xfl
(19)
which is an exact solution. EXAMPLE 2. Duffing equation d’u htL+LL+&=o U(0) = 1, Its correction
(20) U’(0) = 0
functional
can be expressed
u,,+,(t) = u,,(t) +
as follows:
+ u,,(r) + bud:
where U,, is considered as a restricted variation. Making the correction functional (2 1) stationary,
+,
8U,, (f) = 61*,,(t)+ 6
, I [ I,
h
I
(21)
dr
noticing
that &,, = 0,
d’u (7) n + U,,(T) + G;(T) dr dr2 I
= h,,(t) + A(T) &A,‘JT)~,=, - A’(T) ~u,~(T)/_ + ,: (A” + A) Su,, dr = 0 I yields the following
stationary
conditions
A”(7) + A(?-) = 0
(22)
A(r)I = 0 1 A’(r)\,_, = 0 1 The multiplier,
therefore,
A = sin(r - t)
can be determined
as follows: (23)
.I. He I Comput. Methods Appl. Mech. Engrg.
Substitution
61
167 (1998) 57-68
of (23) in (21) results in (24)
Supposing
its initial approximation
has the form
uo(t) = cos cut where cy= a(&) is an unknown function of s with (Y(O)= 1. The initial approximation satisfies all its initial conditions.
By the iteration
formula
The constant cy can be identified by various methods such as method of weighted method, method of collocation, Galerkin method), in this paper we set (Y= \/l + 3&/4 only for simplicity.
(24), we have
residuals
(least square
(27)
The result (26) then leads to the identity
u,(t)=cosLYr+$
I0
sin(T - r) cos 3er dr 6
=
cos at -
(cos 3at - cos t)
4(1-9a’)
with cy defined as (27), which can be approximated cz = Jl Therefore,
+ 3s/4 = 1 + 3&/8 + O(E*) (28) can be approximated
(28) by (29)
as
which is the same as that obtained by perturbation methods. It should be specially pointed out that Eq. (30) is valid only for small 6 while Eq. (28) is uniformly valid for any values of F. The maximal error of the approximate period obtained from Eq. (28) T = 2Am in O<&
4. Variational We consider
iteration formulae the following
dP d2P a*P t +2 +2 = 0, dY
for seepage flow
2-D example (x, y) E f2,
t> 0
(31)
P(t, 0, y) = - 4t + y*
(32a)
-& Pk 0, Y) = 0
(32b)
P(t,x,O)=
(32~)
-4t+x2
-$ P(t, x, 0) = 0
(32d)
P(0, x, y) =x2 + y2
(32e)
62
J. Hr
I Cotnpt.
Methods
Appl.
Med.
Engr,q.
167 (199X) 57-68
the above differential equations are studied by Huang et al. [6] with the Adomian’s Its variational iteration formulae in t-, x- and y-directions can be readily obtained
decomposition
method.
(33)
(35) We start with P,,=P(O,
x, y)=x’
+vl,
by the variational
iteration
formula
in r-direction
(33)
we get (36)
which is an exact solution! We can also assume its initial approximation
has the form
P,#, x, .v) = P(t, 0, J) = - 4t + $ By the variational P,(r,x,~)=
iteration
formula
-4t++2
I
(37) in x-direction
$x)dc=
(34), we obtain
-4t+~‘+x’
which is also an exact solution. It is also very convenient to start with a more general trial-function
with unknowns
like this:
P,,(t, x. v) = C,t + c,x + c,x’ + c&Y + CJ
(38)
where C, (i= 1, 2, 3, 4) are unknowns to be determined. The variational iteration formulae (33) (34) or (35) can be used, or alternately have P,(t, x, y) = c, t + c,x + c,x” + c,y + c,y’ = c,x + c,x’ + c,v + c$ or alternatively
-
used. For example, by (33) we
’ {C, + 2C, + 2C,} dr I0
~ (2C, + 2C,)t
by (35), we get
P,(t,x,~)=C,t+C,x+C,x’+C,~+C,~‘+ = c, t + C&x3
I,;
(s - YHC,+ 2C3 + 2C,l ds
+ c,p - ; (C, + 2C,)$
Imposing the boundary/initial conditions leads also to its exact solution. We can also start with Po(f, x, J) = -4t+_y’+ Ax+&, where A and B are unknowns the iteration formula in x-direction (34), we obtain P, (t, x, y) = - 4t + $ + Ax + (2B - 1)x’
to be determined.
By
(39)
The unknowns can be identified by imposing the boundary condition (32~) and initial condition (32e): A =0 and B = 1, leading also to the exact solution. The approximate solution obtained by Huang et al. in ]6] via Adomian’s decomposition method reads P=
2 2 3+7+.
(
.
2 +7+.
. . (xz++44t) 1
(40)
J. He I Coinput.
Methods
Appl.
Mech.
Engrg.
167 (1998)
57-68
63
though the above series (40) converges to the exact solution, the procedure for evaluation of solutions needs some special techniques and seems to be very cumbersome, while the proposed method is extremely simple and convenient.
5. The fractional
partial differential
equations
of percolation
The above percolation equations (l-4) are deduced under the assumptions of continuity of seepage flow and Darcy law. Generally speaking, these two assumptions are not valid for real seepage flow. The author proposes following modified Darcy law or generalized Darcy law with Riemann-Liouville’s fractional derivatives: 9, = 4
da’P -
axCr’ ’
dW2P qy = Y -
dY”’ ’
O
(414
o
(41b)
O
(4lc)
a*'p q;=&i)z”‘,
in case of (Y1 = cu2= cu3= 1, (4 1a-c) correspond The Riemann-Liouville’s fractional derivative
to Darcy law. is generally defined as [5, lo]
a".44 1 d’ ax"= rtl_ aj dx I o (x -6”4t)dt For example
!!IgxA)= un+1) I‘(h-a+l)X Under the assumption
A-u
,
of continuity
A>-1 of seepage flow, we have following
fractional
differential
equation:
(424 or
(42b) If seepage flow is considered
as rigid body motion, the continuity
equation
can be written down as follows: (43)
with definition ,?‘ulax’=u. Actually the seepage flow is neither continued, be expressed as follows:
nor rigid, so the more general equation
for seepage flow can
(44) where O<,Bl, ,G2, p3<1 It is more exciting that the variational iteration method is also valid for such fractional differential equations, which are very difficult to solve even with numerical simulation. No contribution exists, as far as the author knows, concerning the approximation analytical solution of such problem. Due to the fact that fractional
J. He I Compur. Methods
64
Appl.
Mrch.
Eqrg.
167 (1998)
Z-68
differential equation can excellently describe the natural phenomena, an approximate approach to it has caught much attention by numerous mathematicians and engineers. We will apply the variational iteration method to obtain an analytical solution for a fractional differential equation. Consider first the following system a”u T=f(x+)
(l<(y<2)
(45)
I u(a) = h Though the system is extremely simple, it is remarkably difficult to obtain its approximate analytical solution by traditional nonlinear analytical techniques. According to the variational iteration method, we can construct the following correction functional: u,i + , (xl = q,w
+ I”m)
(46)
where U”(X) is an initial approximation, integrate, defined as follows [5,10]:
F(x) = h[a”u,, li)x” -,f(x, u,,)], and I” is Riemann-Liouville’s
fractional
For example r(A+
I”XA
1)
,i+Ct
zz
r(A+cu+
l)X
,
cu>O,h>
-1
To identify approximately the Lagrange multiplier, we apply restricted variations to nonlinear terms and also to a”ulax” when there exist derivatives with integer orders; But when there exists no derivative with integer order, as far as the author knows, there exists no way to obtain the stationary conditions directly from a functional with Riemann-Liouville’s fractional integrate, failing to determine the multiplier, so, in order to identify approximately the multiplier, some approximation must be made. To do so, we find a minimal integer n = ceil(a) > (Y,or find a maximal integer n = floor(a) < N. For example, in case (Y= 4/3, ceil(u) = 2, floor(a) = 1, so the correction functional (46) can be approximately expressed as follows.
I.4,,+
,(-d = u,,(x) +
1,’{A, [3
-fk
Jj dx
u,,)
(47)
01
The multiplier can be easily determined by the same manipulation as before, substituting Lagrange multipliers, respectively, into (46) results in the following iteration procedures:
U,,+,(X) = U,,(X)+ I”{(B, A, + P*&)[a”%/,xO
-f(x,
where 0, and & are weighted factors with 0, + p, = 1. To illustrate its general process, we consider following EXAMPLE D’l’u
%,I>
the identified
(51)
two examples.
1 II? = x 1/2(u 5-r
+ I),
u(O)=O,O
1
(52)
1. He I Comput. Methods Appl. Mech. Engrg.
its correctional
functional
with Riemand-Liouville’s
fractional
l/2 u,+,(x)
=
u,,(x)
+
D”‘un
Z”2
For ceil( l/2) = 1, its approximate
- & rr
(u, + 1)2
correctional
167 (1998) 57-68
integrate
reads
11
functional
65
(53)
can be expressed
as follows: (54)
and its stationary
conditions
can be readily obtained
d/\(r) 0 a?= 0 1 1 + A(r)I,=, so the multiplier
(55)
can be identified
as: A= - 1, the substitution
of it in (53) results in
l/2
D”‘u-k(u+ Tr
q,+,(x)=u,(x)-Z”2
We start with ~1~(0)=u(O)=O,
1 3 =2x+-x2+Px3 8 while its exact solution u(x)=(l
1)2
by the variational
(56)
> iteration
formula
we have
15 192 [lo] is
-x))“2-
1 =;x++x’+O(x~)
(57)
So, we can see clearly that we can obtain an ideal approximation approximation converges to its exact solution. EXAMPLE
(56)
2. Consider
a more complex
D2u D”2 -+ -+-t++E113=o, Dt2 Dt”’ x(0) = 1, dx(O)ldt = 0,
by its second approximation.
and its nth
system:
O<&<<
1 (58)
The above fractional differential equation can find wide applications in engineering, for example applied to nonlinear oscillator with memory-type damping [lo], or to strong ground vibrations [ll]. Its correction functional can be written down as follows:
it can be
(59) Here, the linear term is also considered evaluation of solutions. Under this condition, as
as a restricted variation in order to simplify the procedures for the natural conditions of above correction functional are expressed
A”( 5) = 0 &$)I,=,
= 0
1 - w315=,
(60) = 0
.I. He I Compur. Mrthd
66
The Lagrange
multiplier
can be identified
Appl.
Mech.
as A= t-r,
I&,( 5) + D”‘u,,W DC’ Supposing
the initial approximation
D5”’
Engrg.
167 (1998)
.57-r%
and the following
iteration
+ u,,( 5) + ~uf( 5)
d5
formula can be obtained:
(61)
has the form
~,~(t) = A e”’ + B e”
(62)
where A, B, a and /3 are unknowns. In view of (61), we have u,(r) = A e”’ + B e” + ,j’([-r)A[(ru2+a”‘+ I + The constant
I
cl’([-t)B[(,B+,f?i”+
l)ePr
l)e”r+,e3’r’]d,$ + E e3”‘] d[ +
’ (5 ~ t){3ABs ec”+B’t (A e”’ + B e”‘)} d[ i0
Q and p can also be identified by the weighted residuals method as mentioned
cyz + (Y”? + I + 3A’Bse’“‘“”
above, herein we set
=0
(63)
and +3,@Ee’“-““=0
@+pr/?+,
(64)
only for simplicity. The above equations (63) and (64) are actually characteristic pseudopolynomial [lo] of linearized (58) when E=O, and can be approximately solved by perturbation methods. We therefore have u,(t) = A emi + Be” =
A
The unknowns A = 1,
e”’
+
B
equation
+ AE ,, ( LJ- t) e3”< dt + BE c: (5 - t) eiBr d< I I
e”’
A and B can be determined
by imposing
the boundary/initial
conditions:
B=++)
Such that we obtain its first approximation without secular terms. It can be found out from (63) and (64) that the natural circular frequency will depend, not only on the amplitude, but also on the time. We now consider 2-D continuous seepage flow which follows the generalized Darcy law (41)
Its correction
functionals
in t-, X- and y-direction
can be written as follows: (67)
(68)
(69) were p,, are considered as restricted variations. The Lagrange multipliers can be identified as follows:
J. He I Comput.
Methods
Appl.
Mech.
Engrg.
167 (1998)
67
57-68
-1
A,=
A2 =
‘-i”? -I
h = 3
ceil(a3)
s-y, 1 -1
floor(a3)
= 1 = 0
Thus, we can obtain its variational
iteration
formulae
in t-, x- and y-direction
P,,+,(t,X,).)=P,~(t,r.y)+jl:{~+~(~)+~(~)}dr
as follows:
(70)
I P,, +*(6 x3 Y) = P,,(L x, Y) +
18,(5-x~-1321{~+~(~)
I0
+$($)}G
(71)
? p,,
+
, 09x> Y> = P,,O> x, Y) +
i 0
(72)
,,,,,-Y,-,,{~+~(~)+~(~)}d~
where p, and /3, are weighted factors with p, + & = 1. As an example we consider the following 1 -D seepage flow !$+A
ax
d’12P ax”2
( )=O
(73)
P(0, x) =x3’?
(74)
3V?i P(t, 0) = - 4 t
(75)
$
(76)
P(t, 0) = 0
its iteration
formula
in t-direction
can be written down as follows: (77)
We start with P,(t, x) = P(0, x)=x~‘~, by the iteration P,(t,x)=x”2-Tt
3J;;;
formula
(77), we get, (78)
which is the exact solution.
6. Conclusion In this paper, the author proposes a very effective and convenient method called variational iteration method to solve (fractional) differential equations, and presents a generalized Darcy law with fractional derivatives and a more general fractional differential equation for seepage flow without the assumption of continuity.
References [I] J.H. He, A new approach to nonlinear partial differential equations, Comm. Nonlinear Sci. Numer. Simul. 2(4) (1997) 230-235, IP address in Internet: http://jiangt.pku.edu.cn/ -nlscomm/, htp://www.pku.edu.cn./docs/xb/xby.html [2] J.H. He, Variational iteration method for delay differential equations, Comm. Nonlinear Sci. Numer. Simul. 2(4) (1997) 235-236. [3] J.H. He, A variational iteration method for nonlinearity and its applications (in Chinese), Mech. Practice 20(I) (1998) 30-32.
68
J. He I Cornput.
Method.\
Appl.
Mrch.
Eqrg.
167 (1998)
57-68
[4] J.H. He, Variational iteration approach to 2.spring system, Mech. Sci. Tech&. 17(2) (1998) 221-223. [51 D. Delbosco, Existence and uniqueness for a nonlinear fractional differential equation, J. Math. Anal. Appl. 204(2) (1996) 609-612. [6] A.X. Huang et al., A new decomposition for solving percolation equations in porous media, 3rd Int. Symp. on Aerothermodynamics of Internal Flows, Beijing, China (1996) 417-420. [7] M. Inokuti et al., General use of the Lagrange multiplier in nonlinear mathematical physics, in: S. Nemat-Nasser, ed., Variational Method in the Mechanics of Solids (Pergamon Press, 1978) 156-162. [8] B.A. Finlayson, The Method of Weighted Residuals and Variational Principles (Academic Press, 1972). [9] G. Adomian, A review of the decomposition method in applied mathematics, J. Math. Anal. Applic. I35 (1988) 510-513. [IO] L.M.B.C. Campos, On the solution of some simple fractional differential equations, Int. J. Math. Math. Sci. l3(3) (1990) 481-496. [I I] R.T. Duarte, A new technique for the analysis of strong ground vibrations and the quantification of earthquake actions, Proc. 7th Europ. Conf. Earth. Eng., Athens, 1983.