Simscape model of latent heat storage in a PCM

Dear fellows,

I'm struggeling to model a latent heat storage component in Simscape for phase change materials (PCM). Hope you have ideas how to takle the problem.

Goal:

Simple PCM model that neglects the volume change, pressure increase or material flux. The only variables shall be the Simscape accross and through variables heat flow q in [J/s] and T in [K]. The parameters shall be:

• Cps... specific heat in solid state [J/(kg K)]
• Cpl... specific heat in liquid state [J/(kg K)]
• m... Mass of PCM material in [kg]
• Tmelt...Melting Temperature of PCM [K]
• L... Specific latent heat in [J/kg]

What I have:

In the change from solid to liquid, the storred energy increases strongly. For the concept, have a look at the PhD Thesis of Simone Arena.

The Formula for the storred heat energy in [J] is:

 Q(T=0 to Tmelt) = m * integral(Cps dT,0,Tmelt)
 Q(T=Tmelt) = m * L
 Q(T=Tmelt to Inf) = m * integral(Cpl dT, Tmelt, Inf)

How would you define the branches and equations for simscape for the heat flow q = dQ/dt?

Thank you so much!

[beta317]

The problem that you are describing is exactly the same problem that I am facing at the moment. Did you find a way to model a latent heat storage component in Simscape for phase change materials (PCM) and how did you accomplish that?

Thanks in advance for your help!

[Anubhav]

This code is not working, can anyone elaborate in simulink

Answer 1

[Andrew Schenk]

It seems this block would be very similar to the Simscape thermal mass block. You can see the Simscape source for this block by running this command:

>> edit foundation.thermal.elements.mass

In the thermal mass block, there is only one equation which defines Q. Based upon the equations you provided, you would modify the single heat flow equation to instead be:

if T < Tmelt
  Q == m * integral(Cps dT,0,Tmelt)
else if T > Tmelt
  Q == m * integral(Cpl dT, Tmelt, Inf)
else 
  Q == m*L;
end

It would be numerically preferable to reformulate the equations to use the Simscape der() operator instead of integ().

[Ralf Reisenauer]

Hi Andrew,

I tried to implement the code you suggested and for me it doesn't work. When I modify the thermal mass block I get the error of an unexpected Matlab operator in your line "Q == m * integral(Cps dT,0,Tmelt)". Do you have any idea why?

Many greetings,

Ralf

[Mohamed FRIKHA]

@Andrew Schenk , i guess it will not work, because Q is on KW and the last part of the equation is on KJ (Q=m*L) and also keep in mind that the temperatuer remain constant for a while in order to do the change phase,

i have a simular problem , i m trying to use the function " delay" i have not succeeded yet if any one have an idea , would you comment ?

[Vishnubharathi thirupathi]

Did you find the answer? How did you finally get to implement PCM block?

[Anubhav]

pls solve this in simulink because it is not working there

Answer 2

[Vaishak]

This code maybe useful:

component PCM_NEW
% Two-port thermal component
nodes
M = foundation.thermal.thermal; % :top
end
nodes(ExternalAccess = none)
N = foundation.thermal.thermal; % :bottom
end
parameters 
mass      = {1,   'kg'};    % Mass
end
parameters    
sp_heat   = {4186, 'J/(kg*K)'};                 % Specific heat
end
parameters    
la_heat   = {336000, 'J/kg'};                 % Latent heat
end
parameters    
T_melt   = {273, 'K'};                 % Melting Temperature
end
parameters (ExternalAccess = none)
diff      = {1,   'K'};    % Mass
end
parameters    
Intial_f   = {1, '1'};                 % Initial fraction of solid (0 to 1) 
end
parameters 
num_ports = foundation.enum.numPorts2.one;     % Number of graphical ports
end
if num_ports == 2
    annotations
    N : ExternalAccess=modify
end
end
variables
% Differential variables
T = {value = {300, 'K'}, priority = priority.high}; % Temperature
Q = {0, 'W'}; % Heat flow rate
f= {1, '1'}; % Initial fraction of solid (0 to 1)
end
branches
Q : M.Q -> *;
end
equations
assert(sp_heat > 0)
assert(la_heat > 0)
assert(T_melt > 0, 'Temperature must be greater than absolute zero')
assert(T > 0, 'Temperature must be greater than absolute zero')
T == M.T;
end
equations(Initial=true)
f==Intial_f;
end
equations
if Q>=0
    if f>=1 && T<(T_melt-diff)
        Q==mass*sp_heat*T.der
        f.der==0;
    elseif f<=0 && T>(T_melt+diff)
        Q==mass*sp_heat*T.der
        f.der==0;
    else
        if f<0
            Q==mass*sp_heat*T.der
            f.der==0
        else
            Q==-mass*la_heat*f.der
            T.der==0;
        end
    end
else
    if f>=1 && T<(T_melt-diff)
        Q==mass*sp_heat*T.der
        f.der==0;
    elseif f<=0 && T>(T_melt+diff)
        Q==mass*sp_heat*T.der
        f.der==0;
    else
        if f>1
            Q==mass*sp_heat*T.der
            f.der==0
        else
            Q==-mass*la_heat*f.der
            T.der==0;
        end
    end
end
end
connections
connect(M,N)
end
end