fork download
copy 
  1. #include<iostream>
  2. #include<cmath>
  3. using namespace std;
  4.  
  5. bool between(double a, double b, double c)
  6. {
  7. 	if(c > b)
  8. 	{
  9. 		if((a>b) && (a<=c))
  10. 			return true;
  11. 		else
  12. 			return false;
  13. 	}
  14. 	else
  15. 	{
  16. 		if((a>b) || (a<c))
  17. 			return true;
  18. 		else
  19. 			return false;
  20. 	}
  21. 	return false;
  22. }
  23.  
  24. int main()
  25. {
  26. 	double RPS, dT, dTH, THio, THic, dAi, Aimax, THeo, THec, dAe, Aemax;
  27. 	double dMi, dMe, M, TH, T, Ai, Ae, Pi, Pe, MWi, MWb, TPi, P, MW, TP;
  28. 	double V, Vmax, Vmin, G, TPh2o, dQb, dMW, s, dQh2o, k, RHOi, Aface;
  29. 	double THbs, THbe, RHO, tmp, Mn, dMio, dMeo, dTP, TPold, Vold, exp, Aio;
  30. 	double Conv = .0174533, sum, Vnt, Gin, Gex, CVin, CVex, CV, Str,Y,C, sum2;
  31. 	bool burn;
  32. 	int Rotors, Rotations, next, sweep, incr, start;
  33.  
  34. 	Rotors = 2;
  35. 	Rotations = 20;
  36. 	sweep = 9000;
  37. 	start = 3000;
  38. 	incr = 300;
  39. 	RPS = 6000/180;
  40. 	dTH = .02;
  41. 	dT = dTH / (360 * RPS);
  42. 	THio = 1;
  43. 	THic = 72;
  44. 	Aio = .01;
  45. 	Aimax = .01310644;
  46. 	dAi = Aimax / 40;
  47. 	THeo = 253;
  48. 	THec = 329;
  49. 	Aemax = .00680644;
  50. 	dAe = Aemax / 30;
  51. 	Pi = 102000;           //approx 2.7 bar , 1.7 boost
  52. 	Pe = 101325;           //atmospheric, no turbo
  53. 	MWi = .02893;          //kg per mole of air
  54. 	MWb = .02832/1.0714;          //final molecular weight after burn
  55. 	TPi = 290;          //temp in
  56. 	Vmax = .654;            //displacement
  57. 	Vmin = Vmax / 10;       //displacement over compression
  58. 	Gin = 1.4;               //gamma of air
  59. 	Gex = 1.24;
  60. 	TPh2o = 350;           //engine water temp, boil = 373
  61. 	dQb = 2125;          //heating rate constant  dTP = dQb * dT / M
  62. 	dMW = (MWb - MWi)/1;  // same except molecular weight change
  63. 	dQh2o = 8.978 / 250;            //dTP = dQ * dT * deltaTP / M
  64. 	THbs = 190;
  65. 	THbe = 240;
  66. 	RHOi = (MWi * Pi) / (8309.28 * TPi);
  67. 	CVin = .71771;
  68. 	CVex = 1.22774;
  69. 	Aface = .0145;
  70. 	Str = .015;
  71. 	C = .74;
  72. 	Y = .863;
  73.  
  74. 	int n = 360/dTH;
  75. 	double Pressure[n], Temp[n], mDoti[n], mDote[n], Press[Rotations], dU[n];
  76.  
  77.  
  78. for(int k = start; k <= sweep; k+=incr)
  79. {
  80. 	RPS = k/180;
  81. 	dT = dTH / (360 * RPS);
  82. 	dMi = 0;
  83. 	dMe = 0;
  84. 	P = (Pe+Pi) / 2;
  85. 	TP = 1000;
  86. 	MW = MWb;
  87. 	RHO = (MW * P) / (8309.28 * TP);
  88. 	M = Vmin * RHO;
  89. 	V = Vmin;
  90. 	Ai = 0;
  91. 	Ae = 0;
  92. 	for(int i = 0; i < Rotations; i++)
  93. 	{
  94. 		for(int j = 0; j < n; j++)
  95. 		{
  96. 			TH = j * dTH;
  97. 			T = dT * j;
  98. 			burn = false;
  99. 			if(between(TH, THio, THic))
  100. 			{
  101. 				Ai+=(dAi * dTH);
  102. 				if(Ai > Aimax)
  103. 					Ai = Aimax;
  104. 			}
  105. 			else
  106. 			{
  107. 				Ai-=(dAi * dTH);
  108. 				if(Ai < 0)
  109. 					Ai = 0;
  110. 			}
  111. 			if(between(TH, THeo, THec))
  112. 			{
  113. 				Ae+=(dAe * dTH);
  114. 				if(Ae > Aemax)
  115. 					Ae = Aemax;
  116. 			}
  117. 			else
  118. 			{
  119. 				Ae-=(dAe * dTH);
  120. 				if(Ae < 0)
  121. 					Ae = 0;
  122. 			}
  123. 			RHO = (MW * P) / (8309.28 * TP);
  124. 			dMio = dMi;
  125. 			dMeo = dMe;
  126. 			//dMi = (Pi - P)*C * dT * Y * Ai *sqrt(2 * RHOi /abs(Pi - P));
  127. 			dMi = (Pi - P) * dT * Ai * sqrt((2 * RHOi)/abs(Pi - P));
  128. 			dMe = (Pe - P) * dT * Ae * sqrt((2 * RHO)/abs(Pe - P));
  129. 			dMio = (dMio + dMi) / 2;
  130. 			dMeo = (dMeo + dMe) / 2;
  131. 			Mn = M + dMeo + dMio;
  132. 			tmp = cos(TH * Conv) * cos(TH * Conv);
  133. 			Vold = V;
  134. 			V = Vmin * tmp + Vmax * (1-tmp);
  135. 			//cout<<MW<<" "<<dMeo<<" "<<M<<" "<<dMio<<" "<<Mn<<endl;
  136. 			MW = (MW * (M + dMeo) + MWi * dMio) / Mn;
  137. 			if(between(TH, THbs, THbe) && (MW > MWb))
  138. 			{
  139. 				MW += (dMW * dT) / Mn;
  140. 				burn = true;
  141. 			}
  142. 			TPold = TP;
  143. 			//cout<<V<<" "<<tmp<<" "<<exp<<" "<<pow(exp, G-1)<<endl;
  144. 			TP = (TP * (M + dMeo) + TPi * dMio) / Mn;
  145. 			G = ((MWi-MW)*Gex + (MW - MWb)*Gin)/(MWi - MWb);
  146. 			exp = V / Vold;
  147. 			CV = ((MWi-MW)*CVex + (MW - MWb)*CVin)/(MWi - MWb);
  148. 			dU[j] = CV * M * TP * (1 - (1/pow(exp, G-1)));
  149. 			Vold *= (1 + dMio / M);
  150. 			tmp = V * (1 - dMeo / M);
  151. 			exp = tmp / Vold;
  152. 			TP /= pow(exp, G-1);
  153.  
  154. 			//dU[j] = CV*(TPold-TP)*M;
  155. 			dTP = 0;
  156. 			if(burn)
  157. 				dTP += (dQb * dT) / Mn;
  158. 			TP += dTP;
  159. 			tmp = dQh2o * dT;
  160. 			TP += tmp * (((TPh2o - TPold) / M) + ((TPh2o - TP) / Mn)) / 2;
  161. 			M = Mn;
  162. 			P = (M * TP * 8309.28) / (MW * V);
  163.  
  164. 			Pressure[j] = P;
  165. 			//cout<<P<<" "<<V<<" "<<TP<<" "<<M<<endl;
  166. 			Temp[j] = TP;
  167. 			mDoti[j] = dMi;
  168. 			mDote[j] = dMe;
  169. 		}
  170. 		Press[i] = P;
  171.  
  172. 	}
  173.  
  174. 	/*sum = 0;
  175. 	for(int i = 0; i < n; i++)
  176. 	{
  177. 		sum+=dU[i];
  178. 	}
  179. 	sum*=(RPS*3*Rotors)/.745;
  180. 	cout<<k<<" "<<sum;
  181. 	sum*=745;
  182. 	sum/=(RPS * 6 * 3.14159);
  183.  
  184. 	cout<<" "<<sum;*/
  185.  
  186. 	sum = 0;
  187. 	for(int i = 0; i < n; i++)
  188. 	{
  189. 		sum+=Pressure[i] * Aface * Str * sin(2 * dTH * i * Conv);
  190. 	}
  191. 	sum/=n;
  192. 	sum*=3;
  193.  
  194. 	cout<<k<<" "<<sum;
  195. 	sum*=(6 * 3.14159 * RPS);
  196. 	sum/=745;
  197.  
  198. 	cout<<" "<<sum;
  199.  
  200. 	sum = 0;
  201. 	sum2=0;
  202. 	for(int i = 0; i < n; i++)
  203. 	{
  204. 		sum+=Temp[i] * (-mDote[i]);
  205. 		sum2 += (-mDote[i]);
  206. 	}
  207. 	sum/=sum2;
  208. 	cout<<" "<<sum<<"K"<<endl;
  209. }	
  210.  
  211. 	return 0;
  212.  
  213. }
Success #stdin #stdout 1.24s 5284KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
3000 89.547 36.2506 691.574K
3300 91.2524 41.5586 738.964K
3600 91.7672 46.4368 779.879K
3900 91.3633 48.544 796.437K
4200 90.9209 52.9098 831.699K
4500 89.8933 56.8607 863.546K
4800 88.9585 58.5201 876.745K
5100 87.1748 61.758 903.819K
5400 84.8151 64.3782 926.794K
5700 83.3837 65.4015 936.296K
6000 80.4309 67.1554 954.186K
6300 77.4942 68.6249 971.525K
6600 75.8638 69.1005 978.499K
6900 72.8641 70.0555 993.887K
7200 69.4396 70.2768 1004.12K
7500 67.9833 70.5229 1011.05K
7800 64.6214 70.3055 1019.33K
8100 61.3883 69.8944 1026.94K
8400 59.894 69.7085 1031.26K
8700 56.928 69.1372 1038.59K
9000 53.8709 68.1504 1042.68K
https://ideone.com/SqwOYW
language:
C++ (gcc 8.3)
created:
5 days ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!