Open Auditore3009 opened 4 months ago
Auditore3009
commented
4 months ago The following code wss used to find the pressure at saturation temperature conditions in the middle of the thermal power plant cycle.
function psat = satpres(input) p1 = input(1); t1 = input(2); s1 = XSteam('s_pT',p1,t1); p = p1; dp = 0.1; error = 5; tol = 1e-2; while error > tol s2 = XSteam('sV_p',p); error = s1 - s2; p = p - dp; end psat = p; end
Auditore3009
commented
4 months ago The following code was used to find the mass of the steam bled to the Reheaters.
%mWhp1 + (m-m1)Whp2 + (m-m1-m2)WIP1 + (m-m1-m2-m3)WIP2 + (m-m1-m2-m3-m4)WIP3 + (m-m1-m2-m3-m4)WLP11 + ((m-m1-m2-m3-m4-m5)/2)WLP12 + %(((m-m1-m2-m3-m4-m5)/2)-m6)WLP13 + (((m-m1-m2-m3-m4-m5)/2)-m6-m8)WLP21 + (((m-m1-m2-m3-m4-m5)/2)-m7)WLP23
N=1000; %number of samples to brute force search
xsamples=linspace(-pi,pi,N);
fsamples=f(xsamples);
[fmax,imax]=max(fsamples);
xmax=xsamples(imax);
close all; clear all; clc;
%% Total Power Output given Pout = 1.096850e+06; %% KW
%% Known Temperature in [C] T1 = 600; T4 = 600;
%% Constant Variables ## Values are based on trial and error for the given data EffHPT1 = 0.75; EffHPT2 = 0.42; EffIPT1 = 0.62; EffIPT2 = 0.73; EffIPT3 = 0.8; EffLPT11 = 0.76; EffLPT12 = 0.58; EffLPT13 = 0.335; EffLPT21 = 0.765; EffLPT22 = 0.34; Eff = [EffHPT1, EffHPT2,EffIPT1,EffIPT2,EffIPT3,EffLPT11,EffLPT12,EffLPT13,EffLPT21,EffLPT22];
%% Known Pressures in [Bar] P1 = 262.5; P2 = 74.94;
P3 = 55.55;
P4 = P3; P5 = 23.64; P6 = 11.65; P7 = 5.89; P810 = 2.41; P820 = 0.63; P811 = 0.25; P910 = 0.0575; P920 = P910;
%% Mass removed for reheater [kg/sec] m1 = 39.98; m2 = 75.82; m3 = 31.28; m4 = 28.20; m5 = 30.74; m6 = 35.60; m7 = 21.27; m8 = 20.96;
%% Calculations for T values for HPT (Points: 1, 2, 3) T2 = ((T1 + 273.15) ((P1 / P2) ^ (-0.4 / 1.4))) -273.15; T21 = (T1 - EffHPT1 (T1 - T2)); T3 = ((T2 + 273.15) ((P2 / P3) ^ (-0.4 / 1.4))) -273.15; T31 = (T21 - EffHPT2 (T21 - T3));
%% HPT Work done [ kJ ] Cp11 = (XSteam('Cp_pT', P1, T1) + XSteam('Cp_pT', P2, T21) ) / 2; Cp12 = (XSteam('Cp_pT', P2, T21) + XSteam('Cp_pT', P3, T31) ) / 2; WHPT1 = Cp11(T1-T21); WHPT2 = Cp12(T21-T31); WHPT = WHPT2 + WHPT1;
%% Calculations for T values for IPT [ points :5, 6, 7] T5 = ((T4 + 273.15) ((P4 / P5) ^ (-0.4 / 1.4))) -273.15; T51 = (T4 - EffIPT1 (T4 - T5)); T6 = ((T51 + 273.15) ((P5 / P6) ^ (-0.4 / 1.4))) -273.15; T61 = (T51 - EffIPT2 (T51 - T6)); T7 = ((T61 + 273.15) ((P6 / P7) ^ (-0.4 / 1.4))) -273.15; T71 = (T61 - EffIPT3 (T61 - T7));
%% IPT Work done [ kJ ] Cp21 = (XSteam('Cp_pT', P4, T4) + XSteam('Cp_pT', P5, T51) ) / 2; Cp22 = (XSteam('Cp_pT', P5, T51) + XSteam('Cp_pT', P6, T61) ) / 2; Cp23 = (XSteam('Cp_pT', P6, T61) + XSteam('Cp_pT', P7, T71) ) / 2; WIPT1 = Cp21(T4-T51); WIPT2 = Cp22(T51-T61); WIPT3 = Cp23*(T61-T71); WIPT = WIPT1 + WIPT2 + WIPT3;
%% Calculations for T values for LPT1 [ points :810, 811, 910] T810 = ((T71 + 273.15) ((P7 / P810) ^ (-0.4 / 1.4))) -273.15; T8101 = (T71 - EffLPT11 (T71 - T810)); T811 = ((T8101 + 273.15) ((P810 / P811) ^ (-0.4 / 1.4))) -273.15; T8111 = (T8101 - EffLPT12 (T8101 - T811)); T910 = ((T8111 + 273.15) ((P811 / P910) ^ (-0.4 / 1.4))) -273.15; T9101 = (T8111 - EffLPT13 (T8111 - T910));
%% LPT1 Work done [ kJ ] Cp31 = (XSteam('Cp_pT', P7, T71) + XSteam('Cp_pT', P810, T8101) ) / 2; Cp32 = (XSteam('Cp_pT', P810, T8101) + XSteam('Cp_pT', P811, T8111) ) / 2; Cp33 = (XSteam('Cp_pT', P811, T8111) + XSteam('Cp_pT', P910, T9101) ) / 2; WLPT11 = Cp31(T71-T8101); WLPT12 = Cp32(T8101-T8111); WLPT13 = Cp33*(T8111-T9101); WLPT1 = WLPT11 + WLPT12 + WLPT13;
%% Calculations for T values for LPT2 [ points :820, 920] T820 = ((T71 + 273.15) ((P7 / P820) ^ (-0.4 / 1.4))) -273.15; T8201 = (T71 - EffLPT21 (T71 - T820)); T920 = ((T8201 + 273.15) ((P820 / P920) ^ (-0.4 / 1.4))) -273.15; T9201 = (T8201 - EffLPT22 (T8201 - T920));
%% LPT2 Work done [ kJ ] Cp41 = (XSteam('Cp_pT', P7, T71) + XSteam('Cp_pT', P820, T8201) ) / 2; Cp42 = (XSteam('Cp_pT', P820, T8201) + XSteam('Cp_pT', P920, T9201) ) / 2; WLPT21 = Cp41(T71-T8201); WLPT22 = Cp42(T8201-T9201); WLPT2 = WLPT21 + WLPT22; %% Total Work Output W = WHPT + WIPT + WLPT1 + WLPT2; Mtot = Pout/W;
%% Calculation of real Value Total Mass Flow Rate syms M1
equation = Pout == M1WHPT1 + (M1-m1)WHPT2...
solutions = solve(equation,M1);
%% Work done by each Turbine M = 620.732496484; Whp1 = MWHPT1 ; Whp2 = (M-m1)WHPT2; Wip1 = (M-m1-m2)WIPT1; Wip2 = (M-m1-m2-m3)WIPT2; Wip3 = (M-m1-m2-m3-m4)WIPT3; Wlp11 = (M-m1-m2-m3-m4)WLPT11; Wlp12 = ((M-m1-m2-m3-m4-m5)/2)WLPT12; Wlp13 = (((M-m1-m2-m3-m4-m5)/2)-m6)WLPT13; Wlp21 = (((M-m1-m2-m3-m4-m5)/2)-m6-m8)WLPT21; Wlp22 = (((M-m1-m2-m3-m4-m5)/2)-m7)WLPT22;
Wtot = Whp1 + Whp2 + Wip1 + Wip2 + Wip3 + Wlp11...
% % % Table with all temperature, pressure, Enthalpy, and mass flow rate values % tableData = [ % T1, P1, XSteam('h_pt', P1, T1), (M1); % T21, P2, XSteam('h_pt', P2, T21), (M1-m1); % T31, P3, XSteam('h_pt', P3, T31), (M1-m1-m2); % T4, P4, XSteam('h_pt', P4, T4), (M1-m1-m2-m3);%IPT1in % T51, P5, XSteam('h_pt', P5, T51), (M1-m1-m2-m3-m4);%IPT1out % T61, P6, XSteam('h_pt', P6, T61), (M1-m1-m2-m3-m4-m5);%IPT2out % T71, P7, XSteam('h_pt', P7, T71), (M1-m1-m2-m3-m4-m5-m6);%IPT3out % T8101, P810, XSteam('h_pt', P810, T8101), (M1-m1-m2-m3-m4-m5-m6-m7);%LPT11,out,actual % T8111, P811, XSteam('h_pt', P811, T8111), (M1-m1-m2-m3-m4-m5-m6-m7-m8);%LPT12,out,actual % T8201, P820, XSteam('h_pt', P820, T8201), (((M-m1-m2-m3-m4-m5)/2)-m7) - m6 - m7 - m8;%LPT21,out,actual % T9101, P910, XSteam('h_pt', P910, T9101), M1; % T9201, P920, XSteam('h_pt', P920, T9201), M2; % ]; % % columnNames = {'Temperature (°C)', 'Pressure (Bar)', 'Enthalpy (kJ/kg)', 'Mass Flow Rate (kg/s)'}; % rowNames = {'T1/P1/h1/M1', 'T21/P2/h3/M2', 'T31/P3/h31/M3',... % 'T4/P4/h4/M4', 'T51/P3/h51/M3', 'T61/P6/h61/M6', ... % 'T71/P7/h71/M7', 'T8101/P810/h8101/M8', 'T8111/P811/h8111/M8',... % 'T8201/P820/h8201/M8', 'T9101/P910/h9101/M9', 'T9201/P920/h9201/M9'}; % % Display the table using the 'uitable' function within a GUI figure % fig = figure; % uitable('Data', tableData, 'ColumnName', columnNames, 'RowName', rowNames, 'Position', [50, 50, 400, 300]);
%%Exergy Analaysis
% Reference state properties (assumed) T_ref = 298; % Reference temperature in Kelvin P_ref = 1; % Reference pressure in bar
%%Exergy Analaysis for HPT
% Assuming you have the specific enthalpies at the inlet and outlet of the Intermediate-pressure turbine h1 = XSteam('h_pt', P1, T1); h21 = XSteam('h_pt', P2, T21); h31 = XSteam('h_pt', P3, T31);
% Calculating specific entropy at the inlet and outlet of the High-pressure turbine s1 = XSteam('s_pT', P1, T1); % Specific entropy at the inlet s21 = XSteam('s_pT', P2, T21); % Specific entropy at the outlet of 1st half of HPT s31 = XSteam('s_pT', P3, T31); % Specific entropy at the outlet of 2nd half of HPT
% Exergy calculation at the inlet and outlet of the high-pressure turbine exin_HPT1 = h1 - T_ref s1; % Exergy at the inlet of 1st half of HPT exout_HPT1 = h21 - T_ref s21; % Exergy at the outlet of 1st half of HPT exout_HPT2 = h31 - T_ref * s31; % Exergy at the outlet of 2nd half of HPT
%Change in Exergy for HPT DE_HPT = (M(exin_HPT1 - exout_HPT1) + (M-m1)(exout_HPT1 - exout_HPT2)) - Whp1 - Whp2;
%The exergy efficiency equation for the HPT Exn_HPT = (exout_HPT1 + exout_HPT2)/(exin_HPT1); %%Exergy Analaysis for IPT
% Calculating the specific enthalpies at the inlet and outlet of the Intermediate-pressure turbine h4 = XSteam('h_pt', P4, T4); h51 = XSteam('h_pt', P5, T51); h61 = XSteam('h_pt', P6, T61); h71 = XSteam('h_pt', P7, T71);
% Calculating the specific entropy at the inlet and outlet of the Intermediate-pressure turbine s4 = XSteam('s_pT', P4, T4); % Specific entropy at the inlet of 1st half of IPT s51 = XSteam('s_pT', P5, T51); % Specific entropy at the outlet of 1st half of IPT s61 = XSteam('s_pT', P6, T61); % Specific entropy at the outlet of 2nd half of IPT s71 = XSteam('s_pT', P7, T71); % Specific entropy at the outlet of 3rd half of IPT
% Exergy calculation at the inlet and outlet of the Intermediate-pressure turbine exin_IPT1 = h4 - T_ref s4; % Exergy at the inlet exout_IPT1 = h51 - T_ref s51; % Exergy at the outlet of 1st half of IPT exout_IPT2 = h61 - T_ref s61; % Exergy at the outlet of 2nd half of HPT exout_IPT3 = h71 - T_ref s71; % Exergy at the outlet of 3rd half of IPT
%Change in Exergy for IPT DE_IPT = ((M-m1-m2)*(exin_IPT1 - exout_IPT1)...
% The exergy efficiency equation for the IPT Exn_IPT = ((exout_IPT1 + exout_IPT2 + exout_IPT3)/(exin_IPT1));
%%Exergy Analaysis for LPT
% Calculating the specific enthalpies at the inlet and outlet of the Low-pressure turbine h8101 = XSteam('h_pt', P810, T8101); h8111 = XSteam('h_pt', P811, T8111); h8201 = XSteam('h_pt', P820, T8201); h9101 = XSteam('h_pt', P910, T9101); h9201 = XSteam('h_pt', P920, T9201);
% Calculating the specific entropy at the inlet and outlet of the Low-pressure turbine %s71 = XSteam('s_pT', P7, T71); % Specific entropy at the outlet of 3rd half of IPT s8101 = XSteam('s_pt', P810, T8101); s8111 = XSteam('s_pt', P811, T8111); s8201 = XSteam('s_pt', P820, T8201); s9101 = XSteam('s_pt', P910, T9101); s9201 = XSteam('s_pt', P920, T9201);
% Exergy calculation at the inlet and outlet of the Intermediate-pressure turbine ex_LPT11_in = h71 - T_ref s71; % Exergy at the inlet ex_LPT11_out = h8101 - T_ref s8101; % Exergy at the outlet of 1st half of IPT ex_LPT12_out = h8111 - T_ref s8111; % Exergy at the outlet of 2nd half of HPT ex_LPT13_out = h9101 - T_ref s9101; % Exergy at the outlet of 3rd half of IPT
%Change in Exergy for LPT1 DE_LPT1 = (M-m1-m2-m3-m4)*(ex_LPT11_in - ex_LPT11_out)...
%The exergy efficiency equation for the LPT1 Exn_LPT1 = (( ex_LPT11_out + ex_LPT12_out + ex_LPT13_out)/(ex_LPT11_in));
% Exergy calculation at the inlet and outlet of the LPT2 ex_LPT21_in = h71 - T_ref s71; % Exergy at the inlet ex_LPT21_out = h8201 - T_ref s8201; % Exergy at the outlet of 1st half of IPT ex_LPT22_out = h9201 - T_ref * s9201; % Exergy at the outlet of 2nd half of HPT
% Wlp21 = (((M-m1-m2-m3-m4-m5)/2)-m6-m8)WLPT21; % Wlp22 = (((M-m1-m2-m3-m4-m5)/2)-m7)WLPT22;
%Change in Exergy for LPT2 DE_LPT2 = (((M-m1-m2-m3-m4-m5)/2)-m6-m8)*(ex_LPT21_in - ex_LPT21_out)...
%The exergy efficiency equation for the LPT2 EXn_LPT1 = ((ex_LPT21_out+ ex_LPT22_out +ex_LPT13_out)/(ex_LPT21_in));