Steady state current or temperature of bare cable

current.m

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