OpenWAM
TCCRamificacion.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 "TCCRamificacion.h"
32 //#include <cmath>
33 #include <iostream>
34 #include "TTubo.h"
35 
36 // ---------------------------------------------------------------------------
37 // ---------------------------------------------------------------------------
38 
39 TCCRamificacion::TCCRamificacion(nmTypeBC TipoCC, int numCC, nmTipoCalculoEspecies SpeciesModel, int numeroespecies,
40  nmCalculoGamma GammaCalculation, bool ThereIsEGR) :
41  TCondicionContorno(TipoCC, numCC, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
42 
43  FTuboExtremo = NULL;
44 
45  FNodoFin = NULL;
46  FIndiceCC = NULL;
47  FEntropia = NULL;
48  FSeccionTubo = NULL;
49  FVelocity = NULL;
50  FDensidad = NULL;
51  FNumeroTubo = NULL;
52 
53  FCC = NULL;
54  FCD = NULL;
55 
56  FMasaEspecie = NULL;
57 
58 }
59 // ---------------------------------------------------------------------------
60 // ---------------------------------------------------------------------------
61 
62 TCCRamificacion::~TCCRamificacion() {
63 
64  if(FTuboExtremo != NULL)
65  delete[] FTuboExtremo;
66  if(FNodoFin != NULL)
67  delete[] FNodoFin;
68  if(FIndiceCC != NULL)
69  delete[] FIndiceCC;
70  if(FEntropia != NULL)
71  delete[] FEntropia;
72  if(FSeccionTubo != NULL)
73  delete[] FSeccionTubo;
74  if(FVelocity != NULL)
75  delete[] FVelocity;
76  if(FDensidad != NULL)
77  delete[] FDensidad;
78  if(FNumeroTubo != NULL)
79  delete[] FNumeroTubo;
80 
81  if(FCC != NULL)
82  delete[] FCC;
83  if(FCD != NULL)
84  delete[] FCD;
85 
86  if(FMasaEspecie != NULL)
87  delete[] FMasaEspecie;
88 
89 }
90 
91 // ---------------------------------------------------------------------------
92 // ---------------------------------------------------------------------------
93 
94 void TCCRamificacion::AsignaTubos(int NumberOfPipes, TTubo **Pipe) {
95  try {
96  int i = 0;
97  int ContadorTubosRamificacion = 0;
98 
99  ContadorTubosRamificacion = 0;
100 
101  for(int i = 0; i < NumberOfPipes; i++) {
102  if(Pipe[i]->getNodoIzq() == FNumeroCC || Pipe[i]->getNodoDer() == FNumeroCC) {
103  ContadorTubosRamificacion++;
104  }
105  }
106 
107  FTuboExtremo = new stTuboExtremo[ContadorTubosRamificacion];
108  FNodoFin = new int[ContadorTubosRamificacion];
109  FIndiceCC = new int[ContadorTubosRamificacion];
110  FCC = new double*[ContadorTubosRamificacion];
111  FCD = new double*[ContadorTubosRamificacion];
112  FEntropia = new double[ContadorTubosRamificacion];
113  FSeccionTubo = new double[ContadorTubosRamificacion];
114  FVelocity = new double[ContadorTubosRamificacion];
115  FDensidad = new double[ContadorTubosRamificacion];
116  FNumeroTubo = new int[ContadorTubosRamificacion];
117 
118  for(int i = 0; i < ContadorTubosRamificacion; i++) {
119  FTuboExtremo[i].Pipe = NULL;
120  FVelocity[i] = 0;
121  }
122 
123  while(FNumeroTubosCC < ContadorTubosRamificacion && i < NumberOfPipes) {
124  if(Pipe[i]->getNodoIzq() == FNumeroCC) {
125  FTuboExtremo[FNumeroTubosCC].Pipe = Pipe[i];
126  FTuboExtremo[FNumeroTubosCC].TipoExtremo = nmLeft;
127  FNodoFin[FNumeroTubosCC] = 0;
128  FIndiceCC[FNumeroTubosCC] = 0;
129  FNumeroTubo[FNumeroTubosCC] = Pipe[i]->getNumeroTubo() - 1;
130  FCC[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Beta);
131  FCD[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Landa);
132  FSeccionTubo[FNumeroTubosCC] = __geom::Circle_area(Pipe[i]->GetDiametro(FNodoFin[FNumeroTubosCC]));
133  FNumeroTubosCC++;
134  }
135  if(Pipe[i]->getNodoDer() == FNumeroCC) {
136  FTuboExtremo[FNumeroTubosCC].Pipe = Pipe[i];
137  FTuboExtremo[FNumeroTubosCC].TipoExtremo = nmRight;
138  FNodoFin[FNumeroTubosCC] = Pipe[i]->getNin() - 1;
139  FIndiceCC[FNumeroTubosCC] = 1;
140  FNumeroTubo[FNumeroTubosCC] = Pipe[i]->getNumeroTubo() - 1;
141  FCC[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Landa);
142  FCD[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Beta);
143  FSeccionTubo[FNumeroTubosCC] = __geom::Circle_area(Pipe[i]->GetDiametro(FNodoFin[FNumeroTubosCC]));
144  FNumeroTubosCC++;
145  }
146  i++;
147  }
148 
149  // Inicializacion del transporte de especies quimicas.
150  FFraccionMasicaEspecie = new double[FNumeroEspecies - FIntEGR];
151  FMasaEspecie = new double[FNumeroEspecies - FIntEGR];
152  for(int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
153  FFraccionMasicaEspecie[i] = FTuboExtremo[0].Pipe->GetFraccionMasicaInicial(i);
154  // Se inicializa con el Pipe 0 de modo arbitrario.
155  }
156 
157  } catch(exception & N) {
158  std::cout << "ERROR: TCCRamificacion::AsignaTubos en la condicion de contorno: " << FNumeroCC << std::endl;
159  std::cout << "Tipo de error: " << N.what() << std::endl;
160  throw Exception(N.what());
161  }
162 }
163 
164 // ---------------------------------------------------------------------------
165 // ---------------------------------------------------------------------------
166 
167 void TCCRamificacion::TuboCalculandose(int TuboActual) {
168  try {
169  FTuboActual = TuboActual;
170  if(FTuboActual == 10000) {
171  FTiempoActual = FTuboExtremo[0].Pipe->getTime1();
172  } else {
173  for(int i = 0; i < FNumeroTubosCC; i++) {
174  if(FNumeroTubo[i] == FTuboActual) {
175  FTiempoActual = FTuboExtremo[i].Pipe->getTime1();
176  }
177  }
178  }
179  } catch(exception & N) {
180  std::cout << "ERROR: TCCRamificacion::TuboCalculandose en la condicion de contorno: " << FNumeroCC << std::endl;
181  std::cout << "Tipo de error: " << N.what() << std::endl;
182  throw Exception(N.what());
183  }
184 }
185 
186 // ---------------------------------------------------------------------------
187 // ---------------------------------------------------------------------------
188 
189 void TCCRamificacion::CalculaCondicionContorno(double Time) {
190  try {
191  double sonido_supuesta_ad, sonido_ant_ad, entropia_entrante, corr_entropia;
192  double suma1 = 0., suma2 = 0., sm1 = 0., sm2 = 0., sm3 = 0.;
193  int TuboCalculado = 0;
194  double DeltaT, MasaTotal = 0., g, m, FraccionMasicaAcum = 0.;
195  // Necesarias para el calculo de especies en la BC.
196 
197  sonido_supuesta_ad = 0.;
198 
199  FTiempoActual = Time;
200  DeltaT = FTiempoActual - FTiempoAnterior;
201  FTiempoAnterior = FTiempoActual;
202 
203  if(FTuboActual == 10000) {
204  TuboCalculado = FTuboActual;
205  FGamma = FTuboExtremo[0].Pipe->GetGamma(FNodoFin[0]);
206  } else {
207  for(int i = 0; i < FNumeroTubosCC; i++) {
208  if(FNumeroTubo[i] == FTuboActual) {
209  TuboCalculado = i;
210  }
211  }
212  FGamma = FTuboExtremo[TuboCalculado].Pipe->GetGamma(FNodoFin[TuboCalculado]);
213  }
214  // FGamma=FTuboExtremo[TuboCalculado].Pipe->GetGamma(FNodoFin[TuboCalculado]);
215  FGamma1 = __Gamma::G1(FGamma);
216  FGamma3 = __Gamma::G3(FGamma);
217  FGamma4 = __Gamma::G4(FGamma);
218 
219  for(int i = 0; i < FNumeroTubosCC; i++) {
220  FEntropia[i] = FTuboExtremo[i].Entropia;
221  }
222 
223  do {
224  /* Determinacion de la velocidad del sonido en la ramificacion. */
225  suma1 = 0.;
226  suma2 = 0.;
227  for(int i = 0; i < FNumeroTubosCC; i++) {
228  suma1 = suma1 + (*FCC[i]) * FSeccionTubo[i] / pow2(FEntropia[i]);
229  suma2 = suma2 + FTuboExtremo[i].Entropia * FSeccionTubo[i] / pow2(FEntropia[i]);
230  }
231  sonido_ant_ad = sonido_supuesta_ad;
232  sonido_supuesta_ad = suma1 / suma2; // Velocity del sonido adimensionalizada (si las variables fuesen dimensionales).
233  // Es una especie de promedio respecto de la entropia de cada tubo.
234 
235  /* Determinacion de la velocidad de cada tubo de la ramificacion. Esta
236  velocidad sera positiva si el flujo sale del tubo (tubo saliente) y
237  negativa si el flujo entra en el tubo (tubo entrante). En realidad se trata
238  de la velocidad solo para extremo derecho. Para el extremo izquierdo,esta multiplicada
239  por un signo negativo. */
240  for(int i = 0; i < FNumeroTubosCC; i++) {
241  FVelocity[i] = (*FCC[i] - sonido_supuesta_ad * FTuboExtremo[i].Entropia) / FGamma3;
242  }
243 
244  /* Calculo de la entropia de los tubos entrantes (el flujo entra en ellos).
245  Esta entropia es igual para todos ellos y se calcula como un balance del
246  flujo que llega de los tubos salientes (de velocidad positiva). */
247  sm1 = 0.;
248  sm2 = 0.;
249  sm3 = 0.;
250  for(int i = 0; i < FNumeroTubosCC; i++) {
251  sm3 = sm3 + FTuboExtremo[i].Entropia;
252  if(FVelocity[i] > 2e-6) {
253  sm1 = sm1 + FVelocity[i] * FSeccionTubo[i] * FEntropia[i];
254  sm2 = sm2 + FVelocity[i] * FSeccionTubo[i];
255  }
256  }
257 
258  if(sm2 < 2e-6) {
259  entropia_entrante = sm3 / FNumeroTubosCC;
260  } else {
261  /* Desde el punto de vista teorico esta es la forma correcta. La
262  formula anterior es para evitar errores de indeterminacion si
263  sm2 es muy pequena,por lo que se acepta como aproximacion. */
264  entropia_entrante = sm1 / sm2;
265  }
266  for(int i = 0; i < FNumeroTubosCC; i++) {
267  FEntropia[i] = FTuboExtremo[i].Entropia;
268  if(FVelocity[i] < 0) { // Flujo entrante al tubo
269  FEntropia[i] = entropia_entrante;
270  }
271  }
272  } while((sonido_supuesta_ad - sonido_ant_ad) / (sonido_ant_ad + 0.01) > 1e-4);
273 
274  /* Calculo de las caracteristicas y la entropia en los extremos del tubo que se
275  esta calculando, una vez resuelta la condicion de contorno */
276  if(TuboCalculado != 10000) {
277  corr_entropia = FTuboExtremo[TuboCalculado].Entropia / FEntropia[TuboCalculado];
278  *FCC[TuboCalculado] = (*FCC[TuboCalculado] + FGamma3 * FVelocity[TuboCalculado] * (corr_entropia - 1)) / corr_entropia;
279  *FCD[TuboCalculado] = *FCC[TuboCalculado] - FGamma1 * FVelocity[TuboCalculado];
280  FTuboExtremo[TuboCalculado].Entropia = FEntropia[TuboCalculado];
281 
282  double ason = (*FCC[TuboCalculado] + *FCD[TuboCalculado]) / 2;
283  double Machx = fabs(FVelocity[TuboCalculado]) / ason;
284  if(Machx > 1) {
285  printf("Sonic condition in boundary: %d\n", FNumeroCC);
286  // double Machy = Machx / fabs(Machx) * sqrt
287  // ((pow(Machx, 2) + 2. / FGamma1) / (FGamma4 * pow(Machx, 2) - 1.));
288  // double asonido = (*FCC[TuboCalculado] + *FCD[TuboCalculado]) / 2;
289  // double Sonidoy = asonido * sqrt
290  // ((FGamma3 * pow(Machx, 2) + 1.) / (FGamma3 * pow(Machy, 2) + 1.));
291  //
292  // double Velocidady = Sonidoy * Machy;
293  ReduceSubsonicFlow(ason, FVelocity[TuboCalculado], FGamma);
294  *FCC[TuboCalculado] = ason + FGamma3 * FVelocity[TuboCalculado];
295  *FCD[TuboCalculado] = ason - FGamma3 * FVelocity[TuboCalculado];
296  }
297  } else {
298  for(int i = 0; i < FNumeroTubosCC; i++) {
299  corr_entropia = FTuboExtremo[i].Entropia / FEntropia[i];
300  *FCC[i] = (*FCC[i] + FGamma3 * FVelocity[i] * (corr_entropia - 1)) / corr_entropia;
301  *FCD[i] = *FCC[i] - FGamma1 * FVelocity[i];
302  FTuboExtremo[i].Entropia = FEntropia[i];
303  double Machx = fabs(*FCC[i] - *FCD[i]) / (*FCC[i] + *FCD[i]) * 2 / FGamma1;
304  if(Machx > 1) {
305  printf("Sonic condition in boundary: %d\n", FNumeroCC);
306  double Machy = Machx / fabs(Machx) * sqrt((pow2(Machx) + 2. / FGamma1) / (FGamma4 * pow2(Machx) - 1.));
307  double asonido = (*FCC[i] + *FCD[i]) / 2;
308  double Sonidoy = asonido * sqrt((FGamma3 * pow2(Machx) + 1.) / (FGamma3 * pow2(Machy) + 1.));
309 
310  double Velocidady = Sonidoy * Machy;
311  *FCC[i] = Sonidoy + FGamma3 * Velocidady;
312  *FCD[i] = Sonidoy - FGamma3 * Velocidady;
313  }
314  }
315  }
316 
317  // Transporte de especies quimicas.
318  for(int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
319  FMasaEspecie[j] = 0.;
320  }
321  for(int i = 0; i < FNumeroTubosCC; i++) {
322  if(FVelocity[i] > 0.) { // Flujo Saliente del tubo
323  FDensidad[i] = pow(((*FCC[i] + *FCD[i]) / 2) / FTuboExtremo[i].Entropia, FGamma4);
324  g = FDensidad[i] * FSeccionTubo[i] * FVelocity[i];
325  m = g * DeltaT;
326  MasaTotal += m;
327  for(int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
328  FMasaEspecie[j] += FTuboExtremo[i].Pipe->GetFraccionMasicaCC(FIndiceCC[i], j) * m;
329  }
330  }
331  }
332 
333  if(MasaTotal != 0) {
334  for(int j = 0; j < FNumeroEspecies - 2; j++) {
335  FFraccionMasicaEspecie[j] = FMasaEspecie[j] / MasaTotal;
336  FraccionMasicaAcum += FFraccionMasicaEspecie[j];
337  }
338  FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
339  if(FHayEGR)
340  FFraccionMasicaEspecie[FNumeroEspecies - 1] = FMasaEspecie[FNumeroEspecies - 1] / MasaTotal;
341  }
342 
343  } catch(exception & N) {
344  std::cout << "ERROR: TCCRamificacion::CalculaCondicionContorno en la condicion de contorno: " << FNumeroCC << std::endl;
345  std::cout << "Tipo de error: " << N.what() << std::endl;
346  throw Exception(N.what());
347  }
348 }
349 
350 // ---------------------------------------------------------------------------
351 // ---------------------------------------------------------------------------
352 
353 #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
TTubo::getTime1
double getTime1() const
Gets the time at the following time-step.
Definition: TTubo.cpp:5492
stTuboExtremo
Definition: Globales.h:730
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
TTubo::getNumeroTubo
int getNumeroTubo() const
Gets the pipe id.
Definition: TTubo.cpp:5464