Description

Gas Turbine System
Deep Dayaramani
Introduction
In today‘s rapidly development of science and technology, countless products and systems have been invented and combined with others to achieve more promising results. The purpose of this project was to examine such an operation, a hybrid solar fossil-fuel gas turbine system (shown in Fig. 1 and 2 below), and to establish a computational method of its performance.

Figure 1. A S-T cycle illustration of various stages in the gas turbine system

Figure 2. A detailed scheme of the gas turbine system.

Initially, it’s known that the atmospheric pressure of air is P1 = 101 kPa at a temperature T1 = 25 °C. The inlet of air to the system had a flow rate of 6.0 kg/s. Moreover, it’s considered as an ideal gas with an effective molecular mass M of 28.97 kg/kmol.

Task 1 (III)
Given that a mixture of propane (C3H 8) and methane (CH4) was injected into and burned in the burner with air, the molar air-to-fuel ratio α and the fuel propane mole fraction γ were defined:
(Eqn.1) Moreover, the mixture had just enough stoichiometric values that converted all the propane and methane to H2O and CO 2 with no oxygen left over. The balanced equation for arbitrary α and γ greater than or equal to the stoichiometric values was also provided below:
(Eqn.2) And the stoichiometric condition for α was given:
(Eqn.3)
The goal of this task was to calculate the molar air-to-fuel ratio α for specific combustor adiabatic flame (exit) temperature Taf and reactant temperature T r. According to definition, T af is the temperature that balances the combustor energy with no additional heat reactions:
(Eqn.4)
The following equation with α, Tr, T p was shown:
(Eqn.5)
where the molar enthalpy can be obtained from the heat of formation and the specific heat:
(Eqn.6)
In the exhaust gas, the water was assumed to be vapor. And α was isolated but only applied when
Taf< Taf,stoich and α>α stoich:
(Eqn.7)

In order to find α, the heat of formation of different species was first found at reference temperature T。= 25°C in the attached table. Various specific heat were also calculated, at Tavg using Eqn. 6, based on the equation provided in the attached table. Hence, the molar enthalpy of different compounds were obtained.
Further, a function of α could be created according to Eqn.7. In this task, Tp was given to be 1900 K, and Tr was 550 K. The fuel propane mole fraction γ depended on the circumstance given: (i) pure propane, (ii) pure methane, and (iii) a 50% mole fraction mixture. The results were calculated and shown below:
Table 1. The values of the molar air-to-fuel ratio α at Tr = 550 K and T p = 1900 K for different mole fractions of mixture.
Pure propane Pure methane A 50% mole fraction mixture
α 39.5022 15.2750 27.3886
Task 2 (IV)
In this part, we examined the exit condition of the compressor and the work required per kilogram of air. First, we guessed the exit temperature T2 was the product temperature given in Task 1, 1900K. Then, an average temperature of T1 and guessed T 2 was calculated and used to find cp (based on the attached table). For air, its molar enthalpy was given to be 0 at reference temperature. The specific heat at constant volume,cv, was found using the constant relationship:
(Eqn.8) For a compressor, its exit temperature is known to be a function of k (the ratio of cp to c v) ,T1, P 2 , and the efficiency ηcomp:
(Eqn.9) In a similar manner, this exit temperature was now taken to be another guess of the outlet temperature. The step was repeated to measure the real exit temperature. And the difference between the successive T2 values were measured. This whole process was iterated until the difference agreed within 0.2 °C.
Moreover, for the compressor, its work was equal to the changed of enthalpy between 2 states. Since the inlet state was 0, work was the exit stage. The molar enthalpy was measured using Eqn. 6, and then converted from kmol to kg of air by dividing 28.97.
Task 3 (V)
This task was first intended to calculate the mole fraction of all exhaust species with respect to the molar air-to-fuel ratio α. Generally, the mole fraction was obtained by dividing the molar coefficient of the particular species by the total number of moles of all species in the inlet or outlet. In our case, we relied on Eqn.2 provided. The total mole of products (n_ex) was found by summing the molar coefficients of N2, O 2, H 2O, and CO 2. The mole fraction of each species was then calculated as mentioned above. For example, the mole fraction of H2O was
(2+2γ)/n_ex. It’s also given that γ was 0.25 in the mixture. The second objective of this task was to find the specific heat of products, cp,prod, as a function of the mole-fraction weight average of the specific heat of all species involved., as shown in the equation below:
(Eqn.9) The method of finding the specific heat was provided on the attached tables, and Tp was given to be 600 K. At the same temperature, cp of pure air was also measured and compared with cp,prod, and should be smaller than it.At 600K, the cp_prod was found to be 33.0904 kJ/kmol K and cp_air = 30.5943 kJ/kmol K.

Task 4 (VI)
Now that we had the ability to find the temperature at step 2, the cp for air at the average temperature of step 1 and 2, and the specific heat of products at step 4, it was time to combine all of the previous tasks to create a whole cycle.
1. First we started at step 1 and used task 2 to find the Temperature at step 2.
2. We then found the mole fractions of products and molar flow rates using the equations below. Equation 10 gives the ability to us to calculate the mass of the products while equation 11 gives us the ability to calculate the moles of products flowing. Equation 13 gives us is the kmol per second of product gases flowing as a result of burning the propane/methane fuel mixture to achieve the specified α for the air flow of 6.0 kg/s (mass flow rate of air).
(Eqn 10)
(Eqn. 11)
3. Next before we could move on to step 2r and 3, we had to find T5, for taking the regenerator into account. For this we guessed T5 and use a similar program as task 2 to reach a final T5. An average temperature of T4 (1600K) and guessed T 5 was calculated and used to find cp using task 3 and αstoich for the first iteration. The specific heat at constant volume,cv , was found using the constant relationship:
(Eqn.14)
For a compressor, its exit temperature is known to be a function of k (the ratio of cp to c v ) ,T4 , P 2 , and the efficiency ηturbine :
T5 = T4 * (1 + ((PP12) (k−k1) − 1) * ηturbine) (Eqn.15)
In a similar manner, this exit temperature was now taken to be another guess of the outlet temperature. The step was repeated to measure the real exit temperature. And the difference between the successive T5 values were measured. This whole process was iterated until the difference agreed within 0.2 °C.
4. After this T5 was found, it was time to incorporate the regenerator into the big picture. For this we determined the cp_air and cp_products at the average temperature of T2 and T5, and then we used these values for calculating T2r using this formula:
(Eqn. 16)
, where (n*cp)min is the minimum of {(n*cp)air , (n*cp)prod }. This gives us T2r.
5. After this, we find T3 using the solar heat input. We use the first law of thermodynamics to find this Temperature. But we need the cp at the average temperature of T3 and T2r to find T3. For this we created a root_finder and using the two equations of finding cp and T3, we solved for T3. These equations are given below
cp_23 == 28.11 + 0.1967 * 10−2 * (0.5 * (T3 + T2_r)) + 0.4802 * 10−5 * (0.5 * (T3 + T2_r))2 − 1.966 * 10−9 * (0.5 * (T3 + T2_r))3
Eqn(17)
T3 == Qs/(n_dot_air * cp_23) + T2_r (Eqn(18))
6. We then used task 1 and T3, T4 to find a new alpha. After finding the new alpha, we then used this alpha to restart the process from step 1 of Task 4. We repeated this for 2 iterations.
7. After finding the right value of alpha, T5, T2, T2r, T3, molar flow rates of products and air, we find the Power output and the heat input to the combustor along with the overall efficiency of the cycle. We do this using the formulas given below.
W net = W t − W c =
Eqn 19
Qnet = Qb + Qsolar = ndotprod * cp, products * (T 4 − T 3) + Qsolar Eqn 20
Eqn 21 8. We then compared this efficiency to a constant cp efficiency whose formula is given below.

