OpenWAM
TCompTubDep.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 #include "TCompTubDep.h"
36 #include "TCCCompresor.h"
37 #include "TDeposito.h"
38 #include "TTubo.h"
39 // ---------------------------------------------------------------------------
40 // ---------------------------------------------------------------------------
41 
42 TCompTubDep::TCompTubDep(int i, nmTipoCalculoEspecies SpeciesModel, int numeroespecies, nmCalculoGamma GammaCalculation,
43  bool ThereIsEGR) :
44  TCompresor(i, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
45 
46  FModeloCompresor = nmCompOriginal;
47  FCorrector = 1;
48  // Mapa = new TMapaComp(FNumeroCompresor);
49 
50  FMass_filt_ant = 0;
51  FMass_ant = 0;
52  RC_ant = 1.;
53  RC_filt_ant = 1.;
54  FDelay = 0.001;
55 
56 }
57 
58 // ---------------------------------------------------------------------------
59 // ---------------------------------------------------------------------------
60 
61 TCompTubDep::~TCompTubDep() {
62 
63 }
64 
65 // ---------------------------------------------------------------------------
66 // ---------------------------------------------------------------------------
67 
68 void TCompTubDep::LeeCompresor(const char *FileWAM, fpos_t &filepos) {
69  try {
70 
71  int format = 0, ac = 0;
72  int InID = 0, OutID = 0, VolID = 0, StaID = 0, RotID = 0;
73 
74  FILE *fich = fopen(FileWAM, "r");
75  fsetpos(fich, &filepos);
76 
77 #ifdef tchtm
78 
79  fscanf(fich, "%d ", &ac);
80  if(ac == 1) {
81  FIsAcoustic = true;
82  fscanf(fich, "%d %d %d %d %d", &InID, &OutID, &VolID, &RotID, &StaID);
83  FAcComp = new TAcousticCompressor(InID, VolID, OutID, RotID, StaID);
84  }
85 
86 #endif
87 
88  fscanf(fich, "%lf", &FDelay);
89 
90  fscanf(fich, "%d ", &format);
91  if(format == 0) {
92  FCompressorMapFormat = nmOldWAMmap;
93  Mapa = new TMapaComp(FNumeroCompresor);
94 
95  } else {
96  FCompressorMapFormat = nmSAMmap;
97  Mapa = new TSAEMap(FNumeroCompresor);
98  }
99 
100  Mapa->LeeMapa(fich);
101 
102  fgetpos(fich, &filepos);
103  fclose(fich);
104  } catch(exception & N) {
105  std::cout << "ERROR: TCompTubDep::LeeCompresor en el compresor: " << FNumeroCompresor << std::endl;
106  std::cout << "Tipo de error: " << N.what() << std::endl;
107  throw Exception("ERROR: LeeCompresor en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
108  }
109 }
110 
111 // ---------------------------------------------------------------------------
112 // ---------------------------------------------------------------------------
113 
114 double TCompTubDep::CalGastoNuevo(double MasaAire) {
115  double ret_val = 0., ac = 0., bc = 0., cc = 0., discr = 0.;
116  try {
117 
118  FRelacionCompresion = Mapa->EvaluaRCHermite(MasaAire);
119  FRendimiento = Mapa->EvaluaRendSplines(MasaAire);
120  if(FRendimiento < 0.01)
121  FRendimiento = 0.01;
122 
123 #ifdef tchtm
124 
125  double Kef = FAcComp->EFCorrector(FCorrector, FRelacionCompresion);
126  FRelacionCompresion = FRelacionCompresion * FCorrector;
127  FRendimiento = FRendimiento * Kef;
128 
129 #endif
130 
131  if(FDeltaTiempo > 0) {
132  RC_filt = ((2 * FDelay - FDeltaTiempo) * RC_filt_ant + FDeltaTiempo * (FRelacionCompresion + RC_ant)) /
133  (2 * FDelay + FDeltaTiempo);
134  FRelacionCompresion = RC_filt;
135 
136  }
137  FPresion20 = FRelacionCompresion * FPresion10;
138  FTemperatura20 = FTemperatura10 + (pow(FRelacionCompresion, FGamma5 * 2.) - 1.) * FTemperatura10 / FRendimiento;
139  FEntropia2 = sqrt(FGamma * FRMezcla * FTemperatura20) / pow(FPresion20 / __cons::PRef, FGamma5) / __cons::ARef;
140 
141  ac = pow2(FEntropia2 * FGamma3 / *FEntro) + FGamma3;
142  bc = -pow2(FEntropia2 / *FEntro) * FGamma1 * *FLanda;
143  cc = pow2(FEntropia2 * *FLanda / *FEntro) - FTemperatura20 * FGamma * FRMezcla / __cons::ARef2;
144  discr = bc * bc - ac * 4 * cc;
145 
146  if(discr < 0.) {
147  FVelocidad2 = -bc / 2. / ac;
148  } else {
149  FVelocidad2 = (-bc - sqrt(discr)) / 2. / ac;
150  }
151 
152  FASonidoSalida = (*FLanda - FGamma3 * FVelocidad2) * FEntropia2 / *FEntro;
153  FDensidad20 = __units::BarToPa(FPresion20) / FRMezcla / FTemperatura20;
154  FTemperatura2 = pow2(FASonidoSalida * __cons::ARef) / FGamma / FRMezcla;
155  FDensidad2 = FDensidad20 * pow(FTemperatura2 / FTemperatura20, 1. / FGamma1);
156  ret_val = -FDensidad2 * FVelocidad2 * FAreaSalComp * __cons::ARef * sqrt(FTemperatura10 / Mapa->getTempRef()) *
157  Mapa->getPresionRef() / __units::BarToPa(FPresion10);
158 
159  return ret_val;
160  } catch(exception & N) {
161  std::cout << "ERROR: CalGastoNuevo en el compresor: " << FNumeroCompresor << std::endl;
162  std::cout << "Tipo de error: " << N.what() << std::endl;
163  throw Exception("ERROR: CalGastoNuevo en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
164  }
165 }
166 
167 // ---------------------------------------------------------------------------
168 // ---------------------------------------------------------------------------
169 
170 double TCompTubDep::RegulaFalsi() {
171  double Masa0, Masa1, MasaX, fMasa0, fMasa1, fMasaX, Masa, MasaXant, GastoNuevo;
172  bool valido;
173  try {
174  valido = false;
175  FCambiaReg = false;
176 
177 #ifdef tchtm
178  FCorrector = FAcComp->CRCorrector() * 0.0001 + FCorrector * 0.9999;
179 #endif
180 
181  Masa0 = 0;
182  Masa = CalGastoNuevo(Masa0);
183  fMasa0 = Masa - Masa0;
184 
185  Masa1 = Mapa->getGastoRelComp1();
186  Masa = CalGastoNuevo(Masa1);
187  fMasa1 = Masa - Masa1;
188 
189  if(fMasa0 * fMasa1 < 0)
190  valido = true;
191 
192  if(!valido) {
193  FCambiaReg = true;
194  if(fMasa0 <= 0) {
195  // FVariacionRegimen=1.05;
196  // std::cout << "WARNING: El compresor: " << FNumeroCompresor <<
197  // " intenta trabajar con flujo inverso" << std::endl;
198  GastoNuevo = 0.;
199  FPresion20 = Mapa->getRelCompBombeo() * FPresion10;
200  FASonidoSalida = *FLanda;
201  FVelocidad2 = 0;
202  FEntropia2 = *FEntro;
203  FRelacionCompresion = Mapa->getRelCompBombeo();
204  if(FDeltaTiempo > 0) {
205  RC_filt = ((2 * FDelay - FDeltaTiempo) * RC_filt_ant + FDeltaTiempo * (FRelacionCompresion + RC_ant)) /
206  (2 * FDelay + FDeltaTiempo);
207  FRelacionCompresion = RC_filt;
208 
209  }
210  FRendimiento = Mapa->EvaluaRendSplines(0);
211  return 0;
212  } else {
213 
214  GastoNuevo = Mapa->getGastoRelComp1() / (sqrt(FTemperatura10 / Mapa->getTempRef()) * Mapa->getPresionRef() /
215  __units::BarToPa(FPresion10));
216  FPresion20 = FPresion10;
217  FTemperatura20 = FTemperatura10;
218 
219  double Adep = sqrt(FGamma * FRMezcla * FTemperatura20) / __cons::ARef;
220  FEntropia2 = Adep / pow(FPresion20 / __cons::PRef, FGamma5);
221 
222  double a = pow2(FGamma3 * FEntropia2 / *FEntro) + FGamma3;
223  double b = __Gamma::G1(FGamma) * *FLanda * pow(FEntropia2 / *FEntro, 2);
224  double c = pow2(FEntropia2 / *FEntro * *FLanda) - pow(Adep, 2);
225 
226  FVelocidad2 = -QuadraticEqP(a, b, c);
227  FASonidoSalida = sqrt(pow2(Adep) - FGamma3 * pow2(FVelocidad2));
228 
229  std::cout << "WARNING: El compresor: " << FNumeroCompresor << " esta trabajando en zona de choque" << std::endl;
230  // FVariacionRegimen=0.95;
231 
232  FDensidad2 = FGamma * __units::BarToPa(pow(FASonidoSalida / FEntropia2, FGamma4) / pow2(FASonidoSalida * __cons::ARef));
233 
234  GastoNuevo = -FVelocidad2 * __cons::ARef * FAreaSalComp * FDensidad2 * (sqrt(FTemperatura10 / Mapa->getTempRef()) /
235  Mapa->getPresionRef() * __units::BarToPa(FPresion10));
236 
237  return GastoNuevo;
238  }
239  // return 0;
240  } else {
241  // Metodo regula-falsi modificado
242  MasaXant = Masa0;
243  int i = 0;
244  do {
245  MasaX = (fMasa0 * Masa1 - fMasa1 * Masa0) / (fMasa0 - fMasa1);
246  Masa = CalGastoNuevo(MasaX);
247  fMasaX = Masa - MasaX;
248  if(fMasa0 * fMasaX < 0) {
249  fMasa1 = fMasaX;
250  Masa1 = MasaX;
251  if(Masa0 == MasaXant) {
252  // Se mantiene 2 veces o mas el extremo izquierdo intervalo
253  fMasa0 = fMasa0 / 2.;
254  }
255  MasaXant = Masa0;
256  } else {
257  fMasa0 = fMasaX;
258  Masa0 = MasaX;
259  if(Masa1 == MasaXant) {
260  // Se mantiene 2 veces o mas el extremo izquierdo intervalo
261  fMasa1 = fMasa1 / 2.;
262  }
263  MasaXant = Masa1;
264  }
265  ++i;
266  } while(fabs(fMasaX) > 1e-6 && i <= 100 && fabs(Masa0 - Masa1) > 1e-6);
267  if(i > 100) {
268  std::cout << "ERROR: The interation method in the compressor does not converge in 100 iterations" << std::endl;
269  throw Exception("ERROR: The interation method in the compressor does not converge in 100 iterations");
270  }
271  GastoNuevo = MasaX;
272  return GastoNuevo;
273  }
274 
275  } catch(exception & N) {
276  std::cout << "ERROR: RegulaFalsi en el compresor: " << FNumeroCompresor << std::endl;
277  std::cout << "Tipo de error: " << N.what() << std::endl;
278  throw Exception("ERROR: RegulaFalsi en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
279  }
280 }
281 
282 // ---------------------------------------------------------------------------
283 // ---------------------------------------------------------------------------
284 
285 void TCompTubDep::CalculaCompresor(double Theta) {
286  double temp = 0., work = 0.;
287  try {
288 
289  temp = __units::BarToPa(pow(*FLanda / *FEntro, FGamma4));
290  //double P20Surge = Mapa->getRelCompBombeo() * FPresion10 * 1e5;
291  double P20Max = Mapa->getMaxCompRatio() * __units::BarToPa(FPresion10);
292 
293 #ifdef tchtm
294  P20Max = P20Max * FCorrector;
295 #endif
296  if(temp >= P20Max) {
297  *FBeta = *FLanda;
298  work = 0.;
299  FGastoCorregido = 0.;
300  FRelacionCompresion = Mapa->getRelCompBombeo();
301  if(FDeltaTiempo > 0) {
302  RC_filt = ((2 * FDelay - FDeltaTiempo) * RC_filt_ant + FDeltaTiempo * (FRelacionCompresion + RC_ant)) /
303  (2 * FDelay + FDeltaTiempo);
304  FRelacionCompresion = RC_filt;
305  }
306  FRendimiento = Mapa->EvaluaRendSplines(0);
307  FCambiaReg = false;
308  FASonidoSalida = *FLanda;
309  FEntropia2 = *FEntro;
310  //FGasto1 = NewDampedSolution(0);
311  FGasto1 = 0;
312  } else {
313  FGastoCorregido = RegulaFalsi();
314 
315  *FEntro = FEntropia2;
316  *FLanda = FASonidoSalida + FGamma3 * FVelocidad2;
317  *FBeta = FASonidoSalida - FGamma3 * FVelocidad2;
318  work = (FTemperatura20 - FTemperatura10) * FRMezcla * FGamma4 / 2.;
319 
320  }
321  RC_ant = FRelacionCompresion;
322  RC_filt_ant = RC_filt;
323 
324  FGasto1 = FGastoCorregido * __units::BarToPa(FPresion10) / Mapa->getPresionRef() / sqrt(
325  FTemperatura10 / Mapa->getTempRef());
326  FGasto0 = FGasto1;
327  FGasto0Correg = FGastoCorregido;
328  FTrabajo = work * FGasto1 * FDeltaTiempo;
329  FGastoCompresor = FGasto1;
330 
331  if(FDeltaTiempo != 0) {
332  FPotencia = FTrabajo / FDeltaTiempo;
333  } else
334  FPotencia = 0.;
335  FTrabajoPaso += FTrabajo;
336  FDeltaTPaso += FDeltaTiempo;
337  // }
338 
339  } catch(exception & N) {
340  std::cout << "ERROR: TCompTubDep::CalculaCompresor en el compresor: " << FNumeroCompresor << std::endl;
341  std::cout << "Tipo de error: " << N.what() << std::endl;
342  throw Exception("ERROR: TCompTubDep::CalculaCompresor en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
343  }
344 }
345 
346 // ---------------------------------------------------------------------------
347 // ---------------------------------------------------------------------------
348 
349 void TCompTubDep::CondicionCompresor(double Theta, stTuboExtremo *TuboExtremo, double AcumulatedTime,
350  int TuboCalculado) {
351  try {
352  double FraccionMasicaAcum = 0.;
353 
354  if(FExtremoSalida == nmRight) {
355  FLanda = &(TuboExtremo[0].Landa);
356  FBeta = &(TuboExtremo[0].Beta);
357  FEntro = &(TuboExtremo[0].Entropia);
358  } else {
359  FBeta = &(TuboExtremo[0].Landa);
360  FLanda = &(TuboExtremo[0].Beta);
361  FEntro = &(TuboExtremo[0].Entropia);
362  }
363  FDeltaTiempo = AcumulatedTime - FTiempo0;
364  FTiempo0 = AcumulatedTime;
365 
366  FCambiaReg = true;
367 
368  CalculaCompresor(Theta);
369 
370  // Transporte de especies quimicas.
371  if(FGasto1 > 0.) {
372  switch(FEntradaCompresor) {
373  case nmPipe:
374  // Transporte de especies quimicas.
375  for(int i = 0; i < FNumeroEspecies - 2; i++) {
376  FFraccionMasicaEspecie[i] = FTuboRotor->GetFraccionMasicaCC(FIndiceCC, i);
377  FraccionMasicaAcum += FFraccionMasicaEspecie[i];
378  }
379  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
380  if(FHayEGR)
381  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboRotor->GetFraccionMasicaCC(FIndiceCC, FNumeroEspecies - 1); // EGR
382  break;
383  case nmPlenum:
384  // Transporte de especies quimicas.
385  for(int i = 0; i < FNumeroEspecies - 2; i++) {
386  FFraccionMasicaEspecie[i] = FDeposito->GetFraccionMasicaEspecie(i);
387  FraccionMasicaAcum += FFraccionMasicaEspecie[i];
388  }
389  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
390  if(FHayEGR)
391  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FDeposito->GetFraccionMasicaEspecie(FNumeroEspecies - 1); // EGR
392  break;
393 
394  // Cuando la entrada sea la atmasfera, la composicion se mantendra siempre constante.
395  }
396  }
397 
398  FRegimenCorregido = Mapa->getRegimenCorregido();
399 
400  } catch(exception & N) {
401  std::cout << "ERROR: CondicionCompresor en el compresor: " << FNumeroCompresor << std::endl;
402  std::cout << "Tipo de error: " << N.what() << std::endl;
403  throw Exception("ERROR: CondicionCompresor en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
404  }
405 }
406 
407 // ---------------------------------------------------------------------------
408 // ---------------------------------------------------------------------------
409 
410 void TCompTubDep::BusquedaEntradaSalida(nmCompressorInlet EntradaCompresor, double AmbientTemperature, int numeroCC,
411  TCondicionContorno **BC, double *AtmosphericComposition) {
412  try {
413  double Cp = 0.;
414 
415  FEntradaCompresor = EntradaCompresor;
416  if(BC[numeroCC - 1]->GetTuboExtremo(0).TipoExtremo == nmLeft) {
417  FAreaSalComp = __geom::Circle_area(BC[numeroCC - 1]->GetTuboExtremo(0).Pipe->GetDiametro(0));
418  FExtremoSalida = nmLeft;
419  FNodoFinTuboSalida = 0;
420  FIndiceCC = 0;
421  } else if(BC[numeroCC - 1]->GetTuboExtremo(0).TipoExtremo == nmRight) {
422  FAreaSalComp = __geom::Circle_area(BC[numeroCC - 1]->GetTuboExtremo(0).Pipe->GetDiametro(
423  BC[numeroCC - 1]->GetTuboExtremo(0).Pipe->getNin() - 1));
424  FExtremoSalida = nmRight;
425  FNodoFinTuboSalida = BC[numeroCC - 1]->GetTuboExtremo(0).Pipe->getNin() - 1;
426  FIndiceCC = 1;
427  }
428 
429  switch(FEntradaCompresor) {
430 
431  case nmAtmosphere:
432  FTemperatura10 = __units::degCToK(AmbientTemperature);
433  // Inicializacion del transporte de especies quimicas.
434  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
435  for(int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
436  FFraccionMasicaEspecie[i] = AtmosphericComposition[i];
437  }
438  break;
439  case nmPipe:
440  FTuboRotor = dynamic_cast<TCCCompresor*>(BC[numeroCC - 1])->getTuboRotor();
441  FExtremoTuboRotor = dynamic_cast<TCCCompresor*>(BC[numeroCC - 1])->getExtremoTuboRotor();
442  break;
443  case nmPlenum:
444  FDeposito = dynamic_cast<TCCCompresor*>(BC[numeroCC - 1])->getPlenum();
445  FDeposito->AsignaCompresor(this, -1); /* Entrada */
446  }
447 
448  } catch(exception & N) {
449  std::cout << "ERROR: BusquedaEntradaSalida en el compresor: " << FNumeroCompresor << std::endl;
450  std::cout << "Tipo de error: " << N.what() << std::endl;
451  throw Exception("ERROR: BusquedaEntradaSalida en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
452  }
453 }
454 
455 void TCompTubDep::Initialize() {
456 
457  double Cp = 0.;
458 
459  switch(FEntradaCompresor) {
460  case nmAtmosphere:
461  if(FCalculoEspecies == nmCalculoCompleto) {
462 
463  FRAtm = CalculoCompletoRMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2], 0,
464  FCalculoGamma, nmMEP);
465  FCpAtm = CalculoCompletoCpMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2], 0,
466  FTemperatura10, FCalculoGamma, nmMEP);
467  FGammaAtm = CalculoCompletoGamma(FRAtm, FCpAtm, FCalculoGamma);
468 
469  } else if(FCalculoEspecies == nmCalculoSimple) {
470 
471  FRAtm = CalculoSimpleRMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FCalculoGamma, nmMEP);
472  FCvAtm = CalculoSimpleCvMezcla(FTemperatura10, FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FCalculoGamma,
473  nmMEP);
474  FGammaAtm = CalculoSimpleGamma(FRAtm, FCvAtm, FCalculoGamma);
475  }
476 
477  break;
478 
479  case nmPipe:
480  if(FExtremoTuboRotor == nmLeft) {
481  FNodoFinEntrada = 0;
482  Cp = (FTuboRotor->GetGamma(0) * FTuboRotor->GetRMezcla(0)) / (FTuboRotor->GetGamma(0) - 1);
483  } else {
484  Cp = (FTuboRotor->GetGamma(FTuboRotor->getNin() - 1) * FTuboRotor->GetRMezcla(FTuboRotor->getNin() - 1)) /
485  (FTuboRotor->GetGamma(FTuboRotor->getNin() - 1) - 1);
486  FNodoFinEntrada = FTuboRotor->getNin() - 1;
487  }
488  FTemperatura10 = __units::degCToK(FTuboRotor->getTemperaturaInicial()) + pow(FTuboRotor->getVelocidadMedia(),
489  2.) / 2. / Cp;
490 
491  // Inicializacion del transporte de especies quimicas.
492  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
493  for(int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
494  FFraccionMasicaEspecie[i] = FTuboRotor->GetFraccionMasicaInicial(i);
495  }
496 
497  break;
498  case nmPlenum:
499  FTemperatura10 = __units::degCToK(FDeposito->getTemperature());
500 
501  // Inicializacion del transporte de especies quimicas.
502  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
503  for(int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
504  FFraccionMasicaEspecie[i] = FDeposito->GetFraccionMasicaEspecie(i);
505  }
506  break;
507  }
508 
509 }
510 
511 // ---------------------------------------------------------------------------
512 // ---------------------------------------------------------------------------
513 
514 void TCompTubDep::DatosEntradaCompresor(double AmbientTemperature, double AmbientPressure, TCondicionContorno *BC) {
515  try {
516  double pentcomp = 0., tentcomp = 0., ventcomp = 0.;
517  double RMezclaEnt = 0., RMezclaSal = 0., GammaEnt = 0., GammaSal = 0., Cp = 0.;
518 
519  switch(FEntradaCompresor) {
520  case nmAtmosphere:
521  // Calculo de Gamma y R en el compresor. Media de valores a la entrada y la salida.
522  RMezclaSal = BC->GetTuboExtremo(0).Pipe->GetRMezcla(FNodoFinTuboSalida);
523  GammaSal = BC->GetTuboExtremo(0).Pipe->GetGamma(FNodoFinTuboSalida);
524 
525  FGamma = (GammaSal + FGammaAtm) / 2.;
526  FRMezcla = (RMezclaSal + FRAtm) / 2.;
527  FGamma1 = __Gamma::G1(FGamma);
528  FGamma3 = __Gamma::G3(FGamma);
529  FGamma4 = __Gamma::G4(FGamma);
530  FGamma5 = __Gamma::G5(FGamma);
531 
532  FPresion10 = AmbientPressure;
533  FTemperatura10 = __units::degCToK(AmbientTemperature);
534  break;
535  case nmPlenum:
536  // Calculo de Gamma y R en el compresor. Media de valores a la entrada y la salida.
537  RMezclaEnt = FDeposito->getR();
538  GammaEnt = FDeposito->getGamma();
539 
540  RMezclaSal = BC->GetTuboExtremo(0).Pipe->GetRMezcla(FNodoFinTuboSalida);
541  GammaSal = BC->GetTuboExtremo(0).Pipe->GetGamma(FNodoFinTuboSalida);
542  if(GammaSal < 1.3) {
543  GammaSal = 1.4;
544  }
545  FGamma = (GammaSal + GammaEnt) / 2.;
546  FRMezcla = (RMezclaSal + RMezclaEnt) / 2.;
547  FGamma1 = __Gamma::G1(FGamma);
548  FGamma3 = __Gamma::G3(FGamma);
549  FGamma4 = __Gamma::G4(FGamma);
550  FGamma5 = __Gamma::G5(FGamma);
551 
552  FPresion10 = FDeposito->getPressure();
553  FTemperatura10 = __units::degCToK(FDeposito->getTemperature());
554 
555  break;
556  case nmPipe:
557 
558  // Calculo de Gamma y R en el compresor. Media de valores a la entrada y la salida.
559  RMezclaEnt = FTuboRotor->GetRMezcla(FNodoFinEntrada);
560  GammaEnt = FTuboRotor->GetGamma(FNodoFinEntrada);
561 
562  RMezclaSal = BC->GetTuboExtremo(0).Pipe->GetRMezcla(FNodoFinTuboSalida);
563  GammaSal = BC->GetTuboExtremo(0).Pipe->GetGamma(FNodoFinTuboSalida);
564 
565  FGamma = (GammaSal + GammaEnt) / 2.;
566  FRMezcla = (RMezclaSal + RMezclaEnt) / 2.;
567  FGamma1 = __Gamma::G1(FGamma);
568  FGamma3 = __Gamma::G3(FGamma);
569  FGamma4 = __Gamma::G4(FGamma);
570  FGamma5 = __Gamma::G5(FGamma);
571 
572  if(FExtremoTuboRotor == nmLeft) {
573  pentcomp = __units::BarToPa(FTuboRotor->GetPresion(0));
574  ventcomp = FTuboRotor->GetVelocidad(0) * __cons::ARef;
575  tentcomp = pow2(FTuboRotor->GetAsonido(0) * __cons::ARef) / FGamma / FRMezcla;
576  Cp = (FTuboRotor->GetGamma(0) * FTuboRotor->GetRMezcla(0)) / (FTuboRotor->GetGamma(0) - 1);
577  } else {
578  pentcomp = __units::BarToPa(FTuboRotor->GetPresion(FTuboRotor->getNin() - 1));
579  ventcomp = FTuboRotor->GetVelocidad(FTuboRotor->getNin() - 1) * __cons::ARef;
580  tentcomp = pow2(FTuboRotor->GetAsonido(FTuboRotor->getNin() - 1) * __cons::ARef) / FGamma / FRMezcla;
581  Cp = (FTuboRotor->GetGamma(FTuboRotor->getNin() - 1) * FTuboRotor->GetRMezcla(FTuboRotor->getNin() - 1)) /
582  (FTuboRotor->GetGamma(FTuboRotor->getNin() - 1) - 1);
583  }
584  FTemperatura10 = tentcomp + pow2(ventcomp) / 2. / Cp;
585  FPresion10 = pentcomp * pow((FTemperatura10 / tentcomp), (FGamma / FGamma1));
586  break;
587  }
588 
589  } catch(exception & N) {
590  std::cout << "ERROR: DatosEntradaCompresor en el compresor: " << FNumeroCompresor << std::endl;
591  std::cout << "Tipo de error: " << N.what() << std::endl;
592  throw Exception("ERROR: DatosEntradaCompresor en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
593  }
594 }
595 
596 double TCompTubDep::NewDampedSolution(double Mass) {
597 
598  double A2, U2, AA2, A20, AA2Old, A2Old, delta = 1, delta0 = 1, deltaA = 1.;
599  double A10 = 0., AA1 = 0.;
600  double Mass_filt = 0.;
601  bool conv;
602 
603  // FDelay = 0.1;
604 
605  if(FDeltaTiempo > 0) {
606  Mass_filt = ((2 * FDelay - FDeltaTiempo) * FMass_filt_ant + FDeltaTiempo * (Mass + FMass_ant)) /
607  (2 * FDelay + FDeltaTiempo);
608  FMass_filt_ant = Mass_filt;
609  FMass_ant = Mass;
610  }
611 
612  // Mass_filt = 0.5 * Mass + 0.5 * FMass_ant;
613  // FMass_ant = Mass;
614 
615  // Mass_filt = Mass;
616 
617  FRendimiento = Mapa->EvaluaRendSplines(Mass_filt) * FAcComp->EFCorrector(FCorrector, FRelacionCompresion);
618 
619  A10 = sqrt(FGamma * FRMezcla * FTemperatura10) / __cons::ARef;
620  AA1 = pow(FPresion10, -FGamma5) * A10;
621 
622  A2 = FASonidoSalida;
623  A2Old = A2;
624  AA2 = FEntropia2;
625  AA2Old = AA2;
626  A20 = sqrt(pow2(A2) + FGamma3 * pow2(FVelocidad2));
627 
628  // double m2 = FDensidad2 * FVelocidad2 * FAreaSalComp * __cons::ARef;
629  int step = 0;
630 
631  while((delta > 1e-5 || deltaA > 1e-5) && step < 100) {
632 
633  U2 = Mass_filt / (FGamma * __units::BarToPa(pow(A2 / AA2, FGamma4)) * FAreaSalComp / pow2(A2) / __cons::ARef);
634 
635  A2 = (*FLanda + FGamma3 * U2) * AA2 / *FEntro;
636 
637  A20 = sqrt(pow2(A2) + FGamma3 * pow2(U2));
638 
639  if(U2 <= 0)
640  AA2 = *FEntro;
641  else
642  AA2 = A20 * AA1 / A10 / sqrt((pow2(A20) - pow2(A10)) / pow2(A10) * FRendimiento + 1);
643 
644  delta = fabs((AA2 - AA2Old) / AA2Old);
645 
646  deltaA = fabs((A2 - A2Old) / A2Old);
647 
648  AA2Old = AA2;
649  A2Old = A2;
650 
651  if(delta > delta0) {
652  // Mass_filt = (Mass_filt + Mass) / 2;
653  FMass_filt_ant = Mass_filt;
654  }
655 
656  delta0 = delta;
657 
658  step++;
659  }
660 
661  FEntropia2 = AA2;
662  FVelocidad2 = -U2;
663  FASonidoSalida = A2;
664  FTemperatura2 = pow(FASonidoSalida * __cons::ARef, 2.0) / FGamma / FRMezcla;
665  FTemperatura20 = pow(A20 * __cons::ARef, 2.0) / FGamma / FRMezcla;
666  // FPresion2 = pow(A2/AA2, FGamma4) * __cons::PRef;
667  FPresion20 = pow(A20 / AA2, FGamma4) * __cons::PRef;
668  FRelacionCompresion = FPresion20 / FPresion10;
669  FDensidad20 = __units::BarToPa(FPresion20) / FRMezcla / FTemperatura20;
670  FDensidad2 = FGamma * __units::BarToPa(pow(FASonidoSalida / FEntropia2, FGamma4)) / pow2(FASonidoSalida * __cons::ARef);
671 
672  printf("%lf %lf %lf\n", Mass_filt, Mass, -FVelocidad2 * __cons::ARef * FDensidad2 * FAreaSalComp);
673 
674  // if(conv)
675  // printf("yes %lf %lf ", FTemperatura2, FPresion20);
676  // else
677  // printf("no %lf %lf ", FTemperatura2, FPresion20);
678 
679  return Mass_filt;
680 
681 }
682 
683 // ---------------------------------------------------------------------------
684 // ---------------------------------------------------------------------------
685 
686 #pragma package(smart_init)
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
TTubo::GetFraccionMasicaCC
double GetFraccionMasicaCC(int j, int i)
Definition: TTubo.h:953
TTubo::GetPresion
double GetPresion(int i) const
Gets the fluid pressure.
Definition: TTubo.cpp:5468
pow2
T pow2(T x)
Returns x to the power of 2.
Definition: Math_wam.h:88
TTubo::GetVelocidad
double GetVelocidad(int i) const
Gets the fluid speed.
Definition: TTubo.cpp:5505
TTubo::GetFraccionMasicaInicial
double GetFraccionMasicaInicial(int i) const
Gets the initial mass fraction of species i.
Definition: TTubo.cpp:5440
Exception
Custom exception class.
Definition: Exception.hpp:39
TTubo::GetRMezcla
double GetRMezcla(int i) const
Gets the gas constant of the mixture at a given cell.
Definition: TTubo.cpp:5476
TTubo::GetAsonido
double GetAsonido(int i) const
Gets the speed of sound.
Definition: TTubo.cpp:5412
TTubo::getVelocidadMedia
double getVelocidadMedia() const
Gets the mean speed.
Definition: TTubo.cpp:5509
Generated on Tue Jun 22 2021 07:09:40 for OpenWAM by   doxygen 1.8.15