OpenWAM
TCCPerdidadePresion.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 "TCCPerdidadePresion.h"
32 
33 //#include <cmath>
34 #ifdef __BORLANDC__
35 #include <vcl.h>
36 #endif
37 #include "TTubo.h"
38 
39 //---------------------------------------------------------------------------
40 //---------------------------------------------------------------------------
41 
42 TCCPerdidadePresion::TCCPerdidadePresion(nmTypeBC TipoCC, int numCC, nmTipoCalculoEspecies SpeciesModel,
43  int numeroespecies, nmCalculoGamma GammaCalculation, bool ThereIsEGR) :
44  TCondicionContorno(TipoCC, numCC, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
45 
46  FTuboExtremo = NULL;
47  FCC = NULL;
48  FCD = NULL;
49  FNodoFin = NULL;
50  FIndiceCC = NULL;
51  FNumeroTubo = NULL;
52 
53  if(TipoCC == nmLinearPressureLoss)
54  FTipoPP = nmPPLineal;
55  else if(TipoCC == nmQuadraticPressureLoss)
56  FTipoPP = nmPPCuadratica;
57  else
58  printf("ERROR en tipo de perdida de presion TCCPerdidadePresion en la condicion de contorno: %d\n", FNumeroCC);
59 
60 }
61 
62 //---------------------------------------------------------------------------
63 //---------------------------------------------------------------------------
64 
65 TCCPerdidadePresion::~TCCPerdidadePresion() {
66  delete[] FTuboExtremo;
67 
68  if(FNodoFin != NULL)
69  delete[] FNodoFin;
70  if(FIndiceCC != NULL)
71  delete[] FIndiceCC;
72  if(*FCC != NULL)
73  delete[] FCC;
74  if(*FCD != NULL)
75  delete[] FCD;
76  if(FNumeroTubo != NULL)
77  delete[] FNumeroTubo;
78 }
79 
80 //---------------------------------------------------------------------------
81 //---------------------------------------------------------------------------
82 
83 void TCCPerdidadePresion::ReadBoundaryData(const char *FileWAM, fpos_t &filepos, int NumberOfPipes, TTubo **Pipe,
84  int nDPF, TDPF **DPF) {
85  try {
86  int i = 0;
87  FK = 0;
88 
89  FTuboExtremo = new stTuboExtremo[2];
90  FNodoFin = new int[2];
91  FIndiceCC = new int[2];
92  FCC = new double*[2];
93  FCD = new double*[2];
94  FNumeroTubo = new int[2];
95 
96  for(int i = 0; i < 2; i++) {
97  FTuboExtremo[i].Pipe = NULL;
98  }
99 
100  while(FNumeroTubosCC < 2 && i < NumberOfPipes) {
101  if(Pipe[i]->getNodoIzq() == FNumeroCC) {
102  FTuboExtremo[FNumeroTubosCC].Pipe = Pipe[i];
103  FTuboExtremo[FNumeroTubosCC].TipoExtremo = nmLeft;
104  FNodoFin[FNumeroTubosCC] = 0;
105  FIndiceCC[FNumeroTubosCC] = 0;
106  FNumeroTubo[FNumeroTubosCC] = Pipe[i]->getNumeroTubo() - 1;
107  FCC[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Beta);
108  FCD[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Landa);
109  FNumeroTubosCC++;
110  }
111  if(Pipe[i]->getNodoDer() == FNumeroCC) {
112  FTuboExtremo[FNumeroTubosCC].Pipe = Pipe[i];
113  FTuboExtremo[FNumeroTubosCC].TipoExtremo = nmRight;
114  FNodoFin[FNumeroTubosCC] = Pipe[i]->getNin() - 1;
115  FIndiceCC[FNumeroTubosCC] = 1;
116  FNumeroTubo[FNumeroTubosCC] = Pipe[i]->getNumeroTubo() - 1;
117  FCC[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Landa);
118  FCD[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Beta);
119  FNumeroTubosCC++;
120  }
121  i++;
122  }
123 
124 // Inicializacion del transporte de especies quimicas.
125  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
126  for(int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
127  // Se elige como composicion inicial la del tubo 0. Es arbitrario.
128  FFraccionMasicaEspecie[i] = FTuboExtremo[0].Pipe->GetFraccionMasicaInicial(i);
129  }
130 
131  FILE *fich = fopen(FileWAM, "r");
132  fsetpos(fich, &filepos);
133 
134  fscanf(fich, "%lf ", &FK); /* Coeficiente de perdidas con signo positivo */
135 
136  fgetpos(fich, &filepos);
137  fclose(fich);
138 
139  } catch(exception &N) {
140  std::cout << "ERROR: TCCPerdidadePresion::LecturaPerdidaPresion en la condicion de contorno: " << FNumeroCC <<
141  std::endl;
142  std::cout << "Tipo de error: " << N.what() << std::endl;
143  throw Exception(N.what());
144  }
145 }
146 
147 //---------------------------------------------------------------------------
148 //---------------------------------------------------------------------------
149 
150 void TCCPerdidadePresion::TuboCalculandose(int TuboActual) {
151  try {
152  FTuboActual = TuboActual;
153  } catch(exception &N) {
154  std::cout << "ERROR: TCCPerdidadePresion::TuboCalculandose en la condicion de contorno: " << FNumeroCC << std::endl;
155  std::cout << "Tipo de error: " << N.what() << std::endl;
156  throw Exception(N.what());
157  }
158 }
159 
160 //---------------------------------------------------------------------------
161 //---------------------------------------------------------------------------
162 
163 void TCCPerdidadePresion::CalculaCondicionContorno(double Time) {
164  try {
165 
166  double vel_sonido_Out = 0., vel_Out = 0., vel_sonido_In = 0., vel_In = 0., xx3 = 0., ei = 0, ed = 0;
167  double flujo, FraccionMasicaAcum = 0.;
168  int TuboCalculado = 0;
169 
170  if(FTuboActual == 10000) {
171  TuboCalculado = FTuboActual;
172  FGamma = FTuboExtremo[0].Pipe->GetGamma(FNodoFin[0]);
173  } else {
174  for(int i = 0; i < FNumeroTubosCC; i++) {
175  if(FNumeroTubo[i] == FTuboActual) {
176  TuboCalculado = i;
177  }
178  }
179  FGamma = FTuboExtremo[TuboCalculado].Pipe->GetGamma(FNodoFin[TuboCalculado]);
180  }
181  FGamma1 = __Gamma::G1(FGamma);
182  FGamma3 = __Gamma::G3(FGamma);
183  FGamma2 = __Gamma::G2(FGamma);
184  FGamma5 = __Gamma::G5(FGamma);
185 
186  flujo = (*FCC[1] / FTuboExtremo[1].Entropia) / (*FCC[0] / FTuboExtremo[0].Entropia);
187 
188  if(flujo < .999995) { /* Flujo de 0 (Saliente) a 1 (Entrante) */
189 
190  FRelacionEntropia = FTuboExtremo[0].Entropia / FTuboExtremo[1].Entropia;
191 
192  /*Acotacion del intervalo de busqueda para A1 (Velocity del sonido a la entrada)*/
193  if((*FCC[0] * 2 / FGamma2) < (*FCC[0] / (FGamma3 / sqrt(FK) + 1) + 1e-6)) {
194  ei = *FCC[0] / (FGamma3 / sqrt(FK) + 1) + 1e-6;
195  } else {
196  ei = *FCC[0] * 2 / FGamma2;
197  }
198  ed = *FCC[0];
199 
200  if(FTipoPP == nmPPLineal) { /* Perdida lineal */
201 
202  ei = QuadraticEqP(FGamma1, 2 * FK, -2 * FK * *FCC[0]) + 1e-6;
203  ed = *FCC[0];
204 
205  stPerdPresAdL PPAL(*FCC[0], *FCC[1], FK, FGamma, FRelacionEntropia, __cons::ARef);
206  vel_sonido_Out = FindRoot(PPAL, ei, ed);
207  vel_sonido_In = PPAL.A2;
208  vel_In = PPAL.U2;
209  vel_Out = PPAL.U1;
210  xx3 = PPAL.xx3;
211 
212  } else if(FTipoPP == nmPPCuadratica) { /* Perdida cuadratica */
213 
214  stPerdPresAd PPA(*FCC[0], *FCC[1], FK, FGamma, FRelacionEntropia);
215  vel_sonido_Out = FindRoot(PPA, ei, ed);
216  vel_sonido_In = PPA.A2;
217  vel_In = PPA.U2;
218  vel_Out = PPA.U1;
219  xx3 = PPA.xx3;
220  /*nuevo if (abs(vel_sonido_Out-vel_Out)<1E-12) {
221  printf ("");
222  } */
223 
224  if(PPA.U1 * PPA.U2 < 0) {
225  PPA.U1 = 0;
226  PPA.U2 = 0;
227  }
228 
229  } else
230  printf("Error en el tipo de perdida de presion en la condicion de contorno: %d\n", FNumeroCC);
231 
232  if(TuboCalculado == 1) {
233  *FCC[1] = vel_sonido_In - FGamma3 * vel_In;
234  *FCD[1] = vel_sonido_In + FGamma3 * vel_In;
235  FTuboExtremo[1].Entropia = FTuboExtremo[1].Entropia * vel_sonido_In / xx3;
236 
237  } else if(TuboCalculado == 0) {
238  *FCC[0] = vel_sonido_Out + FGamma3 * vel_Out;
239  *FCD[0] = vel_sonido_Out - FGamma3 * vel_Out;
240 
241  } else if(TuboCalculado == 10000) {
242  *FCC[1] = vel_sonido_In - FGamma3 * vel_In;
243  *FCD[1] = vel_sonido_In + FGamma3 * vel_In;
244  FTuboExtremo[1].Entropia = FTuboExtremo[1].Entropia * vel_sonido_In / xx3;
245  *FCC[0] = vel_sonido_Out + FGamma3 * vel_Out;
246  *FCD[0] = vel_sonido_Out - FGamma3 * vel_Out;
247 
248  }
249 
250  // Chemical species transport
251  for(int j = 0; j < FNumeroEspecies - 2; j++) {
252  FFraccionMasicaEspecie[j] = FTuboExtremo[0].Pipe->GetFraccionMasicaCC(FIndiceCC[0], j);
253  FraccionMasicaAcum += FFraccionMasicaEspecie[j];
254  }
255  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
256  if(FHayEGR)
257  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboExtremo[0].Pipe->GetFraccionMasicaCC(FIndiceCC[0],
258  FNumeroEspecies - 1);
259  } else if(flujo >= 1.000005) { /* Flujo de 1 (Saliente) a 0 (Entrante) */
260 
261  FRelacionEntropia = FTuboExtremo[1].Entropia / FTuboExtremo[0].Entropia;
262 
263  if((*FCC[1] * 2 / FGamma2) < (*FCC[1] / (FGamma3 / sqrt(FK) + 1) + 1e-6)) {
264  ei = *FCC[1] / (FGamma3 / sqrt(FK) + 1) + 1e-6;
265  } else {
266  ei = *FCC[1] * 2 / FGamma2;
267  }
268  ed = *FCC[1];
269 
270  if(FTipoPP == nmPPLineal) { /* Linear pressure loss */
271 
272  ei = QuadraticEqP(FGamma1, 2 * FK, -2 * FK * *FCC[1]) + 1e-6;
273  ed = *FCC[1];
274 
275  stPerdPresAdL PPAL(*FCC[1], *FCC[0], FK, FGamma, FRelacionEntropia, __cons::ARef);
276  vel_sonido_Out = FindRoot(PPAL, ei, ed);
277  vel_sonido_In = PPAL.A2;
278  vel_In = PPAL.U1;
279  vel_Out = PPAL.U2;
280  xx3 = PPAL.xx3;
281 
282  } else if(FTipoPP == nmPPCuadratica) { /* Quadratic pressure loss */
283  stPerdPresAd PPA(*FCC[1], *FCC[0], FK, FGamma, FRelacionEntropia);
284  vel_sonido_Out = FindRoot(PPA, ei, ed);
285  vel_sonido_In = PPA.A2;
286  vel_In = PPA.U1;
287  vel_Out = PPA.U2;
288  xx3 = PPA.xx3;
289 
290  } else
291  printf("Error en el tipo de perdida de presion en la condicion de contorno: %d\n", FNumeroCC);
292 
293  if(TuboCalculado == 0) {
294  *FCC[0] = vel_sonido_In - FGamma3 * vel_In;
295  *FCD[0] = vel_sonido_In + FGamma3 * vel_In;
296  FTuboExtremo[0].Entropia = FTuboExtremo[0].Entropia * vel_sonido_In / xx3;
297 
298  } else if(TuboCalculado == 1) {
299  *FCC[1] = vel_sonido_Out + FGamma3 * vel_Out;
300  *FCD[1] = vel_sonido_Out - FGamma3 * vel_Out;
301 
302  } else if(TuboCalculado == 10000) {
303  *FCC[0] = vel_sonido_In - FGamma3 * vel_In;
304  *FCD[0] = vel_sonido_In + FGamma3 * vel_In;
305  FTuboExtremo[0].Entropia = FTuboExtremo[0].Entropia * vel_sonido_In / xx3;
306  *FCC[1] = vel_sonido_Out + FGamma3 * vel_Out;
307  *FCD[1] = vel_sonido_Out - FGamma3 * vel_Out;
308  }
309 
310  // Transporte de Especies Quimicas
311  // Se actualiza todos los instantes de calculo.
312  for(int j = 0; j < FNumeroEspecies - 2; j++) {
313  FFraccionMasicaEspecie[j] = FTuboExtremo[1].Pipe->GetFraccionMasicaCC(FIndiceCC[1], j);
314  FraccionMasicaAcum += FFraccionMasicaEspecie[j];
315  }
316  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
317  if(FHayEGR)
318  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboExtremo[1].Pipe->GetFraccionMasicaCC(FIndiceCC[1],
319  FNumeroEspecies - 1);
320 
321  } else {
322  if(TuboCalculado == 1) {
323  *FCD[1] = *FCC[1];
324  } else if(TuboCalculado == 0) {
325  *FCD[0] = *FCC[0];
326  } else if(TuboCalculado == 10000) {
327  *FCD[0] = *FCC[0];
328  *FCD[1] = *FCC[1];
329  }
330 
331  // La composicion se mantiene, al estar el flujo parado.
332  }
333  }
334 
335  catch(exception &N) {
336  std::cout << "ERROR: TCCPerdidadePresion::CalculaCondicionContorno en la condicion de contorno: " << FNumeroCC <<
337  std::endl;
338  std::cout << "Tipo de error: " << N.what() << std::endl;
339  throw Exception(N.what());
340  }
341 }
342 
343 //---------------------------------------------------------------------------
344 //---------------------------------------------------------------------------
345 
346 #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
FindRoot
double FindRoot(T &func, const double x1, const double x2)
Finds the root of a function.
Definition: Math_wam.h:459
stPerdPresAd
Definition: BoundaryFunctions.h:329
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
stPerdPresAdL
Definition: BoundaryFunctions.h:385
Exception
Custom exception class.
Definition: Exception.hpp:39
TTubo.h
TTubo::getNumeroTubo
int getNumeroTubo() const
Gets the pipe id.
Definition: TTubo.cpp:5464