OpenWAM
TIsoSpeedLine.cpp
1 /* --------------------------------------------------------------------------------*\
2 ==========================|
3  \\ /\ /\ // O pen | OpenWAM: The Open Source 1D Gas-Dynamic Code
4  \\ | X | // W ave |
5  \\ \/_\/ // A ction | CMT-Motores Termicos / Universidad Politecnica Valencia
6  \\/ \// M odel |
7  ----------------------------------------------------------------------------------
8  License
9 
10  This file is part of OpenWAM.
11 
12  OpenWAM is free software: you can redistribute it and/or modify
13  it under the terms of the GNU General Public License as published by
14  the Free Software Foundation, either version 3 of the License, or
15  (at your option) any later version.
16 
17  OpenWAM is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  GNU General Public License for more details.
21 
22  You should have received a copy of the GNU General Public License
23  along with OpenWAM. If not, see <http://www.gnu.org/licenses/>.
24 
25 
26  \*-------------------------------------------------------------------------------- */
27 
28 // ---------------------------------------------------------------------------
29 #pragma hdrstop
30 
31 #include "TIsoSpeedLine.h"
32 //#include <cmath>
33 #include "Globales.h"
34 
35 // ---------------------------------------------------------------------------
36 
37 TIsoSpeedLine::TIsoSpeedLine() {
38  iStator = NULL;
39  iRotor = NULL;
40  iEfficiency = NULL;
41 }
42 
43 TIsoSpeedLine::~TIsoSpeedLine() {
44  delete iStator;
45  delete iRotor;
46  delete iEfficiency;
47 }
48 ;
49 
50 void TIsoSpeedLine::ReadIsoSpeed(int points, FILE *Input) {
51  double ER = 0., MF = 0., EF = 0.;
52  FNumDatos = points;
53  for(int i = 0; i < points; i++) {
54  fscanf(Input, "%lf %lf %lf", &MF, &ER, &EF);
55  FReducedAirMassFlow.push_back(MF);
56  FExpansionRatio.push_back(ER);
57  FEfficiency.push_back(EF);
58  }
59 }
60 
61 void TIsoSpeedLine::AsignaValores(double ER, double MF, double EF) {
62  FReducedAirMassFlow.push_back(MF);
63  FExpansionRatio.push_back(ER);
64  FEfficiency.push_back(EF);
65 }
66 
67 void TIsoSpeedLine::EffectiveSection(double Area, bool CalculaGR, double Angle, double Diam1, double Diam2,
68  double Diam3, double n_limit) {
69  double tmp1 = 0., tmp2 = 0.;
70  double FGamma = 1.40, FR = 287;
71  double T00_T0, P00_P0, P2_P0, T2_T0, f_P2_P0, GR, nN, n, P2_P00, T2_T00, T1_T0, gG, g, kK, k, P1_P0, P00_P1, P1_P2,
72  raiZ;
73 
74  FNumDatos = FReducedAirMassFlow.size();
75 
76  for(int i = 0; i < FNumDatos; i++) {
77 
78  tmp1 = 1;
79  tmp2 = 1;
80  do {
81  tmp1 = tmp2;
82  tmp2 = 1 + ((FGamma - 1) / 2) * (FR / FGamma) * (pow2(FReducedAirMassFlow[i] / 1000000) / pow2(Area)) * pow((tmp1),
83  ((FGamma + 1) / (FGamma - 1)));
84  } while(fabs(tmp1 - tmp2) / tmp1 > 1e-12);
85  T00_T0 = tmp2;
86  P00_P0 = pow(T00_T0, FGamma / (FGamma - 1));
87  P2_P0 = (1 / FExpansionRatio[i]) * pow(T00_T0, (FGamma / (FGamma - 1)));
88  T2_T0 = 1 - (FEfficiency[i]) * (1 - pow((1 / FExpansionRatio[i]), ((FGamma - 1) / FGamma))) * T00_T0;
89  f_P2_P0 = FExpansionRatio[i] * (1 / T00_T0 - FEfficiency[i] * (1 - pow((1 / FExpansionRatio[i]),
90  ((FGamma - 1) / FGamma))));
91  if(CalculaGR) {
92  GR = 1 - (((2 * FR * tan(__units::DegToRad(Angle))) / (Diam1 * pow2(Diam2) * pow2(__cons::Pi))) * ((
93  FReducedAirMassFlow[i] / 1000000) / FSpeed) * f_P2_P0);
94  // new code --> if the reaction degree is lower than 0.4 then force it to 0.4 value instead of calculating lower values or even negative values
95  if(GR < 0.5) {
96  GR = 0.5;
97  }
98  } else {
99  GR = 0.5;
100  }
101  nN = log(P2_P0) / log(T2_T0);
102  n = nN / (nN - 1);
103  P2_P00 = P2_P0 / P00_P0;
104  T2_T00 = T2_T0 / T00_T0;
105  T1_T0 = 1 + T00_T0 * (GR - 1) * FEfficiency[i] * (1 - pow(P2_P00, (FGamma - 1) / FGamma));
106  if(n > n_limit) {
107  gG = ((FGamma / (FGamma - 1)) - nN * (log(T2_T00) + log(T00_T0)) / log(T1_T0)) / (1 - (log(T2_T00) + log(T00_T0)) / log(
108  T1_T0));
109  g = (n / (n - n_limit) + (gG / (gG - 1)) / (FGamma - n)) / (1 / (n - n_limit) + 1 / (FGamma - n));
110  } else {
111  gG = 0.;
112  g = (FGamma / (n - 1) + n / (n_limit - n)) / (1 / (n - 1) + 1 / (n_limit - n));
113  }
114  kK = (g / (g - 1)) + (nN - (g / (g - 1))) * (log(T2_T00) + log(T00_T0)) / log(T1_T0);
115  k = kK / (kK - 1);
116  P1_P0 = pow(1 + T00_T0 * (GR - 1) * FEfficiency[i] * (1 - pow(P2_P00, ((FGamma - 1) / FGamma))), kK);
117  P00_P1 = P00_P0 / P1_P0;
118  P1_P2 = P1_P0 / P2_P0;
119  raiZ = (2 / (FGamma - 1)) * (1 - pow(1 / P00_P1, ((FGamma - 1) / FGamma)));
120  if(raiZ < 0)
121  raiZ = 0;
122  else
123  raiZ = sqrt(raiZ);
124  // raiZ = sqrt((2 / (FGamma - 1)) * (1 - pow(1 / P00_P1 , ((FGamma - 1) / FGamma))));
125  StatorEffectiveSection.push_back(FReducedAirMassFlow[i] * 1e-6 * pow(P00_P1,
126  (1 / FGamma)) * (sqrt(FR / FGamma)) / raiZ);
127  // raiZ = sqrt((2 / (FGamma - 1)) * (1 - pow(1 / P1_P2 , ((FGamma - 1) / FGamma))));
128  raiZ = (2 / (FGamma - 1)) * (1 - pow(1 / P1_P2, ((FGamma - 1) / FGamma)));
129  if(raiZ < 0)
130  raiZ = 0;
131  else
132  raiZ = sqrt(raiZ);
133 // RotorEffectiveSection.push_back
134 // (FReducedAirMassFlow[i] * 1e-6 * P00_P1 * pow(P1_P2,
135 // (1 / FGamma)) * (sqrt(FR / FGamma)) / raiZ * sqrt
136 // ((1 / T1_T0) - (T00_T0 / T1_T0) + (T2_T00 * T1_T0 / T00_T0)));
137  RotorEffectiveSection.push_back(FReducedAirMassFlow[i] * 1e-6 * P00_P1 * pow(P1_P2,
138  (1 / FGamma)) * (sqrt(FR / FGamma)) / raiZ);
139  }
140 }
141 
142 void TIsoSpeedLine::CalculatePower(double Tin) {
143 
144  double Cp = 1004;
145  double R = 287;
146  double gamma = Cp / (Cp - R);
147  double gam = (1 - gamma) / gamma;
148  for(int i = 0; i < FNumDatos; i++) {
149  FPower.push_back(FReducedAirMassFlow[i] * Cp * FEfficiency[i] * sqrt(Tin) * FExpansionRatio[i] / 10 * (1 - pow(
150  FExpansionRatio[i], gam)));
151  // printf("%4.2lf\t", FPower[i]);
152  }
153  FPowerMin = FPower.front();
154  FPowerMax = FPower.back();
155  // printf("\n");
156 }
157 
158 void TIsoSpeedLine::PutSpeed(double sp) {
159  FSpeed = sp;
160 }
161 
162 void TIsoSpeedLine::Adimensionaliza() {
163  double DeltaPreMax = FExpansionRatio.back() - FExpansionRatio.front();
164  double m = 0., Rtc = 0.;
165 
166  for(int i = 0; i < FNumDatos; ++i) {
167  FExpansionRatioAdim.push_back((FExpansionRatio[i] - FExpansionRatio.front()) / DeltaPreMax);
168 
169  }
170  FERMin = FExpansionRatio.front();
171  FERMax = FExpansionRatio.back();
172 
173  iStator = new Hermite_interp(FExpansionRatioAdim, StatorEffectiveSection);
174 
175  iRotor = new Hermite_interp(FExpansionRatioAdim, RotorEffectiveSection);
176 
177  iEfficiency = new Hermite_interp(FExpansionRatioAdim, FEfficiency);
178 
179 }
180 
181 void TIsoSpeedLine::GetAdiabaticEfficiency(TTC_HTM *HTM, double TinT, double TinC) {
182 
183 #ifdef tchtm
184  double m = 0., Rtc = 0.;
185 
186  for(int i = 0; i < FNumDatos; ++i) {
187  m = FReducedAirMassFlow[i] / sqrt(TinT) * FExpansionRatio[i] / 10;
188  Rtc = FSpeed * 60 * sqrt(TinT);
189  FEfficiency[i] = HTM->CorrectTurbineMap(m, FExpansionRatio[i],
190  FEfficiency[i], TinC, TinT, Rtc);
191  }
192  if(iEfficiency != NULL)
193  delete iEfficiency;
194  iEfficiency = new Hermite_interp(FExpansionRatioAdim, FEfficiency);
195 #endif
196 }
197 
198 double TIsoSpeedLine::Stator(double ERAdim) {
199  double ret_val = iStator->interp(ERAdim);
200 
201  return ret_val;
202 }
203 
204 double TIsoSpeedLine::Rotor(double ERAdim) {
205  double ret_val = iRotor->interp(ERAdim);
206 
207  return ret_val;
208 }
209 
210 double TIsoSpeedLine::Efficiency(double ERAdim) {
211  double ret_val = iEfficiency->interp(ERAdim);
212 
213  return ret_val;
214 }
215 
216 void TIsoSpeedLine::PrintEffectiveSection(FILE * fich) {
217  for(dVector::size_type i = 0; i < FExpansionRatio.size(); i++) {
218  fprintf(fich, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", FSpeed, FExpansionRatio[i], FReducedAirMassFlow[i], FEfficiency[i],
219  FExpansionRatioAdim[i], StatorEffectiveSection[i],
220  RotorEffectiveSection[i]);
221  }
222 }
223 
224 double TIsoSpeedLine::Speed() {
225  return FSpeed;
226 }
227 
228 double TIsoSpeedLine::ERMax() {
229  return FERMax;
230 }
231 
232 double TIsoSpeedLine::ERMin() {
233  return FERMin;
234 }
235 
236 #pragma package(smart_init)
TTC_HTM
Definition: TTC_HTM.h:307
Hermite_interp
Definition: Math_wam.h:311
pow2
T pow2(T x)
Returns x to the power of 2.
Definition: Math_wam.h:88