31 #include "TCCRamificacion.h"
39 TCCRamificacion::TCCRamificacion(nmTypeBC TipoCC,
int numCC, nmTipoCalculoEspecies SpeciesModel,
int numeroespecies,
40 nmCalculoGamma GammaCalculation,
bool ThereIsEGR) :
41 TCondicionContorno(TipoCC, numCC, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
62 TCCRamificacion::~TCCRamificacion() {
64 if(FTuboExtremo != NULL)
65 delete[] FTuboExtremo;
72 if(FSeccionTubo != NULL)
73 delete[] FSeccionTubo;
78 if(FNumeroTubo != NULL)
86 if(FMasaEspecie != NULL)
87 delete[] FMasaEspecie;
94 void TCCRamificacion::AsignaTubos(
int NumberOfPipes,
TTubo **Pipe) {
97 int ContadorTubosRamificacion = 0;
99 ContadorTubosRamificacion = 0;
101 for(
int i = 0; i < NumberOfPipes; i++) {
102 if(Pipe[i]->getNodoIzq() == FNumeroCC || Pipe[i]->getNodoDer() == FNumeroCC) {
103 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];
118 for(
int i = 0; i < ContadorTubosRamificacion; i++) {
119 FTuboExtremo[i].Pipe = NULL;
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;
130 FCC[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Beta);
131 FCD[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Landa);
132 FSeccionTubo[FNumeroTubosCC] = __geom::Circle_area(Pipe[i]->GetDiametro(FNodoFin[FNumeroTubosCC]));
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;
141 FCC[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Landa);
142 FCD[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Beta);
143 FSeccionTubo[FNumeroTubosCC] = __geom::Circle_area(Pipe[i]->GetDiametro(FNodoFin[FNumeroTubosCC]));
150 FFraccionMasicaEspecie =
new double[FNumeroEspecies - FIntEGR];
151 FMasaEspecie =
new double[FNumeroEspecies - FIntEGR];
152 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
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;
167 void TCCRamificacion::TuboCalculandose(
int TuboActual) {
169 FTuboActual = TuboActual;
170 if(FTuboActual == 10000) {
171 FTiempoActual = FTuboExtremo[0].Pipe->
getTime1();
173 for(
int i = 0; i < FNumeroTubosCC; i++) {
174 if(FNumeroTubo[i] == FTuboActual) {
175 FTiempoActual = FTuboExtremo[i].Pipe->
getTime1();
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;
189 void TCCRamificacion::CalculaCondicionContorno(
double Time) {
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.;
197 sonido_supuesta_ad = 0.;
199 FTiempoActual = Time;
200 DeltaT = FTiempoActual - FTiempoAnterior;
201 FTiempoAnterior = FTiempoActual;
203 if(FTuboActual == 10000) {
204 TuboCalculado = FTuboActual;
205 FGamma = FTuboExtremo[0].Pipe->
GetGamma(FNodoFin[0]);
207 for(
int i = 0; i < FNumeroTubosCC; i++) {
208 if(FNumeroTubo[i] == FTuboActual) {
212 FGamma = FTuboExtremo[TuboCalculado].Pipe->
GetGamma(FNodoFin[TuboCalculado]);
215 FGamma1 = __Gamma::G1(FGamma);
216 FGamma3 = __Gamma::G3(FGamma);
217 FGamma4 = __Gamma::G4(FGamma);
219 for(
int i = 0; i < FNumeroTubosCC; i++) {
220 FEntropia[i] = FTuboExtremo[i].Entropia;
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]);
231 sonido_ant_ad = sonido_supuesta_ad;
232 sonido_supuesta_ad = suma1 / suma2;
240 for(
int i = 0; i < FNumeroTubosCC; i++) {
241 FVelocity[i] = (*FCC[i] - sonido_supuesta_ad * FTuboExtremo[i].Entropia) / FGamma3;
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];
259 entropia_entrante = sm3 / FNumeroTubosCC;
264 entropia_entrante = sm1 / sm2;
266 for(
int i = 0; i < FNumeroTubosCC; i++) {
267 FEntropia[i] = FTuboExtremo[i].Entropia;
268 if(FVelocity[i] < 0) {
269 FEntropia[i] = entropia_entrante;
272 }
while((sonido_supuesta_ad - sonido_ant_ad) / (sonido_ant_ad + 0.01) > 1e-4);
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];
282 double ason = (*FCC[TuboCalculado] + *FCD[TuboCalculado]) / 2;
283 double Machx = fabs(FVelocity[TuboCalculado]) / ason;
285 printf(
"Sonic condition in boundary: %d\n", FNumeroCC);
293 ReduceSubsonicFlow(ason, FVelocity[TuboCalculado], FGamma);
294 *FCC[TuboCalculado] = ason + FGamma3 * FVelocity[TuboCalculado];
295 *FCD[TuboCalculado] = ason - FGamma3 * FVelocity[TuboCalculado];
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;
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.));
310 double Velocidady = Sonidoy * Machy;
311 *FCC[i] = Sonidoy + FGamma3 * Velocidady;
312 *FCD[i] = Sonidoy - FGamma3 * Velocidady;
318 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
319 FMasaEspecie[j] = 0.;
321 for(
int i = 0; i < FNumeroTubosCC; i++) {
322 if(FVelocity[i] > 0.) {
323 FDensidad[i] = pow(((*FCC[i] + *FCD[i]) / 2) / FTuboExtremo[i].Entropia, FGamma4);
324 g = FDensidad[i] * FSeccionTubo[i] * FVelocity[i];
327 for(
int j = 0; j < FNumeroEspecies - FIntEGR; j++) {
334 for(
int j = 0; j < FNumeroEspecies - 2; j++) {
335 FFraccionMasicaEspecie[j] = FMasaEspecie[j] / MasaTotal;
336 FraccionMasicaAcum += FFraccionMasicaEspecie[j];
338 FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
340 FFraccionMasicaEspecie[FNumeroEspecies - 1] = FMasaEspecie[FNumeroEspecies - 1] / MasaTotal;
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;
353 #pragma package(smart_init)