The governing equation is shown as follows:
I first try to employ the NDSolve, but it seems that Mathematica can not handle the fourth boundary condition.
Therefore, I rewrite the code with finite difference method, but still have problem
Mr = 100; Mz = 10; Lr = 10; Lz = 1;
γ = 0.01;
α = 1;
κ = 1;
dr = Lr/Mr;
dz = Lz/dz;
dt = 0.1;
r[i_] := i*dr;
z[j_] := j*dz;
(*i.c.*)
u[i_, j_, 0] := 0;
(*b.c.*)
u[1, j_, k_] := u[2, j, k] + 2;
u[Mr, j_, k_] := u[Mr - 1, j, k];
u[i_, 1, k_] := u[i, 2, k];
u[i_, Mz, k_] :=
If[k == 1, 0,
1/dz^2 (dt dz γ u[i, -1 + Mz, -1 + k] +
dt α γ u[i, -1 + Mz, -1 + k]^2 +
dz^2 u[i, Mz, -1 + k] - dt dz γ u[i, Mz, -1 + k] -
2 dt α γ u[i, -1 + Mz, -1 + k] u[i, Mz, -1 + k] +
dt α γ u[i, Mz, -1 + k]^2)];
u[i_, j_, k_] :=
u[i, j, k] =
1/(dr^2 dz^2 i) (dt dz^2 i u[-1 + i, j, -1 + k] +
dr^2 dt i u[i, -1 + j, -1 + k] - dt dz^2 u[i, j, -1 + k] -
2 dr^2 dt i u[i, j, -1 + k] + dr^2 dz^2 i u[i, j, -1 + k] -
2 dt dz^2 i u[i, j, -1 + k] + dr^2 dt i u[i, 1 + j, -1 + k] +
dt dz^2 u[1 + i, j, -1 + k] + dt dz^2 i u[1 + i, j, -1 + k]);
u[1, Mz, 1]
However, I receive the following error message
$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>
$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>
$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>
General::stop: Further output of $RecursionLimit::reclim will be suppressed during this calculation. >>
The problem may be induced by the imposion of nonlinear and time-derivative and nonlinear boundary condition. Any idea about how to solve my problem? This boundary condition has bothered me for long time.
I re code the provided code as follows
Clear[fdd, pdetoode, tooderule]
fdd[{}, grid_, value_, order_] := value;
fdd[a__] := NDSolve`FiniteDifferenceDerivative@a;
pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]],
rest];
pdetoode[{func__}[var__], rest__] :=
pdetoode[Alternatives[func][var], rest];
pdetoode[rest__, grid_?VectorQ, o_Integer] :=
pdetoode[rest, {grid}, o];
pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid},
pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
MapThread[
Piecewise[{{1, # === #2[[1]]}, {-1, # === #2[[-1]]}},
All] &, {coord, bound}]]},
tooderule@
Flatten@{((u : func) |
Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat,
t_, x2___] :> (Sow@coordtoindex@{x1, x2};
fdd[{dx1, dx2}, {grid},
Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]},
Outer[Slot@i &, grid]]}]]];
tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //.
eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@
Reap[expr /. rule]
Clear@pdetoae;
pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] :=
pdetoae[Alternatives[func][var], rest];
pdetoae[func_[var__], rest__] :=
Module[{t},
Function[
pde, #[pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :>
u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
pdetoode[func[var, t], t, rest]];
\[Gamma] = 100;
\[Alpha] = 0;
\[Kappa] = 1;
R = 10;
Z = 1;
eps = 10^-1;
tend = 1000;
eq = With[{u = u[r, Sqrt[\[Kappa]] z, t]},
Laplacian[u, {r, th, z}, "Cylindrical"] == D[u, t] /.
Sqrt[\[Kappa]] z -> z];
ic = u[r, z, 0] == 0;
bc = With[{u = u[r, z, t]}, {D[u, r] == -2/r /. r -> eps,
u == 0 /. r -> R, D[u, z] == 0 /. z -> 0, D[u, z] == 0 /. z -> Z}]
domain@r = {eps, R};
domain@z = {0, Z};
points@r = 50;
points@z = 50;
difforder = 2;
(grid@# = Array[# &, points@#, domain@#]) & /@ {r, z};
(*Definition of pdetoode isn't included in this post,please find it \
in the link above.*)
ptoofunc = pdetoode[u[r, z, t], t, grid /@ {r, z}, difforder];
delbothside = #[[2 ;; -2]] &;
ode = delbothside /@ delbothside@ptoofunc@eq;
odeic = delbothside /@ delbothside@ptoofunc@ic;
odebc = MapAt[delbothside, ptoofunc@bc, {{1}, {2}}];
sollst = NDSolveValue[{ode, odeic, odebc},
Outer[u, grid@r, grid@z] // Flatten, {t, 0, tend}];
sol = ListInterpolation[
Partition[Developer`ToPackedArray@#["ValuesOnGrid"] & /@ #,
points@z], {grid@r, grid@z, #[[1]]["Coordinates"][[1]]}] &@
sollst; // AbsoluteTiming
Manipulate[
Plot3D[sol[r, z, t], {r, eps, R}, {z, 0, Z},
PlotRange -> {-1, 10}], {t, 0, tend}]
The equation indicate that the upper and under boundary conditions are subjected to no-flow condition. And it can easily be solved and obtain an analytical solution.
To compare the analytical solution (or semi-analytical), the solution can be expressed as
Vi[n_, i_] :=
Vi[n, i] = (-1)^(i + n/2) Sum[
k^(n/2) (2 k)! /( (n/2 - k)! k! (k - 1)! (i - k)! (2 k -
i)! ), { k, Floor[ (i + 1)/2 ], Min[ i, n/2] } ] // N;
Stehfest[F_, s_, t_, n_: 16] :=
If[n > 16, Message[Stehfest::optimalterms, n];
If[ OddQ[n], Message[Stehfest::odd, n];
"Enter an even number of terms",
If[n > 32, Message[Stehfest::terms, n];
" Try a smaller value for n. Maximum allowable n is 32 ",
Log[2]/t Sum[ Vi[n, i]*F /. s -> i Log[2]/t , {i, 1, n} ] ]],
If[ OddQ[n], Message[Stehfest::odd, n];
"Enter an even number of terms",
If[n > 32, Message[Stehfest::terms, n];
" Try a smaller value for n. Maximum allowable n is 32.",
Log[2]/t Sum[
Vi[n, i]*F /. s -> i Log[2]/t , {i, 1, n} ] ]]] // N;
s0[r_, z_, t_] :=
NIntegrate[
Re[Stehfest[(
2 E^(-((Sqrt[a^2 + p] z)/
Sqrt[\[Kappa]])) (-E^((Sqrt[a^2 + p]/Sqrt[\[Kappa]]))
p \[Gamma] -
E^((Sqrt[a^2 + p] (1 + 2 z))/Sqrt[\[Kappa]]) p \[Gamma] +
E^((Sqrt[a^2 + p] z)/
Sqrt[\[Kappa]]) (p \[Gamma] - Sqrt[a^2 + p] Sqrt[\[Kappa]]) +
E^((Sqrt[a^2 + p] (2 + z))/
Sqrt[\[Kappa]]) (p \[Gamma] +
Sqrt[a^2 + p] Sqrt[\[Kappa]])))/(
p (a^2 +
p) ((1 + E^((2 Sqrt[a^2 + p])/
Sqrt[\[Kappa]])) p \[Gamma] + (-1 + E^((2 Sqrt[a^2 + p])/
Sqrt[\[Kappa]])) Sqrt[a^2 + p] Sqrt[\[Kappa]])), p, t, 6]]*
BesselJ[0, a (r)]*a, {a, 0, \[Infinity]},
Method -> {"LevinRule", "LevinFunctions" -> {BesselJ}},
MaxRecursion -> 40];
with \[Gamma] = 0; I however compare the numerical results and analytical one. It seems that the the numerical solution obtain a wrong result. The code and the plot of comparison are shown as follows
\[Gamma] = 0; TT1 = Table[{t, s0[0.9, 0.5, t]}, {t, 0.1, 1, 0.1}];
TT2 = Table[{t, s0[0.9, 0.5, t]}, {t, 1, 10, 1}];
TT3 = Table[{t, s0[0.9, 0.5, t]}, {t, 10, 100, 10}];
TT4 = Table[{t, s0[0.9, 0.5, t]}, {t, 100, 1000, 100}];
C0 = Join[TT1, TT2, TT3, TT4];
TT1 = Table[{t, sol[1, 0.5, t]}, {t, 0.1, 1, 0.1}];
TT2 = Table[{t, sol[1, 0.5, t]}, {t, 1, 10, 1}];
TT3 = Table[{t, sol[1, 0.5, t]}, {t, 10, 100, 10}];
TT4 = Table[{t, sol[1, 0.5, t]}, {t, 100, 1000, 100}];
Cn = Join[TT1, TT2, TT3, TT4];
The figure is u versus t and blue and yellow lines represent the value predicted by numerical and analytical solutions, respectively. I an trying to figure out what happen to the numerical method.
Comments
Post a Comment