OpenWAM
TCCDescargaExtremoAbierto.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 "TCCDescargaExtremoAbierto.h"
32 #include "TTubo.h"
33 
34 // ---------------------------------------------------------------------------
35 // ---------------------------------------------------------------------------
36 
37 TCCDescargaExtremoAbierto::TCCDescargaExtremoAbierto(nmTypeBC TipoCC, int numCC, nmTipoCalculoEspecies SpeciesModel,
38  int numeroespecies, nmCalculoGamma GammaCalculation, bool ThereIsEGR) :
39  TCondicionContorno(TipoCC, numCC, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
40  if(TipoCC == nmOpenEndAtmosphere) {
41  FTipoDescarga = nmDescargaAtmosfera;
42  } else if(TipoCC == nmOpenEndReservoir) {
43  FTipoDescarga = nmDescargaRemanso;
44  } else if(TipoCC == nmOpenEndCalcExtern) {
45  FTipoDescarga = nmDescargaRemansoMatlab;
46  } else
47  printf("ERROR:TCCDescargaExtremoAbierto:Asignacion Tipo BC,en la condicion de contorno: %d\n", FNumeroCC);
48 
49  FComposicion = NULL;
50  FTuboExtremo = NULL;
51  FVelocidadSonidoDep = 0;
52 
53 }
54 
55 // ---------------------------------------------------------------------------
56 // ---------------------------------------------------------------------------
57 
58 TCCDescargaExtremoAbierto::~TCCDescargaExtremoAbierto() {
59  delete[] FTuboExtremo;
60  if(FComposicion != NULL)
61  delete[] FComposicion;
62 }
63 
64 // ---------------------------------------------------------------------------
65 // ---------------------------------------------------------------------------
66 
67 void TCCDescargaExtremoAbierto::AsignAmbientConditions(double Tamb, double Pamb, double *AtmosphericComposition) {
68 
69  double RMezclaDep = 0., CvMezclaDep = 0., CpMezclaDep = 0., GammaDep = 0.;
70  int modeladoescape = 0;
71 
72  FPressure = Pamb;
73  FTemperaturaDep = Tamb;
74  // Inicializacion del transporte de especies quimicas.
75  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
76  FComposicion = new double[FNumeroEspecies - FIntEGR];
77  for(int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
78  FComposicion[i] = AtmosphericComposition[i];
79  FFraccionMasicaEspecie[i] = FTuboExtremo[0].Pipe->GetFraccionMasicaInicial(i);
80  }
81  if(FCalculoEspecies == nmCalculoCompleto) {
82 
83  RMezclaDep = CalculoCompletoRMezcla(FComposicion[0], FComposicion[1], FComposicion[2], 0, FCalculoGamma, nmMEP);
84  CpMezclaDep = CalculoCompletoCpMezcla(FComposicion[0], FComposicion[1], FComposicion[2], 0,
85  __units::degCToK(FTemperaturaDep), FCalculoGamma, nmMEP);
86  GammaDep = CalculoCompletoGamma(RMezclaDep, CpMezclaDep, FCalculoGamma);
87 
88  } else if(FCalculoEspecies == nmCalculoSimple) {
89 
90  RMezclaDep = CalculoSimpleRMezcla(FComposicion[0], 0, FCalculoGamma, nmMEP);
91  CvMezclaDep = CalculoSimpleCvMezcla(__units::degCToK(FTemperaturaDep), FComposicion[0], 0, FCalculoGamma, nmMEP);
92  GammaDep = CalculoSimpleGamma(RMezclaDep, CvMezclaDep, FCalculoGamma);
93 
94  }
95  FVelocidadSonidoDep = sqrt(__units::degCToK(FTemperaturaDep) * GammaDep * RMezclaDep) / __cons::ARef;
96 
97 }
98 
99 // ---------------------------------------------------------------------------
100 // ---------------------------------------------------------------------------
101 
102 void TCCDescargaExtremoAbierto::ReadBoundaryData(const char *FileWAM, fpos_t &filepos, int NumberOfPipes, TTubo **Pipe,
103  int nDPF, TDPF **DPF) {
104  try {
105  int i = 0;
106  double fracciontotal = 0.;
107  double RMezclaDep = 0., CvMezclaDep = 0., CpMezclaDep = 0., GammaDep = 0.;
108  int modeladoescape = 0;
109 
110  FTuboExtremo = new stTuboExtremo[1];
111  FTuboExtremo[0].Pipe = NULL;
112 
113  FPref = 1;
114 
115  while(FNumeroTubosCC < 1 && i < NumberOfPipes) {
116  if(Pipe[i]->getNodoIzq() == FNumeroCC) {
117  FTuboExtremo[FNumeroTubosCC].Pipe = Pipe[i];
118  FTuboExtremo[FNumeroTubosCC].TipoExtremo = nmLeft;
119  FCC = &(FTuboExtremo[FNumeroTubosCC].Beta);
120  FCD = &(FTuboExtremo[FNumeroTubosCC].Landa);
121  FNodoFin = 0;
122  FIndiceCC = 0;
123  FNumeroTubosCC++;
124  }
125  if(Pipe[i]->getNodoDer() == FNumeroCC) {
126  FTuboExtremo[FNumeroTubosCC].Pipe = Pipe[i];
127  FTuboExtremo[FNumeroTubosCC].TipoExtremo = nmRight;
128  FCC = &(FTuboExtremo[FNumeroTubosCC].Landa);
129  FCD = &(FTuboExtremo[FNumeroTubosCC].Beta);
130  FNodoFin = FTuboExtremo[FNumeroTubosCC].Pipe->getNin() - 1;
131  FIndiceCC = 1;
132  FNumeroTubosCC++;
133  }
134  i++;
135  }
136 
137  FILE *fich = fopen(FileWAM, "r");
138  fsetpos(fich, &filepos);
139 
140  if(FTipoDescarga == nmDescargaAtmosfera) { // DESCARGA A LA ATMOSFERA
141  fscanf(fich, "%lf ", &FPerdidaExtremo);
142 
143  } else if(FTipoDescarga == nmDescargaRemanso) { // DESCARGA A UN DEPOSITO DE REMANSO
144  fscanf(fich, "%lf %lf %lf ", &FPressure, &FTemperaturaDep, &FPerdidaExtremo);
145 
146  // Se determina si este remanso modela el escape del motor.
147  fscanf(fich, "%d ", &modeladoescape);
148  modeladoescape == 0 ? FModeladoEscape = false : FModeladoEscape = true;
149 
150  // Inicializacion del transporte de especies quimicas.
151  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
152  FComposicion = new double[FNumeroEspecies - FIntEGR];
153  for(int i = 0; i < FNumeroEspecies - 1; i++) {
154  fscanf(fich, "%lf ", &FComposicion[i]);
155  FFraccionMasicaEspecie[i] = FTuboExtremo[0].Pipe->GetFraccionMasicaInicial(i);
156  fracciontotal += FComposicion[i];
157  }
158  if(FHayEGR) {
159  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboExtremo[0].Pipe->GetFraccionMasicaInicial(FNumeroEspecies - 1);
160  if(FCalculoEspecies == nmCalculoCompleto) {
161  if(FComposicion[0] > 0.20)
162  FComposicion[FNumeroEspecies - 1] = 0.;
163  else
164  FComposicion[FNumeroEspecies - 1] = 1.;
165  } else {
166  if(FComposicion[0] > 0.50)
167  FComposicion[FNumeroEspecies - 1] = 1.;
168  else
169  FComposicion[FNumeroEspecies - 1] = 0.;
170  }
171  }
172  if(fracciontotal < 1. - 1e-10 || fracciontotal > 1. + 1e-10) {
173  std::cout << "ERROR: Total mass fraction must be equal to 1. Check the input data for boundary condition " << FNumeroCC
174  << std::endl;
175  throw Exception(" ");
176  }
177  if(FCalculoEspecies == nmCalculoCompleto) {
178 
179  RMezclaDep = CalculoCompletoRMezcla(FComposicion[0], FComposicion[1], FComposicion[2], 0, FCalculoGamma, nmMEP);
180  CpMezclaDep = CalculoCompletoCpMezcla(FComposicion[0], FComposicion[1], FComposicion[2], 0,
181  __units::degCToK(FTemperaturaDep), FCalculoGamma, nmMEP);
182  GammaDep = CalculoCompletoGamma(RMezclaDep, CpMezclaDep, FCalculoGamma);
183 
184  } else if(FCalculoEspecies == nmCalculoSimple) {
185 
186  RMezclaDep = CalculoSimpleRMezcla(FComposicion[0], 0, FCalculoGamma, nmMEP);
187  CvMezclaDep = CalculoSimpleCvMezcla(__units::degCToK(FTemperaturaDep), FComposicion[0], 0, FCalculoGamma, nmMEP);
188  GammaDep = CalculoSimpleGamma(RMezclaDep, CvMezclaDep, FCalculoGamma);
189 
190  }
191  FVelocidadSonidoDep = sqrt(__units::degCToK(FTemperaturaDep) * GammaDep * RMezclaDep) / __cons::ARef;
192 
193  } else if(FTipoDescarga == nmDescargaRemansoMatlab) {
194  fscanf(fich, "%lf %lf %lf ", &FPressure, &FTemperaturaDep, &FPerdidaExtremo);
195 
196  // Se determina si este remanso modela el escape del motor.
197  fscanf(fich, "%d ", &modeladoescape);
198  modeladoescape == 0 ? FModeladoEscape = false : FModeladoEscape = true;
199 
200  // Inicializacion del transporte de especies quimicas.
201  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
202  FComposicion = new double[FNumeroEspecies - FIntEGR];
203  for(int i = 0; i < FNumeroEspecies - 1; i++) {
204  fscanf(fich, "%lf ", &FComposicion[i]);
205  FFraccionMasicaEspecie[i] = FTuboExtremo[0].Pipe->GetFraccionMasicaInicial(i);
206  fracciontotal += FComposicion[i];
207  }
208  if(FHayEGR) {
209  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboExtremo[0].Pipe->GetFraccionMasicaInicial(FNumeroEspecies - 1);
210  if(FCalculoEspecies == nmCalculoCompleto) {
211  if(FComposicion[0] > 0.2)
212  FComposicion[FNumeroEspecies - 1] = 0.;
213  else
214  FComposicion[FNumeroEspecies - 1] = 1.;
215  } else {
216  if(FComposicion[0] > 0.5)
217  FComposicion[FNumeroEspecies - 1] = 1.;
218  else
219  FComposicion[FNumeroEspecies - 1] = 0.;
220  }
221  }
222  if(fracciontotal != 1.) {
223  std::cout <<
224  "ERROR: La fraccion masica total no puede ser distinta de 1. Repasa la lectura en la condicion de contorno " <<
225  FNumeroCC << std::endl;
226  throw Exception(" ");
227  }
228  if(FCalculoEspecies == nmCalculoCompleto) {
229 
230  RMezclaDep = CalculoCompletoRMezcla(FComposicion[0], FComposicion[1], FComposicion[2], 0, FCalculoGamma, nmMEP);
231  CpMezclaDep = CalculoCompletoCpMezcla(FComposicion[0], FComposicion[1], FComposicion[2], 0,
232  __units::degCToK(FTemperaturaDep), FCalculoGamma, nmMEP);
233  GammaDep = CalculoCompletoGamma(RMezclaDep, CpMezclaDep, FCalculoGamma);
234 
235  } else if(FCalculoEspecies == nmCalculoSimple) {
236 
237  RMezclaDep = CalculoSimpleRMezcla(FComposicion[0], 0, FCalculoGamma, nmMEP);
238  CvMezclaDep = CalculoSimpleCvMezcla(__units::degCToK(FTemperaturaDep), FComposicion[0], 0, FCalculoGamma, nmMEP);
239  GammaDep = CalculoSimpleGamma(RMezclaDep, CvMezclaDep, FCalculoGamma);
240 
241  }
242  FVelocidadSonidoDep = sqrt(__units::degCToK(FTemperaturaDep) * GammaDep * RMezclaDep) / __cons::ARef;
243 
244  } else
245  printf("ERROR:TCCDescargaExtremoAbierto::LeeDescargaExtremoAbierto.Asignacion Tipo BC\n");
246 
247  fgetpos(fich, &filepos);
248  fclose(fich);
249 
250  } catch(exception & N) {
251  std::cout << "ERROR: TCCDescargaExtremoAbierto::LeeDescargaExtremoAbierto en la condicion de contorno: " << FNumeroCC <<
252  std::endl;
253  std::cout << "Tipo de error: " << N.what() << std::endl;
254  throw Exception(N.what());
255  }
256 }
257 
258 // ---------------------------------------------------------------------------
259 // ---------------------------------------------------------------------------
260 
261 //void TCCDescargaExtremoAbierto::PutPresion(double valor) {
262 // try {
263 // FPressure = valor;
264 // }
265 // catch(Exception & N) {
266 // std::cout << "ERROR: TCCDescargaExtremoAbierto::PutPresion en la condicion de contorno: " <<
267 // FNumeroCC << std::endl;
268 // std::cout << "Tipo de error: " << N.what() << std::endl;
269 // throw Exception(N.what());
270 // }
271 //}
272 
273 // ---------------------------------------------------------------------------
274 // ---------------------------------------------------------------------------
275 
276 //void TCCDescargaExtremoAbierto::PutTemperatura(double valor) {
277 // try {
278 // double RMezclaDep = 0., CvMezclaDep = 0., CpMezclaDep = 0., GammaDep = 0.;
279 //
280 // FTemperaturaDep = valor;
281 // if (FCalculoEspecies == nmCalculoCompleto) {
282 //
283 // RMezclaDep = CalculoCompletoRMezcla(FComposicion[0], FComposicion[1], FComposicion[2],
284 // FCalculoGamma);
285 // CpMezclaDep = CalculoCompletoCpMezcla(FComposicion[0], FComposicion[1],
286 // FComposicion[2], FTemperaturaDep, FCalculoGamma, 0);
287 // GammaDep = CalculoCompletoGamma(RMezclaDep, CpMezclaDep, FCalculoGamma, 0);
288 //
289 // }
290 // else if (FCalculoEspecies == nmCalculoSimple) {
291 //
292 // RMezclaDep = CalculoSimpleRMezcla(FComposicion[0], FCalculoGamma, 0);
293 // CvMezclaDep = CalculoSimpleCvMezcla(FTemperaturaDep, FComposicion[0], FCalculoGamma, 0);
294 // GammaDep = CalculoSimpleGamma(RMezclaDep, CvMezclaDep, FCalculoGamma);
295 //
296 // }
297 // FVelocidadSonidoDep = sqrt(FTemperaturaDep * GammaDep * RMezclaDep) / __cons::ARef;
298 // }
299 // catch(Exception & N) {
300 // std::cout << "ERROR: TCCDescargaExtremoAbierto::PutTemperatura en la condicion de contorno: " <<
301 // FNumeroCC << std::endl;
302 // std::cout << "Tipo de error: " << N.what() << std::endl;
303 // throw Exception(N.what());
304 // }
305 //}
306 
307 // ---------------------------------------------------------------------------
308 // ---------------------------------------------------------------------------
309 
310 void TCCDescargaExtremoAbierto::CalculaCondicionContorno(double Time) {
311  try {
312  double rel_CCon_entropia, yyy, pplo, aap, xyx, b, a2, c, u_isen, a_isen, a_real, u_real;
313  double FraccionMasicaAcum = 0.;
314 
315  FGamma = FTuboExtremo[0].Pipe->GetGamma(FNodoFin);
316  FGamma3 = __Gamma::G3(FGamma);
317 
318  rel_CCon_entropia = *FCC / FTuboExtremo[0].Entropia;
319  yyy = pow(FPressure / FPref, __Gamma::G5(FGamma));
320  if(rel_CCon_entropia / yyy >= 1.000005) { // Flujo Saliente del conducto.
321  /* ________ */
322  /* caso > 1 */
323  /* ________ */
324  *FCD = FTuboExtremo[0].Entropia * 2. * yyy - *FCC;
325  pplo = __Gamma::G7(FGamma) * *FCC; /* Condicion flujo supersonico */
326  if(*FCD < pplo) {
327  *FCD = pplo;
328  }
329  // Transporte de especies quimicas.
330  for(int j = 0; j < FNumeroEspecies - 2; j++) {
331  FFraccionMasicaEspecie[j] = FTuboExtremo[0].Pipe->GetFraccionMasicaCC(FIndiceCC, j);
332  FraccionMasicaAcum += FFraccionMasicaEspecie[j];
333  }
334  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
335  if(FHayEGR)
336  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboExtremo[0].Pipe->GetFraccionMasicaCC(FIndiceCC, FNumeroEspecies - 1);
337 
338  } else if(rel_CCon_entropia / yyy < .999995) { // Flujo Entrante al conducto.
339  /* ________ */
340  /* caso < 1 */
341  /* ________ */
342  aap = FVelocidadSonidoDep / yyy;
343  xyx = aap / FTuboExtremo[0].Entropia;
344  b = __Gamma::G1(FGamma) * *FCC * pow2(xyx) * FPerdidaExtremo;
345  a2 = pow2(FGamma3 * xyx * FPerdidaExtremo) + FGamma3;
346  c = pow2(xyx * *FCC) - pow2(FVelocidadSonidoDep);
347  u_isen = (-b + sqrt(pow2(b) - a2 * 4. * c)) / (2. * a2);
348  // Resolucion ecuacion de segundo grado
349  a_isen = sqrt(pow2(FVelocidadSonidoDep) - FGamma3 * pow2(u_isen));
350  u_real = u_isen * FPerdidaExtremo; // Con esta relacion obtenemos la velocidad real.
351  a_real = sqrt(pow2(FVelocidadSonidoDep) - FGamma3 * pow2(u_real));
352  aap = a_real / a_isen * aap;
353  if(fabs(u_real) > a_real) { /* Condicion flujo supersonico */
354  a_real = sqrt(2 / __Gamma::G2(FGamma)) * FVelocidadSonidoDep;
355  u_real = a_real;
356  aap = FTuboExtremo[0].Entropia * a_real / (*FCC + FGamma3 * u_real);
357  }
358  *FCC = a_real - FGamma3 * u_real;
359  *FCD = a_real + FGamma3 * u_real;
360  FTuboExtremo[0].Entropia = aap;
361 
362  if(!FModeladoEscape) {
363  for(int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
364  FFraccionMasicaEspecie[j] = FComposicion[j];
365  }
366  }
367 
368  } else { // Flujo Parado
369  /* ________ */
370  /* caso = 1 */
371  /* ________ */
372  *FCD = *FCC;
373  // La composicion se mantiene, al estar el flujo parado.
374 
375  }
376 
377  } catch(exception & N) {
378  std::cout << "ERROR: TCCDescargaExtremoAbierto::CalculaCondicionesContorno en la condicion de contorno: " << FNumeroCC
379  << std::endl;
380  std::cout << "Tipo de error: " << N.what() << std::endl;
381  throw Exception(N.what());
382  }
383 }
384 
385 // ---------------------------------------------------------------------------
386 // ---------------------------------------------------------------------------
387 
388 #pragma package(smart_init)
TTubo
a Finite differences pipe.
Definition: TTubo.h:116
TTubo::GetGamma
double GetGamma(int i) const
Gets the specific heat capacities ratio at a given cell.
Definition: TTubo.cpp:5444
stTuboExtremo
Definition: Globales.h:730
TDPF
Definition: TDPF.h:45
TTubo::GetFraccionMasicaCC
double GetFraccionMasicaCC(int j, int i)
Definition: TTubo.h:953
TTubo::GetFraccionMasicaInicial
double GetFraccionMasicaInicial(int i) const
Gets the initial mass fraction of species i.
Definition: TTubo.cpp:5440
TCondicionContorno
Definition: TCondicionContorno.h:54
Exception
Custom exception class.
Definition: Exception.hpp:39
TTubo.h
pow2
T pow2(T x)
Returns x to the power of 2.
Definition: Math_wam.h:88