We then ran the computation for γ = 0.25, Q_s = 0, P2 = 500 kPa, T4 = 1600 K, ηcomp = 0.85, ηturb = 0.85, eps_regen = 0.75 and found the results given below.
Power 1.9268e+03 kW
Heat 4.3481e+03 kW
T5_f 1.1826e+03 K
alpha_f 54.6895
W_c 1.2233e+03 kW
W_t 3.1501e+03 kW
n_comp 0.4431
n_constant cp 0.5740
Table 2: Results from running the computation of task 4 for different values.
n_prod was found at a k value of averaged k values of air and product values determined at Tp=0. The efficiency at constant cp, is much higher than the achieved calculating cp at averaged temperatures at every step. This is because the averaged cp just goes from the first step to the last with one step. But if the cp is found at every step it gives a more accurate representation of how the efficiency will work.

Task 5 (VII)
For this Task, we considered a system with ηcomp = 0.85, ηturb = 0.85, ε_regen = 0.75, P2 = 500 kPa, γ = 0.25, and T4 = 1600 K. Because the incident angle of the sun’s rays varies with time over the peak collection hours of the day, Q_s varies from 2.0 to 3.0 MW over the peak collection period of the daytime. We then used the program developed in Task 4 to find the alpha’s and the variations of of the net power output, the heat input in the combustor, and the efficiency of the cycle with varying Q_s. We ran task_4 for 4 different values of Q_s- 0, 2, 2.5, 3 MW. We found the values for alpha, power output, heat input to the combustor and efficiency of given below.

Q_s (kW) 0 2000 2500 3000
alpha 54.6895 114.0972 155.8688 244.8664
Power (kW) 1.9268e+03 1.9063e+03 1.9013e+03 1.8962e+03
Heat (kW) 4.3481e+03 2.2043e+03 1.6673e+03 1.1310e+03
Efficiency 0.4431 0.4534 0.4562 0.4590
Table 3: Results for different values of alpha, Power, Heat and Efficiency for different values of
Qs
This table is summarized with the graphs given below.

Task 5b (VIII)
In this task we asked the question that if the mole fraction of propane (γ ) in the fuel increases to
0.5 or decreases to 0, how will the system performance change? We assumed that the Qs value is 0 MW for this specific case. We ran task 4 for different values of γ- 0, 0.125, 0.25, 0.375, 0.5 and tabulated the results, given below.

γ 0 0.125 0.25 0.375 0.5
Power(kW) 1.9364e+03 1.9308e+03 1.9268e+03 1.9238e+03 1.9214e+03
n_comp 0.4451 0.4440 0.4431 0.4425 0.4420
alpha_f 39.4117 47.0505 54.6895 62.3285 69.9675
Heat (kW) 4.3500e+03 4.3489e+03 4.3481e+03 4.3475e+03 4.3470e+03
Table 4: Results of Power, Efficiency, Heat with different values of gamma This shows that the efficiency, heat to the combustor and power output doesn’t change much with changing gamma while they do decrease. So if the gamma were to decrease to zero, or only methane was being burnt, the power output is the highest as compared to when 50% of propane and methane were burnt. This is because propane combustion requires more amounts of oxygen and hence there is less moles of total gas burned. The heat input to the combustor, decreases as well but not as much as the power input. This results in the efficiency decreasing.

Conclusion
In this project, we analysed the hybrid solar fossil-fuel gas turbine system, to construct a computer simulation of the performance of the system. For this we first developed a program that found the value of the molar air-to-fuel ratio α that would result in a product temperature Tp for specified reactant temperature Tr. We then wrote a computer program to compute the exit condition from the compressor and the work required per kilogram of air from 1 to 2 using a constant cˆ p and cˆ v analysis. Next, using the chemical reaction for this propane and methane mixture with air (O2 + 3.76 N2) with possible products being CO2, H2O, O2 and N2, we derived relations for the mole fraction of each exhaust gas species as a function of molar air to fuel ratio α. Using these relations we set up a computer program to compute cˆ p,prod for the product gas mixture as the mole-fraction weighted average of the cˆ p values for the individual species.

We then combined all these steps into one to recreate a cycle computationally and analyse its overall performance. We then ran the computation for γ = 0.25, Q_s = 0, P2 = 500 kPa, T4 = 1600 K, ηcomp = 0.85, ηturb = 0.85, eps_regen = 0.75 and found the results given in table 2. We then looked at how changing Qs input and gamma change the efficiency, power input, and heat input. We found that the efficiency at constant cp, is much higher than the achieved calculating cp at averaged temperatures at every step. This is because the averaged cp just goes from the first step to the last with one step. But if the cp is found at every step it gives a more accurate representation of how the efficiency will work. As the solar heat input increases, the Power output, heat input to the combustor increase. But the efficiency and alpha increase.
Appendix I

