MATLAB: How to simulate a simple thermodynamical system with simulink

MATLABMATLAB and Simulink Student SuiteSimscapesimulinkthermalthermodynamic

Hello everybody!

I'd like to simulate the thermodynamical model of a Stirling engine using Simulink. The situation is the following:

I would like to know the temperature after a certain time (10 seconds for example) if I place a Heat source under my engine but I have to take care of the conductive heat transfer, the convective heat transfer and the radiative heat transfer. I've done this:

But I don't think that it's ok because Simulink sent me this error message:"The model may not give enough information to make it possible to solve for values of some of its variables. Specific advice is given below.

Cannot solve for one or more variables: 'Conductive_Heat_Transfer.A.T' (Temperature) 'Conductive_Heat_Transfer.B.T' (Temperature) 'Convective_Heat_Transfer.A.T' (Temperature) 'Convective_Heat_Transfer.B.T' (Temperature) 'Ideal_Heat_Flow_Source.B.T' (Temperature) 'Ideal_Heat_Flow_Source.T' (T) 'Ideal_Temperature_Sensor.A.T' (Temperature) 'Ideal_Temperature_Sensor.T' (T) 'Ideal_Temperature_Sensor.Ts' (Ts) 'Radiative_Heat_Transfer.A.T' (Temperature) 'Radiative_Heat_Transfer.B.T' (Temperature)

Equations of one or more components may be dependent or inconsistent. This can cause problems in transient initialization. Here is the set of components involved: 'untitled/Ideal Heat Flow Source' Equation location is: 'C:\Program Files\MATLAB\R2016b\toolbox\physmod\simscape\library\m\+foundation\+thermal\+sources\heat_flow.ssc' (line 34)

'untitled/Ideal Temperature Sensor' Equation location is: 'C:\Program Files\MATLAB\R2016b\toolbox\physmod\simscape\library\m\+foundation\+thermal\+sensors\temperature.ssc' (line 34) Component:Simulink | Category:Modelwarning An error occurred while running the simulation and the simulation was terminated Caused by: Trouble solving algebraic equations in differential-algebraic system. Singular iteration matrix encountered with step size 0.2 at time 0.0. Try tightening the error tolerances. If the problem persists, check the model structure and values of parameters."

Can someone help me ? thanks!

Best Answer

Hi,

I think the issue is occurring because of the way that the system is being modeled. If you add a thermal mass block after the heat transfer blocks, the model will be able to converge to a solution and thus simulate without errors.

Checkout the link below for some examples of modeling Thermal Systems with Simscape: https://www.mathworks.com/help/releases/R2017a/physmod/simscape/examples.html#d0e927

Good luck and have fun!

Related Solutions

MATLAB: How can get the value for each time step for a variable inside of a function

Hi,

you try this code (all in one file as nested function). The values you want to have in every iteration are now added to the solution and also saved in the results table. This should be what you want.

Result = main
function Result = main
%Initial conditions
Wads=0.045;
Tads=80;
Wdes=0.035;
Tdes=27;
Teva=7;
Tcond=27;
Mweva=0.65;         
Initial=[Wads Tads Wdes Tdes Mweva Teva Tcond 0 0 0 0];           
tspan=[1:1:800];
opts = odeset('RelTol',1e-6,'AbsTol',1e-6);
[t,Soln]=ode45(@(t,Soln)ODEdw(t,Soln),tspan,Initial,opts);
colNames = {'Time','Wads','Tads','Wdes','Tdes','Mweva', 'Teva','Tcond',...
    'Tcwaout', 'Thwdout', 'Tchwout', 'Tcwcout'};
Result = array2table([tspan',Soln],'VariableNames',colNames);
function [dYdt] =ODEdw(~,Soln)
Ea=4.2*10^4;         %Activation energy of surface diffusion, (J/mol)
R=8.314;             %Ideal gas constant, (j/mol.K)
Dso=2.54*10^-4;       %Pre-exponential constant, (m^2 S^-1)
Qs=2.8*10^6;         %Heat of adsorption, (J/kg)
%%%%%%%%%%%Operating conditions%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Thwdin=85;            %Inlet temperature of hot water in desorber, (oC)
Tcwain=27;            %Inlet temperature of cooling water in adsorber, (oC)
Tcwcin=27;            %Inlet temperature of cooling water in condenser, (oC)
Tchwin=14;            %Inlet temperature of chilling water in evaporator, (oC)
mcwa=8/60;            %Mass flow rate of cooling water in adsorber, (kg/s)
mcwc=20/60;           %Mass flow rate of cooling water in condenser, (kg/s)
mchw=1.6/60;          %Mass flow rate of chilling water in evaporator, (kg/s)
mhwd=7/60;            %Mass flow rate of hot water in desorber, (kg/s)
%%%%%%%%%Physical and thermal properties of adsorbent%%%%%%%%%%
Cps=924;              %Specific heat capacity of silica gel, (J/kg.K)
Ms=9;                 %Mass of silica, (kg)
Rp=5.6*10^-4;         %Radious of silica gel particle, (m)
ps=700;              %density of silica gel, (kg/m^3)
Por=0.3;              %Porosity of silica gel,
%%%%%%Physical and thermal properties of adsorbate%%%%%%%%%%%%%
Cpw=4200;             %Specific heat capcity of water, (J/kg.k)
Cpv=1900;             %Specific heat capacity of water vapor (J/kg.K)
M=18.02;              %Molar mass of water, (g/mol)
Lv=2500*10^3;         %Lantent heat of vaporization (J/kg) 
%%%%%%%%%Physical and thermal Properties of Condenser%%%%%%%%%%%%%
Acond=5.76;            %Surface area of condenser, (m^2)
Mhexcond=19.87;        %Mass of heat exchanger (Condenser), (kg)
Ucond=486.87;          %Heat transfer coefficient of condenser, (W/m^2.K)
UAcond=Acond*Ucond;
Mwcond=1.0;         %Mass of water in Condenser (kg)
%%%%%%%%%%Physical and thermal properties of adsorber/desorber   %%%%%%%
Mhex=65;             %Mass of heat exchanger (Adsorber/desorber), (kg)
Cphex=386;           %Heat capacity of heat excahgner, (J/kg.K)
Aad=2.11;            %Surface area of adsorber/desorber, (m^2)
Uads=189.60;         %Heat transfer coefficient of adsorber, (W/m^2.K)
Udes=203.99;         %Heat transfer coefficient of desorber, (W/m^2.K)
UAads=Aad*Uads;
UAdes=Aad*Udes;
Vt=0.5;              %volume of adsorber/desorber tank, (m^3)
%%%%%%%%%Physical and thermal Properties of Evaporator%%%%%%%%%%%%%%
Aeva=0.34;           %Surface area of evaporator, (m^2)
Mhexeva=16.1;        %Mass of heat exchanger (Evaporator), (kg)
Ueva=302.58;         %Heat transfer coefficient of evaporator, (W/m^2.K)
UAeva=Aeva*Ueva;
Wads=Soln(1);    
Tads=Soln(2);  
Wdes=Soln(3);    
Tdes=Soln(4);
Mweva=Soln(5);
Teva=Soln(6);
Tcond=Soln(7);
Kads=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tads+273.16)));   %Overall mass tarsnfer coeffcient of adsorber

