3
$\begingroup$

I wish to solve is the heat equation with solution-dependent coefficient.
The equation along with BC's and IC's are as follows:
$\begin{equation} \frac{\partial P[x,t]}{\partial t}-\alpha[x,t]\frac{\partial^2 P[x,t]}{\partial x^2}=0 \end{equation}$,
$\begin{equation}P[x,0]=0,P[L,t]=0 \end{equation}$,
$P[0,t]=Step[t-t1]-0.5Step[t-t2]-0.25Step[t-t3]$,
$L=1$.
The $\alpha[x,t]$ coeffiecnt is determined by the following rules:
$\alpha[x,0]=\alpha_1$,
If $P[x,t]>0.7$ $\rightarrow$ Set $\alpha[x,t]$ to $\alpha_2$
If $P[x,t]<0.3$ $\rightarrow$ Set $\alpha[x,t]$ to $\alpha_1$
If $0.3<P[x,t]<0.7$ $\rightarrow$ Do Nothing.

Those conditions represent the case of bistability of the diffusion coefficient.

I've tried quite a few methods to solve it, including the definition of $\alpha$ as a function of $\frac{\partial P[x,t]}{\partial t}$ and $P[x,t]$ with a combination of Sing[],Boole[],WhenEvent[],ect.

None of them worked at the past few weeks, so any idea may help.
Thanks,
Ofek.

Improve this question
$\endgroup$
5
  • $\begingroup$ If the thermal conductivity is variable, don't you need to express the heat equation as $\rho {{\hat C}_p}(u)\frac{{\partial u}}{{\partial t}} + \nabla \cdot \left( { - k(u)\nabla u} \right) = 0$? $\endgroup$
    – Tim Laska
    Jan 5, 2019 at 13:49
  • $\begingroup$ Yes, the case I've mentioned is a simplified version of the full equation (as you wrote) $\endgroup$
    – Ofek Peretz
    Jan 5, 2019 at 14:14
  • $\begingroup$ In the last application of the WhenEvent Documentation, they have an example of an ON/OFF controller (which is a bistable system with hysteresis) applied to the 2D heat equation. Did you try to adapt that example? $\endgroup$
    – Tim Laska
    Jan 5, 2019 at 15:59
  • $\begingroup$ @OfekPeretz Can you provide data for a numerical model? $\endgroup$
    – Alex Trounev
    Jan 5, 2019 at 16:26
  • 6
    $\begingroup$ It is effectively impossible to diagnose what has gone wrong if the actual code used is not posted. $\endgroup$
    – Daniel Lichtblau
    Jan 5, 2019 at 16:26

2 Answers 2

Reset to default
9
$\begingroup$

I think the OP is trying to demonstrate hysteresis. WhenEvent has an example of an on-off controller, which is an example of a bistable system. Unfortunately, my understanding is, that an event must be a function of time only and cannot have local spatial dependence. The essence of the Numerical Method of Lines is to turn a time-dependent PDE into a system of ODEs depending on time only. This suggests that perhaps we can combine WhenEvent with a roll your own version of the Numerical Method Of Lines. To build this model, it is best to start simple so let us start with a zero dimensional case that demonstrates hysteresis with WhenEvent.

0-D Case

Let us consider a small system with no diffusion limitations driven by a fluctuating heat source. The ODE that we wish to solve is given by:

$${\hat C_p}(u)\frac{du}{dt} = q\left( t \right)$$

We will start with a system that is initially solid. We will set up a simple test case to verify that our WhenEvent construction gives us the correct bistable behavior.

cpl = 2;
cps = 1;
thresholdliquid = Rationalize[0.7];
thresholdsolid = Rationalize[0.3];
q[t] = 5/2 UnitStep[-7 + t] - 5/2 UnitStep[-6 + t] + 
   5/2 UnitStep[-5 + t] - 5/2 UnitStep[-3 + t] + UnitStep[-1 + t];
NDSolve[{cp[t] u'[t] == q[t], cp[0] == cps, u[0] == 0, 
   WhenEvent[{u[t] == thresholdliquid, 
     u[t] == thresholdsolid} , {cp[t] -> 
      If[u[t] == thresholdliquid, cpl, 
       If[u[t] == thresholdsolid, cps, cp[t]]]}]}, {u, cp}, {t, 0, 
   10}, DiscreteVariables -> {cp \[Element] {cps, cpl}}  ];