II
III

function [alpha] = task_1(Tp, Tr, gamma) %the heat of formation from table, [kJ/kmol] h0_c3h8 =-103850; h0_ch4 = -74850; h0_o2 = 0; h0_h2o = -241820; h0_co2 = -393520; h0_n2 = 0;
%specific heat from table, [kJ/kmol*K] T_avg = 0.5 * (Tr+Tp); cp_c3h8 = -4.04 + 30.48*10^(-2)*T_avg – 15.72*10^(-5)*(T_avg)^2 + 31.74*10^(-9)*(T_avg)^3; cp_ch4 = 19.89 + 5.024*10^(-2)*T_avg +1.269*10^(-5)*(T_avg)^2 – 11.01*10^(-9)*(T_avg)^3; cp_o2 = 25.48 + 1.520*10^(-2)*T_avg – 0.7155*10^(-5)*(T_avg)^2 + 1.312*10^(-9)*(T_avg)^3; cp_h2o = 32.24 + 0.1923*10^(-2)*T_avg + 1.055*10^(-5)*(T_avg)^2 – 3.595*10^(-9)*(T_avg)^3; cp_co2 = 22.26 + 5.981*10^(-2)*T_avg – 3.501*10^(-5)*(T_avg)^2 + 7.469*10^(-9)*(T_avg)^3; cp_n2 = 28.9 – 0.1571*10^(-2)*T_avg + 0.8081*10^(-5)*(T_avg)^2 – 2.873*10^(-9)*(T_avg)^3; %enthalpy h_c3h8 =h0_c3h8 + cp_c3h8 *(Tr-298.15); h_ch4 = h0_ch4 + cp_ch4 * (Tr-298.15); h_o2_i = h0_o2 + cp_o2 * (Tr-298.15); h_o2_o = h0_o2 + cp_o2 * (Tp-298.15); h_h2o = h0_h2o + cp_h2o * (Tp-298.15); h_co2 = h0_co2 + cp_co2 * (Tp-298.15); h_n2_i = h0_n2 + cp_n2 * (Tr-298.15); h_n2_o = h0_n2 + cp_n2 * (Tp-298.15); alpha = 4.76*(-gamma*h_c3h8 – (1-gamma)*h_ch4 – (3*gamma+2)*h_o2_o + (2+2*gamma)*h_h2o + (1+2*gamma)*h_co2)/(h_o2_i + 3.76*h_n2_i – 3.76*h_n2_o – h_o2_o); end

function [T2_f, W] =task_2(T1,P1,P2, eta_comp) Appendix IV
%guess T2=T2_g to be average, Tf was taken from task 1
Tf=1900;
T2_g=Tf;
R=8.314;
err=1;
while err >0.2
T_avg = 0.5*(T1+T2_g); cp = 28.11+ 0.1967*10^(-2)*(T_avg) + 0.4802*10^(-5)*(T_avg)^2 – 1.966*10^(-9)*(T_avg)^3; cv = cp-R; k= cp/cv; %first T2 to calculate difference
T2 = T1* (1+((P2/P1)^((k-1)/k) -1)/eta_comp);
%try to get second T2 for difference calculation T_avg_1 = 0.5*(T1+T2);
cp_1 = 28.11+ 0.1967*10^(-2)*( T_avg_1) + 0.4802*10^(-5)*( T_avg_1)^2 – 1.966*10^(-9)*( T_avg_1)^3; cv_1 = cp_1-R; k_1= cp_1/cv_1; T2_1 = T1* (1+((P2/P1)^((k_1-1)/k_1) -1)/eta_comp); err= abs(T2-T2_1); if err <= 0.2 break else
T2_g = T2; end end T2_f = T2_g;
H2= 0 + cp * (T2_f-T1)/(28.97);
%W=H2-H1, and H1=0 W=H2; end
function [y_h20, y_co2, y_n2, y_o2, cp_prod, cp_air] = task_3(Tp, gamma,alpha) Appendix V
%based on the equation given on page1 n_ex = (3.76*alpha/4.76) + (alpha/4.76-3*gamma-2) + (2+2*gamma) + (1+2*gamma); y_h20 = (2+2*gamma) / n_ex; y_co2 = (1+2*gamma)/ n_ex; y_n2 = ((3.76*alpha)/4.76)/ n_ex; y_o2 = ((alpha/4.76)-3*gamma-2) / n_ex;
%cp of each product cp_o2 = 25.48 + 1.520*10^(-2)*Tp – 0.7155*10^(-5)*(Tp)^2 + 1.312*10^(-9)*(Tp)^3; cp_h20 = 32.24 + 0.1923*10^(-2)*Tp + 1.0551*10^(-5)*(Tp)^2 – 3.595*10^(-9)*(Tp)^3; cp_co2 = 22.26 + 5.981*10^(-2)*Tp – 3.501*10^(-5)*(Tp)^2 + 7.469*10^(-9)*(Tp)^3; cp_n2 = 28.9 – 0.1571*10^(-2)*Tp + 0.8081*10^(-5)*(Tp)^2 – 2.873*10^(-9)*(Tp)^3;
%cp of whole product cp_prod = y_h20*cp_h20 + y_co2*cp_co2 + y_n2*cp_n2 + y_o2*cp_o2;
cp_air = 28.11+ 0.1967*10^(-2)*( Tp) + 0.4802*10^(-5)*( Tp)^2 – 1.966*10^(-9)*( Tp)^3; end

