How to discretize a PDE with a Moving Boundary? (a.k.a Stefan Problem, a.k.a Variable Boundary)

[Jax_Xu]

Hi Community! I need your expertise ! I am solving a 1D diffusion equation with a moving boundary. The physical background is the reaction of gas with a solid particle, where the gas diffuses through the product layer but cannot penetrate the dense core so that surface reaction happens. The dense core will be gradually consumed and shrink. In combustion science, this classical model is named shrinking core model.
I tried to explain this model clearly with a diagram and my equations, as shown in the supplementary material. What is tricky here, for me, is the left boundary of the domain – it changes with time. I checked the previous topics in our Julia community, and they do not really help.

My codes are also attached for your perusal, which is problematic (I know~). Feel free to share your ideas!

using MethodOfLines, OrdinaryDiffEq, Plots, DifferentialEquations, Symbolics, Latexify, DomainSets, ModelingToolkit

r0 = 0.01; K = 2e-4; De = 6e-6; mol_rho = 109500
Nx=30
order=2

C=20; k=0.08; x_inf=1 #bulk concentration; reaciton rate const; molar fraction of bulk
@parameters r, t
@variables u(…)
@variables r_c(…)

Dr = Differential(r)
Drr = Differential(r)^2
Dt = Differential(t)

eqs=[Dt(u(r,t)) ~ De*Drr(u(r,t)),
Dt(r_c(t)) ~ (-k)Cu(r_c(t),t)/mol_rho]

bcs=[Dr(u(r_c(t),t)) ~ -ku(r_c(t),t)/r_c(t),
DeDr(u(r0,t)) ~ -K*(u(r0,t)-u_inf/r0),
u(r,0) ~ 0.0, r_c(0) ~ r0]

domains = [r ∈ Interval(0, r0), t ∈ Interval(0.0, 30.0)]

@named scmsys = PDESystem(eqs,bcs,domains,[r,t], [u(r,t), r_c(t)])

discretization = MOLFiniteDifference([r => Nx], t)

@time prob = discretize(scmsys, discretization)

I think the domain should be [0, r_0] or [r_c, r_0]. The former domain discretizes the whole area, and returns zero values to u(r,t) when r<r_c. The latter is a dynamic discretization which I have no idea how to achieve that. Or could this be an interface problem ([0, r_0] + [r_c, r_0])?

Similar cases were found in a Mathematica community here and there. I am trying to understand but really struggling.

Supplementary figures

(Note: the diffusion equation is altered for simplification, which is inaccurate in fact but does not influence the solution of this topic.)

[vettert]

I worked extensively with a PhD student on a similar problem in a different context (spray drying; it’s an evaporating (=shrinking) droplet in which particles precipitate at some point).

Our solution is here: A mechanistic model to predict droplet drying history and particle shell formation in multicomponent systems - ScienceDirect
(In case you don’t have access, pm me; I am sure I can find the preprint somewhere!).

The problem can essentially be solved by changing the coordinate system into one with a fixed boundary; this introduces a virtual reaction term instead of the moving boundary. The problem than becomes a standard one for which many solution approaches work.

If you read the paper, skip the part about particle formation; it’s unimportant in your context.

[Jax_Xu]

Thanks, Thomas! Such a good paper! I am reading already !

[Jax_Xu]

Hi Vettert! Thank you for your sharing! I have successfully discretised the moving boundary by normalising the moving domain into a relatively fixed domain and solve the PDE! It took me some time to understand the Eq.(21) of your paper, which is the transfer of time differentiation. Therefore, I checked your citation (here, for anyone may find this topic helpful), which gave more details.

Again, I am amazed by how smart human can be, especially scientists! Thank you!

(the original publication came up with the normalisation to deal with the moving boundary might be referred to Crank, John. “Free and moving boundary problems.” (No Title) (1984).)