Phase Change Material in Simscape

Dear Sir or Madam,

I would like to model a Phase Change Material in Simscape by editing a normal thermal mass block. I already found this question https://de.mathworks.com/matlabcentral/answers/300839-simscape-model-of-latent-heat-storage-in-a-pcm but it's not working when I put the mentioned code to the thermal mass block. Do you have any idea why it's not working? I'll attach the code which is the normal thermal mass block modified with the code I found on that link. But if I want to introduce a Simcape component block it always says "Line: 37 Column: 25 Unexpected MATLAB operator."

Many greetings,

Ralf

if true
  component mass
% Thermal Mass
% This block models internal energy storage in a thermal network. The rate
% of temperature increase is proportional to the heat flow rate into the
% material and inversely proportioanl to the mass and specific heat of the
% material.

% Copyright 2005-2016 The MathWorks, Inc.

nodes M = foundation.thermal.thermal; % :top end

parameters m = {1, 'kg' }; % Mass Cps = {200, 'J/(kg*K)'}; % Specific heat solid Cpl = {1000, 'J/(kg*K)'}; % Specific heat liquid Tmelt = {310, 'K'}; % Melting Temperature L = {100, 'J/(kg)'}; % Specific Latent Heat end

variables % Differential variables T = {value = {300, 'K'}, priority = priority.high}; % Temperature

    Q = {0, 'W'}; % Heat flow rate
end

branches Q : M.Q -> *; end

equations T == M.T;

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 end end

[Vishnubharathi thirupathi]

Hello Ralf Reisenauer ,

Did you get it running, the PCM block that you created?

[raju karri]

Can you please keep the corrected code

[David John]

It would be easier if you could attach your modified ssc file directly, but I can see an issue:

• "integral" is not a valid Simscape function. Please see https://www.mathworks.com/help/physmod/simscape/lang/integ.html instead. That integral will be over time, however, and not with respect to temperature. You cannot directly evaluate an integral with respect to temperature like you're trying to do here.
• You've also defined Q as a power (or heat flow), but your integral equation defines Q as an energy, so even if the syntax issues were solved, the equations still won't quite make sense.

What you can do, however, is something like the following in the equations section:

let
  specific_heat = if T <= Tmelt, Cps else Cpl end;
in
  Q == m * specific_heat * T.der;
end

This is fairly similar to the Thermal Mass code that you started from...just with an extra let-statement in there to modify the specific heat as appropriate.

Answer 1

[Vaishak]

The following code maybe useful. Separate Cp for liquid and solid can be added.

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