Plot[{Evaluate[{u[t], cp[t] - cps} /. %], thresholdsolid, 
  thresholdliquid}, {t, 0, 10}, PlotTheme -> "Web", 
 PlotStyle -> {Thick, Thick, Dashed, Dashed}]

Bistable 0D case

You can clearly see breaks in the temperature curve that have different thresholds if the temperature is rising or falling. So, it looks as though we have a good prototype for bistability using WhenEvent in the 0D case.

Combining WhenEvent with a Roll-Your-Own Numerical Method of Lines

In the Boundary Conditions section of The Numerical Method Of Lines Tutorial, Shows how you can roll your own numerical method of lines schemes to convert a PDE into a system of ODEs. The following Mathematica code attempts to combine our WhenEvent prototype with the procedure outlined in the numerical method of lines tutorial. I used arbitrary values for the parameters since none were supplied. The thermal diffusivity difference between ice and water is about a factor of 9. However, an ice water system would have such a large latent heat effects that you would need to use the form of the equation I alluded to in the comments.

(* Set Parameters *)
n = 50; Subscript[h, n] = 
 1/n; sf = 1; tmax = 10; alphasolid = 1; alphaliquid = 1/10*alphasolid;
t1 = 1; t2 = 4; t3 = 5 ;
thresholdliquid = Rationalize[0.7];
thresholdsolid = Rationalize[0.3];
(* Create Temperature and Phase State Vector *)
U[t_] = Table[Subscript[u, i][t], {i, 0, n}];
A[t_] = Table[Subscript[a, i][t], {i, 0, n}];
dvbls = Head[#] & /@ A[t];
(* Differentiate BCs *)
bc0 = Subscript[u, 0][t] == 
   UnitStep[t - t1] - 1/2*UnitStep[t - t2] - 1/4*UnitStep[t - t3];
bc0d = Map[D[#, t] + sf # &, bc0];
bc1 = Subscript[u, n][t] == 0;
bc1d = Map[D[#, t] + sf # &, bc1];
(* Set Laplacian *)
laplacianu = 
  NDSolve`FiniteDifferenceDerivative[2, Subscript[h, n] Range[0, n], 
   U[t]];
(* Set up system of ODEs *)
eqns = Thread[D[U[t], t] -  A[t] laplacianu == 0];
(* Replace Boundary Condtions Equations *)
eqns[[1]] = bc0d;
eqns[[-1]] = bc1d;
(* Set Initial Temperature to 0 *)
ics = Thread[U[0] == Table[0, {n + 1}]];
(* Set Initial Phase State to Solid *)
istates = Thread[A[0] == Table[alphasolid, {n + 1}]];
(* Set up Events*)
events = WhenEvent[{Subscript[u, #][t] == tl, 
       Subscript[u, #][t] == ts} , {Subscript[a, #][t] -> 
        If[Subscript[u, #][t] == tl, alphaliquid, 
         If[Subscript[u, #][t] == ts, alphasolid, 
          Subscript[a, #][t]]]}] & /@ 
    Range[0, n] /. {ts -> thresholdsolid, tl -> thresholdliquid};
(* Join System of Equations to be Solved *)
system = Join[eqns, ics, istates, events];
lines = NDSolve[system, Join[U[t], A[t]], {t, 0, tmax}, 
   DiscreteVariables -> dvbls];
(* Create Interpolation Functions for Easier Plotting *)
f = Interpolation[
   Flatten[Table[{{t, i Subscript[h, n]}, 
      First[Subscript[u, i][t] /. lines]}, {i, 0, n}, {t, 0, tmax, 
      tmax/500}], 1], InterpolationOrder -> 2];
fa = Interpolation[
   Flatten[Table[{{t, i Subscript[h, n]}, 
      First[Subscript[a, i][t] /. lines]}, {i, 0, n}, {t, 0, tmax, 
      tmax/500}], 1], InterpolationOrder -> 2];
(* Plot Temperature and Diffusivity *)
ContourPlot[f[t, x], {t, 0.`, tmax}, {x, 0.`, 1}, 
 PlotPoints -> {200, 200}, MaxRecursion -> 4, PlotTheme -> "Web", 
 ColorFunction -> ColorData["DarkBands"]]
ContourPlot[fa[t, x], {t, 0.`, tmax}, {x, 0.`, 1}, 
 PlotPoints -> {200, 200}, MaxRecursion -> 4, PlotTheme -> "Web", 
 ColorFunction -> ColorData["DarkBands"]]

Temperature Contour Plot

Bistable Temperature Plot

Diffusivity State Contour Plot

Bistable Diffusivity Plot

Demonstration of Bistability

We should be able to demonstrate bistability by plotting the temperature and the diffusivity at a value that intersects the phase envelope (e.g., $x=0.1$) as shown with the following Mathematica code.

Plot[{f[t, 0.1], fa[t, 0.1], thresholdsolid, thresholdliquid}, {t, 0, 
  10}, PlotTheme -> "Web", 
 PlotStyle -> {Thick, Thick, Dashed, Dashed}]

Bistable MOL

The system behaves as expected showing transitions at 0.7 going up and 0.3 going down. One can use ParametricPlot to demonstrate the path dependence more clearly as shown below.

p0 = ParametricPlot[{f[t, 0.1], fa[t, 0.1]}, {t, 0, 10}, 
   ColorFunction -> (ColorData["Rainbow"][#3] &), MaxRecursion -> 0, 
   PlotPoints -> 100, PlotTheme -> "Web"];
p0 /. l_Line :> {Arrowheads[ConstantArray[0.05, 30]], Arrow[l]}

Hysteresis Plot

Update To Include No Flux Condition at x=1

In the comments, Ofek requested to also consider the case of a no flux condition at $x$=$1$. It turns out that the MOL Tutorial actually considers this very case. One just needs to adjust bc1 and bc1d accordingly

bc1 = Last[
    NDSolve`FiniteDifferenceDerivative[1, Subscript[h, n] Range[0, n],
      U[t]]] == 0;
bc1d = Map[D[#, t] &, bc1];

Note that the scale factor $sf$ has been set to zero for Neumann type conditions. Without adjustments to the first case, you will obtain plots for temperature and diffusivity state as shown below:

Temperature Contour Plot No Flux on End

Temperature No Flux Case

Diffusivity State Contour Plot No Flux on End

Diffusivity No Flux Case

You will have to adjust parameters to make the plots look more interesting.

Improve this answer
$\endgroup$
8
  • 1
    $\begingroup$ a very good example in MMA! $\endgroup$
    – ABCDEMMM
    Jan 7, 2019 at 7:13
  • $\begingroup$ @TimLaska, Thanks for your detailed solution, this is exactly what I wanted to solve. Just one more question - How can I change the boundary condition at $X=1$ from $P_n(0)==0$ into $\frac{\partial P_n(0)}{\partial x}==0$? $\endgroup$
    – Ofek Peretz
    Jan 14, 2019 at 12:29
  • 1
    $\begingroup$ @OfekPeretz I added a no flux case to the post. $\endgroup$
    – Tim Laska
    Jan 14, 2019 at 13:46
  • $\begingroup$ @TimLaska, Amazing. Thanks! $\endgroup$
    – Ofek Peretz
    Jan 14, 2019 at 15:12
  • $\begingroup$ @OfekPeretz, You are welcome! $\endgroup$
    – Tim Laska
    Jan 14, 2019 at 15:23
7
$\begingroup$

I do not see what the problem is.

t1 = 1; t2 = 2; t3 = 3; L = 10; T = 10;
bc = {P[t, 0] == 0, 
   P[t, L] == 
    UnitStep[t - t1] - UnitStep[t - t2]/2 - UnitStep[t - t3]/4};
a[p_] := Piecewise[{{1, p < .3}, {2, .3 <= p < .5}, {1/2, True}}]
ic = {P[0, x] == 0};
eq = D[P[t, x], t] - a[P[t, x]]*D[P[t, x], x, x] == 0;
sol = NDSolveValue[{eq, ic, bc}, P, {t, 0, T}, {x, 0, L}];
{ContourPlot[sol[t, x], {t, 0, T}, {x, 0, L}, Contours -> 20, 
  ColorFunction -> Hue, PlotLegends -> Automatic, PlotRange -> All, 
  FrameLabel -> {"t", "x"}], 
 ContourPlot[sol[t, x], {t, 0, 5}, {x, L - 2, L}, Contours -> 20, 
  ColorFunction -> Hue, PlotLegends -> Automatic, PlotRange -> All, 
  FrameLabel -> {"t", "x"}]}

fig1

Improve this answer
$\endgroup$

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