the value alpha is not updating with temperature loop(k=1:21).

6 views (last 30 days)
ABHISEK PATI
ABHISEK PATI on 6 Jan 2021
Commented: ABHISEK PATI on 7 Jan 2021
the value alpha and rcmin is not updating with temperature (k=1:21) due to unknown reason .the alpha value depends upon the rc value which i get as asn array after solving the ode45 equation in the j loop .the rc value changes with temperature but rcmin doesnt though it is a function of temperature.please help me if i am doing something wrong with format and sameprobelm happening for the resistance R.
and second problem is i am not able to get the value of Av after typing Av in command window it throws me an error unrecogined function or variable.
i have pasted my while code below please have a look and suggest me changes
%% going to replicate paper 6_Effect of Geometrical Parameters on the Thermal Performance of Ammonia-Based Trapezoidal-Shaped Axial
%% Grooved Heat Pipe
%% SUBROUTINE: Properties:
%it is always advisable to preallocate arrays or speed.
global hfg rowl kl meul sigmal Tv newl rowv meuv newv cpv R twall ks pv T;
R=462;
hfg=zeros(1,200);rowl=zeros(1,200);
kl=zeros(1,200);meul=zeros(1,200);
sigmal=zeros(1,200);newl=zeros(1,200);
rowv=zeros(1,200); meuv=zeros(1,200);
newv=zeros(1,200);pv=zeros(1,200);cpv=zeros(1,200);
Tv =[1:21]; T=[1:41];
Tv=Tv+273;
%declaring properties
for i=1:200
% %% LIQUID PROPERTIES
% %hfg-in j/kg (R=0.9975)
% hfg(i)=(-2.7621*i + 2519.7)*1000;
% %Liqu1 density %in kg/m3 R=0.985
% rowl(i) = -0.7406*i + 1025;
% %Thermal conductivity %in W/mk R=0.997
% kl(i)= -6*10^(-6)*(i^2) + 0.0017*i + 0.5728;
% % viscocity %in pascal sec R=0.9828
% meul(i)= (4*10^(-05)*i^2 - 0.0125*i + 1.1473)*(10^(-3));
% % surface tension %in N/m R=0.999
% sigmal(i)=(-0.0191*i + 7.7407)*(10^(-2));
% %kinematic viscocity %in pascal-sec R=0.999
% newl(i)=meul(i)/rowl(i);
% %% VAPOUR PROPERTIES
% %density %in kg/m3 R=0.992
% rowv(i) = 0.0161*exp(0.0330*i);
% % Vapour viscocity %in pascal-sec R=0.999
% meuv(i) =(0.0038*i + 0.8873)*(10^(-5));
% % kinematic viscocity %in pascal-sec R=0.999
% newv(i)=meuv(i)/rowv(i);
% % vapour pressure %in pascal R=0.987
% pv(i) =(0.3406*(i^2) - 2.1453*(i) + 2.7922)*(10^5);
% %specific heat %in KJ/kgK R=0.987
% cpv(i)=4*10^-05*i^2 - 0.0032*i + 1.9182;
%% LIQUID PROPERTIES
%hfg-in j/kg (R=0.9975)
hfg(i)= ((-523322100)+ (1223.915 - -523322100)/(1 + (i/250210.1)^1.767564))*1000;
%Liqu1 density %in kg/m3 R=0.985
rowl(i) = (-198353100) + (636.5187 - -198353100)/(1 + (i/5530476)^1.274674);
%Thermal conductivity %in W/mk R=0.997
kl(i)= -5109.54 + (0.2978967 - -5109.54)/(1 + (i/1083800)^1.192852);
% viscocity %in pascal sec R=0.9828
meul(i)= (-13431.15 + (0.2479605 - -13431.15)/(1 + (i/4038018)^1.086016))*(10^(-3));
% surface tension %in N/m R=0.999
sigmal(i)= (-3.130892 + (2.457097 - -3.130892)/(1 + (i/146.7721)^1.522479))*(10^(-2));
%kinematic viscocity %in pascal-sec R=0.999
newl(i)=meul(i)/rowl(i);
%% VAPOUR PROPERTIES
%density %in kg/m3 R=0.992
rowv(i) = 36484580 + (5.079686 - 36484580)/(1 + (i/28094.94)^2.398605);
% Vapour viscocity %in pascal-sec R=0.999
meuv(i) =(-2.307407 + (1.152 - -2.307407)/(1 + (i/100.6567)^128.5731))*(10^(-5));
% kinematic viscocity %in pascal-sec R=0.999
newv(i)=meuv(i)/rowv(i);
% vapour pressure %in pascal R=0.987
pv(i) =(32329700 + (4.697606 - 32329700)/(1 + (i/165106.9)^1.786858))*(10^5);
% %specific heat %in KJ/kgK R=0.987
% cpv(i)=4*10^-05*i^2 - 0.0032*i + 1.9182;
end
%%..........................................................................%%
%% Constants defined
global lLt Lt lLe Le lLc Lc La lLa Di Dv rv hg vew w2 w1 N Qin W pi wfin wo gammav phi A2 gamma
lLt=1000;
Lt=1;
lLe=300;
Le=0.3;
lLc=300;
Lc=0.3;
La=0.4;
lLa=400;
Di=10.5*10^(-3);
Dv=7.9*10^(-3);
rv=Dv/2;
hg=1.3;
vew=42*pi/180; %in radians
w2=0.79;%*(10^(-3));
w1=0.5;%*(10^(-3));
gammav=1.33;
N=32; %unknown
Qin=102; %unknown
W=wo/2; %unknown
pi=3.141;
gamma=(33*pi/180);
% w1= (pi*Dv)/(2*N);
wfin=w1/2;
% w2=(2*hg*cot(2*vew))+w1;
wo=w2+(2*wfin);
phi = asin(w1/(Dv*1000));
%.....................................................................
twall=1.1*(10^(-3)); % in mm taken from paper page 1 as Do-Di
ks=167; % taken aluminium
%..............................................................
global A1 fre0 A2;
A1= ((64*W)/((pi^5)*hg))*tanh((pi*hg)/(2*W));
fre0=8*(hg^2)/((W^2)*((1+(hg/W))^2)*((1/3)-A1));
A2= (1-(1.971*exp(-pi*hg/2/W)));
%% finding the value of rc min
%% solving for resistance
R1= twall./(ks*wfin);
R2= hg./(ks*wfin);
R4= twall./(ks*w2);
R7= hg./(ks*wfin);
R8= twall./(ks*wfin);
R10= twall./(ks*w2);
% %to solve
% for k=1:21
% Re(k)= ((R1+R2+R3(rc,k))*(R4+R5(rc)))/((Lc*N)*((R1+R2+R3(rc,k)+(2*R4+2*R5(rc)))));
% Rc(k)= ((R6(k)+R7+R8)*(R9(rc,k)+R10))/((Lc*N)*((R6(k)+R7+R8)+(2*R9(rc,k)+2*R10)));
% end
%% Defining some global parametrs
global mv Qz ml Qmax;m=0;
mv=zeros(1,lLt+1);Qz=zeros(1,lLt+1);ml=zeros(1,lLt+1);
for k=1:21
for j=1:50
for z=1:(lLt+1)
z_m=z/1000;
if (1<=z)&&(z<=lLe+1)
Qz(z)=(z_m/Le)*Qin;
elseif (lLe+1<z)&&(z<=(lLe+lLa+1))
Qz(z)=Qin;
else
Qz(z)=((Lt-z_m)/Lc)*Qin;
end
end
for i=1:(lLt+1)
ml(i)=Qz(i)/hfg(k);
mv(i)=-Qz(i)/hfg(k);
end
rc0=rv*1000;%initial assumption
zspan=[lLt+1 1];
[z,rc]=ode45(@(z,rc) (derpvz(z,rc,k)-derplz(z,rc,k))*(-(rc^2)/sigmal(k)),zspan,rc0);
rc=real(rc);
rcmin(k)=alpha(z,rc,k);
if(rc(41)>=rcmin(k))
Qin=Qin+0.1;
else
Qmax(k)=Qin;
end
end
%finding resistance
Re(k)= real(((R1+R2+R3(rc,k))*(R4+R5(rc)))/((Lc*N)*((R1+R2+R3(rc,k)+(2*R4+2*R5(rc))))));
Rc(k)= real(((R6(k)+R7+R8)*(R9(rc,k)+R10))/((Lc*N)*((R6(k)+R7+R8)+(2*R9(rc,k)+2*R10))));
Rt(k)=(2*Re(k)+Rc(k));
end
%% END OF MAIN %%
% Function definitions for the derivatives
%solving fo the ode of pv for the resistance
% for k=1
% pv0=pv(k);
% pvspan=[lLt 1];
% dpv=ode45(@derpvz(z,rc,k),pvspan,pv0);
% pspan=[1 1000];
% d=deval(dpv,pspan);
% end
%solving derpvz
function dpv= derpvz(z,rc,k)
global mv; global meuv;global rowv;z_i=round(z);
dpv= abs( real(((-2*meuv(k)*frev(z,rc,k)*mv(z_i))/((Dhv(z,rc))^2*Av(rc)*rowv(k)))));
return;
end
%% other function............................................................
%frev
function fv= frev(z,rc,k)
global mv meuv rowv hfg gammav lLt Dv pi av R Tv Qin;
av=Av(rc);Rev=zeros(1,lLt+1);Mav=zeros(1,lLt+1);
for i=1:lLt+1
Rev(i)= (4*mv(i))/(pi*Dv*meuv(k));
Mav(i)= mv(i)/((rowv(k)*av)*((gammav*R*Tv(k))^0.5));
end
z_i=round(z);
if (Rev(z_i)<2300)&&(Mav(z_i)<0.2)
fv=16;
elseif (Rev(z_i)<2300)&&(Mav(z_i)>0.2)
fv=16/(1-((gammav-1)*(Mav(z_i)^2)/2));
elseif (Rev(z_i)>2300)&&(Mav(z_i)<0.2)
fv=0.038*((Dv*Qin/(av*meuv(k)*hfg(k)))^0.75);
else
fv=0.608*((Dv*Qin/(av*meuv(k)*hfg(k)))^0.75)*(1/((1-(0.5*(gammav-1)*(Mav(z_i)^2)))));
end
end
%function Dhv
function dhv = Dhv(z,rc)
global pi rv phi N
psi_f=psi(rc);
dhv = (4*Av(rc))/((2*pi*rv)-(2*N*rv*phi)+(2*N*rc*psi_f));
end
%function Av
function Vapourarea= Av(rc)
global rv pi N w1 phi
psi_v=psi(rc);
Vapourarea=pi*(rv^2)+N*((rc^2*psi_v)-(rv^2)+0.5*(w1*rv*cos(phi))-0.5*(w1*rc*cos(psi_v)));
end
%psi(z)
function p =psi(rc)
global w1
p = asin(w1/2/rc);
end
%% writing all the functions of derplz
function pl=derplz(z,rc,k)
global meul ml rowl N ;z_i=round(z);
pl = abs(real(((-2*meul(k)*frel(z,rc,k)*ml(z_i))/((Dhl(z,rc))^2*Al(rc)*rowl(k))*N)));
end
function fl=frel(z,rc,k)
global fre0 N W pi newl newv A2
fl=fre0*(1+((4*N*(W^3))/(3*pi*(Dhv(z,rc)^3)))*frev(z,rc,k)*(newv(k)/newl(k))*A2);
end
function dl=Dhl(z,rc)
global w1 w2 hg N
dl=(4*Al(rc)/N)/(w2+(2*sqrt(((0.5*(w2-w1))^2+(hg)^2))));
end
function liquidarea=Al(rc)
shi=psi(rc);
global w1 w2 hg N
liquidarea=(N*((hg*0.5*(w1+w2))+(0.5*w1*rc*cos(shi))-((rc^2)*shi)));
end
%% Resistance
function r3=R3(rc,k)
global wfin hg kl
r3= (0.185*wfin)/(kl(k)*(hg-delta_e(rc)));
end
function r5=R5(rc)
global ks w2
r5= delta_e(rc)/(ks*w2);
return
end
function r6=R6(k)
global wfin ks
r6= delta_film(k)/(ks*wfin);
end
function r9=R9(rc,k)
global kl w1
r9= delta_c(rc)/(kl(k)*w1);
return
end
function deltaf= delta_film(k)
global meul wfin Qin rv Lc N sigmal rowl hfg
deltaf= ((64*meul(k)*wfin*Qin*rv)/(Lc*N*sigmal(k)*rowl(k)*hfg(k)))^(0.333);
end
function deltaec1=delta_c(rc)
global hg
deltaec1= hg-(rc1(rc)*(1-(cos(psi(rc1(rc))))));
end
function deltaec=delta_e(rc)
global hg
deltaec= hg-(rc41(rc)*(1-(cos(psi(rc41(rc))))));
end
function rc=rc1(rc)
rc=rc(1)/1000;
end
function rcc=rc41(rc)
rcc= rc(41)/1000;
end
%% Definign the rcmin as a function
function al=alpha(z,rc,k)
rc=[rc(1);rc(41)];
global vew lLt w2 gamma pi w1 ;
alpha=@(rc) pi-real(asin(w1/2./rc))-(2*vew);
alpha_bar1=integral(alpha,1,lLt); % defined the alpha bar
alpha_bar=alpha_bar1/1000;
gamma=33*pi/180;
al=(w2*0.5)/(cos(alpha_bar-gamma)-(tan(-gamma)*(1-sin(alpha_bar-gamma)))); %it is in mm
% psi_rc=abs(real(asin(w1*1000/2./[rc])));
% alpha=pi-psi_rc-2*vew;
% alpha_bar=0.001*trapz(z,alpha);alpha_bar=abs(alpha_bar);
% al=(w2*0.5)/(cosd(alpha_bar-gamma)-(tan(-gamma)*(1-sin(alpha_bar-gamma))));
end
  3 Comments
Show 1 older comment
Cris LaPierre
Cris LaPierre on 6 Jan 2021
Next, your values for rc never change, so you calculation of alpha never changes. If appears you obtain the values of rc from your ode.
[z,rc]=ode45(@(z,rc) (derpvz(z,rc,k)-derplz(z,rc,k))*(-(rc^2)/sigmal(k)),zspan,rc0);
That might be a good place to start your debugging.
ABHISEK PATI
ABHISEK PATI on 7 Jan 2021
thanks for the answer i have checked for myself if you run the code at k=1 and observe the code and again run the code at k=3 and open rc,z curve u can see the rc(1) value will be slightlty upwards giving a higher value than that at k=1.so the rc,z is changing with temperature loop k .
thank you

Sign in to comment.

Sign in to answer this question.

Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

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

Start Hunting!