OpenWAM
TMapaComp2Tub.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 #include <iostream>
30 #ifdef __BORLANDC__
31 #include <vcl.h>
32 #endif
33 //#include <cmath>
34 #pragma hdrstop
35 
36 #include "TMapaComp2Tub.h"
37 
38 //---------------------------------------------------------------------------
39 //---------------------------------------------------------------------------
40 
41 TMapaComp2Tub::TMapaComp2Tub(int i) :
42  TCompressorMap() {
43  FNumeroCompresor = i;
44 
45  FGastoRelComp1 = NULL;
46  FGastoBombeo = NULL;
47  FRelCompBombeo = NULL;
48  FRelComp = NULL;
49  FGastoRend = NULL;
50  FRend = NULL;
51  FSpl = NULL;
52  FOrtp = NULL;
53  FCoefSplBombeo = NULL;
54  FCoefSplRC = NULL;
55  FRegimenCurva = NULL;
56  FNumCurvasRen = NULL;
57  FNumCurvasRenAd = NULL;
58  FCoefbSup = NULL;
59  FCoefbInf = NULL;
60  FCoefcSup = NULL;
61  FCoefcInf = NULL;
62  FCoefdSup = NULL;
63  FCoefdInf = NULL;
64  FGastoInt = NULL;
65  FRelCompInt = NULL;
66  FCoefbX = NULL;
67  FCoefcX = NULL;
68  FCoefdX = NULL;
69  FGastoAdim = NULL;
70  FRendAdim = NULL;
71  FCoefSplRend = NULL;
72 // Se supone un rendimiento para la curva de bombeo igual a 0.4
73  FRendCurvaBombeo = 0.4;
74 // Se supone un rendimiento para el punto de massflow maximo de 0.1
75  FRendGastoMaximo = 0.4;
76  Fnegmas = NULL;
77  Fnegrc = NULL;
78  Fsolneg = NULL;
79  FRelCompNuevo = NULL;
80  FGastoCompNuevo = NULL;
81 
82  FCorrect = true;
83 
84 }
85 
86 //---------------------------------------------------------------------------
87 //---------------------------------------------------------------------------
88 
89 TMapaComp2Tub::~TMapaComp2Tub() {
90  if(FGastoRelComp1 != NULL)
91  delete[] FGastoRelComp1;
92  FGastoRelComp1 = NULL;
93  if(FGastoBombeo != NULL)
94  delete[] FGastoBombeo;
95  FGastoBombeo = NULL;
96  if(FRelCompBombeo != NULL)
97  delete[] FRelCompBombeo;
98  FRelCompBombeo = NULL;
99  if(FRelComp != NULL) {
100  for(int i = 0; i < FNumCurvasReg; i++) {
101  delete[] FRelComp[i];
102  }
103  delete[] FRelComp;
104  }
105  FRelComp = NULL;
106  if(FGastoRend != NULL) {
107  for(int i = 0; i < FNumCurvasReg; i++) {
108  delete[] FGastoRend[i];
109  }
110  delete[] FGastoRend;
111  }
112  FGastoRend = NULL;
113  if(FRend != NULL) {
114  for(int i = 0; i < FNumCurvasReg; i++) {
115  delete[] FRend[i];
116  }
117  delete[] FRend;
118  }
119  FRend = NULL;
120  if(FSpl != NULL)
121  delete[] FSpl;
122  FSpl = NULL;
123  if(FOrtp != NULL)
124  delete[] FOrtp;
125  FOrtp = NULL;
126  if(FCoefSplBombeo != NULL)
127  delete[] FCoefSplBombeo;
128  FCoefSplBombeo = NULL;
129  if(FCoefSplRC != NULL)
130  delete[] FCoefSplRC;
131  FCoefSplRC = NULL;
132  if(FRegimenCurva != NULL)
133  delete[] FRegimenCurva;
134  FRegimenCurva = NULL;
135  if(FNumCurvasRen != NULL)
136  delete[] FNumCurvasRen;
137  FNumCurvasRen = NULL;
138  if(FNumCurvasRenAd != NULL)
139  delete[] FNumCurvasRenAd;
140  FNumCurvasRenAd = NULL;
141  if(FCoefbSup != NULL) {
142  for(int i = 0; i < FNumCurvasReg; i++) {
143  delete[] FCoefbSup[i];
144  }
145  delete[] FCoefbSup;
146  }
147  FCoefbSup = NULL;
148  if(FCoefbInf != NULL) {
149  for(int i = 0; i < FNumCurvasReg; i++) {
150  delete[] FCoefbInf[i];
151  }
152  delete[] FCoefbInf;
153  }
154  FCoefbInf = NULL;
155  if(FCoefcSup != NULL) {
156  for(int i = 0; i < FNumCurvasReg; i++) {
157  delete[] FCoefcSup[i];
158  }
159  delete[] FCoefcSup;
160  }
161  FCoefcSup = NULL;
162  if(FCoefcInf != NULL) {
163  for(int i = 0; i < FNumCurvasReg; i++) {
164  delete[] FCoefcInf[i];
165  }
166  delete[] FCoefcInf;
167  }
168  FCoefcInf = NULL;
169  if(FCoefdSup != NULL) {
170  for(int i = 0; i < FNumCurvasReg; i++) {
171  delete[] FCoefdSup[i];
172  }
173  delete[] FCoefdSup;
174  }
175  FCoefdSup = NULL;
176  if(FCoefdInf != NULL) {
177  for(int i = 0; i < FNumCurvasReg; i++) {
178  delete[] FCoefdInf[i];
179  }
180  delete[] FCoefdInf;
181  }
182  FCoefdInf = NULL;
183  delete[] FGastoInt;
184  FGastoInt = NULL;
185  delete[] FRelCompInt;
186  FRelCompInt = NULL;
187  delete[] FCoefbX;
188  FCoefbX = NULL;
189  delete[] FCoefcX;
190  FCoefcX = NULL;
191  delete[] FCoefdX;
192  FCoefdX = NULL;
193  if(FGastoAdim != NULL) {
194  for(int i = 0; i < FNumCurvasReg; i++) {
195  delete[] FGastoAdim[i];
196  }
197  delete[] FGastoAdim;
198  }
199  FGastoAdim = NULL;
200  if(FRendAdim != NULL) {
201  for(int i = 0; i < FNumCurvasReg; i++) {
202  delete[] FRendAdim[i];
203  }
204  delete[] FRendAdim;
205  }
206  FRendAdim = NULL;
207  if(FCoefSplRend != NULL) {
208  for(int i = 0; i < FNumCurvasReg; i++) {
209  delete[] FCoefSplRend[i];
210  }
211  delete[] FCoefSplRend;
212  }
213  FCoefSplRend = NULL;
214  if(Fnegmas != NULL)
215  delete[] Fnegmas;
216  Fnegmas = NULL;
217  if(Fnegrc != NULL)
218  delete[] Fnegrc;
219  Fnegrc = NULL;
220  if(Fsolneg != NULL)
221  delete[] Fsolneg;
222  Fsolneg = NULL;
223  if(FRelCompNuevo != NULL) {
224  for(int i = 0; i < FNumCurvasReg; i++) {
225  delete[] FRelCompNuevo[i];
226  }
227  delete[] FRelCompNuevo;
228  }
229  FRelCompNuevo = NULL;
230  if(FGastoCompNuevo != NULL)
231  delete[] FGastoCompNuevo;
232  FGastoCompNuevo = NULL;
233 }
234 
235 //---------------------------------------------------------------------------
236 // ---------- LeeMapa -------------------------------------------------------
237 // Lectura del mapa del compresor desde el fichero .wam //
238 // Incializacion de la variables del mapa //
239 // Calculo de los coeficientes de la spline que ajusta la linea de bombeo //
240 // Calculo de los coeficientes ortogonales que ajustan el rendimiento //
241 //---------------------------------------------------------------------------
242 
243 void TMapaComp2Tub::LeeMapa(FILE *fich) {
244  double *W = NULL;
245  try {
246 // Datos de referencia
247  std::cout << "LECTURA MAPA COMPRESOR" << std::endl;
248  std::cout << "______________________" << std::endl;
249  fscanf(fich, "%lf %lf ", &FPresionRef, &FTempRef);
250  FTempRef = __units::degCToK(FTempRef);
251  FPresionRef = __units::BarToPa(FPresionRef);
252 
253  std::cout << "Datos de Referencia:" << std::endl;
254  std::cout << "Pressure: " << FPresionRef << " Pa" << std::endl;
255  std::cout << "Temperature: " << FTempRef << " K" << std::endl;
256 
257 // Rango e incremento de regimenes de giro.
258  fscanf(fich, "%lf %lf %lf ", &FRegMin, &FRegMax, &FIncReg);
259  FNumCurvasReg = floor(((FRegMax - FRegMin) / FIncReg) + 0.5) + 1;
260 
261 // Rango e incremento de gastos masicos.
262  fscanf(fich, "%lf %lf %lf ", &FGastoMin, &FGastoMax, &FIncGasto);
263  FGastoMin *= FMassMultiplier;
264  FGastoMax *= FMassMultiplier;
265  FIncGasto *= FMassMultiplier;
266  FNumPuntosGasto = floor(((FGastoMax - FGastoMin) / FIncGasto) + 0.5) + 1;
267 
268 // Numero maximo de curvas de rendimiento
269  fscanf(fich, "%d ", &FNumCurvasRendMax);
270 
271 // Dimensionado Vectores del mapa
272 
273  FGastoRelComp1 = new double[FNumCurvasReg];
274  FGastoBombeo = new double[FNumCurvasReg + 1];
275  FRelCompBombeo = new double[FNumCurvasReg + 1];
276  FRegimenCurva = new double[FNumCurvasReg];
277  FNumCurvasRen = new int[FNumCurvasReg];
278  FNumCurvasRenAd = new int[FNumCurvasReg];
279  FCoefSplBombeo = new double[FNumCurvasReg + 1];
280  FOrtp = new stOrtoPol2Tub[FNumCurvasRendMax + 1];
281  FRelComp = new double*[FNumCurvasReg];
282  for(int i = 0; i < FNumCurvasReg; i++) {
283  FRelComp[i] = new double[FNumPuntosGasto];
284  }
285 
286  FGastoRend = new double*[FNumCurvasReg];
287  FRend = new double*[FNumCurvasReg];
288  for(int i = 0; i < FNumCurvasReg; i++) {
289  FGastoRend[i] = new double[FNumCurvasRendMax];
290  FRend[i] = new double[FNumCurvasRendMax];
291  }
292 
293 // Lectura del mapa
294 
295  for(int i = 0; i < FNumCurvasReg; i++) {
296  FRegimenCurva[i] = FRegMin + FIncReg * (double) i;
297  fscanf(fich, "%lf ", &FGastoRelComp1[i]);
298  FGastoRelComp1[i] *= FMassMultiplier;
299 
300  fscanf(fich, "%lf %lf ", &FGastoBombeo[i], &FRelCompBombeo[i]);
301  FGastoBombeo[i] *= FMassMultiplier;
302  FRelCompBombeo[i] = (FRelCompBombeo[i] - 1) * FCRMultiplier + 1;
303 
304  for(int j = 0; j < FNumPuntosGasto; j++) {
305  fscanf(fich, "%lf ", &FRelComp[i][j]);
306  FRelComp[i][j] = (FRelComp[i][j] - 1) * FCRMultiplier + 1;
307  }
308  FNumCurvasRen[i] = 0;
309  for(int j = 0; j < FNumCurvasRendMax; j++) {
310  fscanf(fich, "%lf %lf ", &FGastoRend[i][j], &FRend[i][j]);
311  FGastoRend[i][j] = FGastoRend[i][j] * FMassMultiplier;
312  FRend[i][j] = FRend[i][j] * FEffMultiplier;
313  if(FRend[i][j] > 0.) {
314  ++FNumCurvasRen[i];
315  }
316  if(FGastoRend[i][j] > FGastoRelComp1[i]) {
317  std::cout << "WARNING: Existen puntos de rendimiento en la curva de regimen: " << FRegimenCurva[i] << std::endl;
318  std::cout << " con valores de massflow mayores del massflow de relacion de compresion 1" << std::endl;
319  std::cout << " -Valor del massflow del punto de rendimiento: " << FGastoRend[i][j] << std::endl;
320  std::cout << " -Valor del massflow de relacion de compresion 1: " << FGastoRelComp1[i] << std::endl;
321  std::cout << " Revisa el mapa. Esto puede inducir errores importantes de interpolacion" << std::endl;
322  std::cout << " en el compresor numero " << FNumeroCompresor << std::endl << std::endl;
323  }
324  }
325  }
326 
327  Cambio_Mapa(FRadioTip, FRadioHub, FRadioRodete);
328 
329 // Obtencion de los coeficientes de la spline que ajusta la curva de bombeo.
330 
331 //Se desplan todos los valores un lugar para poder incluir el centro de coordenadas massflow=0, relacion de compresion=1
332 
333  for(int i = FNumCurvasReg; i > 0; --i) {
334  FGastoBombeo[i] = FGastoBombeo[i - 1];
335  FRelCompBombeo[i] = FRelCompBombeo[i - 1];
336  }
337  FGastoBombeo[0] = 0.;
338  FRelCompBombeo[0] = 1.;
339 
340  Hermite(FNumCurvasReg + 1, FGastoBombeo, FRelCompBombeo, FCoefSplBombeo);
341 
342 // Obtencion de los splines de las curvas de rendimiento frente a massflow adimensionalizado
343  FGastoAdim = new double*[FNumCurvasReg];
344  FRendAdim = new double*[FNumCurvasReg];
345  FCoefSplRend = new double*[FNumCurvasReg];
346  bool RendEnBombeo = false;
347  for(int i = 0; i < FNumCurvasReg; i++) {
348  if(FGastoRend[i][0] < FGastoBombeo[i + 1] * 1.001 && FGastoRend[i][0] > FGastoBombeo[i + 1] * 0.999) {
349  FNumCurvasRenAd[i] = FNumCurvasRen[i] + 1;
350  RendEnBombeo = true;
351  std::cout << "INFO: Surge limit for " << FRegimenCurva[i] << " rpm iso-speed line has an efficiency value deffined" <<
352  std::endl;
353  } else {
354  FNumCurvasRenAd[i] = FNumCurvasRen[i] + 2;
355  RendEnBombeo = false;
356  std::cout << "INFO: Surge limit for " << FRegimenCurva[i] <<
357  " rpm iso-speed line does not have an efficiency value deffined" << std::endl;
358  }
359  FGastoAdim[i] = new double[FNumCurvasRenAd[i]];
360  FRendAdim[i] = new double[FNumCurvasRenAd[i]];
361  FCoefSplRend[i] = new double[FNumCurvasRenAd[i]];
362  //Punto inicial de la curva. Bombeo.
363  FGastoAdim[i][0] = 0.;
364  if(RendEnBombeo) {
365  FRendAdim[i][0] = FRend[i][0];
366  } else {
367  FRendAdim[i][0] = FRendCurvaBombeo;
368  }
369  //Punto final de la curva. Relacion de compresion 1.
370  FGastoAdim[i][FNumCurvasRenAd[i] - 1] = 1.;
371  FRendAdim[i][FNumCurvasRenAd[i] - 1] = FRendGastoMaximo;
372  for(int j = 1; j < FNumCurvasRenAd[i] - 1; j++) {
373  if(RendEnBombeo) {
374  FGastoAdim[i][j] = (FGastoRend[i][j] - FGastoBombeo[i + 1]) / (FGastoRelComp1[i] - FGastoBombeo[i + 1]);
375  FRendAdim[i][j] = FRend[i][j];
376  } else {
377  FGastoAdim[i][j] = (FGastoRend[i][j - 1] - FGastoBombeo[i + 1]) / (FGastoRelComp1[i] - FGastoBombeo[i + 1]);
378  FRendAdim[i][j] = FRend[i][j - 1];
379  }
380  if(FGastoAdim[i][j] <= FGastoAdim[i][j - 1]) {
381  std::cout << "WARNING: La tabla Massflow-Rendimiento en el regimen de " << FRegimenCurva[i] << " rpm esta desordenada"
382  << std::endl;
383  std::cout << " Si se continua con este mapa pueden aparecer errores en la interpolacion" << std::endl;
384  std::cout << " Ordene los datos de forma creciente en massflow" << std::endl;
385  }
386  }
387  Hermite(FNumCurvasRenAd[i], FGastoAdim[i], FRendAdim[i], FCoefSplRend[i]);
388  }
389  ImprimeMapa();
390  /*if(W!=NULL) delete[] W;*/
391  } catch(exception &N) {
392  std::cout << "ERROR: LeeMapa en el compresor: " << FNumeroCompresor << std::endl;
393  std::cout << "Tipo de error: " << N.what() << std::endl;
394  if(W != NULL)
395  delete[] W;
396  throw Exception("ERROR: LeeMapa en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
397  }
398 }
399 
400 //---------------------------------------------------------------------------
401 //--------- Spline ----------------------------------------------------------
402 // Obtiene los diferentes coeficientes de una curva spline que ajusta los //
403 // pares de puntos (x,y) que tienen un numero n de componentes //
404 //---------------------------------------------------------------------------
405 
406 void TMapaComp2Tub::Spline(int n, double *x, double *y, double *sol) {
407  try {
408  double Espaciado = 0.;
409 //Calculo de diferencias
410  for(int i = 1; i < n; ++i) {
411  FSpl[i].h = x[i] - x[i - 1];
412  FSpl[i].dif = y[i] - y[i - 1];
413  if(FSpl[i].h <= 0.) {
414  std::cout << "ERROR: Error al crear la spline" << std::endl;
415  std::cout << " Los valores en X no estan ordenados o existen 2 puntos situados en el mismo X" << std::endl <<
416  std::endl;
417  throw Exception("ERROR: Error al crear la spline");
418  }
419  }
420 
421 //Elementos matriz,coeficiente y terminos independientes
422 //Solo se almacenan elementos diagonales
423  for(int i = 1; i < n - 1; ++i) {
424  Espaciado = FSpl[i + 1].h / FSpl[i].h;
425  if(Espaciado < 0.1 || Espaciado > 10.) {
426  std::cout << "WARNING: Deberias utilizar una distribucion mas uniforme entre los valores de X" << std::endl;
427  std::cout << " utilizados para ajustar la spline y asi evitar problemas" << std::endl << std::endl;
428  }
429  FSpl[i].d = 2 * (FSpl[i + 1].h + FSpl[i].h);
430  FSpl[i].d1 = FSpl[i + 1].h;
431  FSpl[i].b = (FSpl[i + 1].dif / FSpl[i + 1].h - FSpl[i].dif / FSpl[i].h) * 6.;
432  }
433 //Descomposicion de Cholesky
434  FSpl[1].ud = sqrt(FSpl[1].d);
435  for(int i = 2; i < n - 1; i++) {
436  FSpl[i - 1].ud1 = FSpl[i - 1].d1 / FSpl[i - 1].ud;
437  FSpl[i].ud = sqrt(FSpl[i].d - FSpl[i - 1].ud1 * FSpl[i - 1].ud1);
438  }
439 //Sustitucion directa
440  FSpl[1].yp = FSpl[1].b / FSpl[1].ud;
441  for(int i = 2; i < n - 1; ++i) {
442  FSpl[i].yp = (FSpl[i].b - FSpl[i - 1].ud1 * FSpl[i - 1].yp) / FSpl[i].ud;
443  }
444 //Sustitucion inversa
445  sol[0] = 0.;
446  if(n >= 2) {
447  sol[n - 1] = 0;
448  if(n >= 3) {
449  sol[n - 2] = FSpl[n - 2].yp / FSpl[n - 2].ud;
450  //sol[n-2]=0;
451  if(n >= 4) {
452  for(int i = n - 3; i >= 1; --i) {
453  sol[i] = (FSpl[i].yp - FSpl[i].ud1 * sol[i + 1]) / FSpl[i].ud;
454  //sol[i]=0;
455  }
456  }
457  }
458  }
459  /*for(int i=0;i<n;++i){
460  printf("%d %lf %lf %lf\n",i,x[i],y[i],sol[i]);
461  } */
462  } catch(exception &N) {
463  std::cout << "ERROR: Spline en el compresor: " << FNumeroCompresor << std::endl;
464  std::cout << "Tipo de error: " << N.what() << std::endl;
465  throw Exception("ERROR: Spline en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
466  }
467 }
468 
469 //void TMapaComp2Tub::Hermite(int n,double *x,double *y,double *sol)
470 //{
471 //try
472 //{
473 //
474 //double DeltaK,AlphaK,BetaK,TauK;
475 //
476 //for(int i=1;i<n-1;++i){
477 // sol[i]=(y[i]-y[i-1])/2./(x[i]-x[i-1])+(y[i+1]-y[i])/2./(x[i+1]-x[i]);
478 //}
479 //sol[0]=(y[1]-y[0])/(x[1]-x[0]);
480 //sol[n-1]=(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
481 //
482 //for(int i=0;i<n-1;i++){
483 // DeltaK=(y[i+1]-y[i])/(x[i+1]-x[i]);
484 // if(DeltaK==0){
485 // sol[i]=0;
486 // sol[i+1]=0;
487 // }else{
488 // AlphaK=sol[i]/DeltaK;
489 // BetaK=sol[i+1]/DeltaK;
490 // if(BetaK*BetaK+AlphaK*AlphaK>9){
491 // TauK=3/sqrt(BetaK*BetaK+AlphaK*AlphaK);
492 // sol[i]=TauK*AlphaK*DeltaK;
493 // sol[i+1]=TauK*BetaK*DeltaK;
494 // }
495 // }
496 //}
497 //}
498 //catch(Exception &N)
499 //{
500 //std::cout << "ERROR: Hermite en el compresor: " << FNumeroCompresor << std::endl;
501 //std::cout << "Tipo de error: " << N.what() << std::endl;
502 //throw Exception("ERROR: Spline en el compresor: " +std::to_string(FNumeroCompresor)+N.what());
503 //}
504 //}
505 //
506 //void TMapaComp2Tub::Hermite(int n,std::vector<double> x,std::vector<double> y,std::vector<double> *sol)
507 //{
508 //try
509 //{
510 //
511 //double DeltaK,AlphaK,BetaK,TauK;
512 //
513 //for(int i=1;i<n-1;++i){
514 // (*sol)[i]=(y[i]-y[i-1])/2./(x[i]-x[i-1])+(y[i+1]-y[i])/2./(x[i+1]-x[i]);
515 //}
516 //(*sol)[0]=(y[1]-y[0])/(x[1]-x[0]);
517 //(*sol)[n-1]=(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);
518 //
519 //for(int i=0;i<n-1;i++){
520 // DeltaK=(y[i+1]-y[i])/(x[i+1]-x[i]);
521 // if(DeltaK==0){
522 // (*sol)[i]=0;
523 // (*sol)[i+1]=0;
524 // }else{
525 // AlphaK=(*sol)[i]/DeltaK;
526 // BetaK=(*sol)[i+1]/DeltaK;
527 // if(BetaK*BetaK+AlphaK*AlphaK>9){
528 // TauK=3/sqrt(BetaK*BetaK+AlphaK*AlphaK);
529 // (*sol)[i]=TauK*AlphaK*DeltaK;
530 // (*sol)[i+1]=TauK*BetaK*DeltaK;
531 // }
532 // }
533 //}
534 //}
535 //catch(Exception &N)
536 //{
537 //std::cout << "ERROR: Hermite en el compresor: " << FNumeroCompresor << std::endl;
538 //std::cout << "Tipo de error: " << N.what() << std::endl;
539 //throw Exception("ERROR: Spline en el compresor: " +std::to_string(FNumeroCompresor)+N.what());
540 //}
541 //}
542 
543 void TMapaComp2Tub::SplineVector(int n, std::vector<double> x, std::vector<double> y, std::vector<double> sol) {
544  try {
545  double Espaciado = 0.;
546 
547  std::vector<stSpline2Tub> Spl;
548 
549  Spl.resize(n);
550 
551 //Calculo de diferencias
552  for(int i = 1; i < n; ++i) {
553  Spl[i].h = x[i] - x[i - 1];
554  Spl[i].dif = y[i] - y[i - 1];
555  if(Spl[i].h <= 0.) {
556  std::cout << "ERROR: Error al crear la spline" << std::endl;
557  std::cout << " Los valores en X no estan ordenados o existen 2 puntos situados en el mismo X" << std::endl <<
558  std::endl;
559  throw Exception("ERROR: Error al crear la spline");
560  }
561  }
562 
563 //Elementos matriz,coeficiente y terminos independientes
564 //Solo se almacenan elementos diagonales
565  for(int i = 1; i < n - 1; ++i) {
566  Espaciado = Spl[i + 1].h / Spl[i].h;
567  if(Espaciado < 0.1 || Espaciado > 10.) {
568  std::cout << "WARNING: Deberias utilizar una distribucion mas uniforme entre los valores de X" << std::endl;
569  std::cout << " utilizados para ajustar la spline y asi evitar problemas" << std::endl << std::endl;
570  }
571  Spl[i].d = 2 * (Spl[i + 1].h + Spl[i].h);
572  Spl[i].d1 = Spl[i + 1].h;
573  Spl[i].b = (Spl[i + 1].dif / Spl[i + 1].h - Spl[i].dif / Spl[i].h) * 6.;
574  }
575 //Descomposicion de Cholesky
576  Spl[1].ud = sqrt(Spl[1].d);
577  for(int i = 2; i < n - 1; i++) {
578  Spl[i - 1].ud1 = Spl[i - 1].d1 / Spl[i - 1].ud;
579  Spl[i].ud = sqrt(Spl[i].d - Spl[i - 1].ud1 * Spl[i - 1].ud1);
580  }
581 //Sustitucion directa
582  Spl[1].yp = Spl[1].b / Spl[1].ud;
583  for(int i = 2; i < n - 1; ++i) {
584  Spl[i].yp = (Spl[i].b - Spl[i - 1].ud1 * Spl[i - 1].yp) / Spl[i].ud;
585  }
586 //Sustitucion inversa
587  sol[0] = 0.;
588  if(n >= 2) {
589  sol[n - 1] = 0;
590  if(n >= 3) {
591  sol[n - 2] = Spl[n - 2].yp / Spl[n - 2].ud;
592  if(n >= 4) {
593  for(int i = n - 3; i >= 1; --i) {
594  sol[i] = (Spl[i].yp - Spl[i].ud1 * sol[i + 1]) / Spl[i].ud;
595  }
596  }
597  }
598  }
599 
600  } catch(exception &N) {
601  std::cout << "ERROR: Spline en el compresor: " << FNumeroCompresor << std::endl;
602  std::cout << "Tipo de error: " << N.what() << std::endl;
603  throw Exception("ERROR: Spline en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
604  }
605 }
606 
607 //---------------------------------------------------------------------------
608 //------ EvaluaSpline -------------------------------------------------------
609 // A partir de una spline con coeficientes 'sol' que pasa por los n puntos //
610 // (x,y) se evalua su valor para x=punto //
611 //---------------------------------------------------------------------------
612 
613 double TMapaComp2Tub::EvaluaSpline(double punto, int n, double *x, double *y, double *sol) {
614  try {
615  int k = 0;
616 //Determinacion del indice para evaluacion spline
617  while(x[k] < punto && k < n - 1) {
618  ++k;
619  }
620  k = k - 1;
621  if(k < 0)
622  k = 0;
623  if(k >= n - 1)
624  k = n - 2;
625  double h = x[k + 1] - x[k];
626  double h2 = h * h;
627  double dx1 = x[k + 1] - punto;
628  double dx2 = punto - x[k];
629  double dx13 = pow3(dx1);
630  double dx23 = pow3(dx2);
631 
632  double ret_val = (sol[k] * dx13 + sol[k + 1] * dx23 + (6 * y[k + 1] - h2 * sol[k + 1]) * dx2 +
633  (6 * y[k] - h2 * sol[k]) * dx1) / (6 * h);
634 
635  return ret_val;
636 
637  } catch(exception &N) {
638  std::cout << "ERROR: EvaluaSpline en el compresor: " << FNumeroCompresor << std::endl;
639  std::cout << "Tipo de error: " << N.what() << std::endl;
640  throw Exception("ERROR: EvaluaSpline en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
641  }
642 }
643 
644 //double TMapaComp2Tub::EvaluaHermite(double punto,int n,double *x,double *y,double *sol)
645 //{
646 //try
647 //{
648 //double ret_val,h00,h10,h01,h11,t2,t3,t,h;
649 //int k=0;
650 //
651 //if(punto<=x[0]){
652 // ret_val=y[0];
653 //}else if(punto>=x[n-1]){
654 // ret_val=y[n-1];
655 //}else{
656 // while(x[k]<punto && k<n-1){
657 // ++k;
658 // }
659 // h=(x[k]-x[k-1]);
660 // t=(punto-x[k-1])/h;
661 // t2=t*t;
662 // t3=t2*t;
663 // h00=2*t3-3*t2+1;
664 // h10=t3-2*t2+t;
665 // h01=-2*t3+3*t2;
666 // h11=t3-t2;
667 // ret_val=h00*y[k-1]+h*h10*sol[k-1]+h01*y[k]+h*h11*sol[k];
668 //}
669 //
670 //return ret_val;
671 //
672 //}
673 //catch(Exception &N)
674 //{
675 //std::cout << "ERROR: EvaluaHermite en el compresor: " << FNumeroCompresor << std::endl;
676 //std::cout << "Tipo de error: " << N.what() << std::endl;
677 //throw Exception("ERROR: EvaluaHermite en el compresor: " +std::to_string(FNumeroCompresor)+N.what());
678 //}
679 //}
680 //
681 //double TMapaComp2Tub::EvaluaHermite(double punto,int n,std::vector<double> x,std::vector<double> y,std::vector<double> sol)
682 //{
683 //try
684 //{
685 //double ret_val,h00,h10,h01,h11,t2,t3,t,h;
686 //int k=0;
687 //
688 //if(punto<=x[0]){
689 // ret_val=y[0];
690 //}else if(punto>=x[n-1]){
691 // ret_val=y[n-1];
692 //}else{
693 // while(x[k]<punto && k<n-1){
694 // ++k;
695 // }
696 // h=(x[k]-x[k-1]);
697 // t=(punto-x[k-1])/h;
698 // t2=t*t;
699 // t3=t2*t;
700 // h00=2*t3-3*t2+1;
701 // h10=t3-2*t2+t;
702 // h01=-2*t3+3*t2;
703 // h11=t3-t2;
704 // ret_val=h00*y[k-1]+h*h10*sol[k-1]+h01*y[k]+h*h11*sol[k];
705 //}
706 //
707 //return ret_val;
708 //
709 //}
710 //catch(Exception &N)
711 //{
712 //std::cout << "ERROR: EvaluaHermite en el compresor: " << FNumeroCompresor << std::endl;
713 //std::cout << "Tipo de error: " << N.what() << std::endl;
714 //throw Exception("ERROR: EvaluaHermite en el compresor: " +std::to_string(FNumeroCompresor)+N.what());
715 //}
716 //}
717 
718 double TMapaComp2Tub::InterpolaLineal(double punto, int n, double *x, double *y) {
719  try {
720  int k = 0;
721  double delta = 0.;
722  double ret_val = 0.;
723 //Determinacion del indice para evaluacion spline
724 
725  if(punto < x[0]) {
726  ret_val = y[0];
727  } else if(punto > x[n - 1]) {
728  ret_val = y[n - 1];
729  } else {
730  while(x[k] < punto && k < n - 1) {
731  ++k;
732  }
733  delta = (punto - x[k - 1]) / (x[k] - x[k - 1]);
734  ret_val = y[k] * delta + y[k - 1] * (1 - delta);
735 
736  }
737 
738  return ret_val;
739 
740  } catch(exception &N) {
741  std::cout << "ERROR: EvaluaSpline en el compresor: " << FNumeroCompresor << std::endl;
742  std::cout << "Tipo de error: " << N.what() << std::endl;
743  throw Exception("ERROR: EvaluaSpline en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
744  }
745 }
746 
747 double TMapaComp2Tub::EvaluaSpline(double punto, int n, std::vector<double> x, std::vector<double> y,
748  std::vector<double> sol) {
749  try {
750  int k = 0;
751 //Determinacion del indice para evaluacion spline
752  while(x[k] < punto && k < n - 1) {
753  ++k;
754  }
755  k = k - 1;
756  if(k < 0)
757  k = 0;
758  if(k >= n - 1)
759  k = n - 2;
760  double h = x[k + 1] - x[k];
761  double h2 = h * h;
762  double dx1 = x[k + 1] - punto;
763  double dx2 = punto - x[k];
764  double dx13 = pow3(dx1);
765  double dx23 = pow3(dx2);
766 
767  double ret_val = (sol[k] * dx13 + sol[k + 1] * dx23 + (6 * y[k + 1] - h2 * sol[k + 1]) * dx2 +
768  (6 * y[k] - h2 * sol[k]) * dx1) / (6 * h);
769 
770  return ret_val;
771 
772  } catch(exception &N) {
773  std::cout << "ERROR: EvaluaSpline en el compresor: " << FNumeroCompresor << std::endl;
774  std::cout << "Tipo de error: " << N.what() << std::endl;
775  throw Exception("ERROR: EvaluaSpline en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
776  }
777 }
778 
779 //---------------------------------------------------------------------------
780 // ---------- EvaluaRCSplines ----------- //
781 // Evalua la relacion de compresion de la curva interpolada para un //
782 // determinado 'Massflow' //
783 //---------------------------------------------------------------------------
784 
785 double TMapaComp2Tub::EvaluaRCSplines(double Massflow) {
786  try {
787  double ret_val = 0.;
788  if(Massflow < FGastoInt[0]) {
789  ret_val = FRelCompInt[0];
790  } else if(Massflow >= FGastoInt[FNumPuntos]) {
791  ret_val = 1.;
792  } else {
793  int k = 0;
794  int n = FNumPuntos;
795  //Determinacion del indice para evaluacion spline
796  while(FGastoInt[k] < Massflow && k < n) {
797  ++k;
798  }
799  k = k - 1;
800  if(k < 0)
801  k = 0;
802  if(k >= n)
803  k = n - 1;
804  double h = FGastoInt[k + 1] - FGastoInt[k];
805  double h2 = h * h;
806  double dx1 = FGastoInt[k + 1] - Massflow;
807  double dx2 = Massflow - FGastoInt[k];
808  double dx13 = pow3(dx1);
809  double dx23 = pow3(dx2);
810 
811  ret_val = (FCoefSplRC[k] * dx13 + FCoefSplRC[k + 1] * dx23 + (6 * FRelCompInt[k + 1] - h2 * FCoefSplRC[k + 1]) * dx2 +
812  (6 * FRelCompInt[k] - h2 * FCoefSplRC[k]) * dx1) / (6 * h);
813  }
814  return ret_val;
815 
816  } catch(exception &N) {
817  std::cout << "ERROR: EvaluaRCSpline en el compresor: " << FNumeroCompresor << std::endl;
818  std::cout << "Tipo de error: " << N.what() << std::endl;
819  throw Exception("ERROR: EvaluaRCSpline en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
820  }
821 }
822 
823 double TMapaComp2Tub::EvaluaRCHermite(double Massflow) {
824  try {
825  double ret_val = 0., t2 = 0., t3 = 0., t = 0., h00 = 0., h10 = 0., h01 = 0., h11 = 0., h = 0.;
826 
827  if(Massflow <= FGastoInt[0]) {
828  ret_val = FRelCompInt[0];
829  } else if(Massflow >= FGastoInt[FNumPuntos]) {
830  ret_val = 1.;
831  } else {
832  int k = 0;
833  int n = FNumPuntos;
834  //Determinacion del indice para evaluacion spline
835  while(FGastoInt[k] < Massflow && k < n) {
836  ++k;
837  }
838  h = (FGastoInt[k] - FGastoInt[k - 1]);
839  t = (Massflow - FGastoInt[k - 1]) / h;
840  t2 = t * t;
841  t3 = t2 * t;
842  h00 = 2 * t3 - 3 * t2 + 1;
843  h10 = t3 - 2 * t2 + t;
844  h01 = -2 * t3 + 3 * t2;
845  h11 = t3 - t2;
846  ret_val = h00 * FRelCompInt[k - 1] + h * h10 * FCoefSplRC[k - 1] + h01 * FRelCompInt[k] + h * h11 * FCoefSplRC[k];
847  }
848  return ret_val;
849 
850  } catch(exception &N) {
851  std::cout << "ERROR: EvaluaRCSpline en el compresor: " << FNumeroCompresor << std::endl;
852  std::cout << "Tipo de error: " << N.what() << std::endl;
853  throw Exception("ERROR: EvaluaRCSpline en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
854  }
855 }
856 
857 //---------------------------------------------------------------------------
858 //-------- PolOrtogonal -----------------------------------------------------
859 // Obtiene los coeficiente 'b', 'c' y 'd' del polinomio ortogonal utilizado//
860 // para ajustar el rendimiento de compresor donde 'nterms' es el numero de //
861 // terminos que tienen los coeficientes, 'npoint' es el numero de puntos de//
862 // rendimiento que tiene una curva determinada 'ma' representa el massflow //
863 //'rd' el rendimiento para dicho massflow y 'w' es un vector cuyas componentes//
864 // valen 1 Ajuste por minimos cuadrados //
865 //---------------------------------------------------------------------------
866 
867 void TMapaComp2Tub::PolOrtogonal(int nterms, int npoint, double *ma, double *rd, double *w, double *b, double *c,
868  double *d) {
869  double p = 0;
870  try {
871  for(int j = 0; j < nterms; ++j) {
872  b[j] = 0.;
873  d[j] = 0.;
874  FOrtp[j].s = 0.;
875  }
876  c[0] = 0.;
877  for(int i = 0; i < npoint; ++i) {
878  d[0] += rd[i] * w[i];
879  b[0] += ma[i] * w[i];
880  FOrtp[0].s += w[i];
881  }
882  d[0] /= FOrtp[0].s;
883  for(int i = 0; i < npoint; ++i) {
884  FOrtp[i].error = rd[i] - d[0];
885  }
886  b[0] /= FOrtp[0].s;
887  for(int i = 0; i < npoint; ++i) {
888  FOrtp[i].pjm1 = 1.;
889  FOrtp[i].pj = ma[i] - b[0];
890  }
891  int j = 0;
892  while(j < nterms - 1) {
893  ++j;
894  for(int i = 0; i < npoint; ++i) {
895  p = FOrtp[i].pj * w[i];
896  d[j] += FOrtp[i].error * p;
897  p *= FOrtp[i].pj;
898  b[j] += ma[i] * p;
899  FOrtp[j].s += p;
900  }
901  d[j] /= FOrtp[j].s;
902  for(int i = 0; i < npoint; ++i) {
903  FOrtp[i].error -= d[j] * FOrtp[i].pj;
904  }
905  if(j < nterms) {
906  b[j] /= FOrtp[j].s;
907  c[j] = FOrtp[j].s / FOrtp[j - 1].s;
908  for(int i = 0; i < npoint; ++i) {
909  p = FOrtp[i].pj;
910  FOrtp[i].pj = (ma[i] - b[j]) * FOrtp[i].pj - c[j] * FOrtp[i].pjm1;
911  FOrtp[i].pjm1 = p;
912  }
913  }
914  }
915  } catch(exception &N) {
916  std::cout << "ERROR: PolOrtogonal en el compresor: " << FNumeroCompresor << std::endl;
917  std::cout << "Tipo de error: " << N.what() << std::endl;
918  throw Exception("ERROR: PolOrtogonal en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
919  }
920 }
921 
922 //---------------------------------------------------------------------------
923 //--------InterpolaMapa -----------------------------------------------------
924 // Ajusta los coeficientes de una spline que representa la curva de regimen//
925 // 'rtc' a temperatura ambiente 'AmbientTemperature' e interpola los coeficiente del //
926 // polinomio ortogonal que ajusta el rendimiento //
927 //---------------------------------------------------------------------------
928 
929 void TMapaComp2Tub::InterpolaMapa(double rtc, double AmbientTemperature) {
930 
931  double DeltaN = 0.;
932  int Cent = 0;
933  int index = 0;
934 
935  try {
936 
937 // Se adapta el regimen de giro a interpolar al mapa compresor
938  if(FCorrect)
939  FRegComp = rtc * sqrt(FTempRef / AmbientTemperature);
940  else
941  FRegComp = rtc;
942 
943 // Se halla el regimen de giro inferior al que se quiere interpolar
944 
945  FCurvInf = 0;
946  for(int i = 0; i < FNumCurvasReg; ++i) {
947  if(FRegComp >= FRegimenCurva[i]) {
948  FCurvInf = i;
949  }
950  }
951  DeltaN = (FRegComp - FRegimenCurva[FCurvInf]) / FIncReg;
952 
953 // Se interpola linealmente el massflow de bombeo y el de RC = 1
954 
955  if(FRegComp >= FRegMax) {
956  FGastoBombeoX = FGastoBombeo[FCurvInf + 1] + (FGastoBombeo[FCurvInf + 1] - FGastoBombeo[FCurvInf]) / FIncReg *
957  (FRegComp - FRegMax);
958 
959  FRelCompBombeoX = FRelCompBombeo[FCurvInf + 1] + (FRelCompBombeo[FCurvInf + 1] - FRelCompBombeo[FCurvInf]) / FIncReg *
960  (FRegComp - FRegMax);
961 
962  FGastoRelComp1X = FGastoRelComp1[FCurvInf] + (FGastoRelComp1[FCurvInf] - FGastoRelComp1[FCurvInf - 1]) / FIncReg *
963  (FRegComp - FRegMax);
964 
965  } else if(FRegComp < FRegMin) {
966  FGastoBombeoX = FGastoBombeo[FCurvInf + 1] - (FGastoBombeo[FCurvInf + 2] - FGastoBombeo[FCurvInf + 1]) / FIncReg *
967  (FRegMin - FRegComp);
968 
969  FRelCompBombeoX = EvaluaHermite(FGastoBombeoX, FNumCurvasReg + 1, FGastoBombeo, FRelCompBombeo, FCoefSplBombeo);
970 
971  FGastoRelComp1X = FGastoRelComp1[FCurvInf] - (FGastoRelComp1[FCurvInf + 1] - FGastoRelComp1[FCurvInf]) / FIncReg *
972  (FRegMin - FRegComp);
973  } else {
974  FGastoBombeoX = FGastoBombeo[FCurvInf + 1] + (FGastoBombeo[FCurvInf + 2] - FGastoBombeo[FCurvInf + 1]) * DeltaN;
975 
976  FRelCompBombeoX = EvaluaHermite(FGastoBombeoX, FNumCurvasReg + 1, FGastoBombeo, FRelCompBombeo, FCoefSplBombeo);
977 
978  FGastoRelComp1X = FGastoRelComp1[FCurvInf] + (FGastoRelComp1[FCurvInf + 1] - FGastoRelComp1[FCurvInf]) * DeltaN;
979  }
980  if(FGastoBombeoX < 0)
981  FGastoBombeoX = 0.;
982  if(FGastoRelComp1X < FGastoBombeoX)
983  FGastoRelComp1X = FGastoBombeoX + 1e-6;
984 
985 // Calculo de la curva interpolada Relacion de compresion/Masa corregida
986 
987 // Se haya el ultimo massflow del que se tiene informacion para el regimen de giro inferior
988 
989  if(FRegComp < FRegMin) {
990  FNumPuntos = floor((FGastoRelComp1X - FGastoMin) / FIncGasto) + 1;
991  } else {
992  int j = 0;
993  /*while(FRelCompNuevo[FCurvInf][j]>0. && j<FNumPuntosGastoNuevo-1){
994  ++j;
995  }
996  FNumPuntos=j;*/
997  if(FRegComp >= FRegMax) {
998  index = FCurvInf - 1;
999  } else {
1000  index = FCurvInf;
1001  }
1002  while((FGastoInt[j] = FGastoMin + FIncGasto * j) < FGastoRelComp1[index]) {
1003  ++j;
1004  }
1005  FNumPuntos = j;
1006  }
1007 
1008 // Se interpola linealmente la relacion de compresion entre el regimen de giro inferior
1009 // y superior excepto para el ultimo punto.
1010 
1011  if(FRegComp < FRegMin) {
1012  for(int j = 0; j < FNumPuntos; ++j) {
1013  FGastoInt[j] = FGastoMin + FIncGasto * j;
1014  FRelCompInt[j] = FRelCompNuevo[FCurvInf][j] - (FRelCompNuevo[FCurvInf + 1][j] - FRelCompNuevo[FCurvInf][j]) / FIncReg
1015  * (FRegMin - FRegComp);
1016  if(FRelCompInt[j] < 1.)
1017  FRelCompInt[j] = 1.;
1018  }
1019  } else {
1020  for(int j = 0; j < FNumPuntos; ++j) {
1021  //FGastoInt[j]=FGastoMin+FIncGasto*j;
1022  if(FRegComp >= FRegMax) {
1023  FRelCompInt[j] = FRelCompNuevo[FCurvInf][j] + (FRelCompNuevo[FCurvInf][j] - FRelCompNuevo[FCurvInf - 1][j]) * DeltaN;
1024  } else {
1025  FRelCompInt[j] = FRelCompNuevo[FCurvInf][j] + (FRelCompNuevo[FCurvInf + 1][j] - FRelCompNuevo[FCurvInf][j]) * DeltaN;
1026  }
1027  }
1028  }
1029  FGastoInt[FNumPuntos] = FGastoRelComp1X;
1030  FRelCompInt[FNumPuntos] = 1.;
1031 
1032 // Se obtienen los coeficientes que ajusta la spline de la curva interpolada.
1033 
1034  Hermite(FNumPuntos + 1, FGastoInt, FRelCompInt, FCoefSplRC);
1035 
1036 // Se interpolan linealmente los coeficientes de cada polinomio ortogonal para el
1037 // regimen de giro buscado
1038  /*
1039  FNumTerms=3;
1040  Cent=2;
1041  if(FNumCurvasRen[FCurvInf]<FNumTerms){
1042  FNumTerms=FNumCurvasRen[FCurvInf];
1043  --Cent;
1044  }
1045  if(FNumCurvasRen[FCurvInf+1]<FNumTerms){
1046  FNumTerms=FNumCurvasRen[FCurvInf+1];
1047  --Cent;
1048  }
1049  if(Cent>FNumTerms) FNumTerms=3;
1050 
1051  for(int i=0;i<FNumTerms;++i){
1052  if(FRegComp>FRegMax){
1053  FCoefbX[i]=FCoefbSup[FCurvInf+1][i]+(FCoefbSup[FCurvInf+1][i]-FCoefbInf[FCurvInf][i])/FIncReg*(FRegComp-FRegMax);
1054  FCoefcX[i]=FCoefcSup[FCurvInf+1][i]+(FCoefcSup[FCurvInf+1][i]-FCoefcInf[FCurvInf][i])/FIncReg*(FRegComp-FRegMax);
1055  FCoefdX[i]=FCoefdSup[FCurvInf+1][i]+(FCoefdSup[FCurvInf+1][i]-FCoefdInf[FCurvInf][i])/FIncReg*(FRegComp-FRegMax);
1056  }else if(FRegComp<FRegMin){
1057  FCoefbX[i]=FCoefbInf[FCurvInf][i]-(FCoefbSup[FCurvInf+1][i]-FCoefbInf[FCurvInf][i])/FIncReg*(FRegMin-FRegComp);
1058  FCoefcX[i]=FCoefbInf[FCurvInf][i]-(FCoefcSup[FCurvInf+1][i]-FCoefcInf[FCurvInf][i])/FIncReg*(FRegMin-FRegComp);
1059  FCoefdX[i]=FCoefbInf[FCurvInf][i]-(FCoefdSup[FCurvInf+1][i]-FCoefdInf[FCurvInf][i])/FIncReg*(FRegMin-FRegComp);
1060  }else{
1061  FCoefbX[i]=FCoefbInf[FCurvInf][i]+(FCoefbSup[FCurvInf+1][i]-FCoefbInf[FCurvInf][i])*DeltaN;
1062  FCoefcX[i]=FCoefbInf[FCurvInf][i]+(FCoefcSup[FCurvInf+1][i]-FCoefcInf[FCurvInf][i])*DeltaN;
1063  FCoefdX[i]=FCoefbInf[FCurvInf][i]+(FCoefdSup[FCurvInf+1][i]-FCoefdInf[FCurvInf][i])*DeltaN;
1064  }
1065  } */
1066  } catch(exception &N) {
1067  std::cout << "ERROR: InterpolaMapa en el compresor: " << FNumeroCompresor << std::endl;
1068  std::cout << "Tipo de error: " << N.what() << std::endl;
1069  throw Exception("ERROR: InterpolaMapa en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
1070  }
1071 }
1072 
1073 //---------------------------------------------------------------------------
1074 // ---------EvaluaRendimiento -----------------------------------------------
1075 // Evalua el rendimiento a partir de los coeficientes del polinomio //
1076 // ortogonal de la curva de funcionamiento para un valor de //
1077 // masa de aire 'MasaAire' //
1078 //---------------------------------------------------------------------------
1079 
1080 double TMapaComp2Tub::EvaluaRendimiento(double MasaAire) {
1081  int k = 0;
1082  double prev2 = 0.;
1083  double prev = 0.;
1084  double ret_val = 0.;
1085  try {
1086  k = FNumTerms - 1;
1087  ret_val = FCoefdX[k];
1088  prev = 0.;
1089  --k;
1090  while(k >= 0) {
1091  prev2 = prev;
1092  prev = ret_val;
1093  ret_val = FCoefdX[k] + (MasaAire - FCoefbX[k]) * prev - FCoefcX[k + 1] * prev2;
1094  --k;
1095  }
1096  return ret_val;
1097  } catch(exception &N) {
1098  std::cout << "ERROR: EvaluaRendimiento en el compresor: " << FNumeroCompresor << std::endl;
1099  std::cout << "Tipo de error: " << N.what() << std::endl;
1100  throw Exception("ERROR: EvaluaRendimiento en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
1101  }
1102 }
1103 
1104 //---------------------------------------------------------------------------
1105 //---------------------------------------------------------------------------
1106 
1107 double TMapaComp2Tub::EvaluaRendSplines(double MasaAire) {
1108  double AireAdim = 0.;
1109  double ret_val = 0.;
1110  double RendInf = 0.;
1111  double RendSup = 0.;
1112  double DeltaN = 0.;
1113  int index = 0;
1114  try {
1115  if(FRegComp >= FRegMax) {
1116  index = FCurvInf - 1;
1117  } else {
1118  index = FCurvInf;
1119  }
1120 
1121  if(MasaAire < FGastoBombeoX) {
1122  //En Bombeo se toma como rendimiento el del ultimo punto.
1123  DeltaN = (FRegComp - FRegimenCurva[FCurvInf]) / FIncReg;
1124 
1125  RendInf = FRendAdim[index][0];
1126  RendSup = FRendAdim[index + 1][0];
1127 
1128  if(FRegComp > FRegMax) {
1129  ret_val = RendSup + (RendSup - RendInf) / FIncReg * (FRegComp - FRegMax);
1130  } else if(FRegComp < FRegMin) {
1131  ret_val = RendInf - (RendSup - RendInf) / FIncReg * (FRegMin - FRegComp);
1132  } else {
1133  ret_val = RendInf + (RendSup - RendInf) * DeltaN;
1134  }
1135  } else if(MasaAire > FGastoRelComp1X) {
1136  ret_val = FRendGastoMaximo;
1137  } else {
1138  AireAdim = (MasaAire - FGastoBombeoX) / (FGastoRelComp1X - FGastoBombeoX);
1139 
1140  RendInf = EvaluaHermite(AireAdim, FNumCurvasRenAd[index], FGastoAdim[index], FRendAdim[FCurvInf],
1141  FCoefSplRend[FCurvInf]);
1142 
1143  //RendInf=InterpolaLineal(AireAdim,FNumCurvasRenAd[index],FGastoAdim[index],
1144  // FRendAdim[FCurvInf]);
1145 
1146  RendSup = EvaluaHermite(AireAdim, FNumCurvasRenAd[index + 1], FGastoAdim[index + 1], FRendAdim[index + 1],
1147  FCoefSplRend[index + 1]);
1148 
1149  //RendSup=InterpolaLineal(AireAdim,FNumCurvasRenAd[index+1],FGastoAdim[index+1],
1150  // FRendAdim[index+1]);
1151 
1152  DeltaN = (FRegComp - FRegimenCurva[FCurvInf]) / FIncReg;
1153 
1154  if(FRegComp >= FRegMax) {
1155  ret_val = RendSup + (RendSup - RendInf) / FIncReg * (FRegComp - FRegMax);
1156  } else if(FRegComp < FRegMin) {
1157  ret_val = RendInf - (RendSup - RendInf) / FIncReg * (FRegMin - FRegComp);
1158  } else {
1159  ret_val = RendInf + (RendSup - RendInf) * DeltaN;
1160  }
1161  }
1162  return ret_val;
1163  } catch(exception &N) {
1164  std::cout << "ERROR: EvaluaRendSplines en el compresor: " << FNumeroCompresor << std::endl;
1165  std::cout << "Tipo de error: " << N.what() << std::endl;
1166  throw Exception("ERROR: EvaluaRendSplines en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
1167  }
1168 }
1169 
1170 //---------------------------------------------------------------------------
1171 //---------------------------------------------------------------------------
1172 
1173 double TMapaComp2Tub::GetRelCompInt(int i) {
1174  try {
1175  return FRelCompInt[i];
1176  } catch(exception &N) {
1177  std::cout << "ERROR: GetRelCompInt en el compresor: " << FNumeroCompresor << std::endl;
1178  std::cout << "Tipo de error: " << N.what() << std::endl;
1179  throw Exception("ERROR: GetRelCompInt en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
1180  }
1181 }
1182 
1183 //---------------------------------------------------------------------------
1184 //---------------------------------------------------------------------------
1185 
1186 double TMapaComp2Tub::GetGastoInt(int i) {
1187  try {
1188  return FGastoInt[i];
1189  } catch(exception &N) {
1190  std::cout << "ERROR: GetGastoInt en el compresor: " << FNumeroCompresor << std::endl;
1191  std::cout << "Tipo de error: " << N.what() << std::endl;
1192  throw Exception("ERROR: GetGastoInt en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
1193  }
1194 }
1195 
1196 //---------------------------------------------------------------------------
1197 //---------------------------------------------------------------------------
1198 
1199 double TMapaComp2Tub::BuscaRegimen(double RC, double Massflow, double AmbientTemperature) {
1200  try {
1201  double RCSup = 0.;
1202  double RCInf = 0.;
1203  double ret_val = 0., val1 = 0., val2 = 0.;
1204  double reg = 0.;
1205  int i = 0, sup, inf;
1206  while(RC > RCSup && i < FNumCurvasReg) {
1207  reg = FRegimenCurva[i] * AmbientTemperature / FTempRef;
1208  InterpolaMapa(reg, AmbientTemperature);
1209 
1210  if(Massflow < FGastoRelComp1[i]) {
1211  RCSup = EvaluaRCHermite(Massflow);
1212  } else {
1213  RCSup = 1.;
1214  }
1215  sup = i;
1216  ++i;
1217  }
1218  if(RC > RCSup) {
1219  std::cout << "WARNING: El Punto se sale por arriba del mapa " << std::endl;
1220  //sup=i-1;
1221  }
1222  if(sup > 0) {
1223  inf = sup - 1;
1224  reg = FRegimenCurva[inf] * AmbientTemperature / FTempRef;
1225  InterpolaMapa(reg, AmbientTemperature);
1226 
1227  if(Massflow < FGastoRelComp1[inf]) {
1228  RCInf = EvaluaRCHermite(Massflow);
1229  } else {
1230  std::cout << "WARNING: El siguiente punto puede estar mal calculado:" << std::endl;
1231  std::cout << " Relacionde compresion: " << RC << std::endl;
1232  std::cout << " Massflow masico: " << Massflow << std::endl;
1233  std::cout << " Deberias aumentar el numero de curvas del mapa" << std::endl << std::endl;
1234  RCInf = 1.;
1235  }
1236  ret_val = (RC - RCInf) / (RCSup - RCInf) * (FRegimenCurva[sup] - FRegimenCurva[inf]) + FRegimenCurva[inf];
1237  } else {
1238  inf = sup;
1239  sup = inf + 1;
1240  RCInf = RCSup;
1241  InterpolaMapa(FRegimenCurva[sup], __units::degCToK(AmbientTemperature));
1242  RCSup = EvaluaRCHermite(Massflow);
1243  val1 = (RC - RCInf) / (RCSup - RCInf) * (FRegimenCurva[sup] - FRegimenCurva[inf]) + FRegimenCurva[inf];
1244 
1245  val2 = (RC - 1) / (RCInf - 1) * (FRegimenCurva[inf]);
1246  ret_val = (val1 + val2) / 2;
1247  std::cout << "WARNING: El siguiente punto puede estar mal calculado:" << std::endl;
1248  std::cout << " Relacionde compresion: " << RC << std::endl;
1249  std::cout << " Massflow masico: " << Massflow << std::endl;
1250  std::cout << " Deberias incluir curvas con menor regimen de giro" << std::endl << std::endl;
1251  }
1252  return ret_val;
1253  } catch(exception &N) {
1254  std::cout << "ERROR: BuscaRegimen en el compresor: " << FNumeroCompresor << std::endl;
1255  std::cout << "Tipo de error: " << N.what() << std::endl;
1256  throw Exception("ERROR: BuscaRegimen en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
1257  }
1258 }
1259 
1260 //---------------------------------------------------------------------------
1261 //---------------------------------------------------------------------------
1262 
1263 void TMapaComp2Tub::Cambio_Mapa(double radtip, double radhub, double radrodete) {
1264  try {
1265 //FILE *fplot;
1266  int cont = 0;
1267  double r1 = 0.;
1268  int n = 0;
1269 
1270 //los valores de ra,rb y r2, deberan ser medidos directamente en el compresor y ser leidos desde lee_compresor
1271 
1272  Fnegmas = new double[5];
1273  Fnegrc = new double[5];
1274  Fsolneg = new double[5];
1275  for(int i = 0; i < 5; i++) {
1276  Fsolneg[i] = 0.;
1277  }
1278  n = 0.;
1279  while((double) n * FIncGasto < 0.07) {
1280  ++n;
1281  }
1282  FGastoMin = -n * FIncGasto;
1283 
1284  r1 = sqrt((pow2(radhub) + pow2(radtip)) / 2);
1285 //FGastoMin=-0.07;
1286  FNumPuntosGastoNuevo = floor(((FGastoMax - FGastoMin) / FIncGasto) + 0.5) + 1;
1287 
1288  FSpl = new stSpline2Tub[FNumPuntosGastoNuevo + 1];
1289  FGastoInt = new double[FNumPuntosGastoNuevo + 1];
1290  FRelCompInt = new double[FNumPuntosGastoNuevo + 1];
1291  FCoefSplRC = new double[FNumPuntosGastoNuevo + 1];
1292 
1293  FRelCompNuevo = new double*[FNumCurvasReg];
1294  for(int i = 0; i < FNumCurvasReg; i++) {
1295  FRelCompNuevo[i] = new double[FNumPuntosGastoNuevo];
1296  }
1297 
1298  FGastoCompNuevo = new double[FNumPuntosGastoNuevo];
1299  for(int j = 0; j < FNumPuntosGastoNuevo; ++j) {
1300  FGastoCompNuevo[j] = FGastoMin + j * FIncGasto;
1301  }
1302 
1303  for(int i = 0; i < FNumCurvasReg; ++i) {
1304  Fnegmas[0] = FGastoMin;
1305  Fnegmas[1] = -0.0001;
1306  Fnegmas[2] = 0;
1307  Fnegmas[3] = FGastoBombeo[i] - 0.0001;
1308  Fnegmas[4] = FGastoBombeo[i];
1309  Fnegrc[0] = FRelCompBombeo[i];
1310  Fnegrc[1] = pow(1 + 1 / (__Gamma::Cp_x2 * (FTempRef)) * pow2(__units::RPMToRad_s(FRegMin + i * FIncReg)) * (pow2(
1311  radrodete) - pow2(r1)), __Gamma::G_9);
1312  Fnegrc[2] = pow(1 + 1 / (__Gamma::Cp_x2 * (FTempRef)) * pow2(__units::RPMToRad_s(FRegMin + i * FIncReg)) * (pow2(
1313  radrodete) - pow2(r1)), __Gamma::G_9);
1314  Fnegrc[3] = FRelCompBombeo[i];
1315  Fnegrc[4] = FRelCompBombeo[i];
1316  Hermite(5, Fnegmas, Fnegrc, Fsolneg);
1317  cont = 1;
1318  for(int j = 0; j < FNumPuntosGastoNuevo; ++j) {
1319  if(FGastoMin + j * FIncGasto < FGastoBombeo[i]) {
1320  FRelCompNuevo[i][j] = EvaluaHermite(FGastoMin + j * FIncGasto, 5, Fnegmas, Fnegrc, Fsolneg);
1321  if(FGastoMin + j * FIncGasto > FIncGasto / 10)
1322  ++cont;
1323  } else {
1324  FRelCompNuevo[i][j] = FRelComp[i][cont];
1325  ++cont;
1326  }
1327  //printf("%d %lf %lf\n",i,FGastoMin+j*FIncGasto,FRelCompNuevo[i][j]);
1328  }
1329  }
1330 
1331  /*fplot=fopen ("datos.m","w");
1332  fprintf (fplot, "function [a, NumPuntosGasto, IncGasto, GastoMin, GastoRelComp1]=datos \na=[");
1333  for (int i=0; i<FNumCurvasReg; ++i) {
1334  for(int j=0; j<FNumPuntosGastoNuevo;++j){
1335  if(i!=FNumCurvasReg || j!=FNumPunyosGastoNuevo) fprintf (fplot,"%lf;\n ",Faux[i][j]);
1336  else fprintf (fplot,"%lf ];\n",rcom[i]);
1337 
1338  }
1339  fprintf (fplot,"NumPuntosGasto= %d; \nIncGasto= %lf; \nGastoMin= %lf; \nmarc1= [",FNumPuntosGastoNuevo,FIncGasto,FGastoMin);
1340  for (i=0; i<FNumCurvasReg; i++) fprintf (fplot,"%lf;\n",FGastoRelComp1[i]);
1341  fprintf (fplot,"];");
1342  fclose (fplot); */
1343  } catch(exception &N) {
1344  std::cout << "ERROR: Cambio_Mapa: " << FNumeroCompresor << std::endl;
1345  std::cout << "Tipo de error: " << N.what() << std::endl;
1346  throw Exception(N.what());
1347  }
1348 }
1349 
1350 //---------------------------------------------------------------------------
1351 //---------------------------------------------------------------------------
1352 
1353 void TMapaComp2Tub::ImprimeMapa() {
1354  std::cout << "Printing compressor map .";
1355  FILE *fich;
1356  FILE *fichrd;
1357  fich = fopen("MapaInterp.txt", "w");
1358  fichrd = fopen("MapaRend.txt", "w");
1359  double inc = 0.;
1360  double *massflow;
1361  double **rc;
1362  double **rend;
1363  double *rcbomb;
1364  double *gasbomb;
1365 
1366  massflow = new double[200];
1367 
1368  rcbomb = new double[FNumCurvasReg + 3];
1369 
1370  gasbomb = new double[FNumCurvasReg + 3];
1371 
1372  rc = new double*[200];
1373  rend = new double*[200];
1374  for(int j = 0; j < 200; j++) {
1375  rc[j] = new double[FNumCurvasReg + 2];
1376  rend[j] = new double[FNumCurvasReg + 2];
1377  }
1378  inc = (FGastoMax - FGastoMin) / 199;
1379  for(int j = 0; j < 200; j++) {
1380  massflow[j] = FGastoMin + inc * (double) j;
1381  }
1382  gasbomb[0] = 0.;
1383  rcbomb[0] = 1.;
1384 
1385  for(int i = 0; i < FNumCurvasReg + 1; i++) {
1386  std::cout << ".";
1387  InterpolaMapa(FRegMin + FIncReg * (double) i, FTempRef);
1388  for(int j = 0; j < 200; j++) {
1389  if(massflow[j] > FGastoRelComp1X) {
1390  rc[j][i] = 1.;
1391  rend[j][i] = 0.4;
1392  } else {
1393  rc[j][i] = EvaluaRCHermite(massflow[j]);
1394  if(massflow[j] <= FGastoBombeoX) {
1395  rend[j][i] = EvaluaRendSplines(FGastoBombeoX);
1396  } else {
1397  rend[j][i] = EvaluaRendSplines(massflow[j]);
1398  }
1399  }
1400  }
1401  gasbomb[i + 1] = FGastoBombeoX;
1402  rcbomb[i + 1] = FRelCompBombeoX;
1403 
1404  }
1405 
1406  for(int j = 0; j < 200; j++) {
1407  if(j < FNumCurvasReg + 1 + 2) {
1408  fprintf(fich, "%lf\t%lf\t", gasbomb[j], rcbomb[j]);
1409  } else {
1410  fprintf(fich, "\t\t");
1411  }
1412  fprintf(fich, "%lf\t", massflow[j]);
1413  fprintf(fichrd, "%lf\t", massflow[j]);
1414  for(int i = 0; i < FNumCurvasReg + 1; i++) {
1415  fprintf(fich, "%lf\t", rc[j][i]);
1416  fprintf(fichrd, "%lf\t", rend[j][i]);
1417 
1418  }
1419  fprintf(fich, "\n");
1420  fprintf(fichrd, "\n");
1421  }
1422 
1423  for(int j = 0; j < 200; j++) {
1424  delete[] rc[j];
1425  delete[] rend[j];
1426  }
1427  delete[] massflow;
1428  delete[] rcbomb;
1429  delete[] gasbomb;
1430  delete[] rend;
1431  delete[] rc;
1432 
1433  fclose(fich);
1434  fclose(fichrd);
1435 
1436  FILE *fich2;
1437  fich2 = fopen("MapaReal.txt", "w");
1438 
1439  for(int j = 0; j < FNumPuntosGastoNuevo; ++j) {
1440  fprintf(fich2, "%lf\t", FGastoMin + j * FIncGasto);
1441  for(int i = 0; i < FNumCurvasReg; ++i) {
1442  fprintf(fich2, "%lf\t", FRelCompNuevo[i][j]);
1443  }
1444  fprintf(fich2, "\n");
1445  }
1446 
1447  fclose(fich2);
1448 
1449  std::cout << "Finished." << std::endl;
1450 }
1451 
1452 //---------------------------------------------------------------------------
1453 //---------------------------------------------------------------------------
1454 
1455 void TMapaComp2Tub::ReadGTPowerMap(FILE *fich, int correct) {
1456  try {
1457  double speed = 0., mass = 0., pres = 0., eff = 0.;
1458  int i = 0; //Curva de isoregimen
1459  int j = 0; //Puntos de la curva
1460  if(correct == 1)
1461  FCorrect = true;
1462  else
1463  FCorrect = false;
1464  fscanf(fich, "%lf %lf", &FPresionRef, &FTempRef);
1465  FPresionRef = __units::BarToPa(FPresionRef);
1466  FGTSpeed.resize(i + 1);
1467  FGTMass.resize(i + 1);
1468  FGTPres.resize(i + 1);
1469  FGTEff.resize(i + 1);
1470  while(!feof(fich)) {
1471  fscanf(fich, "%lf %lf %lf %lf", &speed, &mass, &pres, &eff);
1472  if(!feof(fich)) {
1473  if(j > 0) {
1474  if(speed != FGTSpeed[i][j - 1]) {
1475  i++;
1476  j = 0;
1477  FGTSpeed.resize(i + 1);
1478  FGTMass.resize(i + 1);
1479  FGTPres.resize(i + 1);
1480  FGTEff.resize(i + 1);
1481  }
1482  }
1483  FGTSpeed[i].push_back(speed);
1484  FGTMass[i].push_back(mass);
1485  FGTPres[i].push_back(pres);
1486  FGTEff[i].push_back(eff);
1487  j++;
1488  }
1489  }
1490 
1491  } catch(exception &N) {
1492  std::cout << "ERROR: Cambio_Mapa: " << FNumeroCompresor << std::endl;
1493  std::cout << "Tipo de error: " << N.what() << std::endl;
1494  throw Exception(N.what());
1495  }
1496 }
1497 
1498 //---------------------------------------------------------------------------
1499 //---------------------------------------------------------------------------
1500 
1501 void TMapaComp2Tub::RearrangeGTPowerMap(double rtip, double rhub, double rwheel) {
1502  try {
1503  double FlowMass = 0.;
1504 
1505  //FTempRef=303;
1506  //FPresionRef=0.96;
1507  //for(int i=0;i<FGTSpeed.size();i++){
1508  // for(int j=0;j<FGTSpeed[i].size();j++){
1509  // printf("%lf %lf %lf %lf\n",FGTSpeed[i][j],FGTMass[i][j],FGTPres[i][j],FGTEff[i][j]);
1510  // }
1511  //}
1512 
1513  // Regimen maximo y minimo del mapa.
1514  FRegMin = FGTSpeed[0][0];
1515  FRegMax = FGTSpeed[FGTSpeed.size() - 1][0];
1516  FIncReg = (FRegMax - FRegMin) / (FGTSpeed.size() - 1);
1517 
1518  FNumCurvasReg = FGTSpeed.size();
1519  // Massflow maximo y minimo del mapa.
1520  FGastoMin = 0;
1521  FGastoMax = FGTMass[FGTMass.size() - 1][FGTMass[FGTMass.size() - 1].size() - 1];
1522  FIncGasto = (FGastoMax - FGastoMin) / 20;
1523 
1524  FNumPuntosGasto = 21;
1525 
1526  FGTCoefCR.resize(FNumCurvasReg);
1527 
1528  FNumCurvasRendMax = 0;
1529 
1530  for(int i = 0; i < FNumCurvasReg; ++i) {
1531  FGTCoefCR[i].resize(FGTSpeed[i].size());
1532  //int k=FGTSpeed[i].size();
1533  Hermite(FGTSpeed[i].size(), FGTMass[i], FGTPres[i], &FGTCoefCR[i]);
1534  /*for(int j=0;j<FGTSpeed[i].size();j++){
1535  printf("%lf\n",FGTCoefCR[i][j]);
1536  }
1537  printf("\n"); */
1538  //SplineVector(k,FGTMass[i],FGTPres[i],FGTMass[i]);
1539  if(FGTEff[i].size() > (dMatrix::size_type) FNumCurvasRendMax) {
1540  FNumCurvasRendMax = FGTEff[i].size();
1541  }
1542  }
1543 
1544  FGastoRelComp1 = new double[FNumCurvasReg];
1545  FGastoBombeo = new double[FNumCurvasReg + 1];
1546  FRelCompBombeo = new double[FNumCurvasReg + 1];
1547  FRegimenCurva = new double[FNumCurvasReg];
1548  FNumCurvasRen = new int[FNumCurvasReg];
1549  FNumCurvasRenAd = new int[FNumCurvasReg];
1550  FCoefSplBombeo = new double[FNumCurvasReg + 1];
1551  FOrtp = new stOrtoPol2Tub[FNumCurvasRendMax + 1];
1552  FRelComp = new double*[FNumCurvasReg];
1553  for(int i = 0; i < FNumCurvasReg; i++) {
1554  FRelComp[i] = new double[FNumPuntosGasto];
1555  }
1556  FGastoRend = new double*[FNumCurvasReg];
1557  FRend = new double*[FNumCurvasReg];
1558  for(int i = 0; i < FNumCurvasReg; i++) {
1559  FGastoRend[i] = new double[FNumCurvasRendMax];
1560  FRend[i] = new double[FNumCurvasRendMax];
1561  }
1562 
1563  for(int i = 0; i < FNumCurvasReg; i++) {
1564  FRegimenCurva[i] = FRegMin + FIncReg * (double) i;
1565  FGastoRelComp1[i] = FGTMass[i][FGTMass[i].size() - 1];
1566  FGastoBombeo[i] = FGTMass[i][0];
1567  FRelCompBombeo[i] = FGTPres[i][0];
1568  for(int j = 0; j < 21; j++) {
1569  FlowMass = FGastoMin + (double) j * FIncGasto;
1570  if(FlowMass > FGastoBombeo[i]) {
1571  if(FlowMass < FGastoRelComp1[i]) {
1572  FRelComp[i][j] = EvaluaHermite(FlowMass, FGTMass[i].size(), FGTMass[i], FGTPres[i], FGTCoefCR[i]);
1573  } else {
1574  FRelComp[i][j] = 1.;
1575  }
1576 
1577  } else {
1578  FRelComp[i][j] = FRelCompBombeo[i];
1579  }
1580  }
1581  FNumCurvasRen[i] = FGTEff[i].size() - 1;
1582  for(Uint j = 0; j < FGTEff[i].size() - 1; j++) {
1583  FGastoRend[i][j] = FGTMass[i][j];
1584  FRend[i][j] = FGTEff[i][j];
1585  }
1586  }
1587 
1588  Cambio_Mapa(rtip, rhub, rwheel);
1589 
1590  // Interpolacion curva de bombeo.
1591 
1592  for(int i = FNumCurvasReg; i > 0; --i) {
1593  FGastoBombeo[i] = FGastoBombeo[i - 1];
1594  FRelCompBombeo[i] = FRelCompBombeo[i - 1];
1595  }
1596  FGastoBombeo[0] = 0.;
1597  FRelCompBombeo[0] = 1.;
1598 
1599  Hermite(FNumCurvasReg + 1, FGastoBombeo, FRelCompBombeo, FCoefSplBombeo);
1600 
1601  // Obtencion de los splines de las curvas de rendimiento frente a massflow adimensionalizado
1602  FGastoAdim = new double*[FNumCurvasReg];
1603  FRendAdim = new double*[FNumCurvasReg];
1604  FCoefSplRend = new double*[FNumCurvasReg];
1605  bool RendEnBombeo = false;
1606  for(int i = 0; i < FNumCurvasReg; i++) {
1607  if(FGastoRend[i][0] < FGastoBombeo[i + 1] * 1.001 && FGastoRend[i][0] > FGastoBombeo[i + 1] * 0.999) {
1608  FNumCurvasRenAd[i] = FNumCurvasRen[i] + 1;
1609  RendEnBombeo = true;
1610  std::cout << "INFO: El punto de bombeo de la curva de " << FRegimenCurva[i] << " rpm tiene rendimiento" << std::endl;
1611  } else {
1612  FNumCurvasRenAd[i] = FNumCurvasRen[i] + 2;
1613  RendEnBombeo = false;
1614  std::cout << "INFO: El punto de bombeo de la curva de " << FRegimenCurva[i] << " rpm no tiene rendimiento" << std::endl;
1615  }
1616  FGastoAdim[i] = new double[FNumCurvasRenAd[i]];
1617  FRendAdim[i] = new double[FNumCurvasRenAd[i]];
1618  FCoefSplRend[i] = new double[FNumCurvasRenAd[i]];
1619  //Punto inicial de la curva. Bombeo.
1620  FGastoAdim[i][0] = 0.;
1621  if(RendEnBombeo) {
1622  FRendAdim[i][0] = FRend[i][0];
1623  } else {
1624  FRendAdim[i][0] = FRendCurvaBombeo;
1625  }
1626  //Punto final de la curva. Relacion de compresion 1.
1627  FGastoAdim[i][FNumCurvasRenAd[i] - 1] = 1.;
1628  FRendAdim[i][FNumCurvasRenAd[i] - 1] = FRendGastoMaximo;
1629  for(int j = 1; j < FNumCurvasRenAd[i] - 1; j++) {
1630  if(RendEnBombeo) {
1631  FGastoAdim[i][j] = (FGastoRend[i][j] - FGastoBombeo[i + 1]) / (FGastoRelComp1[i] - FGastoBombeo[i + 1]);
1632  FRendAdim[i][j] = FRend[i][j];
1633  } else {
1634  FGastoAdim[i][j] = (FGastoRend[i][j - 1] - FGastoBombeo[i + 1]) / (FGastoRelComp1[i] - FGastoBombeo[i + 1]);
1635  FRendAdim[i][j] = FRend[i][j - 1];
1636  }
1637  if(FGastoAdim[i][j] <= FGastoAdim[i][j - 1]) {
1638  std::cout << "WARNING: La tabla Massflow-Rendimiento en el regimen de " << FRegimenCurva[i] << " rpm esta desordenada"
1639  << std::endl;
1640  std::cout << " Si se continua con este mapa pueden aparecer errores en la interpolacion" << std::endl;
1641  std::cout << " Ordene los datos de forma creciente en massflow" << std::endl;
1642  }
1643  }
1644  Hermite(FNumCurvasRenAd[i], FGastoAdim[i], FRendAdim[i], FCoefSplRend[i]);
1645  }
1646 
1647  ImprimeMapa();
1648 
1649  } catch(exception &N) {
1650  std::cout << "ERROR: RearrangeGTPowerMap: " << FNumeroCompresor << std::endl;
1651  std::cout << "Tipo de error: " << N.what() << std::endl;
1652  throw Exception(N.what());
1653  }
1654 }
1655 
1656 void TMapaComp2Tub::WriteMapForWAM() {
1657  try {
1658  double gast = 0.;
1659 
1660  FILE *fich2;
1661  fich2 = fopen("MapaWAM.cmp", "w");
1662 
1663  fprintf(fich2, "%lf %lf ", FPresionRef, FTempRef);
1664  fprintf(fich2, "%lf %lf %lf ", FRegMin, FRegMax, FIncReg);
1665  fprintf(fich2, "%lf %lf %lf ", FGastoMin, FGastoMax, FIncGasto);
1666  fprintf(fich2, "%d", FNumCurvasRendMax);
1667  fprintf(fich2, "\n");
1668  for(int i = 0; i < FNumCurvasReg; i++) {
1669  fprintf(fich2, "%lf\n", FGastoRelComp1[i]);
1670  fprintf(fich2, "%lf %lf\n", FGastoBombeo[i + 1], FRelCompBombeo[i + 1]);
1671  for(int j = 0; j < FNumPuntosGasto; j++) {
1672  gast = (double) j * FIncGasto;
1673  if(gast > FGastoRelComp1[i]) {
1674  fprintf(fich2, "0.0\n");
1675  } else {
1676  fprintf(fich2, "%lf\n", FRelComp[i][j]);
1677  }
1678  }
1679  for(int j = 0; j < FNumCurvasRendMax; j++) {
1680  if(j < FNumCurvasRen[i]) {
1681  fprintf(fich2, "%lf %lf\n", FGastoRend[i][j], FRend[i][j]);
1682  } else {
1683  fprintf(fich2, "0.0 0.0\n");
1684  }
1685  }
1686 
1687  }
1688 
1689  fclose(fich2);
1690 
1691  } catch(exception &N) {
1692  std::cout << "ERROR: WriteMapForWAM: " << FNumeroCompresor << std::endl;
1693  std::cout << "Tipo de error: " << N.what() << std::endl;
1694  throw Exception(N.what());
1695  }
1696 }
1697 
1698 #pragma package(smart_init)
stOrtoPol2Tub
Definition: Globales.h:949
Uint
unsigned int Uint
Unsigned integer.
Definition: Math_wam.h:69
pow3
T pow3(T x)
Returns x to the power of 3.
Definition: Math_wam.h:101
stSpline2Tub
Definition: Globales.h:935
TCompressorMap
Definition: TCompressorMap.h:12
Exception
Custom exception class.
Definition: Exception.hpp:39
pow2
T pow2(T x)
Returns x to the power of 2.
Definition: Math_wam.h:88