8
$\begingroup$

I came across the paper Solidification dynamics of an impacted drop regarding a heat equations by Thiévenaz et.al and was interested in knowing how they obtained the graphs presented. From what I understand, they solved the PDE presented below, with a moving boundary using the finite differences method. I did solve them in matlab but was wondering how one would do this in mathematica given that the flexibility of matlab with matrices is much greater. The code in matlab may be found here. The relevant equations are: The two heat equations for the temperature field $T(z,t)$ $$\frac{\partial T}{\partial t}=D_s\frac{\partial^2 T}{\partial z^2}, \,\,z\leq0;\qquad \frac{\partial T}{\partial t}=D_i\frac{\partial^2 T}{\partial z^2}, \,\,0\leq z\leq h(t)$$ and the four boundary conditions:

$$T(0^-,t)=T(0^+,t);\quad \lambda_s\frac{\partial T(0^-,t)}{\partial t}=\lambda_i\frac{\partial T(0^+,t)}{\partial t};\quad T(h(t)^-,t)=T_m;$$ $$\lambda_i\frac{\partial T(h(t)^-,t)}{\partial t}=\rho_i L\frac{dh}{dt}$$ where $L$ denotes the latent heat of solidification, $\rho_k$ the thermal density, $\lambda_n$ the thermal conductivity, $D_k=\lambda_k/(\rho_kC_k)$ the diffusion coefficient. Here $k$ denotes $l$ in the liquid phase, $i$ in the solid phase and $s$ for the substrate (they are all constant parameters). It is also imposed a constant temperature $T_s$ deep in the substrate, i.e., $$\lim_{z\to-\infty}T(z,t)=T_s,\quad T(z,0)=T_s\quad z\leq0, h(0)=0$$ which means that at time $0$ the liquid is put into contact with the substrate.

A pictorial description of the model is showed below: enter image description here

Note: The constants mentioned above in my code are chosen arbitrarily (i.e. I chose the values assigned to them)

I did not post the code because it was somewhat long.


As requested in the comments, the scheme I used was:

$$\vec{T^{\,n+1}}=A^{-1}\left(\vec{T^{\,n}}+\vec{B}\right)\implies c\rho\frac{T_i^{n+1}-T_i^n}{\Delta t}=k\frac{T_{i-1}^{n+1}-2T_i^{n+1}+T_{i+1}^n}{\Delta x^2}$$ which is the implicit version. The explicit is similar, and looks like $$c\rho\frac{T_i^{n+1}-T_i^n}{\Delta t}=k\frac{T_{i-1}^{n}-2T_i^{n}+T_{i+1}^n}{\Delta x^2}$$ In my code you will see how I found $A$ and $B$. ($A$ is a matrix and $B$ a vector)

$\endgroup$
  • 1
    $\begingroup$ I'd say Mathematica and MATLAB are about the same when it comes to 95% of the major matrix ops, but maybe you mean element-wise manipulation? Also I think it'd be nice to at least post some code detailing the method used. It's a big ask for us to read, translate, and do all the mental processing without any sort of starting point. $\endgroup$ – b3m2a1 Dec 10 '19 at 6:03
  • $\begingroup$ @b3m2a1 Yes what I meant was element-wise manipulation. And it appears to me it feels more natural in matlab when doing so, or maybe I’m not experienced enough in mathematica. Regarding the method, as I mentioned above, I used the finite difference method which there are a lot of snippets around the web. I will try to post some code, maybe the general idea. Anyways, any help is appreciated. It’s just a curiosity $\endgroup$ – DMH16 Dec 10 '19 at 6:06
  • 2
    $\begingroup$ I'd say most element-wise tasks can be done better as vectorized operations and once you do that there's really no difference between Mathematica and MATLAB. That's really neither here nor there, though. Looks like you just need to use the FDM, from what I glanced at in that image. Lots of questions here pertain to that but maybe xczd will come by and help you. Also how about as a compromise for not having some preliminary code, you post the relevant equations here in LaTeX form? $\endgroup$ – b3m2a1 Dec 10 '19 at 6:09
  • 3
    $\begingroup$ Finally, we already have a few posts about moving b.c., have you read them?: 1. mathematica.stackexchange.com/q/58596/1871 2. mathematica.stackexchange.com/q/130962/1871 3. mathematica.stackexchange.com/q/95358/1871 4. mathematica.stackexchange.com/q/114967/1871 5. mathematica.stackexchange.com/q/184920/1871 $\endgroup$ – xzczd Dec 10 '19 at 6:35
  • 1
    $\begingroup$ The question is much better after the edits. (Still, you need to show us the parameters and i.c.s. ) I'd say the MATLAB code isn't quite related to this question in my view. Is the finite difference scheme used by the paper necessary? If so, it'll be better to add the difference scheme to the question directly, because it's easier for us to understand the difference formula in traditional math notation. If not, then I believe the combination of DChange, pdetoode and NDSolve is a better choice for solving the problem, and again, we don't need to know how the MATLAB code looks like. $\endgroup$ – xzczd Dec 10 '19 at 6:58