function [Power, Heat, T5_f, alpha_f, W_c, W_t, n_comp,Work_Class, Heat_Class, n_class]= task_4 (~)
P1= 101;
T1=298; Appendix VI
P2=500; T4=1600; gamma= .25;
alpha_stoich = 4.76*(2+3*gamma); Qs=0; eff_comp=0.85; eff_turb=0.85; eps_regen=0.75; i = 1; T5_f(1,:) =1000; alpha(1,:) = alpha_stoich; while i<4 m_prod(i,:) = 28.014*3.76*alpha(i,:)/4.76 + 15.9999*2 *(alpha(i,:)/4.76 – 3*gamma -2)+18. 01528*(2+2*gamma)+44.01*(1+2*gamma); m_air = 6; n_dot_air(i,:)= m_air/28.97; n_prod(i,:) = 3.76*alpha(i,:)/4.76 + (alpha(i,:)/4.76 – 3*gamma -2)+(2+2*gamma)+(1+2*gamma); n_dot_prod(i,:)= m_air/(m_prod(i,:)/n_prod(i,:)); [T2_final(i,:),W(i,:),cp_2(i,:)] = task_2(T1,P1,P2,eff_comp,1900);
[T5_f(i+1,:), W_5(i,:),cp_5(i,:)] =task_2a(T4,P1,P2, eff_turb,T5_f(i,:), gamma,alpha(i,:));
T_avg_25= 0.5*(T2_final(i,:)+T5_f(i,:)); [y_h20(1,:), y_co2(i,:), y_n2(i,:), y_o2(i,:), cp_prod(i,:), cp_air(i,:)] = task_3(T_avg_25, gamma,alpha(i,:));
ncp_min = min(cp_prod(i,:)*n_dot_prod(i,:), cp_air(i,:)*n_dot_prod(i,:)); T2_r(i,:) = T2_final(i,:)+ (ncp_min* eps_regen* (T5_f(i,:)- T2_final(i,:)))/(n_dot_air(i,:)
*cp_air(i,:)); [T3, cp_23] = root_finder(Qs,T2_r(i,:),n_dot_air(i,:)); i = i+1; alpha(i,:)= task_1(T4, T3, gamma); end
alpha_f = alpha(3,:); W_t = n_dot_prod(3,:)*cp_5(3,:)*(T4-T5_f(3,:));
W_c = n_dot_air(3,:)*cp_2(3,:)*(T2_final(3,:)-T1);
Power= n_dot_prod(3,:)*cp_5(3,:)*(T4-T5_f(3,:))- n_dot_air(3,:)*cp_2(3,:)*(T2_final(3,:)-T1);
[~, ~, ~, ~, cp_burn, ~] = task_3(.5*(T3+T4), gamma,alpha(3,:)); Heat =(n_dot_prod(3,:)*cp_burn*(T4-T3))+ Qs; k=(28.11/(28.11-8.314)+ (28.3369/(28.3369-8.314)))/2;
Work_Class =(eff_turb*(1-(P1/P2)^((k-1)/k))- (1/eff_comp)*(T1/T4)*((P2/P1)^((k-1)/k)-1));
Heat_Class = (1-eps_regen*(1-eff_turb+eff_turb*((P1/P2)^((k-1)/k)))-(1-eps_regen)*(T1/T4)*(1(1/eff_comp)+(1/eff_comp)*(P2/P1)^((k-1/k)))); n_comp = (n_dot_prod(3,:)*cp_5(3,:)*(T4-T5_f(3,:))- n_dot_air(3,:)*cp_2(3,:)*(T2_final(3,:)-T1))/ (n_dot_prod(3,:)*cp_burn*(T4-T3)); n_class = (eff_turb*(1-(P1/P2)^((k-1)/k))- (1/eff_comp)*(T1/T4)*((P2/P1)^((k-1)/k)-1))/(1 eps_regen*(1-eff_turb+eff_turb*((P1/P2)^((k-1)/k)))-(1-eps_regen)*(T1/T4)*(1-(1/eff_comp)+ (1/eff_comp)*(P2/P1)^((k-1/k)))); end function [T2_f, W,cp] =task_2(T1,P1,P2, eta_comp, T_guess)
%guess T2=T2_g to be average, Tf was taken from task 1
T2_g=T_guess; R=8.314;
err=1; while err >0.2
T_avg = 0.5*(T1+T2_g); cp = 28.11+ 0.1967*10^(-2)*(T_avg) + 0.4802*10^(-5)*(T_avg)^2 – 1.966*10^(-9)*(T_avg)^3; cv = cp-R; k= cp/cv; %first T2 to calculate difference
T2 = T1* (1+((P2/P1)^((k-1)/k) -1)/eta_comp);
%try to get second T2 for difference calculation T_avg_1 = 0.5*(T1+T2);
cp_1 = 28.11+ 0.1967*10^(-2)*( T_avg_1) + 0.4802*10^(-5)*( T_avg_1)^2 – 1.966*10^(-9)*( T_avg_1)^3; cv_1 = cp_1-R; k_1= cp_1/cv_1; T2_1 = T1* (1+((P2/P1)^((k_1-1)/k_1) -1)/eta_comp); err= abs(T2-T2_1); if err <= 0.2 break else
T2_g = T2; end end T2_f = T2_g;
H2= 0 + cp * (T2_f-T1)/(28.97);
%W=H2-H1, and H1=0 W=H2;
end function [T2_f, W,cp] =task_2a(T1,P1,P2, eta_comp, T_guess, gamma,alpha) %guess T2=T2_g to be average, Tf was taken from task 1
T2_g=T_guess; R=8.314;
err=1; while err >0.2 T_avg = 0.5*(T1+T2_g); [~, ~, ~, ~, cp_prod, ~] = task_3(T_avg, gamma,alpha); cp = cp_prod; cv = cp-R; k= cp/cv; %first T2 to calculate difference T2 = T1* (1+((P1/P2)^((k-1)/k) -1)*eta_comp); err= abs(T2-T2_g); if err <= 0.2
break else
T2_g = T2; end end T2_f = T2_g;
H2= 0 + cp * (T1-T2_f)/(28.97);
%W=H2-H1, and H1=0
W=H2;
end
function [y_h20, y_co2, y_n2, y_o2, cp_prod, cp_air] = task_3(Tp, gamma,alpha)
%based on the equation given on page1 n_ex = (3.76*alpha/4.76) + (alpha/4.76-3*gamma-2) + (2+2*gamma) + (1+2*gamma); y_h20 = (2+2*gamma) / n_ex; y_co2 = (1+2*gamma)/ n_ex; y_n2 = (3.76*alpha/4.76)/ n_ex; y_o2 = (alpha/4.76-3*gamma-2) / n_ex;
%cp of each product cp_o2 = 25.48 + 1.520*10^(-2)*Tp – 0.7155*10^(-5)*(Tp)^2 + 1.312*10^(-9)*(Tp)^3; cp_h20 = 32.24 + 0.1923*10^(-2)*Tp + 1.0551*10^(-5)*(Tp)^2 – 3.595*10^(-9)*(Tp)^3; cp_co2 = 22.26 + 5.981*10^(-2)*Tp – 3.501*10^(-5)*(Tp)^2 + 7.469*10^(-9)*(Tp)^3; cp_n2 = 28.9 – 0.1571*10^(-2)*Tp + 0.8081*10^(-5)*(Tp)^2 – 2.873*10^(-9)*(Tp)^3;
%cp of whole product cp_prod = y_h20*cp_h20 + y_co2*cp_co2 + y_n2*cp_n2 + y_o2*cp_o2; cp_air = 28.11+ 0.1967*10^(-2)*( Tp) + 0.4802*10^(-5)*( Tp)^2 – 1.966*10^(-9)*( Tp)^3; end function [alpha] = task_1(Tp, Tr, gamma)
%the heat of formation from table, [kJ/kmol] h0_c3h8 =-103850; h0_ch4 = -74850; h0_o2 = 0; h0_h2o = -241820; h0_co2 = -393520; h0_n2 = 0;
%specific heat from table, [kJ/kmol*K] T_avg = 0.5 * (Tr+Tp); cp_c3h8 = -4.04 + 30.48*10^(-2)*T_avg – 15.72*10^(-5)*(T_avg)^2 + 31.74*10^(-9)*(T_avg)^3; cp_ch4 = 19.89 + 5.024*10^(-2)*T_avg +1.269*10^(-5)*(T_avg)^2 – 11.01*10^(-9)*(T_avg)^3; cp_o2 = 25.48 + 1.520*10^(-2)*T_avg – 0.7155*10^(-5)*(T_avg)^2 + 1.312*10^(-9)*(T_avg)^3; cp_h2o = 32.24 + 0.1923*10^(-2)*T_avg + 1.055*10^(-5)*(T_avg)^2 – 3.595*10^(-9)*(T_avg)^3; cp_co2 = 22.26 + 5.981*10^(-2)*T_avg – 3.501*10^(-5)*(T_avg)^2 + 7.469*10^(-9)*(T_avg)^3; cp_n2 = 28.9 – 0.1571*10^(-2)*T_avg + 0.8081*10^(-5)*(T_avg)^2 – 2.873*10^(-9)*(T_avg)^3;
%enthalpy h_c3h8 =h0_c3h8 + cp_c3h8 *(Tr-298.15); h_ch4 = h0_ch4 + cp_ch4 * (Tr-298.15); h_o2_i = h0_o2 + cp_o2 * (Tr-298.15); h_o2_o = h0_o2 + cp_o2 * (Tp-298.15); h_h2o = h0_h2o + cp_h2o * (Tp-298.15); h_co2 = h0_co2 + cp_co2 * (Tp-298.15); h_n2_i = h0_n2 + cp_n2 * (Tr-298.15); h_n2_o = h0_n2 + cp_n2 * (Tp-298.15);
alpha = 4.76*(-gamma*h_c3h8 – (1-gamma)*h_ch4 – (3*gamma+2)*h_o2_o + (2+2*gamma)*h_h2o +
(1+2*gamma)*h_co2)/(h_o2_i + 3.76*h_n2_i – 3.76*h_n2_o – h_o2_o);
end function [T3_f, cp23] = root_finder(Qs, T2_r, n_dot_air) syms cp_23 T3 eqn = [cp_23 == 28.11+ 0.1967*10^(-2)*(0.5*(T3+T2_r)) + 0.4802*10^(-5)*(0.5*(T3+T2_r))^2 – 1.966 *10^(-9)*(0.5*(T3+T2_r))^3,T3 == Qs/(n_dot_air*cp_23)+T2_r]; v = solve(eqn, [cp_23, T3]); cp_23= vpa(v.cp_23); T3 = vpa(v.T3);
p = find(T3 >800& T3 < 1500); T3_f = double(T3(p)); cp23 = double(cp_23(p)); end

