Question: Heat transfer with phase change and non homogeneous Dirichlet boundary conditions



I'm currently working on the modelling of a thermodynamic process.

Briefly, I cool down a solution (water + polymer) from -5°C to -15°C to induce a phase separation. At the end (and after removing of the water by lyophilisation) I obtain a porous sponge like material.

The process uses a home made cooling system which can be described like this:

- A Peltier module

- An aluminium layer recovered by teflon (And also a layer of ethanol)

- The mold (in stainless steel) in which I pour the (polymer+water) solution

- A layer of air.


The bottom part of the Peltier is kept at -15°C during the all experiment. Then the top part Peltier (x=0) heat to -5°C the full system until thermal equilibrium. At t0, the Peltier lets the temperature in x=0 goes down linearly to -15°C in a controlled way by the adjustement of the cooling rate. A thermocouple measure the temperature in the aluminium part and at ~-10°C (in the aluminium) I observe an increase of the measured temperature due to the phase change of the solution.

I try to model the full process to determine how long is necessary to reach -15°C in the (polymer+water) solution. At the moment I already made it but without the introduction of the phase change of the solution which can be considered as a heat source. Normally I should obtain a decreasing surface with a kind of 2D gaussian in x= polymer solution... But It doesn't work...


At the moment, I use this method:

Definition of the variables

Cooling rate of the Peltier module (in °C.min) :  R1 := .5;  (As T(0,t)=T0-R1*time)
Cooling rate of the air in the cooling chamber (in °C.min) :  R2 := .95*R1; Thid should represent the convection function which is not represented here (It's obviously an approximation) as I don't know how to work with real convection... I need it as boundary conditions.
Initial time of cooling (in min): t0 := 20; (After 20 min I switch on the system and it begins to cool down)
Initial temperature of the system (in °C) : T0 := -5;
Final temperature of the system (in °C) : Tf := -15;
Cooling time in x=0 (in min): tau1 := (T0-Tf)/R1;
Cooling time of air, x= top of the full system (in min): tau2 := (T0-Tf)/R2;
Time of experiment (in min): texp := t0+tau1+50; I wait ~ 50minutes at the end, to reach the equilibrium but this value is completely arbitrariness.

Description of the geometry

1. Thickness (in m):

Aluminium x[1] := 0.6e-2;
Teflon film x[2] := 0.1e-3;
Ethanol x[3] := 0.1e-3;
Stainless steel x[4] := 0.2e-2;
Polymer solution x[5] := 0.2e-2;
Air x[6] := .1;

Size of the system (in m): Ltot := sum(x[n], n = 1 .. 6);
Region of interest (in m): Lint := sum(x[n], n = 1 .. 5);

2. Values of the thermal diffusivity in function of x

Thermal diffusivity (in m2.s-1)
Aluminium alpha1 := 1/(8.418*10^(-5));
Teflon film alpha2 := 1/(.124*10^(-6));
Ethanol alpha3 := 1/(7*10^(-8));
Stainless steel alpha4 := 1/(3.428*10^(-6));
polymer solution (Assimilated to water) alpha5 := 1/(.143*10^(-6));
Air alpha6 := 1/(1.9*10^(-5));

Functions for boundary conditions

1. Thermal diffusivity function (Considered as a constant on this temperature range)

alpha := (x) -> piecewise(0 < x and x

2. Boundary conditions in x = 0
Ti := (t)-> piecewise(0 < t and t

3. Boundary conditions in x = Ltot (The air cools down slowly but we choose to introduce it as a boundary condition)
Tl := (t) -> piecewise(0 < t and t

Solving of the heat equation

1. Initial and boundary conditions

IBC := {T(0, t) = Ti(t), T(x, 0) = T0, T(Ltot, t) = Tl(t)};

2. Numeric solving of the heat equation

PDE := diff(T(x, t), x, x)-alpha(x)*(diff(T(x, t), t)) = 0;

pds := pdsolve(PDE, IBC, numeric);

3. Plotting domain

P1 := pds:-plot3d(x = 0 .. Lint, t = 0 .. texp);

4. Plot

plots[display]({P1}, axes = boxed);

This gives a nice estimation of the temperature evolution in the system, I can change a lot of parameters for optimisation and I see directly the temperature evolution in the system according to the cooling rate in x=0.


HOWEVER, I don't know how to introduce the local temperature increase during the phase change, I tried to add a third boundary condition to solve it by introducing a 2D gaussian curve but the software refuses this and it's not really the best I think. The best would be to add a source term in the heat equation but I don't know how... SOMEBODY TO HELP ME ?


I don't know if I'm clear enough but it would be really nice to help me on this topic which is, of course, fascinating... No ?


Thanks in advance,



Please Wait...