I am trying to use NDSolve to find the solution to a set of coupled diffEQs. They represent 1d (radial) heat balance in a spherical shell.

The goal is to specify a heat flux into the base of the shell, Hb, and have the model output the shell thickness and temperature profile for which heat is balanced.

Within the shell, the temperature T[r,t] is guided by conduction, like

p C ∂T/∂t = ∇∙(A/T ∇T) + Q

where p is density, C specific heat capacity, and Q[r] a source term.
The conductivity is temperature dependent with k=A/T with A a constant.

At the base of the shell, any imbalance in flux into the base (Hb) and out of the base (Fb) is accommodated by phase change which changes the shell thickness, d[t], like

p L ∂d/∂t = -1/3 (Hb - Fb)

where p is density, and L is latent heat.
Flux in (Hb) is an input to the model, but the flux out of the base, Fb, is determined by the temperature gradient at the shell base

Fb = -A/Tb * T'[R-d[t],t]

where Tb = T[R-d[t],t] is the constant temperature at the shell base (rb[t] = R - d[t]).

When there is no source term, this is possible to do by hand so I know it's possible.
With a source, it gets wild though.

As an initial attempt I want to use NDSolve to give the solutions for T and d in the sourceless system (if i can figure it out without a source, I will add that in later).

Here is what I wrote...


Sorry its an image instead of a code block. I tried for some time using the formatting tool, but I use symbols and such that this editor does not seem to like and the code is nearly illegible when I tried. I'll post it here anyway for those who might want to try a copy/paste of it.

Dsolfunc\[LetterSpace]flux[{AI_, \[Rho]I_, cI_, LI_}, {rC_, R_}, {Tb_, Ts_}, Hb_, Qfun_, T0_, d0_] := NDSolve[{
    (* Conduction within shell *)
    \[Rho]I cI (\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\(T[r, t]\)\)) == Div[(AI/T[r, t]) Grad[T[r, t], {r, \[Theta], \[Phi]}, "Spherical"], {r, \[Theta], \[Phi]}, "Spherical"] + Qfun,
    (* Initial temperature *)
    T[r, 0]    == T0,
    (* Temperature at boundaries *)
    T[R, t]    == Ts,
    T[R - d[t], t]  == Tb,
    (* Phase change at base *)
    \[Rho]I LI (\!\(\*SubscriptBox[\(\[PartialD]\), \(t\)]\(d[t]\)\)) == -(1/3) (Hb + (AI/Tb) Derivative[1, 0][T][R - d[t], t]),
    (* Initial shell thickness *)
    d[0] == d0
    }, {T, d}, {r, rC, R}, {t, 0, 10^25}];

When I run the the NDSolve I get the error

enter image description here

I believe this has to do with the fact that NDSolve expects both T and d to be functions of radius, r, and time, t.
My shell thickness, however, is not a function of radius. It just changes in time.

So my question is - is there a way to use Mathematica to model diffEQs where one of the things you solve for is the position of the boundary?

The thickness solution needs the temperature gradient, but the temperature solution needs the thickness to apply the boundary condition, so its certainly coupled.

I viewed many other posts in the hopes I could scrounge together other people's work, but I am having difficulty parsing some of those. In particular, I thought this post an answer to related question was getting close to what I need, but it looks like both of the things they solve for are dependent on the same two variables still, even though it looks like a boundary is moving.