function [Power, Heat, T5_f, alpha, W_c, W_t, n_comp,T3, T2_r,cp_23]= task_5(~) Qs= [0,500,1000,1500,2000,2500,3000]; gamma= 0.25;
P1= 101;
Appendix VII
T1=298;
P2=500; T4=1600; alpha_stoich = 4.76*(2+3*gamma); eff_comp=0.85; eff_turb=0.85; eps_regen=0.75; for l = 1:length(Qs) i = 1; T5_f(1,l) =1000; alpha(1,l) = alpha_stoich; while i<4 m_prod(i,l) = 28.014*3.76*alpha(i,l)/4.76 + 15.9999*2 *(alpha(i,l)/4.76 – 3*gamma -2)+18. 01528*(2+2*gamma)+44.01*(1+2*gamma); m_air = 6; n_dot_air(i,l)= m_air/28.97; n_prod(i,l) = 3.76*alpha(i,l)/4.76 + (alpha(i,l)/4.76 – 3*gamma -2)+(2+2*gamma)+(1+2*gamma); n_dot_prod(i,l)= m_air/(m_prod(i,l)/n_prod(i,l)); [T2_final(i,l),W(i,l),cp_2(i,l)] = task_2(T1,P1,P2,eff_comp,1900);
[T5_f(i+1,l), W_5(i,l),cp_5(i,l)] =task_2a(T4,P1,P2, eff_turb,T5_f(i,l), gamma,alpha(i,l));
T_avg_25= 0.5*(T2_final(i,l)+T5_f(i,l)); [y_h20(i,l), y_co2(i,l), y_n2(i,l), y_o2(i,l), cp_prod(i,l), cp_air(i,l)] = task_3(T_avg_25, gamma,alpha(i,l));
ncp_min = min(cp_prod(i,l)*n_dot_prod(i,l), cp_air(i,l)*n_dot_prod(i,l)); T2_r(i,l) = T2_final(i,l)+ (ncp_min* eps_regen* (T5_f(i,l)- T2_final(i,l)))/(n_dot_air(i,l)
*cp_air(i,l));
%disp(i);
%disp(n_dot_air(i,l)); %[T3(i,l), cp_23(i,l)] = root_finder(Qs(l),T2_r(i,l),n_dot_air(i,l)); cp_23(i,l) = 28.11+ 0.1967*10^(-2)*((T2_r(i,l))) + 0.4802*10^(-5)*((T2_r(i,l)))^2 – 1.966*10^ (-9)*((T2_r(i,l)))^3; T3(i,l) = Qs(l)/(n_dot_air(i,l)*cp_23(i,l))+T2_r(i,l); i = i+1; alpha(i,l)= task_1(T4, T3(i-1,l), gamma); end
alpha_f(1,l) = alpha(3,l); W_t(1,l) = n_dot_prod(3,l)*cp_5(3,l)*(T4-T5_f(3,l));
W_c(1,l) = n_dot_air(3,l)*cp_2(3,l)*(T2_final(3,l)-T1);

