31 #include "TCCPerdidadePresion.h"
42 TCCPerdidadePresion::TCCPerdidadePresion(nmTypeBC TipoCC,
int numCC, nmTipoCalculoEspecies SpeciesModel,
43 int numeroespecies, nmCalculoGamma GammaCalculation,
bool ThereIsEGR) :
44 TCondicionContorno(TipoCC, numCC, SpeciesModel, numeroespecies, GammaCalculation, ThereIsEGR) {
53 if(TipoCC == nmLinearPressureLoss)
55 else if(TipoCC == nmQuadraticPressureLoss)
56 FTipoPP = nmPPCuadratica;
58 printf(
"ERROR en tipo de perdida de presion TCCPerdidadePresion en la condicion de contorno: %d\n", FNumeroCC);
65 TCCPerdidadePresion::~TCCPerdidadePresion() {
66 delete[] FTuboExtremo;
76 if(FNumeroTubo != NULL)
83 void TCCPerdidadePresion::ReadBoundaryData(
const char *FileWAM, fpos_t &filepos,
int NumberOfPipes,
TTubo **Pipe,
84 int nDPF,
TDPF **DPF) {
90 FNodoFin =
new int[2];
91 FIndiceCC =
new int[2];
94 FNumeroTubo =
new int[2];
96 for(
int i = 0; i < 2; i++) {
97 FTuboExtremo[i].Pipe = NULL;
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;
107 FCC[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Beta);
108 FCD[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Landa);
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;
117 FCC[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Landa);
118 FCD[FNumeroTubosCC] = &(FTuboExtremo[FNumeroTubosCC].Beta);
125 FFraccionMasicaEspecie =
new double[FNumeroEspecies - FIntEGR];
126 for(
int i = 0; i < FNumeroEspecies - FIntEGR; i++) {
131 FILE *fich = fopen(FileWAM,
"r");
132 fsetpos(fich, &filepos);
134 fscanf(fich,
"%lf ", &FK);
136 fgetpos(fich, &filepos);
139 }
catch(exception &N) {
140 std::cout <<
"ERROR: TCCPerdidadePresion::LecturaPerdidaPresion en la condicion de contorno: " << FNumeroCC <<
142 std::cout <<
"Tipo de error: " << N.what() << std::endl;
150 void TCCPerdidadePresion::TuboCalculandose(
int TuboActual) {
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;
163 void TCCPerdidadePresion::CalculaCondicionContorno(
double Time) {
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;
170 if(FTuboActual == 10000) {
171 TuboCalculado = FTuboActual;
172 FGamma = FTuboExtremo[0].Pipe->
GetGamma(FNodoFin[0]);
174 for(
int i = 0; i < FNumeroTubosCC; i++) {
175 if(FNumeroTubo[i] == FTuboActual) {
179 FGamma = FTuboExtremo[TuboCalculado].Pipe->
GetGamma(FNodoFin[TuboCalculado]);
181 FGamma1 = __Gamma::G1(FGamma);
182 FGamma3 = __Gamma::G3(FGamma);
183 FGamma2 = __Gamma::G2(FGamma);
184 FGamma5 = __Gamma::G5(FGamma);
186 flujo = (*FCC[1] / FTuboExtremo[1].Entropia) / (*FCC[0] / FTuboExtremo[0].Entropia);
188 if(flujo < .999995) {
190 FRelacionEntropia = FTuboExtremo[0].Entropia / FTuboExtremo[1].Entropia;
193 if((*FCC[0] * 2 / FGamma2) < (*FCC[0] / (FGamma3 / sqrt(FK) + 1) + 1e-6)) {
194 ei = *FCC[0] / (FGamma3 / sqrt(FK) + 1) + 1e-6;
196 ei = *FCC[0] * 2 / FGamma2;
200 if(FTipoPP == nmPPLineal) {
202 ei = QuadraticEqP(FGamma1, 2 * FK, -2 * FK * *FCC[0]) + 1e-6;
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;
212 }
else if(FTipoPP == nmPPCuadratica) {
214 stPerdPresAd PPA(*FCC[0], *FCC[1], FK, FGamma, FRelacionEntropia);
215 vel_sonido_Out =
FindRoot(PPA, ei, ed);
216 vel_sonido_In = PPA.A2;
224 if(PPA.U1 * PPA.U2 < 0) {
230 printf(
"Error en el tipo de perdida de presion en la condicion de contorno: %d\n", FNumeroCC);
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;
237 }
else if(TuboCalculado == 0) {
238 *FCC[0] = vel_sonido_Out + FGamma3 * vel_Out;
239 *FCD[0] = vel_sonido_Out - FGamma3 * vel_Out;
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;
251 for(
int j = 0; j < FNumeroEspecies - 2; j++) {
253 FraccionMasicaAcum += FFraccionMasicaEspecie[j];
255 FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
257 FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboExtremo[0].Pipe->
GetFraccionMasicaCC(FIndiceCC[0],
258 FNumeroEspecies - 1);
259 }
else if(flujo >= 1.000005) {
261 FRelacionEntropia = FTuboExtremo[1].Entropia / FTuboExtremo[0].Entropia;
263 if((*FCC[1] * 2 / FGamma2) < (*FCC[1] / (FGamma3 / sqrt(FK) + 1) + 1e-6)) {
264 ei = *FCC[1] / (FGamma3 / sqrt(FK) + 1) + 1e-6;
266 ei = *FCC[1] * 2 / FGamma2;
270 if(FTipoPP == nmPPLineal) {
272 ei = QuadraticEqP(FGamma1, 2 * FK, -2 * FK * *FCC[1]) + 1e-6;
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;
282 }
else if(FTipoPP == nmPPCuadratica) {
283 stPerdPresAd PPA(*FCC[1], *FCC[0], FK, FGamma, FRelacionEntropia);
284 vel_sonido_Out =
FindRoot(PPA, ei, ed);
285 vel_sonido_In = PPA.A2;
291 printf(
"Error en el tipo de perdida de presion en la condicion de contorno: %d\n", FNumeroCC);
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;
298 }
else if(TuboCalculado == 1) {
299 *FCC[1] = vel_sonido_Out + FGamma3 * vel_Out;
300 *FCD[1] = vel_sonido_Out - FGamma3 * vel_Out;
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;
312 for(
int j = 0; j < FNumeroEspecies - 2; j++) {
314 FraccionMasicaAcum += FFraccionMasicaEspecie[j];
316 FFraccionMasicaEspecie[FNumeroEspecies - 2] = 1. - FraccionMasicaAcum;
318 FFraccionMasicaEspecie[FNumeroEspecies - 1] = FTuboExtremo[1].Pipe->
GetFraccionMasicaCC(FIndiceCC[1],
319 FNumeroEspecies - 1);
322 if(TuboCalculado == 1) {
324 }
else if(TuboCalculado == 0) {
326 }
else if(TuboCalculado == 10000) {
335 catch(exception &N) {
336 std::cout <<
"ERROR: TCCPerdidadePresion::CalculaCondicionContorno en la condicion de contorno: " << FNumeroCC <<
338 std::cout <<
"Tipo de error: " << N.what() << std::endl;
346 #pragma package(smart_init)