1
dC1/dt = -(Fp)(L)/Vp*(dC1/dx) - PSg/Vp(C1/Wp-C2) + Dp*(dC1/dx2);

dC2/dt = PSg/Visfp(C1/Wp-C2) + Disf*(dC2/dx2);

My boundary conditions are:

C1(t,x)|x=0 = Cin(t)

C2(t,x)|x=0 = 0

Cin(t) is a function of a known form e.g. Cin(t)= exp(-t);

dC1(t,x)/dx|x=0 = 0

dC2(t,x)/dx|x=0 = 0

dC1(t,x)/dx|x=L = 0

dC2(t,x)/dx|x=L = 0

The initial conditions are:

C1(t,x)|t=0 = 0

C2(t,x)|t=0 = 0

moreover, the initial condition is piece-wise:

when (t=t.min)  {
C1 = if (x=x.min) 1 else 0
C2 = 0;

I prepared a code in mathematica:

Fp = 0.016666667; Vp = 0.07; Wp = 0.94; Dp = 
 10^-18; PSg = 0.05; Visfp = 0.35; Disf = 10^-18; L = 0.1;

\[Phi][x_] := Piecewise[{{1, x == 0}, {0, x > 0}}]; 
mol[n_, o_: "Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

{sol4, solp} = 
 NDSolve[{Derivative[1, 0][C1][t, 
       x] == -Fp*L/Vp*Derivative[0, 1][C1][t, x] - 
       PSg/Vp*(C1[t, x]/Wp - C2[t, x]) + 
       Dp*Derivative[0, 2][C1][t, x], 
     Derivative[1, 0][C2][t, x] == 
      PSg/Visfp*(C1[t, x]/Wp - C2[t, x]) + 
       Disf*Derivative[0, 2][C2][t, x], C1[t, 0] == Exp[-t], 
     C2[t, 0] == 0, C1[t, L] == 0, C2[t, L] == 0, 
     Derivative[0, 1][C1][t, 0] == 0, Derivative[0, 1][C2][t, 0] == 0,
      Derivative[0, 1][C1][t, L] == 0, 
     Derivative[0, 1][C2][t, L] == 0, C1[0, x] == \[Phi][x], 
     C2[0, x] == 0}, {C1, C2}, {t, 0, 30}, {x, 0, L}, 
    Method -> mol[##]] & @@@ {{33, 4}, {33}}

Plot3D[Evaluate[C1[t, x] /. sol4[[1]]], {t, 0, 30}, {x, 0, L}]

However, no matter what the grid size the solution is wrong. It should look like (C1[x,t]):

1

The problem starts when I use a piece-wise initial condition. if I set C1[0, x] ==1 then the solution is ok (compared to the other software I use). But this piece-wise initial condition is crucial (I compare the results with other software that I know for sure gives correct results)

I'd appreciate any feedback.

Many thanks,

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
user2640106
  • 11
  • 2
  • You might address more people with an insight into this type of problem in the mathematica and scientific computing communities of stackexchange. – Lutz Lehmann Apr 06 '18 at 14:35
  • If scrape and paste your MMA code and I approximate your `Phi` with `Phi[x_]:=E^-(x)` or `Phi[x_]:=E^-(50 x)` then I get "Warning: Boundary and Initial conditions are inconsistent." Do you see that error? Can you correct that explain and correct that error? – Bill Apr 06 '18 at 16:34
  • the inconsistent condition warning is not necessarily a problem. Its not unusual for a well posed problem to have an local inconsistency at the corner. – agentp Apr 06 '18 at 17:51
  • Many thanks for all of your feedback. The problem can be avoided by approximating the stepwise I.C. with C1[0, x] ==Exp[-1000*x] – user2640106 Apr 07 '18 at 18:19
  • To avoid "Warning: Boundary and Initial conditions are inconsistent." you can comment the boundary condtions at x=0 .i.e. remove Derivative[0, 1][C1][t, 0] == 0, Derivative[0, 1][C2][t, 0] == 0 – user2640106 Apr 07 '18 at 18:29

0 Answers0