Power(1,l)= n_dot_prod(3,l)*cp_5(3,l)*(T4-T5_f(3,l))- n_dot_air(3,l)*cp_2(3,l)*(T2_final(3,l)-
T1);
[~, ~, ~, ~, cp_burn, ~] = task_3(.5*(T3(i-1,l)+T4), gamma,alpha(3,l)); Heat(1,l) =(n_dot_prod(3,l)*cp_burn*(T4-T3(i-1,l)));
n_comp(1,l) = (n_dot_prod(3,l)*cp_5(3,l)*(T4-T5_f(3,l))- n_dot_air(3,l)*cp_2(3,l)*(T2_final(3,l)T1))/(Heat(1,l)+Qs(l)); end subplot(2,2,1) plot(Qs, Power ) xlabel(‘Q_{solar} (KW)’) ylabel(‘Power Output(KW)’) subplot(2,2,3) plot(Qs, Heat) xlabel(‘Q_{solar} (KW)’)
ylabel(‘Heat Input in combustor (KW)’) subplot(2,2,2) plot(Qs, n_comp) xlabel(‘Q_{solar} (KW)’) ylabel(‘Overall Efficiency’) subplot(2,2,4) plot(Qs, alpha_f) xlabel(‘Q_{solar} (KW)’) ylabel(‘{alpha}’) end function [T2_f, W,cp] =task_2(T1,P1,P2, eta_comp, T_guess)
%guess T2=T2_g to be average, Tf was taken from task 1
T2_g=T_guess; R=8.314;
err=1; while err >0.2
T_avg = 0.5*(T1+T2_g); cp = 28.11+ 0.1967*10^(-2)*(T_avg) + 0.4802*10^(-5)*(T_avg)^2 – 1.966*10^(-9)*(T_avg)^3; cv = cp-R; k= cp/cv; %first T2 to calculate difference
T2 = T1* (1+((P2/P1)^((k-1)/k) -1)/eta_comp);
%try to get second T2 for difference calculation T_avg_1 = 0.5*(T1+T2);
cp_1 = 28.11+ 0.1967*10^(-2)*( T_avg_1) + 0.4802*10^(-5)*( T_avg_1)^2 – 1.966*10^(-9)*( T_avg_1)^3; cv_1 = cp_1-R; k_1= cp_1/cv_1; T2_1 = T1* (1+((P2/P1)^((k_1-1)/k_1) -1)/eta_comp); err= abs(T2-T2_1); if err <= 0.2 break else
T2_g = T2; end end T2_f = T2_g;
H2= 0 + cp * (T2_f-T1)/(28.97);
%W=H2-H1, and H1=0 W=H2;
end function [T2_f, W,cp] =task_2a(T1,P1,P2, eta_comp, T_guess, gamma,alpha) %guess T2=T2_g to be average, Tf was taken from task 1
T2_g=T_guess; R=8.314;
err=1; while err >0.2 T_avg = 0.5*(T1+T2_g); [~, ~, ~, ~, cp_prod, ~] = task_3(T_avg, gamma,alpha); cp = cp_prod;
cv = cp-R; k= cp/cv; %first T2 to calculate difference T2 = T1* (1+((P1/P2)^((k-1)/k) -1)*eta_comp); err= abs(T2-T2_g); if err <= 0.2
break else T2_g = T2; end end T2_f = T2_g;
H2= 0 + cp * (T1-T2_f)/(28.97);
%W=H2-H1, and H1=0
W=H2; end function [y_h20, y_co2, y_n2, y_o2, cp_prod, cp_air] = task_3(Tp, gamma,alpha) %based on the equation given on page1 n_ex = (3.76*alpha/4.76) + (alpha/4.76-3*gamma-2) + (2+2*gamma) + (1+2*gamma); y_h20 = (2+2*gamma) / n_ex; y_co2 = (1+2*gamma)/ n_ex; y_n2 = (3.76*alpha/4.76)/ n_ex; y_o2 = (alpha/4.76-3*gamma-2) / n_ex; %cp of each product cp_o2 = 25.48 + 1.520*10^(-2)*Tp – 0.7155*10^(-5)*(Tp)^2 + 1.312*10^(-9)*(Tp)^3; cp_h20 = 32.24 + 0.1923*10^(-2)*Tp + 1.0551*10^(-5)*(Tp)^2 – 3.595*10^(-9)*(Tp)^3; cp_co2 = 22.26 + 5.981*10^(-2)*Tp – 3.501*10^(-5)*(Tp)^2 + 7.469*10^(-9)*(Tp)^3; cp_n2 = 28.9 – 0.1571*10^(-2)*Tp + 0.8081*10^(-5)*(Tp)^2 – 2.873*10^(-9)*(Tp)^3; %cp of whole product cp_prod = y_h20*cp_h20 + y_co2*cp_co2 + y_n2*cp_n2 + y_o2*cp_o2; cp_air = 28.11+ 0.1967*10^(-2)*( Tp) + 0.4802*10^(-5)*( Tp)^2 – 1.966*10^(-9)*( Tp)^3; end function [alpha] = task_1(Tp, Tr, gamma)
%the heat of formation from table, [kJ/kmol] h0_c3h8 =-103850; h0_ch4 = -74850; h0_o2 = 0; h0_h2o = -241820; h0_co2 = -393520; h0_n2 = 0;
%specific heat from table, [kJ/kmol*K] T_avg = 0.5 * (Tr+Tp); cp_c3h8 = -4.04 + 30.48*10^(-2)*T_avg – 15.72*10^(-5)*(T_avg)^2 + 31.74*10^(-9)*(T_avg)^3; cp_ch4 = 19.89 + 5.024*10^(-2)*T_avg +1.269*10^(-5)*(T_avg)^2 – 11.01*10^(-9)*(T_avg)^3; cp_o2 = 25.48 + 1.520*10^(-2)*T_avg – 0.7155*10^(-5)*(T_avg)^2 + 1.312*10^(-9)*(T_avg)^3; cp_h2o = 32.24 + 0.1923*10^(-2)*T_avg + 1.055*10^(-5)*(T_avg)^2 – 3.595*10^(-9)*(T_avg)^3; cp_co2 = 22.26 + 5.981*10^(-2)*T_avg – 3.501*10^(-5)*(T_avg)^2 + 7.469*10^(-9)*(T_avg)^3; cp_n2 = 28.9 – 0.1571*10^(-2)*T_avg + 0.8081*10^(-5)*(T_avg)^2 – 2.873*10^(-9)*(T_avg)^3;
%enthalpy h_c3h8 =h0_c3h8 + cp_c3h8 *(Tr-298.15); h_ch4 = h0_ch4 + cp_ch4 * (Tr-298.15); h_o2_i = h0_o2 + cp_o2 * (Tr-298.15); h_o2_o = h0_o2 + cp_o2 * (Tp-298.15); h_h2o = h0_h2o + cp_h2o * (Tp-298.15); h_co2 = h0_co2 + cp_co2 * (Tp-298.15); h_n2_i = h0_n2 + cp_n2 * (Tr-298.15); h_n2_o = h0_n2 + cp_n2 * (Tp-298.15);
alpha = 4.76*(-gamma*h_c3h8 – (1-gamma)*h_ch4 – (3*gamma+2)*h_o2_o + (2+2*gamma)*h_h2o + (1+2*gamma)*h_co2)/(h_o2_i + 3.76*h_n2_i – 3.76*h_n2_o – h_o2_o); end

