My Matlab code stuck in infinite loop, error comes from line 157 to 242 (fourth while loop). I really don't know what is wrong. I would appreciate any help on this!
clear allclear variablesglobal nos noc nof CompMW CompSigma FeedTemp FeedPres Ftotal ... length ir functiontemperrormax=1e-6;maxiterations=500;nos=100;noc=2;z(1)=0.21;z(2)=0.79;Q(2)=1.04e-6;selectivity=5.85;Q(1)=Q(2)*selectivity;FeedPres=150*76/14.7; %psia to cmHg
PermPres=14.7*76/14.7;FeedTemp=25+273.15;CompMW(1)=32; %O2 (g/mol)
CompMW(2)=28.02; %N2
CompMW(3)=44.01; %CO2
CompSigma(1)=3.433; %Angstroms
CompSigma(2)=3.667;CompSigma(3)=3.996;Ptotal=zeros(1,nos); %Pre-allocations for faster code execution
Rtotal=zeros(1,nos);R=zeros(nos,noc);P=zeros(nos,noc);Sweep=zeros(1,noc);R_temp=zeros(nos,noc);P_temp=zeros(nos,noc);nof=1000; %number of fibers
length=50; %cm
or=0.04/2; %cmnow=6;fmin=nof*2*pi*or*length*Q(1)*FeedPres;for n=1:now Ftotal(n,1)=fmin;endPresDropSet=0.04; %psia
PresDropMax=0.2;x(1)=-2.3506;x(2)=-1.33584;x(3)=-0.43607;x(4)=0.43607;x(5)=1.33584;x(6)=2.35060;w(1)=0.00453;w(2)=0.15706;w(3)=0.724629;w(4)=0.724629;w(5)=0.157067;w(6)=0.00453;id_avg=0.021;id_std=0.068; % percent
pcounter=1;while PresDropSet<=PresDropMax error_f=0.1; f_iterations=0; while error_f>1e-4 && f_iterations<100 error_f=0; for n=1:now id=id_avg*(1.4142*id_std*x(n)+1); ir=id/2; %cm TotalArea=2*pi*or*length*nof; A=TotalArea/nos; for m=1:2 for i=1:noc F(i)=z(i)*Ftotal(n,m); end for j=1:nos %Do crossflow calculations
Ptotal(j)=0; Rtotal(j)=0; for i=1:noc %Initialize with permeate vacuum
if j==1 P(j,i)=Q(i)*A*FeedPres*F(i)/Ftotal(n,m); R(j,i)=F(i)-P(j,i); else P(j,i)=Q(i)*A*FeedPres*R(j-1,i)/Rtotal(j-1); R(j,i)=R(j-1,i)-P(j,i); end Ptotal(j)=Ptotal(j)+P(j,i); Rtotal(j)=Rtotal(j)+R(j,i); end c_error=0.1; iterations=0; while c_error>errormax && iterations<maxiterations %Iterate crossflow solution
c_error=0; for i=1:noc if j==1 P_new=(Q(i)*A*FeedPres*F(i)/Rtotal(j))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*FeedPres/Ftotal(n,m)); R_new=F(i)-P_new; else P_new=(Q(i)*A*FeedPres*R(j-1,i)/Rtotal(j))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*FeedPres/Rtotal(j)); R_new=R(j-1,i)-P_new; end c_error=c_error+(abs(R_new-R(j,i))+abs(P_new-P(j,i)))/Ftotal(n,m); iterations=iterations+1; P(j,i)=P_new; R(j,i)=R_new; end Ptotal(j)=0; Rtotal(j)=0; for i=1:noc Ptotal(j)=Ptotal(j)+P(j,i); Rtotal(j)=Rtotal(j)+R(j,i); end end end Sweeptotal=0; %Sum permeates from last stage to first
for i=1:noc Sweep(i)=0; end for j=nos:-1:1 if j==nos Ptotal(j)=Ptotal(j)+Sweeptotal; else Ptotal(j)=Ptotal(j)+Ptotal(j+1); end for i=1:noc if j==nos P(j,i)=P(j,i)+Sweep(i); else P(j,i)=P(j,i)+P(j+1,i); end end end %Permeate flow
if m ==1 for j = 1:nos for n = 1:now for i=1:noc Permf(n,i)=P(j,i); end end end end for j=1:nos RetPres(j)=FeedPres; end d_maxiterations=500; c_error=0.1; %Use direct-substitution method
d_iterations=0; while c_error>errormax && d_iterations<=d_maxiterations c_error=0; functiontemp=Ftotal(n,m); RetPres=RetPresDrop(R,Rtotal,RetPres); for j=1:nos Ptotal(j) = 0; Rtotal(j) = 0; P_sum(j) = 0; for i=1:noc for n = 1:now P(j,i) = P_sum(j) + 1/sqrt(pi)*(w(n)*Permf(n,i)); %Summing up the flows from each fiber type(n)
end Ptotal(j) = Ptotal(j) + P(j,i); Rtotal(j) = Rtotal(j) + R(j,i); end for i=1:noc y(j,i) = P(j,i)/Ptotal(j); end for i=1:noc if j==1 x(j,i) = F(i)/Rtotal(j); %y(j,i) = P(j,i)/Ptotal(j);
J(j,i) = Q(i)*A*(RetPres(j)*x(j,i)-PermPres*y(j,i)); R_new= F(i) - J(j,i); P_new = P(j+1,i) + J(j,i);% P_new=(Q(i)*A*RetPres(j)*F(i)/Rtotal(j)+P(j+1,i)*(1+Q(i)*A*RetPres(j)/Rtotal(j)))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*RetPres(j)/Rtotal(j));
% R_new=F(i)-(P_new-P(j+1,i));
elseif j==nos x(j,i) = R(j,i)/Rtotal(j); % y(j,i) = Sweep(i)/Ptotal(j);
J(j,i) = Q(i)*A*(RetPres(j)*x(j,i)-PermPres*y(j,i)); R_new= R(j-1,i) - J(j,i); P_new = Sweep(i) + J(j,i);% P_new=(Q(i)*A*RetPres(j)*R(j-1,i)/Rtotal(j)+Sweep(i)*(1+Q(i)*A*RetPres(j)/Rtotal(j)))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*RetPres(j)/Rtotal(j));
% R_new=R(j-1,i)-(P_new-Sweep(i));
else x(j,i) = R(j,i)/Rtotal(j); %y(j,i) = P(j,i)/Ptotal(j); J(j,i) = Q(i)*A*(RetPres(j)*x(j,i)-PermPres*y(j,i)); R_new= R(j-1,i) - J(j,i); P_new = P(j+1,i) + J(j,i);% P_new=(Q(i)*A*RetPres(j)*R(j-1,i)/Rtotal(j)+P(j+1,i)*(1+Q(i)*A*RetPres(j)/Rtotal(j)))/(1+Q(i)*A*PermPres/Ptotal(j)+Q(i)*A*RetPres(j)/Rtotal(j));
% R_new=R(j-1,i)-(P_new-P(j+1,i));
end c_error=c_error+(abs(R_new-R(j,i))+abs(P_new-P(j,i)))/Ftotal(n,m); P(j,i)=P_new; R(j,i)=R_new; end Ptotal(j)=0; Rtotal(j)=0; for i=1:noc Ptotal(j)=Ptotal(j)+P(j,i); Rtotal(j)=Rtotal(j)+R(j,i); end if m==1 if j==nos for i=1:noc Retf(n,i)=R(nos,i); end Rtotaltemp=Rtotal(nos); end% if j == 1
% for i=1:noc
% Permf(n,i)=P(j,i);
% end
% Ptotaltemp=Ptotal(1);
% end
if m==1 dimrec(n,j)=0;% dimperm(n,j) = 0;
for i=1:noc dimrec(n,j)=dimrec(n,j)+R(j,i);% dimperm(n,j) = dimperm(n,j) + P(j,i);
dimdist(j)=j/nos; end dimxfrac(n,j)=R(j,1)/dimrec(n,j);% dimyfrac(n,j)=P(j,1)/dimperm(n,j);
dimrec(n,j)=dimrec(n,j)/Ftotal(n,1);% theta(n,j) = dimperm(n,j)/Ftotal(n,1);
end end end d_iterations=d_iterations+1; end PresDrop(m)=(FeedPres-RetPres(nos))*14.7/76; Ftotal(n,2)=1.1*Ftotal(n,1); end if n==3 xfrac_nv(pcounter)=Retf(3,1)/Rtotaltemp; recovery_nv(pcounter)=Rtotaltemp/Ftotal(3,1); end alpha=0.1*Ftotal(n,1)/(PresDrop(2)-PresDrop(1)); error_f=error_f+abs(alpha*(PresDropSet-PresDrop(1))/Ftotal(n,1)); Ftemp(n)=Ftotal(n,1); Ftotal(n,1)=Ftotal(n,1)+alpha*(PresDropSet-PresDrop(1)); end f_iterations=f_iterations+1; end RetTotal=0; % FeedTotal=0;
%PermTotal= 0;
for i=1:noc Ret(i)=0; %Feed(i)=0;
Perm(i)=0; %Permf_nos(i)=0;
for n=1:now Ret(i)=Ret(i)+1/sqrt(pi)*(w(n)*Retf(n,i)); %Perm(i)=Perm(i)+1/sqrt(pi)*(w(n)*Permf(n,i)); %% total permeate at stage 1
%Feed(i)=Feed(i)+1/sqrt(pi)*(w(n)*Ftemp(n));
end RetTotal=RetTotal+Ret(i); %FeedTotal=FeedTotal+Feed(i);
% PermTotal=PermTotal + Perm(i);
end FeedTotal=0; for n=1:now FeedTotal=FeedTotal+1/sqrt(pi)*(w(n)*Ftemp(n)); end % MBCheck = - PermTotal + FeedTotal - RetTotal;
xfrac_v(pcounter)=Ret(1)/RetTotal; recovery_v(pcounter)=RetTotal/FeedTotal; pcounter=pcounter+1; PresDropSet=PresDropSet*1.1; endfunction RetPres = RetPresDrop(R,Rtotal,RetPres)global nos noc nof CompMW CompSigma FeedTemp FeedPres ... length ir functiontempfor j=1:nos RetMu=0; for i=1:noc CompMu=((FeedTemp*CompMW(i))^0.5)/(CompSigma(i)^2); RetMu=RetMu+CompMu*R(j,i)/Rtotal(j); end RetMu=RetMu*2.6693e-6; if j==1 deltaP=(8*RetMu*functiontemp*(76/FeedPres)*(FeedTemp/273.15)* ... (length/nos))/(pi*(ir^4)*nof); deltaP=deltaP*7.5e-4; %Pa to cmHg
RetPres(j)=FeedPres-deltaP; else deltaP=(8*RetMu*Rtotal(j-1)*(76/RetPres(j-1))*(FeedTemp/273.15)* ... (length/nos))/(pi*(ir^4)*nof); deltaP=deltaP*7.5e-4; %Pa to cmHg RetPres(j)=RetPres(j-1)-deltaP; endendend
Best Answer