|
Dear all,
I wrote a small program in matlab using an m-file.
Now I would like to couple the m-file with the Trnsys Studio and call
the m-file every timestep. Therefore I used type 155 and I
followed the recommendations from the example I found in the manual. However I
keep getting an error at "mFileErrorCode"150, so something seems to be wrong
trying to do the necessary calculations. Does someone have an idea what I am
doing wrong? I am sending the m-file and
project-file in attachement.
Thanks a lot in advance,
Best regards,
marijke
ir-arch Marijke Steeman
Universiteit Gent Vakgroep Architectuur en Stedebouw Jozef Plateaustraat 22 9000 Gent tel 09 264 37 52 fax 09 264 41 85 marijke.steeman@ugent.be |
% Data passed from / to TRNSYS
% ----------------------------
%
% 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
% trnInputs
% ---------
%
% trnInputs(1) : ta_inlet
% trnInputs(2) : tf_inlet
% trnInputs(3) : tw_inlet
% trnInputs(4) : RHa_inlet
% trnInputs(5) : RHf_inlet
%
% trnOutputs
% ----------
%
% trnOutputs(1) : ta_outlet
% trnOutputs(2) : tf_outlet
% trnOutputs(3) : tw_outlet
% trnOutputs(4) : RHa_outlet
% trnOutputs(5) : RHf_outlet
mFileErrorCode = 100 % Beginning of the m-file
%--------------------------------------------------------------------------
%----PARAMETERS------------------------------------------------------------
%--------------------------------------------------------------------------
%dimensies warmtewisselaar
L=2; %lengte warmtewisselaar in m
W=1; %breedte warmtewisselaar in m
H=1; %hoogte warmtwisselaar in m
b_a=0.0045; %breedte tussen de platen in m
b_f=0.0045; %breedte tussen de platen in m
debiet=7100; %totaal luchtdebiet in m³/h
ro=1.2; %densiteit van lucht in kg/m³
%definieer constanten
zeta=6.1*10^(-6); %dampcapaciteit in kg/kg/Pa
N=100; %aantal knopen waarin temperaturen worden berekend
deltax=W*L/N; %lengte interval in m*breedte in m
Ca=1008; %specifieke warmte extractielucht J/kgK
Cf=1008; %specifieke warmte pulsielucht J/kgK
hev=2500000; %verdampingswarmte in J/kg
k_air=0.025; %warmtegeleiding van lucht in W/mK
Pr=0.707; %Prandl getal
nu=1.5*10^(-5); %kinematische viscositeit in m²/s
d_p=0.0001; %dikte van de platen in m
k_p=0.040; %warmtegeleidingscoefficient van de wisselaar in W/mK (polypropyleen)
d_w=0.00005; %dikte waterfilm in m
k_w=0.606; %warmtegeleidingscoefficient water in W/mK < uit paper
%berekening overgangscoefficienten en debieten per laag
aantal=H/(b_a+b_f+d_p+d_w); %aantal laagjes resp. in de twee stromen
Q=debiet/aantal; %debiet per laag in m³/h
Ga=(Q/2)*ro/3600; %half debiet per kanaal in kg/s
Gf=(Q/2)*ro/3600; %orde 2.5*10-5 %half debiet per kanaal in kg/s
D_a=2*W*b_a/(W+b_a); %hydraulische diameter van 1 laagje: [4*sectie/natte omtrek] in m
v_a=Q/(3600*b_a*W); % orde 2m/s %luchtsnelheid per luchtlaag in m/s
Re_a=v_a*D_a/nu; % orde 600 %Re berekend ahv hydraulische diam en debiet per laagje
D_f=2*W*b_f/(W+b_f);
v_f=Q/(3600*b_f*W);
Re_f=v_f*D_f/nu;
%Nusselt ook berekenen op basis van hydraulische diameter
if Re_a<2300;
Nu_a=((3.66^3)+(1.61^3)*(Re_a*Pr*D_a/L))^(1/3); %laminair
else zeta_a=(1.82*log(Re_a)-1.64)^(-2);
Nu_a=((zeta_a/8)*(Re_a-1000)*Pr)/(1+12.7*((zeta_a/8)^(0.5))*((Pr^(2/3))-1)); %turbulent
end
if Re_f <2300;
Nu_f=((3.66^3)+(1.61^3)*(Re_f*Pr*D_f/L))^(1/3);
else zeta_f=(1.82*log(Re_f)-1.64)^(-2);
Nu_f=((zeta_f/8)*(Re_f-1000)*Pr)/(1+12.7*((zeta_f/8)^(0.5))*((Pr^(2/3))-1));
end
%boek Heat Transfer ; A. BEJAN. p. 302 eq. (6.43)
%Nu (D) = h*D/k met D daarin de hydraulische diameter
ha=Nu_a*k_air/D_a;
hf=Nu_f*k_air/D_f;
%ha=Nu_a*k_air/(b_a/2); %warmteovergangscoefficient aan zijde extractielucht
%hf=Nu_f*k_air/(b_f/2); %warmteovergangscoefficient aan zijde pulsielucht
alfa_a=ha; %W/m²K
alfa_f=1/((d_p/k_p)+(d_w/k_w)+(1/hf)); %+weerstand warmtewisselaar en waterfilm W/m²K
%berekening dampovergangscoefficient via Lewis - relatie
Pair=101325; %luchtdruk Pa
beta=alfa_a/(Ca*Pair); %in kg/s/Pa/m² %orde 2*10-7
%iteraties
a=0.5; %relaxatie-factor
mFileErrorCode = 110 % After setting parameters
%--------------------------------------------------------------------------
% --- Process Inputs-------------------------------------------------------
% -------------------------------------------------------------------------
nI = trnInfo(3);
nO = trnInfo(6);
ta_inlet = trnInputs(1);%inlaat temperatuur extractielucht in °C
tf_inlet = trnInputs(2);%inlaat temperatuur pulsielucht in °C
tw_inlet = trnInputs(3); %inlaat temperatuur waterfilm in °C
RHa_inlet = trnInputs(4);%inlaat relatieve vochtigheid extractielucht in %
RHf_inlet = trnInputs(5);%inlaat relatieve vochtigheid pulsielucht in %
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)
iCall = 0;
% This is the first time step
iStep = 1;
mFileErrorCode = 130 % After initialization
end
%--------------------------------------------------------------------------
% --- All iterative calls -------------------------------------------------
% -------------------------------------------------------------------------
% --- If this is a first call in the time step, increment counter ---
if ( trnInfo(7) == 0 )
iStep = iStep+1;
end
mFileErrorCode = 150; % After reading inputs
%--------------------------------------------------------------------------
% --- Calculations --------------------------------------------------------
%--------------------------------------------------------------------------
%EERSTE RUN MET BEGINWAARDEN
pa_inlet=(RHa_inlet/100)*exp(65.8094-(7066.27/(273.15+ta_inlet))-5.976*log(273.15+ta_inlet));
pf_inlet=(RHf_inlet/100)*exp(65.8094-(7066.27/(273.15+tf_inlet))-5.976*log(273.15+tf_inlet));
M=zeros(3*N,3*N); %matrix met coefficienten
T=zeros(3*N,1); %matrix met onbekende temperaturen
D=zeros(3*N,1); %matrix met constanten
%coefficienten matrix M
aa=zeros(N,1);
ba=zeros(N,1);
ca=zeros(N,1);
aw=zeros(N,1);
bw=zeros(N,1);
cw=zeros(N,1);
af=zeros(N,1);
bf=zeros(N,1);
cf=zeros(N,1);
%coefficienten matrix D
da=zeros(N,1);
dw=zeros(N,1);
df=zeros(N,1);
%onbekende temperaturen in matrix T
ta=zeros(N,2); %3kolommen; eerste kolom oude waarde, tweede kolom nieuwe
tw=zeros(N,2); %3kolommen; eerste kolom oude waarde, tweede kolom nieuwe
tf=zeros(N,2); %3kolommen; eerste kolom oude waarde, tweede kolom nieuwe
%onbekende dampdrukken
pa=zeros(N,1);
pw=zeros(N,1);
%startwaarden
da(1,1)=-Ga*Ca*ta_inlet; %eerste coefficient matrix D
df(N,1)=-Gf*Cf*tf_inlet; %laatste coefficient matrix D
tw(1,1)=tw_inlet;
ta(1,1)=ta_inlet;
tf(N,1)=tf_inlet;
pw(1,1)=exp(65.8094-(7066.27/(273.15+tw(1,1)))-5.976*log(273.15+tw(1,1)));
pa(1,1)=(Ga*zeta*pa_inlet+beta*deltax*pw(1,1))/(Ga*zeta+beta*deltax);
%VOOR ALLE CONTROLEVOLUMES !!!!
for i=1:N
%schatten temperatuur = inlaatcondities in elke knoop van de warmtewisselaar
ta(i,1)=ta(1,1);
tw(i,1)=tw(1,1);
tf(i,1)=tf(N,1);
%berekening dampdruk wateropp uitgaande van tw_inlet
pw(i,1)=exp(65.8094-(7066.27/(273.15+tw(i,1)))-5.976*log(273.15+tw(i,1)));
end
%berekening dampdruk extractielucht
for i=2:N
pa(i,1)=(Ga*zeta*pa(i-1,1)+beta*deltax*pw(i,1))/(Ga*zeta+beta*deltax);
end
for i=1:N
%%WARMTESYSTEEM
%berekening coefficienten matrix M (veronderstelling deltax constant)
aa(i,1)=Ga*Ca;
ba(i,1)=-(Ga*Ca+deltax*alfa_a);
ca(i,1)=deltax*alfa_a;
aw(i,1)=deltax*alfa_a;
bw(i,1)=-deltax*(alfa_a+alfa_f);
cw(i,1)=deltax*alfa_f;
af(i,1)=deltax*alfa_f;
bf(i,1)=-(Gf*Cf+deltax*alfa_f);
cf(i,1)=Gf*Cf;
dw(i,1)=-hev*beta*deltax*(pa(i,1)-pw(i,1));
%plaatsing elementen in matrix M
M(3*i-2,3*i-2)=ba(i,1);
M(3*i-2,3*i-1)=ca(i,1);
M(3*i-1,3*i-2)=aw(i,1);
M(3*i-1,3*i-1)=bw(i,1);
M(3*i-1,3*i)=cw(i,1);
M(3*i,3*i-1)=af(i,1);
M(3*i,3*i)=bf(i,1);
%plaatsing constanten in matrix D
D(3*i-1,1)=dw(i,1);
end
D(1,1)=da(1,1);
D(3*N,1)=df(N,1);
%aa(1) en cf(N) moet hij niet plaatsen!
for i=2:N
M(3*i-2,3*i-5)=aa(i,1);
end
for i=1:N-1
M(3*i,3*i+3)=cf(i,1);
end
%M*T=D
A=inv(M);
T=A*D;
for i=1:N
%uitlezen matrix T
ta(i,2)=T(3*i-2,1);
tw(i,2)=T(3*i-1,1);
tf(i,2)=T(3*i,1);
end
taout=ta(N,2); %ter controle: waarden na eerste run
twout=tw(N,2);
tfout=tf(1,2);
%--------------------------------------------------------------------------
%ITEREREN SYSTEEM TOT CONVERGENTIEVOORWAARDE VOLDAAN
%--------------------------------------------------------------------------
q=0;
%eerste kolom = oude waarde
while ((abs(ta(N,2)-ta(N,1)))> 0.01 | (abs(tf(1,2)-tf(1,1)))> 0.01 | (abs(tw(N,2)-tw(N,1)))> 0.01 | (q<=50))
%((abs(ta(N,2)-ta(N,1))> 0.1)| (q<=10) ) %oude - nieuwe waarde => convergentievoorwaarde !!
%schuif "oude nieuwe" waarde naar "oude" waarde: shift kolom 2 ->1
q=q+1
for k=1:N
ta(k,1)=ta(k,1)+a*(ta(k,2)-ta(k,1));
tw(k,1)=tw(k,1)+a*(tw(k,2)-tw(k,1));
tf(k,1)=tf(k,1)+a*(tf(k,2)-tf(k,1));
pw(k,1)=exp(65.8094-(7066.27/(273.15+tw(k,1)))-5.976*log(273.15+tw(k,1)));
end
pa(1,1)=(Ga*zeta*pa_inlet+beta*deltax*pw(1,1))/(Ga*zeta+beta*deltax);
for k=2:N
pa(k,1)=(Ga*zeta*pa(k-1,1)+beta*deltax*pw(k,1))/(Ga*zeta+beta*deltax);
end
D=zeros(3*N,1);
M=zeros(3*N,3*N);
T=zeros(3*N,1);
for k=1:N
%%WARMTESYSTEEM
%berekening coefficienten
aa(k,1)=Ga*Ca;
ba(k,1)=-(Ga*Ca+deltax*alfa_a);
ca(k,1)=deltax*alfa_a;
aw(k,1)=deltax*alfa_a;
bw(k,1)=-deltax*(alfa_a+alfa_f);
cw(k,1)=deltax*alfa_f;
af(k,1)=deltax*alfa_f;
bf(k,1)=-(Gf*Cf+deltax*alfa_f);
cf(k,1)=Gf*Cf;
dw(k,1)=-hev*beta*deltax*(pa(k,1)-pw(k,1));
%plaatsing elementen in matrix M
M(3*k-2,3*k-2)=ba(k,1);
M(3*k-2,3*k-1)=ca(k,1);
M(3*k-1,3*k-2)=aw(k,1);
M(3*k-1,3*k-1)=bw(k,1);
M(3*k-1,3*k)=cw(k,1);
M(3*k,3*k-1)=af(k,1);
M(3*k,3*k)=bf(k,1);
%plaatsing constanten in matrix D
D(3*k-1,1)=dw(k,1);
end
D(1,1)=da(1,1);
D(3*N,1)=df(N,1);
%aa(1) en cf(N) moet hij niet plaatsen!
for k=2:N
M(3*k-2,3*k-5)=aa(k,1);
end
for k=1:N-1
M(3*k,3*k+3)=cf(k,1);
end
%M*T=D
B=inv(M);
T=B*D;
%uitlezen matrix T :: nieuwe waarden
for k=1:N
ta(k,2)=T(3*k-2,1);
tw(k,2)=T(3*k-1,1);
tf(k,2)=T(3*k,1);
end
ta(N,2);
tw(N,2);
tf(1,2);
end
%outputs
ta_outlet=ta(N,2);
tf_outlet=tf(1,2);
tw_oulet=tw(N,2);
pa_outlet=pa(N,1);
%berekening RH outputs
pa_sat=exp(65.8094-(7066.27/(273.15+ta_outlet))-5.976*log(273.15+ta_outlet));
RHa_outlet=100*pa_outlet/pa_sat
pf_sat=exp(65.8094-(7066.27/(273.15+tf_outlet))-5.976*log(273.15+tf_outlet));
RHf_outlet=100*pf_inlet/pf_sat
%--------------------------------------------------------------------------
% --- Set outputs----------------------------------------------------------
%--------------------------------------------------------------------------
trnOutputs(1) = ta_outlet;
trnOutputs(2) = tf_outlet;
trnOutputs(3) = tw_outlet;
trnOutputs(4) = RHaf_outlet;
trnOutputs(5) = RHf_outlet;
mFileErrorCode = 0; % Tell TRNSYS that we reached the end of the m-file without errors
return
Attachment:
koppeling trnsys - matlab.TPF
Description: Binary data