ABSTRACT
Ω = R ′
D∗im (p + λ)
Then, it is required to solve for:
∂z2 −Ωc′ = 0
Applying boundary conditions (23.8), the result is:
c′ = c [
exp ( Ω1/2z
)+ exp (−Ω1/2z)
exp ( Ω1/2B
)+ exp (−Ω1/2B) ]
or:
c′ = c cos h ( Ω1/2z
)
cos h ( Ω1/2B
) (23.9)
The total diffusive exchange flux, q, across the interfacial area separating the mobile and immobile regions for any slab, per unit volume of immobile zone, is obtained differentiating equation (23.9) and making use of Fick’ first law [2]:
q = c g = ( θimD
∂z
∣∣∣∣ z=B
)( 2λd
2Bλd
) = θimD
B
∂c′
∂z
∣∣∣∣ z=B
(23.10)
∂c′
∂z = c cos h
( Ω1/2B
) sin h
( Ω1/2z
) Ω1/2 − 0
cos h2 ( Ω1/2B
) = cΩ1/2 sin h ( Ω1/2z
)
cos h ( Ω1/2B
)
Thus:
∂c′
∂z
∣ ∣ ∣ ∣ z=B
= cΩ1/2 sin h ( Ω1/2B
)
cos h ( Ω1/2B
) = cΩ1/2 tan h (Ω1/2B)
Substituting in (23.10):
c g = θimD ∗ im
B cΩ1/2 tan h
( Ω1/2B
) (23.11)
where the minus sign conventionally used in Fick’s law is neglected. Laplace transforming (23.1) gives:
pc + ν′ ∂c ∂x
− D′x ∂2c
∂2c
∂2c
∂z2 + λc +Φcg = 0
Substituting (23.11):
ν′ ∂c ∂x
− D′x ∂2c
∂2c
∂2c
∂z2 + [ p + λ + ΦθimD
B Ω1/2 tan h
( Ω1/2B
)] c = 0
Let:
Λ = p + λ + ΦθimD ∗ im
B Ω1/2 tan h
( Ω1/2B
)
So:
ν′ ∂c ∂x
− D′x ∂2c
∂2c
∂2c
∂z2 +Λc = 0 (23.12)
Laplace transforming the boundary conditions (23.3):
c(0, y, z, p) = Co p + γ δ(y − y
′)δ(z − z′) c(∞, y, z, p) = 0 c(x,+∞, z, p) = 0 c(x,−∞, z, p) = 0 ∂c
∂z (x, y, 0, p) = 0
∂c
∂z (x, y, T , p) = 0
(23.13)
c = ∞∫
−∞ c(x,α, z, p) exp(−iαy)dy
and apply it to (23.12) and (23.13). This leads to:
ν′ ∂c ∂x
− D′x ∂2c
∂x2 + D′yα2c − D′z
∂2c
∂z2 +Λc = 0 (23.14)
Transforming boundary conditions (23.13):
c(0,α, z, p) = Co p + γ exp(−iαy
′)δ(z − z′) c(∞,α, z, p) = 0
(23.15)
∂c
∂z (x,α, 0, p) = 0
∂c
∂z (x,α, T , p) = 0
(23.16)
Finally, define the finite Fourier cosine transform in z:
c = T∫
c cos (nπz
T
) dz n = 0, 1, 2, . . . ,∞
and apply it to (23.14); thus:
ν′ ∂c ∂x
− D′x ∂2c
∂x2 + D′yα2c + D′z
( n2π2
T
) c +Λc = 0
or:
∂2c
∂x2 − ν
∂x − (
D′x + D
( n2π2
T
) + Λ
)
c = 0
which can be written as:
∂2c
∂x2 − v
∂c
∂x − (
Dx + Dz
( n2π2
T
) + R
)
c = 0 (23.16)
Boundary conditions, from (23.15a):
c(0,α, n, p) = Co p + γ exp(−iαy
′) cos (
T
)
c(∞,α, n, p) = 0 (23.17)
k1 = v Dx
; and k2 = Dyα 2
Dx + Dz
( n2π2
T
) + ΛR
In order to obtain:
∂2c
∂x2 − k1∂c
∂x − k2c = 0
This is a second order linear O.D.E., whose solution, applying boundary conditions (23.17), is:
c = Co p + γ cos
( nπz′
T
) exp
(−iαy′) exp [ x 2 (k1 −
√ k12 + 4k2)
]
Substituting k1, k2, Λ, Ω, Φ and ϕ, we get:
c = Co p + γ cos
( nπz′
T
) exp
(−iαy′)
exp
⎧ ⎪⎪⎨
⎪⎪⎩
x
⎡
⎢⎢ ⎣
ν
Dx −
√√√√√√
) + 4RDx (p + λ)
+ 4θim √
√ p + λ tan h
{ B √
[p + λ] }
⎤
⎥⎥ ⎦
⎫ ⎪⎪⎬
⎪⎪⎭ (23.18)
In order to invert the cosine Fourier transform, let:
ω1 = Co p + γ exp
(−iαy′) ; ω2 = νx2Dx
ω3 = ν 2
Dx + 4Dyα
Dx + 4R
Dx (p + λ) + 4θim
θmDx (b + B) √
p + λ tan h {
B
√ R′
D∗im [p + λ]
}
ω4 = 4Dz Dx
( π2
T
)
Thus:
c = ω1 cos (
T
) exp
[ ω2 − x2
√ ω3 + ω4n2
]
Define the inverse cosine Fourier transform in z:
F−1c [ c (n)
] = c (n = 0)
T + 2
T
n=1 c (n) cos
(nπz T
)
Then,
c = ω1 T
exp ( ω2 − x2
√ ω3
) + 2ω1
T
( ω2 − x2
√ ω3 + 4n2
) cos
(nπz T
) cos
( nπz′
T
)
1, z2) z2 > z1, leads to:
c = ω1 T
exp ( ω2 − x2
√ ω3
) + 2ω1
π
n exp
( ω2 − x2
√ ω3 + ω4n2
)
× cos (nπz
T
) [ sin (nπz2
T
) − sin
(nπz1 T
)]
Substituting ωi-values:
c = Co T (p + γ ) exp
( vx
) exp
(−iαy′) exp
× ⎛
⎝− x 2
√√√√ v 2
Dx + 4Dyα
Dx + 4R
Dx (p + λ) + 4θim
θmDx (b + B) √
p + λ tan h {
B
√ R′
D∗im [p + λ]
}⎞
⎠
+ 2Co π (p + γ ) exp
( νx
) ∞∑
[ 1
n cos
(nπz T
) [ sin (nπz2
T
) − sin
(nπz1 T
)] exp
(−iαy′)
× exp
⎛
⎜⎜ ⎝−
x
√√√√√√
Dx + 4Dyα2Dx + 4RDx (p + λ)
p + λ tan h { B √
[p + λ] }
⎞
⎟⎟ ⎠
⎤
⎥⎥ ⎦
Now, let:
ε1 = Co (p + γ ) exp
( vx
) ; ε2 = 4Dy
Dx ; ε3 = 1
n cos
(nπz T
) [ sin (nπz2
T
) − sin
(nπz1 T
)]
ε4 = ν 2
Dx + 4R
Dx (p + λ) + 4θim
θmDx (b + B) √
p + λ tan h {
B
√ R′
D∗im [p + λ]
}
;
ε5 = 4Dzπ 2n2
DxT + ε4
Hence:
c = ε1 T
exp (−iαy′) exp
⎛
⎝−2 √
( ε4 + ε2α2
) ⎞
⎠
+ 2ε1 π
⎡
⎣ε3 exp (−iαy′) exp
⎛
⎝−2 √
( ε5 + ε2α2
) ⎞
⎠
⎤
⎦
exp ( −a1ς2 − a2
) dς = 1
√ π
a1 exp(−2√a1a2)
→ exp (−2√a1a2 ) = 2
√ a1 π
exp ( −a1ς2 − a2
) dς
Leads to:
c = xε1 2T
√ π
exp ( − x
16 ς2 )
exp ( − ε4
) exp(−iαy′) exp
( −ε2α
) dς
+ xε1 π √
π
⎡
⎣ε3
exp ( − x
16 ς2 )
exp ( − ε5
) exp(−iαy′) exp
( −ε2α
) dς
⎤
⎦
In order to invert the exponential Fourier transform in y, the following inverse is needed [3]:
F−1e [exp(−cα2)] = 1
2 √
cπ exp
( − y
4c
) c > 0
As well as the shift theorem [3]:
F−1e [exp(−iαy′)f (α)] = f (y − y′) Hence,
c = xε1 4Tπ
√ ε2
ς exp { − [
16 + (y − y
] ς2 − ε4
} dς
+ xε1 2π2
√ ε2
⎡
⎣ε3
ς exp ( − [
16 + (y − y
] ς2 − ε5
) dς
⎤
⎦ (23.19)
Let u = ζ 2 to convert equation (23.19) into:
c = xε1 8Tπ
√ ε2
exp { − [ ε2x2 + 4(y − y′)2
] u − ε4
u
} du
+ xε1 4π2
√ ε2
⎡
⎣ε3
exp ( − [ ε2x2 + 4(y − y′)2
] u − ε5
u
) du
⎤
⎦ (23.20)
Define the following identity:
exp ( −Ψ1u − Ψ24u
) du =
√ Ψ2
Ψ1 K1 [√ Ψ1Ψ2
]
c = xε1 √
Tπ √
ε2x2 + 4 (y − y′)2 K1
⎧ ⎨
⎩ 1
√ ε4 [ ε2x2 + 4 (y − y′)2
]
⎫ ⎬
⎭
+ 2xε1 π2 √
ε2x2 + 4 (y − y′)2 ∞∑
⎡
⎣ε3 √
⎧ ⎨
⎩ 1
√ ε5 [ ε2x2 + 4 (y − y′)2
]
⎫ ⎬
⎭
⎤
⎦
Extend the source in the y-direction by integrating y′ from −y0 to y0. This leads to:
c = xε1 √
Tπ
√ ε2x2 + 4 (y − y′)2 K1
⎧ ⎨
⎩ 1
√ ε4 [ ε2x2 + 4 (y − y′)2
]
⎫ ⎬
⎭ dy′
+ 2xε1 π2
⎡
⎢ ⎣ε3
√ ε5
√ ε2x2 + 4 (y − y′)2 K1
⎧ ⎨
⎩ 1
√ ε5 [ ε2x2 + 4 (y − y′)2
]
⎫ ⎬
⎭ dy′
⎤
⎥ ⎦
Substituting εi values:
c(x, y, z, p) = xCo π(p + γ ) exp
( vx
)
⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎪⎩
√
v2 + 4R (p + λ) + 4θim √
√ p + λ tan h
{ B √
[p + λ] } y0∫
⎛
⎜ ⎝
1 √
Dyx2 + Dx (y − y′)2
×K1 ⎛
[ ν2 + 4R(p + λ) + 4θim
p + λ tan h { B √
[p + λ] }] ⎞
⎠
⎞
⎟ ⎠dy′
+ 1 π
⎡
⎢⎢ ⎣
n cos
(nπz T
) [ sin (nπz2
T
) − sin
(nπz1 T
)]
× √√√√4Dzπ
T + v2 + 4R (p + λ) + 4θim
θm (b + B) √
p + λ tan h {
B
√ R′
D∗im [p + λ]
}
× y0∫
1 √
Dyx2 + Dx(y − y′)2 K1
⎛
⎜⎜ ⎝
√√√√√√
T + ν2 + 4Rλ + 4θim
θm(b+B) tan h { B √
λ }]
⎞
⎟⎟ ⎠ dy
⎤
⎥⎥ ⎦
⎫ ⎪⎪⎬
⎪⎪⎭
(23.21)
The Laplace transform of this solution cannot be easily inverted analytically; a numerical method must be used for inversion [4].