Pwa=133.32*exp(18.3-(3820/((Teva+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)

Psga=133.32*exp(18.3-(3820/((Tads+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)

Weqa=0.346*(Pwa/Psga)^0.625;
Aads=(Kads*(Weqa-Wads));
TLads=(Ms*(Cpw+Cpw*Wads)+Mhex*Cphex); %Left hand side of adsorber equation
Tcwaout=Tads+(Tcwain-Tads)*exp((-UAads)/(mcwa*Cpw));
   
Kdes=((15*Dso)/(Rp).^2)*exp(-Ea/(R*(Tdes+273.16)));   %Overall mass tarsnfer coeffcient of adsorber
Pwd=133.32*exp(18.3-(3820/((Tcond+273.16)-46.1))); % Saturated vapour pressure of waterm (vapour) at adsorbent temperature, (Pa)
Psgd=133.32*exp(18.3-(3820/((Tdes+273.16)-46.1))); % Saturated vapour pressure of water at adsorbent temperature, (Pa)
Weqd=0.346*(Pwd/Psgd)^0.625;
Ades=(Kdes*(Weqd-Wdes));
TLdes=(Ms*(Cps+Cpw*Wdes)+Mhex*Cphex); %Left hand side of desorber equation
Thwdout=Tdes+(Thwdin-Tdes)*exp((-UAdes)/(mhwd*Cpw));
  
TLeva=(Mweva*Cpw+Mhexeva*Cphex); %Left hand side of evaporator equation

Tchwout=Teva+(Tchwin-Teva)*exp((-UAeva)/(mchw*Cpw));
  
TLcond=(Mwcond*Cpw+Mhexcond*Cphex); %Left hand side of evaporator equation
Tcwcout=Tcond+(Tcwcin-Tcond)*exp((-UAcond)/(mcwc*Cpw));
dYdt=[Kads*(Weqa-Wads);((Qs*Ms*Aads)+(Ms*Cpv*Aads*(Teva-Tads))+mcwa*Cpw*(Tcwain-Tcwaout))/TLads;
    Kdes*(Weqd-Wdes); ((Qs*Ms*Ades)+(mcwa*Cpw*(Thwdin-Thwdout)))/TLdes; 
    -Ms*((Kads*(Weqa-Wads))+(Kdes*(Weqd-Wdes)));((-Lv*Ms*Aads)+Ades*Ms*Cpw*(Tcond-Teva)+mchw*Cpw*(Tchwin-Tchwout))/TLeva;
    (Ades*Ms*Cpv*(Tdes-Tcond)+(Lv*Ms*Ades)+mcwc*Cpw*(Tcwcin-Tcwcout))/TLcond; Tcwaout; Thwdout; Tchwout; Tcwcout];
end
end

Best regards

Stephan

MATLAB: Simscape block parameter, change dynamically during simulation

To do this, create a custom block by modifying the shipping Conductive Heat Transfer block to make the "Thermal conductivity" parameter an input.

Open the source code for the Conductive Heat Transfer block
Save it to a new .ssc file.
Add an inputs section (you can look at the source code for the temperature source block to see what it should look like -- essentially a line with "inputs" and another with "end").
Move the line "th_cond = { 401, 'W/(m*K)' }; % Thermal conductivity" from the "parameters" section of the file to the newly created "inputs" section you just added.
Get a Simscape Components block from the Foundation Utilities Library and point it at that new .ssc file you just created.

Congratulations! You have just created your first custom Simscape component.

--Steve

Best Answer

Related Solutions

MATLAB: How can get the value for each time step for a variable inside of a function

MATLAB: Simscape block parameter, change dynamically during simulation

Related Question