OpenWAM
TDepVolCte.cpp
1 /* --------------------------------------------------------------------------------*\
2 |==========================|
3  |\\ /\ /\ // O pen | OpenWAM: The Open Source 1D Gas-Dynamic Code
4  | \\ | X | // W ave |
5  | \\ \/_\/ // A ction | CMT-Motores Termicos / Universidad Politecnica Valencia
6  | \\/ \// M odel |
7  ----------------------------------------------------------------------------------
8  License
9 
10  This file is part of OpenWAM.
11 
12  OpenWAM is free software: you can redistribute it and/or modify
13  it under the terms of the GNU General Public License as published by
14  the Free Software Foundation, either version 3 of the License, or
15  (at your option) any later version.
16 
17  OpenWAM is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  GNU General Public License for more details.
21 
22  You should have received a copy of the GNU General Public License
23  along with OpenWAM. If not, see <http://www.gnu.org/licenses/>.
24 
25 
26  \*-------------------------------------------------------------------------------- */
27 
28 // ---------------------------------------------------------------------------
29 #pragma hdrstop
30 
31 #include "TDepVolCte.h"
32 
33 #include "TCCDeposito.h"
34 #include "TCCUnionEntreDepositos.h"
35 #include "TTubo.h"
36 #include "TCompresor.h"
37 #include "TDPF.h"
38 #include "TCanalDPF.h"
39 
40 // ---------------------------------------------------------------------------
41 // ---------------------------------------------------------------------------
42 
43 TDepVolCte::TDepVolCte(int i, nmTipoCalculoEspecies SpeciesModel, int numeroespecies, nmCalculoGamma GammaCalculation,
44  bool ThereIsEGR) :
45  TDepVolCteBase(i, nmDepVolCte, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
46 }
47 
48 // ---------------------------------------------------------------------------
49 // ---------------------------------------------------------------------------
50 
51 TDepVolCte::~TDepVolCte() {
52 
53 }
54 
55 // ---------------------------------------------------------------------------
56 // ---------------------------------------------------------------------------
57 
58 void TDepVolCte::ActualizaPropiedades(double TimeCalculo) {
59 
60  double H = 0.; // Entalpia de entrada
61  double Energia = 0.;
62  double MasaEntrante, FraccionMasicaAcum = 0.;
63  double DeltaT = 0.;
64  double g = 0., v = 0., a = 0., m = 0., g1 = 0.;
65  int SignoFlujo = 1; // Inicializado por si el flujo esta parado.
66  int SignoFlujoED = 1; // Inicializado por si el flujo esta parado.
67 
68  try {
69  FMasa0 = FMasa;
70  MasaEntrante = 0.;
71  H = 0.;
72  DeltaT = TimeCalculo - FTime;
73 
74  if(FCalculoEspecies == nmCalculoCompleto) {
75 
76  FRMezcla = CalculoCompletoRMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2], 0,
77  FCalculoGamma, nmMEP);
78  FCpMezcla = CalculoCompletoCpMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2], 0,
79  __units::degCToK(FTemperature), FCalculoGamma, nmMEP);
80  FGamma = CalculoCompletoGamma(FRMezcla, FCpMezcla, FCalculoGamma);
81 
82  } else if(FCalculoEspecies == nmCalculoSimple) {
83 
84  FRMezcla = CalculoSimpleRMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FCalculoGamma, nmMEP);
85  FCvMezcla = CalculoSimpleCvMezcla(__units::degCToK(FTemperature), FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1],
86  FCalculoGamma, nmMEP);
87  FGamma = CalculoSimpleGamma(FRMezcla, FCvMezcla, FCalculoGamma);
88 
89  }
90 
91  bool Converge = false;
92  bool FirstStep = true;
93  double H0 = 0.;
94  double Asonido0 = FAsonido;
95  double Asonido1 = FAsonido;
96  double Error = 0.;
97  double Diff = 0.;
98  double Heat = 0.;
99  double MTemp = FGamma / (pow2(FAsonido * __cons::ARef) * FMasa0);
100 
101  while(!Converge) {
102  H = 0.;
103  for(int i = 0; i < FNumeroUniones; i++) {
104  if(FCCDeposito[i]->getTipoCC() == nmPipeToPlenumConnection) {
105  if(dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getSentidoFlujo() == nmEntrante) {
106  SignoFlujo = 1;
107  } else if(dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getSentidoFlujo() == nmSaliente) {
108  SignoFlujo = -1;
109  }
110  g = (double) - dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getMassflow();
111  if(!FCCDeposito[i]->getUnionDPF()) {
112  m = g * DeltaT * FCCDeposito[i]->GetTuboExtremo(0).Pipe->getNumeroConductos();
113  } else {
114 #ifdef ParticulateFilter
115  int NumeroCanales = 0;
116  int NumeroHaz = FCCDeposito[i]->GetTuboExtremo(0).NumeroHaz;
117  int TipoCanal = FCCDeposito[i]->GetTuboExtremo(0).TipoCanal;
118  NumeroCanales = FCCDeposito[i]->GetTuboExtremo(0).DPF->GetCanal(NumeroHaz, TipoCanal)->getNumeroCanales();
119  m = g * DeltaT * NumeroCanales;
120 #endif
121  }
122  v = (double) SignoFlujo * dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getVelocity();
123  a = dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getSpeedSound();
124  if(FirstStep) {
125  MasaEntrante += m;
126  for(int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
127  FMasaEspecie[j] += FCCDeposito[i]->GetFraccionMasicaEspecie(j) * m;
128  }
129  }
130  if(v > 0) {
131  H += EntalpiaEntrada(a, v, m, Asonido1, FMasa, FCCDeposito[i]->getGamma());
132  }
133  }
134  }
135  for(int i = 0; i < FNumeroUnionesED; i++) {
136 
137  if(FCCUnionEntreDep[i]->getTipoCC() == nmUnionEntreDepositos) {
138  if(FNumeroDeposito == dynamic_cast<TCCUnionEntreDepositos*>(FCCUnionEntreDep[i])->getNumeroDeposito1()) {
139  SignoFlujoED = dynamic_cast<TCCUnionEntreDepositos*>(FCCUnionEntreDep[i])->getSentidoFlujoED1();
140  } else if(FNumeroDeposito == dynamic_cast<TCCUnionEntreDepositos*>(FCCUnionEntreDep[i])->getNumeroDeposito2()) {
141  SignoFlujoED = dynamic_cast<TCCUnionEntreDepositos*>(FCCUnionEntreDep[i])->getSentidoFlujoED2();
142  } else {
143  printf("ERROR:TDepVolCte::ActualizaPropiedades en el deposito %d, union entre depositos %d\n", FNumeroDeposito, i);
144  }
145  g = (double) SignoFlujoED * dynamic_cast<TCCUnionEntreDepositos*>(FCCUnionEntreDep[i])->getMassflow();
146  m = g * DeltaT;
147  a = (double) dynamic_cast<TCCUnionEntreDepositos*>(FCCUnionEntreDep[i])->getSpeedSound();
148  v = (double) SignoFlujoED * dynamic_cast<TCCUnionEntreDepositos*>(FCCUnionEntreDep[i])->getVelocity() / __cons::ARef;
149  if(FirstStep) {
150  MasaEntrante += m;
151  for(int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
152  FMasaEspecie[j] += FCCUnionEntreDep[i]->GetFraccionMasicaEspecie(j) * m;
153  }
154  }
155  if(g > 0) {
156  H += EntalpiaEntrada(a, 0., m, Asonido1, FMasa, FCCUnionEntreDep[i]->getGamma());
157 
158  }
159 
160  }
161  }
162 
163  if(FHayCompresor) {
164  g = (double) FCompresorSentido * FCompresor->getMassflow();
165  m = g * DeltaT;
166  a = FCompresor->getSpeedSound();
167  if(FirstStep) {
168  MasaEntrante += m;
169  for(int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
170  FMasaEspecie[j] += FCompresor->GetFraccionMasicaEspecie(j) * m;
171  }
172  }
173  if(g > 0) {
174  H += EntalpiaEntrada(a, 0, m, Asonido1, FMasa, FCompresor->getGamma());
175  }
176  }
177  if(FirstStep) {
178  FMasa = FMasa0 + MasaEntrante;
179  for(int j = 0; j < FNumeroEspecies - 2; j++) {
180  FFraccionMasicaEspecie[j] = FMasaEspecie[j] / FMasa;
181  FraccionMasicaAcum += FFraccionMasicaEspecie[j];
182  }
183  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
184  if(FHayEGR)
185  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FMasaEspecie[FNumeroEspecies - 1] / FMasa;
186  FirstStep = false;
187  H0 = H;
188  }
189  Heat = FHeatPower * DeltaT * (MTemp + FGamma / (pow2(Asonido1 * __cons::ARef) * FMasa)) / 2.;
190  Energia = pow(FMasa / FMasa0 * exp((H0 + H) / 2 - Heat), __Gamma::G1(FGamma));
191 
192  Asonido1 = FAsonido * sqrt(Energia);
193  Error = (Diff = Asonido1 - Asonido0, fabs(Diff)) / Asonido1;
194  if(Error > 1e-6) {
195  Asonido0 = Asonido1;
196  } else {
197  Converge = true;
198  FAsonido = Asonido1;
199  }
200  }
201  double A2 = pow2(__cons::ARef * FAsonido);
202  FPressure = __units::PaToBar(A2 / FGamma / FVolumen * FMasa);
203  FPresionIsen = pow(FPressure / FPresRef, __Gamma::G5(FGamma));
204  FTemperature = __units::KTodegC(A2 / FGamma / FRMezcla);
205  FTime = TimeCalculo;
206  } catch(exception & N) {
207  std::cout << "ERROR: TDepVolCte::ActualizaPropiedades en el deposito: " << FNumeroDeposito << std::endl;
208  std::cout << "Tipo de error: " << N.what() << std::endl;
209  throw Exception(N.what());
210  }
211 }
212 
213 void TDepVolCte::UpdateProperties0DModel(double TimeCalculo) {
214 
215  ActualizaPropiedades(TimeCalculo);
216 
217  AcumulaResultadosMedios(TimeCalculo);
218 
219 }
220 
221 // ---------------------------------------------------------------------------
222 // ---------------------------------------------------------------------------
223 
224 #pragma package(smart_init)
TCCUnionEntreDepositos
Definition: TCCUnionEntreDepositos.h:41
TDepVolCteBase
Definition: TDepVolCteBase.h:33
Exception
Custom exception class.
Definition: Exception.hpp:39
TTubo.h
TCCDeposito
Definition: TCCDeposito.h:40
pow2
T pow2(T x)
Returns x to the power of 2.
Definition: Math_wam.h:88