1
$\begingroup$

The formulation of the problem:

enter image description here

I tried to solve it with MOL and the method of V.R. Voller:

enter image description here

enter image description here

The script:

Clear["Global`*"]
n = 100; (*Number of Elements*)
H = 0.05; (*Length*)
dH = H/n; (*Length of an Element*)

\[Rho]s = 917;
\[Rho]l = 1/(1.003*^-3);
cs = 0.185*^3;
cl = 4.179*^3;
ks = 1.16*1.91;
kl = 613*^-3;

\[Alpha]l = kl/(\[Rho]l*cl);
\[Alpha]s = ks/(\[Rho]s*cs);

L = 333;

q = 40000;
\[Theta]i = 271;
\[Theta]m = 273;
Z = Table[Subscript[z, i], {i, 0, n}];
For[i = 0, i <= n, i++, Subscript[z, i] = i*dH];

tmax = 10;

(*Defining Tables *)
\[CapitalTheta][t_] = Table[Subscript[\[Theta], i][t], {i, 0, n}];
\[CapitalAlpha][t_] = Table[Subscript[\[Alpha], i][t], {i, 0, n}];
F[t_] = Table[HeavisideTheta[Subscript[\[Theta], i][t]], {i, 0, n}];

dd\[CapitalTheta] = 
  NDSolve`FiniteDifferenceDerivative[2, 
   dH*Range[0, n], \[CapitalTheta][t]];
d\[CapitalTheta] = 
  NDSolve`FiniteDifferenceDerivative[1, 
   dH*Range[0, n], \[CapitalTheta][t]];

