Solidification Heat Transfer Model using PDE Toolbox

Hello,

I would like to model the solidifcation of a material at initial temperature that is cooled down in outer space up to by radiation transfer on the surface. Due to the range of temperature, a phase change happens inside my domain. The material I use is basalt and I know from documentation those values :

• Specific heat of solid basalt:
• Latent heat of crystallization of basaltic magma:
• Crystallization temperature of basalt: ,

I have already implemented the following Finite Element Model but would now like to take into account the phase change.

%% Parameters
rho=1000;           % (kg/m**3) Density
cp=1000;            % (J/kg/K) Specific heat
T0=2000;            % (K) Initial temperature
T_out=300;          % (K) outer space temperature
eps=1;              % Emissivity
dt=200;             % (s) time-step
day=3600*24;
tmax=6*day;
tlist = [0:dt:tmax];
lambda = @(location,state) (0.46+0.95*exp(-2.3e-3*state.u)); 

%% Model 
thermalModel = createpde('thermal','transient');
gm = multisphere(1);
thermalModel.Geometry=gm;
generateMesh(thermalModel,'Hmax',0.2,"GeometricOrder","quadratic");
thermalModel.StefanBoltzmannConstant = 5.670373E-8;
thermalIC(thermalModel,T0);
thermalProperties(thermalModel,'ThermalConductivity',lambda,'MassDensity',rho,'SpecificHeat',cp);
thermalBC(thermalModel,"Face",1,"Emissivity",@(region,state) eps,"AmbientTemperature",T_out, "Vectorized","on");
thermalResults = solve(thermalModel,tlist);

Does anyone happen to know how to model the solidifcation inside the domain when it's cooled down over time ? Knowing that this phenomenon would start to occur at crystallization temperature of basalt. I thought about the internalHeatSource function but can't figure how to properly implement it.

Many thanks for any help you can give me !

Regards,

Tom

[Umar]

To model solidification with a phase change in Matlab, you can utilize the internalHeatSource function to account for the latent heat of crystallization. This function allows you to introduce a heat source term that represents the energy released or absorbed during the phase change process.

Here's a general outline of how you can implement the internalHeatSource function in your existing Finite Element Model:

Define the latent heat of crystallization value for basalt:

L = {latent heat value}; % Latent heat of crystallization of basaltic magma Modify your thermal model to include the internal heat source term:

internalHeatSource(thermalModel, L, 'Face', {face index where phase change occurs}); By incorporating the internalHeatSource function with the appropriate latent heat value and specifying the face where the phase change occurs, you can simulate the solidification process with the phase change effect in your model.

Feel free to adjust the parameters and functions based on your specific requirements and the details of your Finite Element Model. This approach should help you simulate the solidification process accurately while considering the phase change behavior of the material.

[Tom]

However, the heat source on a 3D model is on a cell. And I must take into account the crystallization temperature. I know I shoud get a cooling curve as follows :

[Yifeng Tang]

I was reading the documentation, and it appears that the specific heat can be made a function "that depends on space, time, or temperature." I've modeled phase-changing material in Simscape using custom equations, but the same idea may work here. The idea is to make cp a function of temperature, and around the melting point, cp will have a dramatic peak over a small T range (e.g. +/-1K). This is usually good enough for engineering purpose. The width, height and overall shape of the cp peak needs to ensure that the integral of (cp*dT) reflect the latent heat.

Below are the key equations in Simscape (not MATLAB)

T == M.T;
Q == mass * sp_heat_pcm * T.der;
sp_heat_pcm == sp_heat + h_trans/pi/(((T-T_trans)*R_trans)^2+1)*R_trans;

T_trans is the melting temperature, R_trans is the inverse of the temperature range (1/K), and h_trans is the latent heat.

I believe the integral will given you h_trans over the neighborhood of T_trans, but please verify before using.

[Tom]

Hi, thank you for your response. Once I modeled the peak for the specific heat as a function of temperature where the area reflect the latent heat, do I need to do anything else to take into account the phase change ?

Because now I have a function called in the material properties :

cp = @(location, state) ... ; 
thermalProperties(thermalModel,'ThermalConductivity',lambda,'MassDensity',rho,'SpecificHeat',cp);

Thank you

[Yifeng Tang]

If the material stays in place (doesn't flow away when in liquid phase) and the density change is neglible (so volume stays the same and you don't have to remap the coordinates), I feel like a cp function is sufficient to model the thermal aspects of the melting/solidification. You may be able to consider other properties, like conductivity, as functions of temperature, if you believe they do change significantly over the temperature range in question. I don't have specific domain knowledge here to recommend one way or the other.