6
+500
$\begingroup$

As pointed out in the comment, OP's MATLAB implementation doesn't seem to be correct, so I'll just give a solution to OP's problem.

Since the heat flux continuity involves in, there's indeed something new in this problem compared to previous ones. We first handle the continuity condition. Notice that, strictly speaking, 1D heat conduction equation is not

$$\frac{\partial T}{\partial t}=\frac{\lambda}{\rho C}\frac{\partial^2 T}{\partial z^2}$$

but

$${\rho C}\frac{\partial T}{\partial t}=\frac{\partial}{\partial z}\left(\lambda\frac{\partial T}{\partial z}\right)$$

These 2 equations are equivalent only if $\lambda$ is constant, and piecewise constant is not constant. "So what? Why you mention this? " Because the latter equation owns a good property: when it's discretized directly, the heat continuity condition will be automatically satisfied.

Then we need to deal with the moving b.c.. The "change the variable -> discretize" method can be used in principle, but to make this answer simpler and more interesting, I'd like to try something I've never tried: I'll directly discretize the PDE in $z$ direction on a moving grid. This can be easily done with the help of pdetoode.

Given that OP hasn't shown the parameters, they'll be casually chosen as

ρL = 1; eps = 10^-2;
inf = 2; Tm = 1; Ts = 2; tend = 10;

SetAttributes[#, Listable] & /@ {ρc, λ, icfunc};
ρc[a_ /; ! FreeQ[a, h]] := 2
ρc[a_ /; a <= 0] := 1
λ[a_ /; ! FreeQ[a, h]] := 3
λ[a_ /; a <= 0] := 1
icfunc[a_ /; ! FreeQ[a, h]] := Ts
icfunc[a_ /; a <= 0] := Tm

$\rho C$ (ρc), $\lambda$ (λ) and initial condition (icfunc) are defined in a seemingly strange way, so that e.g. ρc[{-1/10, h[t]/10}] will evaluate to {2, 1}. This will make subsequent programming more convenient.

Next we interpret the PDE, i.c. and b.c. to Mathematica code:

With[{T = T[z, t], dT = dT[z, t]},
 eq = ρc[z] D[T, t] == D[dT, z];
 eqd = dT == λ[z] D[T, z];
 bc = {{λ[z] D[T, z] == ρL h'[t], T == Ts} /. z -> h[t], 
   T == Tm /. z -> -inf};
 icextra = h[0] == eps;
 ic = T == icfunc[z] /. t -> 0;
 ]

Create the moving mesh and discretize:

domain@1 = {-inf, 0};
domain@2 = {0, h[t]};
points = 400;
difforder = 2;
gridfunc = Array[# &, points, #] &;
grid = gridfunc@domain@1~Join~Rest@gridfunc@domain@2;

(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = pdetoode[{T, dT}[z, t], t, grid, difforder];
del = #[[2 ;; -2]] &;
ode = Block[{dT}, Set @@ ptoofunc /@ eqd;
   del@ptoofunc@eq];
odeic = ptoofunc@ic;
odebc = ptoofunc@bc;

The resulting system is a DAE system and NDSolve can't handle it well, so we create a toode function to transform the system to an ODE system (For more info check this and this post) and solve:

toode = With[{sf = 10}, sf # + D[#, t] & /@ # &];
{funch, funcTlist} = 
  NDSolveValue[
   MapAt[toode, {ode, icextra, odeic, odebc} /. T[i_] :> T@ToString@i // 
     Flatten, {{-1}, {-2}}], {h, T@ToString@# & /@ grid}, {t, 0, tend}, 
   MaxSteps -> Infinity];

funcT[1] = rebuild[funcTlist[[;; points]], grid[[;; points]], 2];
funcT[2] = rebuild[funcTlist[[-points ;;]], grid[[-points ;;]]/h[t], 2];
func = Function[{z, t}, 
   Piecewise[{{funcT[1][z, t], z <= 0}, {Ts, z >= funch[t]}}, funcT[2][z/funch[t], t]]];

Finally, just some illustration:

ListLinePlot[funch, PlotRange -> All]

enter image description here

Here funch is plotted using an undocumented syntax.

Manipulate[Plot[func[z, t], {z, -inf, funch[tend]}, PlotRange -> {0, 2}, 
  AxesLabel -> {"z", "T"}], {t, 0, tend}]

enter image description here

Animate[DensityPlot[func[z, t], {x, 0, 1}, {z, -inf, funch[tend]}, PlotRange -> {0, 2}, 
  PlotPoints -> 50, PerformanceGoal -> "Quality", ColorFunction -> "DeepSeaColors"], {t, 
  0., tend, 0.05}]

enter image description here

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ The process is very clear. I will have to look up the definitions to understand the code. Anyways, thank you for putting the time $\endgroup$ – DMH16 Dec 13 '19 at 23:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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