dvbls = Head[#] & /@ \[CapitalAlpha][t];

(*eqn*)
\[Theta]eqn = 
  Thread[D[\[CapitalTheta][t], t] - \[CapitalAlpha][t]*
      dd\[CapitalTheta] + \[Rho]l*L*D[F[t], t] == 0];

(*BC*)
bc0 = d\[CapitalTheta][[1]] == -\[Rho]l*cl*q/kl;
bcH = Subscript [\[Theta], n][t] == \[Rho]s*cs*(\[Theta]i - \[Theta]m);

\[Theta]eqn[[1]] = bc0;
\[Theta]eqn[[-1]] = bcH;

(*IC*)
\[Theta]ic = 
  Thread[\[CapitalTheta][0] == cs*\[Rho]s*(\[Theta]i - \[Theta]m)];
\[Alpha]ic = Thread[\[CapitalAlpha][0] == \[Alpha]s];

(*Events*)
event1 = WhenEvent[Subscript[\[Theta], #][t] <= 0, 
     Subscript[\[Alpha], #][t] -> \[Alpha]s] & /@ Range[0, n];
event2 = WhenEvent[Subscript[\[Theta], #][t] > 0, 
     Subscript[\[Alpha], #][t] -> \[Alpha]l] & /@ Range[0, n];

system = Join[\[Theta]eqn, \[Theta]ic, \[Alpha]ic, event1, event2];

lines = NDSolveValue[system, 
  Join[\[CapitalTheta][t], \[CapitalAlpha][t]], {t, 0, tmax}, 
  DiscreteVariables -> dvbls, 
  StepMonitor :> (laststep = thisstep; thisstep = t; 
    stepsize = thisstep - laststep;)]

If I try to solve it without WhenEvent - assuming constant diffusing coefficient, I get excellent results. But if I include the events, I get a lot of errors:

enter image description here

I didn't find any good documentation to those errors. It would be very helpful if someone had a reference to that kind of error, or knows where my code went wrong.

Thanks!

$\endgroup$
3
  • 1
    $\begingroup$ Did you seen this post mathematica.stackexchange.com/questions/184920/… ? $\endgroup$ Mar 15, 2021 at 20:44
  • $\begingroup$ @AlexTrounev I had - it's a really good one so thanks. But I also try to understand why my script does not work, as it applies a different method. $\endgroup$
    – Ymir
    Mar 16, 2021 at 6:25
  • $\begingroup$ You code is fine, but it not for numerical computation. You can compare it with Tim Laska approach to compute solidification problem with enthalpy method. Did you seen this paper emerald.com/insight/publication/issn/0961-5539 $\endgroup$ Mar 16, 2021 at 15:20

1 Answer 1

0
$\begingroup$

With a little modification this code works also with the message

NDSolveValue::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.

Nevertheless we have output and it looks reasonable for alpha

Clear["Global`*"]
n = 50;(*Number of Elements*)H = 0.05;(*Length*)dH = 
 1/n;(*Length of an Element*)\[Rho]s = 917;
\[Rho]l = 1/(1.003*^-3);
cs = 0.185*^3;
cl = 4.179*^3;
ks = 1.16*1.91;
kl = 613*^-3;

\[Alpha]l = kl/(\[Rho]l*cl);
\[Alpha]s = ks/(\[Rho]s*cs);

L = 333;

q = 40000;
\[Theta]i = 271;
\[Theta]m = 273;
Z = Table[Subscript[z, i], {i, 0, n}];
For[i = 0, i <= n, i++, Subscript[z, i] = i*dH];

tmax = 10;

(*Defining Tables*)
\[CapitalTheta][t_] = Table[Subscript[\[Theta], i][t], {i, 0, n}];
\[CapitalAlpha][t_] = Table[Subscript[\[Alpha], i][t], {i, 0, n}];
F[t_] = Table[
   ArcTan[10^2 Subscript[\[Theta], i][t]]/Pi + 1/2, {i, 0, n}];

dd\[CapitalTheta] = 
  NDSolve`FiniteDifferenceDerivative[2, 
   dH*Range[0, n], \[CapitalTheta][t]];
d\[CapitalTheta] = 
  NDSolve`FiniteDifferenceDerivative[1, 
   dH*Range[0, n], \[CapitalTheta][t]];

dvbls = Head[#] & /@ \[CapitalAlpha][t];

(*eqn*)
\[Theta]eqn = 
  Thread[D[\[CapitalTheta][t], t] - \[CapitalAlpha][t]*
      dd\[CapitalTheta]/H^2 + \[Rho]l*L*D[F[t], t] == 0];

(*BC*)
bc0 = d\[CapitalTheta][[1]] == -H \[Rho]l*cl*q/kl;
bcH = Subscript[\[Theta], n][t] == \[Rho]s*cs*(\[Theta]i - \[Theta]m);

\[Theta]eqn[[1]] = bc0;
\[Theta]eqn[[-1]] = bcH;

(*IC*)
\[Theta]ic = 
  Thread[\[CapitalTheta][0] == cs*\[Rho]s*(\[Theta]i - \[Theta]m)];
\[Alpha]ic = Thread[\[CapitalAlpha][0] == \[Alpha]s];

(*Events*)

event1 = WhenEvent[Subscript[\[Theta], #][t] <= 0, 
     Subscript[\[Alpha], #][t] -> \[Alpha]s] & /@ Range[0, n];
event2 = WhenEvent[Subscript[\[Theta], #][t] > 0, 
     Subscript[\[Alpha], #][t] -> \[Alpha]l] & /@ Range[0, n];


system = Join[\[Theta]eqn, \[Theta]ic, \[Alpha]ic, event1, event2];


lines = NDSolveValue[system, 
   Join[\[CapitalTheta][t], \[CapitalAlpha][t]], {t, 0, tmax}, 
   DiscreteVariables -> dvbls];

Note that we map solution on the unit interval and use continues function for F. Visualization

{ListLinePlot[
  Table[Drop[lines, -n - 1][[1 ;; 20]]/(\[Rho]s*cs), {t, 1, 10, 1}], 
  PlotRange -> All, PlotLegends -> Automatic], 
 ListLinePlot[Table[Drop[lines, -n - 1]/(\[Rho]s*cs), {t, 1, 10, 1}], 
  PlotRange -> {-10, 5}, Frame -> True, PlotLabel -> T], 
 ListLinePlot[Table[Drop[lines, n + 1], {t, 1, 10, 1}], 
  PlotRange -> All, PlotLabel -> "\[Alpha]"]}

Figure 1

$\endgroup$

Your Answer

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

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