57 #include "TBloqueMotor.h" 63 nmTipoCalculoEspecies SpeciesModel, nmCalculoGamma GammaCalculation,
bool ThereIsEGR) {
66 FAnguloTotalCiclo = Engine[0]->getAngTotalCiclo();
67 }
else if(Engine == NULL) {
68 FAnguloTotalCiclo = 720.;
69 FRegimenFicticio = 720. / 6. / SimulationDuration;
72 #ifdef ParticulateFilter 73 FDPFEntradaTubo = NULL;
74 FDPFSalidaTubo = NULL;
75 FHayDPFNodoIzq =
false;
76 FHayDPFNodoDer =
false;
89 FNumeroEspecies = SpeciesNumber;
90 FCalculoEspecies = SpeciesModel;
91 FCalculoGamma = GammaCalculation;
92 FNumEcuaciones = 3 + (FNumeroEspecies - 1 - FIntEGR);
94 FComposicionInicial = NULL;
95 FFraccionMasicaEspecie = NULL;
96 FFraccionMasicaCC = NULL;
110 FDiametroTubo = NULL;
116 FDerLinArea12 = NULL;
122 FVelocidadMedia = NULL;
123 FAsonidoMedia = NULL;
124 FPresionMedia = NULL;
147 FCoefTurbulencia = NULL;
148 FSUMTPTuboPro = NULL;
151 ResultadosMedios = NULL;
152 ResultInstantaneos = NULL;
155 FNumDistSensores = 0;
156 FIntercooler =
false;
157 FConcentrico =
false;
158 FMod.Modelo = nmLaxWendroff;
159 FMod.SubModelo = nmNinguno;
160 FMod.OpcionSubModelo = nmNinguna;
161 FMod.FormulacionLeyes = nmSinArea;
172 FCourantLocal = NULL;
196 FResistRadInt = NULL;
197 FResistRadExt = NULL;
198 FResistAxiAnt = NULL;
199 FResistAxiPos = NULL;
204 FVelocidadDim = NULL;
216 DestruyeVector(FVelocidadDim);
217 DestruyeVector(FAsonidoDim);
218 DestruyeVector(FTemperature);
219 DestruyeVector(FFlowMass);
221 DestruyeVector(FComposicionInicial);
224 if(FFraccionMasicaEspecie != NULL) {
225 for(
int i = 0; i < FNin; i++)
226 delete[] FFraccionMasicaEspecie[i];
227 delete[] FFraccionMasicaEspecie;
229 if(FFraccionMasicaCC != NULL) {
230 for(
int i = 0; i < 2; i++)
231 delete[] FFraccionMasicaCC[i];
232 delete[] FFraccionMasicaCC;
234 if(FVelocidadCC != NULL)
235 delete[] FVelocidadCC;
236 if(FDensidadCC != NULL)
237 delete[] FDensidadCC;
255 if(FDExtTramo != NULL)
259 if(FDiametroTubo != NULL)
260 delete[] FDiametroTubo;
261 if(FDiametroD12 != NULL)
262 delete[] FDiametroD12;
263 if(FDiametroS12 != NULL)
264 delete[] FDiametroS12;
267 if(FDerLin12 != NULL)
269 if(FDerLinArea != NULL)
270 delete[] FDerLinArea;
271 if(FDerLinArea12 != NULL)
272 delete[] FDerLinArea12;
277 if(FPresion0 != NULL)
279 if(FAsonido0 != NULL)
281 if(FVelocidad0 != NULL)
282 delete[] FVelocidad0;
283 if(FVelocidadMedia != NULL)
284 delete[] FVelocidadMedia;
285 if(FAsonidoMedia != NULL)
286 delete[] FAsonidoMedia;
287 if(FPresionMedia != NULL)
288 delete[] FPresionMedia;
289 if(FPresion1 != NULL)
291 if(FAsonido1 != NULL)
293 if(FVelocidad1 != NULL)
294 delete[] FVelocidad1;
305 for(
int i = 0; i < 3; i++)
310 for(
int i = 0; i < FNumEcuaciones; i++)
315 for(
int i = 0; i < FNumEcuaciones; i++)
319 if(FU0Medio != NULL) {
320 for(
int i = 0; i < FNumEcuaciones; i++)
321 delete[] FU0Medio[i];
325 for(
int i = 0; i < FNumEcuaciones; i++)
330 for(
int i = 0; i < FNumEcuaciones; i++)
335 for(
int i = 0; i < FNumEcuaciones; i++)
340 for(
int i = 0; i < FNumEcuaciones; i++)
345 for(
int i = 0; i < FNumEcuaciones; i++)
350 for(
int i = 0; i < FNumEcuaciones; i++)
355 for(
int i = 0; i < FNumEcuaciones; i++)
360 for(
int i = 0; i < FNumEcuaciones; i++)
364 if(FUfctad != NULL) {
365 for(
int i = 0; i < FNumEcuaciones; i++)
370 for(
int i = 0; i < FNumEcuaciones; i++)
375 for(
int i = 0; i < FNumEcuaciones; i++)
379 if(FDeltaFCTd != NULL) {
380 for(
int i = 0; i < FNumEcuaciones; i++)
381 delete[] FDeltaFCTd[i];
385 for(
int i = 0; i < FNumEcuaciones; i++)
390 for(
int i = 0; i < FNumEcuaciones; i++)
396 if(FCoefTurbulencia != NULL)
397 delete[] FCoefTurbulencia;
398 if(FSUMTPTuboPro != NULL) {
399 for(
int i = 0; i < 2; i++) {
400 for(
int j = 0; j < 3; j++) {
401 delete[] FSUMTPTuboPro[i][j];
403 delete[] FSUMTPTuboPro[i];
405 delete[] FSUMTPTuboPro;
407 if(FTPTubo != NULL) {
408 for(
int i = 0; i < 3; i++)
412 if(FTParedAnt != NULL) {
413 for(
int i = 0; i < 3; i++)
414 delete[] FTParedAnt[i];
418 for(
int i = 0; i < FNumResMedios; i++) {
419 if(ResultadosMedios[i].FraccionSUM != NULL)
420 delete[] ResultadosMedios[i].FraccionSUM;
421 if(ResultadosMedios[i].FraccionMED != NULL)
422 delete[] ResultadosMedios[i].FraccionMED;
424 if(ResultadosMedios != NULL)
425 delete[] ResultadosMedios;
426 for(
int i = 0; i < FNumResInstant; i++) {
427 if(ResultInstantaneos[i].FraccionINS != NULL)
428 delete[] ResultInstantaneos[i].FraccionINS;
430 if(ResultInstantaneos != NULL)
431 delete[] ResultInstantaneos;
433 if(FCourantLocal != NULL)
434 delete[] FCourantLocal;
436 for(
int i = 0; i < FNumEcuaciones; i++) {
437 for(
int k = 0; k < FNumEcuaciones; k++) {
438 if(FTVD.Pmatrix != NULL)
439 delete[] FTVD.Pmatrix[i][k];
440 if(FTVD.Qmatrix != NULL)
441 delete[] FTVD.Qmatrix[i][k];
443 if(FTVD.Pmatrix != NULL)
444 delete[] FTVD.Pmatrix[i];
445 if(FTVD.Qmatrix != NULL)
446 delete[] FTVD.Qmatrix[i];
448 for(
int i = 0; i < FNumEcuaciones; i++) {
449 if(FTVD.Bmas != NULL)
450 delete[] FTVD.Bmas[i];
451 if(FTVD.Bmen != NULL)
452 delete[] FTVD.Bmen[i];
453 if(FTVD.Bvector != NULL)
454 delete[] FTVD.Bvector[i];
455 if(FTVD.gflux != NULL)
456 delete[] FTVD.gflux[i];
457 if(FTVD.Alpha != NULL)
458 delete[] FTVD.Alpha[i];
459 if(FTVD.Beta != NULL)
460 delete[] FTVD.Beta[i];
461 if(FTVD.DeltaU != NULL)
462 delete[] FTVD.DeltaU[i];
463 if(FTVD.DeltaB != NULL)
464 delete[] FTVD.DeltaB[i];
465 if(FTVD.DeltaW != NULL)
466 delete[] FTVD.DeltaW[i];
467 if(FTVD.hLandaD != NULL)
468 delete[] FTVD.hLandaD[i];
469 if(FTVD.LandaD != NULL)
470 delete[] FTVD.LandaD[i];
472 delete[] FTVD.Phi[i];
479 if(FTVD.Pmatrix != NULL)
480 delete[] FTVD.Pmatrix;
481 if(FTVD.Qmatrix != NULL)
482 delete[] FTVD.Qmatrix;
483 if(FTVD.Bmas != NULL)
485 if(FTVD.Bmen != NULL)
487 if(FTVD.Bvector != NULL)
488 delete[] FTVD.Bvector;
489 if(FTVD.gflux != NULL)
491 if(FTVD.Alpha != NULL)
493 if(FTVD.Beta != NULL)
495 if(FTVD.DeltaU != NULL)
496 delete[] FTVD.DeltaU;
497 if(FTVD.DeltaB != NULL)
498 delete[] FTVD.DeltaB;
499 if(FTVD.DeltaW != NULL)
500 delete[] FTVD.DeltaW;
501 if(FTVD.hLandaD != NULL)
502 delete[] FTVD.hLandaD;
503 if(FTVD.LandaD != NULL)
504 delete[] FTVD.LandaD;
515 if(FResistRadInt != NULL)
516 delete[] FResistRadInt;
517 if(FResistRadExt != NULL)
518 delete[] FResistRadExt;
519 if(FResistAxiAnt != NULL)
520 delete[] FResistAxiAnt;
521 if(FResistAxiPos != NULL)
522 delete[] FResistAxiPos;
539 int TipTC = 0, Concentrico;
541 double fracciontotal = 0.;
543 FILE *fich = fopen(FileWAM,
"r");
544 fsetpos(fich, &filepos);
546 fscanf(fich,
"%d %d %d %d", &FNodoIzq, &FNodoDer, &FNTramos, &FNumeroConductos);
547 fscanf(fich,
"%lf ", &FFriccion);
548 fscanf(fich,
"%lf %lf %lf %lf ", &FTIniParedTub, &FTini, &FPini, &FVelMedia);
549 fscanf(fich,
"%d %lf %lf ", &TipTC, &FCoefAjusTC, &FCoefAjusFric);
553 FTipoTransCal = nmTuboAdmision;
556 FTipoTransCal = nmTuboEscape;
559 FTipoTransCal = nmPipaEscape;
562 FTipoTransCal = nmPipaAdmision;
566 FComposicionInicial =
new double[FNumeroEspecies - FIntEGR];
567 for(
int i = 0; i < FNumeroEspecies - 1; i++) {
568 fscanf(fich,
"%lf ", &FComposicionInicial[i]);
569 fracciontotal += FComposicionInicial[i];
572 if(FCalculoEspecies == nmCalculoCompleto) {
573 if(FComposicionInicial[0] > 0.2)
574 FComposicionInicial[FNumeroEspecies - 1] = 0.;
576 FComposicionInicial[FNumeroEspecies - 1] = 1.;
578 if(FComposicionInicial[0] > 0.5)
579 FComposicionInicial[FNumeroEspecies - 1] = 1.;
581 FComposicionInicial[FNumeroEspecies - 1] = 0.;
585 if(fracciontotal > 1 + 1.e-10 && fracciontotal < 1 - 1e-10) {
586 std::cout <<
"ERROR: Total mass fraction must be equal to 1. Check input data for pipe " << FNumeroTubo << std::endl;
590 fscanf(fich,
"%lf %d ", &FMallado, &FTctpt);
594 FTipoCalcTempPared = nmVariableConInerciaTermica;
597 FTipoCalcTempPared = nmVariableSinInerciaTermica;
600 FTipoCalcTempPared = nmTempConstante;
604 fscanf(fich,
"%d ", &metodo[0]);
607 FMod.Modelo = nmLaxWendroff;
608 fscanf(fich,
"%d ", &metodo[1]);
614 FMod.FormulacionLeyes = nmSinArea;
617 FMod.FormulacionLeyes = nmConArea;
624 FMod.FormulacionLeyes = nmConArea;
625 FMod.SubModelo = nmFCT;
626 fscanf(fich,
"%d ", &metodo[2]);
629 FMod.OpcionSubModelo = nmDDNAD;
630 FMod.Difusion = nmDamping;
631 FMod.Antidifusion = nmNaive;
634 FMod.OpcionSubModelo = nmDDPAD;
635 FMod.Difusion = nmDamping;
636 FMod.Antidifusion = nmPhoenical;
639 FMod.OpcionSubModelo = nmDDEAD;
640 FMod.Difusion = nmDamping;
641 FMod.Antidifusion = nmExplicit;
644 FMod.OpcionSubModelo = nmDSNAD;
645 FMod.Difusion = nmSmoothing;
646 FMod.Antidifusion = nmNaive;
649 FMod.OpcionSubModelo = nmDSPAD;
650 FMod.Difusion = nmSmoothing;
651 FMod.Antidifusion = nmPhoenical;
654 FMod.OpcionSubModelo = nmDSEAD;
655 FMod.Difusion = nmSmoothing;
656 FMod.Antidifusion = nmPhoenical;
660 }
else if(metodo[0] == 2) {
662 FMod.FormulacionLeyes = nmConArea;
665 fscanf(fich,
"%lf ", &FCourant);
667 fgetpos(fich, &filepos);
670 }
catch(exception & N) {
671 std::cout <<
"ERROR: TTubo::LeeDatosGeneralesTubo en el tubo: " << FNumeroTubo << std::endl;
672 std::cout <<
"Tipo de error: " << N.what() << std::endl;
683 double EspesorPrin = 0.;
684 int EsPrincipal = 0, refrigerante = 0, EsFluida = 0;
687 FDExtTramo =
new double[FNTramos + 1];
688 FLTramo =
new double[FNTramos + 1];
693 FILE *fich = fopen(FileWAM,
"r");
694 fsetpos(fich, &filepos);
696 switch(tipomallado) {
698 FTipoMallado = nmDistancia;
701 FTipoMallado = nmAngular;
705 if(FTipoMallado == nmAngular && ene < 0.) {
706 std::cout <<
"ERROR: El mallado no puede ser angular al no existir motor.Pipe: " << FNumeroTubo << std::endl;
710 fscanf(fich,
"%lf ", &FDExtTramo[0]);
711 for(
int i = 1; i <= FNTramos; i++) {
712 fscanf(fich,
"%lf %lf ", &FLTramo[i], &FDExtTramo[i]);
714 CalculoPuntosMalla(ene);
716 if(FTipoCalcTempPared != nmTempConstante) {
717 FResistRadInt =
new double[FNin];
718 FResistRadExt =
new double[FNin];
719 FResistAxiAnt =
new double[FNin];
720 FResistAxiPos =
new double[FNin];
721 FCapInt =
new double[FNin];
722 FCapMed =
new double[FNin];
723 FCapExt =
new double[FNin];
726 fscanf(fich,
"%lf %d ", &FDuracionCiclo, &FNumCiclosSinInerciaTermica);
729 if(FTipoTransCal != nmPipaEscape && FTipoTransCal != nmPipaAdmision) {
730 fscanf(fich,
"%lf %lf %d ", &FCoefExt, &FEmisividad, &refrigerante);
731 switch(refrigerante) {
739 if(FTipRefrig == nmAgua)
740 fscanf(fich,
"%lf ", &FTRefrigerante);
744 fscanf(fich,
"%d", &datoWAMer);
747 fscanf(fich,
"%d ", &FNumCapas);
748 FCapa =
new stCapa[FNumCapas];
750 for(
int i = 0; i < FNumCapas; i++) {
751 fscanf(fich,
"%d %d ", &EsPrincipal, &EsFluida);
752 EsPrincipal == 0 ? FCapa[i].EsPrincipal = false : FCapa[i].EsPrincipal =
true;
753 EsFluida == 0 ? FCapa[i].EsFluida = false : FCapa[i].EsFluida ==
true;
754 if(FCapa[i].EsFluida) {
755 fscanf(fich,
"%lf %lf ", &FCapa[i].EmisividadInterior, &FCapa[i].EmisividadExterior);
757 fscanf(fich,
"%lf %lf %lf %lf ", &FCapa[i].Density, &FCapa[i].CalorEspecifico, &FCapa[i].Conductividad,
759 fscanf(fich,
"%lf %lf ", &FCapa[i].EmisividadInterior, &FCapa[i].EmisividadExterior);
761 if(FCapa[i].EsPrincipal) {
762 FEspesorIntPrin = EspesorPrin;
763 FEspesorPrin = FCapa[i].Espesor;
764 FDensidadPrin = FCapa[i].Density;
765 FCalEspPrin = FCapa[i].CalorEspecifico;
766 FConductPrin = FCapa[i].Conductividad;
769 EspesorPrin += FCapa[i].Espesor;
773 FEspesorExtPrin = EspesorPrin;
774 }
else if(FTipoTransCal == nmPipaEscape || FTipoTransCal == nmPipaAdmision) {
775 fscanf(fich,
"%d", &datoWAMer);
778 fgetpos(fich, &filepos);
781 }
catch(exception & N) {
782 std::cout <<
"ERROR: TTubo::LeeDatosGeometricoTubo en el tubo: " << FNumeroTubo << std::endl;
783 std::cout <<
"Tipo de error: " << N.what() << std::endl;
792 void TTubo::CalculoPuntosMalla(
double ene) {
796 double *FLTotalTramo;
799 FLTotalTramo =
new double[FNTramos + 1];
802 FLTotalTramo[0] = FLTramo[0];
804 for(
int i = 1; i <= FNTramos; i++)
805 FLTotalTramo[i] = FLTotalTramo[i - 1] + FLTramo[i];
806 FLongitudTotal = FLTotalTramo[FNTramos];
808 if(FTipoMallado == nmDistancia)
809 xx = FLongitudTotal / FMallado + 0.5;
810 if(FTipoMallado == nmAngular)
811 xx = FLongitudTotal / 650. / FMallado * 6. * ene + 0.5;
812 if(FTipoMallado != nmDistancia && FTipoMallado != nmAngular) {
813 std::cout <<
"WARNING: No se ha definido una forma correcta de mallado en el tubo: " << FNumeroTubo << std::endl;
820 FXref = FLongitudTotal / (double)(FNin - 1);
823 FDiametroTubo =
new double[FNin];
824 FArea =
new double[FNin];
825 for(
int i = 1; i < FNin - 1; i++) {
828 }
while(FLTotalTramo[ii] < i * FXref);
830 double r1 = i * FXref - FLTotalTramo[ii];
831 FDiametroTubo[i] = Interpola(FDExtTramo[ii], FDExtTramo[ii + 1], FLTramo[ii + 1], r1);
832 FArea[i] = __geom::Circle_area(FDiametroTubo[i]);
834 FDiametroTubo[0] = FDExtTramo[0];
835 FArea[0] = __geom::Circle_area(FDiametroTubo[0]);
836 FDiametroTubo[FNin - 1] = FDExtTramo[FNTramos];
837 FArea[FNin - 1] = __geom::Circle_area(FDiametroTubo[FNin - 1]);
839 FDiametroD12 =
new double[FNin];
840 FDiametroS12 =
new double[FNin];
841 FArea12 =
new double[FNin];
842 for(
int i = 0; i < FNin - 1; i++) {
843 FDiametroD12[i] = (FDiametroTubo[i + 1] + FDiametroTubo[i]) / 2.;
844 FDiametroS12[i] = sqrt((
pow2(FDiametroTubo[i + 1]) +
pow2(FDiametroTubo[i])) / 2.);
845 FArea12[i] = (FArea[i + 1] + FArea[i]) / 2.;
848 delete[] FLTotalTramo;
850 }
catch(exception & N) {
851 std::cout <<
"ERROR: TTubo::CalculoPuntosMalla en el tubo: " << FNumeroTubo << std::endl;
852 std::cout <<
"Tipo de error: " << N.what() << std::endl;
866 for(
int i = 0; i < BC[FNodoIzq - 1]->getNumeroTubosCC(); i++) {
867 if(FNumeroTubo == BC[FNodoIzq - 1]->GetTuboExtremo(i).Pipe->getNumeroTubo()) {
872 for(
int i = 0; i < BC[FNodoDer - 1]->getNumeroTubosCC(); i++) {
873 if(FNumeroTubo == BC[FNodoDer - 1]->GetTuboExtremo(i).Pipe->getNumeroTubo()) {
880 catch(exception & N) {
881 std::cout <<
"ERROR: TTubo::ComunicacionTubo_CC en el tubo: " << FNumeroTubo << std::endl;
882 std::cout <<
"Tipo de error: " << N.what() << std::endl;
891 #ifdef ParticulateFilter 896 bool PrimeraVez =
false;
898 if(FTipoCalcTempPared != nmTempConstante) {
899 FTipoCanal =
new int[2];
900 if(CC[FNodoIzq - 1]->getTipoCC() == nmPipeToPlenumConnection) {
901 numDeposito = dynamic_cast<TCCDeposito*>(CC[FNodoIzq - 1])
902 ->getNumeroDeposito();
903 for(
int k = 0; k < Deposito[numDeposito - 1]->getNUniones();
905 if(Deposito[numDeposito - 1]->GetCCDeposito(k)->getUnionDPF
907 FHayDPFNodoIzq =
true;
908 FTipoCanal[0] = Deposito[numDeposito - 1]->GetCCDeposito
909 (k)->GetTuboExtremo(0).TipoCanal;
910 FDPFEntradaTubo = Deposito[numDeposito - 1]
911 ->GetCCDeposito(k)->GetTuboExtremo(0).DPF;
912 if(Deposito[numDeposito - 1]->GetCCDeposito(k)
913 ->GetTuboExtremo(0).TipoExtremo == nmLeft) {
916 FNodoDPFEntrada = FDPFEntradaTubo->GetCanal(0, 0)
924 if(CC[FNodoDer - 1]->getTipoCC() == nmPipeToPlenumConnection) {
925 numDeposito = dynamic_cast<TCCDeposito*>(CC[FNodoDer - 1])
926 ->getNumeroDeposito();
927 for(
int k = 0; k < Deposito[numDeposito - 1]->getNUniones();
929 if(Deposito[numDeposito - 1]->GetCCDeposito(k)->getUnionDPF
931 FHayDPFNodoDer =
true;
932 FTipoCanal[1] = Deposito[numDeposito - 1]->GetCCDeposito
933 (k)->GetTuboExtremo(0).TipoCanal;
934 FDPFSalidaTubo = Deposito[numDeposito - 1]
935 ->GetCCDeposito(k)->GetTuboExtremo(0).DPF;
936 if(Deposito[numDeposito - 1]->GetCCDeposito(k)
937 ->GetTuboExtremo(0).TipoExtremo == nmLeft) {
940 FNodoDPFSalida = FDPFSalidaTubo->GetCanal(0, 0)
950 std::cout <<
"ERROR: TTubo::ComunicacionDPF en el Tubo " <<
951 FNumeroTubo << std::endl;
952 std::cout <<
"Tipo de error: " << N.what() << std::endl;
965 double RMezclaIni = 0., CpMezclaIni = 0., CvMezclaIni = 0., GammaIni = 0.;
967 FPresion0 =
new double[FNin];
968 FAsonido0 =
new double[FNin];
969 FVelocidad0 =
new double[FNin];
970 FPresion1 =
new double[FNin];
971 FAsonido1 =
new double[FNin];
972 FVelocidadMedia =
new double[FNin];
973 FPresionMedia =
new double[FNin];
974 FAsonidoMedia =
new double[FNin];
975 FVelocidad1 =
new double[FNin];
977 FVelocidadDim =
new double[FNin];
978 FAsonidoDim =
new double[FNin];
979 FTemperature =
new double[FNin];
980 FFlowMass =
new double[FNin];
982 FFraccionMasicaEspecie =
new double*[FNin];
983 for(
int i = 0; i < FNin; i++)
984 FFraccionMasicaEspecie[i] =
new double[FNumeroEspecies - FIntEGR];
985 FFraccionMasicaCC =
new double*[2];
986 for(
int i = 0; i < 2; i++)
987 FFraccionMasicaCC[i] =
new double[FNumeroEspecies - FIntEGR];
988 FVelocidadCC =
new double[2];
989 FDensidadCC =
new double[2];
990 FAreaCC =
new double[2];
991 FGamma =
new double[FNin];
992 FGamma1 =
new double[FNin];
993 FGamma3 =
new double[FNin];
994 FGamma4 =
new double[FNin];
995 FGamma5 =
new double[FNin];
996 FGamma6 =
new double[FNin];
997 FRMezcla =
new double[FNin];
999 FUt =
new double*[3];
1000 for(
int i = 0; i < 3; i++)
1001 FUt[i] =
new double[FNin];
1002 FU0 =
new double*[FNumEcuaciones];
1003 for(
int i = 0; i < FNumEcuaciones; i++)
1004 FU0[i] =
new double[FNin];
1005 FU0Sum =
new double*[FNumEcuaciones];
1006 for(
int i = 0; i < FNumEcuaciones; i++)
1007 FU0Sum[i] =
new double[FNin];
1008 for(
int i = 0; i < FNin; i++) {
1009 for(
int k = 0; k < FNumEcuaciones; k++) {
1013 FU0Medio =
new double*[FNumEcuaciones];
1014 for(
int i = 0; i < FNumEcuaciones; i++)
1015 FU0Medio[i] =
new double[FNin];
1016 FU1 =
new double*[FNumEcuaciones];
1017 for(
int i = 0; i < FNumEcuaciones; i++)
1018 FU1[i] =
new double[FNin];
1019 FU12 =
new double*[FNumEcuaciones];
1020 for(
int i = 0; i < FNumEcuaciones; i++)
1021 FU12[i] =
new double[FNin];
1022 FW =
new double*[FNumEcuaciones];
1023 for(
int i = 0; i < FNumEcuaciones; i++)
1024 FW[i] =
new double[FNin];
1025 FV1 =
new double*[FNumEcuaciones];
1026 for(
int i = 0; i < FNumEcuaciones; i++)
1027 FV1[i] =
new double[FNin];
1028 FV2 =
new double*[FNumEcuaciones];
1029 for(
int i = 0; i < FNumEcuaciones; i++)
1030 FV2[i] =
new double[FNin];
1031 if(FMod.SubModelo == nmFCT) {
1032 FUfct0 =
new double*[FNumEcuaciones];
1033 for(
int i = 0; i < FNumEcuaciones; i++)
1034 FUfct0[i] =
new double[FNin + 1];
1035 FUfct1 =
new double*[FNumEcuaciones];
1036 for(
int i = 0; i < FNumEcuaciones; i++)
1037 FUfct1[i] =
new double[FNin + 1];
1038 FUfctd =
new double*[FNumEcuaciones];
1039 for(
int i = 0; i < FNumEcuaciones; i++)
1040 FUfctd[i] =
new double[FNin + 1];
1041 FUfctad =
new double*[FNumEcuaciones];
1042 for(
int i = 0; i < FNumEcuaciones; i++)
1043 FUfctad[i] =
new double[FNin + 1];
1044 Ffl =
new double*[FNumEcuaciones];
1045 for(
int i = 0; i < FNumEcuaciones; i++)
1046 Ffl[i] =
new double[FNin + 1];
1047 FdU =
new double*[FNumEcuaciones];
1048 for(
int i = 0; i < FNumEcuaciones; i++)
1049 FdU[i] =
new double[FNin + 1];
1050 FDeltaFCTd =
new double*[FNumEcuaciones];
1051 for(
int i = 0; i < FNumEcuaciones; i++)
1052 FDeltaFCTd[i] =
new double[FNin + 1];
1053 FflU =
new double*[FNumEcuaciones];
1054 for(
int i = 0; i < FNumEcuaciones; i++)
1055 FflU[i] =
new double[FNin + 1];
1056 FaU =
new double*[FNumEcuaciones];
1057 for(
int i = 0; i < FNumEcuaciones; i++)
1058 FaU[i] =
new double[FNin + 1];
1061 Frho =
new double[FNin];
1062 FRe =
new double[FNin];
1064 if(FMod.Modelo == nmTVD)
1067 FDerLin =
new double[FNin - 1];
1068 FDerLin12 =
new double[FNin - 1];
1070 FDerLinArea =
new double[FNin - 1];
1071 FDerLinArea12 =
new double[FNin - 1];
1077 for(
int i = 0; i < FNin - 1; i++) {
1078 FDerLin[i] = DerLinF(FDiametroTubo[i], FDiametroTubo[i + 1], FXref);
1079 FDerLinArea[i] = DerLinFArea(FArea[i], FArea[i + 1], FXref);
1082 for(
int i = 1; i < FNin - 1; i++) {
1083 FDerLin12[i] = DerLinF(FDiametroD12[i], FDiametroD12[i - 1], FXref);
1084 FDerLinArea12[i] = DerLinFArea(FArea12[i - 1], FArea12[i], FXref);
1088 if(FCalculoEspecies == nmCalculoCompleto) {
1090 RMezclaIni = CalculoCompletoRMezcla(FComposicionInicial[0], FComposicionInicial[1], FComposicionInicial[2], 0,
1091 FCalculoGamma, nmMEP);
1092 CpMezclaIni = CalculoCompletoCpMezcla(FComposicionInicial[0], FComposicionInicial[1], FComposicionInicial[2], 0,
1093 __units::degCToK(FTini), FCalculoGamma, nmMEP);
1094 GammaIni = CalculoCompletoGamma(RMezclaIni, CpMezclaIni, FCalculoGamma);
1096 }
else if(FCalculoEspecies == nmCalculoSimple) {
1098 RMezclaIni = CalculoSimpleRMezcla(FComposicionInicial[0], FComposicionInicial[1], FCalculoGamma, nmMEP);
1099 CvMezclaIni = CalculoSimpleCvMezcla(__units::degCToK(FTini), FComposicionInicial[0], FComposicionInicial[1],
1100 FCalculoGamma, nmMEP);
1101 GammaIni = CalculoSimpleGamma(RMezclaIni, CvMezclaIni, FCalculoGamma);
1105 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
1106 FFraccionMasicaCC[0][j] = FComposicionInicial[j];
1107 FFraccionMasicaCC[1][j] = FComposicionInicial[j];
1110 double viscgas = 1.4615e-6 *
pow150(__units::degCToK(FTini)) / (__units::degCToK(FTini) + 110.4);
1112 for(
int i = 0; i < FNin; i++) {
1114 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
1115 FFraccionMasicaEspecie[i][j] = FComposicionInicial[j];
1117 FRMezcla[i] = RMezclaIni;
1118 FGamma[i] = GammaIni;
1119 FGamma1[i] = __Gamma::G1(FGamma[i]);
1120 FGamma3[i] = __Gamma::G3(FGamma[i]);
1121 FGamma4[i] = __Gamma::G4(FGamma[i]);
1122 FGamma5[i] = __Gamma::G5(FGamma[i]);
1123 FGamma6[i] = __Gamma::G6(FGamma[i]);
1125 FPresion0[i] = FPini;
1126 FTemperature[i] = __units::degCToK(FTini);
1127 FAsonidoDim[i] = sqrt(__units::degCToK(FTini) * FGamma[i] * FRMezcla[i]);
1128 FAsonido0[i] = FAsonidoDim[i] / __cons::ARef;
1129 FVelocidadDim[i] = FVelMedia;
1130 FVelocidad0[i] = FVelMedia / __cons::ARef;
1131 FPresion1[i] = FPini;
1132 FAsonido1[i] = FAsonido0[i];
1133 FVelocidad1[i] = FVelocidad0[i];
1134 FFlowMass[i] = FArea[i] * FVelocidadDim[i] * __units::BarToPa(FPresion0[i]) / FRMezcla[i] / FTemperature[i];
1136 if(FMod.FormulacionLeyes == nmSinArea) {
1137 Transforma1(FVelocidad0[i], FAsonido0[i], FPresion0[i], FU0, FGamma[i], FGamma1[i], FFraccionMasicaEspecie[i], i);
1139 Transforma1(FVelocidad1[i], FAsonido1[i], FPresion1[i], FU1, FGamma[i], FGamma1[i], FFraccionMasicaEspecie[i], i);
1141 Frho[i] = FU0[0][i];
1142 }
else if(FMod.FormulacionLeyes == nmConArea) {
1143 Transforma1Area(FVelocidad0[i], FAsonido0[i], FPresion0[i], FU0, FArea[i], FGamma[i], FGamma1[i],
1144 FFraccionMasicaEspecie[i], i);
1146 Transforma1Area(FVelocidad1[i], FAsonido1[i], FPresion1[i], FU1, FArea[i], FGamma[i], FGamma1[i],
1147 FFraccionMasicaEspecie[i], i);
1149 Frho[i] = FU0[0][i] / FArea[i];
1153 std::cout <<
"ERROR: El tipo de formulacion de las leyes de conservacion no esta bien definido" << std::endl;
1156 FRe[i] = Frho[i] * FVelocidadDim[i] * FDiametroTubo[i] / viscgas;
1159 }
catch(exception & N) {
1160 std::cout <<
"ERROR: TTubo::IniciaVariableFundamentalesTubo en el tubo: " << FNumeroTubo << std::endl;
1161 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1174 for(
int i = 0; i < FNin; i++) {
1175 FTemperature[i] =
pow2(FAsonidoDim[i]) / (FGamma[i] * FRMezcla[i]);
1178 if(FCalculoEspecies == nmCalculoSimple) {
1180 FRMezcla[i] = CalculoSimpleRMezcla(FFraccionMasicaEspecie[i][0], 0, FCalculoGamma, nmMEP);
1181 double CvMezcla = CalculoSimpleCvMezcla(FTemperature[i], FFraccionMasicaEspecie[i][0], 0, FCalculoGamma, nmMEP);
1182 FGammaN = CalculoSimpleGamma(FRMezcla[i], CvMezcla, FCalculoGamma);
1183 if(abs(FGammaN - FGamma[i]) > 0.025) {
1186 FGamma[i] = 0.9995 * FGamma[i] + 0.0005 * FGammaN;
1190 FRMezcla[i] = CalculoCompletoRMezcla(FFraccionMasicaEspecie[i][0], FFraccionMasicaEspecie[i][1],
1191 FFraccionMasicaEspecie[i][2], 0, FCalculoGamma, nmMEP);
1192 double CpMezcla = CalculoCompletoCpMezcla(FFraccionMasicaEspecie[i][0], FFraccionMasicaEspecie[i][1],
1193 FFraccionMasicaEspecie[i][2], 0, FTemperature[i], FCalculoGamma, nmMEP);
1194 FGammaN = CalculoCompletoGamma(FRMezcla[i], CpMezcla, FCalculoGamma);
1195 if(abs(FGammaN - FGamma[i]) > 0.025) {
1198 FGamma[i] = 0.9995 * FGamma[i] + 0.0005 * FGammaN;
1203 FGamma1[i] = __Gamma::G1(FGamma[i]);
1204 FGamma3[i] = __Gamma::G3(FGamma[i]);
1205 FGamma4[i] = __Gamma::G4(FGamma[i]);
1206 FGamma5[i] = __Gamma::G5(FGamma[i]);
1207 FGamma6[i] = __Gamma::G6(FGamma[i]);
1209 Frho[i] = __units::BarToPa(FPresion0[i]) / FRMezcla[i] / FTemperature[i];
1212 }
catch(exception & N) {
1213 std::cout <<
"ERROR: TTubo::CalculoPropiedadesGas en el tubo: " << FNumeroTubo << std::endl;
1214 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1223 void TTubo::Transforma1(
const double& v,
const double& a,
const double& p,
double **U,
const double& Gamma,
1224 const double& Gamma1,
double *Yespecie,
const int& i) {
1228 double V = v * __cons::ARef;
1229 double Ppa = __units::BarToPa(p);
1231 U[0][i] = Gamma * Ppa /
pow2(a * __cons::ARef);
1232 U[1][i] = U[0][i] * V;
1233 U[2][i] = Ppa / Gamma1 + U[1][i] * V / 2.0;
1234 for(
int j = 3; j < 3 + (FNumeroEspecies - 2); j++) {
1235 U[j][i] = U[0][i] * Yespecie[j - 3];
1238 U[3 + (FNumeroEspecies - 2)][i] = U[0][i] * Yespecie[FNumeroEspecies - 1];
1240 }
catch(exception & N) {
1241 std::cout <<
"ERROR: TTubo::Transforma1 en el tubo: " << FNumeroTubo << std::endl;
1242 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1251 void TTubo::Transforma1Area(
const double& v,
const double& a,
const double& p,
double **U,
const double& area,
1252 const double& Gamma,
const double& Gamma1,
double *Yespecie,
const int& i) {
1256 double APpa = area * __units::BarToPa(p);
1257 double V = v * __cons::ARef;
1259 U[0][i] = Gamma * APpa /
pow2(a * __cons::ARef);
1260 U[1][i] = U[0][i] * V;
1261 U[2][i] = APpa / Gamma1 + U[1][i] * V / 2.0;
1262 for(
int j = 3; j < 3 + (FNumeroEspecies - 2); j++) {
1263 U[j][i] = U[0][i] * Yespecie[j - 3];
1266 U[3 + (FNumeroEspecies - 2)][i] = U[0][i] * Yespecie[FNumeroEspecies - 1];
1271 catch(exception & N) {
1272 std::cout <<
"ERROR: TTubo::Transforma1Area en el tubo: " << FNumeroTubo << std::endl;
1273 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1282 void TTubo::Transforma2(
double& v,
double& a,
double& p,
double **U,
const double& Gamma,
const double& Gamma1,
1283 double *Yespecie,
const int& i) {
1287 double fraccionmasicaacum = 0.;
1288 double V = U[1][i] / U[0][i];
1289 v = V / __cons::ARef;
1290 double P = (U[2][i] - U[1][i] * V / 2.0) * Gamma1;
1291 p = __units::PaToBar(P);
1292 a = sqrt(Gamma * P / U[0][i] / __cons::ARef2);
1293 if(v > 1e200 || v < -1e200) {
1294 std::cout <<
"ERROR: Valor de velocidad no valido" << std::endl;
1297 if(p > 1e200 || p < 0) {
1298 std::cout <<
"ERROR: Valor de presion no valido en el tubo " << FNumeroTubo <<
" nodo " << i << std::endl;
1301 if(a > 1e200 || a < 0) {
1302 std::cout <<
"ERROR: Valor de velocidad del sonido no valido" << std::endl;
1303 throw Exception(
"Error velocidad del sonido");
1307 for(
int j = 0; j < FNumeroEspecies - 2; j++) {
1308 Yespecie[j] = U[j + 3][i] / U[0][i];
1309 fraccionmasicaacum += Yespecie[j];
1311 Yespecie[FNumeroEspecies - 2] = 1. - fraccionmasicaacum;
1313 Yespecie[FNumeroEspecies - 1] = U[FNumeroEspecies - 2 + 3][i] / U[0][i];
1315 }
catch(exception & N) {
1316 std::cout <<
"ERROR: TTubo::Transforma2 en el tubo: " << FNumeroTubo << std::endl;
1317 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1326 void TTubo::Transforma2Area(
double& v,
double& a,
double& p,
double **U,
const double& area,
const double& Gamma,
1327 const double& Gamma1,
double *Yespecie,
const int& i) {
1331 double fraccionmasicaacum = 0.;
1333 std::cout <<
"ERROR: Calculation in pipe " << FNumeroTubo <<
" is unstable" << std::endl;
1334 if(FMod.Modelo == nmLaxWendroff)
1335 std::cout <<
" Try to use TVD scheme for this pipe" << std::endl;
1337 std::cout <<
" Check the input data" << std::endl;
1338 throw Exception(
"ERROR: the pipe calculation is unstable");
1341 double V = U[1][i] / U[0][i];
1342 double P = (U[2][i] - U[1][i] * V / 2.0) * Gamma1;
1344 v = V / __cons::ARef;
1345 p = __units::PaToBar(P) / area;
1346 a = sqrt(Gamma * P / U[0][i] / __cons::ARef2);
1348 for(
int j = 0; j < FNumeroEspecies - 2; j++) {
1349 Yespecie[j] = U[j + 3][i] / U[0][i];
1350 fraccionmasicaacum += Yespecie[j];
1352 Yespecie[FNumeroEspecies - 2] = 1. - fraccionmasicaacum;
1354 Yespecie[FNumeroEspecies - 1] = U[FNumeroEspecies - 2 + 3][i] / U[0][i];
1356 }
catch(exception & N) {
1357 std::cout <<
"ERROR: TTubo::Transforma2Area en el tubo: " << FNumeroTubo << std::endl;
1358 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1367 void TTubo::Transforma3Area(
double **Ufct,
double **U,
double Area,
double Gamma,
double Gamma1,
double Gamma6,
int i) {
1372 Ufct[0][i] = U[1][i];
1373 Ufct[1][i] = Gamma * U[2][i] / U[0][i] -
pow2(U[1][i]) * Gamma1 / 2. / U[0][i] / U[0][i];
1374 Ufct[2][i] = ((U[2][i] -
pow2(U[1][i]) / U[0][i] / 2.) * Gamma1) * pow(
1375 (1. + (Gamma1 / 2.) *
pow2(U[1][i]) / U[0][i] / U[0][i] / (Gamma * ((U[2][i] -
pow2(
1376 U[1][i]) / U[0][i] / 2.) * Gamma1) / U[0][i])), Gamma * Gamma6) / Area;
1378 for(
int j = 3; j < FNumEcuaciones; j++) {
1379 Ufct[j][i] = U[j][i] * U[1][i] / U[0][i];
1383 }
catch(exception & N) {
1384 std::cout <<
"ERROR: TTubo::Transforma3Area en el tubo: " << FNumeroTubo << std::endl;
1385 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1394 void TTubo::Transforma4Area(
double **U1,
double **Ufctd,
double Area,
double Gamma,
double Gamma1,
double Gamma3,
1395 double Gamma4,
double Gamma6,
int i) {
1396 double error = 0., fu = 0., dfu = 0., vel = 0., vel1 = 0.;
1404 Y =
new double[FNumeroEspecies - 1];
1406 vel = Ufctd[0][i] / (U1[0][i]);
1410 while(error > 0.00000001) {
1412 fu = vel - Ufctd[0][i] / (Area * Ufctd[2][i] * Gamma4) * pow(2. * Ufctd[1][i],
1413 (Gamma / Gamma1)) * pow((2. * Ufctd[1][i] -
pow2(vel)), -Gamma6);
1415 dfu = 1. - Ufctd[0][i] * vel / (Area * Ufctd[2][i] * Gamma) * pow(2. * Ufctd[1][i],
1416 (Gamma / Gamma1)) * pow((2. * Ufctd[1][i] -
pow2(vel)), -Gamma / Gamma1);
1418 vel1 = vel - fu / dfu;
1419 error = fabs(vel1 - vel);
1424 v = vel / __cons::ARef;
1425 a = sqrt(Gamma1 * (Ufctd[1][i] - (vel * vel) / 2.)) / __cons::ARef;
1426 p = Ufctd[2][i] / pow((1 + Gamma3 *
pow2(v / a)), Gamma / Gamma1) / 1.e5;
1427 for(
int j = 0; j < FNumeroEspecies - 1 - FIntEGR; j++) {
1428 if(Ufctd[0][i] != 0.) {
1429 Y[j] = Ufctd[j + 3][i] / Ufctd[0][i];
1432 Y[j] = FFraccionMasicaEspecie[i][j];
1436 Transforma1Area(v, a, p, U1, Area, Gamma, Gamma1, Y, i);
1441 }
catch(exception & N) {
1442 std::cout <<
"ERROR: TTubo::Transforma4Area en el tubo: " << FNumeroTubo << std::endl;
1443 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1453 double dist1 = 0., dist2 = 0.;
1458 Fhi =
new double[FNin];
1459 Fhe =
new double[FNin];
1460 FVelPro =
new double[FNin];
1461 FCoefTurbulencia =
new double[FNin];
1462 FSUMTPTuboPro =
new double**[2];
1463 for(
int i = 0; i < 2; i++)
1464 FSUMTPTuboPro[i] =
new double*[3];
1465 for(
int i = 0; i < 2; i++) {
1466 for(
int j = 0; j < 3; j++) {
1467 FSUMTPTuboPro[i][j] =
new double[FNin];
1470 FTPTubo =
new double*[3];
1471 for(
int i = 0; i < 3; i++)
1472 FTPTubo[i] =
new double[FNin];
1473 FTParedAnt =
new double*[3];
1474 for(
int i = 0; i < 3; i++)
1475 FTParedAnt[i] =
new double[FNin];
1477 for(
int i = 0; i < FNin; i++) {
1478 dist1 = FXref * (double) i + BC[FNodoIzq - 1]->getPosicionNodo();
1479 dist2 = FXref * (double)(FNin - 1 - i) + BC[FNodoDer - 1]->getPosicionNodo();
1480 if(Minimo(dist1, dist2) >= 10000.)
1481 FCoefTurbulencia[i] = 1.;
1483 FCoefTurbulencia[i] = 1. + 3. * exp(-Minimo(dist1, dist2) / (4 * FDiametroTubo[i]));
1484 for(
int j = 0; j < 3; j++) {
1485 FTPTubo[j][i] = FTIniParedTub;
1486 FTParedAnt[j][i] = FTIniParedTub;
1488 for(
int k = 0; k < 2; k++) {
1489 for(
int j = 0; j < 3; j++) {
1490 FSUMTPTuboPro[k][j][i] = 0.;
1498 if(FTipoCalcTempPared != nmTempConstante) {
1499 if(FTipoTransCal == nmPipaAdmision || FTipoTransCal == nmPipaEscape) {
1500 FTExt = __units::degCToK(Engine[0]->getTempRefrigerante());
1502 if(FTipRefrig == nmAgua)
1503 FTExt = __units::degCToK(FTRefrigerante);
1505 FTExt = __units::degCToK(AmbientTemperature);
1509 }
catch(exception & N) {
1510 std::cout <<
"ERROR: TTubo::IniciaVariablesTransmisionCalor en el tubo: " << FNumeroTubo << std::endl;
1511 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1520 double VTotalMax = 0., VTotalNodo = 0.;
1525 VTotalMax = (FAsonidoDim[0] + fabs(FVelocidadDim[0]));
1526 for(
int i = 0; i < FNin; i++) {
1527 VTotalNodo = (FAsonidoDim[i] + fabs(FVelocidadDim[i]));
1528 if(VTotalNodo > VTotalMax)
1529 VTotalMax = VTotalNodo;
1533 if(FMod.Modelo == nmTVD) {
1540 }
catch(exception & N) {
1541 std::cout <<
"ERROR: TTubo::EstabilidadMetodoCalculo en el tubo: " << FNumeroTubo << std::endl;
1542 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1554 if(FMod.Modelo == nmLaxWendroff && FMod.FormulacionLeyes == nmSinArea)
1556 else if(FMod.Modelo == nmLaxWendroff && FMod.FormulacionLeyes == nmConArea) {
1558 if(FMod.SubModelo == nmFCT) {
1560 FluxCorrectedTransport();
1562 }
else if(FMod.Modelo == nmTVD) {
1565 std::cout <<
"ERROR: Metodo de calculo no implementado" << std::endl;
1569 }
catch(exception & N) {
1570 std::cout <<
"ERROR: TTubo::CalculaVariablesFundamentales en el tubo: " << FNumeroTubo << std::endl;
1571 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1580 void TTubo::LaxWendroff() {
1585 double x1, x2, x3, x4, *hi12, *rho12, *Re12, *TPTubo12, *Gamma12, *Rmezcla12, *Gamma1_12;
1587 hi12 =
new double[FNin - 1];
1588 rho12 =
new double[FNin - 1];
1589 Re12 =
new double[FNin - 1];
1590 TPTubo12 =
new double[FNin - 1];
1591 Gamma12 =
new double[FNin - 1];
1592 Rmezcla12 =
new double[FNin - 1];
1593 Gamma1_12 =
new double[FNin - 1];
1601 CalculaFlujo(FU0, FW, FGamma, FGamma1, Nodos);
1603 CalculaFuente1(FU0, FV1, FGamma, FGamma1, Nodos);
1605 CalculaFuente2(FU0, FV2, FDiametroTubo, Fhi, Frho, FRe, FTPTubo[0], FGamma, FRMezcla, FGamma1, Nodos);
1607 for(
int i = 0; i < FNin - 1; i++) {
1608 for(
int j = 0; j < FNumEcuaciones; j++) {
1609 x1 = FU0[j][i] + FU0[j][i + 1];
1610 x2 = -dtdx * (FW[j][i + 1] - FW[j][i]);
1611 x3 = -dt2 * (FV1[j][i + 1] + FV1[j][i]) * FDerLin[i];
1612 x4 = -dt2 * (FV2[j][i + 1] + FV2[j][i]);
1613 FU12[j][i] = (x1 + x2 + x3 + x4) / 2;
1618 for(
int i = 0; i < FNin - 1; i++) {
1619 hi12[i] = (Fhi[i] + Fhi[i + 1]) / 2.;
1620 rho12[i] = (Frho[i] + Frho[i + 1]) / 2.;
1621 Re12[i] = (FRe[i] + FRe[i + 1]) / 2.;
1622 TPTubo12[i] = (FTPTubo[0][i] + FTPTubo[0][i + 1]) / 2.;
1623 Gamma12[i] = (FGamma[i] + FGamma[i + 1]) / 2.;
1624 Rmezcla12[i] = (FRMezcla[i] + FRMezcla[i + 1]) / 2.;
1625 Gamma1_12[i] = (FGamma1[i] + FGamma1[i + 1]) / 2.;
1628 CalculaFlujo(FU12, FW, Gamma12, Gamma1_12, Nodos);
1630 CalculaFuente1(FU12, FV1, Gamma12, Gamma1_12, Nodos);
1632 CalculaFuente2(FU12, FV2, FDiametroD12, hi12, rho12, Re12, TPTubo12, Gamma12, Rmezcla12, Gamma1_12, Nodos);
1634 for(
int i = 1; i < FNin - 1; i++) {
1635 for(
int j = 0; j < FNumEcuaciones; j++) {
1637 x2 = -dtdx * (FW[j][i] - FW[j][i - 1]);
1638 x3 = -dt2 * (FV1[j][i] + FV1[j][i - 1]) * FDerLin12[i];
1639 x4 = -dt2 * (FV2[j][i] + FV2[j][i - 1]);
1640 FU1[j][i] = x1 + x2 + x3 + x4;
1653 }
catch(exception & N) {
1654 std::cout <<
"ERROR: TTubo::LaxWendrof en el tubo: " << FNumeroTubo << std::endl;
1655 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1664 void TTubo::FluxCorrectedTransport() {
1665 double c1 = 0., c2 = 0., c3 = 0., c4 = 0., sign = 0.;
1670 for(
int k = 0; k < FNumEcuaciones; k++) {
1671 FU1[k][0] = FU0[k][0];
1672 FU1[k][FNin - 1] = FU0[k][FNin - 1];
1675 for(
int i = 0; i < FNin; i++) {
1676 Transforma3Area(FUfct0, FU0, FArea[i], FGamma[i], FGamma1[i], FGamma6[i], i);
1678 Transforma3Area(FUfct1, FU1, FArea[i], FGamma[i], FGamma1[i], FGamma6[i], i);
1683 if(FMod.Difusion == nmDamping) {
1684 for(
int i = 0; i < FNin - 1; i++) {
1685 for(
int k = 0; k < FNumEcuaciones; k++) {
1686 Ffl[k][i] = (FUfct0[k][i + 1] - FUfct0[k][i]) / 8.;
1689 }
else if(FMod.Difusion == nmSmoothing) {
1690 for(
int i = 0; i < FNin - 1; i++) {
1691 for(
int k = 0; k < FNumEcuaciones; k++) {
1692 Ffl[k][i] = (FUfct1[k][i + 1] - FUfct1[k][i]) / 8.;
1696 std::cout <<
"ERROR: Metodo de Difusion mal definido para el FCT en el tubo n. " << FNumeroTubo << std::endl;
1700 for(
int i = 1; i < FNin - 1; i++) {
1701 for(
int k = 0; k < FNumEcuaciones; k++) {
1702 FdU[k][i] = Ffl[k][i] - Ffl[k][i - 1];
1703 FUfctd[k][i] = FUfct1[k][i] + FdU[k][i];
1707 for(
int k = 0; k < FNumEcuaciones; k++) {
1708 FUfctd[k][0] = FUfct1[k][0];
1709 FUfctd[k][FNin - 1] = FUfct1[k][FNin - 1];
1713 if(FMod.Antidifusion == nmExplicit) {
1714 for(
int i = 1; i < FNin; i++) {
1715 for(
int k = 0; k < FNumEcuaciones; k++) {
1716 FDeltaFCTd[k][i] = FUfctd[k][i] - FUfctd[k][i - 1];
1719 }
else if(FMod.Antidifusion == nmNaive) {
1720 for(
int i = 1; i < FNin; i++) {
1721 for(
int k = 0; k < FNumEcuaciones; k++) {
1722 FDeltaFCTd[k][i] = FUfct0[k][i] - FUfct0[k][i - 1];
1725 }
else if(FMod.Antidifusion == nmPhoenical) {
1726 for(
int i = 1; i < FNin; i++) {
1727 for(
int k = 0; k < FNumEcuaciones; k++) {
1728 FDeltaFCTd[k][i] = FUfct1[k][i] - FUfct1[k][i - 1];
1732 std::cout <<
"ERROR: Metodo de Antidifusion mal definido para el FCT en el tubo n. " << FNumeroTubo << std::endl;
1736 for(
int k = 0; k < FNumEcuaciones; k++) {
1737 FDeltaFCTd[k][0] = FDeltaFCTd[k][1];
1738 FDeltaFCTd[k][FNin] = FDeltaFCTd[k][FNin - 1];
1740 for(
int i = 1; i < FNin - 2; i++) {
1741 for(
int k = 0; k < FNumEcuaciones; k++) {
1742 if(FDeltaFCTd[k][i + 1] >= 0.) {
1747 c1 = 5. * sign * FDeltaFCTd[k][i] / 8.;
1748 c2 = fabs(FDeltaFCTd[k][i + 1]) / 8.;
1749 c3 = 5. * sign * FDeltaFCTd[k][i + 2] / 8.;
1757 FflU[k][i] = sign * c4;
1764 for(
int k = 0; k < FNumEcuaciones; k++) {
1765 if(FDeltaFCTd[k][1] >= 0.) {
1770 c2 = fabs(FDeltaFCTd[k][1]) / 8.;
1771 c3 = 5. * sign * FDeltaFCTd[k][2] / 8.;
1777 FflU[k][0] = sign * c4;
1782 for(
int k = 0; k < FNumEcuaciones; k++) {
1783 if(FDeltaFCTd[k][FNin - 2] >= 0.) {
1788 c1 = 5. * sign * FDeltaFCTd[k][FNin - 2] / 8.;
1789 c2 = fabs(FDeltaFCTd[k][FNin - 1]) / 8.;
1795 FflU[k][FNin - 2] = sign * c4;
1797 FflU[k][FNin - 2] = 0.;
1800 for(
int i = 1; i < FNin - 1; i++) {
1801 for(
int k = 0; k < FNumEcuaciones; k++) {
1802 FaU[k][i] = -FflU[k][i] + FflU[k][i - 1];
1803 FUfctad[k][i] = FUfctd[k][i] + FaU[k][i];
1805 Transforma4Area(FU1, FUfctad, FArea[i], FGamma[i], FGamma1[i], FGamma3[i], FGamma4[i], FGamma6[i], i);
1810 }
catch(exception & N) {
1811 std::cout <<
"ERROR: TTubo::FluxCorrectedTransport en el tubo: " << FNumeroTubo << std::endl;
1812 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1821 void TTubo::LaxWendroffArea() {
1826 double x1, x2, x3, x4, *hi12, *rho12, *Re12, *TPTubo12, *Gamma12, *Rmezcla12, *Gamma1_12;
1828 hi12 =
new double[FNin - 1];
1829 rho12 =
new double[FNin - 1];
1830 Re12 =
new double[FNin - 1];
1831 TPTubo12 =
new double[FNin - 1];
1832 Gamma12 =
new double[FNin - 1];
1833 Rmezcla12 =
new double[FNin - 1];
1834 Gamma1_12 =
new double[FNin - 1];
1842 CalculaFlujo(FU0, FW, FGamma, FGamma1, Nodos);
1844 CalculaFuente1Area(FU0, FV1, FArea, FGamma1, Nodos);
1846 CalculaFuente2Area(FU0, FV2, FArea, Fhi, Frho, FRe, FTPTubo[0], FGamma, FRMezcla, FGamma1, Nodos);
1848 for(
int i = 0; i < FNin - 1; i++) {
1849 for(
int j = 0; j < FNumEcuaciones; j++) {
1850 x1 = (FU0[j][i] + FU0[j][i + 1]);
1851 x2 = -dtdx * (FW[j][i + 1] - FW[j][i]);
1852 x3 = -dt2 * (FV1[j][i + 1] + FV1[j][i]) * FDerLinArea[i];
1853 x4 = -dt2 * (FV2[j][i + 1] + FV2[j][i]);
1854 FU12[j][i] = (x1 + x2 + x3 + x4) / 2.;
1859 for(
int i = 0; i < FNin - 1; i++) {
1860 hi12[i] = (Fhi[i] + Fhi[i + 1]) / 2.;
1861 rho12[i] = (Frho[i] + Frho[i + 1]) / 2.;
1862 Re12[i] = (FRe[i] + FRe[i + 1]) / 2.;
1863 TPTubo12[i] = (FTPTubo[0][i] + FTPTubo[0][i + 1]) / 2.;
1864 Gamma12[i] = (FGamma[i] + FGamma[i + 1]) / 2.;
1865 Rmezcla12[i] = (FRMezcla[i] + FRMezcla[i + 1]) / 2.;
1866 Gamma1_12[i] = (FGamma1[i] + FGamma1[i + 1]) / 2.;
1869 CalculaFlujo(FU12, FW, Gamma12, Gamma1_12, Nodos);
1871 CalculaFuente1Area(FU12, FV1, FArea12, Gamma1_12, Nodos);
1873 CalculaFuente2Area(FU12, FV2, FArea12, hi12, rho12, Re12, TPTubo12, Gamma12, Rmezcla12, Gamma1_12, Nodos);
1875 for(
int i = 1; i < FNin - 1; i++) {
1876 for(
int j = 0; j < FNumEcuaciones; j++) {
1878 x2 = -dtdx * (FW[j][i] - FW[j][i - 1]);
1879 x3 = -dt2 * (FV1[j][i] + FV1[j][i - 1]) * FDerLinArea12[i];
1880 x4 = -dt2 * (FV2[j][i] + FV2[j][i - 1]);
1881 FU1[j][i] = x1 + x2 + x3 + x4;
1894 }
catch(exception & N) {
1895 std::cout <<
"ERROR: TTubo::LaxWendroffArea en el tubo: " << FNumeroTubo << std::endl;
1896 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1905 void TTubo::CalculaFlujo(
double **U,
double **W,
double *Gamma,
double *Gamma1,
int Nodos) {
1913 for(
int i = 0; i < Nodos; i++) {
1914 if(!DoubEqZero(U[1][i])) {
1915 U1U0 = U[1][i] / U[0][i];
1918 W[1][i] = U[2][i] * Gamma1[i] - (Gamma[i] - 3.0) * U[1][i] * U1U0 / 2.;
1919 W[2][i] = Gamma[i] * U[2][i] * U1U0 - Gamma1[i] * U[1][i] *
pow2(U1U0) / 2.;
1920 for(
int j = 3; j < FNumEcuaciones; j++) {
1921 W[j][i] = U[j][i] * U1U0;
1925 W[1][i] = U[2][i] * Gamma1[i];
1927 for(
int j = 3; j < FNumEcuaciones; j++) {
1933 }
catch(exception & N) {
1934 std::cout <<
"ERROR: TTubo::CalculaFlujo en el tubo: " << FNumeroTubo << std::endl;
1935 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1944 void TTubo::CalculaFuente1(
double **U,
double **V1,
double *Gamma,
double *Gamma1,
int Nodos) {
1949 for(
int i = 0; i < Nodos; i++) {
1950 U1U0 = U[1][i] / U[0][i];
1953 V1[1][i] = U[1][i] * U1U0;
1954 V1[2][i] = Gamma[i] * U[2][i] * U1U0 - Gamma1[i] * U[1][i] *
pow2(U1U0) / 2.;
1955 for(
int j = 3; j < FNumEcuaciones; j++) {
1956 V1[j][i] = U[j][i] * U1U0;
1960 }
catch(exception & N) {
1961 std::cout <<
"ERROR: TTubo::CalculaFuente1 en el tubo: " << FNumeroTubo << std::endl;
1962 std::cout <<
"Tipo de error: " << N.what() << std::endl;
1971 void TTubo::CalculaFuente1Area(
double **U,
double **V1,
double *Area,
double *Gamma1,
int Nodos) {
1976 for(
int i = 0; i < Nodos; i++) {
1977 if(DoubEqZero(U[1][i])) {
1978 p = (U[2][i]) * Gamma1[i] / Area[i];
1980 p = (U[2][i] -
pow2(U[1][i]) / U[0][i] / 2.0) * Gamma1[i] / Area[i];
1986 for(
int j = 3; j < FNumEcuaciones; j++) {
1991 }
catch(exception & N) {
1992 std::cout <<
"ERROR: TTubo::CalculaFuente1Area en el tubo: " << FNumeroTubo << std::endl;
1993 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2002 void TTubo::CalculaFuente2(
double **U,
double **V2,
double *diame,
double *hi,
double *rho,
double *Re,
2003 double *TempParedTubo,
double *Gamma,
double *Rmezcla,
double *Gamma1,
int Nodos) {
2004 double v = 0, a = 0., p = 0., tgas = 0., g = 0., q = 0., f = 0.;
2008 for(
int i = 0; i < Nodos; i++) {
2010 v = U[1][i] / U[0][i];
2011 p = (U[2][i] - U[1][i] * v / 2.0) * Gamma1[i];
2012 a = sqrt(Gamma[i] * p / U[0][i]);
2013 if(v > 1e200 || v < -1e200) {
2014 std::cout <<
"ERROR: Valor de velocidad no valido" << std::endl;
2017 if(p > 1e200 || p < 0) {
2018 std::cout <<
"ERROR: Valor de presion no valido en el tubo " << FNumeroTubo <<
" nodo " << i << std::endl;
2021 if(a > 1e200 || a < 0) {
2022 std::cout <<
"ERROR: Valor de velocidad del sonido no valido" << std::endl;
2023 throw Exception(
"Error velocidad del sonido");
2027 if(DoubEqZero(v) || DoubEqZero(FCoefAjusFric)) {
2030 Colebrook(FFriccion, diame[i], f, Re[i]);
2031 g = f * v * v * v / fabs(v) * 2 / diame[i] * FCoefAjusFric;
2035 if(DoubEqZero(FCoefAjusTC)) {
2038 tgas = a * a / (Gamma[i] * Rmezcla[i]);
2039 TransmisionCalor(tgas, diame[i], q, hi[i], rho[i], TempParedTubo[i]);
2040 q = q * FCoefAjusTC;
2045 V2[1][i] = U[0][i] * g;
2046 V2[2][i] = -U[0][i] * q;
2047 for(
int j = 3; j < FNumEcuaciones; j++) {
2052 }
catch(exception & N) {
2053 std::cout <<
"ERROR: TTubo::CalculaFuente1 en el tubo: " << FNumeroTubo << std::endl;
2054 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2063 void TTubo::CalculaFuente2Area(
double **U,
double **V2,
double *Area,
double *hi,
double *rho,
double *Re,
2064 double *TempParedTubo,
double *Gamma,
double *Rmezcla,
double *Gamma1,
int Nodos) {
2065 double v = 0., a = 0., pA = 0., tgas = 0., g = 0., q = 0., f = 0.;
2070 for(
int i = 0; i < Nodos; i++) {
2072 v = U[1][i] / U[0][i];
2073 pA = (U[2][i] - U[1][i] * v / 2.0) * Gamma1[i];
2074 a = sqrt(Gamma[i] * pA / U[0][i]);
2076 diame = sqrt(Area[i] * __cons::_4_Pi);
2078 if(DoubEqZero(v) || DoubEqZero(FCoefAjusFric)) {
2081 Colebrook(FFriccion, diame, f, Re[i]);
2082 g = f * v * v * v / fabs(v) * 2 / diame * FCoefAjusFric;
2086 if(DoubEqZero(FCoefAjusTC)) {
2089 tgas = a * a / (Gamma[i] * Rmezcla[i]);
2090 TransmisionCalor(tgas, diame, q, hi[i], rho[i], TempParedTubo[i]);
2091 q = q * FCoefAjusTC;
2096 V2[1][i] = U[0][i] * g;
2097 V2[2][i] = -U[0][i] * q;
2098 for(
int j = 3; j < FNumEcuaciones; j++) {
2103 }
catch(exception & N) {
2104 std::cout <<
"ERROR: TTubo::CalculaFuente1 en el tubo: " << FNumeroTubo << std::endl;
2105 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2114 void TTubo::Colebrook(
double rug,
double dia,
double& f,
double Re) {
2115 double temp = 0., temp2 = 0.;
2117 double pow_approx = 0.;
2121 temp = rug / (3700.0 * dia) + 5.74 / pow(Re, 0.9);
2122 }
else if(Re > 50000.) {
2125 pow_approx = 7.04596786e-13 *
pow3(Re) - 3.55878012e-07 *
pow2(Re) + 3.34536053e-01 * Re + 1.02365893e+03;
2126 temp = rug / (3700.0 * dia) + 5.74 / pow_approx;
2127 }
else if(Re > 20000.) {
2130 pow_approx = 5.09028031e-12 *
pow3(Re) - 1.00134262e-06 *
pow2(Re) + 3.67476796e-01 * Re + 4.40171892e+02;
2131 temp = rug / (3700.0 * dia) + 5.74 / pow_approx;
2132 }
else if(Re > 10000.) {
2135 pow_approx = 2.92727184e-11 *
pow3(Re) - 2.48693415e-06 *
pow2(Re) + 3.98901290e-01 * Re + 2.11628209e+02;
2136 temp = rug / (3700.0 * dia) + 5.74 / pow_approx;
2137 }
else if(Re > 4000.) {
2140 pow_approx = 1.49478492e-10 *
pow3(Re) - 5.88098006e-06 *
pow2(Re) + 4.31645206e-01 * Re + 1.03406848e+02;
2141 temp = rug / (3700.0 * dia) + 5.74 / pow_approx;
2143 temp = rug / (3700.0 * dia) + 0.00328895476345;
2146 temp2 =
pow2(log10(temp));
2147 }
else if(temp > 0.019) {
2150 temp2 = 6.76495011e+05 *
pow4(temp) - 1.07370937e+05 *
pow3(temp) + 6.95435632e+03 *
pow2(
2151 temp) - 2.44427511e+02 * temp + 5.74410691e+00;
2152 }
else if(temp > 0.008) {
2155 temp2 = 2.08516955e+07 *
pow4(temp) - 1.46552875e+06 *
pow3(temp) + 4.20051994e+04 *
pow2(
2156 temp) - 6.58067557e+02 * temp + 7.63733377e+00;
2157 }
else if(temp > 0.00329) {
2160 temp2 = 7.84991530e+08 *
pow4(temp) - 2.30952512e+07 *
pow3(temp) + 2.77321091e+05 *
pow2(
2161 temp) - 1.83296385e+03 * temp + 9.92239603e+00;
2162 }
else if(temp > 0.0013) {
2165 temp2 = 3.29050277e+10 *
pow4(temp) - 3.93727194e+08 *
pow3(temp) + 1.92219159e+06 *
pow2(
2166 temp) - 5.18844300e+03 * temp + 1.25952908e+01;
2167 }
else if(temp > 0.0005) {
2170 temp2 = 1.57146781e+12 *
pow4(temp) - 7.37699295e+09 *
pow3(temp) + 1.41288296e+07 *
pow2(
2171 temp) - 1.50182432e+04 * temp + 1.56956366e+01;
2172 }
else if(temp > 0.00019) {
2175 temp2 = 8.09890844e+13 *
pow4(temp) - 1.45825030e+11 *
pow3(temp) + 1.07169184e+08 *
pow2(
2176 temp) - 4.38671762e+04 * temp + 1.92059658e+01;
2177 }
else if(temp > 0.00007) {
2180 temp2 = 4.47104917e+15 *
pow4(temp) - 3.03334424e+12 *
pow3(temp) + 8.39454836e+08 *
pow2(
2181 temp) - 1.29630353e+05 * temp + 2.31540981e+01;
2183 temp2 =
pow2(log10(temp));
2196 void TTubo::TransmisionCalor(
double tgas,
double diametro,
double& q,
double hi,
double rho,
double Tw) {
2202 q = 4. * hi * (__units::degCToK(Tw) - tgas) / rho / diametro;
2206 }
catch(exception & N) {
2207 std::cout <<
"ERROR: TTubo::TransmisionCalor en el tubo: " << FNumeroTubo << std::endl;
2208 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2217 inline double TTubo::DerLinF(
double d1,
double d2,
double xref) {
2223 return 2. * (d2 - d1) / ((d1 + d2) / 2.0) / xref;
2226 }
catch(exception & N) {
2227 std::cout <<
"ERROR: TTubo::DerLinF en el tubo: " << FNumeroTubo << std::endl;
2228 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2237 inline double TTubo::DerLinFArea(
double area1,
double area2,
double xref) {
2242 return (area2 - area1) / xref;
2245 }
catch(exception & N) {
2246 std::cout <<
"ERROR: TTubo::DerLinFArea en el tubo: " << FNumeroTubo << std::endl;
2247 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2261 double a = 0., v = 0., p = 0.;
2262 double LandaIzq = 0., BetaIzq = 0., EntropiaIzq = 0.;
2263 double LandaDer = 0., BetaDer = 0., EntropiaDer = 0.;
2264 double *YIzq, *YDer;
2266 YIzq =
new double[FNumeroEspecies - FIntEGR];
2267 YDer =
new double[FNumeroEspecies - FIntEGR];
2269 LandaIzq = BC[FNodoIzq - 1]->GetTuboExtremo(FTuboCCNodoIzq).Landa;
2270 BetaIzq = BC[FNodoIzq - 1]->GetTuboExtremo(FTuboCCNodoIzq).Beta;
2271 EntropiaIzq = BC[FNodoIzq - 1]->GetTuboExtremo(FTuboCCNodoIzq).Entropia;
2273 TransformaContorno(LandaIzq, BetaIzq, EntropiaIzq, a, v, p, 1, FGamma1[0], FGamma3[0], FGamma4[0], FGamma5[0]);
2274 if(BC[FNodoIzq - 1]->getTipoCC() == nmBranch) {
2276 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
2277 YIzq[i] = FFraccionMasicaCC[0][i];
2280 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
2281 YIzq[i] = BC[FNodoIzq - 1]->GetFraccionMasicaEspecie(i);
2285 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
2286 YIzq[i] = BC[FNodoIzq - 1]->GetFraccionMasicaEspecie(i);
2290 if(FMod.FormulacionLeyes == nmSinArea)
2291 Transforma1(v, a, p, FU1, FGamma[0], FGamma1[0], YIzq, 0);
2292 else if(FMod.FormulacionLeyes == nmConArea)
2293 Transforma1Area(v, a, p, FU1, FArea[0], FGamma[0], FGamma1[0], YIzq, 0);
2295 LandaDer = BC[FNodoDer - 1]->GetTuboExtremo(FTuboCCNodoDer).Landa;
2296 BetaDer = BC[FNodoDer - 1]->GetTuboExtremo(FTuboCCNodoDer).Beta;
2297 EntropiaDer = BC[FNodoDer - 1]->GetTuboExtremo(FTuboCCNodoDer).Entropia;
2299 TransformaContorno(LandaDer, BetaDer, EntropiaDer, a, v, p, 1, FGamma1[FNin - 1], FGamma3[FNin - 1], FGamma4[FNin - 1],
2301 if(BC[FNodoDer - 1]->getTipoCC() == nmBranch) {
2303 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
2304 YDer[i] = FFraccionMasicaCC[1][i];
2307 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
2308 YDer[i] = BC[FNodoDer - 1]->GetFraccionMasicaEspecie(i);
2312 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
2313 YDer[i] = BC[FNodoDer - 1]->GetFraccionMasicaEspecie(i);
2317 if(FMod.FormulacionLeyes == nmSinArea)
2318 Transforma1(v, a, p, FU1, FGamma[FNin - 1], FGamma1[FNin - 1], YDer, FNin - 1);
2319 else if(FMod.FormulacionLeyes == nmConArea)
2320 Transforma1Area(v, a, p, FU1, FArea[FNin - 1], FGamma[FNin - 1], FGamma1[FNin - 1], YDer, FNin - 1);
2322 if(FMod.FormulacionLeyes == nmSinArea) {
2323 for(
int i = 0; i < FNin; i++) {
2324 Transforma2(FVelocidad0[i], FAsonido0[i], FPresion0[i], FU1, FGamma[i], FGamma1[i], FFraccionMasicaEspecie[i], i);
2325 for(
int k = 0; k < FNumEcuaciones; k++) {
2326 FU0[k][i] = FU1[k][i];
2328 FVelocidadDim[i] = FVelocidad0[i] * __cons::ARef;
2329 FAsonidoDim[i] = FAsonido0[i] * __cons::ARef;
2330 FFlowMass[i] = FU0[1][i] * FArea[i];
2332 }
else if(FMod.FormulacionLeyes == nmConArea) {
2333 for(
int i = 0; i < FNin; i++) {
2334 Transforma2Area(FVelocidad0[i], FAsonido0[i], FPresion0[i], FU1, FArea[i], FGamma[i], FGamma1[i],
2335 FFraccionMasicaEspecie[i], i);
2336 for(
int k = 0; k < FNumEcuaciones; k++) {
2337 FU0[k][i] = FU1[k][i];
2339 FVelocidadDim[i] = FVelocidad0[i] * __cons::ARef;
2340 FAsonidoDim[i] = FAsonido0[i] * __cons::ARef;
2341 FFlowMass[i] = FU0[1][i];
2342 if(FVelocidadDim[i] > FAsonidoDim[i] + 1.0e-10) {
2343 printf(
"Supersonic flow in pipe: %d node: %d, Mach = %lf\n", FNumeroTubo, i, FVelocidadDim[i] / FAsonidoDim[i]);
2354 catch(exception & N) {
2355 std::cout <<
"ERROR: TTubo::ValoresDeContorno tubo: " << FNumeroTubo << std::endl;
2356 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2364 void TTubo::TransformaContorno(
double& L,
double& B,
double& E,
double& a,
double& v,
double& p,
const int& modo,
2365 const double& Gamma1,
const double& Gamma3,
const double& Gamma4,
const double& Gamma5) {
2370 L = (a + Gamma3 * v);
2371 B = (a - Gamma3 * v);
2372 E = a / pow(p, Gamma5);
2375 v = (L - B) / Gamma1;
2376 p = pow(a / E, Gamma4);
2379 }
catch(exception & N) {
2380 std::cout <<
"ERROR: TTubo::TransformaContorno en el tubo: " << FNumeroTubo << std::endl;
2381 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2391 double Machx = 0., Machy = 0., Velocidady = 0., Sonidoy = 0.;
2396 for(
int i = 0; i < FNin; i++) {
2397 Machx = FVelocidad0[i] / FAsonido0[i];
2398 if(-1. >= Machx || Machx > 1.) {
2399 Machy = Machx / fabs(Machx) * sqrt((
pow2(Machx) + 2. / FGamma1[i]) / (FGamma4[i] *
pow2(Machx) - 1.));
2400 Sonidoy = FAsonido0[i] * sqrt((FGamma1[i] / 2. *
pow2(Machx) + 1.) / (FGamma1[i] / 2. *
pow2(Machy) + 1.));
2402 Velocidady = Sonidoy * Machy;
2403 FAsonido0[i] = Sonidoy;
2404 FVelocidad0[i] = Velocidady;
2408 }
catch(exception & N) {
2409 std::cout <<
"ERROR: TTubo::ReduccionFlujoSubsonico en el tubo: " << FNumeroTubo << std::endl;
2410 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2420 double Machx = 0., Machy = 0., Velocidady = 0., Sonidoy = 0.;
2421 double velocidad = 0., asonido = 0., presion = 0.;
2426 for(
int i = 1; i < FNin - 1; i++) {
2427 velocidad = FU1[1][i] / FU1[0][i] / __cons::ARef;
2428 presion = __units::PaToBar((FU1[2][i] - FU1[1][i] * velocidad * __cons::ARef / 2.0) * FGamma1[i] / FArea[i]);
2429 asonido = sqrt(FGamma[i] * __units::BarToPa(presion) * FArea[i] / FU1[0][i] / __cons::ARef2);
2430 Machx = velocidad / asonido;
2431 if(-1. >= Machx || Machx > 1.) {
2432 Machy = Machx / fabs(Machx) * sqrt((
pow2(Machx) + 2. / FGamma1[i]) / (FGamma4[i] *
pow2(Machx) - 1.));
2433 Sonidoy = asonido * sqrt((FGamma1[i] / 2. *
pow2(Machx) + 1.) / (FGamma1[i] / 2. *
pow2(Machy) + 1.));
2435 Velocidady = Sonidoy * Machy;
2437 velocidad = Velocidady;
2438 FU1[0][i] = FGamma[i] * __units::BarToPa(presion) /
pow2(asonido * __cons::ARef) * FArea[i];
2439 FU1[1][i] = FU1[0][i] * velocidad * __cons::ARef;
2440 FU1[2][i] = FArea[i] * __units::BarToPa(presion) / FGamma1[i] + FU1[1][i] * velocidad * __cons::ARef / 2.0;
2445 }
catch(exception & N) {
2446 std::cout <<
"ERROR: TTubo::ReduccionFlujoSubsonicoFCT en el tubo: " << FNumeroTubo << std::endl;
2447 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2457 int NumVars = 0, TipoVar = 0;
2462 FILE *fich = fopen(FileWAM,
"r");
2463 fsetpos(fich, &filepos);
2465 fscanf(fich,
"%d ", &FNumResMedios);
2470 for(
int i = 0; i < FNumResMedios; i++) {
2471 ResultadosMedios[i].TemperaturaGas =
false;
2472 ResultadosMedios[i].TemperaturaGasSUM = 0.;
2473 ResultadosMedios[i].TemperaturaGasMED = 0;
2474 ResultadosMedios[i].Pressure =
false;
2475 ResultadosMedios[i].PresionSUM = 0.;
2476 ResultadosMedios[i].PresionMED = 0.;
2477 ResultadosMedios[i].Velocity =
false;
2478 ResultadosMedios[i].VelocidadSUM = 0.;
2479 ResultadosMedios[i].VelocidadMED = 0.;
2480 ResultadosMedios[i].Massflow =
false;
2481 ResultadosMedios[i].GastoSUM = 0.;
2482 ResultadosMedios[i].GastoMED = 0.;
2483 ResultadosMedios[i].TemperaturaInternaPared =
false;
2484 ResultadosMedios[i].TemperaturaInternaParedSUM = 0.;
2485 ResultadosMedios[i].TemperaturaInternaParedMED = 0.;
2486 ResultadosMedios[i].TemperaturaIntermediaPared =
false;
2487 ResultadosMedios[i].TemperaturaIntermediaParedSUM = 0.;
2488 ResultadosMedios[i].TemperaturaIntermediaParedMED = 0.;
2489 ResultadosMedios[i].TemperaturaExternaPared =
false;
2490 ResultadosMedios[i].TemperaturaExternaParedSUM = 0.;
2491 ResultadosMedios[i].TemperaturaExternaParedMED = 0.;
2492 ResultadosMedios[i].NITmedio =
false;
2493 ResultadosMedios[i].NITmedioSUM = 0;
2494 ResultadosMedios[i].NITmedioMED = 0;
2495 ResultadosMedios[i].CoefPelInterior =
false;
2496 ResultadosMedios[i].CoefPelInteriorSUM = 0.;
2497 ResultadosMedios[i].CoefPelInteriorMED = 0.;
2498 ResultadosMedios[i].FraccionMasicaEspecies =
false;
2499 ResultadosMedios[i].FraccionSUM =
new double[FNumeroEspecies - FIntEGR];
2500 ResultadosMedios[i].FraccionMED =
new double[FNumeroEspecies - FIntEGR];
2501 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
2502 ResultadosMedios[i].FraccionSUM[j] = 0.;
2503 ResultadosMedios[i].FraccionMED[j] = 0.;
2505 ResultadosMedios[i].PonderacionSUM = 0.;
2506 ResultadosMedios[i].GastoPonderacionSUM = 0.;
2508 fscanf(fich,
"%lf %d ", &ResultadosMedios[i].Distancia, &NumVars);
2510 for(
int j = 0; j < NumVars; j++) {
2511 fscanf(fich,
"%d ", &TipoVar);
2514 ResultadosMedios[i].TemperaturaGas =
true;
2517 ResultadosMedios[i].Pressure =
true;
2520 ResultadosMedios[i].Velocity =
true;
2523 ResultadosMedios[i].Massflow =
true;
2526 ResultadosMedios[i].TemperaturaInternaPared =
true;
2529 ResultadosMedios[i].TemperaturaIntermediaPared =
true;
2532 ResultadosMedios[i].TemperaturaExternaPared =
true;
2535 ResultadosMedios[i].NITmedio =
true;
2537 std::cout <<
"ERROR: No puedes pedir el NIT como resultado si no hay motor" << std::endl;
2542 ResultadosMedios[i].CoefPelInterior =
true;
2545 ResultadosMedios[i].FraccionMasicaEspecies =
true;
2548 std::cout <<
"WARNING: El tipo de variable seleccionada para la salida de resultados medios no es valida" << std::endl;
2555 fgetpos(fich, &filepos);
2558 }
catch(exception & N) {
2559 std::cout <<
"ERROR: TTubo::ReadAverageResults en el tubo: " << FNumeroTubo << std::endl;
2560 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2575 std::ostringstream TextDist;
2576 TextDist.precision(8);
2578 for(
int i = 0; i < FNumResMedios; i++) {
2579 TextDist << ResultadosMedios[i].Distancia;
2581 if(ResultadosMedios[i].TemperaturaGas) {
2584 medoutput << Label.c_str();
2586 if(ResultadosMedios[i].Pressure) {
2589 medoutput << Label.c_str();
2591 if(ResultadosMedios[i].Velocity) {
2594 medoutput << Label.c_str();
2596 if(ResultadosMedios[i].Massflow) {
2599 medoutput << Label.c_str();
2601 if(ResultadosMedios[i].TemperaturaInternaPared) {
2604 medoutput << Label.c_str();
2606 if(ResultadosMedios[i].TemperaturaIntermediaPared) {
2609 medoutput << Label.c_str();
2611 if(ResultadosMedios[i].TemperaturaExternaPared) {
2614 medoutput << Label.c_str();
2616 if(ResultadosMedios[i].NITmedio) {
2619 medoutput << Label.c_str();
2621 if(ResultadosMedios[i].CoefPelInterior) {
2624 medoutput << Label.c_str();
2626 if(ResultadosMedios[i].FraccionMasicaEspecies) {
2627 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
2630 medoutput << Label.c_str();
2636 }
catch(exception & N) {
2637 std::cout <<
"ERROR: TTubo::HeaderAverageResults en el tubo: " << FNumeroTubo << std::endl;
2638 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2654 for(
int i = 0; i < FNumResMedios; i++) {
2655 if(ResultadosMedios[i].TemperaturaGas)
2656 medoutput <<
"\t" << __units::KTodegC(ResultadosMedios[i].TemperaturaGasMED);
2657 if(ResultadosMedios[i].Pressure)
2658 medoutput <<
"\t" << ResultadosMedios[i].PresionMED;
2659 if(ResultadosMedios[i].Velocity)
2660 medoutput <<
"\t" << ResultadosMedios[i].VelocidadMED;
2661 if(ResultadosMedios[i].Massflow)
2662 medoutput <<
"\t" << ResultadosMedios[i].GastoMED;
2663 if(ResultadosMedios[i].TemperaturaInternaPared)
2664 medoutput <<
"\t" << ResultadosMedios[i].TemperaturaInternaParedMED;
2665 if(ResultadosMedios[i].TemperaturaIntermediaPared)
2666 medoutput <<
"\t" << ResultadosMedios[i].TemperaturaIntermediaParedMED;
2667 if(ResultadosMedios[i].TemperaturaExternaPared)
2668 medoutput <<
"\t" << ResultadosMedios[i].TemperaturaExternaParedMED;
2669 if(ResultadosMedios[i].NITmedio)
2670 medoutput <<
"\t" << ResultadosMedios[i].NITmedioMED;
2671 if(ResultadosMedios[i].CoefPelInterior)
2672 medoutput <<
"\t" << ResultadosMedios[i].CoefPelInteriorMED;
2673 if(ResultadosMedios[i].FraccionMasicaEspecies) {
2674 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
2675 medoutput <<
"\t" << ResultadosMedios[i].FraccionMED[j];
2681 }
catch(exception & N) {
2682 std::cout <<
"ERROR: TTubo::ResultadosMedios en el tubo: " << FNumeroTubo << std::endl;
2683 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2692 int NumVars = 0, TipoVar = 0;
2697 FILE *fich = fopen(FileWAM,
"r");
2698 fsetpos(fich, &filepos);
2700 fscanf(fich,
"%d ", &FNumResInstant);
2703 for(
int i = 0; i < FNumResInstant; i++) {
2704 ResultInstantaneos[i].Pressure =
false;
2705 ResultInstantaneos[i].PresionINS = 0.;
2706 ResultInstantaneos[i].Velocity =
false;
2707 ResultInstantaneos[i].VelocidadINS = 0.;
2708 ResultInstantaneos[i].TemperaturaGas =
false;
2709 ResultInstantaneos[i].TemperaturaGasINS = 0.;
2710 ResultInstantaneos[i].FlujoMasico =
false;
2711 ResultInstantaneos[i].FlujoMasicoINS = 0.;
2712 ResultInstantaneos[i].VelocidadDerecha =
false;
2713 ResultInstantaneos[i].VelocidadDerechaINS = 0.;
2714 ResultInstantaneos[i].VelocidadIzquierda =
false;
2715 ResultInstantaneos[i].VelocidadIzquierdaINS = 0.;
2716 ResultInstantaneos[i].PresionDerecha =
false;
2717 ResultInstantaneos[i].PresionDerechaINS = 0.;
2718 ResultInstantaneos[i].PresionIzquierda =
false;
2719 ResultInstantaneos[i].PresionIzquierdaINS = 0.;
2720 ResultInstantaneos[i].NIT =
false;
2721 ResultInstantaneos[i].NITINS = 0.;
2722 ResultInstantaneos[i].TemperaturaInternaPared =
false;
2723 ResultInstantaneos[i].TemperaturaInternaParedINS = 0.;
2724 ResultInstantaneos[i].TemperaturaIntermediaPared =
false;
2725 ResultInstantaneos[i].TemperaturaIntermediaParedINS = 0.;
2726 ResultInstantaneos[i].TemperaturaExternaPared =
false;
2727 ResultInstantaneos[i].TemperaturaExternaParedINS = 0.;
2728 ResultInstantaneos[i].CoefPelInterior =
false;
2729 ResultInstantaneos[i].CoefPelInteriorINS = 0.;
2730 ResultInstantaneos[i].FraccionMasicaEspecies =
false;
2731 ResultInstantaneos[i].FraccionINS =
new double[FNumeroEspecies - FIntEGR];
2732 ResultInstantaneos[i].Gamma =
false;
2733 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
2734 ResultInstantaneos[i].FraccionINS[j] = 0.;
2737 fscanf(fich,
"%lf %d ", &ResultInstantaneos[i].Distancia, &NumVars);
2739 for(
int j = 0; j < NumVars; j++) {
2740 fscanf(fich,
"%d ", &TipoVar);
2743 ResultInstantaneos[i].Pressure =
true;
2746 ResultInstantaneos[i].Velocity =
true;
2749 ResultInstantaneos[i].TemperaturaGas =
true;
2752 ResultInstantaneos[i].FlujoMasico =
true;
2755 ResultInstantaneos[i].VelocidadDerecha =
true;
2758 ResultInstantaneos[i].VelocidadIzquierda =
true;
2761 ResultInstantaneos[i].PresionDerecha =
true;
2764 ResultInstantaneos[i].PresionIzquierda =
true;
2767 ResultInstantaneos[i].NIT =
true;
2769 std::cout <<
"ERROR: No puedes pedir el NIT como resultado si no hay motor" << std::endl;
2774 ResultInstantaneos[i].TemperaturaInternaPared =
true;
2777 ResultInstantaneos[i].TemperaturaIntermediaPared =
true;
2780 ResultInstantaneos[i].TemperaturaExternaPared =
true;
2783 ResultInstantaneos[i].CoefPelInterior =
true;
2786 ResultInstantaneos[i].FraccionMasicaEspecies =
true;
2789 ResultInstantaneos[i].Gamma =
true;
2792 std::cout <<
"WARNING: El tipo de variable seleccionada para la salida de resultados instantaneos no es valida" <<
2798 fgetpos(fich, &filepos);
2801 }
catch(exception & N) {
2802 std::cout <<
"ERROR: TTubo::ReadInstantaneousResults en el tubo: " << FNumeroTubo << std::endl;
2803 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2818 std::ostringstream TextDist;
2819 TextDist.precision(8);
2821 for(
int i = 0; i < FNumResInstant; i++) {
2822 TextDist << ResultInstantaneos[i].Distancia;
2823 if(ResultInstantaneos[i].Pressure) {
2826 insoutput << Label.c_str();
2828 if(ResultInstantaneos[i].Velocity) {
2831 insoutput << Label.c_str();
2833 if(ResultInstantaneos[i].TemperaturaGas) {
2836 insoutput << Label.c_str();
2838 if(ResultInstantaneos[i].FlujoMasico) {
2841 insoutput << Label.c_str();
2843 if(ResultInstantaneos[i].VelocidadDerecha) {
2846 insoutput << Label.c_str();
2848 if(ResultInstantaneos[i].VelocidadIzquierda) {
2851 insoutput << Label.c_str();
2853 if(ResultInstantaneos[i].PresionDerecha) {
2856 insoutput << Label.c_str();
2858 if(ResultInstantaneos[i].PresionIzquierda) {
2861 insoutput << Label.c_str();
2863 if(ResultInstantaneos[i].NIT) {
2866 insoutput << Label.c_str();
2868 if(ResultInstantaneos[i].TemperaturaInternaPared) {
2871 insoutput << Label.c_str();
2873 if(ResultInstantaneos[i].TemperaturaIntermediaPared) {
2876 insoutput << Label.c_str();
2878 if(ResultInstantaneos[i].TemperaturaExternaPared) {
2881 insoutput << Label.c_str();
2883 if(ResultInstantaneos[i].CoefPelInterior) {
2886 insoutput << Label.c_str();
2888 if(ResultInstantaneos[i].FraccionMasicaEspecies) {
2889 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
2892 insoutput << Label.c_str();
2895 if(ResultInstantaneos[i].Gamma) {
2898 insoutput << Label.c_str();
2903 }
catch(exception & N) {
2904 std::cout <<
"ERROR: TTubo::HeaderInstantaneousResults en el tubo nº: " << FNumeroTubo << std::endl;
2905 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2919 for(
int i = 0; i < FNumResInstant; i++) {
2920 if(ResultInstantaneos[i].Pressure)
2921 insoutput <<
"\t" << ResultInstantaneos[i].PresionINS;
2922 if(ResultInstantaneos[i].Velocity)
2923 insoutput <<
"\t" << ResultInstantaneos[i].VelocidadINS;
2924 if(ResultInstantaneos[i].TemperaturaGas)
2925 insoutput <<
"\t" << ResultInstantaneos[i].TemperaturaGasINS;
2926 if(ResultInstantaneos[i].FlujoMasico)
2927 insoutput <<
"\t" << ResultInstantaneos[i].FlujoMasicoINS;
2928 if(ResultInstantaneos[i].VelocidadDerecha)
2929 insoutput <<
"\t" << ResultInstantaneos[i].VelocidadDerechaINS;
2930 if(ResultInstantaneos[i].VelocidadIzquierda)
2931 insoutput <<
"\t" << ResultInstantaneos[i].VelocidadIzquierdaINS;
2932 if(ResultInstantaneos[i].PresionDerecha)
2933 insoutput <<
"\t" << ResultInstantaneos[i].PresionDerechaINS;
2934 if(ResultInstantaneos[i].PresionIzquierda)
2935 insoutput <<
"\t" << ResultInstantaneos[i].PresionIzquierdaINS;
2936 if(ResultInstantaneos[i].NIT)
2937 insoutput <<
"\t" << ResultInstantaneos[i].NITINS;
2938 if(ResultInstantaneos[i].TemperaturaInternaPared)
2939 insoutput <<
"\t" << ResultInstantaneos[i].TemperaturaInternaParedINS;
2940 if(ResultInstantaneos[i].TemperaturaIntermediaPared)
2941 insoutput <<
"\t" << ResultInstantaneos[i].TemperaturaIntermediaParedINS;
2942 if(ResultInstantaneos[i].TemperaturaExternaPared)
2943 insoutput <<
"\t" << ResultInstantaneos[i].TemperaturaExternaParedINS;
2944 if(ResultInstantaneos[i].CoefPelInterior)
2945 insoutput <<
"\t" << ResultInstantaneos[i].CoefPelInteriorINS;
2946 if(ResultInstantaneos[i].FraccionMasicaEspecies) {
2947 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
2948 insoutput <<
"\t" << ResultInstantaneos[i].FraccionINS[j];
2951 if(ResultInstantaneos[i].Gamma)
2952 insoutput <<
"\t" << ResultInstantaneos[i].GammaINS;
2957 }
catch(exception & N) {
2958 std::cout <<
"ERROR: TTubo::ResultadosInstantaneos en el tubo nº: " << FNumeroTubo << std::endl;
2959 std::cout <<
"Tipo de error: " << N.what() << std::endl;
2969 double dist = 0., Vble = 0., d = 0., Rmezcla = 0., Gamma = 0., GastoPonderacion = 0.;
2977 for(
int i = 0; i < FNumResMedios; i++) {
2978 dist = ResultadosMedios[i].Distancia / FXref;
2979 n1 = (int) floor(dist);
2980 if(n1 >= FNin - 1) {
2981 if(ResultadosMedios[i].TemperaturaGas || ResultadosMedios[i].FraccionMasicaEspecies) {
2982 GastoPonderacion = fabs(FFlowMass[FNin - 1]);
2983 ResultadosMedios[i].PonderacionSUM += FFlowMass[FNin - 1] *
FDeltaTime;
2984 ResultadosMedios[i].GastoPonderacionSUM += GastoPonderacion;
2986 if(ResultadosMedios[i].TemperaturaGas)
2987 ResultadosMedios[i].TemperaturaGasSUM += (FTemperature[FNin - 1]) * GastoPonderacion;
2988 if(ResultadosMedios[i].Pressure)
2989 ResultadosMedios[i].PresionSUM += FPresion0[FNin - 1] *
FDeltaTime;
2990 if(ResultadosMedios[i].Velocity)
2991 ResultadosMedios[i].VelocidadSUM += FVelocidadDim[FNin - 1] *
FDeltaTime;
2992 if(ResultadosMedios[i].Massflow) {
2993 ResultadosMedios[i].GastoSUM += FFlowMass[i] *
FDeltaTime;
2995 if(ResultadosMedios[i].TemperaturaInternaPared)
2996 ResultadosMedios[i].TemperaturaInternaParedSUM += FTPTubo[0][FNin - 1] *
FDeltaTime;
2997 if(ResultadosMedios[i].TemperaturaIntermediaPared)
2998 ResultadosMedios[i].TemperaturaIntermediaParedSUM += FTPTubo[1][FNin - 1] *
FDeltaTime;
2999 if(ResultadosMedios[i].TemperaturaExternaPared)
3000 ResultadosMedios[i].TemperaturaExternaParedSUM += FTPTubo[2][FNin - 1] *
FDeltaTime;
3001 if(ResultadosMedios[i].NITmedio) {
3002 double nit = CalculaNIT(FAsonido0[FNin - 1], FVelocidad0[FNin - 1], FPresion0[FNin - 1], FDiametroTubo[FNin - 1],
3003 FGamma[FNin - 1], FRMezcla[FNin - 1]);
3004 ResultadosMedios[i].NITmedioSUM += nit *
FDeltaTime;
3006 if(ResultadosMedios[i].CoefPelInterior) {
3007 ResultadosMedios[i].CoefPelInteriorSUM += Fhi[FNin - 1] *
FDeltaTime;
3009 if(ResultadosMedios[i].FraccionMasicaEspecies) {
3010 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
3011 ResultadosMedios[i].FraccionSUM[j] += FFraccionMasicaEspecie[FNin - 1][j] * GastoPonderacion *
FDeltaTime;
3016 d = dist - (double) n1;
3018 if(ResultadosMedios[i].TemperaturaGas || ResultadosMedios[i].FraccionMasicaEspecies) {
3019 GastoPonderacion = fabs(Interpola(FFlowMass[n1], FFlowMass[n2], 1., d));
3020 ResultadosMedios[i].GastoPonderacionSUM += GastoPonderacion;
3021 ResultadosMedios[i].PonderacionSUM += GastoPonderacion *
FDeltaTime;
3024 if(ResultadosMedios[i].TemperaturaGas)
3025 Vble = Interpola(FTemperature[n1], FTemperature[n2], 1., d);
3026 ResultadosMedios[i].TemperaturaGasSUM += Vble * GastoPonderacion;
3027 if(ResultadosMedios[i].Pressure) {
3028 Vble = Interpola(FPresion0[n1], FPresion0[n2], 1., d);
3029 ResultadosMedios[i].PresionSUM += Vble *
FDeltaTime;
3031 if(ResultadosMedios[i].Velocity) {
3032 Vble = Interpola(FVelocidad0[n1], FVelocidad0[n2], 1., d);
3033 ResultadosMedios[i].VelocidadSUM += Vble * __cons::ARef *
FDeltaTime;
3035 if(ResultadosMedios[i].Massflow) {
3036 Vble = Interpola(FFlowMass[n1], FFlowMass[n2], 1., d);
3037 ResultadosMedios[i].GastoSUM += Vble *
FDeltaTime;
3039 if(ResultadosMedios[i].TemperaturaInternaPared) {
3040 Vble = Interpola(FTPTubo[0][n1], FTPTubo[0][n2], 2., d);
3041 ResultadosMedios[i].TemperaturaInternaParedSUM += Vble *
FDeltaTime;
3043 if(ResultadosMedios[i].TemperaturaIntermediaPared) {
3044 Vble = Interpola(FTPTubo[1][n1], FTPTubo[1][n2], 2., d);
3045 ResultadosMedios[i].TemperaturaIntermediaParedSUM += Vble *
FDeltaTime;
3047 if(ResultadosMedios[i].TemperaturaExternaPared) {
3048 Vble = Interpola(FTPTubo[2][n1], FTPTubo[2][n2], 2., d);
3049 ResultadosMedios[i].TemperaturaExternaParedSUM += Vble *
FDeltaTime;
3051 if(ResultadosMedios[i].NITmedio) {
3052 double a = Interpola(FAsonido0[n1], FAsonido0[n2], 1., d);
3053 double v = Interpola(FVelocidad0[n1], FVelocidad0[n2], 1., d);
3054 double p = Interpola(FPresion0[n1], FPresion0[n2], 1., d);
3055 double diam = Interpola(FDiametroTubo[n1], FDiametroTubo[n2], 1., d);
3056 Gamma = Interpola(FGamma[n1], FGamma[n2], 1., d);
3057 Rmezcla = Interpola(FRMezcla[n1], FRMezcla[n2], 1., d);
3058 double nit = CalculaNIT(a, v, p, diam, Gamma, Rmezcla);
3059 ResultadosMedios[i].NITmedioSUM += nit *
FDeltaTime;
3061 if(ResultadosMedios[i].CoefPelInterior) {
3062 Vble = Interpola(Fhi[n1], Fhi[n2], 1., d);
3063 ResultadosMedios[i].CoefPelInteriorSUM += Vble *
FDeltaTime;
3065 if(ResultadosMedios[i].FraccionMasicaEspecies) {
3066 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
3067 Vble = Interpola(FFraccionMasicaEspecie[n1][j], FFraccionMasicaEspecie[n2][j], 1., d);
3068 ResultadosMedios[i].FraccionSUM[j] += Vble * GastoPonderacion *
FDeltaTime;
3075 if(Theta > FControlResMed * FAnguloTotalCiclo) {
3077 for(
int i = 0; i < FNumResMedios; i++) {
3078 if(ResultadosMedios[i].Pressure) {
3079 ResultadosMedios[i].PresionMED = ResultadosMedios[i].PresionSUM / FTiempoMedSUM;
3080 ResultadosMedios[i].PresionSUM = 0.;
3082 if(ResultadosMedios[i].TemperaturaGas) {
3083 if(ResultadosMedios[i].GastoPonderacionSUM != 0)
3084 ResultadosMedios[i].TemperaturaGasMED = ResultadosMedios[i].TemperaturaGasSUM / ResultadosMedios[i].GastoPonderacionSUM;
3086 ResultadosMedios[i].TemperaturaGasMED = 0;
3087 ResultadosMedios[i].TemperaturaGasSUM = 0.;
3088 ResultadosMedios[i].GastoPonderacionSUM = 0.;
3090 if(ResultadosMedios[i].Velocity) {
3091 ResultadosMedios[i].VelocidadMED = ResultadosMedios[i].VelocidadSUM / FTiempoMedSUM;
3092 ResultadosMedios[i].VelocidadSUM = 0.;
3094 if(ResultadosMedios[i].Massflow) {
3095 ResultadosMedios[i].GastoMED = ResultadosMedios[i].GastoSUM / FTiempoMedSUM;
3096 ResultadosMedios[i].GastoSUM = 0.;
3098 if(ResultadosMedios[i].TemperaturaInternaPared) {
3099 ResultadosMedios[i].TemperaturaInternaParedMED = ResultadosMedios[i].TemperaturaInternaParedSUM / FTiempoMedSUM;
3100 ResultadosMedios[i].TemperaturaInternaParedSUM = 0.;
3102 if(ResultadosMedios[i].TemperaturaIntermediaPared) {
3103 ResultadosMedios[i].TemperaturaIntermediaParedMED = ResultadosMedios[i].TemperaturaIntermediaParedSUM / FTiempoMedSUM;
3104 ResultadosMedios[i].TemperaturaIntermediaParedSUM = 0.;
3106 if(ResultadosMedios[i].TemperaturaExternaPared) {
3107 ResultadosMedios[i].TemperaturaExternaParedMED = ResultadosMedios[i].TemperaturaExternaParedSUM / FTiempoMedSUM;
3108 ResultadosMedios[i].TemperaturaExternaParedSUM = 0.;
3110 if(ResultadosMedios[i].NITmedio) {
3111 ResultadosMedios[i].NITmedioMED = ResultadosMedios[i].NITmedioSUM / FTiempoMedSUM;
3112 ResultadosMedios[i].NITmedioSUM = 0.;
3114 if(ResultadosMedios[i].CoefPelInterior) {
3115 ResultadosMedios[i].CoefPelInteriorMED = ResultadosMedios[i].CoefPelInteriorSUM / FTiempoMedSUM;
3116 ResultadosMedios[i].CoefPelInteriorSUM = 0.;
3118 if(ResultadosMedios[i].FraccionMasicaEspecies) {
3119 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
3120 if(DoubEqZero(ResultadosMedios[i].PonderacionSUM)) {
3121 ResultadosMedios[i].FraccionMED[j] = 0.;
3123 ResultadosMedios[i].FraccionMED[j] = ResultadosMedios[i].FraccionSUM[j] / ResultadosMedios[i].PonderacionSUM;
3125 ResultadosMedios[i].FraccionSUM[j] = 0.;
3127 ResultadosMedios[i].PonderacionSUM = 0.;
3131 FControlResMed = FControlResMed + 1.;
3134 }
catch(exception & N) {
3135 std::cout <<
"ERROR: TTubo::CalculaResultadosMedios en el tubo: " << FNumeroTubo << std::endl;
3136 std::cout <<
"Tipo de error: " << N.what() << std::endl;
3146 double dist = 0., d = 0., ason = 0., vel = 0., Aa = 0., ason1 = 0., vel1 = 0., Aa1 = 0., ason2 = 0., vel2 = 0.,
3153 if(FNumResInstant != 0) {
3154 for(
int i = 0; i < FNumResInstant; i++) {
3155 dist = ResultInstantaneos[i].Distancia / FXref;
3156 n1 = (int) floor(dist);
3157 if(n1 >= FNin - 1) {
3158 if(ResultInstantaneos[i].Pressure)
3159 ResultInstantaneos[i].PresionINS = FPresion0[FNin - 1];
3160 if(ResultInstantaneos[i].Velocity)
3161 ResultInstantaneos[i].VelocidadINS = FVelocidad0[FNin - 1] * __cons::ARef;
3162 if(ResultInstantaneos[i].TemperaturaGas) {
3163 double temp = __units::KTodegC(
pow2(FAsonido0[FNin - 1] * __cons::ARef) / (FGamma[FNin - 1] * FRMezcla[FNin - 1]));
3164 ResultInstantaneos[i].TemperaturaGasINS = temp;
3166 if(ResultInstantaneos[i].FlujoMasico) {
3167 ResultInstantaneos[i].FlujoMasicoINS = FFlowMass[FNin - 1];
3169 if(ResultInstantaneos[i].VelocidadDerecha || ResultInstantaneos[i].VelocidadIzquierda
3170 || ResultInstantaneos[i].PresionDerecha || ResultInstantaneos[i].PresionIzquierda) {
3171 ason = FAsonidoDim[FNin - 1];
3172 vel = FGamma1[FNin - 1] / 2 * FVelocidadDim[FNin - 1];
3173 Aa = ason / pow(FPresion0[FNin - 1], FGamma5[FNin - 1]);
3175 if(ResultInstantaneos[i].VelocidadDerecha) {
3176 double VelDer = ((ason + vel) - Aa) / FGamma1[FNin - 1];
3177 ResultInstantaneos[i].VelocidadDerechaINS = VelDer;
3179 if(ResultInstantaneos[i].VelocidadIzquierda) {
3180 double VelIzq = -((ason - vel) - Aa) / FGamma1[FNin - 1];
3181 ResultInstantaneos[i].VelocidadIzquierdaINS = VelIzq;
3183 if(ResultInstantaneos[i].PresionDerecha) {
3184 double PreDer = pow(((ason + vel) / Aa + 1) / 2., FGamma4[FNin - 1]);
3185 ResultInstantaneos[i].PresionDerechaINS = PreDer;
3187 if(ResultInstantaneos[i].PresionIzquierda) {
3188 double PreIzq = pow(((ason - vel) / Aa + 1) / 2., FGamma4[FNin - 1]);
3189 ResultInstantaneos[i].PresionIzquierdaINS = PreIzq;
3191 if(ResultInstantaneos[i].NIT) {
3192 double nit = CalculaNIT(FAsonido0[FNin - 1], FVelocidad0[FNin - 1], FPresion0[FNin - 1], FDiametroTubo[FNin - 1],
3193 FGamma[FNin - 1], FRMezcla[FNin - 1]);
3194 ResultInstantaneos[i].NITINS = nit;
3196 if(ResultInstantaneos[i].TemperaturaInternaPared)
3197 ResultInstantaneos[i].TemperaturaInternaParedINS = FTPTubo[0][FNin - 1];
3198 if(ResultInstantaneos[i].TemperaturaIntermediaPared)
3199 ResultInstantaneos[i].TemperaturaIntermediaParedINS = FTPTubo[1][FNin - 1];
3200 if(ResultInstantaneos[i].TemperaturaExternaPared)
3201 ResultInstantaneos[i].TemperaturaExternaParedINS = FTPTubo[2][FNin - 1];
3202 if(ResultInstantaneos[i].CoefPelInterior)
3203 ResultInstantaneos[i].CoefPelInteriorINS = Fhi[FNin - 1];
3204 if(ResultInstantaneos[i].FraccionMasicaEspecies) {
3205 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
3206 ResultInstantaneos[i].FraccionINS[j] = FFraccionMasicaEspecie[FNin - 1][j];
3209 if(ResultInstantaneos[i].Gamma)
3210 ResultInstantaneos[i].GammaINS = FGamma[FNin - 1];
3213 d = dist - (double) n1;
3214 if(ResultInstantaneos[i].Pressure) {
3215 double pres = Interpola(FPresion0[n1], FPresion0[n2], 1., d);
3216 ResultInstantaneos[i].PresionINS = pres;
3218 if(ResultInstantaneos[i].Velocity) {
3219 double vel = Interpola(FVelocidadDim[n1], FVelocidadDim[n2], 1., d);
3220 ResultInstantaneos[i].VelocidadINS = vel;
3222 if(ResultInstantaneos[i].TemperaturaGas) {
3223 double temp1 = FTemperature[n1];
3224 double temp2 = FTemperature[n2];
3225 double temp = Interpola(temp1, temp2, 1., d);
3226 ResultInstantaneos[i].TemperaturaGasINS = __units::KTodegC(temp);
3228 if(ResultInstantaneos[i].FlujoMasico) {
3231 double gto1 = FFlowMass[n1];
3234 double gto2 = FFlowMass[n2];
3235 double gto = Interpola(gto1, gto2, 1., d);
3236 ResultInstantaneos[i].FlujoMasicoINS = gto;
3238 if(ResultInstantaneos[i].VelocidadDerecha || ResultInstantaneos[i].VelocidadIzquierda
3239 || ResultInstantaneos[i].PresionDerecha || ResultInstantaneos[i].PresionIzquierda) {
3240 ason1 = FAsonidoDim[n1];
3241 vel1 = FGamma1[n1] / 2 * FVelocidadDim[n1];
3242 Aa1 = ason1 / pow(FPresion0[n1], FGamma5[n1]);
3243 ason2 = FAsonidoDim[n2];
3244 vel2 = FGamma1[n2] / 2 * FVelocidadDim[n2];
3245 Aa2 = ason2 / pow(FPresion0[n2], FGamma5[n2]);
3247 if(ResultInstantaneos[i].VelocidadDerecha) {
3248 double VelDer1 = ((ason1 + vel1) - Aa1) / FGamma1[n1];
3249 double VelDer2 = ((ason2 + vel2) - Aa2) / FGamma1[n2];
3250 double VelDer = Interpola(VelDer1, VelDer2, 1., d);
3251 ResultInstantaneos[i].VelocidadDerechaINS = VelDer;
3253 if(ResultInstantaneos[i].VelocidadIzquierda) {
3254 double VelIzq1 = -((ason1 - vel1) - Aa1) / FGamma1[n1];
3255 double VelIzq2 = -((ason2 - vel2) - Aa2) / FGamma1[n2];
3256 double VelIzq = Interpola(VelIzq1, VelIzq2, 1., d);
3257 ResultInstantaneos[i].VelocidadIzquierdaINS = VelIzq;
3259 if(ResultInstantaneos[i].PresionDerecha) {
3260 double PreDer1 = pow(((ason1 + vel1) / Aa1 + 1) / 2., FGamma4[n1]);
3261 double PreDer2 = pow(((ason2 + vel2) / Aa2 + 1) / 2., FGamma4[n2]);
3262 double PreDer = Interpola(PreDer1, PreDer2, 1., d);
3263 ResultInstantaneos[i].PresionDerechaINS = PreDer;
3265 if(ResultInstantaneos[i].PresionIzquierda) {
3266 double PreIzq1 = pow(((ason1 - vel1) / Aa1 + 1) / 2., FGamma4[n1]);
3267 double PreIzq2 = pow(((ason2 - vel2) / Aa2 + 1) / 2., FGamma4[n2]);
3268 double PreIzq = Interpola(PreIzq1, PreIzq2, 1., d);
3269 ResultInstantaneos[i].PresionIzquierdaINS = PreIzq;
3271 if(ResultInstantaneos[i].NIT) {
3272 double nit1 = CalculaNIT(FAsonido0[n1], FVelocidad0[n1], FPresion0[n1], FDiametroTubo[n1], FGamma[n1], FRMezcla[n1]);
3273 double nit2 = CalculaNIT(FAsonido0[n2], FVelocidad0[n2], FPresion0[n2], FDiametroTubo[n2], FGamma[n2], FRMezcla[n2]);
3274 double nit = Interpola(nit1, nit2, 1., d);
3275 ResultInstantaneos[i].NITINS = nit;
3277 if(ResultInstantaneos[i].TemperaturaInternaPared) {
3278 double TP = Interpola(FTPTubo[0][n1], FTPTubo[0][n2], 1., d);
3279 ResultInstantaneos[i].TemperaturaInternaParedINS = TP;
3281 if(ResultInstantaneos[i].TemperaturaIntermediaPared) {
3282 double TP = Interpola(FTPTubo[1][n1], FTPTubo[1][n2], 1., d);
3283 ResultInstantaneos[i].TemperaturaIntermediaParedINS = TP;
3285 if(ResultInstantaneos[i].TemperaturaExternaPared) {
3286 double TP = Interpola(FTPTubo[2][n1], FTPTubo[2][n2], 1., d);
3287 ResultInstantaneos[i].TemperaturaExternaParedINS = TP;
3289 if(ResultInstantaneos[i].CoefPelInterior) {
3290 double hi = Interpola(Fhi[n1], Fhi[n2], 1., d);
3291 ResultInstantaneos[i].CoefPelInteriorINS = hi;
3293 if(ResultInstantaneos[i].FraccionMasicaEspecies) {
3294 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
3295 double Fraccion = Interpola(FFraccionMasicaEspecie[n1][j], FFraccionMasicaEspecie[n2][j], 1., d);
3296 ResultInstantaneos[i].FraccionINS[j] = Fraccion;
3299 if(ResultInstantaneos[i].Gamma) {
3300 double gamma = Interpola(FGamma[n1], FGamma[n2], 1., d);
3301 ResultInstantaneos[i].GammaINS = gamma;
3308 }
catch(exception & N) {
3309 std::cout <<
"ERROR: TTubo::CalculaResultadosInstantaneos en el tubo: " << FNumeroTubo << std::endl;
3310 std::cout <<
"Tipo de error: " << N.what() << std::endl;
3319 double TTubo::CalculaNIT(
double a,
double v,
double p,
double d,
double Gamma,
double R) {
3324 double tem0 = 0., pre0 = 0., area = 0., gto = 0., nit = 0., V = 0., A = 0., V2 = 0., A2 = 0.;
3326 kp = Gamma * R / (Gamma - 1);
3327 V = v * __cons::ARef;
3329 A = a * __cons::ARef;
3331 tem0 = A2 / (Gamma * R) + V2 / 2. / kp;
3332 pre0 = __units::BarToPa(p) * pow((1 + V2 / 2. / kp / A2), (kp / 287.));
3333 area = __geom::Circle_area(d);
3334 gto = Gamma * __units::BarToPa(p) * area * V / A2;
3335 nit = gto * kp * tem0 * (1 - pow(pre0 / 100000., (-287. / kp)));
3338 }
catch(exception & N) {
3339 std::cout <<
"ERROR: TTubo::CalculaNIT en el tubo: " << FNumeroTubo << std::endl;
3340 std::cout <<
"Tipo de error: " << N.what() << std::endl;
3350 double AmbientTemperature) {
3355 double dtem, temed, rhog, viscext, Re, Pr, cond, viscpared, vel, n1, n2, L;
3357 if(FTipoCalcTempPared != nmTempConstante && FCoefAjusTC != 0) {
3359 if(FTipoTransCal == nmTuboAdmision || FTipoTransCal == nmTuboEscape) {
3360 for(
int i = 0; i < FNin; i++) {
3361 dtem = fabs(__units::degCToK(FTPTubo[2][i]) - FTExt);
3362 temed = (__units::degCToK(FTPTubo[2][i]) + FTExt) / 2.;
3367 switch(FTipRefrig) {
3370 rhog = __units::BarToPa(AmbientPressure) / __R::Air / temed;
3372 viscext = 1.4615e-6 *
pow150(temed) / (temed + 110.4);
3374 cond = (-8.39061e-09 * temed + 7.05256e-05) * temed + 6.51528e-03;
3375 Re = rhog * 1. * (FDiametroTubo[i] + 2 * (FEspesorIntPrin + FEspesorPrin + FEspesorExtPrin)) / viscext;
3385 viscext = ((-2.632351E-09 * temed + 2.737629E-06) * temed - 9.530709E-04) * temed + 1.114642E-01;
3386 Pr = ((-2.022269E-05 * temed + 2.106518E-02) * temed - 7.340298E+00) * temed + 8.581110E+02;
3387 cond = ((9.496332E-09 * temed - 1.707697E-05) * temed + 9.183462E-03) * temed - 8.626578E-01;
3389 Re = rhog * 1. * (FDiametroTubo[i] + 2 * (FEspesorIntPrin + FEspesorPrin + FEspesorExtPrin)) / viscext;
3392 std::cout <<
"WARNING: Tipo de refrigeracion mal definida en el tubo: " << FNumeroTubo << std::endl;
3395 if((2e4 < Re) && (Re < 4e5)) {
3402 Fhe[i] = 0.3 + 0.62 * sqrt(Re) * cbrt(Pr) /
pow025(1 + pow(0.4 / Pr, 0.666666)) * pow(1 + pow(Re / 282000, n1),
3403 n2) * cond / FDiametroTubo[i];
3407 Fhe[i] = Fhe[i] + __cons::Sigma * FEmisividad * (
pow4(__units::degCToK(FTPTubo[2][i])) -
pow4(FTExt)) / dtem;
3409 Fhe[i] = Fhe[i] * FCoefExt;
3412 temed = __units::degCToK(Engine[0]->getTempRefrigerante());
3414 viscext = ((-2.632351E-09 * temed + 2.737629E-06) * temed - 9.530709E-04) * temed + 1.114642E-01;
3415 L = 1.5 * Engine[0]->getGeometria().Diametro;
3417 vel = 5.64268e-7 * Engine[0]->getRegimen() / 60. * Engine[0]->getParPotMax() / Engine[0]->getGeometria().NCilin / pow(L,
3419 Re = rhog * fabs(vel) * L / 2.3 / viscext;
3420 Pr = ((-2.022269E-05 * temed + 2.106518E-02) * temed - 7.340298E+00) * temed + 8.581110E+02;
3421 cond = ((9.496332E-09 * temed - 1.707697E-05) * temed + 9.183462E-03) * temed - 8.626578E-01;
3423 for(
int i = 0; i < FNin; i++) {
3427 if(FTPTubo[2][i] > 100.) {
3428 viscpared = 0.000282;
3430 double Tp2 = __units::degCToK(FTPTubo[2][i]);
3431 viscpared = ((-2.632351E-09 * Tp2 + 2.737629E-06) * Tp2 - 9.530709E-04) * Tp2 + 1.114642E-01;
3434 Fhe[i] = 0.027 * (1 + 24.2 / pow(2.3, 0.7) /
pow025(Re)) * pow(Re, 0.8) * cbrt(Pr) * pow(viscext / viscpared,
3435 0.14) * cond / (L / 2.3);
3440 }
catch(exception & N) {
3441 std::cout <<
"ERROR: TTubo::CalculaCoeficientePeliculaExterior en el tubo: " << FNumeroTubo << std::endl;
3442 std::cout <<
"Tipo de error: " << N.what() << std::endl;
3455 double Dext, Dint, DIntPrin , Text, Tint, Rcond, Rrad, UnionEspes, UnionConduct, Cap;
3458 if(FTipoCalcTempPared != nmTempConstante && FCoefAjusTC != 0 && !FConcentrico) {
3459 for(
int i = 0; i < FNin; i++) {
3462 Dint = FDiametroTubo[i];
3463 FResistRadInt[i] = 0;
3464 FResistRadExt[i] = 0;
3465 for(
int j = 0; j < FNumCapas; j++) {
3466 Dext = Dint + 2 * FCapa[j].Espesor;
3467 if(FCapa[j].EsPrincipal ==
false) {
3468 Rcond = log(Dext / Dint) / __cons::Pi_x_2 / FCapa[j].Conductividad / FXref;
3471 if(FCapa[j].EsFluida) {
3473 Text = __units::degCToK(FTPTubo[1][i]);
3474 Tint = __units::degCToK(FTPTubo[0][i]);
3477 Rrad = (1 / FCapa[j].EmisividadInterior + Dint / Dext * (1 / FCapa[j].EmisividadExterior - 1)) /
3478 (__cons::Sigma * __cons::Pi * Dint * FXref) * (Tint - Text) / (
pow4(
3479 Tint) -
pow4(Text));
3480 FResistRadInt[i] += 1 / (1 / Rcond + 1 / Rrad);
3482 FResistRadInt[i] += Rcond;
3486 FResistRadInt[i] += Rcond;
3491 if(FCapa[j].EsFluida) {
3493 Text = __units::degCToK(FTPTubo[2][i]);
3494 Tint = __units::degCToK(FTPTubo[1][i]);
3496 Rrad = (1 / FCapa[j].EmisividadInterior + Dint / Dext * (1 / FCapa[j].EmisividadExterior - 1)) /
3497 (__cons::Sigma * __cons::Pi * Dint * FXref) * (Tint - Text) / (
pow4(
3498 Tint) -
pow4(Text));
3499 FResistRadExt[i] += 1 / (1 / Rcond + 1 / Rrad);
3501 FResistRadExt[i] += Rcond;
3505 FResistRadExt[i] += Rcond;
3510 FResistRadInt[i] += log((Dint + FCapa[j].Espesor) / Dint) / __cons::Pi_x_2 / FCapa[j].Conductividad / FXref;
3511 FResistRadExt[i] += log(Dext / (Dint + FCapa[j].Espesor)) / __cons::Pi_x_2 / FCapa[j].Conductividad / FXref;
3518 FResistAxiAnt[i] = 0;
3519 FResistAxiPos[i] = 0;
3520 DIntPrin = FDiametroTubo[i] + 2 * FEspesorIntPrin;
3522 if(BC[FNodoIzq - 1]->getTipoCC() == nmPipesConnection) {
3523 UnionEspes = dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoIzq - 1])->getEspesor();
3524 UnionConduct = dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoIzq - 1])->getConductividad();
3525 if(UnionConduct > 0) {
3526 FResistAxiAnt[i] = UnionEspes / UnionConduct / (__cons::Pi * (DIntPrin + FEspesorPrin) * FEspesorPrin);
3529 FResistAxiPos[i] = FXref / FConductPrin / (__cons::Pi * (DIntPrin + FEspesorPrin) * FEspesorPrin);
3531 #ifdef ParticulateFilter 3532 else if(BC[FNodoIzq - 1]->getTipoCC()
3533 == nmPipeToPlenumConnection) {
3534 if(FHayDPFNodoIzq) {
3535 if(FTipoCanal[0] == 0) {
3537 FResistAxiAnt[i] = FDPFEntradaTubo->getAjustRAxAnt
3539 (Pi * FDPFEntradaTubo->getConductividadMetal
3540 () * (FDPFEntradaTubo->getDiametroEfect()
3541 + 2 * FDPFEntradaTubo->getEspesorAislante
3542 () + 2 * FDPFEntradaTubo->getEspesorAire()
3543 + FDPFEntradaTubo->getEspesorMetal())
3544 * FDPFEntradaTubo->getEspesorMetal());
3546 FResistAxiAnt[i] = FDPFEntradaTubo->getAjustRAxPos
3548 (Pi * FDPFEntradaTubo->getConductividadMetal
3549 () * (FDPFEntradaTubo->getDiametroEfect()
3550 + 2 * FDPFEntradaTubo->getEspesorAislante
3551 () + 2 * FDPFEntradaTubo->getEspesorAire()
3552 + FDPFEntradaTubo->getEspesorMetal())
3553 * FDPFEntradaTubo->getEspesorMetal());
3558 else if(i == FNin - 1) {
3559 FResistAxiAnt[i] = FXref / FConductPrin / (__cons::Pi * (DIntPrin + FEspesorPrin) * FEspesorPrin);
3560 if(BC[FNodoDer - 1]->getTipoCC() == nmPipesConnection) {
3561 UnionEspes = dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoDer - 1])->getEspesor();
3562 UnionConduct = dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoDer - 1])->getConductividad();
3563 if(UnionConduct > 0) {
3564 FResistAxiPos[i] = UnionEspes / UnionConduct / (__cons::Pi * (DIntPrin + FEspesorPrin) * FEspesorPrin);
3567 #ifdef ParticulateFilter 3568 else if(BC[FNodoDer - 1]->getTipoCC()
3569 == nmPipeToPlenumConnection) {
3570 if(FHayDPFNodoDer) {
3571 if(FTipoCanal[1] == 0) {
3574 = FDPFSalidaTubo->getAjustRAxAnt() /
3575 (Pi * FDPFSalidaTubo->getConductividadMetal
3576 () * (FDPFSalidaTubo->getDiametroEfect()
3578 FDPFSalidaTubo->getEspesorAislante()
3579 + 2 * FDPFSalidaTubo->getEspesorAire()
3580 + FDPFSalidaTubo->getEspesorMetal())
3581 * FDPFSalidaTubo->getEspesorMetal());
3584 = FDPFSalidaTubo->getAjustRAxPos() /
3585 (Pi * FDPFSalidaTubo->getConductividadMetal
3586 () * (FDPFSalidaTubo->getDiametroEfect()
3588 FDPFSalidaTubo->getEspesorAislante()
3589 + 2 * FDPFSalidaTubo->getEspesorAire()
3590 + FDPFSalidaTubo->getEspesorMetal())
3591 * FDPFSalidaTubo->getEspesorMetal());
3597 FResistAxiAnt[i] = FXref / FConductPrin / (__cons::Pi * (DIntPrin + FEspesorPrin) * FEspesorPrin);
3598 FResistAxiPos[i] = FXref / FConductPrin / (__cons::Pi * (DIntPrin + FEspesorPrin) * FEspesorPrin);
3603 Dint = FDiametroTubo[i];
3607 for(
int j = 0; j < FNumCapas; j++) {
3608 Dext = Dint + 2 * FCapa[j].Espesor;
3609 if(FCapa[j].EsPrincipal ==
false) {
3612 if(!FCapa[j].EsFluida) {
3614 Cap = FCapa[j].Density * FCapa[j].CalorEspecifico * __geom::Ring_area(Dint, Dext) * FXref;
3619 if(!FCapa[j].EsFluida) {
3621 Cap = FCapa[j].Density * FCapa[j].CalorEspecifico * __geom::Ring_area(Dint, Dext) * FXref;
3627 FCapInt[i] += FDensidadPrin * FCalEspPrin * __geom::Ring_area(DIntPrin, DIntPrin + 0.5 * FEspesorPrin) * FXref;
3628 FCapMed[i] = FDensidadPrin * FCalEspPrin * __geom::Ring_area(DIntPrin + 0.5 * FEspesorPrin,
3629 DIntPrin + 1.5 * FEspesorPrin) * FXref;
3630 FCapExt[i] += FDensidadPrin * FCalEspPrin * __geom::Ring_area(DIntPrin + 1.5 * FEspesorPrin,
3631 DIntPrin + 2 * FEspesorPrin) * FXref;
3639 }
catch(exception & N) {
3640 std::cout <<
"ERROR: TTubo::CalculaResistenciasdePared en el tubo: " << FNumeroTubo << std::endl;
3641 std::cout <<
"Tipo de error: " << N.what() << std::endl;
3654 double Tg = 0., cesp = 0., viscgas = 0., cond = 0., viscpared = 0.;
3656 if(FCoefAjusTC != 0) {
3657 for(
int i = 0; i < FNin; i++) {
3659 Tg = FTemperature[i];
3660 cesp = FGamma[i] * FRMezcla[i] / (FGamma[i] - 1);
3661 viscgas = 1.4615e-6 *
pow150(Tg) / (Tg + 110.4);
3662 cond = viscgas * cesp / 0.709;
3663 if(DoubEqZero(FVelPro[i])) {
3666 FRe[i] = Frho[i] * FVelPro[i] * FDiametroTubo[i] / viscgas;
3670 switch(FTipoTransCal) {
3671 case nmTuboAdmision:
3673 Fhi[i] = 0.0694 *
pow075(FRe[i]) * cond / FDiametroTubo[i];
3675 case nmPipaAdmision:
3677 Fhi[i] = 0.0694 *
pow075(FRe[i]) * cond / FDiametroTubo[i];
3681 Fhi[i] = 1.6 * pow(FRe[i], 0.4) * cond / FDiametroTubo[i] * FCoefTurbulencia[i];
3685 viscpared = 1.4615e-6 *
pow150(__units::degCToK(FTPTubo[0][i])) / (__units::degCToK(FTPTubo[0][i]) + 110.4);
3686 Fhi[i] = 0.1 * pow(FRe[i], 0.8) * 0.709 * pow(viscgas / viscpared, 0.14) * cond / FDiametroTubo[i];
3689 std::cout <<
"WARNING: Transmision de calor mal definida en el tubo: " << FNumeroTubo << std::endl;
3692 }
else if(FCoefAjusFric != 0) {
3693 for(
int i = 0; i < FNin; i++) {
3694 Tg = FTemperature[i];
3695 viscgas = 1.4615e-6 *
pow150(Tg) / (Tg + 110.4);
3696 FRe[i] = Frho[i] * FVelPro[i] * FDiametroTubo[i] / viscgas;
3700 }
catch(exception & N) {
3701 std::cout <<
"ERROR: TTubo::CalculaCoeficientePeliculaInterior en el tubo: " << FNumeroTubo << std::endl;
3702 std::cout <<
"Tipo de error: " << N.what() << std::endl;
3713 double zzz = 0., czz = 0., cz1 = 0., uq1 = 0.;
3714 double DeltaTTPared = 0.;
3715 double Tpant0 = 0., Tpant1 = 0., Tpant2 = 0., Tpantant = 0., Tpantpos = 0., Text = 0., Ri = 0., Re = 0., ErrorTp = 0.;
3717 int extremo = 0, nodo = 0;
3724 for(
int i = 0; i < FNin; i++) {
3725 FTParedAnt[0][i] = FTPTubo[0][i];
3726 FTParedAnt[1][i] = FTPTubo[1][i];
3727 FTParedAnt[2][i] = FTPTubo[2][i];
3729 zzz = 0.013 / DeltaTTPared;
3730 czz = 2 / (zzz + 1);
3731 uq1 = fabs(FVelocidadDim[i]);
3732 if(FTipoTransCal == nmTuboEscape)
3736 FVelPro[i] = cz1 * uq1 + (1 - cz1) * FVelPro[i];
3739 if(FTipoCalcTempPared != nmTempConstante && FCoefAjusTC != 0) {
3741 if(Theta > FAnguloTotalCiclo) {
3742 FSUMTime += DeltaTTPared;
3744 for(
int i = 0; i < FNin; i++) {
3746 Tg = FTemperature[i];
3749 if(FTipoTransCal == nmTuboAdmision || FTipoTransCal == nmTuboEscape) {
3752 Text = __units::degCToK(Engine[0]->getTempRefrigerante());
3760 if(BC[FNodoIzq - 1]->getTipoCC() == nmPipesConnection) {
3761 if(dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoIzq - 1])->getConductividad() > 0) {
3762 if(BC[FNodoIzq - 1]->GetTuboExtremo(0).Pipe->getNumeroTubo() == FNumeroTubo) {
3764 if(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
3767 nodo = BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
3771 if(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
3774 nodo = BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
3777 Tpantant = __units::degCToK(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->GetTPTuboAnt(1, nodo));
3780 #ifdef ParticulateFilter 3781 else if(BC[FNodoIzq - 1]->getTipoCC() == nmPipeToPlenumConnection) {
3782 if(FHayDPFNodoIzq) {
3783 Tpantant = __units::degCToK(FDPFEntradaTubo->GetTSuperficie
3784 (FNodoDPFEntrada, 2));
3788 Tpantpos = __units::degCToK(FTParedAnt[1][i + 1]);
3789 }
else if(i == FNin - 1) {
3790 Tpantant = __units::degCToK(FTParedAnt[1][i - 1]);
3791 if(BC[FNodoDer - 1]->getTipoCC() == nmPipesConnection) {
3792 if(dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoDer - 1])->getConductividad() > 0) {
3793 if(BC[
getNodoDer() - 1]->GetTuboExtremo(0).Pipe->getNumeroTubo() == FNumeroTubo) {
3795 if(BC[FNodoDer - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
3798 nodo = BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
3802 if(BC[FNodoDer - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
3805 nodo = BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
3808 Tpantpos = __units::degCToK(BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->GetTPTuboAnt(1, nodo));
3811 #ifdef ParticulaFilter 3812 else if(BC[FNodoDer - 1]->getTipoCC() == nmPlenum) {
3813 if(FHayDPFNodoDer) {
3814 Tpantpos = __units::degCToK(FDPFSalidaTubo->GetTSuperficie
3815 (FNodoDPFSalida, 2));
3820 Tpantant = __units::degCToK(FTParedAnt[1][i - 1]);
3821 Tpantpos = __units::degCToK(FTParedAnt[1][i + 1]);
3824 Tpant0 = __units::degCToK(FTParedAnt[0][i]);
3825 Tpant1 = __units::degCToK(FTParedAnt[1][i]);
3826 Tpant2 = __units::degCToK(FTParedAnt[2][i]);
3829 if(DoubEqZero(Fhi[i])) {
3832 Ri = 1 / __cons::Pi / FDiametroTubo[i] / Fhi[i] / FCoefAjusTC / FXref;
3834 if(DoubEqZero(Fhe[i])) {
3837 Re = 1 / __cons::Pi / (FDiametroTubo[i] + 2 * (FEspesorIntPrin + FEspesorPrin + FEspesorExtPrin)) / Fhe[i] / FXref;
3840 FTPTubo[2][i] = DeltaTTPared / FCapExt[i] * (1 / FResistRadExt[i] * (Tpant1 - Tpant2) + 1 / Re *
3841 (Text - Tpant2)) + Tpant2;
3842 if(FResistAxiAnt[i] > 0. && FResistAxiPos[i] > 0.) {
3843 FTPTubo[1][i] = DeltaTTPared / FCapMed[i] * (1 / FResistRadInt[i] * (Tpant0 - Tpant1) + 1 / FResistRadExt[i] *
3844 (Tpant2 - Tpant1) + 1 / FResistAxiAnt[i] * (Tpantant - Tpant1) + 1 / FResistAxiPos[i] * (Tpantpos - Tpant1)) + Tpant1;
3845 }
else if(FResistAxiAnt[i] > 0.) {
3846 FTPTubo[1][i] = DeltaTTPared / FCapMed[i] * (1 / FResistRadInt[i] * (Tpant0 - Tpant1) + 1 / FResistRadExt[i] *
3847 (Tpant2 - Tpant1) + 1 / FResistAxiAnt[i] * (Tpantant - Tpant1)) + Tpant1;
3848 }
else if(FResistAxiPos[i] > 0.) {
3849 FTPTubo[1][i] = DeltaTTPared / FCapMed[i] * (1 / FResistRadInt[i] * (Tpant0 - Tpant1) + 1 / FResistRadExt[i] *
3850 (Tpant2 - Tpant1) + 1 / FResistAxiPos[i] * (Tpantpos - Tpant1)) + Tpant1;
3852 FTPTubo[1][i] = DeltaTTPared / FCapMed[i] * (1 / FResistRadInt[i] * (Tpant0 - Tpant1) + 1 / FResistRadExt[i] *
3853 (Tpant2 - Tpant1)) + Tpant1;
3855 FTPTubo[0][i] = DeltaTTPared / FCapInt[i] * (1 / Ri * (Tg - Tpant0) + 1 / FResistRadInt[i] *
3856 (Tpant1 - Tpant0)) + Tpant0;
3857 for(
int k = 0; k < 3; k++) {
3858 FTPTubo[k][i] = __units::KTodegC(FTPTubo[k][i]);
3862 if(FTipoCalcTempPared == nmVariableSinInerciaTermica
3863 || Theta / FAnguloTotalCiclo <= Engine[0]->getNumCiclosSinInerciaTermica()) {
3864 if(Fhi[i] != 0. && Theta > FAnguloTotalCiclo) {
3866 FSUMTPTuboPro[0][1][i] += 1 / (Ri + FResistRadInt[i]) * Tg * DeltaTTPared;
3868 FSUMTPTuboPro[1][1][i] += 1 / (Ri + FResistRadInt[i]) * DeltaTTPared;
3870 FSUMTPTuboPro[0][0][i] += 1 / Ri * Tg * DeltaTTPared;
3872 FSUMTPTuboPro[1][0][i] += 1 / Ri * DeltaTTPared;
3879 if(FCicloTubo != Engine[0]->getCiclo() && FSUMTime > 0.) {
3881 if((FTipoCalcTempPared == nmVariableSinInerciaTermica
3882 || Theta / FAnguloTotalCiclo <= Engine[0]->getNumCiclosSinInerciaTermica()) && Theta > FAnguloTotalCiclo + 1) {
3884 EsPrimeraVez =
true;
3888 for(
int i = 0; i < FNin; i++) {
3891 Tpantant = __units::degCToK(FTPTubo[1][i]);
3892 Tpantpos = __units::degCToK(FTPTubo[1][i]);
3894 if(BC[FNodoIzq - 1]->getTipoCC() == nmPipesConnection) {
3895 if(dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoIzq - 1])->getConductividad() > 0) {
3896 if(BC[FNodoIzq - 1]->GetTuboExtremo(0).Pipe->getNumeroTubo() == FNumeroTubo) {
3898 if(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
3901 nodo = BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
3905 if(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
3908 nodo = BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
3911 Tpantant = __units::degCToK(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->GetTPTuboAnt(1, nodo));
3914 #ifdef ParticulateFilter 3916 (BC[FNodoIzq - 1]->getTipoCC()
3917 == nmPipeToPlenumConnection) {
3918 if(FHayDPFNodoIzq) {
3920 __units::degCToK(FDPFEntradaTubo->GetTSuperficie
3921 (FNodoDPFEntrada, 2));
3925 Tpantpos = __units::degCToK(FTPTubo[1][i + 1]);
3926 }
else if(i == FNin - 1) {
3927 Tpantant = __units::degCToK(FTPTubo[1][i - 1]);
3928 if(BC[FNodoDer - 1]->getTipoCC() == nmPipesConnection) {
3929 if(dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoDer - 1])->getConductividad() > 0) {
3930 if(BC[FNodoDer - 1]->GetTuboExtremo(0).Pipe->getNumeroTubo() == FNumeroTubo) {
3932 if(BC[FNodoDer - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
3935 nodo = BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
3939 if(BC[FNodoDer - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
3942 nodo = BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
3945 Tpantpos = __units::degCToK(BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->GetTPTuboAnt(1, nodo));
3948 #ifdef ParticulateFilter 3950 (BC[FNodoDer - 1]->getTipoCC()
3951 == nmPipeToPlenumConnection) {
3952 if(FHayDPFNodoDer) {
3954 __units::degCToK(FDPFSalidaTubo->GetTSuperficie
3955 (FNodoDPFSalida, 2));
3960 Tpantant = __units::degCToK(FTPTubo[1][i - 1]);
3961 Tpantpos = __units::degCToK(FTPTubo[1][i + 1]);
3965 Tpant0 = __units::degCToK(FTPTubo[0][i]);
3966 Tpant1 = __units::degCToK(FTPTubo[1][i]);
3967 Tpant2 = __units::degCToK(FTPTubo[2][i]);
3970 FTPTubo[1][i] = (FSUMTime * 1 / (Re + FResistRadExt[i]) * Text + FSUMTPTuboPro[0][1][i]) / (FSUMTime * 1 /
3971 (Re + FResistRadExt[i]) + FSUMTPTuboPro[1][1][i]);
3973 if(FResistAxiAnt[i] > 0. && FResistAxiPos[i] > 0.) {
3974 FTPTubo[1][i] = (FSUMTime * (1 / (Re + FResistRadExt[i]) * Text + 1 / FResistAxiAnt[i] * Tpantant + 1 / FResistAxiPos[i]
3975 * Tpantpos) + FSUMTPTuboPro[0][1][i]) / (FSUMTime * (1 / (Re + FResistRadExt[i]) + 1 / FResistAxiAnt[i] + 1 /
3976 FResistAxiPos[i]) + FSUMTPTuboPro[1][1][i]);
3977 }
else if(FResistAxiAnt[i] > 0.) {
3978 FTPTubo[1][i] = (FSUMTime * (1 / (Re + FResistRadExt[i]) * Text + 1 / FResistAxiAnt[i] * Tpantant) +
3979 FSUMTPTuboPro[0][1][i]) / (FSUMTime * (1 / (Re + FResistRadExt[i]) + 1 / FResistAxiAnt[i]) + FSUMTPTuboPro[1][1][i]);
3980 }
else if(FResistAxiPos[i] > 0.) {
3981 FTPTubo[1][i] = (FSUMTime * (1 / (Re + FResistRadExt[i]) * Text + 1 / FResistAxiPos[i] * Tpantpos) +
3982 FSUMTPTuboPro[0][1][i]) / (FSUMTime * (1 / (Re + FResistRadExt[i]) + 1 / FResistAxiPos[i]) + FSUMTPTuboPro[1][1][i]);
3984 FTPTubo[1][i] = (FSUMTime * (1 / (Re + FResistRadExt[i]) * Text) + FSUMTPTuboPro[0][1][i]) / (FSUMTime * (1 /
3985 (Re + FResistRadExt[i])) + FSUMTPTuboPro[1][1][i]);
3988 FTPTubo[0][i] = ((FSUMTime * 1 / FResistRadInt[i] * Tpant1 + FSUMTPTuboPro[0][0][i]) /
3989 (FSUMTime * 1 / FResistRadInt[i] + FSUMTPTuboPro[1][0][i]));
3990 FTPTubo[2][i] = (FSUMTime * (1 / Re * Text + 1 / FResistRadExt[i] * Tpant1)) / (FSUMTime *
3991 (1 / Re + 1 / FResistRadExt[i]));
3992 if(ErrorTp < fabs(Tpant1 - FTPTubo[1][i])) {
3993 ErrorTp = fabs(Tpant1 - FTPTubo[1][i]);
3995 for(
int k = 0; k < 3; k++) {
3996 FTPTubo[k][i] = __units::KTodegC(FTPTubo[k][i]);
3999 EsPrimeraVez =
false;
4002 for(
int i = 0; i < FNin; i++) {
4003 for(
int j = 0; j < 2; j++) {
4004 for(
int k = 0; k < 3; k++) {
4005 FSUMTPTuboPro[j][k][i] = 0.;
4012 if(FTipoCalcTempPared != nmTempConstante && FCoefAjusTC != 0) {
4013 if(FCicloTubo != Engine[0]->getCiclo()) {
4015 FCicloTubo = Engine[0]->getCiclo();
4019 }
catch(exception & N) {
4020 std::cout <<
"ERROR: TTubo::CalculaTemperaturaPared en el tubo: " << FNumeroTubo << std::endl;
4021 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4032 double zzz = 0., czz = 0., cz1 = 0., uq1 = 0.;
4033 double DeltaTTPared = 0.;
4034 double Tpant0 = 0., Tpant1 = 0., Tpant2 = 0., Tpantant = 0., Tpantpos = 0., Text = 0., Ri = 0., Re = 0., ErrorTp = 0.;
4036 int extremo = 0, nodo = 0;
4044 for(
int i = 0; i < FNin; i++) {
4045 FTParedAnt[0][i] = FTPTubo[0][i];
4046 FTParedAnt[1][i] = FTPTubo[1][i];
4047 FTParedAnt[2][i] = FTPTubo[2][i];
4049 zzz = 0.013 / DeltaTTPared;
4050 czz = 2 / (zzz + 1);
4051 uq1 = fabs(FVelocidadDim[i]);
4052 if(FTipoTransCal == nmTuboEscape)
4056 FVelPro[i] = cz1 * uq1 + (1 - cz1) * FVelPro[i];
4059 if(FTipoCalcTempPared != nmTempConstante && FCoefAjusTC != 0) {
4060 FCicloActual =
FTime1 / FDuracionCiclo;
4061 if(
FTime1 > FDuracionCiclo) {
4062 FSUMTime += DeltaTTPared;
4064 for(
int i = 0; i < FNin; i++) {
4066 Tg = FTemperature[i];
4077 if(BC[FNodoIzq - 1]->getTipoCC() == nmPipesConnection) {
4078 if(dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoIzq - 1])->getConductividad() > 0) {
4079 if(BC[FNodoIzq - 1]->GetTuboExtremo(0).Pipe->getNumeroTubo() == FNumeroTubo) {
4081 if(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
4084 nodo = BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
4088 if(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
4091 nodo = BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
4094 Tpantant = __units::degCToK(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->GetTPTuboAnt(1, nodo));
4097 #ifdef ParticulateFilter 4098 else if(BC[FNodoIzq - 1]->getTipoCC() == nmPipeToPlenumConnection) {
4099 if(FHayDPFNodoIzq) {
4100 Tpantant = __units::degCToK(FDPFEntradaTubo->GetTSuperficie
4101 (FNodoDPFEntrada, 2));
4105 Tpantpos = __units::degCToK(FTParedAnt[1][i + 1]);
4106 }
else if(i == FNin - 1) {
4107 Tpantant = __units::degCToK(FTParedAnt[1][i - 1]);
4108 if(BC[FNodoDer - 1]->getTipoCC() == nmPipesConnection) {
4109 if(dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoDer - 1])->getConductividad() > 0) {
4110 if(BC[
getNodoDer() - 1]->GetTuboExtremo(0).Pipe->getNumeroTubo() == FNumeroTubo) {
4112 if(BC[FNodoDer - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
4115 nodo = BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
4119 if(BC[FNodoDer - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
4122 nodo = BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
4125 Tpantpos = __units::degCToK(BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->GetTPTuboAnt(1, nodo));
4128 #ifdef ParticulateFilter 4129 else if(BC[FNodoDer - 1]->getTipoCC() == nmPipeToPlenumConnection) {
4130 if(FHayDPFNodoDer) {
4131 Tpantpos = __units::degCToK(FDPFSalidaTubo->GetTSuperficie
4132 (FNodoDPFSalida, 2));
4137 Tpantant = __units::degCToK(FTParedAnt[1][i - 1]);
4138 Tpantpos = __units::degCToK(FTParedAnt[1][i + 1]);
4141 Tpant0 = __units::degCToK(FTParedAnt[0][i]);
4142 Tpant1 = __units::degCToK(FTParedAnt[1][i]);
4143 Tpant2 = __units::degCToK(FTParedAnt[2][i]);
4146 if(DoubEqZero(Fhi[i])) {
4149 Ri = 1 / __cons::Pi / FDiametroTubo[i] / Fhi[i] / FCoefAjusTC / FXref;
4151 if(DoubEqZero(Fhe[i])) {
4154 Re = 1 / __cons::Pi / (FDiametroTubo[i] + 2 * (FEspesorIntPrin + FEspesorPrin + FEspesorExtPrin)) / Fhe[i] / FXref;
4157 FTPTubo[2][i] = DeltaTTPared / FCapExt[i] * (1 / FResistRadExt[i] * (Tpant1 - Tpant2) + 1 / Re *
4158 (Text - Tpant2)) + Tpant2;
4159 if(FResistAxiAnt[i] > 0. && FResistAxiPos[i] > 0.) {
4160 FTPTubo[1][i] = DeltaTTPared / FCapMed[i] * (1 / FResistRadInt[i] * (Tpant0 - Tpant1) + 1 / FResistRadExt[i] *
4161 (Tpant2 - Tpant1) + 1 / FResistAxiAnt[i] * (Tpantant - Tpant1) + 1 / FResistAxiPos[i] * (Tpantpos - Tpant1)) + Tpant1;
4162 }
else if(FResistAxiAnt[i] > 0.) {
4163 FTPTubo[1][i] = DeltaTTPared / FCapMed[i] * (1 / FResistRadInt[i] * (Tpant0 - Tpant1) + 1 / FResistRadExt[i] *
4164 (Tpant2 - Tpant1) + 1 / FResistAxiAnt[i] * (Tpantant - Tpant1)) + Tpant1;
4165 }
else if(FResistAxiPos[i] > 0.) {
4166 FTPTubo[1][i] = DeltaTTPared / FCapMed[i] * (1 / FResistRadInt[i] * (Tpant0 - Tpant1) + 1 / FResistRadExt[i] *
4167 (Tpant2 - Tpant1) + 1 / FResistAxiPos[i] * (Tpantpos - Tpant1)) + Tpant1;
4169 FTPTubo[1][i] = DeltaTTPared / FCapMed[i] * (1 / FResistRadInt[i] * (Tpant0 - Tpant1) + 1 / FResistRadExt[i] *
4170 (Tpant2 - Tpant1)) + Tpant1;
4172 FTPTubo[0][i] = DeltaTTPared / FCapInt[i] * (1 / Ri * (Tg - Tpant0) + 1 / FResistRadInt[i] *
4173 (Tpant1 - Tpant0)) + Tpant0;
4174 for(
int k = 0; k < 3; k++) {
4175 FTPTubo[k][i] = __units::KTodegC(FTPTubo[k][i]);
4179 if(FTipoCalcTempPared == nmVariableSinInerciaTermica || FCicloActual <= FNumCiclosSinInerciaTermica) {
4180 if(
FTime1 > FDuracionCiclo) {
4182 FSUMTPTuboPro[0][1][i] += 1 / (Ri + FResistRadInt[i]) * Tg * DeltaTTPared;
4184 FSUMTPTuboPro[1][1][i] += 1 / (Ri + FResistRadInt[i]) * DeltaTTPared;
4186 FSUMTPTuboPro[0][0][i] += 1 / Ri * Tg * DeltaTTPared;
4188 FSUMTPTuboPro[1][0][i] += 1 / Ri * DeltaTTPared;
4195 if(FCicloTubo != FCicloActual && FSUMTime > 0. && FCicloActual > 1) {
4197 if((FTipoCalcTempPared == nmVariableSinInerciaTermica || FCicloActual <= FNumCiclosSinInerciaTermica)
4198 &&
FTime1 > FDuracionCiclo) {
4200 EsPrimeraVez =
true;
4204 for(
int i = 0; i < FNin; i++) {
4207 Tpantant = __units::degCToK(FTPTubo[1][i]);
4208 Tpantpos = __units::degCToK(FTPTubo[1][i]);
4210 if(BC[FNodoIzq - 1]->getTipoCC() == nmPipesConnection) {
4211 if(dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoIzq - 1])->getConductividad() > 0) {
4212 if(BC[FNodoIzq - 1]->GetTuboExtremo(0).Pipe->getNumeroTubo() == FNumeroTubo) {
4214 if(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
4217 nodo = BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
4221 if(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
4224 nodo = BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
4227 Tpantant = __units::degCToK(BC[FNodoIzq - 1]->GetTuboExtremo(extremo).Pipe->GetTPTuboAnt(1, nodo));
4230 #ifdef ParticulateFilter 4232 (BC[FNodoIzq - 1]->getTipoCC()
4233 == nmPipeToPlenumConnection) {
4234 if(FHayDPFNodoIzq) {
4236 __units::degCToK(FDPFEntradaTubo->GetTSuperficie
4237 (FNodoDPFEntrada, 2));
4241 Tpantpos = __units::degCToK(FTPTubo[1][i + 1]);
4242 }
else if(i == FNin - 1) {
4243 Tpantant = __units::degCToK(FTPTubo[1][i - 1]);
4244 if(BC[FNodoDer - 1]->getTipoCC() == nmPipesConnection) {
4245 if(dynamic_cast<TCCUnionEntreTubos*>(BC[FNodoDer - 1])->getConductividad() > 0) {
4246 if(BC[FNodoDer - 1]->GetTuboExtremo(0).Pipe->getNumeroTubo() == FNumeroTubo) {
4248 if(BC[FNodoDer - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
4251 nodo = BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
4255 if(BC[FNodoDer - 1]->GetTuboExtremo(extremo).TipoExtremo == nmLeft) {
4258 nodo = BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->getNin() - 1;
4261 Tpantpos = __units::degCToK(BC[FNodoDer - 1]->GetTuboExtremo(extremo).Pipe->GetTPTuboAnt(1, nodo));
4264 #ifdef ParticulateFilter 4266 (BC[FNodoDer - 1]->getTipoCC()
4267 == nmPipeToPlenumConnection) {
4268 if(FHayDPFNodoDer) {
4270 __units::degCToK(FDPFSalidaTubo->GetTSuperficie
4271 (FNodoDPFSalida, 2));
4276 Tpantant = __units::degCToK(FTPTubo[1][i - 1]);
4277 Tpantpos = __units::degCToK(FTPTubo[1][i + 1]);
4281 Tpant0 = __units::degCToK(FTPTubo[0][i]);
4282 Tpant1 = __units::degCToK(FTPTubo[1][i]);
4283 Tpant2 = __units::degCToK(FTPTubo[2][i]);
4286 FTPTubo[1][i] = (FSUMTime / (Re + FResistRadExt[i]) * Text + FSUMTPTuboPro[0][1][i]) / (FSUMTime /
4287 (Re + FResistRadExt[i]) + FSUMTPTuboPro[1][1][i]);
4289 if(FResistAxiAnt[i] > 0. && FResistAxiPos[i] > 0.) {
4290 FTPTubo[1][i] = (FSUMTime * (1 / (Re + FResistRadExt[i]) * Text + 1 / FResistAxiAnt[i] * Tpantant + 1 / FResistAxiPos[i]
4291 * Tpantpos) + FSUMTPTuboPro[0][1][i]) / (FSUMTime * (1 / (Re + FResistRadExt[i]) + 1 / FResistAxiAnt[i] + 1 /
4292 FResistAxiPos[i]) + FSUMTPTuboPro[1][1][i]);
4293 }
else if(FResistAxiAnt[i] > 0.) {
4294 FTPTubo[1][i] = (FSUMTime * (1 / (Re + FResistRadExt[i]) * Text + 1 / FResistAxiAnt[i] * Tpantant) +
4295 FSUMTPTuboPro[0][1][i]) / (FSUMTime * (1 / (Re + FResistRadExt[i]) + 1 / FResistAxiAnt[i]) + FSUMTPTuboPro[1][1][i]);
4296 }
else if(FResistAxiPos[i] > 0.) {
4297 FTPTubo[1][i] = (FSUMTime * (1 / (Re + FResistRadExt[i]) * Text + 1 / FResistAxiPos[i] * Tpantpos) +
4298 FSUMTPTuboPro[0][1][i]) / (FSUMTime * (1 / (Re + FResistRadExt[i]) + 1 / FResistAxiPos[i]) + FSUMTPTuboPro[1][1][i]);
4300 FTPTubo[1][i] = (FSUMTime * (1 / (Re + FResistRadExt[i]) * Text) + FSUMTPTuboPro[0][1][i]) / (FSUMTime * (1 /
4301 (Re + FResistRadExt[i])) + FSUMTPTuboPro[1][1][i]);
4304 FTPTubo[0][i] = ((FSUMTime / FResistRadInt[i] * Tpant1 + FSUMTPTuboPro[0][0][i]) / (FSUMTime / FResistRadInt[i] +
4305 FSUMTPTuboPro[1][0][i]));
4306 FTPTubo[2][i] = (FSUMTime * (1 / Re * Text + 1 / FResistRadExt[i] * Tpant1)) / (FSUMTime *
4307 (1 / Re + 1 / FResistRadExt[i]));
4308 if(ErrorTp < fabs(Tpant1 - FTPTubo[1][i])) {
4309 ErrorTp = fabs(Tpant1 - FTPTubo[1][i]);
4311 for(
int k = 0; k < 3; k++) {
4312 FTPTubo[k][i] = __units::KTodegC(FTPTubo[k][i]);
4315 EsPrimeraVez =
false;
4318 for(
int i = 0; i < FNin; i++) {
4319 for(
int j = 0; j < 2; j++) {
4320 for(
int k = 0; k < 3; k++) {
4321 FSUMTPTuboPro[j][k][i] = 0.;
4326 if(FCicloTubo != FCicloActual) {
4327 FCicloTubo = FCicloActual;
4332 }
catch(exception & N) {
4333 std::cout <<
"ERROR: TTubo::CalculaTemperaturaParedSinMotor en el tubo: " << FNumeroTubo << std::endl;
4334 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4347 if(FNumResMedios > 0) {
4348 printf(
"__________________________________\n");
4349 printf(
" AVERAGE RESULTS IN PIPES \n");
4350 printf(
"__________________________________\n\n\n");
4351 for(
int i = 0; i < FNumResMedios; i++) {
4352 std::cout <<
"Average results in point " << i <<
" place in pipe " << FNumeroTubo;
4353 std::cout <<
" at " << ResultadosMedios[i].Distancia <<
" metres from left end" << std::endl;
4354 if(ResultadosMedios[i].TemperaturaGas)
4355 std::cout <<
"Average Temperature = " << ResultadosMedios[i].TemperaturaGasMED <<
" C" << std::endl;
4356 if(ResultadosMedios[i].Pressure)
4357 std::cout <<
"Average Pressure = " << ResultadosMedios[i].PresionMED <<
" bar" << std::endl;
4358 if(ResultadosMedios[i].Velocity)
4359 std::cout <<
"Average Velocity = " << ResultadosMedios[i].VelocidadMED <<
" m/s" << std::endl;
4360 if(ResultadosMedios[i].Massflow)
4361 std::cout <<
"Average Massflow = " << ResultadosMedios[i].GastoMED <<
" kg/s" << std::endl;
4362 if(ResultadosMedios[i].NITmedio)
4363 std::cout <<
"NIT Medio = " << ResultadosMedios[i].NITmedioMED <<
" Watts" << std::endl << std::endl;
4368 }
catch(exception & N) {
4369 std::cout <<
"ERROR: TTubo::SalidaGeneralTubos en el tubo nº: " << FNumeroTubo << std::endl;
4370 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4379 inline double TTubo::Maximo(
double x,
double y) {
4380 return x > y ? x : y;
4386 inline double TTubo::Minimo(
double x,
double y) {
4387 return x > y ? y : x;
4400 }
catch(exception & N) {
4401 std::cout <<
"ERROR: TTubo::AjustaPaso en el tubo: " << FNumeroTubo << std::endl;
4402 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4416 if(FVelocidad0[0] <= 0) {
4417 BC[FNodoIzq - 1]->PutEntropia(FTuboCCNodoIzq,
4418 Interpola_Entropia(BC[FNodoIzq - 1]->GetTuboExtremo(FTuboCCNodoIzq).TipoExtremo, DeltaTiempo));
4420 BC[FNodoIzq - 1]->PutBeta(FTuboCCNodoIzq,
4423 if(FVelocidad0[FNin - 1] >= 0) {
4424 BC[FNodoDer - 1]->PutEntropia(FTuboCCNodoDer,
4425 Interpola_Entropia(BC[FNodoDer - 1]->GetTuboExtremo(FTuboCCNodoDer).TipoExtremo, DeltaTiempo));
4427 BC[FNodoDer - 1]->PutLanda(FTuboCCNodoDer,
4428 Interpola_Caracteristica(BC[FNodoDer - 1]->GetTuboExtremo(FTuboCCNodoDer).Entropia, -1, getNin() - 1, DeltaTiempo));
4430 }
catch(exception & N) {
4431 std::cout <<
"ERROR: TTubo::CalculaCaracteristicasExtremos en el tubo: " << FNumeroTubo << std::endl;
4432 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4450 if(TipoExtremoTubo == nmRight) {
4455 double dtdx = DeltaTiempo / FXref;
4457 double entropia = 0.;
4458 double velocidadp = 0.;
4459 double asonidop = 0.;
4461 if(DeltaTiempo < 1e-15 || DoubEqZero(FVelocidadDim[extremo])) {
4463 Calculo_Entropia(entropia, velocidadp, extremo, 0., signo, DeltaTiempo, indiceCC);
4466 int ind1 = ind + signo;
4468 stPathOrigin PathOrigin(FU0[0][ind], FU0[1][ind], FU0[0][ind1], FU0[1][ind1], dtdx, signo);
4470 double dist =
zbrent(PathOrigin, 0., 1., 1e-5);
4472 Calculo_Entropia(entropia, velocidadp, ind, dist, signo, DeltaTiempo, indiceCC);
4476 return entropia / __cons::ARef;
4478 }
catch(exception & N) {
4479 std::cout <<
"ERROR: TTubo::Interpola_Entropia: " << std::endl;
4480 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4489 void TTubo::Calculo_Entropia(
double& entropia,
double& velocidadp,
int ind,
double dist,
int signo,
double DeltaTiempo,
4495 int ind1 = ind + signo;
4497 double w0 = FU0[0][ind];
4498 double w1 = FU0[1][ind];
4499 double w2 = FU0[2][ind];
4500 double gammap = FGamma[ind];
4501 double Rmezclap = FRMezcla[ind];
4502 double diamep = FDiametroTubo[ind];
4504 double tptubop = FTPTubo[0][ind];
4505 double hip = Fhi[ind];
4506 double rhop = Frho[ind];
4507 double Rep = FRe[ind];
4509 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
4510 FFraccionMasicaCC[indiceCC][j] = FFraccionMasicaEspecie[ind][j];
4513 if(dist > 0. || dist < 1.) {
4514 w0 = Interpola(FU0[0][ind], FU0[0][ind1], 1., dist);
4515 w1 = Interpola(FU0[1][ind], FU0[1][ind1], 1., dist);
4516 w2 = Interpola(FU0[2][ind], FU0[2][ind1], 1., dist);
4517 gammap = Interpola(FGamma[ind], FGamma[ind1], 1., dist);
4518 Rmezclap = Interpola(FRMezcla[ind], FRMezcla[ind1], 1., dist);
4519 diamep = Interpola(FDiametroTubo[ind], FDiametroTubo[ind1], 1., dist);
4520 if(FCoefAjusTC != 0 || FCoefAjusFric != 0) {
4521 tptubop = Interpola(FTPTubo[0][ind], FTPTubo[0][ind1], 1., dist);
4522 hip = Interpola(Fhi[ind], Fhi[ind1], 1., dist);
4523 Rep = Interpola(FRe[ind], FRe[ind1], 1., dist);
4526 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
4527 FFraccionMasicaCC[indiceCC][j] = Interpola(FFraccionMasicaEspecie[ind][j], FFraccionMasicaEspecie[ind1][j], 1., dist);
4531 double gamma1p = __Gamma::G1(gammap);
4532 double gamma3p = __Gamma::G3(gammap);
4533 double gamma5p = __Gamma::G5(gammap);
4534 velocidadp = w1 / w0;
4535 rhop = w0 / __geom::Circle_area(diamep);
4536 double asonidop = sqrt(gammap * gamma1p * (w2 / w0 -
pow2(velocidadp) / 2));
4537 double presionp = __units::PaToBar((w2 -
pow2(w1) / 2. / w0) * gamma1p / __geom::Circle_area(diamep));
4538 double entropiap = asonidop / pow(presionp, gamma5p);
4539 entropia = entropiap;
4541 FAreaCC[indiceCC] = FArea[ind];
4542 FDensidadCC[indiceCC] = rhop;
4543 FVelocidadCC[indiceCC] = velocidadp;
4547 if(!DoubEqZero(FCoefAjusTC)) {
4550 double tgasp =
pow2(asonidop) / (gammap * Rmezclap);
4552 TransmisionCalor(tgasp, diamep, q, hip, rhop, tptubop);
4555 double dacal = gamma3p * entropiap * q * FCoefAjusTC * DeltaTiempo /
pow2(asonidop);
4562 if(!DoubEqZero(velocidadp) && !DoubEqZero(FCoefAjusFric)) {
4565 double velabs = fabs(velocidadp);
4566 Colebrook(FFriccion, diamep, f, Rep);
4567 double dafric = gamma1p * FCoefAjusFric * f * entropiap *
pow3(velabs) * DeltaTiempo / (diamep *
pow2(asonidop));
4573 catch(exception & N) {
4574 std::cout <<
"ERROR: TTubo::Calculo_Entropia " << FNumeroTubo << std::endl;
4575 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4589 double dtdx = DeltaTiempo / FXref;
4591 double caracteristica = 0.;
4592 double velocidadp = 0.;
4593 double asonidop = 0.;
4595 if(DeltaTiempo < 1e-15) {
4596 Calculo_Caracteristica(caracteristica, velocidadp, asonidop, ind, 0., signo, entropia, DeltaTiempo);
4599 dtdx = DeltaTiempo / FXref;
4601 int ind1 = ind + signo;
4603 stCharOrigin CharOrigin(FU0[0][ind], FU0[1][ind], FU0[2][ind], FU0[0][ind1], FU0[1][ind1], FU0[2][ind1], FGamma[ind],
4604 FGamma[ind1], dtdx, signo);
4606 double dist =
zbrent(CharOrigin, 0., 1., 1e-5);
4608 Calculo_Caracteristica(caracteristica, velocidadp, asonidop, ind, dist, signo, entropia, DeltaTiempo);
4610 return caracteristica / __cons::ARef;
4612 }
catch(exception & N) {
4613 std::cout <<
"ERROR: TTubo::Interpola_Caracteristica " << FNumeroTubo << std::endl;
4614 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4623 void TTubo::Calculo_Caracteristica(
double& caracteristica,
double& velocidadp,
double& asonidop,
int ind,
double dist,
4624 int signo,
double entropia,
double DeltaTiempo) {
4629 int ind1 = ind + signo;
4631 double w0 = FU0[0][ind];
4632 double w1 = FU0[1][ind];
4633 double w2 = FU0[2][ind];
4634 double gammap = FGamma[ind];
4635 double Rmezclap = FRMezcla[ind];
4636 double diamep = FDiametroTubo[ind];
4638 double tptubop = FTPTubo[0][ind];
4639 double hip = Fhi[ind];
4640 double rhop = Frho[ind];
4641 double Rep = FRe[ind];
4643 if(dist < 1. && dist > 0) {
4645 w0 = Interpola(FU0[0][ind], FU0[0][ind1], 1., dist);
4646 w1 = Interpola(FU0[1][ind], FU0[1][ind1], 1., dist);
4647 w2 = Interpola(FU0[2][ind], FU0[2][ind1], 1., dist);
4648 gammap = Interpola(FGamma[ind], FGamma[ind1], 1., dist);
4649 Rmezclap = Interpola(FRMezcla[ind], FRMezcla[ind1], 1., dist);
4650 diamep = Interpola(FDiametroTubo[ind], FDiametroTubo[ind1], 1., dist);
4651 if(FCoefAjusTC != 0 || FCoefAjusFric != 0) {
4652 tptubop = Interpola(FTPTubo[0][ind], FTPTubo[0][ind1], 1., dist);
4653 hip = Interpola(Fhi[ind], Fhi[ind1], 1., dist);
4654 rhop = Interpola(Frho[ind], Frho[ind1], 1., dist);
4655 Rep = Interpola(FRe[ind], FRe[ind1], 1., dist);
4659 double gamma1p = __Gamma::G1(gammap);
4660 double gamma3p = __Gamma::G3(gammap);
4661 double gamma5p = __Gamma::G5(gammap);
4662 velocidadp = w1 / w0;
4663 asonidop = sqrt(gammap * gamma1p * (w2 / w0 -
pow2(velocidadp) / 2));
4664 caracteristica = asonidop - signo * gamma3p * velocidadp;
4670 if(FCoefAjusTC != 0) {
4674 double tgasp =
pow2(asonidop) / (gammap * Rmezclap);
4676 TransmisionCalor(tgasp, diamep, q, hip, rhop, tptubop);
4678 double dacal = gamma3p * gamma1p * DeltaTiempo * q * FCoefAjusTC / asonidop;
4680 caracteristica += dacal;
4686 double presionp = __units::PaToBar((w2 -
pow2(w1) / 2. / w0) * gamma1p / __geom::Circle_area(diamep));
4687 double entropiap = asonidop / pow(presionp, gamma5p);
4688 double increentropia = entropia * __cons::ARef - entropiap;
4689 double daen = asonidop * increentropia / entropiap;
4690 caracteristica += daen;
4694 if(FDiametroTubo[ind1] != FDiametroTubo[ind]) {
4697 daar = -gamma3p * asonidop * velocidadp * 2 * (FDiametroTubo[ind1] - FDiametroTubo[ind]) * DeltaTiempo /
4699 }
else if(signo == -1) {
4700 daar = -gamma3p * asonidop * velocidadp * 2 * (FDiametroTubo[ind] - FDiametroTubo[ind1]) * DeltaTiempo /
4703 caracteristica += daar;
4709 if(velocidadp != 0. && FCoefAjusFric != 0.) {
4711 double velabs = fabs(velocidadp);
4714 Colebrook(FFriccion, diamep, f, Rep);
4716 double dafric = signo * gamma1p * (1. + signo * gamma1p * velocidadp / asonidop) * f * FCoefAjusFric *
pow3(
4717 velocidadp) * DeltaTiempo / (diamep * velabs);
4719 caracteristica += dafric;
4723 }
catch(exception & N) {
4724 std::cout <<
"ERROR: TTubo::Calculo_Caracteristica " << FNumeroTubo << std::endl;
4725 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4739 BC[FNodoIzq - 1]->PutLanda(FTuboCCNodoIzq, FAsonido0[0] + FGamma3[0] * FVelocidad0[0]);
4740 BC[FNodoIzq - 1]->PutBeta(FTuboCCNodoIzq, FAsonido0[0] - FGamma3[0] * FVelocidad0[0]);
4741 BC[FNodoIzq - 1]->PutEntropia(FTuboCCNodoIzq, FAsonido0[0] / pow(FPresion0[0], FGamma5[0]));
4743 BC[FNodoDer - 1]->PutLanda(FTuboCCNodoDer, FAsonido0[FNin - 1] + FGamma3[FNin - 1] * FVelocidad0[FNin - 1]);
4744 BC[FNodoDer - 1]->PutBeta(FTuboCCNodoDer, FAsonido0[FNin - 1] - FGamma3[FNin - 1] * FVelocidad0[FNin - 1]);
4745 BC[FNodoDer - 1]->PutEntropia(FTuboCCNodoDer, FAsonido0[FNin - 1] / pow(FPresion0[FNin - 1], FGamma5[FNin - 1]));
4747 }
catch(exception & N) {
4748 std::cout <<
"ERROR: TTubo::InicializaCaracteristicas tubo:" << FNumeroTubo << std::endl;
4749 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4758 void TTubo::CalculaB() {
4762 double v = 0., p = 0., f = 0., tgas = 0., g = 0., q = 0., diamemed = 0., Rm = 0., Rm1 = 0.;
4763 double Remed, gamma, gamma1, Rmed, himed, rhomed, twallmed, Vmed, H1, H2, Hmed, Amed, rhoAmed;
4765 for(
int i = 0; i < FNin - 1; i++) {
4767 FTVD.Bvector[0][i] = 0.;
4769 FTVD.Bvector[1][i] = 0.;
4771 FTVD.Bvector[2][i] = 0.;
4773 if(FArea[i] != FArea[i + 1] || FCoefAjusFric != 0 || FCoefAjusTC != 0) {
4775 Rm = sqrtRhoA[i + 1] / sqrtRhoA[i + 1];
4777 gamma = (Rm * FGamma[i + 1] + FGamma[i]) / Rm1;
4781 Vmed = (Rm * FVelocidadDim[i + 1] + FVelocidadDim[i]) / Rm1;
4783 if(FArea[i] != FArea[i + 1] || FCoefAjusTC != 0) {
4785 H1 = 0.5 * FVelocidadDim[i] * FVelocidadDim[i] + FAsonidoDim[i] * FAsonidoDim[i] / gamma1;
4786 H2 = 0.5 * FVelocidadDim[i + 1] * FVelocidadDim[i + 1] + FAsonidoDim[i + 1] * FAsonidoDim[i + 1] / gamma1;
4787 Hmed = (Rm * H2 + H1) / Rm1;
4788 Amed = sqrt(gamma1 * (Hmed - 0.5 * Vmed * Vmed));
4789 rhomed = sqrt(Frho[i] * Frho[i + 1]);
4792 if(FCoefAjusFric != 0 || FCoefAjusTC != 0) {
4793 rhoAmed = sqrtRhoA[i + 1] * sqrtRhoA[i + 1];
4796 if(FArea[i] != FArea[i + 1]) {
4797 FTVD.Bvector[1][i] += rhomed * Amed * Amed / gamma * (FArea[i] - FArea[i + 1]);
4800 if(FCoefAjusFric != 0) {
4801 Remed = 0.5 * (FRe[i] + FRe[i + 1]);
4803 Colebrook(FFriccion, FDiametroD12[i], f, Remed);
4805 g = f * Vmed * Vmed * 2 / FDiametroD12[i] * FCoefAjusFric;
4807 g = -f * Vmed * Vmed * 2 / FDiametroD12[i] * FCoefAjusFric;
4809 FTVD.Bvector[1][i] += FXref * g * rhoAmed;
4813 if(FCoefAjusTC != 0) {
4814 Rmed = 0.5 * (FRMezcla[i] + FRMezcla[i + 1]);
4815 himed = 0.5 * (Fhi[i] + Fhi[i + 1]);
4816 twallmed = 0.5 * (FTPTubo[0][i] + FTPTubo[0][i + 1]);
4818 tgas = Amed * Amed / gamma / Rmed;
4820 TransmisionCalor(tgas, FDiametroD12[i], q, himed, rhomed, twallmed);
4822 q = q * FCoefAjusTC;
4824 FTVD.Bvector[2][i] = -FXref * q * rhoAmed;
4832 }
catch(exception & N) {
4833 std::cout <<
"ERROR: TTubo::CalculaB tubo:" << FNumeroTubo << std::endl;
4834 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4843 void TTubo::CalculaBmen() {
4847 double v, p, f, tgas, g, q, diamemed, B, Rm, Rm1, gamma, gamma1, Vmed, Amed, H1, H2, Hmed, rhoAmed;
4848 double Remed = 0., Rmed = 0., himed = 0., rhomed = 0., twallmed = 0.;
4850 for(
int i = 1; i < FNin; i++) {
4852 FTVD.Bmen[0][i] = 0.;
4854 FTVD.Bmen[1][i] = 0.;
4856 FTVD.Bmen[2][i] = 0.;
4858 if(FArea[i] != FArea12[i - 1] || FCoefAjusFric != 0 || FCoefAjusTC != 0) {
4860 rhoAmed = sqrt(sqrtRhoA[i - 1] *
pow3(sqrtRhoA[i]));
4861 B = FU0[0][i] + rhoAmed + sqrtRhoA[i] * sqrtRhoA[i - 1];
4862 Rm = B / sqrt(
pow3(sqrtRhoA[i - 1]) * sqrtRhoA[i]);
4864 gamma = (Rm * FGamma[i] + FGamma[i - 1]) / Rm1;
4868 Vmed = (Rm * FVelocidadDim[i] + FVelocidadDim[i - 1]) / Rm1;
4870 if(FArea[i] != FArea12[i - 1] || FCoefAjusTC != 0) {
4872 H1 = 0.5 * FVelocidadDim[i] * FVelocidadDim[i] + FAsonidoDim[i] * FAsonidoDim[i] / gamma1;
4873 H2 = 0.5 * FVelocidadDim[i - 1] * FVelocidadDim[i - 1] + FAsonidoDim[i - 1] * FAsonidoDim[i - 1] / gamma1;
4874 Hmed = (Rm * H1 + H2) / Rm1;
4875 Amed = sqrt(gamma1 * (Hmed - 0.5 * Vmed * Vmed));
4876 rhomed = sqrt(Frho[i] * sqrt(Frho[i] * Frho[i - 1]));
4879 if(FCoefAjusFric != 0 || FCoefAjusTC != 0) {
4881 diamemed = 0.75 * FDiametroTubo[i] + 0.25 * FDiametroTubo[i - 1];
4884 if(FArea[i] != FArea12[i - 1]) {
4885 FTVD.Bmen[1][i] += rhomed * Amed * Amed / gamma * (FArea12[i - 1] - FArea[i]);
4888 if(FCoefAjusFric != 0) {
4889 Remed = (Rm * FRe[i] + FRe[i - 1]) / Rm1;
4891 Colebrook(FFriccion, diamemed, f, Remed);
4893 g = f * Vmed * Vmed * 2 / diamemed * FCoefAjusFric;
4895 g = -f * Vmed * Vmed * 2 / diamemed * FCoefAjusFric;
4897 FTVD.Bmen[1][i] += FXref * g * rhoAmed;
4901 if(FCoefAjusTC != 0) {
4902 Rmed = (Rm * FRMezcla[i] + FRMezcla[i - 1]) / Rm1;
4903 himed = (Rm * Fhi[i] + Fhi[i - 1]) / Rm1;
4904 twallmed = (Rm * FTPTubo[0][i] + FTPTubo[0][i - 1]) / Rm1;
4906 tgas = Amed * Amed / gamma / Rmed;
4908 TransmisionCalor(tgas, diamemed, q, himed, rhomed, twallmed);
4910 q = q * FCoefAjusTC;
4912 FTVD.Bmen[2][i] = -FXref * q * rhoAmed;
4920 }
catch(exception & N) {
4921 std::cout <<
"ERROR: TTubo::CalculaBmen tubo:" << FNumeroTubo << std::endl;
4922 std::cout <<
"Tipo de error: " << N.what() << std::endl;
4931 void TTubo::CalculaBmas() {
4935 double v, p, f, tgas, g, q, diamemed, B, Rm, Rm1, gamma, gamma1, Vmed, Amed, H1, H2, Hmed, rhoAmed;
4936 double Remed = 0., Rmed = 0., himed = 0., rhomed = 0., twallmed = 0.;
4938 for(
int i = 0; i < FNin - 1; i++) {
4940 FTVD.Bmas[0][i] = 0.;
4942 FTVD.Bmas[1][i] = 0.;
4944 FTVD.Bmas[2][i] = 0.;
4946 if(FArea[i] != FArea12[i] || FCoefAjusFric != 0 || FCoefAjusTC != 0) {
4948 rhoAmed = sqrt(sqrtRhoA[i + 1] *
pow3(sqrtRhoA[i]));
4949 B = FU0[0][i] + rhoAmed + sqrtRhoA[i] * sqrtRhoA[i + 1];
4950 Rm = B / sqrt(
pow3(sqrtRhoA[i + 1]) * sqrtRhoA[i]);
4952 gamma = (Rm * FGamma[i] + FGamma[i + 1]) / Rm1;
4956 Vmed = (Rm * FVelocidadDim[i] + FVelocidadDim[i + 1]) / Rm1;
4958 if(FArea[i] != FArea12[i] || FCoefAjusTC != 0) {
4960 H1 = 0.5 * FVelocidadDim[i] * FVelocidadDim[i] + FAsonidoDim[i] * FAsonidoDim[i] / gamma1;
4961 H2 = 0.5 * FVelocidadDim[i + 1] * FVelocidadDim[i + 1] + FAsonidoDim[i + 1] * FAsonidoDim[i + 1] / gamma1;
4962 Hmed = (Rm * H1 + H2) / Rm1;
4963 Amed = sqrt(gamma1 * (Hmed - 0.5 * Vmed * Vmed));
4964 rhomed = sqrt(Frho[i] * sqrt(Frho[i] * Frho[i + 1]));
4967 if(FCoefAjusFric != 0 || FCoefAjusTC != 0) {
4969 diamemed = 0.75 * FDiametroTubo[i] + 0.25 * FDiametroTubo[i + 1];
4972 if(FArea[i] != FArea12[i]) {
4973 FTVD.Bmas[1][i] += rhomed * Amed * Amed / gamma * (FArea[i] - FArea12[i]);
4976 if(FCoefAjusFric != 0) {
4977 Remed = (Rm * FRe[i] + FRe[i + 1]) / Rm1;
4979 Colebrook(FFriccion, diamemed, f, Remed);
4981 g = f * Vmed * Vmed * 2 / diamemed * FCoefAjusFric;
4983 g = -f * Vmed * Vmed * 2 / diamemed * FCoefAjusFric;
4985 FTVD.Bmas[1][i] += FXref * g * rhoAmed;
4989 if(FCoefAjusTC != 0) {
4990 Rmed = (Rm * FRMezcla[i] + FRMezcla[i + 1]) / Rm1;
4991 himed = (Rm * Fhi[i] + Fhi[i + 1]) / Rm1;
4992 twallmed = (Rm * FTPTubo[0][i] + FTPTubo[0][i + 1]) / Rm1;
4994 tgas = Amed * Amed / gamma / Rmed;
4996 TransmisionCalor(tgas, diamemed, q, himed, rhomed, twallmed);
4998 q = q * FCoefAjusTC;
5000 FTVD.Bmas[2][i] = -FXref * q * rhoAmed;
5008 }
catch(exception & N) {
5009 std::cout <<
"ERROR: TTubo::CalculaBmas tubo:" << FNumeroTubo << std::endl;
5010 std::cout <<
"Tipo de error: " << N.what() << std::endl;
5019 void TTubo::CalculaMatrizJacobiana() {
5023 double Rmed = 0., Vmed = 0., Vmed2 = 0., Hmed = 0., Amed = 0., Amed2 = 0., gamma = 0., H1 = 0., H2 = 0., Rmed1 = 0.;
5024 double *Ymed, gamma1, gamma2, gaAmed2;
5026 Ymed =
new double[FNumeroEspecies - 1 - FIntEGR];
5028 for(
int i = 0; i < FNin - 1; i++) {
5033 Rmed = sqrtRhoA[i + 1] / sqrtRhoA[i + 1];
5035 gamma = (Rmed * FGamma[i + 1] + FGamma[i]) / Rmed1;
5037 gamma2 = gamma1 / 2;
5038 H1 = 0.5 * FVelocidadDim[i] * FVelocidadDim[i] + FAsonidoDim[i] * FAsonidoDim[i] / gamma1;
5039 H2 = 0.5 * FVelocidadDim[i + 1] * FVelocidadDim[i + 1] + FAsonidoDim[i + 1] * FAsonidoDim[i + 1] / gamma1;
5040 Vmed = (Rmed * FVelocidadDim[i + 1] + FVelocidadDim[i]) / Rmed1;
5041 Vmed2 = Vmed * Vmed;
5042 Hmed = (Rmed * H2 + H1) / Rmed1;
5043 Amed2 = gamma1 * (Hmed - 0.5 * Vmed2);
5045 for(
int j = 0; j < FNumeroEspecies - 1 - FIntEGR; j++) {
5046 Ymed[j] = (Rmed * FFraccionMasicaEspecie[i + 1][j] + FFraccionMasicaEspecie[i][j]) / Rmed1;
5049 FTVD.Pmatrix[1][0][i] = Vmed - Amed;
5050 FTVD.Pmatrix[1][1][i] = Vmed;
5051 FTVD.Pmatrix[1][2][i] = Vmed + Amed;
5053 FTVD.Pmatrix[2][0][i] = Hmed - Vmed * Amed;
5054 FTVD.Pmatrix[2][1][i] = Vmed2 / 2.;
5055 FTVD.Pmatrix[2][2][i] = Hmed + Vmed * Amed;
5057 gaAmed2 = gamma2 / Amed2;
5059 FTVD.Qmatrix[0][0][i] = (Vmed / Amed + gaAmed2 * Vmed2) / 2.;
5060 FTVD.Qmatrix[0][1][i] = -0.5 / Amed - gaAmed2 * Vmed;
5061 FTVD.Qmatrix[0][2][i] = gaAmed2;
5063 FTVD.Qmatrix[1][0][i] = 1. - gaAmed2 * Vmed2;
5064 FTVD.Qmatrix[1][1][i] = gamma1 * Vmed / Amed2;
5065 FTVD.Qmatrix[1][2][i] = -gamma1 / Amed2;
5067 FTVD.Qmatrix[2][0][i] = -(Vmed / Amed + gaAmed2 * Vmed2) / 2.;
5068 FTVD.Qmatrix[2][1][i] = 0.5 / Amed - gaAmed2 * Vmed;
5069 FTVD.Qmatrix[2][2][i] = gaAmed2;
5072 FTVD.Alpha[0][i] = FTVD.Pmatrix[1][0][i];
5073 FTVD.Alpha[1][i] = Vmed;
5074 FTVD.Alpha[2][i] = FTVD.Pmatrix[1][2][i];
5076 for(
int j = 3; j < FNumEcuaciones; j++) {
5077 FTVD.Alpha[j][i] = (Vmed);
5083 }
catch(exception & N) {
5084 std::cout <<
"ERROR: TTubo::CalculaMatrizJacobiana tubo:" << FNumeroTubo << std::endl;
5085 std::cout <<
"Tipo de error: " << N.what() << std::endl;
5094 void TTubo::TVD_Estabilidad() {
5098 double VTotalMax = 0.;
5099 double VLocal = 0., DeltaT_tvd = 0.;
5101 CalculaFlujo(FU0, FTVD.W, FGamma, FGamma1, FNin);
5103 CalculaMatrizJacobiana();
5107 for(
int i = 0; i < FNin - 1; i++) {
5108 for(
int k = 0; k < 3; ++k) {
5109 FTVD.DeltaU[k][i] = FTVD.Qmatrix[k][0][i] * (FU0[0][i + 1] - FU0[0][i]) + FTVD.Qmatrix[k][1][i] *
5110 (FU0[1][i + 1] - FU0[1][i]) + FTVD.Qmatrix[k][2][i] * (FU0[2][i + 1] - FU0[2][i]);
5111 FTVD.DeltaB[k][i] = FTVD.Qmatrix[k][0][i] * FTVD.Bvector[0][i] + FTVD.Qmatrix[k][1][i] * FTVD.Bvector[1][i] +
5112 FTVD.Qmatrix[k][2][i] * FTVD.Bvector[2][i];
5113 FTVD.DeltaW[k][i] = FTVD.Qmatrix[k][0][i] * (FTVD.W[0][i + 1] - FTVD.W[0][i] + FTVD.Bvector[0][i]) +
5114 FTVD.Qmatrix[k][1][i] * (FTVD.W[1][i + 1] - FTVD.W[1][i] + FTVD.Bvector[1][i]) + FTVD.Qmatrix[k][2][i] *
5115 (FTVD.W[2][i + 1] - FTVD.W[2][i] + FTVD.Bvector[2][i]);
5116 if(fabs(FTVD.DeltaU[k][i]) < 1e-3) {
5117 FTVD.Beta[k][i] = FTVD.DeltaB[k][i] / 1e-3;
5119 FTVD.Beta[k][i] = FTVD.DeltaB[k][i] / FTVD.DeltaU[k][i];
5128 if(FTVD.Alpha[k][i] + FTVD.Beta[k][i] != 0) {
5129 if((VLocal = fabs(FTVD.Alpha[k][i]) + fabs(FTVD.Beta[k][i])) > VTotalMax) {
5134 for(
int k = 3; k < FNumEcuaciones; k++) {
5135 FTVD.Beta[k][i] = 0.;
5138 DeltaT_tvd = FCourant * FXref / VTotalMax;
5142 }
catch(exception & N) {
5143 std::cout <<
"ERROR: TTubo::TVD_Estabilidad tubo:" << FNumeroTubo << std::endl;
5144 std::cout <<
"Tipo de error: " << N.what() << std::endl;
5153 void TTubo::TVD_Limitador() {
5163 for(
int i = 0; i < FNin - 1; ++i) {
5165 for(
int k = 0; k < FNumEcuaciones; k++) {
5166 FTVD.LandaD[k][i] = dtdx * (FTVD.Alpha[k][i] + FTVD.Beta[k][i]);
5167 if(FTVD.LandaD[k][i] >= 0.) {
5168 FTVD.hLandaD[k][i] = 1;
5170 FTVD.hLandaD[k][i] = -1;
5174 for(
int i = 1; i < FNin - 2; ++i) {
5175 for(
int k = 0; k < 3; k++) {
5176 double den = (((double) FTVD.hLandaD[k][i] - FTVD.LandaD[k][i]) * FTVD.DeltaW[k][i]);
5177 double num = ((double) FTVD.hLandaD[k][i - FTVD.hLandaD[k][i]] - FTVD.LandaD[k][i - FTVD.hLandaD[k][i]]) *
5178 (FTVD.DeltaW[k][i - FTVD.hLandaD[k][i]]);
5179 if(fabs(den) > 1e-10)
5180 FTVD.R[k][i] = num / den;
5185 for(
int k = 3; k < FNumEcuaciones; k++) {
5186 double num = ((double) FTVD.hLandaD[k][i - FTVD.hLandaD[k][i]] - FTVD.LandaD[k][i - FTVD.hLandaD[k][i]]) *
5187 (FTVD.W[k][i + 1 - FTVD.hLandaD[k][i]] - FTVD.W[k][i - FTVD.hLandaD[k][i]]);
5188 double den = ((double) FTVD.hLandaD[k][i] - FTVD.LandaD[k][i]) * (FTVD.W[k][i + 1] - FTVD.W[k][i]);
5189 if(fabs(den) > 1e-10)
5190 FTVD.R[k][i] = num / den;
5195 for(
int k = 0; k < FNumEcuaciones; k++) {
5196 FTVD.R[k][0] = FTVD.R[k][1];
5197 FTVD.R[k][FNin - 2] = FTVD.R[k][FNin - 3];
5199 for(
int i = 0; i < FNin - 1; ++i) {
5200 for(
int k = 0; k < 3; k++) {
5201 FTVD.Phi[k][i] = (double) FTVD.hLandaD[k][i] - Limita(FTVD.R[k][i]) * ((double) FTVD.hLandaD[k][i] - FTVD.LandaD[k][i]);
5205 for(
int i = 0; i < FNin - 1; ++i) {
5206 for(
int k = 0; k < 3; ++k) {
5207 FTVD.gflux[k][i] = 0.5 * (FTVD.W[k][i] + FTVD.W[k][i + 1] - FTVD.Bmas[k][i] + FTVD.Bmen[k][i + 1] -
5208 (FTVD.Pmatrix[k][0][i] * FTVD.Phi[0][i] * FTVD.DeltaW[0][i] + FTVD.Pmatrix[k][1][i] * FTVD.Phi[1][i] * FTVD.DeltaW[1][i]
5209 + FTVD.Pmatrix[k][2][i] * FTVD.Phi[2][i] * FTVD.DeltaW[2][i]));
5211 for(
int k = 3; k < FNumEcuaciones; k++) {
5212 FTVD.gflux[k][i] = 0.5 * (FTVD.W[k][i] + FTVD.W[k][i + 1] - (double) FTVD.hLandaD[k][i] *
5213 (FTVD.W[k][i + 1] - FTVD.W[k][i])) + 0.5 * Limita(FTVD.R[k][i]) * ((double) FTVD.hLandaD[k][i] - FTVD.LandaD[k][i]) *
5214 (FTVD.W[k][i + 1] - FTVD.W[k][i]);
5218 for(
int i = 1; i < FNin - 1; ++i) {
5219 for(
int k = 0; k < FNumEcuaciones; k++) {
5220 FU1[k][i] = FU0[k][i] - dtdx * ((FTVD.gflux[k][i] - FTVD.gflux[k][i - 1]) + (FTVD.Bmen[k][i] + FTVD.Bmas[k][i]));
5224 }
catch(exception & N) {
5225 std::cout <<
"ERROR: TTubo::TVD_Limitador tubo:" << FNumeroTubo << std::endl;
5226 std::cout <<
"Tipo de error: " << N.what() << std::endl;
5235 inline double TTubo::Limita(
double r) {
5236 double ret_val = 0.;
5246 ret_val = (fabs(r) + r) / (1. + fabs(r));
5248 if(ret_val != ret_val)
5257 ret_val = (fabs(r) + r) / (1. + fabs(r));
5259 if(std::isnan(ret_val)) {
5318 void TTubo::DimensionaTVD() {
5323 FTVD.Bmas =
new double*[FNumEcuaciones];
5324 FTVD.Bvector =
new double*[FNumEcuaciones];
5325 FTVD.Bmen =
new double*[FNumEcuaciones];
5326 FTVD.gflux =
new double*[FNumEcuaciones];
5327 FTVD.Alpha =
new double*[FNumEcuaciones];
5328 FTVD.Beta =
new double*[FNumEcuaciones];
5329 FTVD.DeltaU =
new double*[FNumEcuaciones];
5330 FTVD.DeltaB =
new double*[FNumEcuaciones];
5331 FTVD.DeltaW =
new double*[FNumEcuaciones];
5332 FTVD.hLandaD =
new int*[FNumEcuaciones];
5333 FTVD.LandaD =
new double*[FNumEcuaciones];
5334 FTVD.Phi =
new double*[FNumEcuaciones];
5335 FTVD.R =
new double*[FNumEcuaciones];
5336 FTVD.W =
new double*[FNumEcuaciones];
5337 FTVD.Qmatrix =
new double**[FNumEcuaciones];
5338 FTVD.Pmatrix =
new double**[FNumEcuaciones];
5340 sqrtRhoA =
new double[FNin];
5342 for(
int k = 0; k < FNumEcuaciones; ++k) {
5343 FTVD.Bmas[k] =
new double[FNin];
5344 FTVD.Bvector[k] =
new double[FNin];
5345 FTVD.Bmen[k] =
new double[FNin];
5346 FTVD.gflux[k] =
new double[FNin];
5347 FTVD.Alpha[k] =
new double[FNin];
5348 FTVD.Beta[k] =
new double[FNin];
5349 FTVD.DeltaU[k] =
new double[FNin];
5350 FTVD.DeltaB[k] =
new double[FNin];
5351 FTVD.DeltaW[k] =
new double[FNin];
5352 FTVD.hLandaD[k] =
new int[FNin];
5353 FTVD.LandaD[k] =
new double[FNin];
5354 FTVD.Phi[k] =
new double[FNin];
5355 FTVD.R[k] =
new double[FNin];
5356 FTVD.W[k] =
new double[FNin];
5359 for(
int k = 0; k < FNumEcuaciones; ++k) {
5360 FTVD.Qmatrix[k] =
new double*[FNumEcuaciones];
5361 FTVD.Pmatrix[k] =
new double*[FNumEcuaciones];
5363 for(
int i = 0; i < FNumEcuaciones; i++) {
5364 FTVD.Qmatrix[k][i] =
new double[FNin];
5365 FTVD.Pmatrix[k][i] =
new double[FNin];
5368 for(
int i = 0; i < FNin; i++) {
5370 for(
int k = 3; k < FNumEcuaciones; k++) {
5371 for(
int j = 3; j < FNumEcuaciones; j++) {
5372 FTVD.Pmatrix[k][j][i] = 0.;
5375 FTVD.Pmatrix[0][0][i] = 1.;
5376 FTVD.Pmatrix[0][1][i] = 1.;
5377 FTVD.Pmatrix[0][2][i] = 1.;
5378 for(
int j = 3; j < FNumEcuaciones; j++) {
5379 FTVD.Pmatrix[j][j][i] = 1.;
5382 for(
int k = 3; k < FNumEcuaciones; k++) {
5383 for(
int j = 3; j < FNumEcuaciones; j++) {
5384 FTVD.Qmatrix[k][j][i] = 0.;
5387 for(
int j = 3; j < FNumEcuaciones; j++) {
5388 FTVD.Qmatrix[j][j][i] = 1.;
5392 }
catch(exception & N) {
5393 std::cout <<
"ERROR: TTubo::TVD_Limitador tubo:" << FNumeroTubo << std::endl;
5394 std::cout <<
"Tipo de error: " << N.what() << std::endl;
5400 void TTubo::RoeConstants() {
5402 sqrtRhoA[0] = sqrt(FU1[0][0]);
5403 for(
int i = 0; i < FNin - 1; i++) {
5404 sqrtRhoA[i + 1] = sqrt(FU1[0][i + 1]);
5413 return FAsonido0[i];
5417 return FCoefTurbulencia[i];
5421 return FCpMezcla[i];
5425 return FCvMezcla[i];
5437 return FDiametroTubo[i];
5441 return FComposicionInicial[i];
5449 return FLongitudTotal;
5469 return FPresion0[i];
5485 return FTIniParedTub;
5498 return FTPTubo[j][i];
5502 return FTParedAnt[j][i];
5506 return FVelocidad0[i];
5530 FTPTubo[k][i] = valor;
5537 #pragma package(smart_init) double Interpola_Entropia(nmPipeEnd TipoExtremoTubo, double DeltaTiempo)
void HeaderInstantaneousResults(std::stringstream &insoutput, stEspecies *DatosEspecies) const
void PutTime0(double valor)
Sets the current time.
double GetArea(int i) const
Gets the cross section at a given cell.
void ActualizaValoresNuevos(TCondicionContorno **BC)
T pow4(T x)
Returns x to the power of 4.
double getTime1() const
Gets the time at the following time-step.
int getNodoIzq() const
Gets the left-hand side node.
TTubo(int SpeciesNumber, int j, double SimulationDuration, TBloqueMotor **Engine, nmTipoCalculoEspecies SpeciesModel, nmCalculoGamma GammaCalculation, bool ThereIsEGR)
void LeeDatosGeometricosTubo(const char *FileWAM, fpos_t &filepos, double ene, int tipomallado, TBloqueMotor **Engine)
double GetDiametro(int i) const
Gets the cell diameter.
double FTime1
Time at following time step.
void ImprimeResultadosInstantaneos(std::stringstream &insoutput) const
double GetTPTuboAnt(int j, int) const
Gets the previous wall temperature at a given cell and node.
int getNodoDer() const
Gets the right-hand side node.
double Interpola_Caracteristica(double entropia, int signo, int extremo, double DeltaTiempo)
void ComunicacionDPF(TCondicionContorno **CC, TDeposito **Deposito)
double GetCoefTurbulencia(int i) const
Gets the turbulence coefficient.
void CalculaVariablesFundamentales()
T pow150(T x)
Returns x to the power of 1.5.
void CalculaCaracteristicasExtremos(TCondicionContorno **BC, double DeltaTiempo)
double GetTPTubo(int j, int i) const
Gets the wall temperature at a given cell and node.
double getTemperaturaInicial() const
Gets the initial temperature.
double zbrent(T &func, const double &x1, const double &x2, const double &tol)
Solves a function using Brent's method.
double GetGamma(int i) const
Gets the specific heat capacities ratio at a given cell.
std::string PutLabel(int idx)
Returns an integer.
void ImprimeResultadosMedios(std::stringstream &medoutput) const
double FDeltaTime
Time step.
double getTempWallIni() const
Gets the initial temperature of the wall.
void InicializaCaracteristicas(TCondicionContorno **BC)
double GetPresion(int i) const
Gets the fluid pressure.
double GetVelPro(int i) const
Gets the integrated gas velocity.
T pow075(T x)
Returns x to the power of 0.75.
double GetCpMezcla(int i) const
Gets the specific heat capacity at constant pressure at a given cell.
void PutDeltaTime(double valor)
Sets the time step.
T pow2(T x)
Returns x to the power of 2.
void HeaderAverageResults(std::stringstream &medoutput, stEspecies *DatosEspecies) const
void AjustaPaso(double Intervalo)
T pow3(T x)
Returns x to the power of 3.
double getDeltaTime() const
Gets the time step.
void CalculaResistenciasdePared(TCondicionContorno **BC)
double GetVelocidad(int i) const
Gets the fluid speed.
void PutTime1(double valor)
Sets the time after the following time-step.
void PutTPTubo(int k, int i, double valor)
Sets the wall temperature at a given cell and node.
void ReduccionFlujoSubsonico()
void CalculaTemperaturaPared(TBloqueMotor **Engine, double Theta, double CrankAngle, TCondicionContorno **BC)
double GetFraccionMasicaInicial(int i) const
Gets the initial mass fraction of species i.
void LeeDatosGeneralesTubo(const char *FileWAM, fpos_t &filepos)
double GetCvMezcla(int i) const
Gets the specific heat capacity at constant volume at a given cell.
void ComunicacionTubo_CC(TCondicionContorno **BC)
void SalidaGeneralTubos(stEspecies *DatosEspecies) const
void ReadInstantaneousResultsTubo(const char *FileWAM, fpos_t &filepos, bool HayMotor)
T pow025(T x)
Returns x to the power of 0.25.
double FTime0
Time at current time step.
double getPresionInicial() const
Gets the initial pressure.
double getMallado() const
Gets the mesh size.
double GetDensidad(int i) const
Gets the density.
void PutVelPro(int i, double valor)
Sets the integrated gas velocity.
void EstabilidadMetodoCalculo()
double GetRMezcla(int i) const
Gets the gas constant of the mixture at a given cell.
int getNumeroTubo() const
Gets the pipe id.
void CalculaCoeficientePeliculaInterior(TCondicionContorno **BC)
double getTime0() const
Gets the current time.
void ReduccionFlujoSubsonicoFCT()
void IniciaVariablesFundamentalesTubo()
double GetAsonido(int i) const
Gets the speed of sound.
double getVelocidadMedia() const
Gets the mean speed.
void CalculaResultadosInstantaneos()
void ActualizaPropiedadesGas()
void IniciaVariablesTransmisionCalor(TCondicionContorno **BC, TBloqueMotor **Engine, double AmbientTemperature)
void CalculaResultadosMedios(double Theta)
void CalculaCoeficientePeliculaExterior(TBloqueMotor **Engine, double AmbientPressure, double AmbientTemperature)
void ReadAverageResultsTubo(const char *FileWAM, fpos_t &filepos, bool HayMotor)
double getLongitudTotal() const
Gets the total length of the pipe.