OpenWAM
TCCUnionEntreTubos.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 "TCCUnionEntreTubos.h"
32 //#include <cmath>
33 #include <iostream>
34 #include "TTubo.h"
35 
36 // ---------------------------------------------------------------------------
37 // ---------------------------------------------------------------------------
38 
39 TCCUnionEntreTubos::TCCUnionEntreTubos(nmTypeBC TipoCC, int numCC, nmTipoCalculoEspecies SpeciesModel,
40  int numeroespecies, nmCalculoGamma GammaCalculation, bool ThereIsEGR) :
41  TCondicionContorno(TipoCC, numCC, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
42 
43  FTuboExtremo = NULL;
44  FNodoFin = NULL;
45  FIndiceCC = NULL;
46  FNumeroTubo = NULL;
47  FTubo = NULL;
48  FCC = NULL;
49  FCD = NULL;
50 
51 }
52 // ---------------------------------------------------------------------------
53 // ---------------------------------------------------------------------------
54 
55 TCCUnionEntreTubos::~TCCUnionEntreTubos() {
56 
57  if(FTuboExtremo != NULL)
58  delete[] FTuboExtremo;
59  if(FNodoFin != NULL)
60  delete[] FNodoFin;
61  if(FIndiceCC != NULL)
62  delete[] FIndiceCC;
63  if(FNumeroTubo != NULL)
64  delete[] FNumeroTubo;
65  if(FTubo != NULL)
66  delete[] FTubo;
67  if(FCC != NULL)
68  delete[] FCC;
69  if(FCD != NULL)
70  delete[] FCD;
71 
72 }
73 
74 // ---------------------------------------------------------------------------
75 // ---------------------------------------------------------------------------
76 
77 void TCCUnionEntreTubos::ReadBoundaryData(const char *FileWAM, fpos_t &filepos, int NumberOfPipes, TTubo **Pipe,
78  int nDPF, TDPF **DPF) {
79  try {
80  int i = 0;
81 
82  FEspesor = 0.;
83  FConductividad = 0.;
84 
85  FTuboExtremo = new stTuboExtremo[2];
86  FNodoFin = new int[2];
87  FIndiceCC = new int[2];
88  FCC = new double*[2];
89  FCD = new double*[2];
90  FNumeroTubo = new int[2];
91  FTubo = new int[2];
92 
93  for(int i = 0; i < 2; i++) {
94  FTuboExtremo[i].Pipe = NULL;
95  }
96 
97  while(FNumeroTubosCC < 2 && i < NumberOfPipes) {
98  if(Pipe[i]->getNodoIzq() == FNumeroCC || Pipe[i]->getNodoDer() == FNumeroCC) {
99  FTubo[FNumeroTubosCC] = i;
100  if(Pipe[FTubo[FNumeroTubosCC]]->getNodoIzq() == FNumeroCC) {
101  FNodoFin[FNumeroTubosCC] = 0;
102  }
103  if(Pipe[FTubo[FNumeroTubosCC]]->getNodoDer() == FNumeroCC) {
104  FNodoFin[FNumeroTubosCC] = Pipe[FTubo[FNumeroTubosCC]]->getNin() - 1;
105  }
106  FNumeroTubosCC++;
107  }
108  i++;
109  }
110 
111  /* Ahora al tubo de mayor diametro se le asignara la posicion 1 de los vectores
112  y al de menor diametro la posicion 0 */
113  if(Pipe[FTubo[0]]->GetDiametro(FNodoFin[0]) >= Pipe[FTubo[1]]->GetDiametro(FNodoFin[1])) {
114  if(Pipe[FTubo[0]]->getNodoIzq() == FNumeroCC) {
115  FTuboExtremo[1].Pipe = Pipe[FTubo[0]];
116  FTuboExtremo[1].TipoExtremo = nmLeft;
117  FNodoFin[1] = 0;
118  FIndiceCC[1] = 0;
119  FNumeroTubo[1] = FTubo[0];
120  FCC[1] = &(FTuboExtremo[1].Beta);
121  FCD[1] = &(FTuboExtremo[1].Landa);
122  }
123  if(Pipe[FTubo[0]]->getNodoDer() == FNumeroCC) {
124  FTuboExtremo[1].Pipe = Pipe[FTubo[0]];
125  FTuboExtremo[1].TipoExtremo = nmRight;
126  FNodoFin[1] = Pipe[FTubo[0]]->getNin() - 1;
127  FIndiceCC[1] = 1;
128  FNumeroTubo[1] = FTubo[0];
129  FCC[1] = &(FTuboExtremo[1].Landa);
130  FCD[1] = &(FTuboExtremo[1].Beta);
131  }
132  if(Pipe[FTubo[1]]->getNodoIzq() == FNumeroCC) {
133  FTuboExtremo[0].Pipe = Pipe[FTubo[1]];
134  FTuboExtremo[0].TipoExtremo = nmLeft;
135  FNodoFin[0] = 0;
136  FIndiceCC[0] = 0;
137  FNumeroTubo[0] = FTubo[1];
138  FCC[0] = &(FTuboExtremo[0].Beta);
139  FCD[0] = &(FTuboExtremo[0].Landa);
140  }
141  if(Pipe[FTubo[1]]->getNodoDer() == FNumeroCC) {
142  FTuboExtremo[0].Pipe = Pipe[FTubo[1]];
143  FTuboExtremo[0].TipoExtremo = nmRight;
144  FNodoFin[0] = Pipe[FTubo[1]]->getNin() - 1;
145  FIndiceCC[0] = 1;
146  FNumeroTubo[0] = FTubo[1];
147  FCC[0] = &(FTuboExtremo[0].Landa);
148  FCD[0] = &(FTuboExtremo[0].Beta);
149  }
150 
151  } else {
152  if(Pipe[FTubo[1]]->getNodoIzq() == FNumeroCC) {
153  FTuboExtremo[1].Pipe = Pipe[FTubo[1]];
154  FTuboExtremo[1].TipoExtremo = nmLeft;
155  FNodoFin[1] = 0;
156  FIndiceCC[1] = 0;
157  FNumeroTubo[1] = FTubo[1];
158  FCC[1] = &(FTuboExtremo[1].Beta);
159  FCD[1] = &(FTuboExtremo[1].Landa);
160  }
161  if(Pipe[FTubo[1]]->getNodoDer() == FNumeroCC) {
162  FTuboExtremo[1].Pipe = Pipe[FTubo[1]];
163  FTuboExtremo[1].TipoExtremo = nmRight;
164  FNodoFin[1] = Pipe[FTubo[1]]->getNin() - 1;
165  FIndiceCC[1] = 1;
166  FNumeroTubo[1] = FTubo[1];
167  FCC[1] = &(FTuboExtremo[1].Landa);
168  FCD[1] = &(FTuboExtremo[1].Beta);
169  }
170  if(Pipe[FTubo[0]]->getNodoIzq() == FNumeroCC) {
171  FTuboExtremo[0].Pipe = Pipe[FTubo[0]];
172  FTuboExtremo[0].TipoExtremo = nmLeft;
173  FNodoFin[0] = 0;
174  FIndiceCC[0] = 0;
175  FNumeroTubo[0] = FTubo[0];
176  FCC[0] = &(FTuboExtremo[0].Beta);
177  FCD[0] = &(FTuboExtremo[0].Landa);
178  }
179  if(Pipe[FTubo[0]]->getNodoDer() == FNumeroCC) {
180  FTuboExtremo[0].Pipe = Pipe[FTubo[0]];
181  FTuboExtremo[0].TipoExtremo = nmRight;
182  FNodoFin[0] = Pipe[FTubo[0]]->getNin() - 1;
183  FIndiceCC[0] = 1;
184  FNumeroTubo[0] = FTubo[0];
185  FCC[0] = &(FTuboExtremo[0].Landa);
186  FCD[0] = &(FTuboExtremo[0].Beta);
187  }
188  }
189 
190  // Inicializacion del transporte de especies quimicas.
191  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
192  for(int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
193  // Se elige como composicion inicial la del tubo 0. Es arbitrario.
194  FFraccionMasicaEspecie[i] = FTuboExtremo[0].Pipe->GetFraccionMasicaInicial(i);
195  }
196 
197  FILE *fich = fopen(FileWAM, "r");
198  fsetpos(fich, &filepos);
199 
200  fscanf(fich, "%lf %lf", &FEspesor, &FConductividad);
201  /* Coeficiente de perdidas con signo positivo */
202 
203  fgetpos(fich, &filepos);
204  fclose(fich);
205 
206  } catch(exception & N) {
207  std::cout << "ERROR: TCCUnionEntreTubos::LeeUnionEntreTubos en la condicion de contorno: " << FNumeroCC << std::endl;
208  std::cout << "Tipo de error: " << N.what() << std::endl;
209  throw Exception(N.what());
210  }
211 }
212 
213 // ---------------------------------------------------------------------------
214 // ---------------------------------------------------------------------------
215 
216 void TCCUnionEntreTubos::TuboCalculandose(int TuboActual) {
217  try {
218  FTuboActual = TuboActual;
219  } catch(exception & N) {
220  std::cout << "ERROR: TCCUnionEntreTubos::TuboCalculandose en la condicion de contorno: " << FNumeroCC << std::endl;
221  std::cout << "Tipo de error: " << N.what() << std::endl;
222  throw Exception(N.what());
223  }
224 }
225 
226 // ---------------------------------------------------------------------------
227 // ---------------------------------------------------------------------------
228 
229 void TCCUnionEntreTubos::CalculaCondicionContorno(double Time) {
230  try {
231  double rel_entropia = 0., rel_area = 0., vel_sonido_In = 0., vel_sonido_Out = 0., vel_In = 0., vel_Out = 0.,
232  correc_sonido_In = 0.;
233  double flujo, FraccionMasicaAcum = 0., exd, exi;
234  int TuboCalculado = 0;
235 
236  if(FTuboActual == 10000) {
237  TuboCalculado = FTuboActual;
238  FGamma = FTuboExtremo[0].Pipe->GetGamma(FNodoFin[0]);
239  } else {
240  for(int i = 0; i < FNumeroTubosCC; i++) {
241  if(FNumeroTubo[i] == FTuboActual) {
242  TuboCalculado = i;
243  }
244  }
245  FGamma = FTuboExtremo[TuboCalculado].Pipe->GetGamma(FNodoFin[TuboCalculado]);
246  }
247 
248  FGamma3 = __Gamma::G3(FGamma);
249  FGamma2 = __Gamma::G2(FGamma);
250  FGamma1 = __Gamma::G1(FGamma);
251 
252  /* Criterio para determinar el sentido el flujo */
253  flujo = (*FCC[1] / FTuboExtremo[1].Entropia) / (*FCC[0] / FTuboExtremo[0].Entropia);
254  if(flujo < 0.999995) { /* Sentido del flujo: de 0(saliente (out)) a 1(entrante (in)) */
255  rel_entropia = FTuboExtremo[0].Entropia / FTuboExtremo[1].Entropia;
256  rel_area = pow2(FTuboExtremo[1].Pipe->GetDiametro(FNodoFin[1]) / FTuboExtremo[0].Pipe->GetDiametro(FNodoFin[0]));
257 
258  int cont = 0;
259  /* Intervalo de acotacion de A1 */
260  exd = *FCC[0];
261  exi = *FCC[0] * 2. / FGamma2;
262 
263  /* LLAMADA A LA ESTRUCTURA-ensanchamiento */
264  stExpansion EnsA1(*FCC[0], *FCC[1], rel_entropia, rel_area, FGamma);
265  vel_sonido_Out = FindRoot(EnsA1, exi, exd);
266  vel_In = EnsA1.U2;
267  vel_Out = EnsA1.U1;
268  vel_sonido_In = EnsA1.A2;
269  correc_sonido_In = EnsA1.xx3;
270  /* nuevo if (abs(vel_sonido_Out-vel_Out)<1E-12) {
271  printf ("");
272  } */
273  if(TuboCalculado == 1) {
274  *FCD[1] = vel_sonido_In + FGamma3 * vel_In;
275  *FCC[1] = vel_sonido_In - FGamma3 * vel_In;
276  FTuboExtremo[1].Entropia = vel_sonido_In * FTuboExtremo[1].Entropia / correc_sonido_In;
277 
278  } else if(TuboCalculado == 0) {
279  *FCD[0] = vel_sonido_Out - FGamma3 * vel_Out;
280 
281  } else if(TuboCalculado == 10000) {
282  *FCD[1] = vel_sonido_In + FGamma3 * vel_In;
283  *FCC[1] = vel_sonido_In - FGamma3 * vel_In;
284  FTuboExtremo[1].Entropia = vel_sonido_In * FTuboExtremo[1].Entropia / correc_sonido_In;
285  *FCD[0] = vel_sonido_Out - FGamma3 * vel_Out;
286  }
287 
288  // Transporte de Especies Quimicas
289  // Se actualiza todos los instantes de calculo.
290  for(int j = 0; j < FNumeroEspecies - 2; j++) {
291  FFraccionMasicaEspecie[j] = FTuboExtremo[0].Pipe->GetFraccionMasicaCC(FIndiceCC[0], j);
292  FraccionMasicaAcum += FFraccionMasicaEspecie[j];
293  }
294  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
295  if(FHayEGR)
296  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboExtremo[0].Pipe->GetFraccionMasicaCC(FIndiceCC[0],
297  FNumeroEspecies - 1);
298 
299  } else if(flujo > 1.000005) { /* Sentido del flujo: de 1(saliente) a 0(entrante) */
300 
301  rel_entropia = FTuboExtremo[0].Entropia / FTuboExtremo[1].Entropia;
302  rel_area = pow2(FTuboExtremo[0].Pipe->GetDiametro(FNodoFin[0]) / FTuboExtremo[1].Pipe->GetDiametro(FNodoFin[1]));
303 
304  /* LLAMADA A LA ESTRUCTURA-estrechamiento */
305 
306  /* Intervalo de acotacion de U1 */
307  exi = 0;
308  exd = *FCC[1] * 2. / FGamma2;
309 
310  stContraction EstU1(*FCC[1], *FCC[0], rel_entropia, rel_area, FGamma);
311  vel_Out = FindRoot(EstU1, exi, exd);
312  vel_sonido_In = EstU1.A2;
313  vel_sonido_Out = EstU1.A1;
314  vel_In = EstU1.U2;
315 
316  if(TuboCalculado == 0) {
317  *FCD[0] = vel_sonido_In + FGamma3 * vel_In;
318  *FCC[0] = vel_sonido_In - FGamma3 * vel_In;
319  FTuboExtremo[0].Entropia = FTuboExtremo[1].Entropia;
320 
321  } else if(TuboCalculado == 1) {
322  *FCD[1] = vel_sonido_Out - FGamma3 * vel_Out;
323 
324  } else if(TuboCalculado == 10000) {
325  *FCD[0] = vel_sonido_In + FGamma3 * vel_In;
326  *FCC[0] = vel_sonido_In - FGamma3 * vel_In;
327  FTuboExtremo[0].Entropia = FTuboExtremo[1].Entropia;
328  *FCD[1] = vel_sonido_Out - FGamma3 * vel_Out;
329  }
330 
331  // Transporte de Especies Quimicas
332  // Se actualiza todos los instantes de calculo (al igual que la temperatura en la BC).
333  for(int j = 0; j < FNumeroEspecies - 2; j++) {
334  FFraccionMasicaEspecie[j] = FTuboExtremo[1].Pipe->GetFraccionMasicaCC(FIndiceCC[1], j);
335  FraccionMasicaAcum += FFraccionMasicaEspecie[j];
336  }
337  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
338  if(FHayEGR)
339  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboExtremo[1].Pipe->GetFraccionMasicaCC(FIndiceCC[1],
340  FNumeroEspecies - 1);
341 
342  } else { /* Flujo Parado */
343 
344  if(TuboCalculado == 0) {
345  *FCD[0] = *FCC[0];
346  } else if(TuboCalculado == 1) {
347  *FCD[1] = *FCC[1];
348  } else if(TuboCalculado == 10000) {
349  *FCD[0] = *FCC[0];
350  *FCD[1] = *FCC[1];
351  }
352  // La composicion se mantiene, al estar el flujo parado.
353 
354  }
355  } catch(exception & N) {
356  std::cout << "ERROR: TCCUnionEntreTubos::CalculaCondicionContorno en la condicion de contorno: " << FNumeroCC <<
357  std::endl;
358  std::cout << "Tipo de error: " << N.what() << std::endl;
359  throw Exception(N.what());
360  }
361 }
362 
363 // ---------------------------------------------------------------------------
364 // ---------------------------------------------------------------------------
365 
366 /* void TCCUnionEntreTubos::Estrechamiento(double CCS,double CCE,double rel_entropia,
367  double rel_area,double smag,double *xa1,double *xa2,double *xu2,double *xu1)
368  {
369  try
370  {
371  double xx, xx1, xx2, xx3, exd, exi,ex1,ex2,xxx, xu1p, ytty, xxxx;
372 
373 
374 
375  }
376  catch(Exception &N)
377  {
378  std::cout << "ERROR: TCCUnionEntreTubos::Estrechamiento en la condicion de contorno: " << FNumeroCC << std::endl;
379  std::cout << "Tipo de error: " << N.what() << std::endl;
380  throw Exception(N.what());
381  }
382  }
383 
384  //---------------------------------------------------------------------------
385  //---------------------------------------------------------------------------
386 
387  void TCCUnionEntreTubos::Ensanchamiento(double CCS,double CCE,double rel_entropia,
388  double rel_area,double )
389  {
390  try
391  {
392 
393 
394  }
395  catch(Exception &N)
396  {
397  std::cout << "ERROR: TCCUnionEntreTubos::Ensanchamiento en la condicion de contorno: " << FNumeroCC << std::endl;
398  std::cout << "Tipo de error: " << N.what() << std::endl;
399  throw Exception(N.what());
400  }
401  }
402 
403  //void TCCUnionEntreTubos::CalculaCaracteristicas(double Time)
404  //{
405  //int signo = 0;
406  //
407  // for(int i=0;i<FNumeroTubosCC;i++){
408  // signo=1.;
409  // if(FTuboExtremo[i].TipoExtremo==nmRight) signo=-1;
410  // signo=FTuboExtremo[i].Pipe->getNumeroTubo();
411  // //if(FTuboExtremo[i].Pipe->GetVelocidad(0)*signo<=0){
412  // //FTuboExtremo[i].Entropia=FTuboExtremo[i].Pipe
413  // //}
414  // }
415  //}
416 
417 
418  //---------------------------------------------------------------------------
419  //---------------------------------------------------------------------------
420 
421 
422  #pragma package(smart_init)
423  */
stExpansion
Definition: BoundaryFunctions.h:203
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
stTuboExtremo
Definition: Globales.h:730
TDPF
Definition: TDPF.h:45
TTubo::GetDiametro
double GetDiametro(int i) const
Gets the cell diameter.
Definition: TTubo.cpp:5436
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
stContraction
Definition: BoundaryFunctions.h:271
TTubo.h
pow2
T pow2(T x)
Returns x to the power of 2.
Definition: Math_wam.h:88