function [Power, Heat, T5_f, alpha, W_c, W_t, n_comp,T3, T2_r,cp_23]= task_5b(~)
Qs= 0;
Appendix VIII gamma= [0,0.125,0.25,0.375,0.5]; P1= 101;
T1=298;
P2=500; T4=1600;
eff_comp=0.85; eff_turb=0.85; eps_regen=0.75; for l = 1:length(gamma) i = 1; T5_f(1,l) =1000; alpha_stoich = 4.76*(2+3*gamma(l)); alpha(i,l) = alpha_stoich; while i<4 m_prod(i,l) = 28.014*3.76*alpha(i,l)/4.76 + 15.9999*2 *(alpha(i,l)/4.76 – 3*gamma(l) -2)+18. 01528*(2+2*gamma(l))+44.01*(1+2*gamma(l)); m_air = 6; n_dot_air(i,l)= m_air/28.97; n_prod(i,l) = 3.76*alpha(i,l)/4.76 + (alpha(i,l)/4.76 – 3*gamma(l) -2)+(2+2*gamma(l))+ (1+2*gamma(l)); n_dot_prod(i,l)= m_air/(m_prod(i,l)/n_prod(i,l)); [T2_final(i,l),W(i,l),cp_2(i,l)] = task_2(T1,P1,P2,eff_comp,1900); [T5_f(i+1,l), W_5(i,l),cp_5(i,l)] =task_2a(T4,P1,P2, eff_turb,T5_f(i,l), gamma(l),alpha(i, l)); T_avg_25= 0.5*(T2_final(i,l)+T5_f(i,l)); [y_h20(i,l), y_co2(i,l), y_n2(i,l), y_o2(i,l), cp_prod(i,l), cp_air(i,l)] = task_3(T_avg_25, gamma(l),alpha(i,l));
ncp_min = min(cp_prod(i,l)*n_dot_prod(i,l), cp_air(i,l)*n_dot_prod(i,l)); T2_r(i,l) = T2_final(i,l)+ (ncp_min* eps_regen* (T5_f(i,l)- T2_final(i,l)))/(n_dot_air(i,l)
*cp_air(i,l));
%disp(i);
%disp(n_dot_air(i,l));
[T3(i,l), cp_23(i,l)] = root_finder(Qs,T2_r(i,l),n_dot_air(i,l));
%cp_23(i,l) = 28.11+ 0.1967*10^(-2)*((T2_r(i,l))) + 0.4802*10^(-5)*((T2_r(i,l)))^2 – 1.966
*10^(-9)*((T2_r(i,l)))^3; %T3(i,l) = Qs/(n_dot_air(i,l)*cp_23(i,l))+T2_r(i,l); i = i+1; alpha(i,l)= task_1(T4, T3(i-1,l), gamma(l)); end
alpha_f(1,l) = alpha(3,l); W_t(1,l) = n_dot_prod(3,l)*cp_5(3,l)*(T4-T5_f(3,l));
W_c(1,l) = n_dot_air(3,l)*cp_2(3,l)*(T2_final(3,l)-T1);

