36 #include "TCompTubos.h"
37 #include "TCondicionContorno.h"
43 TCompTubos::TCompTubos(
int i, nmTipoCalculoEspecies SpeciesModel,
int numeroespecies, nmCalculoGamma GammaCalculation,
45 TCompresor(i, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
46 FModeloCompresor = nmCompPipes;
49 FCoefPresiones = 0.10344;
58 TCompTubos::~TCompTubos() {
65 void TCompTubos::CondicionCompresor(
double Theta,
stTuboExtremo *TuboExtremo,
double TiempoActual,
int TuboCalculado) {
67 double ErrorVelMax = 0., ErrorVelMin = 0., ErrorCheck = 0., RC = 0., RD = 0., CP = 0.;
69 double CoefPresiones_anterior, Rendimiento_ant = 0.;
70 double FraccionMasicaAcum = 0.;
71 double vout = 0., aaout = 0.;
75 CoefPresiones_anterior = FCoefPresiones;
76 Rendimiento_ant = FRendimiento;
78 FGamma = (FTuboRot->
GetGamma(FNodoFinRotor) + FTuboEst->
GetGamma(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);
88 FDeltaTiempo = TiempoActual - FTiempo0;
89 FTiempo0 = TiempoActual;
93 double vR = FSignoRotor * (TuboExtremo[FPosicionRotor].Landa - TuboExtremo[FPosicionRotor].Beta) / FGamma3 *
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));
99 double aS = (TuboExtremo[FPosicionEstator].Landa + TuboExtremo[FPosicionEstator].Beta) / 2. * __cons::ARef;
100 double vS = FSignoStator * (TuboExtremo[FPosicionEstator].Landa - TuboExtremo[FPosicionEstator].Beta) / FGamma3 *
102 double pS = __units::BarToPa(pow(aS / TuboExtremo[FPosicionEstator].Entropia / __cons::ARef, FGamma4));
111 FPresionDep = __units::PaToBar(pS);
113 FTempDep = aS * aS / FGamma / FRMezcla;
114 FVolumen = FLongitudCaract * FAreaStator;
115 FMasaDep = pS / FRMezcla / FTempDep * FVolumen;
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;
122 if(FExtremoRotor == nmLeft) {
123 FCarCIn = &(TuboExtremo[FPosicionRotor].Beta);
124 FAaIn = &(TuboExtremo[FPosicionRotor].Entropia);
125 FCarDIn = &(TuboExtremo[FPosicionRotor].Landa);
127 FCarCIn = &(TuboExtremo[FPosicionRotor].Landa);
128 FAaIn = &(TuboExtremo[FPosicionRotor].Entropia);
129 FCarDIn = &(TuboExtremo[FPosicionRotor].Beta);
131 if(FExtremoStator == nmRight) {
132 FCarCOut = &(TuboExtremo[FPosicionEstator].Landa);
133 FAaOut = &(TuboExtremo[FPosicionEstator].Entropia);
134 FCarDOut = &(TuboExtremo[FPosicionEstator].Beta);
136 FCarCOut = &(TuboExtremo[FPosicionEstator].Beta);
137 FAaOut = &(TuboExtremo[FPosicionEstator].Entropia);
138 FCarDOut = &(TuboExtremo[FPosicionEstator].Landa);
140 FAreaIn = FAreaRotor;
141 FAreaOut = FAreaStator;
142 FPrimerCiclo =
false;
143 FCoefPresiones =
pow2(*FAaOut / *FAaIn) * pow(FRelacionCompresion, FGamma1 / FGamma) - 1;
146 if(TuboCalculado == 10000) {
148 SolveInletBoundary(FA1in, FU1in, FA1out, FU1out);
153 AA1fin = FAsonidoDep / __cons::ARef * pow(FPresionDep, -FGamma5) * pow(FRelacionCompresion,
154 FGamma5) / sqrt(1 + FCoefPresiones);
158 FFlowIn = FGamma * __units::BarToPa(pow(FA1in * __cons::ARef / pow(AA1fin * __cons::ARef, FGamma),
159 1 / FGamma3)) * FU1in * __cons::ARef * FAreaIn;
163 FAsonidoIn = FA1in * __cons::ARef;
164 FVelocidadIn = FU1in * __cons::ARef;
166 SolveOutletBoundary(FA2, FU2);
170 AA2fin = FAsonidoDep / __cons::ARef * pow(FPresionDep, -FGamma5);
175 NewPropertiesInTheVolume();
176 FAsonidoOut = FA2 * __cons::ARef;
177 FVelocidadOut = FU2 * __cons::ARef;
180 *FCarDIn = (FAsonidoIn - FGamma1 / 2. * FVelocidadIn) / __cons::ARef;
181 if(FVelocidadIn < 0) {
184 *FCarDOut = (FAsonidoOut + FGamma1 / 2. * FVelocidadOut) / __cons::ARef;
185 if(FVelocidadOut > 0) {
190 }
else if(TuboCalculado == FPosicionRotor) {
192 SolveInletBoundary(FA1in, FU1in, FA1out, FU1out);
197 AA1fin = FAsonidoDep / __cons::ARef * pow(FPresionDep, -FGamma5) * pow(FRelacionCompresion,
198 FGamma5) / sqrt(1 + FCoefPresiones);
202 FFlowIn = FGamma * __units::BarToPa(pow(FA1in * __cons::ARef / pow(AA1fin * __cons::ARef, FGamma),
203 1 / FGamma3)) * FU1in * __cons::ARef * FAreaIn;
205 NewPropertiesInTheVolume();
207 FAsonidoIn = FA1in * __cons::ARef;
208 FVelocidadIn = FU1in * __cons::ARef;
211 *FCarDIn = (FAsonidoIn - FGamma1 / 2. * FVelocidadIn) / __cons::ARef;
212 if(FVelocidadIn < 0) {
217 SolveOutletBoundary(FA2, FU2);
221 AA2fin = FAsonidoDep / __cons::ARef * pow(FPresionDep, -FGamma5);
226 NewPropertiesInTheVolume();
227 FAsonidoOut = FA2 * __cons::ARef;
228 FVelocidadOut = FU2 * __cons::ARef;
231 *FCarDOut = (FAsonidoOut + FGamma1 / 2. * FVelocidadOut) / __cons::ARef;
232 if(FVelocidadOut > 0) {
238 FGasto1 = FGamma * __units::BarToPa(pow(FAsonidoIn / pow(*FAaIn * __cons::ARef, FGamma),
239 1 / FGamma3)) * FVelocidadIn * FAreaIn;
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);
250 FGastoCorregido = FGasto1 * sqrt(FTempTotalIn / Mapa2T->getTempRef()) / (FPreTotalIn / Mapa2T->getPresionRef());
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());
265 FRelacionCompresion = FRelacionCompresion + (1 / FLongitudCaract) * (RC - FRelacionCompresion) * 0.5 * fabs(
266 FVelocidadIn + FVelocidadOut) * FDeltaTiempo;
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());
284 FCoefPresiones = (pow(FRelacionCompresion, FGamma1 / FGamma) - 1) / FRendimiento;
286 if(FDeltaTiempo > 0) {
288 FTrabajo = (FTempTotalOut - FTempTotalIn) * FRMezcla * FGamma4 / 2. * fabs(FGasto1) * FDeltaTiempo;
289 FPotencia = FTrabajo / FDeltaTiempo;
291 FTrabajoPaso += FTrabajo;
292 FDeltaTPaso += FDeltaTiempo;
294 FGastoCompresor = FGasto1;
295 FTemperatura10 = FTempTotalIn;
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());
309 void TCompTubos::Biseccion(
double *VelIn,
double *VelOut,
double *AIn,
double *AOut,
double CarIn,
double AaIn,
310 double CarOut,
int sig) {
314 FVelInMax = (double) sig * 2 * __cons::ARef * CarIn / FGamma2;
316 FErrorVelIn = 10000.;
320 while(fabs(FErrorVelIn) > FErrorCarIn && FCuentaVelIn < 200) {
321 *VelIn = (FVelInMin + FVelInMax) / 2.;
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;
331 FPreTotalOut = FPreTotalIn * pow(FRelacionCompresion, FSentidoFlujo);
334 FTempTotalOut = FTempTotalIn * (1 + FCoefPresiones);
335 }
else if(!FPrimerCiclo) {
336 FTempTotalOut = FTempTotalInAnt;
337 FTempTotalIn = FTempTotalOut;
338 }
else if(FPrimerCiclo) {
339 FTempTotalOut = FTempTotalInAnt;
340 FTempTotalIn = FTempTotalOutAnt;
341 FPrimerCiclo =
false;
345 FVelOutMax = (double) sig * 2 * __cons::ARef * CarOut / (3 - FGamma);
347 if(FVelOutMax < -sqrt(2 * FCpMezcla * FTempTotalOut)) {
348 FVelOutMax = -sqrt(2 * FCpMezcla * FTempTotalOut);
350 if(FVelOutMax < -sqrt(2 * FRMezcla * FGamma * FTempTotalOut / FGamma2)) {
351 FVelOutMax = -sqrt(2 * FRMezcla * FGamma * FTempTotalOut / FGamma2);
354 FErrorVelOut = 1000.;
356 while(fabs(FErrorVelOut) > FErrorCarOut && FCuentaVelOut < 1050.) {
357 *VelOut = (FVelOutMin + FVelOutMax) / 2.;
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);
364 FVelOutMax = *VelOut;
366 FVelOutMin = *VelOut;
368 *AOut = sqrt(FGamma * FRMezcla * FTempOut);
369 FErrorVelIn = CarOut - (*AOut - (double) sig * *VelOut * FGamma1 / 2) / __cons::ARef;
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());
389 void TCompTubos::LeeCompresor(
const char *FileWAM, fpos_t &filepos) {
393 FILE *fich = fopen(FileWAM,
"r");
394 fsetpos(fich, &filepos);
396 fscanf(fich,
"%d %d ", &FTuboRotor, &FTuboStator);
397 fscanf(fich,
"%lf %lf %lf %lf", &FRadioTip, &FRadioHub, &FRadioRodete, &FLongitudCaract);
400 fscanf(fich,
"%d ", &format);
402 std::cout <<
"SAE Compressor map format is not compatible with two pipes compressor model" << endl;
404 cout <<
"Press anykey to exit.";
409 Mapa2T->PutRadioHub(FRadioHub);
410 Mapa2T->PutRadioRadioTip(FRadioTip);
411 Mapa2T->PutRadioRodete(FRadioRodete);
413 Mapa2T->LeeMapa(fich);
415 fgetpos(fich, &filepos);
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());
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;
435 }
else if(FTuboStator == BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe->getNumeroTubo()) {
436 FTuboEst = BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe;
437 FPosicionEstator = i;
442 FExtremoRotor = nmLeft;
446 FAreaRotor = __geom::Circle_area(FTuboRot->
GetDiametro(0));
449 }
else if(FTuboRot->
getNodoDer() == NumeroCC) {
450 FExtremoRotor = nmRight;
451 FNodoFinRotor = FTuboRot->getNin() - 1;
454 FAreaRotor = __geom::Circle_area(FTuboRot->
GetDiametro(FTuboRot->getNin() - 1));
457 std::cout <<
"ERROR: Asigncion tubo entrada compresor " << FNumeroCompresor << std::endl;
460 FExtremoStator = nmLeft;
464 FAreaStator = __geom::Circle_area(FTuboEst->
GetDiametro(0));
466 }
else if(FTuboEst->
getNodoDer() == NumeroCC) {
467 FExtremoStator = nmRight;
468 FNodoFinStator = FTuboEst->getNin() - 1;
471 FAreaStator = __geom::Circle_area(FTuboEst->
GetDiametro(FTuboEst->getNin() - 1));
474 std::cout <<
"ERROR: Asignacion tubo salida compresor " << FNumeroCompresor << std::endl;
476 if(FExtremoRotor == nmLeft) {
480 FTuboRot->getNin() - 1) - 1);
485 FFraccionMasicaEspecie =
new double[FNumeroEspecies - FIntEGR];
486 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
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());
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) {
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;
509 double a1_local = 0., a2_local = 0., u1_local = 0., u2_local = 0.;
519 Lim_u1 = (double) sig * 2 * __cons::ARef * cc1 / FGamma2;
520 Lim_u2 = (double) sig * 2 * __cons::ARef * cc2 / FGamma2;
527 if(fabs(*u1 / *a1) < 0.7 && fabs(*u2 / *a2) < 0.7) {
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);
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);
542 DetJ = J11 * J22 - J12 * J21;
544 if(fabs(DetJ) < 1e-10) {
548 da1 = (J22 * f1 - J12 * f2) / DetJ;
549 da2 = (-J21 * f1 + J11 * f2) / DetJ;
555 *u1 = (double) sig * (cc1 * __cons::ARef - *a1) / FGamma3;
556 *u2 = -(double) sig * (cc2 * __cons::ARef - *a2) / FGamma3;
559 *u1 = *a1 * *u1 / fabs(*u1);
561 *u2 = *a2 * *u2 / fabs(*u2);
563 Error = sqrt(
pow2(da1 / (*a1 + da1)) +
pow2(da2 / (*a2 + da2)));
566 }
while(Error > 1e-12 && cont < 5000);
568 if(cont == 5000 || *u1 * *u2 < 0) {
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);
585 Biseccion(&u1_local, &u2_local, &a1_local, &a2_local, cc1, aa1, cc2, sig);
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());
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) {
603 stCompSolverA1 FunA1(cc1, cc2, s1, s2, FGamma, aa1, aa2, FRelacionCompresion, FRendimiento, FLongitudCaract,
604 FCoefPresiones, FDeltaTiempo);
606 double X1 = cc1 * __cons::ARef - (FGamma - 1) / 2 * (*u1 + 5);
607 double X2 = cc1 * __cons::ARef - (FGamma - 1) / 2 * (*u1 - 5);
609 double X1lim = cc1 * __cons::ARef * 2 / (FGamma + 1);
610 double X2lim = cc1 * __cons::ARef * 2 / (3 - FGamma);
612 if(zbrac2(FunA1, X1, X2, X1lim, X2lim)) {
617 FCoefPresiones = FunA1.CPfin;
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());
629 void TCompTubos::ExtremoCerrado() {
636 *FCarDOut = *FCarCOut;
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());
648 void TCompTubos::SolveInletBoundary(
double &A,
double &U,
double &Ao,
double &Uo) {
649 double K1 = 0., K2 = 0., K3 = 0.;
651 double A_dep = pow(FPresionDep, FGamma5);
652 double A_pipe = *FCarCIn / *FAaIn * pow(FRelacionCompresion, FGamma5);
654 K1 = FGamma2 / FGamma1;
655 K2 = -4 / FGamma1 * *FCarCIn;
657 if(A_dep / A_pipe > 1 + 1e-10) {
658 FFlujoDirecto =
false;
660 K3 = 2 *
pow2(*FCarCIn) / FGamma1 -
pow2(FAsonidoDep / __cons::ARef) / (1 + FCoefPresiones);
662 }
else if(A_dep / A_pipe < 1 - 1e-10) {
663 FFlujoDirecto =
true;
665 K3 = 2 *
pow2(*FCarCIn) / FGamma1 - pow(FPresionDep / FRelacionCompresion, FGamma1 / FGamma) *
pow2(*FAaIn);
667 FFlujoDirecto =
true;
674 A = QuadraticEqP(K1, K2, K3);
675 U = (*FCarCIn - A) / FGamma3;
679 stCompOut FCompOut(A, U, FRelacionCompresion, FCoefPresiones, FGamma);
683 double Lim1 = Uo / fabs(Uo) * 1e-10;
684 double Lim2 = Uo + Uo / fabs(Uo) * 1e-3;
689 if(zbrac(FCompOut, Lim1, Lim2)) {
690 Uo =
FindRoot(FCompOut, Lim1, Lim2);
692 Ao = pow(FCompOut.A1outA, FGamma3);
696 void TCompTubos::SolveOutletBoundary(
double &A,
double &U) {
699 double Ad = pow(FPresionDep / Pref, 1. / FGamma4);
700 double rel_CCon_Entropia = *FCarCOut / *FAaOut;
702 if(rel_CCon_Entropia / Ad > 1 + 1e-10) {
705 double xaa2 = pow(*FAaOut, FGamma4);
706 FFlowOut = -FGamma * FAreaOut * __units::BarToPa(pow(A, 1 / FGamma3)) * U / (__cons::ARef * xaa2);
709 }
else if(rel_CCon_Entropia / Ad < 1 - 1e-10) {
713 FFlowOut = -FAreaOut * FGamma * pow(Ad * __cons::ARef / FAsonidoDep, FGamma4) * U * __units::BarToPa(pow(A,
714 2. / FGamma1)) / __cons::ARef;
724 void TCompTubos::NewPropertiesInTheVolume() {
727 double MasaIn = FFlowIn * FDeltaTiempo;
728 double MasaOut = FFlowOut * FDeltaTiempo;
731 H += EntalpiaEntrada(FA1out, FU1out, MasaIn, FAsonidoDep / __cons::ARef, FMasaDep, FGamma);
734 H += EntalpiaEntrada(FA2, FU2, MasaOut, FAsonidoDep / __cons::ARef, FMasaDep, FGamma);
736 double MasaDep0 = FMasaDep;
737 FMasaDep = FMasaDep + MasaIn + MasaOut;
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);
747 void TCompTubos::InFlow(
double Ad,
double &A,
double &U) {
748 double AThroat = *FAaOut * Ad;
749 nmCaso Caso = nmFlujoEntranteSaltoSubcritico;
752 U = -(*FCarCOut - A) / FGamma3;
756 if(UThroat > AThroat) {
763 void TCompTubos::OutFlow(
double Ad,
double &A,
double &U) {
766 double U2cr = FAsonidoDep * sqrt(2. / FGamma2) * (sqrt(1 + FGamma1 * FGamma2) - 1) / FGamma1;
767 double A2cr = sqrt(
pow2(FAsonidoDep) - FGamma3 *
pow2(U2cr));
769 stFSSub FSA2(*FAaOut, Ad, FGamma, 1, *FCarCOut, FAsonidoDep / __cons::ARef);
771 double error = FSA2(A2cr / __cons::ARef);
775 A =
FindRoot(FSA2, A2cr / __cons::ARef, FAsonidoDep / __cons::ARef);
781 double TCompTubos::EntalpiaEntrada(
double ASonidoE,
double VelocidadE,
double MasaE,
double ASonidoD,
double MasaD,
784 double xx = 0., yy = 0., ret_val = 0.;
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;
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;
803 }
else if(FTuboStator == BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe->getNumeroTubo()) {
804 FTuboEst = BC[NumeroCC - 1]->GetTuboExtremo(i).Pipe;
805 FPosicionEstator = i;
810 FExtremoRotor = nmLeft;
814 FAreaRotor = __geom::Circle_area(FTuboRot->
GetDiametro(0));
817 }
else if(FTuboRot->
getNodoDer() == NumeroCC) {
818 FExtremoRotor = nmRight;
819 FNodoFinRotor = FTuboRot->getNin() - 1;
822 FAreaRotor = __geom::Circle_area(FTuboRot->
GetDiametro(FTuboRot->getNin() - 1));
825 std::cout <<
"ERROR: Asigncion tubo entrada compresor " << FNumeroCompresor << std::endl;
828 FExtremoStator = nmLeft;
832 FAreaStator = __geom::Circle_area(FTuboEst->
GetDiametro(0));
834 }
else if(FTuboEst->
getNodoDer() == NumeroCC) {
835 FExtremoStator = nmRight;
836 FNodoFinStator = FTuboEst->getNin() - 1;
839 FAreaStator = __geom::Circle_area(FTuboEst->
GetDiametro(FTuboEst->getNin() - 1));
842 std::cout <<
"ERROR: Asignacion tubo salida compresor " << FNumeroCompresor << std::endl;
847 void TCompTubos::Initialize() {
851 if(FExtremoRotor == nmLeft) {
855 FTuboRot->getNin() - 1) - 1);
860 FFraccionMasicaEspecie =
new double[FNumeroEspecies - FIntEGR];
861 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
867 #pragma package(smart_init)