Heat Equation with Mobile Boundary

Question

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

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.

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

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