Alright if anyone has some insight I'd appreciate it!! Thanks much.

  • 1
    $\begingroup$ Here are instructions to post readable code from MMA: How to copy code from Mathematica so it looks good on this site. $\endgroup$
    – MarcoB
    Commented Jun 15, 2022 at 15:21
  • 1
    $\begingroup$ We already have quite a number of posts related to moving b.c. e.g. mathematica.stackexchange.com/q/211080/1871 (Don't miss the links in the comments. ) Have you read them? $\endgroup$
    – xzczd
    Commented Jun 15, 2022 at 15:43
  • $\begingroup$ I have spent the last two weeks of work dedicated only to looking through posts on this site to help. It's a problem Ive been trying to deal with for years, however, so I've come back to these posts in the hopes of finding something, on and off, for quite some time. However, I did not specifically find the one you mention xzczd, so I appreciate you sharing it I will take a look. $\endgroup$
    – Jeruwaho
    Commented Jun 15, 2022 at 15:49
  • $\begingroup$ Im at three weeks now staring at this problem only, and the links nested within the link that the user above provided. I still can not tell what about my problem is not solvable to Mathematica, especially considering it is a very easy problem by hand. My confusion is certainly related to being a novice at Mathematica. I see that xzczd wants me to use the functions they repeatedly promote in all these links, but I just can not tell what their function does, and how that helps me. Specifically the #, &, @ syntax boggles me. Ill just write my own integrator i guess, instead of NDSolve $\endgroup$
    – Jeruwaho
    Commented Jun 22, 2022 at 15:48

1 Answer 1


Since OP is having difficulty in understanding previous related posts, I'll solve the problem as a favor. But I won't explain much because there's nothing new in this problem compared with previous ones. The only thing I want to emphasize is, if you've really spent weeks on the problem but made no progress, you're probably in wrong way. Trying to read or even modify existing code without a basic understanding for the core language of Mathematica is not effective. To understand those seemingly horrible #, &, @, etc. you may start from this post. Books mentioned here are also good sources. Do remember #, &, @, etc. are unavoidable for any non-trivial work.

Definition of pdetoode and DChange are not included in this answer, please find them in the links.

{{AI, ρI, cI, LI}, {rC, R}, {Tb, Ts}, Hb, Qfun, T0, 
   d0} = {{567, 1000, 2000, 334000}, {1266/10, 1982/10} 10^3, {273, 80}, 10 10^-3, 0,
    273, 25000};

{eq, ic, bc} = {{ρI cI (D[T[r, t], t]) == 
     Div[(AI/T[r, t]) Grad[T[r, t], {r, θ, ϕ}, 
         "Spherical"], {r, θ, ϕ}, "Spherical"] + Qfun, ρI LI (D[d[t], 
        t]) == -(1/3) (Hb + (AI/Tb) Derivative[1, 0][T][R - d[t], t])}, {T[r, 0] == 
     T0, d[0] == d0}, {T[R, t] == Ts, T[R - d[t], t] == Tb}};

domain = {0, 1};
{neweq, newicmid, newbc} = 
 DChange[{eq, ic, bc}, ρ == Rescale[r, {R - d[t], R}, domain], r, ρ, T[r, t]]

newic = {newicmid[[1]] /. t -> 0 /. Rule @@ newicmid[[2]], newicmid[[2]]}

With[{eq = neweq, ic = newic, bc = newbc},
 points = 100; grid = Array[# &, points, domain]; difforder = 2;
 ptoofunc = pdetoode[T[ρ, t], t, grid, difforder];
 del = #[[2 ;; -2]] &;
 ode = {del@ptoofunc@eq[[1]], ptoofunc@eq[[2]]};
 odebc = ptoofunc@bc;
 odeic = {ic[[1]] // ptoofunc // del, ic[[2]]} /. t -> 0;
 tend = 10^25;]

{dsol, Tsollst} = 
   NDSolveValue[{ode, odeic, odebc}, {d, T /@ grid}, {t, 0, tend}]; // AbsoluteTiming

Tsolmid = rebuild[Tsollst, grid, 2]

tendreal = dsol["Domain"][[1, -1]]
(* 3.15985*10^16 *)

The solver stops half-way. This is reasonable, because d[t] == R at t == 3.15985*10^16:

Plot[dsol[t], {t, 0, tendreal}, GridLines -> {Automatic, {R}}]

enter image description here

 Plot[Tsolmid[(r - R)/dsol[t], t], {r, R - dsol[t], R}, 
  PlotRange -> {{R - dsol[tendreal], R}, {0, 300}}], {t, 0, 0.999 tendreal}]

enter image description here

  • $\begingroup$ I do appreciate your time for this problem. Very kind of you. It does look like it works! I have not had the chance to mess around with it and figure out how you got there, but it does seem relatively straightforward. I am a little surprised to see d[t] = R as time goes on because that means the shell thickens, while I know from the analytic solutions that it in fact thins. It might just be some flipped definitions I was not clear enough in my description, and I bet I can tease out. Ill look at this next week and reply with what I see (and mark this as answered as well.) $\endgroup$
    – Jeruwaho
    Commented Jun 29, 2022 at 15:06
  • $\begingroup$ And thank you for pointing to some documentation on the shorthand symbols. I anticipate trying to avoid those, in favor of explicit function notion and such, but perhaps that limits what I can do. Ill check that out too $\endgroup$
    – Jeruwaho
    Commented Jun 29, 2022 at 15:07
  • $\begingroup$ Well after a lot of playing with this, I think it does work in some circumstances! I'll give it the answered check, and thank you again. Whether it can work for the full problem I want it for, and not this simple example, remains to be seen. Specifically, I have a temperature dependent source term, which Ive not even attempted to re-include yet. I also noticed though, that the stuff does break down when I use a radially dependent initial temperature condition, as opposed to it being constant. Specifically T0 = Ts(Tb/Ts)^{(R-r)(R-d)/(r d)}. It works for a fixed constant initial temp tho $\endgroup$
    – Jeruwaho
    Commented Jul 20, 2022 at 16:35
  • 1
    $\begingroup$ @Jeruwaho I'm not sure how you code it, but modifying T0 and odeic to T0 = Ts (Tb/Ts)^((R - r) (R - d0)/(r d0)); and odeic = {ic[[1]] // ptoofunc // del, ic[[2]]} /. freeze /. d[t] -> d0; should work. $\endgroup$
    – xzczd
    Commented Jul 21, 2022 at 2:40
  • $\begingroup$ great thanks! ill give it a go. ill continue to post my updates back in here as new answers as I include more of the whole problem back in, if thats allowed, in case it helps some others with similar problems $\endgroup$
    – Jeruwaho
    Commented Jul 26, 2022 at 17:29

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.