35 #include "TCompTubDep.h"
36 #include "TCCCompresor.h"
37 #include "TDeposito.h"
42 TCompTubDep::TCompTubDep(
int i, nmTipoCalculoEspecies SpeciesModel,
int numeroespecies, nmCalculoGamma GammaCalculation,
44 TCompresor(i, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
46 FModeloCompresor = nmCompOriginal;
61 TCompTubDep::~TCompTubDep() {
68 void TCompTubDep::LeeCompresor(
const char *FileWAM, fpos_t &filepos) {
71 int format = 0, ac = 0;
72 int InID = 0, OutID = 0, VolID = 0, StaID = 0, RotID = 0;
74 FILE *fich = fopen(FileWAM,
"r");
75 fsetpos(fich, &filepos);
79 fscanf(fich,
"%d ", &ac);
82 fscanf(fich,
"%d %d %d %d %d", &InID, &OutID, &VolID, &RotID, &StaID);
88 fscanf(fich,
"%lf", &FDelay);
90 fscanf(fich,
"%d ", &format);
92 FCompressorMapFormat = nmOldWAMmap;
96 FCompressorMapFormat = nmSAMmap;
97 Mapa =
new TSAEMap(FNumeroCompresor);
102 fgetpos(fich, &filepos);
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());
114 double TCompTubDep::CalGastoNuevo(
double MasaAire) {
115 double ret_val = 0., ac = 0., bc = 0., cc = 0., discr = 0.;
118 FRelacionCompresion = Mapa->EvaluaRCHermite(MasaAire);
119 FRendimiento = Mapa->EvaluaRendSplines(MasaAire);
120 if(FRendimiento < 0.01)
125 double Kef = FAcComp->EFCorrector(FCorrector, FRelacionCompresion);
126 FRelacionCompresion = FRelacionCompresion * FCorrector;
127 FRendimiento = FRendimiento * Kef;
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;
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;
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;
147 FVelocidad2 = -bc / 2. / ac;
149 FVelocidad2 = (-bc - sqrt(discr)) / 2. / ac;
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);
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());
170 double TCompTubDep::RegulaFalsi() {
171 double Masa0, Masa1, MasaX, fMasa0, fMasa1, fMasaX, Masa, MasaXant, GastoNuevo;
178 FCorrector = FAcComp->CRCorrector() * 0.0001 + FCorrector * 0.9999;
182 Masa = CalGastoNuevo(Masa0);
183 fMasa0 = Masa - Masa0;
185 Masa1 = Mapa->getGastoRelComp1();
186 Masa = CalGastoNuevo(Masa1);
187 fMasa1 = Masa - Masa1;
189 if(fMasa0 * fMasa1 < 0)
199 FPresion20 = Mapa->getRelCompBombeo() * FPresion10;
200 FASonidoSalida = *FLanda;
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;
210 FRendimiento = Mapa->EvaluaRendSplines(0);
214 GastoNuevo = Mapa->getGastoRelComp1() / (sqrt(FTemperatura10 / Mapa->getTempRef()) * Mapa->getPresionRef() /
215 __units::BarToPa(FPresion10));
216 FPresion20 = FPresion10;
217 FTemperatura20 = FTemperatura10;
219 double Adep = sqrt(FGamma * FRMezcla * FTemperatura20) / __cons::ARef;
220 FEntropia2 = Adep / pow(FPresion20 / __cons::PRef, FGamma5);
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);
226 FVelocidad2 = -QuadraticEqP(a, b, c);
227 FASonidoSalida = sqrt(
pow2(Adep) - FGamma3 *
pow2(FVelocidad2));
229 std::cout <<
"WARNING: El compresor: " << FNumeroCompresor <<
" esta trabajando en zona de choque" << std::endl;
232 FDensidad2 = FGamma * __units::BarToPa(pow(FASonidoSalida / FEntropia2, FGamma4) /
pow2(FASonidoSalida * __cons::ARef));
234 GastoNuevo = -FVelocidad2 * __cons::ARef * FAreaSalComp * FDensidad2 * (sqrt(FTemperatura10 / Mapa->getTempRef()) /
235 Mapa->getPresionRef() * __units::BarToPa(FPresion10));
245 MasaX = (fMasa0 * Masa1 - fMasa1 * Masa0) / (fMasa0 - fMasa1);
246 Masa = CalGastoNuevo(MasaX);
247 fMasaX = Masa - MasaX;
248 if(fMasa0 * fMasaX < 0) {
251 if(Masa0 == MasaXant) {
253 fMasa0 = fMasa0 / 2.;
259 if(Masa1 == MasaXant) {
261 fMasa1 = fMasa1 / 2.;
266 }
while(fabs(fMasaX) > 1e-6 && i <= 100 && fabs(Masa0 - Masa1) > 1e-6);
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");
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());
285 void TCompTubDep::CalculaCompresor(
double Theta) {
286 double temp = 0., work = 0.;
289 temp = __units::BarToPa(pow(*FLanda / *FEntro, FGamma4));
291 double P20Max = Mapa->getMaxCompRatio() * __units::BarToPa(FPresion10);
294 P20Max = P20Max * FCorrector;
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;
306 FRendimiento = Mapa->EvaluaRendSplines(0);
308 FASonidoSalida = *FLanda;
309 FEntropia2 = *FEntro;
313 FGastoCorregido = RegulaFalsi();
315 *FEntro = FEntropia2;
316 *FLanda = FASonidoSalida + FGamma3 * FVelocidad2;
317 *FBeta = FASonidoSalida - FGamma3 * FVelocidad2;
318 work = (FTemperatura20 - FTemperatura10) * FRMezcla * FGamma4 / 2.;
321 RC_ant = FRelacionCompresion;
322 RC_filt_ant = RC_filt;
324 FGasto1 = FGastoCorregido * __units::BarToPa(FPresion10) / Mapa->getPresionRef() / sqrt(
325 FTemperatura10 / Mapa->getTempRef());
327 FGasto0Correg = FGastoCorregido;
328 FTrabajo = work * FGasto1 * FDeltaTiempo;
329 FGastoCompresor = FGasto1;
331 if(FDeltaTiempo != 0) {
332 FPotencia = FTrabajo / FDeltaTiempo;
335 FTrabajoPaso += FTrabajo;
336 FDeltaTPaso += FDeltaTiempo;
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());
349 void TCompTubDep::CondicionCompresor(
double Theta,
stTuboExtremo *TuboExtremo,
double AcumulatedTime,
352 double FraccionMasicaAcum = 0.;
354 if(FExtremoSalida == nmRight) {
355 FLanda = &(TuboExtremo[0].Landa);
356 FBeta = &(TuboExtremo[0].Beta);
357 FEntro = &(TuboExtremo[0].Entropia);
359 FBeta = &(TuboExtremo[0].Landa);
360 FLanda = &(TuboExtremo[0].Beta);
361 FEntro = &(TuboExtremo[0].Entropia);
363 FDeltaTiempo = AcumulatedTime - FTiempo0;
364 FTiempo0 = AcumulatedTime;
368 CalculaCompresor(Theta);
372 switch(FEntradaCompresor) {
375 for(
int i = 0; i < FNumeroEspecies - 2; i++) {
377 FraccionMasicaAcum += FFraccionMasicaEspecie[i];
379 FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
381 FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboRotor->
GetFraccionMasicaCC(FIndiceCC, FNumeroEspecies - 1);
385 for(
int i = 0; i < FNumeroEspecies - 2; i++) {
386 FFraccionMasicaEspecie[i] = FDeposito->GetFraccionMasicaEspecie(i);
387 FraccionMasicaAcum += FFraccionMasicaEspecie[i];
389 FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
391 FFraccionMasicaEspecie[FNumeroEspecies - 1] = FDeposito->GetFraccionMasicaEspecie(FNumeroEspecies - 1);
398 FRegimenCorregido = Mapa->getRegimenCorregido();
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());
410 void TCompTubDep::BusquedaEntradaSalida(nmCompressorInlet EntradaCompresor,
double AmbientTemperature,
int numeroCC,
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;
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;
429 switch(FEntradaCompresor) {
432 FTemperatura10 = __units::degCToK(AmbientTemperature);
434 FFraccionMasicaEspecie =
new double[FNumeroEspecies - FIntEGR];
435 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
436 FFraccionMasicaEspecie[i] = AtmosphericComposition[i];
440 FTuboRotor =
dynamic_cast<TCCCompresor*
>(BC[numeroCC - 1])->getTuboRotor();
441 FExtremoTuboRotor =
dynamic_cast<TCCCompresor*
>(BC[numeroCC - 1])->getExtremoTuboRotor();
444 FDeposito =
dynamic_cast<TCCCompresor*
>(BC[numeroCC - 1])->getPlenum();
445 FDeposito->AsignaCompresor(
this, -1);
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());
455 void TCompTubDep::Initialize() {
459 switch(FEntradaCompresor) {
461 if(FCalculoEspecies == nmCalculoCompleto) {
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);
469 }
else if(FCalculoEspecies == nmCalculoSimple) {
471 FRAtm = CalculoSimpleRMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FCalculoGamma, nmMEP);
472 FCvAtm = CalculoSimpleCvMezcla(FTemperatura10, FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FCalculoGamma,
474 FGammaAtm = CalculoSimpleGamma(FRAtm, FCvAtm, FCalculoGamma);
480 if(FExtremoTuboRotor == nmLeft) {
484 Cp = (FTuboRotor->
GetGamma(FTuboRotor->getNin() - 1) * FTuboRotor->
GetRMezcla(FTuboRotor->getNin() - 1)) /
485 (FTuboRotor->
GetGamma(FTuboRotor->getNin() - 1) - 1);
486 FNodoFinEntrada = FTuboRotor->getNin() - 1;
492 FFraccionMasicaEspecie =
new double[FNumeroEspecies - FIntEGR];
493 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
499 FTemperatura10 = __units::degCToK(FDeposito->getTemperature());
502 FFraccionMasicaEspecie =
new double[FNumeroEspecies - FIntEGR];
503 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
504 FFraccionMasicaEspecie[i] = FDeposito->GetFraccionMasicaEspecie(i);
514 void TCompTubDep::DatosEntradaCompresor(
double AmbientTemperature,
double AmbientPressure,
TCondicionContorno *BC) {
516 double pentcomp = 0., tentcomp = 0., ventcomp = 0.;
517 double RMezclaEnt = 0., RMezclaSal = 0., GammaEnt = 0., GammaSal = 0., Cp = 0.;
519 switch(FEntradaCompresor) {
522 RMezclaSal = BC->GetTuboExtremo(0).Pipe->
GetRMezcla(FNodoFinTuboSalida);
523 GammaSal = BC->GetTuboExtremo(0).Pipe->
GetGamma(FNodoFinTuboSalida);
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);
532 FPresion10 = AmbientPressure;
533 FTemperatura10 = __units::degCToK(AmbientTemperature);
537 RMezclaEnt = FDeposito->getR();
538 GammaEnt = FDeposito->getGamma();
540 RMezclaSal = BC->GetTuboExtremo(0).Pipe->
GetRMezcla(FNodoFinTuboSalida);
541 GammaSal = BC->GetTuboExtremo(0).Pipe->
GetGamma(FNodoFinTuboSalida);
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);
552 FPresion10 = FDeposito->getPressure();
553 FTemperatura10 = __units::degCToK(FDeposito->getTemperature());
559 RMezclaEnt = FTuboRotor->
GetRMezcla(FNodoFinEntrada);
560 GammaEnt = FTuboRotor->
GetGamma(FNodoFinEntrada);
562 RMezclaSal = BC->GetTuboExtremo(0).Pipe->
GetRMezcla(FNodoFinTuboSalida);
563 GammaSal = BC->GetTuboExtremo(0).Pipe->
GetGamma(FNodoFinTuboSalida);
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);
572 if(FExtremoTuboRotor == nmLeft) {
573 pentcomp = __units::BarToPa(FTuboRotor->
GetPresion(0));
575 tentcomp =
pow2(FTuboRotor->
GetAsonido(0) * __cons::ARef) / FGamma / FRMezcla;
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);
584 FTemperatura10 = tentcomp +
pow2(ventcomp) / 2. / Cp;
585 FPresion10 = pentcomp * pow((FTemperatura10 / tentcomp), (FGamma / FGamma1));
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());
596 double TCompTubDep::NewDampedSolution(
double Mass) {
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.;
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;
617 FRendimiento = Mapa->EvaluaRendSplines(Mass_filt) * FAcComp->EFCorrector(FCorrector, FRelacionCompresion);
619 A10 = sqrt(FGamma * FRMezcla * FTemperatura10) / __cons::ARef;
620 AA1 = pow(FPresion10, -FGamma5) * A10;
626 A20 = sqrt(
pow2(A2) + FGamma3 *
pow2(FVelocidad2));
631 while((delta > 1e-5 || deltaA > 1e-5) && step < 100) {
633 U2 = Mass_filt / (FGamma * __units::BarToPa(pow(A2 / AA2, FGamma4)) * FAreaSalComp /
pow2(A2) / __cons::ARef);
635 A2 = (*FLanda + FGamma3 * U2) * AA2 / *FEntro;
637 A20 = sqrt(
pow2(A2) + FGamma3 *
pow2(U2));
642 AA2 = A20 * AA1 / A10 / sqrt((
pow2(A20) -
pow2(A10)) /
pow2(A10) * FRendimiento + 1);
644 delta = fabs((AA2 - AA2Old) / AA2Old);
646 deltaA = fabs((A2 - A2Old) / A2Old);
653 FMass_filt_ant = Mass_filt;
664 FTemperatura2 = pow(FASonidoSalida * __cons::ARef, 2.0) / FGamma / FRMezcla;
665 FTemperatura20 = pow(A20 * __cons::ARef, 2.0) / FGamma / FRMezcla;
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);
672 printf(
"%lf %lf %lf\n", Mass_filt, Mass, -FVelocidad2 * __cons::ARef * FDensidad2 * FAreaSalComp);
686 #pragma package(smart_init)