# 2-step diffusion equation with Neumann boundary conditions

I want to resolve a PDE model of 2-step 1D heat diffusion equation.

1st step: If 0< t < 60, it has a constant heat flax at x = 0, and has a Neumann boundary conditions at x = 6000.

2nd step: If t >= 60, it as a Neumann boundary conditions at both edge (x = 0, and x = 6000).

The key problem is that I have some trouble when I use the 1st step results as 2nd step initial value. Any suggestions how to fix it? Consider the following code:

[1st step]

Needs["NDSolveFEM"]

SoD = 100;
Dif = 69340;

deqN = D[u[t, x], t] - Dif*D[u[t, x], {x, 2}] - If[x == 0 , SoD, 0] ==
NeumannValue[-SoD, x == 0] + NeumannValue[0, x == 6000];

ic = u[0, x] == 0;

sol = NDSolveValue[{deqN, ic}, u, {t, 0, 60}, {x, 0, 6000}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.1}}}}];
Plot3D[sol[t, x], {t, 0, 60}, {x, 0, 6000}, PlotRange -> Full,
PlotStyle -> Automatic] I sucessfully got a 1st step result, but I cannot use this result in following 2nd step.

[2nd step]

deqN2 = D[u2[t, x], t] - Dif*D[u2[t, x], {x, 2}] ==
NeumannValue[0, x == 0] + NeumannValue[0, x == h];
ic2 = u2[0, x] == Evaluate[sol[t, x] /. t -> 0];
sol2 = NDSolveValue[{deqN2, ic2}, u2, {t, 0, 100}, {x, 0, 60000}, Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.1}}}}];


You don't need to separate the problem in two steps :

SoD = 100;
Dif = 69340;

deqN = D[u[t, x], t] - Dif*D[u[t, x], {x, 2}]  ==
NeumannValue[If[t<60 ,100, 0], x == 0] + NeumannValue[0, x == 6000];

ic = u[0, x] == 0;

sol = NDSolveValue[{deqN, ic}, u, {t, 0, 160}, {x, 0, 6000},
Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> {"Length" -> 0.1}}}}];

Plot3D[sol[t, x], {t, 0, 160}, {x, 0, 6000}, PlotRange -> Full, PlotStyle -> Automatic] Note

NeumannValue[If[t<60 ,SoD, 0], x == 0] doesn't work, though it should be equivalent to NeumannValue[If[t<60 ,100, 0], x == 0], so far I know.

EDIT

If you want to conserve SoD in the expression If[t<60 ,SoD, 0], you can do :

With[{SoD=SoD},D[u[t, x], t] - Dif*D[u[t, x], {x, 2}] == NeumannValue[If[t<60 ,SoD, 0], x == 0] + NeumannValue[0, x == 6000]]

• (+1) While this is what OP requested, I'd like to point out that the hard jump If[t<60 ,100, 0] for t=60 is not 'natural' - such a jump is not observed in nature. One could fade out the 100 to 0 starting at t=60 and have it decay at some rate. Just a thought. – user21 Apr 5 '17 at 8:28
• @user21 Physically, it's easy to make relatively hard jumps with a (electrical) resistor. – andre314 Apr 5 '17 at 13:06
• Yes, that's what I am trying to say. A more accurate model would use a relative hard jump, not an absolute jump. Maybe I'll add a message to NDSolve that warns about this. It's a bit like inconsistent initial and boundary conditions. – user21 Apr 5 '17 at 13:25