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.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
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]]

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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