How to specify a NeumannValue in a system of non-linear partial differential equations?

Question

I am trying to solve a nonlinear system of partial differential equations. If I assume DirichletConditions it all works well.

tmax = 100;
zmax = 1;
pDGL0=D[c[z, t], t] ==Sqrt[T[z, t]] D[c[z, t], z, z] - 100 D[c[z, t],z];
pDGL1 = D[T[z, t], t] == 0.001 D[T[z, t], z, z] - 0.01 D[T[z, t], z];
aW0 = c[z, 0] == If[0 < z <= zmax, 0, 1];
aW1 = T[z, 0] == If[0 < z < zmax, 300, 400];
rW0 = DirichletCondition[c[z, t] == 0, z == zmax]; 
rW1 = DirichletCondition[c[z, t] == 1, z == 0] ;
rW2 = DirichletCondition[T[z, t] == 400, z == zmax];
rW3 = DirichletCondition[T[z, t] == 400, z == 0];
dgln = {pDGL0, pDGL1, aW0, aW1, rW0, rW1, rW2, rW3};
sol = NDSolve[dgln, {c, T}, {t, 0, tmax}, {z, 0, zmax}, MaxSteps-> {50,Infinity}]

This leads to a solution like:

cLsg[z_, t_] := c[z, t] /. sol[[1]]; 
TLsg[z_, t_] := T[z, t] /. sol[[1]];
Plot3D[cLsg[z, t], {z, 0, zmax}, {t, 0, tmax}]
Plot3D[TLsg[z, t] - 273.15, {z, 0, zmax}, {t, 0, tmax}]

However when I try to set a NeumannValue

Derivative[1, 0][c][0, t] == 1;

at one end of the boundary namely z==0 I get an error message like

CoefficientArrays::poly: c 23536+100 c23537-c23538Sqrt[T]-NeumannValue[0.2,z==0] is not a polynomial.

and also

NDSolve::femper: PDE parsing error of {c23536+100\ c23537-c23538\ Sqrt[T]-NeumannValue[0.2,z==0],T23539+0.1\ T23540-0.1\ T23541}. Inconsistent equation dimensions.

The code creating the error message is:

tmax = 2;
zmax = 1;
pDGL0 = D[c[z, t], t]==Sqrt[T[z, t]] D[c[z, t], z, z]-100 D[c[z,t],z]+NeumannValue[1, z == 0];
pDGL1 =  D[T[z, t], t] == 0.1 D[T[z, t], z, z] - 0.1 D[T[z, t], z];
aW0 = c[z, 0] == If[0 < z <= zmax, 0, 1];
aW1 = T[z, 0] == If[0 < z < zmax, 300, 400];
rW0 = DirichletCondition[c[z, t] == 0, z == zmax] ;
rW1 = DirichletCondition[T[z, t] == 400, z == zmax];
rW2 = DirichletCondition[T[z, t] == 400, z == 0];
dgln = {pDGL0, pDGL1, aW0, aW1, rW0, rW1, rW2};
sol = NDSolve[dgln, {c, T}, {t, 0, tmax}, {z, 0, zmax}]

The problem might be that there is a coupled system of non-linear pdes. When I tried to give the NeumannValue explicitly by

Derivative[1, 0][c][0, t] == 1;

this often gave the error message that it is consistant with its inital value or it said that the number of boundary conditions is insufficient. How can I fix this problem? Is there a better way to formulate the NeumannValues?

Answer 1

You may do it by MethodOfLines. Try this:

 tmax = 100;
zmax = 1;
pDGL0 = D[c[z, t], t] == 
   Sqrt[T[z, t]] D[c[z, t], z, z] - 100 D[c[z, t], z];
pDGL1 = D[T[z, t], t] == 0.001 D[T[z, t], z, z] - 0.01 D[T[z, t], z];
aW0 = c[z, 0] == If[0 < z <= zmax, 0, 1];
aW1 = T[z, 0] == If[0 < z < zmax, 300, 400];
rW0 = c[zmax, t] == 0;
rW1 = c[0, t] == 1;
rW2 = T[zmax, t] == 400;
rW3 = T[0, t] == 400;
nV = (D[c[z, t], z] /. z -> 0) == 1;
dgln = {pDGL0, pDGL1, aW0, aW1, rW0, rW1, rW2, rW3, nV};
sol = NDSolve[dgln, {c, T}, {t, 0, tmax}, {z, 0, zmax}, 
  Method -> "MethodOfLines"]

There will be several warnings, but the solution, nevertheless, looks reasonably:

    cLsg[z_, t_] := c[z, t] /. sol[[1]];
TLsg[z_, t_] := T[z, t] /. sol[[1]];

Row[{Plot3D[cLsg[z, t], {z, 0, zmax}, {t, 0, tmax}, 
   AxesLabel -> {"z", "t", "c      "}, ImageSize -> 250], Spacer[20],
  Plot3D[TLsg[z, t], {z, 0, zmax}, {t, 0, tmax}, 
   AxesLabel -> {"z", "t", "T   "}, PlotRange -> All, 
   ImageSize -> 250]}]

giving this:

I mean that the boundary conditions are fulfilled in the solution. You may then further try to improve the MethodOfLines by methods you may find in this tutorial.