Chapter 2 examples in Deen’s book

3.2-1

The equation describes diffusion and reaction in the liquid layer.  The Thiele modulus, ϕ = Da^(1/2), is a convenient form of ratio of the reaction timescale to the diffusion timescale in the entire layer.

ss=DSolve[{θ''[η]- ϕ^2θ[η] ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

θ[η]/.ss[[1]]

ExpToTrig[%]

ans =FullSimplify[%]

TrigToExp[%]

ans1 = Cosh[η ϕ]/Cosh[ϕ]

answerplot =Plot[ans1/.ϕ->{1,10,100},{η,0,1},AxesLabel->{"η","θ"},PlotStyle->{Dashing[{.07,.04}]}]

Some other plot examples

Plot[ans1/.ϕ->{1},{η,0,1},AxesLabel->{"η","θ"},PlotStyle->Red,PlotLegends->{"ϕ=1"},PlotRange->All]

Plot[ans1/.ϕ->{1},{η,0,1},AxesLabel->{"η","θ"},PlotStyle->Dashing[{.02,.02}],PlotLegends->{"ϕ=1"}]

Plot[ans1/.ϕ->{10},{η,0,1},AxesLabel->{"η","θ"},PlotStyle->Green,PlotLegends->{"ϕ=10"},PlotRange->All]

Plot[ans1/.ϕ->{100},{η,0,1},AxesLabel->{"η","θ"},PlotStyle->Blue,PlotLegends->{"ϕ=100"},PlotRange->All]

Show[%,%%,%%%,PlotRange->{0,1}]

We should look at solution for the effect of the reaction close to the interface.  If ϕ is large, then the concentration gradient is increased and hence the flux is increased.  Further, the distance until the solute is depleted is reduced.

TrigToExp[ans1]

What is occurring at large ϕ?  (I could not make this come out from some clever limit, so I just picked the dominate pieces of the numerator and denomenator.)

Exp[ϕ η]/Exp[ϕ]

The exponential is ~0 if the argument gets up to -4.  What value of η is this?

Solve[(η-1) ϕ == -4,η]

Check the green line.  It works and the other values will also work.

1-4/10.

(-4 + ϕ)/ϕ

Apart[%]

So the concentration layer is confined to this thickness.

Now examine the flux

How is the flux increased with chemical reaction?  Many chemical processes use the “trick” of a reaction to either enhance the interfacial transfer, to increase the capacity of the solvent or to shift the equilibrium in a favorable direction (which probably enhances the flux and surely increases capacity.)

Start with the dimensionless flux:

flux =D[ans1,η]

Let’s do our analysis first

TrigToExp[flux]

The second term dominates

%[[2]]

We take the most important piece

ϕ Exp[ϕ (η-1)]

Expand this around η=1

Series[%,{η,1,2}]

So the dimensionless flux will be as large as ϕ.  We can check numbers.

flux/.{η->1,ϕ->100.}

flux/.{η->1,ϕ->10.}

flux/.{η->1,ϕ->5.}

Here are some plots for different values.  You can see that because of the reaction, the flux for a “faster” reaction condition may drop and cross the flux of a slower reaction.  Does this crossing tell us anything interesting?

Plot[flux/.ϕ->{1,2,5},{η,0,1},AxesLabel->{"η","θ"},PlotRange->All]

Is the crossing interesting?  I don’t know.

(flux/.ϕ->ϕ1)-(flux/.ϕ->ϕ2)

Series[%58,{η,0,3}]

%90/.{ϕ1->2.8,ϕ2->1.5}

Solve[Normal[%] ==0,η]

%90/.{ϕ2->ϕ1+ϵ}

Series[Normal[%],{ϵ,0,3}]

FindRoot[(%58/.{ϕ1->2.8,ϕ2->1.5})==0,{η,.7}]

Numerical solution

zz=NDSolve[{θ''[η]- 1^2θ[η] ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot1 =Plot[θ[η]/. zz[[1]],{η,0,1},PlotLegends->{"ϕ=1"},PlotStyle->Green,PlotRange->{0,1}]

zz10=NDSolve[{θ''[η]- 10^2θ[η] ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot10 =Plot[θ[η]/. zz10[[1]],{η,0,1},PlotLegends->{"ϕ=10"},PlotStyle->Red,PlotRange->{0,1}]

zz100=NDSolve[{θ''[η]- 35^2θ[η] ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot100 = Plot[θ[η]/. {zz100[[1]]},{η,0,1},PlotRange->{0,1},PlotLegends->{"ϕ=35"},PlotStyle->Blue,PlotRange->{0,1}]

Show[plot1,plot10,plot100,PlotRange->{0,1}]

answerplot

They almost match

Show[%,%%]

Suppose 2nd order

zz102=NDSolve[{θ''[η]- 10^2θ[η]^2 ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot102 = Plot[θ[η]/. {zz102[[1]]},{η,0,1},PlotRange->{0,1},PlotLegends->{"ϕ=10, 2nd order"},PlotStyle->Orange,PlotRange->{0,1}]

zz10=NDSolve[{θ''[η]- 10^2θ[η] ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot10 = Plot[θ[η]/. {zz10[[1]]},{η,0,1},PlotRange->{0,1},PlotLegends->{"ϕ=10, 1st order"},PlotStyle->Red,PlotRange->{0,1}]

Show[plot10,plot102]

zz1=NDSolve[{θ''[η]- 1θ[η] ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot11 = Plot[θ[η]/. {zz1[[1]]},{η,0,1},PlotRange->{0,1},PlotLegends->{"ϕ=1, 1st order"},PlotStyle->Green,PlotRange->{0,1}]

zz12=NDSolve[{θ''[η]- 1θ[η]^2 ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot12 = Plot[θ[η]/. {zz12[[1]]},{η,0,1},PlotRange->{0,1},PlotLegends->{"ϕ=1, 2nd order"},PlotStyle->Red,PlotRange->{0,1}]

Show[plot11,plot12]

Let’s try a different boundary condition

3.2-1

Set the flux

This introduces the parameter Bi, the Biot number.  It is the ratio of outside transport to inside transport

ss=DSolve[{θ''[η]- ϕ^2θ[η] ==0,θ'[1] == -Bi(θ[1] - 1),θ'[0] == 0},θ[η],η]

θ[η]/.ss[[1]]

ExpToTrig[%]

ans =FullSimplify[%]

answerplot =Plot[ans/.{ϕ->{1,10,100},Bi->1000},{η,0,1},AxesLabel->{"η","θ"}]

answerplot =Plot[ans/.{ϕ->{1,10,100},Bi->10000},{η,0,1},AxesLabel->{"η","θ"}]

answerplot =Plot[ans/.{ϕ->{1,10,100},θθ->1,Bi->1000},{η,.9,1},AxesLabel->{"η","θ"},PlotRange->{.7,1}]

answerplot =Plot[ans/.{ϕ->{1,10,100},θθ->1,Bi->500},{η,.9,1},AxesLabel->{"η","θ"},PlotRange->{.7,1}]

answerplot =Plot[ans/.{ϕ->{1,10,100},θθ->1,Bi->100},{η,.7,1},AxesLabel->{"η","θ"},PlotRange->{0,1}]

answerplot =Plot[ans/.{ϕ->{1,10,100},θθ->1,Bi->10},{η,.7,1},AxesLabel->{"η","θ"},PlotRange->{0,1}]

Now what is the flux?

flux =D[ans,η]

ss1 = TrigToExp[flux]

ss2=Apart[ss1]

ss3 =%[[2]]

ss4=Limit[ss3,Bi->Infinity]

Expand this around η=1

Series[ss4,{η,1,2}]

So we see if Bi->Infinty, then the flux is still ϕ

Suppose the Bi>>1, but not infinity,  will the flux still be ϕ?

ss3

Series[ss3,{η,1,2}]

Bi Exp[2 ϕ] ϕ/(Bi +ϕ)/ Exp[2 ϕ]

Limit[%,ϕ->Infinity]

So as expected, the flux is determined by the Biot number!

answerplot =Plot[flux/.{ϕ->{1,10,100},θθ->1,Bi->10},{η,.7,1},AxesLabel->{"η","θ"},PlotRange->{0,100}]

So the dimensionless flux will be as large as ϕ.  We can check numbers.

flux/.{η->1,ϕ->100.}

flux/.{η->1,ϕ->10.}

flux/.{η->1,ϕ->5.}

Here are some plots for different values.  You can see that because of the reaction, the flux for a “faster” reaction condition may drop and cross the flux of a slower reaction.  Does this crossing tell us anything interesting?

Plot[flux/.ϕ->{1,2,5},{η,0,1},AxesLabel->{"η","θ"},PlotRange->All]

Is the crossing interesting?  I don’t know.

(flux/.ϕ->ϕ1)-(flux/.ϕ->ϕ2)

Series[%58,{η,0,3}]

%90/.{ϕ1->2.8,ϕ2->1.5}

Solve[Normal[%] ==0,η]

%90/.{ϕ2->ϕ1+ϵ}

Series[Normal[%],{ϵ,0,3}]

FindRoot[(%58/.{ϕ1->2.8,ϕ2->1.5})==0,{η,.7}]

Numerical solution

zz=NDSolve[{θ''[η]- 1^2θ[η] ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot1 =Plot[θ[η]/. zz[[1]],{η,0,1},PlotLegends->{"ϕ=1"},PlotStyle->Green,PlotRange->{0,1}]

zz10=NDSolve[{θ''[η]- 10^2θ[η] ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot10 =Plot[θ[η]/. zz10[[1]],{η,0,1},PlotLegends->{"ϕ=10"},PlotStyle->Red,PlotRange->{0,1}]

zz100=NDSolve[{θ''[η]- 35^2θ[η] ==0,θ[1] == 1,θ'[0] == 0},θ[η],η]

plot100 = Plot[θ[η]/. {zz100[[1]]},{η,0,1},PlotRange->{0,1},PlotLegends->{"ϕ=35"},PlotStyle->Blue,PlotRange->{0,1}]

Show[plot1,plot10,plot100,PlotRange->{0,1}]

answerplot

They almost match

Show[%,%%]

3.2-2

This one considers the case when the reaction is reversible.

DSolve[{Da Ca''[y]- k1 Ca[y]  + k2 Cb[y]  ==0,
              Db Cb''[y]+ k1 Ca[y]  - k2 Cb[y]  ==0,       
           Ca[0] == Ca0 , Cb[0] == Cb0,Ca'[L] == 0,Cb'[L] == 0},
             {Ca[y], Cb[y]},y]

cafunc =Ca[y]/.%[[1]]

cbfunc = Cb[y]/.%%[[1]]

Plot[{cafunc,cbfunc}/.{k1->.2,k2->.2,Ca0->1,Cb0->0,L->1,Da->.001,Db->.001},{y,0,1},PlotRange->{0,1}]

Plot[{cafunc,cbfunc}/.{k1->.2,k2->.2,Ca0->1,Cb0->0,L->1,Da->.01,Db->.01},{y,0,1},PlotRange->{0,1}]

3.2-3

DSolve[{D[θ[η],η] ==- Da ϕ^n ( 1- x0 (1-m) θ[η]),θ[0] == 1},θ[η],η]

thetaans = θ[η]/.%[[1]]

Limit[thetaans,m->1]

DSolve[θ'[η]==-Da (ϕ^n) (1-(1-m) x0 θ[η]),{θ[η]},{η}]

thetaans/.x0->1

thetaans/.{Da->.5, m->2,n->2,x0->1}

Plot[thetaans/.{Da->1, m->2,n->1,x0->1,ϕ->1},{η,0,1}]

3.2-5

DSolve[{1/r D[ r D[T[r],r],r] + Hv/k == 0, T'[0] == 0,T'[R] == -h/k (T[R] - Tinf)},T[r],r]

T[r]/.%[[1]]

FullSimplify[%]

Apart[%34]

Apart[%[[1]]]

R^2 Cancel[Together[%38[[1]]/R^2 + %38[[3]]/R^2] ]+ %38[[2]] + Tinf

ans325 = (2 + Bi (1-η^2))/(2 + Bi)

Plot[ans325/.{Bi->{0,.5,2.,10,1000}},{η,0,1}]

3.2-7

DSolve[{1/η D[ η D[θ[η],η],η] - Da θ[η] == 0,θ'[0] == 0,θ[1] ==1},θ[η],η]

thetah = θ[η]/.%[[1]]

cylindrical = Plot[thetah/.{Da->{0.1,1,10,100}},{η,0,1}]

3.2-8

DSolve[{1/η^2 D[ η^2 D[θ[η],η],η] - Da θ[η] == 0,θ'[0] == 0,θ[1] ==1},θ[η],η]

thetas = θ[η]/.%[[1]]

spherical = Plot[thetas/.{Da->{0.1,1,10,100}},{η,0,1},PlotStyle->Dashing[{.1,.05}]]

Show[spherical,cylindrical]

3.2-6

DSolve[{U D[C[y],y] - DA D[C[y],{y,2}] == 0,C'[0] == U (1-K)/DA C[0],C[-δ] ==Cinf},C[y],y]

3.2-9

ss=DSolve[{θ''[η]- Da θ[η]/(κ+ θ[η]) ==0,θ[0] == 1,θ'[1] == 0},θ[η],η]

zzz=NDSolve[({θ''[η]- Da θ[η]/(κ+ θ[η]) ==0,θ[0] == 1,θ'[1] == 0})/.{κ->100,Da->1000},θ[η],η]

plotzzz= Plot[θ[η]/. {zzz[[1]]},{η,0,1},PlotRange->{0,1},PlotLegends->{"κ =100"},PlotStyle->Blue,PlotRange->{0,1}]

zzz=NDSolve[({θ''[η]- Da θ[η]/(κ+ θ[η]) ==0,θ[0] == 1,θ'[1] == 0})/.{κ->10,Da->1000},θ[η],η]

plotzzz= Plot[θ[η]/. {zzz[[1]]},{η,0,1},PlotRange->{0,1},PlotLegends->{"κ =10"},PlotStyle->Blue,PlotRange->{0,1}]

zzz=NDSolve[({θ''[η]- Da θ[η]/(κ+ θ[η]) ==0,θ[0] == 1,θ'[1] == 0})/.{κ->0,Da->10},θ[η],η]

plotzzz= Plot[θ[η]/. {zzz[[1]]},{η,0,1},PlotRange->{0,1},PlotLegends->{"κ =0"},PlotStyle->Blue,PlotRange->{0,1}]

zzz=NDSolve[({θ''[η]- Da θ[η]/(κ+ θ[η]) ==0,θ[0] == 1,θ'[1] == 0})/.{κ->10,Da->1000},θ[η],{η,0,.1}]

plotzzz= Plot[θ[η]/. {zzz[[1]]},{η,0,.1},PlotRange->{0,1},PlotLegends->{"κ =10"},PlotStyle->Blue,PlotRange->{0,1}]

Sqrt[2 Da (1-κ Log[(κ-1)/κ])]

Plot[%207/.Da->{10,100,1000},{κ,0,5},PlotRange->All]