OpenWAM
TUnionDireccional.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 "TUnionDireccional.h"
32 
33 #include "TCCDeposito.h"
34 #include "TCondicionContorno.h"
35 #include "TTubo.h"
36 #include "TTipoValvula.h"
37 
38 // ---------------------------------------------------------------------------
39 // ---------------------------------------------------------------------------
40 
41 TUnionDireccional::TUnionDireccional(int i, int NUnionDireccional, nmTipoCalculoEspecies SpeciesModel,
42  int numeroespecies, nmCalculoGamma GammaCalculation, bool ThereIsEGR) :
43  TDepVolCteBase(i, nmUnionDireccional, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
44 
45  FNumUnionDireccional = NUnionDireccional;
46  FCCEntrada = NULL;
47 
48  FNodoEntrada = NULL;
49  FSentidoEntrada = NULL;
50  FCDSalidaInicial = NULL;
51  FVelocidadCorte = NULL;
52  FVelocidadFin = NULL;
53  FVelocity = NULL;
54  FCoefA = NULL;
55  FCoefB = NULL;
56 
57 }
58 
59 // ---------------------------------------------------------------------------
60 // ---------------------------------------------------------------------------
61 
62 TUnionDireccional::~TUnionDireccional() {
63 
64  if(FNodoEntrada != NULL)
65  delete[] FNodoEntrada;
66  if(FSentidoEntrada != NULL)
67  delete[] FSentidoEntrada;
68  if(FCDSalidaInicial != NULL)
69  delete[] FCDSalidaInicial;
70  if(FVelocidadCorte != NULL)
71  delete[] FVelocidadCorte;
72  if(FVelocidadFin != NULL)
73  delete[] FVelocidadFin;
74  if(FVelocity != NULL)
75  delete[] FVelocity;
76  if(FCoefA != NULL)
77  delete[] FCoefA;
78  if(FCoefB != NULL)
79  delete[] FCoefB;
80 
81  if(FCCEntrada != NULL)
82  delete[] FCCEntrada;
83 
84 }
85 
86 // ---------------------------------------------------------------------------
87 // ---------------------------------------------------------------------------
88 
89 void TUnionDireccional::LeeDatosUnionDireccional(const char *FileWAM, fpos_t &filepos) {
90  try {
91  int numid = 0; // Dato para Wamer
92 
93  FNodoEntrada = new int[2];
94  FSentidoEntrada = new int[2];
95  FCDSalidaInicial = new double[2];
96  FVelocidadCorte = new double[2];
97  FVelocidadFin = new double[2];
98  FVelocity = new double[2];
99  FCoefA = new double[2];
100  FCoefB = new double[2];
101 
102  FCCEntrada = new TCondicionContorno*[2];
103 
104  FILE *fich = fopen(FileWAM, "r");
105  fsetpos(fich, &filepos);
106 
107  fscanf(fich, "%d ", &numid); /* DATO PARA WAMER */
108  fscanf(fich, "%d %d %d ", &FNodoEntrada[0], &FNodoEntrada[1], &FNodoSalida);
109  /* Lectura de informacion para el calculo del coeficiente de descarga de salida
110  para las uniones de entrada al deposito de union direccional */
111  fscanf(fich, "%lf %lf %lf ", &FCDSalidaInicial[0], &FVelocidadCorte[0], &FVelocidadFin[0]);
112  fscanf(fich, "%lf %lf %lf ", &FCDSalidaInicial[1], &FVelocidadCorte[1], &FVelocidadFin[1]);
113 
114  fgetpos(fich, &filepos);
115  fclose(fich);
116 
117  FCoefA[0] = -FCDSalidaInicial[0] * FVelocidadFin[0] / (FVelocidadCorte[0] - FVelocidadFin[0]);
118  FCoefB[0] = FCDSalidaInicial[0] / (FVelocidadCorte[0] - FVelocidadFin[0]);
119  FCoefA[1] = -FCDSalidaInicial[1] * FVelocidadFin[1] / (FVelocidadCorte[1] - FVelocidadFin[1]);
120  FCoefB[1] = FCDSalidaInicial[1] / (FVelocidadCorte[1] - FVelocidadFin[1]);
121 
122  } catch(exception & N) {
123  std::cout << "ERROR: TUnionDireccional::LeeDatosUnionDireccional en el deposito: " << FNumeroDeposito << std::endl;
124  std::cout << "Tipo de error: " << N.what() << std::endl;
125  throw Exception(N.what());
126  }
127 }
128 
129 // ---------------------------------------------------------------------------
130 // ---------------------------------------------------------------------------
131 
132 void TUnionDireccional::AsignaCCUnionDireccional() {
133  try {
134 
135  if(FTipoDeposito == nmUnionDireccional) {
136  for(int i = 0; i < FNumeroUniones; i++) {
137  if(FCCDeposito[i]->getNumeroCC() == FNodoEntrada[0]) {
138  FCCEntrada[0] = FCCDeposito[i];
139  } else if(FCCDeposito[i]->getNumeroCC() == FNodoEntrada[1]) {
140  FCCEntrada[1] = FCCDeposito[i];
141  } else if(FCCDeposito[i]->getNumeroCC() == FNodoSalida) {
142  FCCSalida = FCCDeposito[i];
143  }
144  }
145  }
146 
147  } catch(exception & N) {
148  std::cout << "ERROR: TUnionDireccional::AsignaCCUnionDireccional en la union direccional " << FNumUnionDireccional <<
149  std::endl;
150  std::cout << "Tipo de error: " << N.what() << std::endl;
151  throw Exception(N.what());
152  }
153 }
154 
155 // ---------------------------------------------------------------------------
156 // ---------------------------------------------------------------------------
157 
158 void TUnionDireccional::ActualizaPropiedades(double TimeCalculo) {
159  try {
160  double H = 0.; // Entalpia de entrada
161  double Energia = 0.;
162  double MasaEntrante, FraccionMasicaAcum = 0.;
163  double DeltaT = 0.;
164  double g = 0., v = 0., a = 0., m = 0.;
165  int SignoFlujo = 1;
166 
167  FMasa0 = FMasa;
168  MasaEntrante = 0.;
169  H = 0.;
170  DeltaT = TimeCalculo - FTime;
171 
172  if(FCalculoEspecies == nmCalculoCompleto) {
173 
174  FRMezcla = CalculoCompletoRMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2], 0,
175  FCalculoGamma, nmMEP);
176  FCpMezcla = CalculoCompletoCpMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2], 0,
177  __units::degCToK(FTemperature), FCalculoGamma, nmMEP);
178  FGamma = CalculoCompletoGamma(FRMezcla, FCpMezcla, FCalculoGamma);
179 
180  } else if(FCalculoEspecies == nmCalculoSimple) {
181 
182  FRMezcla = CalculoSimpleRMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FCalculoGamma, nmMEP);
183  FCvMezcla = CalculoSimpleCvMezcla(__units::degCToK(FTemperature), FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1],
184  FCalculoGamma, nmMEP);
185  FGamma = CalculoSimpleGamma(FRMezcla, FCvMezcla, FCalculoGamma);
186 
187  }
188 
189  bool Converge = false;
190  bool FirstStep = true;
191  double H0 = 0.;
192  double Asonido0 = FAsonido;
193  double Asonido1 = FAsonido;
194  double Error = 0.;
195  double Diff = 0.;
196 
197  while(!Converge) {
198  H = 0.;
199  for(int i = 0; i < FNumeroUniones; i++) {
200  if(dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getSentidoFlujo() == nmEntrante) {
201  SignoFlujo = 1;
202  } else if(dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getSentidoFlujo() == nmSaliente) {
203  SignoFlujo = -1;
204  }
205  g = (double) - dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getMassflow();
206  v = (double) SignoFlujo * dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getVelocity();
207  a = dynamic_cast<TCCDeposito*>(FCCDeposito[i])->getSpeedSound();
208  m = g * DeltaT * FCCDeposito[i]->GetTuboExtremo(0).Pipe->getNumeroConductos();
209  if(FirstStep) {
210  MasaEntrante += m;
211  for(int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
212  FMasaEspecie[j] += FCCDeposito[i]->GetFraccionMasicaEspecie(j) * m;
213  }
214  }
215  if(v > 0) {
216  H += EntalpiaEntrada(a, v, m, FAsonido, FMasa, FCCDeposito[i]->getGamma());
217  }
218 
219  }
220  if(FirstStep) {
221  FMasa = FMasa0 + MasaEntrante;
222  for(int j = 0; j < FNumeroEspecies - 2; j++) {
223  FFraccionMasicaEspecie[j] = FMasaEspecie[j] / FMasa;
224  FraccionMasicaAcum += FFraccionMasicaEspecie[j];
225  }
226  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
227  if(FHayEGR)
228  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FMasaEspecie[FNumeroEspecies - 1] / FMasa;
229  FirstStep = false;
230  H0 = H;
231  }
232 
233  Energia = pow(FMasa / FMasa0 * exp((H + H0) / 2), __Gamma::G1(FGamma));
234  Asonido1 = FAsonido * sqrt(Energia);
235  Error = (Diff = Asonido1 - Asonido0, fabs(Diff)) / Asonido1;
236  if(Error > 1e-6) {
237  Asonido0 = Asonido1;
238  } else {
239  FAsonido = Asonido1;
240  Converge = true;
241  }
242  }
243  FTemperature = __units::KTodegC(pow2(FAsonido * __cons::ARef) / (FGamma * FRMezcla));
244  FPressure = __units::PaToBar(pow2(__cons::ARef * FAsonido) / FGamma / FVolumen * FMasa);
245  FPresionIsen = pow(FPressure / FPresRef, __Gamma::G5(FGamma));
246  FTime = TimeCalculo;
247  } catch(exception & N) {
248  std::cout << "ERROR: TUnionDireccional::ActualizaPropiedades en la union direccional: " << FNumUnionDireccional <<
249  std::endl;
250  std::cout << "Tipo de error: " << N.what() << std::endl;
251  throw Exception(N.what());
252  }
253 }
254 
255 // ---------------------------------------------------------------------------
256 // ---------------------------------------------------------------------------
257 
258 void TUnionDireccional::CalculoUnionDireccional() {
259  try {
260 
261  double coefA = 0., coefB = 0.;
262  /* Parametro independiente y pendiente de la recta para el calculo del Coeficiende de Descarga */
263 
264  for(int i = 0; i < 2; i++) {
265  if(dynamic_cast<TCCDeposito*>(FCCEntrada[i])->getSentidoFlujo() == nmEntrante) {
266  FSentidoEntrada[i] = 1;
267  } else if(dynamic_cast<TCCDeposito*>(FCCEntrada[i])->getSentidoFlujo() == nmSaliente) {
268  FSentidoEntrada[i] = -1;
269  } else
270  FSentidoEntrada[i] = 0; // Flujo parado
271  FVelocity[i] = FSentidoEntrada[i] * dynamic_cast<TCCDeposito*>(FCCEntrada[i])->getVelocity() * __cons::ARef;
272  }
273 
274  /* Calculo del coeficiente de descarga de salida en el Pipe de Entrada 0 */
275  if(FVelocity[1] <= FVelocidadCorte[0]) {
276  dynamic_cast<TCCDeposito*>(FCCEntrada[0])->PutCDSalida(FCDSalidaInicial[0]);
277  } else if(FVelocity[1] >= FVelocidadFin[0]) {
278  dynamic_cast<TCCDeposito*>(FCCEntrada[0])->PutCDSalida(0);
279  } else {
280  dynamic_cast<TCCDeposito*>(FCCEntrada[0])->PutCDSalida(FCoefA[0] + FCoefB[0] * FVelocity[1]);
281 
282  }
283 
284  /* Calculo del coeficiente de descarga de salida en el Pipe de Entrada 1 */
285  if(FVelocity[0] <= FVelocidadCorte[1]) {
286  dynamic_cast<TCCDeposito*>(FCCEntrada[1])->PutCDSalida(FCDSalidaInicial[1]);
287  } else if(FVelocity[0] >= FVelocidadFin[0]) {
288  dynamic_cast<TCCDeposito*>(FCCEntrada[1])->PutCDSalida(0);
289  } else {
290  dynamic_cast<TCCDeposito*>(FCCEntrada[1])->PutCDSalida(FCoefA[1] + FCoefB[1] * FVelocity[0]);
291 
292  }
293 
294  } catch(exception & N) {
295  std::cout << "ERROR: TUnionDireccional::CalculoUnionDireccional en la union direccional: " << FNumUnionDireccional <<
296  std::endl;
297  std::cout << "Tipo de error: " << N.what() << std::endl;
298  throw Exception(N.what());
299  }
300 }
301 
302 void TUnionDireccional::UpdateProperties0DModel(double TimeCalculo) {
303  ActualizaPropiedades(TimeCalculo);
304 
305  CalculoUnionDireccional();
306 
307  AcumulaResultadosMedios(TimeCalculo);
308 
309 }
310 
311 // ---------------------------------------------------------------------------
312 // ---------------------------------------------------------------------------
313 
314 #pragma package(smart_init)
TDepVolCteBase
Definition: TDepVolCteBase.h:33
TCondicionContorno
Definition: TCondicionContorno.h:54
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