OpenWAM
TCompTubos.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 #include <iostream>
31 #ifdef __BORLANDC__
32 #include <vcl.h>
33 #endif
34 //#include <cmath>
35 
36 #include "TCompTubos.h"
37 #include "TCondicionContorno.h"
38 #include "TTubo.h"
39 
40 // ---------------------------------------------------------------------------
41 // ---------------------------------------------------------------------------
42 
43 TCompTubos::TCompTubos(int i, nmTipoCalculoEspecies SpeciesModel, int numeroespecies, nmCalculoGamma GammaCalculation,
44  bool ThereIsEGR) :
45  TCompresor(i, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
46  FModeloCompresor = nmCompPipes;
47  Mapa2T = new TMapaComp2Tub(FNumeroCompresor);
48  FPrimerCiclo = true;
49  FCoefPresiones = 0.10344;
50  FGastoOut = 0.;
51  FFlowIn = 0;
52  FFlowOut = 0;
53 }
54 
55 // ---------------------------------------------------------------------------
56 // ---------------------------------------------------------------------------
57 
58 TCompTubos::~TCompTubos() {
59 
60 }
61 
62 // ---------------------------------------------------------------------------
63 // ---------------------------------------------------------------------------
64 
65 void TCompTubos::CondicionCompresor(double Theta, stTuboExtremo *TuboExtremo, double TiempoActual, int TuboCalculado) {
66 
67  double ErrorVelMax = 0., ErrorVelMin = 0., ErrorCheck = 0., RC = 0., RD = 0., CP = 0.;
68  double temp = 0.;
69  double CoefPresiones_anterior, Rendimiento_ant = 0.;
70  double FraccionMasicaAcum = 0.;
71  double vout = 0., aaout = 0.;
72 
73  try {
74  // FContadorCheckSentido=0;
75  CoefPresiones_anterior = FCoefPresiones;
76  Rendimiento_ant = FRendimiento;
77 
78  FGamma = (FTuboRot->GetGamma(FNodoFinRotor) + FTuboEst->GetGamma(FNodoFinStator)) / 2.;
79  // Se ha de hacer una media de entrada y salida
80  FRMezcla = (FTuboRot->GetRMezcla(FNodoFinRotor) + FTuboEst->GetRMezcla(FNodoFinStator)) / 2.;
81  FCpMezcla = FGamma * FRMezcla / (FGamma - 1);
82  FGamma1 = __Gamma::G1(FGamma);
83  FGamma2 = __Gamma::G2(FGamma);
84  FGamma3 = __Gamma::G3(FGamma);
85  FGamma4 = __Gamma::G4(FGamma);
86  FGamma5 = __Gamma::G5(FGamma);
87 
88  FDeltaTiempo = TiempoActual - FTiempo0;
89  FTiempo0 = TiempoActual;
90 
91  if(FPrimerCiclo) {
92 
93  double vR = FSignoRotor * (TuboExtremo[FPosicionRotor].Landa - TuboExtremo[FPosicionRotor].Beta) / FGamma3 *
94  __cons::ARef;
95  double aR = (TuboExtremo[FPosicionRotor].Landa + TuboExtremo[FPosicionRotor].Beta) / 2. * __cons::ARef;
96  double pR = __units::BarToPa(pow(aR / TuboExtremo[FPosicionRotor].Entropia / __cons::ARef, FGamma4));
97  // double dR=FGamma*pR/aR/aR;
98 
99  double aS = (TuboExtremo[FPosicionEstator].Landa + TuboExtremo[FPosicionEstator].Beta) / 2. * __cons::ARef;
100  double vS = FSignoStator * (TuboExtremo[FPosicionEstator].Landa - TuboExtremo[FPosicionEstator].Beta) / FGamma3 *
101  __cons::ARef;
102  double pS = __units::BarToPa(pow(aS / TuboExtremo[FPosicionEstator].Entropia / __cons::ARef, FGamma4));
103 
104  FAsonidoIn = aR;
105  FAsonidoOut = aS;
106  FVelocidadIn = vR;
107  FVelocidadOut = vS;
108  FAout = aS;
109  FPreTotalOut = pS;
110 
111  FPresionDep = __units::PaToBar(pS);
112  FAsonidoDep = aS;
113  FTempDep = aS * aS / FGamma / FRMezcla;
114  FVolumen = FLongitudCaract * FAreaStator;
115  FMasaDep = pS / FRMezcla / FTempDep * FVolumen;
116 
117  FRendimiento = Mapa2T->EvaluaRendSplines(0.);
118  FRelacionCompresion = pS * pow(1 + FGamma3 * vS * vS / aS / aS,
119  FGamma / FGamma1) / pR / pow(1 + FGamma3 * vR * vR / aR / aR, FGamma / FGamma1);
120  FCoefPresiones = (pow(FRelacionCompresion, FGamma1 / FGamma) - 1) / FRendimiento;
121 
122  if(FExtremoRotor == nmLeft) {
123  FCarCIn = &(TuboExtremo[FPosicionRotor].Beta);
124  FAaIn = &(TuboExtremo[FPosicionRotor].Entropia);
125  FCarDIn = &(TuboExtremo[FPosicionRotor].Landa);
126  } else {
127  FCarCIn = &(TuboExtremo[FPosicionRotor].Landa);
128  FAaIn = &(TuboExtremo[FPosicionRotor].Entropia);
129  FCarDIn = &(TuboExtremo[FPosicionRotor].Beta);
130  }
131  if(FExtremoStator == nmRight) {
132  FCarCOut = &(TuboExtremo[FPosicionEstator].Landa);
133  FAaOut = &(TuboExtremo[FPosicionEstator].Entropia);
134  FCarDOut = &(TuboExtremo[FPosicionEstator].Beta);
135  } else {
136  FCarCOut = &(TuboExtremo[FPosicionEstator].Beta);
137  FAaOut = &(TuboExtremo[FPosicionEstator].Entropia);
138  FCarDOut = &(TuboExtremo[FPosicionEstator].Landa);
139  }
140  FAreaIn = FAreaRotor;
141  FAreaOut = FAreaStator;
142  FPrimerCiclo = false;
143  FCoefPresiones = pow2(*FAaOut / *FAaIn) * pow(FRelacionCompresion, FGamma1 / FGamma) - 1;
144  }
145 
146  if(TuboCalculado == 10000) {
147 
148  SolveInletBoundary(FA1in, FU1in, FA1out, FU1out);
149 
150  double AA1fin = 0.;
151 
152  if(!FFlujoDirecto) {
153  AA1fin = FAsonidoDep / __cons::ARef * pow(FPresionDep, -FGamma5) * pow(FRelacionCompresion,
154  FGamma5) / sqrt(1 + FCoefPresiones);
155  } else {
156  AA1fin = *FAaIn;
157  }
158  FFlowIn = FGamma * __units::BarToPa(pow(FA1in * __cons::ARef / pow(AA1fin * __cons::ARef, FGamma),
159  1 / FGamma3)) * FU1in * __cons::ARef * FAreaIn;
160 
161  // NewPropertiesInTheVolume();
162 
163  FAsonidoIn = FA1in * __cons::ARef;
164  FVelocidadIn = FU1in * __cons::ARef;
165 
166  SolveOutletBoundary(FA2, FU2);
167 
168  double AA2fin = 0.;
169  if(FFlowOut < 0) {
170  AA2fin = FAsonidoDep / __cons::ARef * pow(FPresionDep, -FGamma5);
171  } else {
172  AA2fin = *FAaOut;
173  }
174 
175  NewPropertiesInTheVolume();
176  FAsonidoOut = FA2 * __cons::ARef;
177  FVelocidadOut = FU2 * __cons::ARef;
178  // <--- NUEVO MODELO
179 
180  *FCarDIn = (FAsonidoIn - FGamma1 / 2. * FVelocidadIn) / __cons::ARef;
181  if(FVelocidadIn < 0) {
182  *FAaIn = AA1fin;
183  }
184  *FCarDOut = (FAsonidoOut + FGamma1 / 2. * FVelocidadOut) / __cons::ARef;
185  if(FVelocidadOut > 0) {
186  *FAaOut = AA2fin;
187 
188  }
189 
190  } else if(TuboCalculado == FPosicionRotor) {
191  // NUEVO MODELO ---->
192  SolveInletBoundary(FA1in, FU1in, FA1out, FU1out);
193 
194  double AA1fin = 0.;
195 
196  if(!FFlujoDirecto) {
197  AA1fin = FAsonidoDep / __cons::ARef * pow(FPresionDep, -FGamma5) * pow(FRelacionCompresion,
198  FGamma5) / sqrt(1 + FCoefPresiones);
199  } else {
200  AA1fin = *FAaIn;
201  }
202  FFlowIn = FGamma * __units::BarToPa(pow(FA1in * __cons::ARef / pow(AA1fin * __cons::ARef, FGamma),
203  1 / FGamma3)) * FU1in * __cons::ARef * FAreaIn;
204 
205  NewPropertiesInTheVolume();
206 
207  FAsonidoIn = FA1in * __cons::ARef;
208  FVelocidadIn = FU1in * __cons::ARef;
209  // <--- NUEVO MODELO
210 
211  *FCarDIn = (FAsonidoIn - FGamma1 / 2. * FVelocidadIn) / __cons::ARef;
212  if(FVelocidadIn < 0) {
213  *FAaIn = AA1fin;
214  }
215  } else {
216  // NUEVO MODELO ---->
217  SolveOutletBoundary(FA2, FU2);
218 
219  double AA2fin = 0.;
220  if(FFlowOut < 0) {
221  AA2fin = FAsonidoDep / __cons::ARef * pow(FPresionDep, -FGamma5);
222  } else {
223  AA2fin = *FAaOut;
224  }
225 
226  NewPropertiesInTheVolume();
227  FAsonidoOut = FA2 * __cons::ARef;
228  FVelocidadOut = FU2 * __cons::ARef;
229  // <--- NUEVO MODELO
230 
231  *FCarDOut = (FAsonidoOut + FGamma1 / 2. * FVelocidadOut) / __cons::ARef;
232  if(FVelocidadOut > 0) {
233  *FAaOut = AA2fin;
234  }
235 
236  }
237 
238  FGasto1 = FGamma * __units::BarToPa(pow(FAsonidoIn / pow(*FAaIn * __cons::ARef, FGamma),
239  1 / FGamma3)) * FVelocidadIn * FAreaIn;
240 
241  FTempIn = FAsonidoIn * FAsonidoIn / (FGamma * FRMezcla);
242  FTempTotalIn = FTempIn + FVelocidadIn * FVelocidadIn / 2 / FCpMezcla;
243  FPresionIn = __units::BarToPa(pow(*FAaIn * __cons::ARef / FAsonidoIn, -FGamma4));
244  FPreTotalIn = FPresionIn * pow(FTempTotalIn / FTempIn, FGamma / FGamma1);
245  FTempOut = FAsonidoOut * FAsonidoOut / (FGamma * FRMezcla);
246  FTempTotalOut = FTempOut + FVelocidadOut * FVelocidadOut / 2 / FCpMezcla;
247  FPresionOut = __units::BarToPa(pow(*FAaOut * __cons::ARef / FAsonidoOut, -FGamma4));
248  FPreTotalOut = FPresionOut * pow(FTempTotalOut / FTempOut, FGamma / FGamma1);
249 
250  FGastoCorregido = FGasto1 * sqrt(FTempTotalIn / Mapa2T->getTempRef()) / (FPreTotalIn / Mapa2T->getPresionRef());
251 
252  // Busqueda de relacion de compresion en el compresor dado por el mapa
253  if(FGastoCorregido < -0.07) {
254  RC = Mapa2T->EvaluaRCHermite(Mapa2T->getGastoBombeo());
255  } else if(FGastoCorregido >= -0.07 && FGastoCorregido < Mapa2T->getGastoBombeo()) {
256  RC = Mapa2T->EvaluaRCHermite(FGastoCorregido);
257  } else if(FGastoCorregido >= Mapa2T->getGastoBombeo() && FGastoCorregido <= Mapa2T->getGastoRelComp1()) {
258  RC = Mapa2T->EvaluaRCHermite(FGastoCorregido);
259  } else if(FGastoCorregido > Mapa2T->getGastoRelComp1()) {
260  RC = Mapa2T->EvaluaRCHermite(Mapa2T->getGastoRelComp1());
261  }
262 
263  if(RC < 1.)
264  RC = 1.;
265  FRelacionCompresion = FRelacionCompresion + (1 / FLongitudCaract) * (RC - FRelacionCompresion) * 0.5 * fabs(
266  FVelocidadIn + FVelocidadOut) * FDeltaTiempo;
267 
268  // Busqueda del rendimiento en el compresor dado por el mapa
269  if(FGastoCorregido < -0.07) {
270  RD = Mapa2T->EvaluaRendSplines(Mapa2T->getGastoBombeo());
271  } else if(FGastoCorregido >= -0.07 && FGastoCorregido < Mapa2T->getGastoBombeo()) {
272  RD = Mapa2T->EvaluaRendSplines(Mapa2T->getGastoBombeo());
273  } else if(FGastoCorregido >= Mapa2T->getGastoBombeo() && FGastoCorregido <= Mapa2T->getGastoRelComp1()) {
274  RD = Mapa2T->EvaluaRendSplines(FGastoCorregido);
275  } else if(FGastoCorregido > Mapa2T->getGastoRelComp1()) {
276  RD = Mapa2T->EvaluaRendSplines(Mapa2T->getGastoRelComp1());
277  }
278 
279  if(RD < 0.4)
280  RD = 0.4;
281 
282  FRendimiento = RD;
283 
284  FCoefPresiones = (pow(FRelacionCompresion, FGamma1 / FGamma) - 1) / FRendimiento;
285 
286  if(FDeltaTiempo > 0) {
287 
288  FTrabajo = (FTempTotalOut - FTempTotalIn) * FRMezcla * FGamma4 / 2. * fabs(FGasto1) * FDeltaTiempo;
289  FPotencia = FTrabajo / FDeltaTiempo;
290 
291  FTrabajoPaso += FTrabajo;
292  FDeltaTPaso += FDeltaTiempo;
293 
294  FGastoCompresor = FGasto1;
295  FTemperatura10 = FTempTotalIn;
296  }
297 
298  FRegimenCorregido = Mapa2T->getRegimenCorregido();
299  } catch(exception & N) {
300  std::cout << "ERROR: CalculaCondicionContorno en el compresor: " << FNumeroCompresor << std::endl;
301  std::cout << "Tipo de error: " << N.what() << std::endl;
302  throw Exception("ERROR: CalculaCondicionContorno en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
303  }
304 }
305 
306 // ---------------------------------------------------------------------------
307 // ---------------------------------------------------------------------------
308 
309 void TCompTubos::Biseccion(double *VelIn, double *VelOut, double *AIn, double *AOut, double CarIn, double AaIn,
310  double CarOut, int sig) {
311  try {
312  bool Primera = true;
313  FVelInMin = 0;
314  FVelInMax = (double) sig * 2 * __cons::ARef * CarIn / FGamma2;
315 
316  FErrorVelIn = 10000.;
317  FErrorCarIn = 1e-6;
318  FCuentaVelIn = 0.;
319 
320  while(fabs(FErrorVelIn) > FErrorCarIn && FCuentaVelIn < 200) {
321  *VelIn = (FVelInMin + FVelInMax) / 2.;
322 
323  *AIn = __cons::ARef * CarIn - (double) sig * 0.5 * FGamma1 * *VelIn;
324  FPresionIn = __units::BarToPa(pow(AaIn / *AIn, -FGamma4));
325  FTempIn = pow2(*AIn) / FGamma / FRMezcla;
326  FDensidadIn = FPresionIn / (FRMezcla * FTempIn);
327  FTempTotalIn = FTempIn + pow2(*VelIn) / 2. / FCpMezcla;
328  FPreTotalIn = FPresionIn * pow(FTempTotalIn / FTempIn, FGamma / FGamma1);
329  FGasto1 = FSentidoFlujo * FAreaIn * FDensidadIn * *VelIn;
330 
331  FPreTotalOut = FPreTotalIn * pow(FRelacionCompresion, FSentidoFlujo);
332 
333  if(FFlujoDirecto) {
334  FTempTotalOut = FTempTotalIn * (1 + FCoefPresiones);
335  } else if(!FPrimerCiclo) {
336  FTempTotalOut = FTempTotalInAnt; // andres
337  FTempTotalIn = FTempTotalOut;
338  } else if(FPrimerCiclo) {
339  FTempTotalOut = FTempTotalInAnt;
340  FTempTotalIn = FTempTotalOutAnt;
341  FPrimerCiclo = false;
342  }
343 
344  FVelOutMin = 0.;
345  FVelOutMax = (double) sig * 2 * __cons::ARef * CarOut / (3 - FGamma);
346 
347  if(FVelOutMax < -sqrt(2 * FCpMezcla * FTempTotalOut)) {
348  FVelOutMax = -sqrt(2 * FCpMezcla * FTempTotalOut);
349  }
350  if(FVelOutMax < -sqrt(2 * FRMezcla * FGamma * FTempTotalOut / FGamma2)) {
351  FVelOutMax = -sqrt(2 * FRMezcla * FGamma * FTempTotalOut / FGamma2);
352  }
353  FCuentaVelOut = 0;
354  FErrorVelOut = 1000.;
355 
356  while(fabs(FErrorVelOut) > FErrorCarOut && FCuentaVelOut < 1050.) {
357  *VelOut = (FVelOutMin + FVelOutMax) / 2.;
358  ++FCuentaVelOut;
359  FDensidadOut = fabs(FGasto1 / (FAreaOut * *VelOut));
360  FTempOut = FTempTotalOut - pow2(*VelOut) / 2. / FCpMezcla;
361  FPresionOut = FDensidadOut * FTempOut * FRMezcla;
362  FErrorVelOut = FPreTotalOut - FPresionOut * pow(FTempTotalOut / FTempOut, FGamma / FGamma1);
363  if(FErrorVelOut > 0)
364  FVelOutMax = *VelOut;
365  else
366  FVelOutMin = *VelOut;
367  }
368  *AOut = sqrt(FGamma * FRMezcla * FTempOut);
369  FErrorVelIn = CarOut - (*AOut - (double) sig * *VelOut * FGamma1 / 2) / __cons::ARef;
370 
371  if(FErrorVelIn > 0)
372  FVelInMax = *VelIn;
373  else
374  FVelInMin = *VelIn;
375 
376  ++FCuentaVelIn;
377  }
378 
379  } catch(exception & N) {
380  std::cout << "ERROR: BuscaErrorCaracteristica en el compresor: " << FNumeroCompresor << std::endl;
381  std::cout << "Tipo de error: " << N.what() << std::endl;
382  throw Exception("ERROR: BuscaErrorCaracteristica en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
383  }
384 }
385 
386 // ---------------------------------------------------------------------------
387 // ---------------------------------------------------------------------------
388 
389 void TCompTubos::LeeCompresor(const char *FileWAM, fpos_t &filepos) {
390  int tipo = 0;
391  try {
392 
393  FILE *fich = fopen(FileWAM, "r");
394  fsetpos(fich, &filepos);
395 
396  fscanf(fich, "%d %d ", &FTuboRotor, &FTuboStator);
397  fscanf(fich, "%lf %lf %lf %lf", &FRadioTip, &FRadioHub, &FRadioRodete, &FLongitudCaract);
398 
399  int format = 0;
400  fscanf(fich, "%d ", &format);
401  if(format == 1) {
402  std::cout << "SAE Compressor map format is not compatible with two pipes compressor model" << endl;
403 
404  cout << "Press anykey to exit.";
405  cin.ignore();
406 
407  exit(0);
408  }
409  Mapa2T->PutRadioHub(FRadioHub);
410  Mapa2T->PutRadioRadioTip(FRadioTip);
411  Mapa2T->PutRadioRodete(FRadioRodete);
412 
413  Mapa2T->LeeMapa(fich);
414 
415  fgetpos(fich, &filepos);
416  fclose(fich);
417  } catch(exception & N) {
418  std::cout << "ERROR: LeeCompresor en el compresor: " << FNumeroCompresor << std::endl;
419  std::cout << "Tipo de error: " << N.what() << std::endl;
420  throw Exception("ERROR: LeeCompresor en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
421  }
422 }
423 
424 // ---------------------------------------------------------------------------
425 // ---------------------------------------------------------------------------
426 
427 void TCompTubos::RelacionTubos(TCondicionContorno **BC, int NumeroCC) {
428  try {
429  double Cp = 0.;
430 
431  for(int i = 0; i < 2; i++) {
432  if(FTuboRotor == BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe->getNumeroTubo()) {
433  FTuboRot = BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe;
434  FPosicionRotor = i;
435  } else if(FTuboStator == BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe->getNumeroTubo()) {
436  FTuboEst = BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe;
437  FPosicionEstator = i;
438  }
439  }
440 
441  if(FTuboRot->getNodoIzq() == NumeroCC) {
442  FExtremoRotor = nmLeft;
443  FSignoRotor = -1;
444  FNodoFinRotor = 0;
445  // FIndExtRotor=0;
446  FAreaRotor = __geom::Circle_area(FTuboRot->GetDiametro(0));
447  FIndiceCCRotor = 0;
448 
449  } else if(FTuboRot->getNodoDer() == NumeroCC) {
450  FExtremoRotor = nmRight;
451  FNodoFinRotor = FTuboRot->getNin() - 1;
452  FSignoRotor = 1;
453  // FIndExtRotor=1;
454  FAreaRotor = __geom::Circle_area(FTuboRot->GetDiametro(FTuboRot->getNin() - 1));
455  FIndiceCCRotor = 1;
456  } else {
457  std::cout << "ERROR: Asigncion tubo entrada compresor " << FNumeroCompresor << std::endl;
458  }
459  if(FTuboEst->getNodoIzq() == NumeroCC) {
460  FExtremoStator = nmLeft;
461  FNodoFinStator = 0;
462  FSignoStator = -1;
463  // FIndExtStator=0;
464  FAreaStator = __geom::Circle_area(FTuboEst->GetDiametro(0));
465  FIndiceCCStator = 0;
466  } else if(FTuboEst->getNodoDer() == NumeroCC) {
467  FExtremoStator = nmRight;
468  FNodoFinStator = FTuboEst->getNin() - 1;
469  FSignoStator = 1;
470  // FIndExtStator=1;
471  FAreaStator = __geom::Circle_area(FTuboEst->GetDiametro(FTuboEst->getNin() - 1));
472  FIndiceCCStator = 1;
473  } else {
474  std::cout << "ERROR: Asignacion tubo salida compresor " << FNumeroCompresor << std::endl;
475  }
476  if(FExtremoRotor == nmLeft) {
477  Cp = (FTuboRot->GetGamma(0) * FTuboRot->GetRMezcla(0)) / (FTuboRot->GetGamma(0) - 1);
478  } else
479  Cp = (FTuboRot->GetGamma(FTuboRot->getNin() - 1) * FTuboRot->GetRMezcla(FTuboRot->getNin() - 1)) / (FTuboRot->GetGamma(
480  FTuboRot->getNin() - 1) - 1);
481 
482  FTemperatura10 = __units::degCToK(FTuboRot->getTemperaturaInicial()) + pow2(FTuboRot->getVelocidadMedia()) / 2. / Cp;
483 
484  // Inicializacion del transporte de especies quimicas.
485  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
486  for(int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
487  FFraccionMasicaEspecie[i] = FTuboRot->GetFraccionMasicaInicial(i);
488  }
489 
490  } catch(exception & N) {
491  std::cout << "ERROR: TCompTubos::RelacionTubos en el compresor: " << FNumeroCompresor << std::endl;
492  std::cout << "Tipo de error: " << N.what() << std::endl;
493  throw Exception("ERROR: LeeCompresor en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
494  }
495 }
496 
497 // ---------------------------------------------------------------------------
498 // ---------------------------------------------------------------------------
499 
500 void TCompTubos::MetodoNewton2D(double *a1, double *a2, double *u1, double *u2, double aa1, double aa2, double cc1,
501  double cc2, double s1, double s2, double k, int sig) {
502  try {
503  double f1 = 0., f2 = 0., J11 = 0., J12 = 0., J22 = 0., J21 = 0., DetJ = 0., da1 = 0., da2 = 0., Error = 0.;
504  double a10 = 0., a20 = 0., v10 = 0., v20 = 0., Lim_u1 = 0., Lim_u2 = 0., aa2new = 0.;
505  bool biseccion = false;
506  int cont = 0;
507 
508  /* Pedro */
509  double a1_local = 0., a2_local = 0., u1_local = 0., u2_local = 0.;
510 
511  a1_local = *a1;
512  a2_local = *a2;
513  u1_local = *u1;
514  u2_local = *u2;
515 
516  a10 = *a1;
517  a20 = *a2;
518 
519  Lim_u1 = (double) sig * 2 * __cons::ARef * cc1 / FGamma2;
520  Lim_u2 = (double) sig * 2 * __cons::ARef * cc2 / FGamma2;
521 
522  // aa2new=sqrt(1+FCoefPresiones)* aa1/pow(FRelacionCompresion,sig*FGamma5);
523 
524  aa2new = aa2;
525  // k_local=k;
526 
527  if(fabs(*u1 / *a1) < 0.7 && fabs(*u2 / *a2) < 0.7) {
528 
529  do {
530  f1 = FGamma * __units::BarToPa(pow(*a1 / pow(aa1, FGamma), 1 / FGamma3) * *u1 * s1 - pow(*a2 / pow(aa2new, FGamma),
531  1 / FGamma3) * *u2 * s2);
532  f2 = (*a2 * *a2 + FGamma3 * *u2 * *u2) / aa2 / aa2 - (*a1 * *a1 + FGamma3 * *u1 * *u1) / aa1 / aa1 * pow(
533  FRelacionCompresion, sig * FGamma1 / FGamma);
534 
535  J11 = FGamma * 2 * s1 / FGamma1 * __units::BarToPa(pow(*a1 / pow(aa1, FGamma),
536  1 / FGamma3)) * (*u1 / *a1 - (double) sig);
537  J12 = -FGamma * 2 * s2 / FGamma1 * __units::BarToPa(pow(*a2 / pow(aa2new, FGamma),
538  1 / FGamma3)) * (*u2 / *a1 + (double) sig);
539  J21 = -2 * pow(FRelacionCompresion, (sig * FGamma1 / FGamma)) / aa1 / aa1 * (*a1 - (double) sig * *u1);
540  J22 = 2 / aa2 / aa2 * (*a2 - (double) sig * *u2);
541 
542  DetJ = J11 * J22 - J12 * J21;
543 
544  if(fabs(DetJ) < 1e-10) {
545  da1 = 0;
546  da2 = *a2 - *a1;
547  } else {
548  da1 = (J22 * f1 - J12 * f2) / DetJ;
549  da2 = (-J21 * f1 + J11 * f2) / DetJ;
550  }
551 
552  *a1 = *a1 - da1;
553  *a2 = *a2 - da2;
554 
555  *u1 = (double) sig * (cc1 * __cons::ARef - *a1) / FGamma3;
556  *u2 = -(double) sig * (cc2 * __cons::ARef - *a2) / FGamma3;
557 
558  if(fabs(*u1) > *a1)
559  *u1 = *a1 * *u1 / fabs(*u1);
560  if(fabs(*u2) > *a2)
561  *u2 = *a2 * *u2 / fabs(*u2);
562 
563  Error = sqrt(pow2(da1 / (*a1 + da1)) + pow2(da2 / (*a2 + da2)));
564 
565  cont++;
566  } while(Error > 1e-12 && cont < 5000);
567 
568  if(cont == 5000 || *u1 * *u2 < 0) {
569  biseccion = true;
570  } else {
571  // FPresionOut=1e5*pow(aa2new/ *a2,-FGamma4);
572  FTempIn = *a1 * *a1 / FGamma / FRMezcla;
573  FTempTotalIn = FTempIn + pow2(*u1) / 2. / FCpMezcla;
574  FPresionIn = __units::BarToPa(pow(aa1 / *a1, -FGamma4));
575  FPreTotalIn = FPresionIn * pow(FTempTotalIn / FTempIn, FGamma / FGamma1);
576  FPreTotalOut = FPreTotalIn * pow(FRelacionCompresion, FSentidoFlujo);
577  FPresionOut = FPreTotalOut / pow(1 + FGamma3 * pow2(*u2 / *a2), FGamma / FGamma1);
578  FTempTotalOut = 1 / FGamma / FRMezcla * (*a2 * *a2 + FGamma3 * *u2 * *u2);
579  }
580  } else {
581  biseccion = true;
582  }
583 
584  if(biseccion) {
585  Biseccion(&u1_local, &u2_local, &a1_local, &a2_local, cc1, aa1, cc2, sig);
586  *u1 = u1_local;
587  *a1 = a1_local;
588  *u2 = u2_local;
589  *a2 = a2_local;
590  }
591 
592  } catch(exception & N) {
593  std::cout << "ERROR: MetodoNewton2D en el compresor: " << FNumeroCompresor << std::endl;
594  std::cout << "Tipo de error: " << N.what() << std::endl;
595  throw Exception("ERROR: LeeCompresor en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
596  }
597 }
598 
599 void TCompTubos::Solver(double *a1, double *a2, double *u1, double *u2, double aa1, double aa2, double cc1, double cc2,
600  double s1, double s2, double k, int sig) {
601  try {
602 
603  stCompSolverA1 FunA1(cc1, cc2, s1, s2, FGamma, aa1, aa2, FRelacionCompresion, FRendimiento, FLongitudCaract,
604  FCoefPresiones, FDeltaTiempo);
605 
606  double X1 = cc1 * __cons::ARef - (FGamma - 1) / 2 * (*u1 + 5);
607  double X2 = cc1 * __cons::ARef - (FGamma - 1) / 2 * (*u1 - 5);
608 
609  double X1lim = cc1 * __cons::ARef * 2 / (FGamma + 1);
610  double X2lim = cc1 * __cons::ARef * 2 / (3 - FGamma);
611 
612  if(zbrac2(FunA1, X1, X2, X1lim, X2lim)) {
613  *a1 = FindRoot(FunA1, X1, X2);
614  *a2 = FunA1.A2;
615  *u1 = FunA1.V1;
616  *u2 = FunA1.V2;
617  FCoefPresiones = FunA1.CPfin;
618  }
619  } catch(exception & N) {
620  std::cout << "ERROR: MetodoNewton2D en el compresor: " << FNumeroCompresor << std::endl;
621  std::cout << "Tipo de error: " << N.what() << std::endl;
622  throw Exception("ERROR: LeeCompresor en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
623  }
624 }
625 
626 // ---------------------------------------------------------------------------
627 // ---------------------------------------------------------------------------
628 
629 void TCompTubos::ExtremoCerrado() {
630  try {
631 
632  FVelocidadIn = 0.;
633  FVelocidadOut = 0.;
634  FGasto1 = 0.;
635  *FCarDIn = *FCarCIn;
636  *FCarDOut = *FCarCOut;
637 
638  } catch(exception & N) {
639  std::cout << "ERROR: TCompTubos::ExtremoCerrado en el compresor: " << FNumeroCompresor << std::endl;
640  std::cout << "Tipo de error: " << N.what() << std::endl;
641  throw Exception("ERROR: TCompTubos::ExtremoCerrado en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
642  }
643 }
644 
645 // ---------------------------------------------------------------------------
646 // ---------------------------------------------------------------------------
647 
648 void TCompTubos::SolveInletBoundary(double &A, double &U, double &Ao, double &Uo) {
649  double K1 = 0., K2 = 0., K3 = 0.;
650 
651  double A_dep = pow(FPresionDep, FGamma5);
652  double A_pipe = *FCarCIn / *FAaIn * pow(FRelacionCompresion, FGamma5);
653 
654  K1 = FGamma2 / FGamma1;
655  K2 = -4 / FGamma1 * *FCarCIn;
656 
657  if(A_dep / A_pipe > 1 + 1e-10) {
658  FFlujoDirecto = false;
659 
660  K3 = 2 * pow2(*FCarCIn) / FGamma1 - pow2(FAsonidoDep / __cons::ARef) / (1 + FCoefPresiones);
661 
662  } else if(A_dep / A_pipe < 1 - 1e-10) {
663  FFlujoDirecto = true;
664 
665  K3 = 2 * pow2(*FCarCIn) / FGamma1 - pow(FPresionDep / FRelacionCompresion, FGamma1 / FGamma) * pow2(*FAaIn);
666  } else {
667  FFlujoDirecto = true;
668 
669  U = 0;
670  A = *FCarCIn;
671  return;
672  }
673 
674  A = QuadraticEqP(K1, K2, K3);
675  U = (*FCarCIn - A) / FGamma3;
676 
677  if(FFlujoDirecto) {
678 
679  stCompOut FCompOut(A, U, FRelacionCompresion, FCoefPresiones, FGamma);
680 
681  double Uo = U;
682 
683  double Lim1 = Uo / fabs(Uo) * 1e-10;
684  double Lim2 = Uo + Uo / fabs(Uo) * 1e-3;
685 
686  if(Lim1 * Lim2 < 0)
687  printf("cosa rara");
688 
689  if(zbrac(FCompOut, Lim1, Lim2)) {
690  Uo = FindRoot(FCompOut, Lim1, Lim2);
691  }
692  Ao = pow(FCompOut.A1outA, FGamma3);
693  }
694 }
695 
696 void TCompTubos::SolveOutletBoundary(double &A, double &U) {
697  double Pref = 1;
698 
699  double Ad = pow(FPresionDep / Pref, 1. / FGamma4);
700  double rel_CCon_Entropia = *FCarCOut / *FAaOut;
701 
702  if(rel_CCon_Entropia / Ad > 1 + 1e-10) { // Flujo entrante al deposito
703  InFlow(Ad, A, U);
704 
705  double xaa2 = pow(*FAaOut, FGamma4);
706  FFlowOut = -FGamma * FAreaOut * __units::BarToPa(pow(A, 1 / FGamma3)) * U / (__cons::ARef * xaa2);
707  // Massflow entrante al deposito negativo
708 
709  } else if(rel_CCon_Entropia / Ad < 1 - 1e-10) {
710  // Flujo saliente del deposito
711  OutFlow(Ad, A, U);
712 
713  FFlowOut = -FAreaOut * FGamma * pow(Ad * __cons::ARef / FAsonidoDep, FGamma4) * U * __units::BarToPa(pow(A,
714  2. / FGamma1)) / __cons::ARef;
715 
716  } else { // Flujo Parado
717  U = 0;
718  A = *FCarCOut;
719  FFlowOut = 0.;
720  }
721 
722 }
723 
724 void TCompTubos::NewPropertiesInTheVolume() {
725  double H = 0.;
726 
727  double MasaIn = FFlowIn * FDeltaTiempo;
728  double MasaOut = FFlowOut * FDeltaTiempo;
729 
730  if(FFlowIn > 0) {
731  H += EntalpiaEntrada(FA1out, FU1out, MasaIn, FAsonidoDep / __cons::ARef, FMasaDep, FGamma);
732  }
733  if(FFlowOut > 0) {
734  H += EntalpiaEntrada(FA2, FU2, MasaOut, FAsonidoDep / __cons::ARef, FMasaDep, FGamma);
735  }
736  double MasaDep0 = FMasaDep;
737  FMasaDep = FMasaDep + MasaIn + MasaOut;
738 
739  double Energia = pow(FMasaDep / MasaDep0 * exp(H), FGamma1);
740  FAsonidoDep *= sqrt(Energia);
741  double A2 = FAsonidoDep * FAsonidoDep;
742  FPresionDep = __units::PaToBar(A2 / FGamma / FVolumen * FMasaDep);
743  FTempDep = __units::KTodegC(A2 / FGamma / FRMezcla);
744 
745 }
746 
747 void TCompTubos::InFlow(double Ad, double &A, double &U) {
748  double AThroat = *FAaOut * Ad;
749  nmCaso Caso = nmFlujoEntranteSaltoSubcritico;
750 
751  A = AThroat;
752  U = -(*FCarCOut - A) / FGamma3;
753 
754  double UThroat = U;
755 
756  if(UThroat > AThroat) {
757 
758  } else {
759 
760  }
761 }
762 
763 void TCompTubos::OutFlow(double Ad, double &A, double &U) {
764  // Critical conditions
765 
766  double U2cr = FAsonidoDep * sqrt(2. / FGamma2) * (sqrt(1 + FGamma1 * FGamma2) - 1) / FGamma1;
767  double A2cr = sqrt(pow2(FAsonidoDep) - FGamma3 * pow2(U2cr));
768 
769  stFSSub FSA2(*FAaOut, Ad, FGamma, 1, *FCarCOut, FAsonidoDep / __cons::ARef);
770 
771  double error = FSA2(A2cr / __cons::ARef);
772 
773  if(error < 0.) { // Salto de presiones supercritico.
774  } else {
775  A = FindRoot(FSA2, A2cr / __cons::ARef, FAsonidoDep / __cons::ARef);
776  U = FSA2.U2;
777  }
778 
779 }
780 
781 double TCompTubos::EntalpiaEntrada(double ASonidoE, double VelocidadE, double MasaE, double ASonidoD, double MasaD,
782  double Gamma) {
783 
784  double xx = 0., yy = 0., ret_val = 0.;
785 
786  if(fabs(MasaE) != 0.) {
787  xx = (pow2(ASonidoE / ASonidoD) - 1.) / __Gamma::G1(Gamma);
788  yy = pow2(VelocidadE / ASonidoD) / 2.;
789  ret_val = Gamma * MasaE * (xx + yy) / MasaD;
790  } else {
791  ret_val = 0.;
792  }
793  return ret_val;
794 
795 }
796 
797 void TCompTubos::AsignPipes(TCondicionContorno **BC, int NumeroCC) {
798 
799  for(int i = 0; i < 2; i++) {
800  if(FTuboRotor == BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe->getNumeroTubo()) {
801  FTuboRot = BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe;
802  FPosicionRotor = i;
803  } else if(FTuboStator == BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe->getNumeroTubo()) {
804  FTuboEst = BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe;
805  FPosicionEstator = i;
806  }
807  }
808 
809  if(FTuboRot->getNodoIzq() == NumeroCC) {
810  FExtremoRotor = nmLeft;
811  FSignoRotor = -1;
812  FNodoFinRotor = 0;
813  // FIndExtRotor=0;
814  FAreaRotor = __geom::Circle_area(FTuboRot->GetDiametro(0));
815  FIndiceCCRotor = 0;
816 
817  } else if(FTuboRot->getNodoDer() == NumeroCC) {
818  FExtremoRotor = nmRight;
819  FNodoFinRotor = FTuboRot->getNin() - 1;
820  FSignoRotor = 1;
821  // FIndExtRotor=1;
822  FAreaRotor = __geom::Circle_area(FTuboRot->GetDiametro(FTuboRot->getNin() - 1));
823  FIndiceCCRotor = 1;
824  } else {
825  std::cout << "ERROR: Asigncion tubo entrada compresor " << FNumeroCompresor << std::endl;
826  }
827  if(FTuboEst->getNodoIzq() == NumeroCC) {
828  FExtremoStator = nmLeft;
829  FNodoFinStator = 0;
830  FSignoStator = -1;
831  // FIndExtStator=0;
832  FAreaStator = __geom::Circle_area(FTuboEst->GetDiametro(0));
833  FIndiceCCStator = 0;
834  } else if(FTuboEst->getNodoDer() == NumeroCC) {
835  FExtremoStator = nmRight;
836  FNodoFinStator = FTuboEst->getNin() - 1;
837  FSignoStator = 1;
838  // FIndExtStator=1;
839  FAreaStator = __geom::Circle_area(FTuboEst->GetDiametro(FTuboEst->getNin() - 1));
840  FIndiceCCStator = 1;
841  } else {
842  std::cout << "ERROR: Asignacion tubo salida compresor " << FNumeroCompresor << std::endl;
843  }
844 
845 }
846 
847 void TCompTubos::Initialize() {
848 
849  double Cp = 0.;
850 
851  if(FExtremoRotor == nmLeft) {
852  Cp = (FTuboRot->GetGamma(0) * FTuboRot->GetRMezcla(0)) / (FTuboRot->GetGamma(0) - 1);
853  } else
854  Cp = (FTuboRot->GetGamma(FTuboRot->getNin() - 1) * FTuboRot->GetRMezcla(FTuboRot->getNin() - 1)) / (FTuboRot->GetGamma(
855  FTuboRot->getNin() - 1) - 1);
856 
857  FTemperatura10 = __units::degCToK(FTuboRot->getTemperaturaInicial()) + pow2(FTuboRot->getVelocidadMedia()) / 2. / Cp;
858 
859  // Inicializacion del transporte de especies quimicas.
860  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
861  for(int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
862  FFraccionMasicaEspecie[i] = FTuboRot->GetFraccionMasicaInicial(i);
863  }
864 
865 }
866 
867 #pragma package(smart_init)
TTubo::getNodoDer
int getNodoDer() const
Gets the right-hand side node.
Definition: TTubo.cpp:5456
TMapaComp2Tub
Definition: TMapaComp2Tub.h:41
TTubo::getTemperaturaInicial
double getTemperaturaInicial() const
Gets the initial temperature.
Definition: TTubo.cpp:5480
TTubo::GetGamma
double GetGamma(int i) const
Gets the specific heat capacities ratio at a given cell.
Definition: TTubo.cpp:5444
FindRoot
double FindRoot(T &func, const double x1, const double x2)
Finds the root of a function.
Definition: Math_wam.h:459
stTuboExtremo
Definition: Globales.h:730
TTubo::GetDiametro
double GetDiametro(int i) const
Gets the cell diameter.
Definition: TTubo.cpp:5436
stCompSolverA1
Definition: NewCompSolver.h:78
TTubo::GetFraccionMasicaInicial
double GetFraccionMasicaInicial(int i) const
Gets the initial mass fraction of species i.
Definition: TTubo.cpp:5440
TCompresor
Definition: TCompresor.h:47
TCondicionContorno
Definition: TCondicionContorno.h:54
stFSSub
Definition: BoundaryFunctions.h:93
Exception
Custom exception class.
Definition: Exception.hpp:39
stCompOut
Definition: TCompTubos.h:44
TTubo.h
pow2
T pow2(T x)
Returns x to the power of 2.
Definition: Math_wam.h:88
TTubo::GetRMezcla
double GetRMezcla(int i) const
Gets the gas constant of the mixture at a given cell.
Definition: TTubo.cpp:5476
TTubo::getVelocidadMedia
double getVelocidadMedia() const
Gets the mean speed.
Definition: TTubo.cpp:5509
TTubo::getNodoIzq
int getNodoIzq() const
Gets the left-hand side node.
Definition: TTubo.cpp:5460