close all;clear;clc;umatrix = zeros(200, 1);Cmatrix = zeros(200, 1);j = 1; % (1) Get t and di
Pp = 15.5; % Primary pressure
do = 0.01748; % Outer tube diameter
S = 160.65; % Allowable stress
y = 0.4; % Coeffcient given by ASME boiler code
t = Pp * do/(2 * S + 2 * y * Pp); % Wall thickness
di = do - 2 * t; % Inner tube diameter
% (2) With guessed up, get the total number of tube Nsg
mp = (1.38 * 10^4) / 3; % Mass flow rate into one SG
rhop = 706.917; % Primary density, T = 309.055 C and P = 155 bar
for up = 0.1:0.1:20 % To find the optimal velocity
Nsg = 1.2 * mp / (rhop * up * di^2 * pi / 4); % Number of SG tubes
% (3) Get Tpout, Rpmf, q''(z=0), q''(z=Lsg)
% To calculate Rpmf
% Calcuation Rp
kp = 0.54471; % Thermal conductivity of primary system
Ap = pi * di^2 / 4; % Flow area of one tube of primary system
mup = 84.9498 * 10^(-6); % Viscosity of primary system
cp = 5.71027 * 10^3; % Heat capacity of primary system
hp = (kp / di) * 0.023 * ((mp * di) / (Ap * Nsg * mup))^(0.8) * (cp * mup / kp)^(0.33); % HTC of primary system
Rp = (do / di) * (1 / hp); % Primary heat transfer resistance
% Calculation Rw
kw = 12.552; % Wall thermal conductivity of Inconel 690
Rw = (do * log(do / di)) / (2 * kw); % Tube wall resistance
% Calculation Rf
Rf = 1.057 * 10^(-5); % It is given
% Calculation Rs
hf = 1249.38 * 1000; % Enthalpy of fluid
hg = 2777.06 * 1000; % Enthalpy of gas
hfg = hg - hf; % Unit : J/kg
Prf = 0.855144 % Prandtl number of fluid
cf = 5.33358 * 1000; % Heat capacity of fluid, unit : J/kg
Cws = 0.006; % It is given muf = 92.5798 * 10^(-6); % Dynamic liquid viscosity
sigma = 0.0184373; % Surface tension
g = 9.80665; % Gravitational acceleration
rhof = 746.02; % Density of liquid
rhog = 34.4984; % Density of gas
kf = 0.577425; % Thermal conductivity of the liquid
F = (hfg / cf) * Cws * (1 / (muf * hfg) * ((sigma / (g * (rhof - rhog))^(1/2)))^(1/3)) * ((cf * muf / kf)^(1.7)); % Rohsenow correlation
Tp1 = 326.61; % Tin
Tp2 = 291.50; % Tout
Ts = 282.379; % Saturation temperature at P = 6.65MPa
Rc = Rp + Rw + Rf; % q''(z=0), q''(z=Lsg)
syms Q1 Q2 Q1 = solve((Tp1 - Ts) == Rc * Q1 + F * Q1^(1/3), Q1); % Q1 = q''(z=0)
Q2 = solve((Tp2 - Ts) == Rc * Q2 + F * Q2^(1/3), Q2); % Q2 = q''(z=Lsg)
% (4) Get Asg and Lsg
Asg = mp * cp * (Rc * log(Q1 / Q2) + (1/2) * F * (Q2^(-2/3) - Q1^(-2/3))); % SG wall area
Lsg = (Asg / (Nsg * do * pi)); % Length of a tube
% (5) Get pumping power
epsilon = 1.01 * 10^(-5); % SG tube surface roughness
K = 3; % Total form loss coefficient
mp = Ap * Nsg * rhop * up; % mp = rhoAg
f = 0.0055 * (1 + ((2 * 10^4 * epsilon / di) + (10^6 * Ap * mup / ((mp / Nsg) * di)))^(0.333)); % Rough condition
Ploss = (((K + ((f * Lsg) / di)) * up^2 * rhop) / 2); Q = mp / rhop; E = Ploss * Q / Asg; % Pumping power per unit area
% (6) Get yearly cost during lifetime
Cao = 945.86; % Cost per unit area(based on do) of reference SG
Asgo = 4.645 * 10^3; % Heat exchanger area(based on do) of reference SG
Cp = 5.33 * 10^(-9); % Cost per unit energy
Ca = Cao * (Asg / Asgo)^(0.6); % Six-tenth rule, cost per unit area
i = 0.05; % Interest rate
n = 30; % SG lifetime
R = i * (1 + i)^(n) / ((1 + i)^(n) - 1); % Capital recovery factor
CF = 0.88; % Capacity factor
theta = 365 * 24 * 3600 * CF; % Pump operation time
. C = Asg * (Ca * R + E * theta * Cp); % yearly cost during lifetime
umatrix(j, 1) = up; Cmatrix(j, 1) = C; j = j + 1;end plot(umatrix, Cmatrix);[mincost, I] = min(Cmatrix);disp(mincost);disp(umatrix(I, 1));I made the MATLAB code described above, but it gives Invalid use of operator error and doesn't progress the for loops. How can I solve this problem?
Best Answer