1
$\begingroup$

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["NDSolve`FEM`"]

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]

enter image description here

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}}}}];

Thank you for your cooperation.

$\endgroup$
5
$\begingroup$

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]  

enter image description here

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]]

$\endgroup$
  • $\begingroup$ (+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. $\endgroup$ – user21 Apr 5 '17 at 8:28
  • 1
    $\begingroup$ @user21 Physically, it's easy to make relatively hard jumps with a (electrical) resistor. $\endgroup$ – andre314 Apr 5 '17 at 13:06
  • 1
    $\begingroup$ 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. $\endgroup$ – user21 Apr 5 '17 at 13:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.