API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Mar 21st, 2019
57
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 16.86 KB | None | 0 0
raw download clone embed print report
  1. close all; clear all; clc
  2. %% To make sure that matlab will find the functions. You must change it to your situation
  3. abspath_to_generalfolder='My folder; % absolute reference to General folder
  4. addpath(abspath_to_generalfolder);
  5. %% Load Nasadatabase
  6. TdataBase=fullfile('General','Nasa','NasaThermalDatabase');
  7. load(TdataBase);
  8. fuel = 'Gasoline';
  9. addpath('General\Nasa')
  10. iSp = myfind({Sp.Name},{fuel,'C2H5OH','O2','CO2','H2O','N2'}); % Find indexes of these species
  11. SpS=Sp(iSp); % Subselection of the database in the order according to {'Gasoline','O2','CO2','H2O','N2'}
  12. NSp = length(SpS);
  13. Mi = [SpS.Mass];
  14. %% Values should not be changed. Ready for use
  15. global Runiv Pref
  16. Runiv=8.314472;
  17. [T0, a0, P0, rho0] = atmosisa(0); % Official ambient pressure at ground level
  18. %% Start composition
  19. V_gas = 1; % volume fraction gasoline
  20. V_eth = 1 - V_gas; % volume fraction ethanol
  21.  
  22. m_fuel = 7.22e-06; % Mass of fuel for one cycle
  23.  
  24. N_gas = V_gas*m_fuel/Mi(1); % moles fuel
  25. N_eth = V_eth*m_fuel/Mi(2); % moles ethanol
  26. % C2H5OH + 3(O2+ 3.76N2) = 2CO2 + 3H2O + 11.28N2
  27. Neth1=[0,N_eth, 3*N_eth, 0, 0, (0.79*3*N_eth)/0.21];
  28. % C8H18+ 12.5(O2+ 3.76N2) = 8CO2+ 9H2O+ 47N2 (SSA Victor)
  29. Ngas1=[N_gas,0, 12.5*N_gas, 0, 0, 12.5*3.67*N_gas];
  30. N1tot=Neth1+Ngas1; % Moles at start
  31.  
  32. m_gas(1:659) = N_gas*Mi(1); % Initial mass of gasoline
  33. m_eth(1:659) = N_eth*Mi(2); % Initial mass of ethanol
  34. m_O2(1:659) = (Neth1(3)+Ngas1(3))/Mi(3); % Initial mass of oxygen
  35.  
  36. Mtot1 = N1tot*Mi'; % Total Mass
  37. X1_tot = N1tot/sum(N1tot); % X (mole) fraction
  38. Y1_tot = N1tot.*Mi/Mtot1; % Y (mass) fraction at start
  39. %% After combustion
  40. % C2H5OH + 3(O2+ 3.76N2) = 2CO2 + 3H2O + 11.28N2
  41. Neth2=[0,0,0, 2*N_eth,3*N_eth,(0.79*3*N_eth)/0.21];
  42. % C8H18+ 12.5(O2+ 3.76N2) = 8CO2+ 9H2O+ 47N2 (SSA Victor)
  43. Ngas2=[0,0,0,8*N_gas,9*N_gas, 12.5*3.67*N_gas];
  44. N2tot=Neth2+Ngas2; % molefractions of the mixture after combustion
  45.  
  46. Mtot2 = N2tot*Mi'; % Mass of components after combustion of 1 mole fuel and 1 mole ethanol
  47. X2_tot = N2tot/sum(N2tot); % X (mole) fraction
  48. Y2_tot = N2tot.*Mi/Mtot2; % Massfractions of the mixture after combustion
  49. %% Begin conditions
  50. Mair = 28.96e-3; % Molair mass air
  51. Rg = Runiv/Mair; % specific gas constant of air
  52. p(1)=P0;T(1)=T0; % Initial pressure and temperature values
  53. Ca(1)=0.0;V(1)=Vcyl(Ca(1)); % Initial crank angle and volume
  54. m1(1)= p(1)*V(1)/(Rg*T(1)); % Initial mass
  55. k = 1.4; % Simple k for the Poisson relations (from air)(only used for motorized mass)
  56. Cv = 718; % Cv air
  57. dCa=0.5; % Stepsize
  58. Tw = 647; % Temperature of the wall(?)
  59. V_d = 2.2281e-04-2.9300e-05; % Displaced volume (maximum volume change)
  60. r_c = 2.2281e-04/2.9300e-05; % Compression ratio
  61. B = 0.0678; % Bore
  62. S = 0.0536; % Stroke
  63. S_p = 2*S*3000/60; % Average piston speed
  64. c = 0.0064; % Height of cylinder head
  65. E = sqrt(1-((c^2)/((0.5*B)^2))); % Eccencitry
  66. A = pi*B*S + pi * (0.5*B)^2 + 2*pi *(0.5*B)^2 + pi * (c^2)/E * log((1+E)/(1-E)); %Area
  67. dt = dCa/(3000*360/60); % Time difference for every stepsize (used for temperature change)
  68. Qlhv_eth = 26.952*10^6; % LHV of ethanol [J/kg]
  69. Qlhv_gas = 43.448*10^6; % LHV of gasoline [J/kg]
  70. Comb_rate_gas = m_gas(1)/55; % Combustion rate of gasoline [kg/deg]
  71. Comb_rate_eth = m_eth(1)/55; % Combustion rate of ethanol [kg/deg]
  72. %% Fuel intake
  73.  
  74. dm = (Mtot1-m1(1))/(2.2281e-04-(V(1))); % Mass change per unit of volume
  75.  
  76. for i=2:360
  77. Ca(i) = Ca(i-1)+dCa; % Compute new angle
  78. V(i) = Vcyl(Ca(i)); % Cylinder volume for current crank-angle
  79. p(i) = p(i-1); % Pressure is constant
  80. m1(i) = m1(i-1) + dm*(V(i)-V(i-1)); % Compute mass for new volume
  81. T(i) = T(i-1); % Temperature is constant
  82. end
  83.  
  84. %% Compression
  85.  
  86. M_av1 = sum(X1_tot.*Mi); % Molar mass of gas mixture
  87.  
  88. for i = 361:720
  89. if i>= 660; % Initial point of combustion
  90. C1 = 2.28; % Constant for velocity w
  91. C2 = 3.24e-03; % Constant for velocity w
  92. Ca(i)=Ca(i-1)+dCa; % Compute new angle
  93. V(i)=Vcyl(Ca(i)); % Compute new volume
  94. Qc(i) = Qlhv_gas*Comb_rate_gas*dCa + Qlhv_eth*Comb_rate_eth*dCa; % Combustion heat
  95. pm = (((r_c*V_d)/(r_c - 1))^k * P0)/(V(660)^k); % Motorized pressure
  96. w = C1*S_p + C2*((V_d*(T(660-1)))/(p(660-1)*V(660)))*(p(660-1)-pm); % Average cylinder gas velocity
  97. Hc(i) = 3.26*B^-0.2*p(i-1)^0.8*(T(i-1))^-0.55*w^0.8; % Compute heat transfer coefficient for current temperature and pressure
  98. dQ(i) = Hc(i)*A*((T(i-1))-Tw)*dt; % Heat loss
  99. m_added_CO2(i) = ((Comb_rate_gas*dCa*Mi(1)*8) + (Comb_rate_eth*dCa*Mi(2)*2))/Mi(4);
  100. % Change in mass of CO2 by combustion
  101. m_CO2(i) = sum(m_added_CO2); % Total mass of CO2
  102. m_added_H2O(i) = ((Comb_rate_gas*dCa*Mi(1)*9) + (Comb_rate_eth*dCa*Mi(2)*3))/Mi(5);
  103. % Change in mass of H2O by combustion
  104. m_H2O(i) = sum(m_added_H2O); % Total mass of H2O
  105. m_added_O2(i) = -((Comb_rate_gas*dCa*Mi(1)*12.5) + (Comb_rate_eth*dCa*Mi(2)*3))/Mi(3);
  106. % Change in mass of O2 by combustion
  107. m_O2(i) = m_O2(i-1) + sum(m_added_O2(i)); % Total mass of O2
  108. m_gas(i) = m_gas(i-1) - Comb_rate_gas*dCa; % Total mass of gasoline
  109. m_eth(i) = m_eth(i-1) - Comb_rate_eth*dCa; % Total mass of ethanol
  110. m_N2(i) = Y1_tot(6)*m1(1); % Total mass of N2
  111. m1(i) = m1(i-1)-Comb_rate_gas*dCa - Comb_rate_eth*dCa + m_added_CO2(i) + m_added_H2O(i) + m_added_O2(i);
  112. % Compute new total mass
  113. Cv(i) = (m_gas(i)/m1(i))*CvNasa((T(i-1)),SpS(1)) + (m_eth(i)/m1(i))*CvNasa((T(i-1)),SpS(2)) + (m_O2(i)/m1(i))*CvNasa((T(i-1)),SpS(3)) + (m_CO2(i)/m1(i))*CvNasa((T(i-1)),SpS(4)) + (m_H2O(i)/m1(i))*CvNasa((T(i-1)),SpS(5)) + (m_N2(i)/m1(i))*CvNasa((T(i-1)),SpS(6));
  114. % Specific heat capacity for current temperature
  115. T(i) = (T(i-1)) + (-dQ(i) + Qc(i) - p(i-1)*(V(i)-V(i-1)))/(Cv(i)*m1(i));% Compute new temperature using 1st law of thermodynamics
  116. p(i) = m1(i)*(Runiv/M_av1)*T(i)/V(i); % Compute new pressure using ideal gas law
  117. else % Loop before combustion
  118. C1 = 2.28; % Constant for velocity w
  119. C2 = 0; % Constant for velocity w
  120. Ca(i)=Ca(i-1)+dCa; % Compute new angle
  121. V(i)=Vcyl(Ca(i)); % Compute new volume
  122. pm = (((r_c*V_d)/(r_c - 1))^k * P0)/(V(360)^k); % Motorized pressure
  123. w = C1*S_p + C2*((V_d*(T(360)))/(p(360)*V(360)))*(p(360)-pm); % Average cylinder gas velocity
  124. Hc(i) = 3.26*B^-0.2*p(i-1)^0.8*(T(i-1))^-0.55*w^0.8; % Compute heat transfer coefficient for current temperature and pressure
  125. dQ(i) = Hc(i)*A*((T(i-1))-Tw)*dt; % Heat loss
  126. Cv(i) = Y1_tot(1)*CvNasa((T(i-1)),SpS(1)) + Y1_tot(2)*CvNasa((T(i-1)),SpS(2)) + Y1_tot(3)*CvNasa((T(i-1)),SpS(3)) + Y1_tot(4)*CvNasa((T(i-1)),SpS(4)) + Y1_tot(5)*CvNasa((T(i-1)),SpS(5)) + Y1_tot(6)*CvNasa((T(i-1)),SpS(6));
  127. % Specific heat capacity for current temperature
  128. m1(i) = m1(i-1); % Compute new mass
  129. T(i) = (T(i-1)) + (-dQ(i) - p(i-1)*(V(i)-V(i-1)))/(Cv(i)*m1(i)); % Compute new temperature using 1st law of thermodynamics
  130. p(i) = m1(i)*(Runiv/M_av1)*T(i)/V(i); % Compute new pressure using idael gas law
  131. end
  132. end
  133.  
  134. %% Combustion (OLD VERSION)
  135.  
  136. % pm(2) = (((r_c*V_d)/(r_c - 1))^k * P0)/(V(719)^k); % Compute motorized pressure for new phase
  137. %
  138. % C1 = 2.28; % Constant for w (velocity)
  139. % C2 = 3.24e-03; % Constant for w (velocity)
  140. % w = C1*S_p + C2*(((V(720)-V(719))*(T(719)))/(p(719)*V(719)))*(p(719)-pm(2));% Average cylinder velocity
  141. % Hc(720) = 3.26*B^-0.2*p(720-1)^0.8*(T(720-1))^-0.55*w^0.8; % Heat transfer coefficient
  142. % dQ(720) = Hc(720)*A*((T(720-1))-Tw)*dt; % Heat losses
  143. %
  144. % Qh = Y1_tot(1)*Qlhv_gas + Y1_tot(2)*Qlhv_eth -dQ(720); % Combustion heat
  145. %
  146. % Cv(720) = CvNasa((T(719)),SpS(1)); % Specific heat capacity for current temperature
  147. % m1(720) = Mtot2; % Current mass
  148. % M_av2 = sum(X2_tot.*Mi); % Molucular mass of gas mixture
  149. % T(720) = (Qh)/(sum(Y2_tot)*Cv(720))+T(719); % Compute current temperature
  150. % p(720) = m1(720)*(Runiv/M_av2)*T(720)/V(720); % Current pressure
  151.  
  152.  
  153. %% Expansion
  154. M_av2 = sum(X2_tot.*Mi); % Molucular mass of gas mixture
  155. C1 = 2.28; % Constant for w (velocity)
  156. C2 = 3.24e-03; % Constant for w (velocity)
  157.  
  158.  
  159. for i = 721:1080
  160. if i<= 770 % Last part of combustion
  161. Ca(i) = Ca(i-1)+dCa; % Compute new angle
  162. V(i) = Vcyl(Ca(i)); % Compute new volume
  163. Qc(i) = Qlhv_gas*Comb_rate_gas*dCa + Qlhv_eth*Comb_rate_eth*dCa; % Combustion heat
  164. pm = (((r_c*V_d)/(r_c - 1))^k * P0)/(V(721)^k); % Motorized pressure
  165. w = C1*S_p + C2*((V_d*(T(721-1)))/(p(721-1)*V(721)))*(p(721-1)-pm); % Average cylinder gas velocity
  166. Hc(i) = 3.26*B^-0.2*p(i-1)^0.8*T(i-1)^-0.55*w^0.8; % Compute heat transfer coefficient for current temperature and pressure
  167. dQ(i) = Hc(i)*A*(T(i-1)-Tw)*dt; % Heat losses
  168. m_added_CO2(i) = ((Comb_rate_gas*dCa*Mi(1)*8) + (Comb_rate_eth*dCa*Mi(2)*2))/Mi(4);
  169. % Change in mass of CO2 by combustion
  170. m_CO2(i) = sum(m_added_CO2); % Total mass of CO2
  171. m_added_H2O(i) = ((Comb_rate_gas*dCa*Mi(1)*9) + (Comb_rate_eth*dCa*Mi(2)*3))/Mi(5);
  172. % Change in mass of H2O by combustion
  173. m_H2O(i) = sum(m_added_H2O); % Total mass of H2O
  174. m_added_O2(i) = -((Comb_rate_gas*dCa*Mi(1)*12.5) + (Comb_rate_eth*dCa*Mi(2)*3))/Mi(3);
  175. % Change in mass of O2 by combustion
  176. m_O2(i) = m_O2(i-1) + sum(m_added_O2(i)); % Total mass of O2
  177. m_gas(i) = m_gas(i-1) - Comb_rate_gas*dCa; % Total mass of gasoline
  178. m_eth(i) = m_eth(i-1) - Comb_rate_eth*dCa; % Total mass of ethanol
  179. m_N2(i) = Y1_tot(6)*m1(1); % Total mass of N2
  180. m1(i) = m1(i-1)-Comb_rate_gas*dCa - Comb_rate_eth*dCa + m_added_CO2(i) + m_added_H2O(i) + m_added_O2(i);
  181. % Compute new mass
  182. Cv(i) = (m_gas(i)/m1(i))*CvNasa((T(i-1)),SpS(1)) + (m_eth(i)/m1(i))*CvNasa((T(i-1)),SpS(2)) + (m_O2(i)/m1(i))*CvNasa((T(i-1)),SpS(3)) + (m_CO2(i)/m1(i))*CvNasa((T(i-1)),SpS(4)) + (m_H2O(i)/m1(i))*CvNasa((T(i-1)),SpS(5)) + (m_N2(i)/m1(i))*CvNasa((T(i-1)),SpS(6));
  183. % Specific heat capacity for current temperature
  184. T(i) = T(i-1) + (-dQ(i) + Qc(i) - p(i-1)*(V(i)-V(i-1)))/(Cv(i)*m1(i)); % Compute new temperature
  185. p(i) = m1(i)*(Runiv/M_av2)*T(i)/V(i); % Compute new pressure
  186. else
  187. Ca(i)=Ca(i-1)+dCa; % Compute new angle
  188. V(i)=Vcyl(Ca(i)); % Compute new volume
  189. pm = (((r_c*V_d)/(r_c - 1))^k * P0)/(V(720)^k); % Motorized pressure
  190. w = C1*S_p + C2*(((V(720)-V(719))*T(720))/(p(720)*V(720)))*(p(720)-pm); % Average cylinder gas velocity
  191. Hc(i) = 3.26*B^-0.2*p(i-1)^0.8*T(i-1)^-0.55*w^0.8; % Compute heat transfer coefficient for current temperature and pressure
  192. dQ(i) = Hc(i)*A*(T(i-1)-Tw)*dt; % Heat lossed
  193. Cv(i) = Y2_tot(1)*CvNasa((T(i-1)),SpS(1)) + Y2_tot(2)*CvNasa((T(i-1)),SpS(2)) + Y2_tot(3)*CvNasa((T(i-1)),SpS(3)) + Y2_tot(4)*CvNasa((T(i-1)),SpS(4)) + Y2_tot(5)*CvNasa((T(i-1)),SpS(5)) + Y2_tot(6)*CvNasa((T(i-1)),SpS(6));
  194. % Specific heat capacity for current temperature
  195. m1(i) = m1(i-1); % Compute new mass
  196. T(i) = T(i-1) + (-dQ(i) - p(i-1)*(V(i)-V(i-1)))/(Cv(i)*m1(i)); % Compute new temperature
  197. p(i) = m1(i)*(Runiv/M_av2)*T(i)/V(i); % Compute new pressure
  198. end
  199. end
  200.  
  201. %% Exhaust
  202. for i=1080:1440
  203. Ca(i)=Ca(i-1)+dCa; % Compute new angle
  204. V(i)=Vcyl(Ca(i)); % Compute new volume
  205. T(i)=T0; % Compute new temperature
  206. p(i)=P0; % Compute new pressure
  207. end
  208.  
  209. %% Plot
  210. figure(1) % Create figure
  211. plot(V,p) % Create plot
  212. xlabel('Volume [m^3]') % xLabel figure
  213. ylabel('Pressure [Pa]') % yLabel figure
  214. title('p-V diagram')
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!