Pond model, calling Matlab from TRNSYS errorcode 150, info(7), info(13)

I am trying to simulate heat transfers in a water pond model and couple with Trnsys Type 155. I get an error in TRNSYS mFileErrorCode 150, at info(7) = 0 , info(13) = 0. What can be the problem?
mFileErrorCode = 100; % Beginning of the m-file
% Pond_test.m
% -----------Pond parameters------------------------------------------------------------------------------------------------
%Upper layer thickness (m) (ul)
Xul=0.3;
%Middle layer elements
n=800;
%Middle layer thickness (m) (ml)
Xml=8;
%Delta x midle layer
dxw=0.01;
%Bottom layer thickness (m) (bl)
Xbl=0.7;
%Concrete layer thickness (m) (cl)
Xcl=1;
%Concrete layer elements
ncl=100;
%Delta x concrete layer
dxcl=Xcl/ncl;
%Concrete parameters
%%Density (RO) (kg/m3)
ROco=2300;
%%Specific heat capacity (Cp) (J/kg*K)
Cpco=1000;
%%Thermal conductivity (k) (W/(m*K))
kco=0.8;
% Thermal diffustivity concrete (alfa) [m2/s]
acl= 0.000000347;
%Water parameters
%%Density (RO) (kg/m3)
ROwa=1000;
%%Specific heat capacity (Cp) (J/kg*K)
Cpwa=4190;
%%Thermal conductivity (k) (W/(m*K))
kwa=0.6;
% Thermal diffustivity water (alfa) [m2/s]
aw= 0.000000143;
%Area pond [m2]
Apond=50000;
% trnTime (1x1) : simulation time
% trnInfo (15x1) : TRNSYS info array
% trnInputs (nIx1) : TRNSYS inputs
% trnStartTime (1x1) : TRNSYS Simulation Start time
% trnStopTime (1x1) : TRNSYS Simulation Stop time
% trnTimeStep (1x1) : TRNSYS Simulation time step
% mFileErrorCode (1x1) : Error code for this m-file. It is set to 1 by TRNSYS and the m-file should set it to 0 at the
% end to indicate that the call was successful. Any non-zero value will stop the simulation
% trnOutputs (nOx1) : TRNSYS outputs
mFileErrorCode = 110; % After setting parameters
% --- Process Inputs from TRNSYS ---------------------------------------------------------------------------------------------------
% ----------------------------------------------------------------------------------------------------------------------
Theap = trnInputs(1); % °C
mdotheap = trnInputs(2); % kg/h
Tamb = trnInputs(3); % °C
Hhor = trnInputs(4); % kJ/hr-m2
Wind = trnInputs(5); % m/s
RH = trnInputs(6); % %
Tsky = trnInputs(7); % °C
AI = trnInputs(8); % °
mFileErrorCode = 120; % After processing inputs
% --- First call of the simulation: initial time step (no iterations) --------------------------------------------------
% ----------------------------------------------------------------------------------------------------------------------
% (note that Matlab is initialized before this at the info(7) = -1 call, but the m-file is not called)
if ( (trnInfo(7) == 0) && (trnTime-trnStartTime < 1e-6) )
% This is the first call (Counter will be incremented later for this very first call)
nCall = 0;
% This is the first time step
nStep = 1;
% Initialize history of the variables for plotting at the end of the simulation
nTimeSteps = (trnStopTime-trnStartTime)/trnTimeStep + 1;
history.Theap = zeros(nTimeSteps,1);
history.mdotheap = zeros(nTimeSteps,1);
history.Tamb = zeros(nTimeSteps,1);
history.Hhor = zeros(nTimeSteps,1);
history.Wind = zeros(nTimeSteps,1);
history.RH = zeros(nTimeSteps,1);
history.Tsky = zeros(nTimeSteps,1);
history.AI = zeros(nTimeSteps,1);
history.To = zeros(nTimeSteps,1)+15;
history.Tu = zeros(nTimeSteps,1)+15;
history.TN2 = zeros(nTimeSteps,1)+15;
history.Tcl1 = zeros(nTimeSteps,1)+15;
history.Tcl2 = zeros(nTimeSteps,1)+15;
history.Tcli = zeros(nTimeSteps,1)+15;
history.Tb = zeros(nTimeSteps,1)+15;
history.ql = zeros(nTimeSteps,1);
history.Ti = zeros(nTimeSteps,1)+15;
% No return, calculate the pls pond performance during
% this call1
mFileErrorCode = 130; % After initialization
end
% --- Very last call of the simulation (after the user clicks "OK"): Do nothing ----------------------------------------
% ----------------------------------------------------------------------------------------------------------------------
if ( trnInfo(8) == -1 )
mFileErrorCode = 1000;
mFileErrorCode = 0; % Tell TRNSYS that we reached the end of the m-file without errors
return
end
% --- Post convergence calls: store values -----------------------------------------------------------------------------
% ----------------------------------------------------------------------------------------------------------------------
if (trnInfo(13) == 1)
mFileErrorCode = 140; % Beginning of a post-convergence call
history.Theap(nStep) = Theap;
history.mdotheap(nStep) = mdotheap;
history.Tamb(nStep) = Tamb;
history.Hhor(nStep) = Hhor;
history.Wind(nStep) = Wind;
history.RH(nStep) = RH;
history.Tsky(nStep) = Tsky;
history.AI(nStep) = AI;
history.To(nStep+1) = To;
history.Tu(nStep+1) = Tu;
history.TN2(nStep+1) = TN2;
history.Ti(nStep) = Ti;
history.Tb(nStep) = Tb;
history.Tcl1(nStep+1) = Tcl1;
history.Tcl2(nStep+1) = Tcl2;
history.Tcli(nStep+1) = Tcli;
history.qloss(nStep) = qloss;
mFileErrorCode = 0; % Tell TRNSYS that we reached the end of the m-file without errors
return % Do not update outputs at this call
end
% --- All iterative calls ----------------------------------------------------------------------------------------------
% ----------------------------------------------------------------------------------------------------------------------
% --- If this is a first call in the time step, increment counter ---
if (trnInfo(7) == 0)
nStep=nStep+1;
end
% --- Get TRNSYS Inputs ---
nI = trnInfo(3); % For bookkeeping
nO = trnInfo(6); % For bookkeeping
Theap = trnInputs(1);
mdotheap = trnInputs(2);
Tamb = trnInputs(3);
Hhor = trnInputs(4);
Wind = trnInputs(5);
RH = trnInputs(6);
Tsky = trnInputs(7);
AI = trnInputs(8);
mFileErrorCode = 150; % After reading inputs
% --- Calculate PLS pond performance ---
%elements thickness or space steps
%dt=300
dt=trnTimeStep/3600;
%Fourier concrete
wc =acl*dt/(dxcl*dxcl);
%Fourier water
ww =aw*dt/(dxw*dxw);
%Refraction angle (RA) and Snells Law
RA = asind((sind(AI))/(1.333));
%solar radiation in pond
for i=1:n
x(i) =((i-1)*dxw)+Xul;
H(i) =(Hhor/3.6)*(0.36-(0.08*(log(x(i)/cosd(RA)))));
end
Hb =(Hhor/3.6)*(0.36-(0.08*(log(Xul+Xml+Xbl)/cosd(RA))));
%Initial condition
for i=1:n
T(i,1) =15;
end
for i=1:ng
Tg(i,1)=15;
end
Tb(1) =15;
Tg(1) =15;
j =nStep;
%------------------------------------------------------------
%------Heat loss pond surface--------------------------------
tsu(j) = history.Tu(nstep);
Patm = 101325;
La = 2260;
Cs = 1.001;
hc = 5.7+(3.8*Wind);
P1(j) = exp(18.403-(3885/(tsu(j)+230)));
Pa = RH*exp(18.403-(3885/((Tamb)+230)));
STbc = 0.0000000567;
AbW = 0.85;
%Evaporation
qe(j) = (La*hc*(P1(j)*Pa))/(1.6*Cs*Patm);
%Convection
qc(j) = hc*(tsu(j)-Tamb);
%Radiation
qr(j) = AbW*STbc*(((tsu(j)+273)^4)-((Tsky+273)^4));
%Tot heat loss
qloss(j) = qe(j)+qc(j)+qr(j);
%---------------------------------------------------
%-------------Boundary conditions-------------------
%------------------Upper layer----------------------
T(1,j+1)=(history.Tu(nStep))+((dt/(ROwa*Cpwa*Xul))*(((Hhor/3.6)-H(1))+((kwa/((dxw/2)))*((history*TN2(nStep))-(history.Tu(nStep))))-(qloss(j))));
%-----------Internal nodes middle layer---------------------
for i=2:n-1
T(i,j+1) = (T(i,j))+(ww(T(i+1,j)-2*T(i,j)+T(i-1,j)))+((ww*(dxw/kwa))*(H(i-1)-H(i)));
end
%From heap mdot[kg/s] Cpwa [kj/kg*K] T [K]
qinheap = ((mdotheap/(3600))*(Cpwa/1000)*(history.Tb(nStep)-Theap))/50000; %pond area 50000m2
%------------Lower zone-------------------------------------
T(n,j+1)=((history.To(nStep))+((dt/(ROwa*Cpwa*Xbl)))*((H(n-1))-((kwa/(dxw/2))*((history.To(nStep))-(T(n-1,j))))-(qinheap)-((kco/dxcl)*((history.Tcl(nStep))-(history.Tcl(nStep))))));
%------------Concrete heat loss-----------------------------
%-------------First point layer-----------------------------
Tcl(1,j+1)=(history.Tcl1(nStep))+((wc)*(history.Tb(nStep)-(2*(history.Tcl1(nStep)))+(history.Tcl2(nStep))));
%-------Internal node concrete layer------------------------
for i=2:ncl-1
Tcl(i,j+1)=(Tcl(i,j))+((wc)*(Tcl(i+1,j)-(2*(Tcl(i,j)))+(Tcl(i-1,j))));
end
%----Temp ground=Temp concrete lower layer------------------
Tcl(ncl,j+1)=25;
%-----------------Bottom pond-------------------------------
hbottom = 200;
Tb(j+1)=(((Hb)+((kco/dxcl)*(Tg(1,j+1)))+(hbottom*(T(n,j+1))))/((kco/dxcl)+hbottom));
%--------------------Node temp and qloss--------------------
To = T(n,j+1);
Tu = T(1,j+1);
Ti = T(i,j+1);
TN2 = T(2,j+1);
ql = qloss(j);
Tcl1 = Tcl(1,j+1);
Tcl2 = Tc2(2,j+1);
Tcli = Tcli(i,j+1);
Tb = Tb(j+1);
%-------------------------Outputs---------------------------
trnOutputs(1) = To;
trnOutputs(2) = Tu;
trnOutputs(3) = Hhor;
trnOutputs(4) = Ti;
trnOutputs(5) = TN2;
trnOutputs(6) = ql;
trnOutputs(7) = Tb;
trnOutputs(8) = Tcl1;
trnOutputs(9) = mdot;
trnOutputs(10) = qinheap;
mFileErrorCode = 0; %Tell trnsys no errors occurred

Answers (1)

Categories

Asked:

on 21 May 2020

Answered:

on 12 Mar 2021

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!