Problems when solving a nonlinear PDE system with NDSolve

The nonlinear PDE system is actually extracted from a research paper published in 2000 Here is the paper link

The authors solved the system by using an ordinary differential equation integrator in FortranVariable-coefficient Ordinary Differential Equation solver, which suggests this is based on the "Method of Lines" available in Mathematica's NDSolve for transient nonlinear PDE systems.

The Following Codes have been updated

However, when I am trying to put the system into NDSolve, different error messages occur:

(*
Dear challenger:

If you've tried to solve this set of equations but failed,
please increment the following counter:

totalHoursSacrificedHere = 8 +6/60
*)
ClearAll["Global`*"]
(*Predefined Constants*)
ho=-57*10^5;hp=57*10^4;ha=-25*10^6;

Here I convert the units into SI, so these are different from original authors'.

no1=410;nc1=210;nc2=240;no3=1650;na3=30;(*converted from grams to kilograms *)
ng1=1000+no1-nc1;ng2=1000-nc2;ng3=1000+no3-na3;(*Derived constants do not appear in the table *)

(*Try Leach's original data settings *)
(*ho=ho/1000;hp=hp/1000;ha=ha/1000;*)
no1=no1/1000;nc1=nc1/1000;nc2=nc2/1000;no3=no3/1000;na3=na3/1000;
ng1=ng1/1000;ng2=ng2/1000;ng3=ng3/1000;

I think there are typo's in the original article. The unit of pre-exponential factors should be $s^{-1}$, not $L/s$. So converting from $L/s$ into $m^3/s$ is not necessary. This can be verified by the unit of the same constants used in table 2 of the same article. Such constants were determined by the cited reference  of the article, which further verifies this issue.

Ao=569*10^9;Ap=2*10^17;Aa=5*10^8;
EoR=19245;EpR=26500;EaR=19244;
f=13/10;g=18/10;h=78/100;m=1/2;
kc=42/1000;kf=63/1000;kg=258/10000;
ρc=10;ρf=53/2;
Cc=1100;Cf=1700;
df=5*10^-5;dc=13/10000;
ϕ=98/100;
Do2m=453*10^-7;
(*Constants for boundary and initial conditions*)
Tw=600;Tgin=300;Tin=300;Yo2in=23/100;uin=53/10000;
xl=0;xr=15/100;

(*Constants not appeared apparently in literature of Leach et al*)
σ=56704*10^-12;R=83144/10000;
patm=101325;M=28956*10^-6;Pr=7/10;
A=8000/5;cg=1005;μg=2*10^-5;ν=144*10^-7;

(*Functions defined*)
ωo[x_, t_] := (1 - yc[x, t] - ya[x, t])^f *ρf*Ao*(Ybo[x, t])^m *Exp[-EoR/T[x, t]]
ωp[x_, t_] := (1 - yc[x, t] - ya[x, t])^g *ρf*Ap*Exp[-EpR/T[x, t]]
ωa[x_, t_] :=yc[x, t]*ρc*Aa*(Ybo[x, t])^h*Exp[-EaR/T[x, t]]
ρg[x_, t_] := patm*M/R/Tg[x, t]

The $V_{O_2}$ has been given in the article and already defined similarly as a function. There do be a series of physical constants, e.g., $\sigma$ Stefan-Boltzmann constant, universal gas constant $R$, Prantl constant $Pr\approx 1 \text{ or } 0.7$, which should be determined based on the chemical engineering background knowledge and I have included them into the equations and definitions implicitly without further clarification.

Vo2[x_, t_] := -Do2m*D[Ybo[x, t], x]
hmA3p[x_,t_] := (4*10^-5/df + 42/100*ug[x, t]*(2*10^-5)^(2/3)*ν^(-1/6)) A;
hA3p[x_, t_] := A*(2 + 11/10*(7/10)^(1/3)*((df*ρg[x, t]*ug[x, t])/μg)^(3/5))*kg/df;

(*Nonlinear PDEs*)
eq1 = ((yc[x, t] + ya[x, t])*ρc + ρf*(1 - yc[x, t] -ya[x, t]))*((yc[x, t] + ya[x, t])*Cc + Cf*(1 - yc[x, t] - ya[x, t]))*D[T[x, t], t] == D[((yc[x, t] + ya[x, t])*kc +kf*(1 - yc[x, t] - ya[x, t]) + (16*σ*df*T[x, t]^3)/3)*
D[T[x, t], x], x] - ωo[x, t]*ho - ωp[x, t]* hp - ωa[x, t]*ha + hA3p[x, t] (Tg[x, t] - T[x, t]);

eq2 = 10*D[yc[x, t], t] ==nc1* ωo[x, t] + nc2* ωp[x, t] - ωa[x, t];
eq3 = 10*D[ya[x, t], t] == na3* ωa[x, t];

eq4 = ϕ*ρg[x, t]*cg*(D[Tg[x, t], t] + ug[x, t]*D[Tg[x, t], x]) == ϕ*kg*D[Tg[x, t], {x, 2}] +
hA3p[x, t] (T[x, t] - Tg[x, t]);

In the same article, all the porosity ratios $\phi'$s with different subscripts have been considered as approximately the same constant 0.98. So I followed this convention. The ϕ in equation 5 can be completely cancelled actually.

Equations 6 and 7 have been combined together with $y^s_{O_2}$ eliminated. Also by confirmation with cited , the $y_{O_2}$ in the initial and boundary conditions section and in equations 2 and 3, should be $y^b_{O_2}$, i.e., oxygen concentration in the bulk of the gas.

eq5 = ϕ (D[ρg[x, t], t] + ϕ*D[ρg[x, t]*ug[x, t], x]) == ϕ (ng2 *ωp[x,t] + (ng1 - no1) ωo[x, t] + (ng3 - no3)*ωa[x, t]);
eq6 = ϕ (D[ρg[x, t]*Ybo[x, t], t] + D[ρg[x, t]*Ybo[x, t]*ug[x, t], x] + D[ρg[x, t] Ybo[x, t] Vo2[x, t], x]) == -(no1* ωo[x, t] + no3* ωa[x, t]);
sys = {eq1, eq2, eq3, eq4, eq5, eq6};

tr = 2000; trs = Min[1800, tr];
(*Calculate ya[xr=15/100,t] and yc[xr=15/100,t] boundary conditions*)

convertRule = {T[x, t] -> Tw, Ybo[x, t] -> Yo2in,
Derivative[0, 1][yc][x, t] -> fyc'[t], yc[x, t] -> fyc[t],
Derivative[0, 1][ya][x, t] -> fya'[t], ya[x, t] -> fya[t]};
ode1 = eq2 /. convertRule // Simplify;
ode2 = eq3 /. convertRule // Simplify;

sol = NDSolve[{ode1, ode2, fya == 0, fyc == 0}, {fya, fyc}, {t,
0, tr}, Compiled -> False, PrecisionGoal -> 20, WorkingPrecision -> 100];
{Plot[fyc[t] /. sol, {t, 0, trs}, ImageSize -> 300,
PlotPoints -> Floor[trs]] Plot[fya[t] /. sol, {t, 0, trs},
ImageSize -> 300, PlotPoints -> Floor[trs]]}

(*Initial and boundary conditions*)
ibcs={Derivative[1,0][T][xl,t]==0,T[xr,t]==Tw,T[x,0]==(UnitStep[x-xr]+1)Tgin,Derivative[1,0][Tg][xl,t]==0,Tg[xr,t]==Tgin,Tg[x,0]==Tgin,
Derivative[1,0][Ybo][xl,t]==0,Ybo[xr,t]==Yo2in,Ybo[x,0]==Yo2in,yc[x,0]==0,yc[xr,t]==fyc[t]/.(sol[[1,2]]),ya[x,0]==0,ya[xr,t]==fya[t]/.(sol[[1,1]]),(*Derivative[1,0][ug][0,t]\[Equal]0,*)ug[x,0]==uin,ug[xr,t]==uin};

fsol = NDSolve[sys~Join~ibcs, {Tg, T, yc, ya, Ybo, ug}, {x, xl, xr}, {t, 0, tr},Method -> {"PDEDiscretization" -> {"MethodOfLines","SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 50}}}, MaxSteps -> 10^3, AccuracyGoal -> 8(*,MaxStepFraction\[Rule]1/100*)(*,(*Compiled\[Rule]False,*)WorkingPrecision\[Rule]50*)
]

During evaluation of In:= NDSolve::pdord: Some of the functions have zero differential order, so the equations will be solved as a system of differential-algebraic equations. >>

Does this mean the Method of Lines cannot be used?

During evaluation of In:= NDSolve::bcart: Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable x. Artificial boundary effects may be present in the solution. >>

For Method of lines algorithm, such initial and boundary conditions should already be enough.

During evaluation of In:= NDSolve::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. >>

Out= {}

When I change the options Compiled->False and WorkingPrecision->50, muc more memory is needed than my PC has.

Can such precision issue be solved by just increasing the working precison? or I will have to use some NDSolve Options which are dedicated to such very stiff cases?

Using IndexReduction method gives strange results do not match the real cases. (It seems this method can only be immediately used for differential-algebraic equations with ordinary differential equations only. For those cases with "epllitic-parabolic" PDEs, tutorial indicates difference matrices should be constructed first. and I have another question posted on the tutorial actually)

fsol = NDSolve[sys~Join~ibcs, {Tg, T, yc, ya, Ybo, ug}, {x, xl, xr}, {t, 0, tr},Method -> {"PDEDiscretization" -> {"MethodOfLines","SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 50}}}, MaxSteps -> 10^3, AccuracyGoal -> 8(*,MaxStepFraction\[Rule]1/100*)(*,(*Compiled\[Rule]False,*)WorkingPrecision\[Rule]50*)
]

Further Update

I have tried to rewrite the code per the tutorial mentioned here, first construct difference matrices; then the difference equations with boundary and initial conditions; but still obtain warning message as below, and the results are different from those obtained by Leach et al(2000) article though I already added one more initial condition $y_a(x,0)=0$ than the 2000 article:

NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended. >>

• I suggest you to place those parameters and equations in the order they appeared in the paper, currently it's really frustrating to check if there's any simple mistake. – xzczd Jan 21 '16 at 5:59
• You can consider uploading it to a net-disk. – xzczd Jan 21 '16 at 11:25
• As mentioned above, you'd better not to change the order of parameters, just place them in the order they appeared in the paper so it's easier to check if they're correct. Also, it's not a good idea to include those value of parameters (101325*29 etc.) in the equations directly, it's a potential source of mistake. – xzczd Jan 22 '16 at 5:59
• One way is simply to copy the symbols from here or elsewhere. Halirutin's program can be used for more extensive copying of special characters. I too have a package, and I intend to write it up soon. – bbgodfrey Jan 30 '16 at 4:16
• It would be best if Wolfram can add one more option to the "Copy As" menu command to solve it immediately. I look forward to your package. Once finished, do let me know. thank you. – LCFactorization Jan 30 '16 at 4:33

I guess the culprit behind is the not-strong-enough DAE solver. Using the 3rd approach in this answer i.e. by descretizing the spatial derivatives manually, I managed to solve your equation set without warning, but the results seems to be different from those in the paper, I'm not sure about the reason. (For such a complicated problem there're too many possibilities! ) Anyway, I'd like to show you my solution.

Equations and parameters are the same as what you gave in the question so they're omitted in this answer.

pdetoode is a general purpose function for discretizing PDE to ODE, its definition can be found here.

First, some preparation:

lb = xl; rb = xr;

xdifforder = 4;

points = 50;
grid = Array[# &, points, {lb, rb}];
var = {Tg, T, yc, ya, Ybo, ug};

ptoo = pdetoode[var[x, t], t, grid, xdifforder];

odevar = # /@ grid & /@ var;

OK, let's discretize ibcs:

odeibc = ibcs // ptoo;

In the next part of code, sys is also discretized, then the obtained equations at the "endpoints" that are duplicates of b.c.s are removed. This is the most tricky part of the solution, because equations in sys are coupled so it's not that easy to identify the redundant ODEs. There should be cleverer ways, but here I simply found the right ones to remove by trial and error:

eqngroup = {#, Complement[sys, #]} & /@ Subsets[sys, {3}];
i = 8;(*Found by trial and error*)
odeq1 = (eqngroup[[i, 1]] // ptoo)[[All, 2 ;; -2]];
odeq2 = (eqngroup[[i, 2]] // ptoo)[[All, 1 ;; -2]];
odesol = NDSolveValue[{odeq1, odeq2, odeibc}, odevar, {t, 0, tr}];

sollst = rebuild[#, grid, -1] & /@ odesol

Partition[Plot3D[#[x, t], {t, 0, tr}, {x, lb, rb}, PlotRange -> All] & /@ sollst,
3] // Grid • thank you very much. I got similar result as yours. I contacted the corresponding author with my questions and results half a year ago but got no response yet. – LCFactorization Oct 8 '16 at 9:04
• Can you please post a ready- to- work version of your code for verification so that I can try it easily? – LCFactorization Oct 8 '16 at 12:40
• @LCFactorization Just execute those code lines before fsol = NDSolve[sys~Join~ibcs,…… in your question first. – xzczd Oct 8 '16 at 12:56
• The solution this way is easy for me to understand. I will spend time in understanding you code. How long did you spend in solving the PDE system this way I am just curious? – LCFactorization Oct 8 '16 at 17:08
• @LCFactorization The spirit of your solution is the same as mine, the only difference is I've automate the discretization with pdetoode etc.. I've added some explanation for pdetoode in my answer, have a look. As to the timing for solving your system, I should say I'm not sure, but the upper limit is 4 hours. – xzczd Oct 9 '16 at 2:25