Applying NDSolve to a multi-layer problem (three coupled PDEs)

Question

I am trying to solve the following governing equations via NDSolve.

The system consists of three layers. The following shows how I constructed a numerical model my system with given parameter values:

rw = 0.01;
er = 0.5;
ez = er;
D2 = 1;
Dp = 0.1;
bp = 1;
b1 = 5;
b2 = 5;
R1 = 1;
R2 = 1;
Rp = 1;
λ1 = 0;
λ2 = λ1;
λp = λ1;
Q = 1;
θ1 = 0.5;
θ2 = 0.5;
θp = 0.5;
A = Q/(2*Pi*θ1*b1);
C0 = 10;
vp = 0.1;
v2 = 0.1;
R = 10;
NR[z_] := If[z <= b1, R1, If[b1 < z < b1 + bp, Rp, R2]];
Ner[z_] := If[z <= b1, er, Exp[-10000 z]];
Nez[z_] := If[z <= b1, ez, If[b1 < z < b1 + bp, (Dp*vp)/A, (D2*v2)/A]];
Nv[z_] := If[z <= b1, Exp[-10000 z], If[b1 < z < b1 + bp, vp, v2]]; 
NA[z_] := If[z <= b1, A, Exp[-10000 z]]; 
sys = 
  {NR[z]*Derivative[0, 0, 1][NC][r, z, t] == 
    (A*Ner[z]*Derivative[2, 0, 0][NC][r, z, t] + 
      (A*Ner[z] + NA[z])/r*Derivative[1, 0, 0][NC][r, z, t]) + 
      A*Nez[z]*Derivative[0, 2, 0][NC][r, z, t] - 
      Nv[z]*Derivative[0, 1, 0][NC][r, z, t],
  NC[rw, z, t] == C0,
  NC[R, z, t] == Exp[-10000 z],
  Derivative[0, 1, 0][NC][r, 0, t] == Exp[-10000 r],
  Derivative[0, 1, 0][NC][r, b1 + bp + b2, t] == Exp[-10000 r],
  NC[r, z, 0] == Exp[-10000 r]}; 
sol = 
  NDSolve[
    sys, NC, {r, rw, R}, {z, 0, b1 + bp + b2}, {t, 0, 100}, 
    Method -> "BDF", AccuracyGoal -> 24, PrecisionGoal -> 10]

The coefficients are a function of z; thus, the governing equations are suitable for each governing equation of the three layers.

To validate the numerical results, I derived a semi-analytical solution as a reference. The code producing the reference is longer than I am permitted to post, so I can not include its code here.

The numerical results seem to be unstable and give the totally difference from reference solution.

Show[
  Plot[Evaluate[NC[r, 2.5, 100] /. sol], {r, rw, R}, PlotRange -> All], 
  ListPlot[RR, PlotRange -> All]]

The solid line is the numerical result and the dotted line is the reference solution,

Can you see any problems with the way I am handling the numerical model for the three layers? Do you have any suggestion about dealing with this kind problem?

Answer 1

It's because you've made at least one simple mistake when coding the equation. Check the A Ner[z] + NA[z] part carefully. It's a very bad idea to change the form of the equation when coding it, you should make the equations in the code and in the original text as similar as possible.

A trivial issue is those Exp[10000 …] doesn't make sense in your code, if you're trying to omit the ibcinc warning, you still need to adjust their coefficients carefully. In your case the inconsistency doesn't cause trouble, so let it be.

Another trivial issue is the adjustion for AccuracyGoal and PrecisionGoal doesn't make sense here. For more informatin, check this post.

Here's my attempt for coding your equation:

pw[part1_, part2_, part3_] = 
  Piecewise[{{part1, 0 < z < b1}, {part2, b1 < z < b1 + bp}}, part3];
coedt[z_] = pw[R1, Rp, R2];
coedr2[z_] = pw[er A, 0, 0];
coedr[z_] = pw[(er A) (1 - 1/er)/r, 0, 0];
coedz2[z_] = pw[(er A) ez/er, Dp, D2];
coedz[z_] = pw[0, vp, v2];

sys = With[{NC = NC[r, z, t]},
           {coedt[z] D[NC, t] == coedr2[z] D[NC, r, r] + coedr[z] D[NC, r] + 
                                 coedz2[z] D[NC, z, z] - coedz[z] D[NC, z],
            NC == C0 /. r -> rw, NC == 0 /. r -> R, 
            D[NC, z] == 0 /. z -> 0, D[NC, z] == 0 /. z -> b1 + bp + b2, 
            NC == 0 /. t -> 0}];

sol = NDSolve[sys, NC, {r, rw, R}, {z, 0, b1 + bp + b2}, {t, 0, 100}]; // AbsoluteTiming

Plot[Evaluate[NC[r, 2.5, 100] /. sol], {r, rw, R}, PlotRange -> All]