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

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
    

[Cris LaPierre]

I get an error the first time I try to run your code because W is defined using w0, which is not defined until 6 lines later.

[Cris LaPierre]

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]

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