Power(1,l)= n_dot_prod(3,l)*cp_5(3,l)*(T4-T5_f(3,l))- n_dot_air(3,l)*cp_2(3,l)*(T2_final(3,l)-
T1);
[~, ~, ~, ~, cp_burn, ~] = task_3(.5*(T3(i-1,l)+T4), gamma(l),alpha(3,l)); Heat(1,l) =(n_dot_prod(3,l)*cp_burn*(T4-T3(i-1,l)));
n_comp(1,l) = (n_dot_prod(3,l)*cp_5(3,l)*(T4-T5_f(3,l))- n_dot_air(3,l)*cp_2(3,l)*(T2_final(3,l)T1))/(Heat(1,l)); end end function [T2_f, W,cp] =task_2(T1,P1,P2, eta_comp, T_guess)
%guess T2=T2_g to be average, Tf was taken from task 1
T2_g=T_guess; R=8.314;
err=1; while err >0.2
T_avg = 0.5*(T1+T2_g); cp = 28.11+ 0.1967*10^(-2)*(T_avg) + 0.4802*10^(-5)*(T_avg)^2 – 1.966*10^(-9)*(T_avg)^3; cv = cp-R; k= cp/cv; %first T2 to calculate difference
T2 = T1* (1+((P2/P1)^((k-1)/k) -1)/eta_comp);
%try to get second T2 for difference calculation T_avg_1 = 0.5*(T1+T2);
cp_1 = 28.11+ 0.1967*10^(-2)*( T_avg_1) + 0.4802*10^(-5)*( T_avg_1)^2 – 1.966*10^(-9)*( T_avg_1)^3; cv_1 = cp_1-R; k_1= cp_1/cv_1; T2_1 = T1* (1+((P2/P1)^((k_1-1)/k_1) -1)/eta_comp); err= abs(T2-T2_1); if err <= 0.2 break else
T2_g = T2; end end T2_f = T2_g;
H2= 0 + cp * (T2_f-T1)/(28.97);
%W=H2-H1, and H1=0 W=H2;
end function [T2_f, W,cp] =task_2a(T1,P1,P2, eta_comp, T_guess, gamma,alpha) %guess T2=T2_g to be average, Tf was taken from task 1
T2_g=T_guess; R=8.314;
err=1; while err >0.2 T_avg = 0.5*(T1+T2_g); [~, ~, ~, ~, cp_prod, ~] = task_3(T_avg, gamma,alpha); cp = cp_prod; cv = cp-R; k= cp/cv; %first T2 to calculate difference T2 = T1* (1+((P1/P2)^((k-1)/k) -1)*eta_comp); err= abs(T2-T2_g); if err <= 0.2
break else
T2_g = T2; end end T2_f = T2_g;
H2= 0 + cp * (T1-T2_f)/(28.97);
%W=H2-H1, and H1=0
W=H2; end
function [y_h20, y_co2, y_n2, y_o2, cp_prod, cp_air] = task_3(Tp, gamma,alpha)
%based on the equation given on page1 n_ex = (3.76*alpha/4.76) + (alpha/4.76-3*gamma-2) + (2+2*gamma) + (1+2*gamma); y_h20 = (2+2*gamma) / n_ex; y_co2 = (1+2*gamma)/ n_ex; y_n2 = (3.76*alpha/4.76)/ n_ex; y_o2 = (alpha/4.76-3*gamma-2) / n_ex;
%cp of each product cp_o2 = 25.48 + 1.520*10^(-2)*Tp – 0.7155*10^(-5)*(Tp)^2 + 1.312*10^(-9)*(Tp)^3; cp_h20 = 32.24 + 0.1923*10^(-2)*Tp + 1.0551*10^(-5)*(Tp)^2 – 3.595*10^(-9)*(Tp)^3; cp_co2 = 22.26 + 5.981*10^(-2)*Tp – 3.501*10^(-5)*(Tp)^2 + 7.469*10^(-9)*(Tp)^3; cp_n2 = 28.9 – 0.1571*10^(-2)*Tp + 0.8081*10^(-5)*(Tp)^2 – 2.873*10^(-9)*(Tp)^3;
%cp of whole product cp_prod = y_h20*cp_h20 + y_co2*cp_co2 + y_n2*cp_n2 + y_o2*cp_o2; cp_air = 28.11+ 0.1967*10^(-2)*( Tp) + 0.4802*10^(-5)*( Tp)^2 – 1.966*10^(-9)*( Tp)^3; end function [alpha] = task_1(Tp, Tr, gamma)
%the heat of formation from table, [kJ/kmol] h0_c3h8 =-103850; h0_ch4 = -74850; h0_o2 = 0; h0_h2o = -241820; h0_co2 = -393520; h0_n2 = 0;
%specific heat from table, [kJ/kmol*K] T_avg = 0.5 * (Tr+Tp); cp_c3h8 = -4.04 + 30.48*10^(-2)*T_avg – 15.72*10^(-5)*(T_avg)^2 + 31.74*10^(-9)*(T_avg)^3; cp_ch4 = 19.89 + 5.024*10^(-2)*T_avg +1.269*10^(-5)*(T_avg)^2 – 11.01*10^(-9)*(T_avg)^3; cp_o2 = 25.48 + 1.520*10^(-2)*T_avg – 0.7155*10^(-5)*(T_avg)^2 + 1.312*10^(-9)*(T_avg)^3; cp_h2o = 32.24 + 0.1923*10^(-2)*T_avg + 1.055*10^(-5)*(T_avg)^2 – 3.595*10^(-9)*(T_avg)^3; cp_co2 = 22.26 + 5.981*10^(-2)*T_avg – 3.501*10^(-5)*(T_avg)^2 + 7.469*10^(-9)*(T_avg)^3; cp_n2 = 28.9 – 0.1571*10^(-2)*T_avg + 0.8081*10^(-5)*(T_avg)^2 – 2.873*10^(-9)*(T_avg)^3;
%enthalpy h_c3h8 =h0_c3h8 + cp_c3h8 *(Tr-298.15); h_ch4 = h0_ch4 + cp_ch4 * (Tr-298.15); h_o2_i = h0_o2 + cp_o2 * (Tr-298.15); h_o2_o = h0_o2 + cp_o2 * (Tp-298.15); h_h2o = h0_h2o + cp_h2o * (Tp-298.15); h_co2 = h0_co2 + cp_co2 * (Tp-298.15); h_n2_i = h0_n2 + cp_n2 * (Tr-298.15); h_n2_o = h0_n2 + cp_n2 * (Tp-298.15); alpha = 4.76*(-gamma*h_c3h8 – (1-gamma)*h_ch4 – (3*gamma+2)*h_o2_o + (2+2*gamma)*h_h2o +

(1+2*gamma)*h_co2)/(h_o2_i + 3.76*h_n2_i – 3.76*h_n2_o – h_o2_o); end
11/23/18 6:35 PM C:UsersprincDocumentsFall … oot_finder.m 1 of 1
function [T3_f, cp23] = root_finder(Qs, T2_r, n_dot_air) syms cp_23 T3 eqn = [cp_23 == 28.11+ 0.1967*10^(-2)*(0.5*(T3+T2_r)) + 0.4802*10^(-5)*(0.5*(T3+T2_r))^2 – 1.966*10^(-9)*(0.5*(T3+T2_r))^3,T3 == Qs/(n_dot_air*cp_23)+T2_r]; v = solve(eqn, [cp_23, T3]); cp_23= vpa(v.cp_23);
Appendix IX
T3 = vpa(v.T3);
p = find(T3 >800& T3 < 1500); T3_f = double(T3(p)); cp23 = double(cp_23(p)); end