Ocean Thermal Energy Conversion
Programming Language: MatLab
Dependencies:
- XSteam
- Cooldrop--> Python
%givens
T_SW =27; % Celsius (Temperature of the warm surface water)
P_SW = 1; % atm (Pressure of warm surface water)
T6 = 4; % Celsius (Temperature of the cold deep water)
x1 = 0; % dimensionless (quality)
TTD = 3; % Celsius (Terminal Temperature Difference)
P_net = 20; % MW (total net power output of the OTEC power cycle system)
C_p = 4.18; % (kJ/kg*K)
zeta = 0.9; % Dimensionless (Efficiency of the work generator)
rho3 = 1000; % kg/m^3 (Density of Water)
% Asssumed Givens:
T1 = 28; % Celsius
T2 = 20; % Celsius (Temperature of the Flash Evaporator)
T3 = T2;
T4 = T2;
T_CW = 2; % Celsius
T7 = T6 + TTD; % Celsius
p_eff = 0.64; % Dimensionless unit (Assuming the pump efficiency of all pumps to be
equivalent)
%givens for piping system
pi = 3.141592653589793238462643; % pi constant
rho6 = 999.97; % (kg/m^3) (Density of Water at 4 C from engineeringtoolbox.com)
rho1 = 996.53; % (kg/m^3) (Density of Water at 27 C from engineeringtoolbox.com)
g = 9.81; % (Gravitational Constant in meters per seconds squared)
D_p1 = 12; % inches (Diameter of cold deep ocean water pipe)(Assumed value from engineering toolbox.com)
D_p2 = 8; % inches (Diameter of pipe warm surface water pipe)(Assumed value from engineering toolbox.com)
D_p3 = 6; % inches (Diameter of the fresh water pipe)(Assumed value from engineering toolbox.com)
z2 = 600;% Meters (Pipe Length of the cold water deep ocean)
Z2 = 5; % Meters (Pipe Length of the warm surface water)
Mu= 1.5705*(10^-3); % (Pa*s) (Dynamic Viscosity of Seawater at 4 degrees Celsius from engineeringtoolbox.com)
R_p = 0.002; % (mm) (Relative Roughness Factor from
Vi = 0; % m/s
z1 = 0; % m
L_pipes1 = Z2;
L_pipes2 = z2;
K_l = 0;
% Givens for the Industrial Cooling load:
T_load = 0; % Celsius (Required Cooling Load Temperature)
L_annual = 25; % (kWh/ft^2) (Annual Base Cooling Electric Load)
A_tot = 500000; % ft^2 (Total Area of the cooling load)
zeta_c = 0.85; % dimensionless (Compressor Isentropic Efficiency)
T_ambient = 30; % Celsius (Ambient Temperature)
hrs = 24 * 365; % hrs/year
% Financial Givens:
Price_electricity = 0.47; % $/kWh (Price of electricity)
Price_desalinated = 2.00; % $/m^3 (Price of desalinated water)
Cost_initial = 10000; % $/kWnet (Powerplant's initial costs)
M_pipe = 49.55; % lbf/ft (T4-12 ratio of weight per foot of Schedule 40 Stainless Steel Pipe from metalsdepot.com)
Cost_steel = 36.86; % $/m (Cost of steel pipe depending on pipe diameter)
%Solution:
%Conversions:
W_net = P_net * 1000; % kW -MW to kW (Total Work of the OTEC Plant)
D_pipes1 = D_p1 / 39.3700787; % m -inches to meters (Diameter of the pipe for Cold water pump)
D_pipes2 = D_p2 / 39.3700787; % m -inches to meters (Diameter of the pipes for Warm surface water pump)
D_pipes3 = D_p3 / 39.3700787; % m -inches to meters (Diameter of the pipe for the freshwater)
R_pipes = R_p / 1000; % m - millimeters to meters (Relative Roughness of stainless steel pipe)
Mp = (M_pipe * 2.2) / (3.2808399); % kg/m -pounds to kilograms, -feet to meters (Weight of the T4-12 per foot)
% Solution to find all of the OTEC's Enthalpy:
hSW = XSteam('hL_T', T_SW);
h1 = XSteam('hL_T', T1);
h2 = h1;
h2f = XSteam('hL_T', T2);
P2= XSteam('psat_T',T2);
P3 = P2;
h4 = h2f;
h2g = XSteam('hV_T', T2);
h3 = h2g;
x2= (h2-h2f)/(h2g-h2f);
T5 = T7 + TTD; % Celsius
T8 = T5;
h5f = XSteam('hL_T', T5);
h5g = XSteam('hV_T', T5);
s3 = XSteam('sV_T', T3);
s5s = s3;
P5 = XSteam('psat_T',T5);
s5f = XSteam('sL_T', T5);
s5g = XSteam('sV_T', T5);
x5 = (s5s-s5f)/(s5g-s5f);
h5s = h5f + x5*(h5g-h5f);
h5 = zeta * h5g;
h8 = h5f;
hCW = XSteam('hL_T', T_CW);
h6 = XSteam('hL_T', T6);
s8 = XSteam('sL_T', T8);
s9 = s8;
h9 = XSteam('h_ps',P5, s9);
h4 = XSteam('hL_p', P5);
s2f = XSteam('sL_T', T4);
s4 = s2f;
h10 = XSteam('h_ps',P5, s4);
h7 = XSteam('hL_T', T7);
% Solving for Turbine efficiency:
zeta_t = (h3 - h5s) / (h3 - h5); % dimensionless
% Solving for m_dot_3:
m_dot_3 = ((W_net) / (zeta*zeta_t*(h3 - h5s)));
% Solving for the remaining m_dots except m_dot_6:
m_dot_SW = m_dot_3 / x2; % (kg/s)
m_dot_4 = m_dot_SW - m_dot_3; % (kg/s)
m_dot_5 = m_dot_3; % (kg/s)
m_dot_8 = m_dot_5; % (kg/s)
m_dot_9 = m_dot_8; % (kg/s)
% Solving for m_dot_6 and m_dot_7:
W_turb = m_dot_3*(h3-h5s);
% Q_condenser = W_turb --> substitute in the equation below
% Q_condenser = m_dot_6*Cp*(T7 - T6) --> rearranging the equation to solve
% for m_dot_6
% m_dot_6 = (Q_condenser / (Cp*(T7-T6)))
Q_condenser = m_dot_3*(h3 - h5); % kW
m_dot_6 = (Q_condenser / (C_p*(T7-T6))); % (kg/s)
m_dot_7 = m_dot_6; % (kg/s)
% Solving for Work of Turbine:
W_turb = m_dot_3*(h3-h5s); % kW
% Solving for the Pressure Difference of the pipe:
% Tabulating for area of the pipe:
% A = (pi/4)*(D^2)
A1 = (pi/4)*(D_pipes1^2); % m^2 (Area of Pipe for the deep cold water pump)
A2 = (pi/4)*(D_pipes2^2); % m^2 (Area of Pipe for the surface water pump)
A3 = (pi/4)*(D_pipes3^2); % m^2 (Area of Pipe for the Freshwater)
% Tabulating the Volumetric Flow rate of Pumps 1 and 3:
% m_dot = V_dot*rho
V_dot_1 = m_dot_SW / rho1;% (m^3/s) (Volumetric Flow Rate for the surface water pump)
V_dot_6 = m_dot_6 / rho6; % (m^3/s) (Volumetric Flow Rate for the cold deep water pump)
V_dot_8 = m_dot_8 / rho3; % (m^3/s) (Volumetric Flow Rate for the Freshwater pump)
V_dot_9 = V_dot_8
% Tabulating the velocity of Pumps 1 and 3:
% V_dot = A*V
V1 = V_dot_1/A2; % m/s (Velocity for the surface water pump)
V6 = V_dot_6/A1; % m/s (Velocity for the cold deep water pump)
% Calculation of Reynold's Number:
% Re = ((rho*V*D)/Mu)
Re1 = ((rho1*V1*D_pipes2)/Mu); % Dimensionless Unit
Re6 = ((rho6*V6*D_pipes1)/Mu); % Dimensionless Unit
% Finding the Relative Roughness:
releps1 = R_pipes/D_pipes2; % dimensionless unit
releps6 = R_pipes/D_pipes1; % dimensionless unit
% Finding the Friction Factor:
f1=Haaland(releps1, Re1);
f2=Haaland(releps6, Re6);
% Finding the Major Head Losses of the Pipes:
% H_major_loss=((f*L*(V^2))/(D*2*g)
H_major_loss1 = (f1*L_pipes1*(V1^2) / (D_pipes1 * 2 *g)); % m
H_major_loss2 = (f2*L_pipes2*(V6^2) / (D_pipes2 * 2 *g)); % m
% Finding the Pressure difference:
% deltaP = P1 - P2
deltaP1 = P_SW - P2; % Pa (Pressure difference of the surface water pump)
deltaP3 = P_SW; % Pa (Pressure difference for the cold deep surface water pump)
% Finding the Head Losses:
% (P1/(rho*g))+((V1^2)/(2*g)+z1=(P2/(rho*g))+((V2^2)/(2*g))+z2+H_major_loss+H_minor_loss
;z1=0,H_major_loss=((f*L*(V1^2))/(D_pipes*2*g),H_minor_loss=K_L*.5*rho*(V1^2)=0
% Solve for H
% H = ((V1^2)-(V2^2)/(2g)) + ((deltaP)/(rho*g)) + (z1 - z2)H_loss1 = (((V1^2)-(Vi^2))
/ (2*g)) + ((deltaP1)/(rho1*g)) + (z2 - z1); % m (Head loss of the pipe of the surface water)
H_loss2 = (((V6^2)-(Vi^2))/ (2*g)) + ((deltaP3)/(rho6*g)) + (Z2 - z1)
; % m (Head loss of the pipe of the surface water)
% Finding C for the equation H=C*V_dot^2
% C =(Sum(K_l) + f*(L/d)*)*(2*g*(A^2))Q^2
C_warm = ((K_l) + f1*(L_pipes1/D_pipes1))*(2*g*(A2^2));
C_cold = ((K_l) + f2*(L_pipes2/D_pipes2))*(2*g*(A1^2));
%Graphing H-Q Curve for Warm Surface Pump
Q = 0:1:25;
H = C_warm*Q.^2;
figure
plot(Q,H);, xlabel('Flow Rate(Q) (m^3/s)'), ylabel('Head(H) (m)'), title('Warm
Surface Pump H-Q Curve'), grid on
%Graphing H-Q Curve for Cold-Deep Water Pump
Q = 0:1:25;
H = C_cold*Q.^2;
figure
plot(Q,H);, xlabel('Flow Rate(Q) (m^3/s)'), ylabel('Head(H) (m)'), title('Cold Deep
Pump H-Q Curve'), grid on
% Solving for Work of Turbine:
W_turb = m_dot_3*(h3-h5s) % kW
% Solving for the Work of all pumps:
W_pump1 = p_eff * ((deltaP1 * V_dot_1) * 1000) % kW
W_pump2 = p_eff * (m_dot_4*((P_SW/rho3) - (P2/rho3))) % kw
W_pump3 = p_eff * ((deltaP3 * V_dot_6) * 1000) % kW
W_pump4 = p_eff * (m_dot_8*(P_SW/rho3) - (P5/rho3)) % kW
% Verification of the OTEC Power Plant power output is 20MW:
Wnet = W_turb - (W_pump1 + W_pump2 + W_pump3 + W_pump4) % kW
% Finding the Compressor Power:
W_compressor = ((L_annual * A_tot) / hrs); % kW
% Mechanical Vapor Cycle Configuration 1: Air Cooled
TTD1 = 1; % Celsius
gas = 'Air';
TH = T_ambient + 273.15; % K (Conversion from Celsius to Kelvin)
TC = T_load + 273.15; % K (Conversion from Celsius to Kelvin)
% Carnot Coefficient of Performance;
% COP = ((TC) / (TH-TC)) Dimensionless unit
COP1 = ((TC) / (TH-TC));
L_load1 = COP1 * (zeta_c *W_compressor); % kW
% Mechanical Vapor Cycle Configuration 2: Water Cooled before the condensor of the OTEC
gas = 'Water';
TH = T6 + TTD; % Celsius
T1= TH +273.15; % Kelvin
P1=py.CoolProp.CoolProp.PropsSI('P','T',T1,'Q',1,gas);
h1=py.CoolProp.CoolProp.PropsSI('H','T',T1,'Q',1,gas);
s1=py.CoolProp.CoolProp.PropsSI('S','T',T1,'Q',1,gas);
P2=10e6; % Pa
beta=P2/P1;
s2is=s1;
h2is=py.CoolProp.CoolProp.PropsSI('H','P',P2,'S',s2is,gas);
eta_is=0.8;
h2=(h2is-h1)/eta_is+h1;
T2=py.CoolProp.CoolProp.PropsSI('T','H',h2,'P',P2,gas);
s2=py.CoolProp.CoolProp.PropsSI('S','H',h2,'P',P2,gas);
T3=(T_ambient + TTD) +273.15;
P3=P2;
h3=py.CoolProp.CoolProp.PropsSI('H','T',T3,'P',P3,gas);
s3=py.CoolProp.CoolProp.PropsSI('S','T',T3,'P',P3,gas);
s4=s3;
P4 = P1;
h4=py.CoolProp.CoolProp.PropsSI('H','S',s4,'P',P1,gas);
T4=py.CoolProp.CoolProp.PropsSI('T','S',s4,'P',P1,gas);
Q4=py.CoolProp.CoolProp.PropsSI('Q','S',s4,'P',P1,gas);
COP2 = (h2-h3)/(h2-h1)
L_load2 = COP2 * (zeta_c *W_compressor) % kW
% Plot P-h diagram
figure;
plot([h1, h2, h3, h4, h1], [P1, P2, P3, P4, P1], '-o');
xlabel('Enthalpy (h)');
ylabel('Pressure (P)');
title('P-h Diagram');
legend('Cycle', 'Location', 'Best');
grid on;
% Mechanical Vapor Cycle Configuration 2: Ammonia Cooled after the condensor of the OTEC
gas = 'Ammonia';
TH = T7 + TTD; % Celsius
T1= TH +273.15; % Kelvin
P1=py.CoolProp.CoolProp.PropsSI('P','T',T1,'Q',1,gas);
h1=py.CoolProp.CoolProp.PropsSI('H','T',T1,'Q',1,gas);
s1=py.CoolProp.CoolProp.PropsSI('S','T',T1,'Q',1,gas);
P2=10e6; % Pa
beta=P2/P1;
s2is=s1;
h2is=py.CoolProp.CoolProp.PropsSI('H','P',P2,'S',s2is,gas);
eta_is=0.8;
h2=(h2is-h1)/eta_is+h1;
T2=py.CoolProp.CoolProp.PropsSI('T','H',h2,'P',P2,gas);
s2=py.CoolProp.CoolProp.PropsSI('S','H',h2,'P',P2,gas);
T3=(T_ambient + TTD) +273.15;
P3=P2;
h3=py.CoolProp.CoolProp.PropsSI('H','T',T3,'P',P3,gas);
s3=py.CoolProp.CoolProp.PropsSI('S','T',T3,'P',P3,gas);
s4=s3;
P4 = P1;
h4=py.CoolProp.CoolProp.PropsSI('H','S',s4,'P',P1,gas);
T4=py.CoolProp.CoolProp.PropsSI('T','S',s4,'P',P1,gas);
Q4=py.CoolProp.CoolProp.PropsSI('Q','S',s4,'P',P1,gas);
COP3 = (h2-h3)/(h2-h1)
L_load3 = COP3 * (zeta_c *W_compressor) % kW
% Plot P-h diagram
figure;
plot([h1, h2, h3, h4, h1], [P1, P2, P3, P4, P1], '-o');
xlabel('Enthalpy (h)');
ylabel('Pressure (P)');
title('P-h Diagram');
legend('Cycle', 'Location', 'Best');
grid on;
% Economic Analysis:
% Income generated from the OTEC Power Plant
Income1 = Price_electricity * Wnet; % $/hour
Income_1 = Income1 * hrs % $/year (Income of the electricity per Year)
Income2 = Price_desalinated * V_dot_9 % $/s (Income of Desalinated water per second)
Income_2 = Income2 * (3600 * 24 * 365) % $/year (Income of Desalinated water per year)
% Cost from Each MVC configurations:
Load_base_cost1 = (L_load1 * Price_electricity) * hrs % $/year
Load_base_cost2 = (L_load2 * Price_electricity) * hrs % $/year
Load_base_cost3 = (L_load3 * Price_electricity) * hrs % $/year
% Savings comparing each MVC configuration:
Savings2 = Load_base_cost2 - Load_base_cost1 % $/year
Savings3 = Load_base_cost3 - Load_base_cost1 % $/year
% Total Income and Saving of the OTEC power plant per year:
Total_Revenue1 = Income_1 + Income_2 % $/year
Total_Revenue2 = Income_1 + Income_2 + Savings2 % $/year
Total_Revenue3 = Income_1 + Income_2 + Savings3; % $/year
% Total Cost of the OTEC Power Plant:
% Total Cost of Material
Wt_total= Mp*(L_pipes1 + L_pipes2); % kg
Total_cost= Wt_total*Cost_steel; % $
Operation_cost = Cost_initial * 10000 * 20; % $/year
Total_expenditures1 = Operation_cost + Total_cost + Load_base_cost1; % $/year
Total_expenditures2 = Operation_cost + Total_cost + Load_base_cost2; % $/year
Total_expenditures3 = Operation_cost + Total_cost + Load_base_cost3; % $/year
% Total profit of otec
Profit1 = Total_Revenue1 - Total_expenditures1 % $/year
Profit2 = Total_Revenue2 - Total_expenditures2 % $/year
Profit3 = Total_Revenue3 - Total_expenditures3 % $/year
%plot income vs expenditure graph
x = 0:1:25
y1 = Total_Revenue1*x - Total_expenditures1
y2 = Total_Revenue2*x - Total_expenditures2
y3 = Total_Revenue3*x - Total_expenditures3
y4 = Total_expenditures1
figure
plot(x, y1, x, y2, x, y3, x, y4), xlabel('Years'), ylabel('Dollars'), title('Income per Year')
% Determining the most economical pipe diameter
function f=Haaland(releps, Re)
answer=-1.8*log10((releps/3.7)^1.11+6.9/Re);
f= (1/answer)^2;
end