Skip to content

Instantly share code, notes, and snippets.

@c1b3rh4ck c1b3rh4ck/current.m
Created Nov 10, 2013

Embed
What would you like to do?
Steady state current or temperature of bare cable
function itcabledesnudo( )
% Program for computing steady state current or temperature of bare cable
% IEEE Standard 738 - 2006
% Carlos J. Zapata
%Modified by : Hector Jimenez S.
% Check the last updates.
% Last update
% Oct 01 - 2012: A bug in qsun function was removed
% May 06 - 2012: The progran was released
clear
% Input data --------------------------------------------------------------
caso=1; % Case to analyse 1-->I given T, 2-->T given I
Y=100; % T in [centigrade] for case 1, I in [A] for case 2
D=28.14; % Diameter of cable [mm]
T1=25; % Temperature for ac resistence 1 [centigrades]
RT1=7.283E-05; % AC cable resistence at temperature 1 [Ohm/m]
T2=75; % Temperature for ac resistence 2 [centigrades]
RT2=8.688E-05; % AC cable resistence at temperature 1 [Ohm/m]
e=0.50; % Cable emissivity [adimentional] 0.23 to 0.91
a=0.50; % Cable solar absortivity [adimentional) 0.23 to 0.91
fi=90; % Angle between wind and cable axis [degrees]
He=100; % Average altitude of conductor above sea level [m]
Zl=90; % Azimuth of line [degrees]
Lat=30; % Latitude where the cable is located [degrees]
N=161; % Day of year, 1 to 365
hour=14; % Hour of the day, 1 to 24
Vw=0.61; % Speed of wind [m/s]
Ta=40; % Ambient air temperature [centigrade]
atmosphere=1; % 1-->clear, 2-->industrial
% ------------------------ MAIN PROGRAM -----------------------------------
clc
disp('PROGRAM FOR COMPUTING BARE CABLE CURRENT-TEMPERATURE')
disp(' ')
disp('STANDARD IEEE 738 - 2006')
disp(' ')
switch caso
case 1
printindata()
[X]=it(Y);
disp(' ')
disp(['I=',num2str(X),' [Amperes] for a conductor temperature of ',num2str(Y),' [centigrades]'])
disp(' ')
case 2
printindata()
[X]=ti(Y);
disp(' ')
disp(['Tc=',num2str(X),' [centigrades] for a conductor current of ',num2str(Y),' [Amperes]'])
disp(' ')
otherwise
disp('ERROR: case of study no valid')
end
% ------------------------- EMBEBED FUNCTIONS -----------------------------
function [x]=it(Y) %-------------------------------------------------------
qc=convected(Y);
qr=radiated(Y);
[qs]=qfromsun( );
RTC=resistance(Y);
I=sqrt((qc+qr-qs)/RTC);
x=I;
salida0=[{'qc'} {num2str(qc)} {'[W/m]'}
{'qr'} {num2str(qr)} {'[W/m]'}
{'qs'} {num2str(qs)} {'[W/m]'}
{'RTC'} {num2str(RTC)} {'[Ohm/m]'}
];
disp(' ')
disp(salida0)
end %----------------------------------------------------------------------
function [x]=ti(y) %-------------------------------------------------------
end %----------------------------------------------------------------------
function [qc]=convected(Tc) %----------------------------------------------
Kangle=1.194-cos(fi*pi/180)+0.194*cos(2*fi*pi/180)+0.368*sin(2*fi*pi/180);
Tfilm=(Tc+Ta)/2;
kf=2.424E-02+7.477E-05*Tfilm-4.407E-09*Tfilm^2;
mhuf=1.458E-06*(Tfilm+273)^1.5/(Tfilm+383.4);
rhof=(1.293-1.525E-04*He+6.379E-09*He^2)/(1+0.00367*Tfilm);
qcn=0.0205*rhof^0.5*D^0.75*(Tc-Ta)^1.25;
qc1=(1.01+0.0372*(D*rhof*Vw/mhuf)^0.52)*kf*Kangle*(Tc-Ta);
qc2=0.0119*(D*rhof*Vw/mhuf)^0.60*kf*Kangle*(Tc-Ta);
if Vw==0 | fi==0
qc=qcn;
else
qc=max(qc1,qc2);
end
salida1=[{'Convection'} {' '} {' '}
{'Kangle'} {num2str(Kangle)} {' '}
{'Tfilm'} {num2str(Tfilm)} {'[centigrades]'}
{'kf'} {num2str(kf)} {'[W/(m-centigrade)]'}
{'mhuf'} {num2str(mhuf)} {'[Pa-s]'}
{'rhof'} {num2str(rhof)} {'[kg/m3]'}
{'qcn'} {num2str(qcn)} {'[W/m]'}
{'qc1'} {num2str(qc1)} {'[W/m]'}
{'qc2'} {num2str(qc2)} {'[W/m]'}
];
disp(' ')
disp(salida1)
end %----------------------------------------------------------------------
function [qr]=radiated(Tc) %-----------------------------------------------
qr=0.0178*D*e*(((Tc+273)/100)^4-((Ta+273)/100)^4);
end %----------------------------------------------------------------------
function [qs]=qfromsun( ) %------------------------------------------------
delta=23.4583*sin((284+N)/365*360*pi/180);
if hour==12
w=0;
else
w=+15*(hour-12);
end
Hc=asin(cos(Lat*pi/180)*cos(delta*pi/180)*cos(w*pi/180)+sin(Lat*pi/180)*sin(delta*pi/180));
chi=sin(w*pi/180)/(sin(Lat*pi/180)*cos(w*pi/180)-cos(Lat*pi/180)*tan(delta*pi/180));
if w<0 & w>=-180
if chi>=0
C=0;
else
C=180;
end
end
if w>=0 & w<=180
if chi>=0
C=180;
else
C=360;
end
end
Zc=C*pi/180+atan(chi);
tetha=acos(cos(Hc)*cos(Zc-Zl*pi/180));
[KA KB KC KD KE KF KG]=sunsky(atmosphere);
Qs=KA+KB*(Hc*180/pi)+KC*(Hc*180/pi)^2+KD*(Hc*180/pi)^3+KE*(Hc*180/pi)^4+KF*(Hc*180/pi)^5+KG*(Hc*180/pi)^6;
Ksolar=1+1.148E-04*He-1.108E-08*He^2;
Qse=Ksolar*Qs;
Ap=D/1000;
qs=a*Qse*sin(tetha)*Ap;
salida2=[{'Solar'} {' '} {' '}
{'delta'} {num2str(delta)} {'[degrees]'}
{'w'} {num2str(w)} {'[hour degree]'}
{'Hc'} {num2str(Hc*180/pi)} {'[degrees]'}
{'Chi'} {num2str(chi)} {' '}
{'C'} {num2str(C)} {'[degrees]'}
{'Zc'} {num2str(Zc*180/pi)} {'[degrees]'}
{'tetha'} {num2str(tetha*180/pi)} {'[degrees]'}
{'Ksolar'} {num2str(Ksolar)} {' '}
{'Qs'} {num2str(Qs)} {'[W/m]'}
{'Ksolar'} {num2str(Ksolar)} {' '}
{'Qse'} {num2str(Qse)} {'[W/m]'}
{'Ap'} {num2str(Ap)} {'[m2/m]'}
];
disp(salida2)
end %----------------------------------------------------------------------
function [rtc]=resistance(Tc) %--------------------------------------------
rtc=(RT2-RT1)/(T2-T1)*(Tc-T1)+RT1;
end %----------------------------------------------------------------------
function [A B C D E F G]=sunsky(atmosphere) %------------------------------
switch atmosphere
case 1 %Clear atmosphere
A=-42.2391;
B=+63.8044;
C=-1.9220;
D=+3.46921E-02;
E=-3.61118E-04;
F=+1.94318E-06;
G=-4.07608E-09;
case 2 %Industrial atmosphere
A=+53.1821;
B=+14.2110;
C=+6.6138E-01;
D=-3.1658E-02;
E=+5.464E-04;
F=-4.3446E-06;
G=+1.3236E-08;
end
end %----------------------------------------------------------------------
function printindata() %---------------------------------------------------
disp('Input Data:')
switch atmosphere
case 1 %Clear atmosphere
descrip='clear';
case 2 %Industrial atmosphere
desc='industrial';
end
salida0=[{'Cable'} {'-------'} {'--------'}
{'D'} {num2str(D)} {'[mm'}
{'Rac1'} {num2str(RT1)} {'[Ohm/m]'}
{'T1'} {num2str(T1)} {'[centigrade]'}
{'Rac2'} {num2str(RT2)} {'[Ohm/m]'}
{'T2'} {num2str(T2)} {'[centigrade]'}
{'e'} {num2str(e)} {' '}
{'alpha'} {num2str(a)} {' '}
{'Line'} {'-------'} {'--------'}
{'fi'} {num2str(fi)} {'[degrees]'}
{'He'} {num2str(He)} {'[m] '}
{'Zl'} {num2str(Zl)} {'[degrees]'}
{'Lat'} {num2str(Lat)} {'[degrees]'}
{'Weather'} {'-------'} {'--------'}
{'Day N'} {num2str(N)} {' '}
{'Hour'} {num2str(hour)} {' '}
{'Vw'} {num2str(Vw)} {'[m/s]'}
{'Ta'} {num2str(Ta)} {'[centigrade]'}
{'atmosphere'} {descrip} {' '}
];
disp(salida0)
end %----------------------------------------------------------------------
end %----------------------------------------------------------------------
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.