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) {
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) {
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) {
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) {
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) {
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) {
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) {
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) {
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) {
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) {
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)