my main program is
Vd=40; %volume of drum m^3
Vr=37; %volume of riser m^3
Vdc=11; %volume of downcomer m^3
Vt=Vd+Vr+Vdc;%total volume
Ad=20; %drum area at normal operating level m^2
md=95000;%mass of drum in kg
mt=300000; %total metal mass kg
mr=160000; %total riser mass kg
k=25;%friction coefficent in downcomer riser loop
B=0.3;%parameter
Td=12; % residence time of steam min drum s
Cp=0.42;%specific heat of steel Kj/kgdegC
n=100;t=zeros(1,n);t(1,1)=0;qf=50; %feed flow rate kg/s(to be specified by user)
p=zeros(1,n);p(1,1)=8.5; %steam pressue MPa(to be specified by user)
Vwt=zeros(1,n);Vwt(1,1)=57.2;%initial total volume of water in m^3
Vst=zeros(1,n);Vst(1,1)=Vt-Vwt(1,1);pf=9.142235;%pressure at which feed is entering Mpa(to be specified by user)
hf=1026.33;%-0.0003*pf^6+0.019*pf^5-0.47*pf^4+6.1*pf^3-46*pf^2+2.4e+002*pf+5.7e+002;%enthalpy of saturated feed kj/kg
hs=zeros(1,n);hs(1,1)=0.059*p(1,1)^3-2.3*p(1,1)^2+10*p(1,1)+2.8e+003; %specific enthalpy of steam (kJ/kg)
hw=zeros(1,n);hw(1,1)=0.24*p(1,1)^3-7.5*p(1,1)^2+1.2e+002*p(1,1)+7.1e+002;%specific enthalpy of water (kJ/kg)
Ps=zeros(1,n);Ps(1,1)=-0.0029*p(1,1)^6+0.11*p(1,1)^5-1.8*p(1,1)^4+15*p(1,1)^3-64*p(1,1)^2+1.5e+002*p(1,1)-1.3e+002; % density of steam(kg/m^3)
Pw=zeros(1,n);Pw(1,1)=-0.056*p(1,1)^3+1.6*p(1,1)^2-32*p(1,1)+9e+002; %density of water(kg/m^3)
U=zeros(1,n);U(1,1)=-(6*0.0029)*p(1,1)^5+(0.11*5)*p(1,1)^4-(1.8*4)*p(1,1)^3+(15*3)*p(1,1)^2-(64*2)*p(1,1)+002; %dPs/dp
V=zeros(1,n);V(1,1)=-(0.056*3)*p(1,1)^2+(1.6*2)*p(1,1)-32; %dPw/dp
W=zeros(1,n);W(1,1)=(0.059*3)*p(1,1)^2-(2.3*2)*p(1,1)+10; %dhs/dp
X=zeros(1,n);X(1,1)=(0.24*3)*p(1,1)^2-15*p(1,1)+002; %dhw/dp
ts=zeros(1,n);ts(1,1)=0.06*p(1,1)^3-1.9*p(1,1)^2+27*p(1,1)+1.7e+002; %temperature of steam(degC)
O=zeros(1,n);O(1,1)=(0.06*3)*p(1,1)^2-(1.9*2)*p(1,1)+27; %dts/dp
g=9.81; %acceleration due to gravity m/s^2
%first system is at equilibrium
qs=qf;Q=qs*hs(1,1)-qf*hf;%in kW
Vsd0=5.95569;hc=zeros(1,n);hc(1,1)=hs(1,1)-hw(1,1);%enthalpy of condensation in kJ/kg
Vsd=zeros(1,n);Vsd(1,1)=Vsd0-((Td*(hw(1,1)-hf)*qf)/(Ps(1,1)*hc(1,1)));%volume of steam in the drum below water level
Lr=11.739;%length of riser in m
Ldc=11.739;%length of downcomer in m
Adc=Vdc/Ldc;% area of downcomer in m^2
Ar=Vr/Lr;% area of riser in m^2
ar=zeros(1,n);ar(1,1)=0.051;aar=zeros(1,n);aar(1,1)=0.1262;neta=zeros(1,n);neta(1,1)=ar(1,1)*(Pw(1,1)-Ps(1,1))/Ps(1,1);daarp=zeros(1,n);daarp(1,1)=(1/(Pw(1,1)-Ps(1,1))^2)*(Pw(1,1)*U(1,1)-Ps(1,1)*V(1,1))*(1+(Pw(1,1)/(Ps(1,1)*(1+neta(1,1)))))-(((Ps(1,1)+Pw(1,1))*log(1+neta(1,1)))/(neta(1,1)*Ps(1,1)));daarar=zeros(1,n);daarar(1,1)=(Pw(1,1)/(Ps(1,1)*neta(1,1)))*((log(1+neta(1,1))/neta(1,1))-(1/1+neta(1,1)));Vwd=zeros(1,n);Vwd(1,1)=Vwt(1,1)-Vdc-((1-aar(1,1))*Vr);lw=zeros(1,n);lw(1,1)=Vwd(1,1)/Ad;ls=zeros(1,n);ls(1,1)=Vsd(1,1)/Ad;l=zeros(1,n);l(1,1)=lw(1,1)+ls(1,1);qsd=zeros(1,n);qsd(1,1)=(Ps(1,1)*Vsd0/Td);qdc=zeros(1,n);qdc(1,1)=((2*Pw(1,1)*Adc*(Pw(1,1)-Ps(1,1))*g*aar(1,1)*Vr)/k)^0.5;qct=zeros(1,n);qct(1,1)=(hw(1,1)-hf(1,1))*qf/hc(1,1);qr=zeros(1,n);qr(1,1)=qdc(1,1);e11=zeros(1,n);e11(1,1)=Pw(1,1)-Ps(1,1);e12=zeros(1,n);e12(1,1)=(Vwt(1,1)*V(1,1))+(Vst(1,1)*U(1,1));e21=zeros(1,n);e21(1,1)=(Pw(1,1)*hw(1,1))-(Ps(1,1)*hs(1,1));e22=zeros(1,n);e22(1,1)=(Vwt(1,1)*(hw(1,1)*V(1,1)+Pw(1,1)*X(1,1)))+(Vst(1,1)*(hs(1,1)*U(1,1)+Ps(1,1)*W(1,1)))-Vt+(mt*Cp*O(1,1));e32=zeros(1,n);e32(1,1)=((Pw(1,1)*X(1,1)-ar(1,1)*hc(1,1)*V(1,1))*(1-aar(1,1))*Vr)+((((1-ar(1,1))*hc(1,1)*U(1,1))+Ps(1,1)*W(1,1))*aar(1,1)*Vr)+((Ps(1,1)+(Pw(1,1)-Ps(1,1))*ar(1,1))*hc(1,1)*Vr*daarp(1,1))-Vr+(mr*Cp*O(1,1));e33=zeros(1,n);e33(1,1)=((1-ar(1,1))*Ps(1,1)+ar(1,1)*Pw(1,1))*hc(1,1)*Vr*daarar(1,1);e42=zeros(1,n);e42(1,1)=(Vsd(1,1)*U(1,1))+((1/hc(1,1))*(Ps(1,1)*Vsd(1,1)*W(1,1)+Pw(1,1)*Vwd(1,1)*X(1,1)-Vsd(1,1)-Vwd(1,1)+md*Cp*O(1,1)))+(ar(1,1)*(1+B)*Vr*(aar(1,1)*U(1,1)+(1-aar(1,1))*V(1,1)+(Ps(1,1)-Pw(1,1))*daarp(1,1)));e43=zeros(1,n);e43(1,1)=ar(1,1)*(1+B)*(Ps(1,1)-Pw(1,1))*Vr*daarar(1,1);e44=zeros(1,n);e44(1,1)=Ps(1,1);%give step input
qs=qs+10;%Q=Q+10000;
%solving the dynamic equations
for i=1:n; dt=10; param=[e11(1,i),Vwt(1,i),e12(1,i),p(1,i),qf,qs,dt,e21(1,i),e22(1,i),Q,hf,hs(1,i),e32(1,i),e33(1,i),ar(1,i),... hc(1,i),qdc(1,i),e42(1,i),e43(1,i),e44(1,i),Vsd(1,i),Ps(1,i),Vsd0,Td,hw(1,i)]; x0=[Vwt(1,i); p(1,i);ar(1,i); Vsd(1,i)]; f=@(x) myfun(x,param); x=fsolve(f,x0); x(1,2)=p(1,i+1); x(1,3)=ar(1,i+1); x(1,4)=Vsd(1,i+1); t(1,i+1)=t(1,i)+dt; Vst(1,i+1)=Vt-Vwt(1,i+1); hs(1,i+1)=0.059*p(1,i+1)^3-2.3*p(1,i+1)^2+10*p(1,i+1)+2.8e+003; %specific enthalpy of steam (kJ/kg) hw(1,i+1)=0.24*p(1,i+1)^3-7.5*p(1,i+1)^2+1.2e+002*p(1,i+1)+7.1e+002;%specific enthalpy of water (kJ/kg) Ps(1,i+1)=-0.0029*p(1,i+1)^6+0.11*p(1,i+1)^5-1.8*p(1,i+1)^4+15*p(1,i+1)^3-64*p(1,i+1)^2+1.5e+002*p(1,i+1)-1.3e+002; % density of steam(kg/m^3) Pw(1,i+1)=-0.056*p(1,i+1)^3+1.6*p(1,i+1)^2-32*p(1,i+1)+9e+002; %density of water(kg/m^3) U(1,i+1)=-(6*0.0029)*p(1,i+1)^5+(0.11*5)*p(1,i+1)^4-(1.8*4)*p(1,i+1)^3+(15*3)*p(1,i+1)^2-(64*2)*p(1,i+1)+002; %dPs/dp V(1,i+1)=-(0.056*3)*p(1,i+1)^2+(1.6*2)*p(1,i+1)-32; %dPw/dp W(1,i+1)=(0.059*3)*p(1,i+1)^2-(2.3*2)*p(1,i+1)+10; %dhs/dp X(1,i+1)=(0.24*3)*p(1,i+1)^2-15*p(1,i+1)+002; %dhw/dp ts(1,i+1)=0.06*p(1,i+1)^3-1.9*p(1,i+1)^2+27*p(1,i+1)+1.7e+002; %temperature of steam(degC) O(1,i+1)=(0.06*3)*p(1,i+1)^2-(1.9*2)*p(1,i+1)+27; hc(1,i+1)=hs(1,i+1)-hw(1,i+1);%enthalpy of condensation in kJ/kg Vsd(1,i+1)=Vsd0-((Td*(hw(1,i+1)-hf)*qf)/(Ps(1,i+1)*hc(1,i+1))); aar(1,i+1)=(Pw(1,i+1)/(Pw(1,i+1)-Ps(1,i+1)))*(1-(Ps(1,i+1)/((Pw(1,i+1)-Ps(1,i+1))*ar(1,i+1)))*log(1+((Pw(1,i+1)-Ps(1,i+1))*ar(1,i+1))/Ps(1,i+1))); neta(1,i+1)=ar(1,i+1)*(Pw(1,i+1)-Ps(1,i+1))/Ps(1,i+1); daarp(1,i+1)=(1/(Pw(1,i+1)-Ps(1,i+1))^2)*(Pw(1,i+1)*U(1,i+1)-Ps(1,i+1)*V(1,i+1))*(1+(Pw(1,i+1)/(Ps(1,i+1)*(1+neta(1,i+1)))))-(((Ps(1,i+1)+Pw(1,i+1))*log(1+neta(1,i+1)))/(neta(1,i+1)*Ps(1,i+1))); daarar(1,i+1)=(Pw(1,i+1)/(Ps(1,i+1)*neta(1,i+1)))*((log(1+neta(1,i+1))/neta(1,i+1))-(1/1+neta(1,i+1))); Vwd(1,i+1)=Vwt(1,i+1)-Vdc-((1-aar(1,i+1))*Vr); lw(1,i+1)=Vwd(1,i+1)/Ad; ls(1,i+1)=Vsd(1,i+1)/Ad; l(1,i+1)=lw(1,i+1)+ls(1,i+1); qsd(1,i+1)=(Ps(1,i+1)*Vsd0/Td); qdc(1,i+1)=((2*Pw(1,i+1)*Adc*(Pw(1,i+1)-Ps(1,i+1))*g*aar(1,i+1)*Vr)/k)^0.5; qct(1,i+1)=(hw(1,i+1)-hf(1,i+1))*qf/hc(1,i+1)+(1/hc(1,i+1)*dt)*(Ps(1,i+1)*Vst(1,i+1)*W+Pw(1,i+1)*Vwt(1,i+1)*X(1,i+1)-Vt+mt*Cp*O(1,i+1))*(p(1,i+1)-p(1,i)); qr(1,1)=qdc(1,1)-(Vr/dt)*(aar*U+(1-aar)*V+(Pw-Ps)*daarp)*(p(1,i+1)-p(1,i))+((Pw-Ps)*Vr*daarar*(ar(1,i+1)-ar(1,i+1))/dt); e11(1,i+1)=Pw(1,i+1)-Ps(1,i+1); e12(1,i+1)=(Vwt(1,i+1)*V(1,i+1))+(Vst(1,i+1)*U(1,i+1)); e21(1,i+1)=(Pw(1,i+1)*hw(1,i+1))-(Ps(1,i+1)*hs(1,i+1)); e22(1,i+1)=(Vwt(1,i+1)*(hw(1,i+1)*V(1,i+1)+Pw(1,i+1)*X(1,i+1)))+(Vst(1,i+1)*(hs(1,i+1)*U(1,i+1)+Ps(1,i+1)*W(1,i+1)))-Vt+(mt*Cp*O(1,i+1)); e32(1,i+1)=((Pw(1,i+1)*X(1,i+1)-ar(1,i+1)*hc(1,i+1)*V(1,i+1))*(1-aar(1,i+1))*Vr)+((((1-ar(1,i+1))*hc(1,i+1)*U(1,i+1))+Ps(1,i+1)*W(1,i+1))*aar(1,i+1)*Vr)+((Ps(1,i+1)+(Pw(1,i+1)-Ps(1,i+1))*ar(1,i+1))*hc(1,i+1)*Vr*daarp(1,i+1))-Vr+(mr*Cp*O(1,i+1)); e33(1,i+1)=((1-ar(1,i+1))*Ps(1,i+1)+ar(1,i+1)*Pw(1,i+1))*hc(1,i+1)*Vr*daarar(1,i+1); e42(1,i+1)=(Vsd(1,i+1)*U(1,i+1))+((1/hc(1,i+1))*(Ps(1,i+1)*Vsd(1,i+1)*W(1,i+1)+Pw(1,i+1)*Vwd(1,i+1)*X(1,i+1)-Vsd(1,i+1)-Vwd(1,i+1)+md*Cp*O(1,i+1)))+(ar(1,i+1)*(1+B)*Vr*(aar(1,i+1)*U(1,i+1)+(1-aar(1,i+1))*V(1,i+1)+(Ps(1,i+1)-Pw(1,i+1))*daarp(1,i+1))); e43(1,i+1)=ar(1,i+1)*(1+B)*(Ps(1,i+1)-Pw(1,i+1))*Vr*daarar(1,i+1); e44(1,i+1)=Ps(1,i+1);end
my function for fsolve is myfun.m
function F=myfun(x,param)e11(1,i)=param(1);Vwt(1,i)=param(2);e12(1,i)=param(3);p(1,i)=param(4);qf=param(5);qs=param(6);dt=param(7);e21(1,i)=param(8);e22(1,i)=param(9);Q=param(10);hf=param(11);hs(1,i)=param(12);e32(1,i)=param(13);e33(1,i)=param(14);ar(1,i)=param(15);hc(1,i)=param(16);qdc(1,i)=param(17);e42(1,i)=param(18);e43(1,i)=param(19);e44(1,i)=param(20);Vsd(1,i)=param(21);Ps(1,i)=param(22);Vsd0=param(23);Td=param(24);hw(1,i)=param(25);F=[e11(1,i)*(x(1)-Vwt(1,i))+e12(1,i)*(x(2)-p(1,i))-(qf-qs)*dt;... e21(1,i)*(x(1)-Vwt(1,i))+e22(1,i)*(x(2)-p(1,i))-(Q+qf*hf-qs*hs(1,i))*dt;... e32(1,i)*(x(2)-p(1,i))+e33(1,i)*(x(3)-ar(1,i))-(Q-ar(1,i)*hc(1,i)*qdc(1,i))*dt;... e42(1,i)*(x(2)-p(1,i))+e43(1,i)*(x(3)-ar(1,i))+e44(1,i)*(x(4)-Vsd(1,i))-(Ps(1,i)*(Vsd0-Vsd(1,i))/Td)*dt-((hf-hw(1,i))/hc(1,i))*dt];endthe error shown is??? Subscript indices must either be real positive integers or logicals.
Error in ==> myfun at 2 e11(1,i)=param(1);
Error in ==> @(x)myfun(x,param)
Error in ==> fsolve at 254 fuser = feval(funfcn{3},x,varargin{:});
Error in ==> drumboiler at 120 x=fsolve(f,x0);
Caused by: Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.
Please suggest a solution
Best Answer