OpenWAM
TMapaComp.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 "TMapaComp.h"
37 
38 //---------------------------------------------------------------------------
39 //---------------------------------------------------------------------------
40 
41 TMapaComp::TMapaComp(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  FRegimenCurva = NULL;
55  FNumCurvasRen = NULL;
56  FNumCurvasRenAd = NULL;
57  FCoefbSup = NULL;
58  FCoefbInf = NULL;
59  FCoefcSup = NULL;
60  FCoefcInf = NULL;
61  FCoefdSup = NULL;
62  FCoefdInf = NULL;
63  FGastoInt = NULL;
64  FRelCompInt = NULL;
65  FCoefSplRC = 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.1;
76 
77 }
78 
79 //---------------------------------------------------------------------------
80 //---------------------------------------------------------------------------
81 
82 TMapaComp::~TMapaComp() {
83  if(FGastoRelComp1 != NULL)
84  delete[] FGastoRelComp1;
85  if(FGastoBombeo != NULL)
86  delete[] FGastoBombeo;
87  if(FRelCompBombeo != NULL)
88  delete[] FRelCompBombeo;
89  if(FRelComp != NULL) {
90  for(int i = 0; i < FNumCurvasReg; i++) {
91  delete[] FRelComp[i];
92  }
93  delete[] FRelComp;
94  }
95  if(FGastoRend != NULL) {
96  for(int i = 0; i < FNumCurvasReg; i++) {
97  delete[] FGastoRend[i];
98  }
99  delete[] FGastoRend;
100  }
101  if(FRend != NULL) {
102  for(int i = 0; i < FNumCurvasReg; i++) {
103  delete[] FRend[i];
104  }
105  delete[] FRend;
106  }
107  if(FSpl != NULL)
108  delete[] FSpl;
109  if(FOrtp != NULL)
110  delete[] FOrtp;
111  if(FCoefSplBombeo != NULL)
112  delete[] FCoefSplBombeo;
113  if(FRegimenCurva != NULL)
114  delete[] FRegimenCurva;
115  if(FNumCurvasRen != NULL)
116  delete[] FNumCurvasRen;
117  if(FNumCurvasRenAd != NULL)
118  delete[] FNumCurvasRenAd;
119  if(FCoefbSup != NULL) {
120  for(int i = 0; i < FNumCurvasReg; i++) {
121  delete[] FCoefbSup[i];
122  }
123  delete[] FCoefbSup;
124  }
125  if(FCoefbInf != NULL) {
126  for(int i = 0; i < FNumCurvasReg; i++) {
127  delete[] FCoefbInf[i];
128  }
129  delete[] FCoefbInf;
130  }
131  if(FCoefcSup != NULL) {
132  for(int i = 0; i < FNumCurvasReg; i++) {
133  delete[] FCoefcSup[i];
134  }
135  delete[] FCoefcSup;
136  }
137  if(FCoefcInf != NULL) {
138  for(int i = 0; i < FNumCurvasReg; i++) {
139  delete[] FCoefcInf[i];
140  }
141  delete[] FCoefcInf;
142  }
143  if(FCoefdSup != NULL) {
144  for(int i = 0; i < FNumCurvasReg; i++) {
145  delete[] FCoefdSup[i];
146  }
147  delete[] FCoefdSup;
148  }
149  if(FCoefdInf != NULL) {
150  for(int i = 0; i < FNumCurvasReg; i++) {
151  delete[] FCoefdInf[i];
152  }
153  delete[] FCoefdInf;
154  }
155  delete[] FGastoInt;
156  delete[] FRelCompInt;
157  delete[] FCoefSplRC;
158  delete[] FCoefbX;
159  delete[] FCoefcX;
160  delete[] FCoefdX;
161  if(FGastoAdim != NULL) {
162  for(int i = 0; i < FNumCurvasReg; i++) {
163  delete[] FGastoAdim[i];
164  }
165  delete[] FGastoAdim;
166  }
167  if(FRendAdim != NULL) {
168  for(int i = 0; i < FNumCurvasReg; i++) {
169  delete[] FRendAdim[i];
170  }
171  delete[] FRendAdim;
172  }
173  if(FCoefSplRend != NULL) {
174  for(int i = 0; i < FNumCurvasReg; i++) {
175  delete[] FCoefSplRend[i];
176  }
177  delete[] FCoefSplRend;
178  }
179 }
180 
181 //---------------------------------------------------------------------------
182 // ---------- LeeMapa -------------------------------------------------------
183 // Lectura del mapa del compresor desde el fichero .wam //
184 // Incializacion de la variables del mapa //
185 // Calculo de los coeficientes de la spline que ajusta la linea de bombeo //
186 // Calculo de los coeficientes ortogonales que ajustan el rendimiento //
187 //---------------------------------------------------------------------------
188 
189 void TMapaComp::LeeMapa(FILE *fich) {
190  double *W = NULL;
191  try {
192 // Datos de referencia
193  std::cout << "LECTURA MAPA COMPRESOR" << std::endl;
194  std::cout << "______________________" << std::endl;
195  fscanf(fich, "%lf %lf ", &FPresionRef, &FTempRef);
196  FTempRef = __units::degCToK(FTempRef);
197  ;
198  FPresionRef = __units::BarToPa(FPresionRef);
199  fscanf(fich, "%lf %lf %lf ", &FMassMultiplier, &FCRMultiplier, &FEffMultiplier);
200 
201  std::cout << "Datos de Referencia:" << std::endl;
202  std::cout << "Pressure: " << FPresionRef << " Pa" << std::endl;
203  std::cout << "Temperature: " << FTempRef << " K" << std::endl;
204 
205 // Rango e incremento de regimenes de giro.
206  fscanf(fich, "%lf %lf %lf ", &FRegMin, &FRegMax, &FIncReg);
207  FNumCurvasReg = floor(((FRegMax - FRegMin) / FIncReg) + 0.5) + 1;
208 
209 // Rango e incremento de gastos masicos.
210  fscanf(fich, "%lf %lf %lf ", &FGastoMin, &FGastoMax, &FIncGasto);
211  FGastoMin *= FMassMultiplier;
212  FGastoMax *= FMassMultiplier;
213  FIncGasto *= FMassMultiplier;
214  FNumPuntosGasto = floor(((FGastoMax - FGastoMin) / FIncGasto) + 0.5) + 1;
215 
216 // Numero maximo de curvas de rendimiento
217  fscanf(fich, "%d ", &FNumCurvasRendMax);
218 
219 // Dimensionado Vectores del mapa
220 
221  FGastoRelComp1 = new double[FNumCurvasReg];
222  FGastoBombeo = new double[FNumCurvasReg + 1];
223  FRelCompBombeo = new double[FNumCurvasReg + 1];
224  FRegimenCurva = new double[FNumCurvasReg];
225  FNumCurvasRen = new int[FNumCurvasReg];
226  FNumCurvasRenAd = new int[FNumCurvasReg];
227  FCoefSplBombeo = new double[FNumCurvasReg + 1];
228  FSpl = new stSpline[FNumPuntosGasto + 1];
229  FOrtp = new stOrtoPol[FNumCurvasRendMax + 1];
230  FRelComp = new double*[FNumCurvasReg];
231  for(int i = 0; i < FNumCurvasReg; i++) {
232  FRelComp[i] = new double[FNumPuntosGasto];
233  }
234  FGastoInt = new double[FNumPuntosGasto + 1];
235  FRelCompInt = new double[FNumPuntosGasto + 1];
236  FCoefSplRC = new double[FNumPuntosGasto + 1];
237 
238  FGastoRend = new double*[FNumCurvasReg];
239  FRend = new double*[FNumCurvasReg];
240  for(int i = 0; i < FNumCurvasReg; i++) {
241  FGastoRend[i] = new double[FNumCurvasRendMax];
242  FRend[i] = new double[FNumCurvasRendMax];
243  }
244 
245 // Lectura del mapa
246 
247  for(int i = 0; i < FNumCurvasReg; i++) {
248  FRegimenCurva[i] = FRegMin + FIncReg * (double) i;
249  fscanf(fich, "%lf ", &FGastoRelComp1[i]);
250  FGastoRelComp1[i] *= FMassMultiplier;
251 
252  fscanf(fich, "%lf %lf ", &FGastoBombeo[i], &FRelCompBombeo[i]);
253  FGastoBombeo[i] = FGastoBombeo[i] * FMassMultiplier;
254  FRelCompBombeo[i] = (FRelCompBombeo[i] - 1) * FCRMultiplier + 1;
255 
256  for(int j = 0; j < FNumPuntosGasto; j++) {
257  fscanf(fich, "%lf ", &FRelComp[i][j]);
258  FRelComp[i][j] = (FRelComp[i][j] - 1) * FCRMultiplier + 1;
259  }
260  FNumCurvasRen[i] = 0;
261  for(int j = 0; j < FNumCurvasRendMax; j++) {
262  fscanf(fich, "%lf %lf ", &FGastoRend[i][j], &FRend[i][j]);
263  FGastoRend[i][j] = FGastoRend[i][j] * FMassMultiplier;
264  FRend[i][j] = FRend[i][j] * FEffMultiplier;
265  if(FRend[i][j] > 0.) {
266  ++FNumCurvasRen[i];
267  }
268  if(FGastoRend[i][j] > FGastoRelComp1[i]) {
269  std::cout << "WARNING: Existen puntos de rendimiento en la curva de regimen: " << FRegimenCurva[i] << std::endl;
270  std::cout << " con valores de massflow mayores del massflow de relacion de compresion 1" << std::endl;
271  std::cout << " -Valor del massflow del punto de rendimiento: " << FGastoRend[i][j] << std::endl;
272  std::cout << " -Valor del massflow de relacion de compresion 1: " << FGastoRelComp1[i] << std::endl;
273  std::cout << " Revisa el mapa. Esto puede inducir errores importantes de interpolacion" << std::endl;
274  std::cout << " en el compresor numero " << FNumeroCompresor << std::endl << std::endl;
275  }
276  }
277  }
278 // Obtencion de los coeficientes de la spline que ajusta la curva de bombeo.
279 
280 //Se desplan todos los valores un lugar para poder incluir el centro de coordenadas massflow=0, relacion de compresion=1
281 
282  for(int i = FNumCurvasReg; i > 0; --i) {
283  FGastoBombeo[i] = FGastoBombeo[i - 1];
284  FRelCompBombeo[i] = FRelCompBombeo[i - 1];
285  }
286  FGastoBombeo[0] = 0.;
287  FRelCompBombeo[0] = 1.;
288 
289  Hermite(FNumCurvasReg + 1, FGastoBombeo, FRelCompBombeo, FCoefSplBombeo);
290 
291 // Obtencion de los coeficientes polinomios ortogonales.
292  /*double NumTerms = 0.;
293 
294  W=new double[FNumCurvasRendMax];
295 
296  FCoefbInf=new double*[FNumCurvasReg];
297  FCoefbSup=new double*[FNumCurvasReg];
298  FCoefcInf=new double*[FNumCurvasReg];
299  FCoefcSup=new double*[FNumCurvasReg];
300  FCoefdInf=new double*[FNumCurvasReg];
301  FCoefdSup=new double*[FNumCurvasReg];
302 
303  for(int i=0;i<FNumCurvasReg;i++){
304  FCoefbInf[i]=new double[3];
305  FCoefbSup[i]=new double[3];
306  FCoefcInf[i]=new double[3];
307  FCoefcSup[i]=new double[3];
308  FCoefdInf[i]=new double[3];
309  FCoefdSup[i]=new double[3];
310  }
311 
312  FCoefbX=new double[3];
313  FCoefcX=new double[3];
314  FCoefdX=new double[3];
315 
316  for(int i=0;i<FNumCurvasReg-1;i++){
317  NumTerms=3;
318  if(FNumCurvasRen[i]<NumTerms){
319  NumTerms=FNumCurvasRen[i];
320  }
321  for(int j=0;j<FNumCurvasRen[i];++j){
322  W[j]=1.;
323  }
324 
325  PolOrtogonal(NumTerms,FNumCurvasRen[i],FGastoRend[i],FRend[i],W,FCoefbInf[i],FCoefcInf[i],
326  FCoefdInf[i]);
327 
328  if(FNumCurvasRen[i+1]<NumTerms){
329  NumTerms=FNumCurvasRen[i+1];
330  }
331  for(int j=0;j<FNumCurvasRen[i+1];++j){
332  W[j]=1.;
333  }
334 
335  PolOrtogonal(NumTerms,FNumCurvasRen[i+1],FGastoRend[i+1],FRend[i+1],W,FCoefbSup[i+1],FCoefcSup[i+1],
336  FCoefdSup[i+1]);
337  } */
338 // Obtencion de los splines de las curvas de rendimiento frente a massflow adimensionalizado
339  FGastoAdim = new double*[FNumCurvasReg];
340  FRendAdim = new double*[FNumCurvasReg];
341  FCoefSplRend = new double*[FNumCurvasReg];
342  bool RendEnBombeo = false;
343  for(int i = 0; i < FNumCurvasReg; i++) {
344  if(FGastoRend[i][0] < FGastoBombeo[i + 1] * 1.001 && FGastoRend[i][0] > FGastoBombeo[i + 1] * 0.999) {
345  FNumCurvasRenAd[i] = FNumCurvasRen[i] + 1;
346  RendEnBombeo = true;
347  std::cout << "INFO: El punto de bombeo de la curva de " << FRegimenCurva[i] << " rpm tiene rendimiento" << std::endl;
348  } else {
349  FNumCurvasRenAd[i] = FNumCurvasRen[i] + 2;
350  RendEnBombeo = false;
351  std::cout << "INFO: El punto de bombeo de la curva de " << FRegimenCurva[i] << " rpm no tiene rendimiento" << std::endl;
352  }
353  FGastoAdim[i] = new double[FNumCurvasRenAd[i]];
354  FRendAdim[i] = new double[FNumCurvasRenAd[i]];
355  FCoefSplRend[i] = new double[FNumCurvasRenAd[i]];
356  //Punto inicial de la curva. Bombeo.
357  FGastoAdim[i][0] = 0.;
358  if(RendEnBombeo) {
359  FRendAdim[i][0] = FRend[i][0];
360  } else {
361  FRendAdim[i][0] = FRendCurvaBombeo;
362  }
363  //Punto final de la curva. Relacion de compresion 1.
364  FGastoAdim[i][FNumCurvasRenAd[i] - 1] = 1.;
365  FRendAdim[i][FNumCurvasRenAd[i] - 1] = FRendGastoMaximo;
366  for(int j = 1; j < FNumCurvasRenAd[i] - 1; j++) {
367  if(RendEnBombeo) {
368  FGastoAdim[i][j] = (FGastoRend[i][j] - FGastoBombeo[i + 1]) / (FGastoRelComp1[i] - FGastoBombeo[i + 1]);
369  FRendAdim[i][j] = FRend[i][j];
370  } else {
371  FGastoAdim[i][j] = (FGastoRend[i][j - 1] - FGastoBombeo[i + 1]) / (FGastoRelComp1[i] - FGastoBombeo[i + 1]);
372  FRendAdim[i][j] = FRend[i][j - 1];
373  }
374  if(FGastoAdim[i][j] <= FGastoAdim[i][j - 1]) {
375  std::cout << "WARNING: La tabla Massflow-Rendimiento en el regimen de " << FRegimenCurva[i] << " rpm esta desordenada"
376  << std::endl;
377  std::cout << " Si se continua con este mapa pueden aparecer errores en la interpolacion" << std::endl;
378  std::cout << " Ordene los datos de forma creciente en massflow" << std::endl;
379  }
380  }
381  Hermite(FNumCurvasRenAd[i], FGastoAdim[i], FRendAdim[i], FCoefSplRend[i]);
382  }
383 //ImprimeMapa();
384  /*if(W!=NULL) delete[] W;*/
385  } catch(exception &N) {
386  std::cout << "ERROR: LeeMapa en el compresor: " << FNumeroCompresor << std::endl;
387  std::cout << "Tipo de error: " << N.what() << std::endl;
388  if(W != NULL)
389  delete[] W;
390  throw Exception("ERROR: LeeMapa en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
391  }
392 }
393 
394 //---------------------------------------------------------------------------
395 //--------- Spline ----------------------------------------------------------
396 // Obtiene los diferentes coeficientes de una curva spline que ajusta los //
397 // pares de puntos (x,y) que tienen un numero n de componentes //
398 //---------------------------------------------------------------------------
399 
400 void TMapaComp::Spline(int n, double *x, double *y, double *sol) {
401  try {
402  double Espaciado = 0.;
403 //Calculo de diferencias
404  for(int i = 1; i < n; ++i) {
405  FSpl[i].h = x[i] - x[i - 1];
406  FSpl[i].dif = y[i] - y[i - 1];
407  if(FSpl[i].h <= 0.) {
408  std::cout << "ERROR: Error al crear la spline" << std::endl;
409  std::cout << " Los valores en X no estan ordenados o existen 2 puntos situados en el mismo X" << std::endl <<
410  std::endl;
411  throw Exception("ERROR: Error al crear la spline");
412  }
413  }
414 
415 //Elementos matriz,coeficiente y terminos independientes
416 //Solo se almacenan elementos diagonales
417  for(int i = 1; i < n - 1; ++i) {
418  Espaciado = FSpl[i + 1].h / FSpl[i].h;
419  if(Espaciado < 0.1 || Espaciado > 10.) {
420  std::cout << "WARNING: Deberias utilizar una distribucion mas uniforme entre los valores de X" << std::endl;
421  std::cout << " utilizados para ajustar la spline y asi evitar problemas" << std::endl << std::endl;
422  }
423  FSpl[i].d = 2 * (FSpl[i + 1].h + FSpl[i].h);
424  FSpl[i].d1 = FSpl[i + 1].h;
425  FSpl[i].b = (FSpl[i + 1].dif / FSpl[i + 1].h - FSpl[i].dif / FSpl[i].h) * 6.;
426  }
427 //Descomposicion de Cholesky
428  FSpl[1].ud = sqrt(FSpl[1].d);
429  for(int i = 2; i < n - 1; i++) {
430  FSpl[i - 1].ud1 = FSpl[i - 1].d1 / FSpl[i - 1].ud;
431  FSpl[i].ud = sqrt(FSpl[i].d - FSpl[i - 1].ud1 * FSpl[i - 1].ud1);
432  }
433 //Sustitucion directa
434  FSpl[1].yp = FSpl[1].b / FSpl[1].ud;
435  for(int i = 2; i < n - 1; ++i) {
436  FSpl[i].yp = (FSpl[i].b - FSpl[i - 1].ud1 * FSpl[i - 1].yp) / FSpl[i].ud;
437  }
438 //Sustitucion inversa
439  sol[0] = 0.;
440  if(n >= 2) {
441  sol[n - 1] = 0;
442  if(n >= 3) {
443  sol[n - 2] = FSpl[n - 2].yp / FSpl[n - 2].ud;
444  if(n >= 4) {
445  for(int i = n - 3; i >= 1; --i) {
446  sol[i] = (FSpl[i].yp - FSpl[i].ud1 * sol[i + 1]) / FSpl[i].ud;
447  }
448  }
449  }
450  }
451 
452  } catch(exception &N) {
453  std::cout << "ERROR: Spline en el compresor: " << FNumeroCompresor << std::endl;
454  std::cout << "Tipo de error: " << N.what() << std::endl;
455  throw Exception("ERROR: Spline en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
456  }
457 }
458 
459 //---------------------------------------------------------------------------
460 //------ EvaluaSpline -------------------------------------------------------
461 // A partir de una spline con coeficientes 'sol' que pasa por los n puntos //
462 // (x,y) se evalua su valor para x=punto //
463 //---------------------------------------------------------------------------
464 
465 double TMapaComp::EvaluaSpline(double punto, int n, double *x, double *y, double *sol) {
466  try {
467  int k = 0;
468 //Determinacion del indice para evaluacion spline
469  while(x[k] < punto && k < n - 1) {
470  ++k;
471  }
472  k = k - 1;
473  if(k < 0)
474  k = 0;
475  if(k >= n - 1)
476  k = n - 2;
477  double h = x[k + 1] - x[k];
478  double h2 = h * h;
479  double dx1 = x[k + 1] - punto;
480  double dx2 = punto - x[k];
481  double dx13 = pow3(dx1);
482  double dx23 = pow3(dx2);
483 
484  double ret_val = (sol[k] * dx13 + sol[k + 1] * dx23 + (6 * y[k + 1] - h2 * sol[k + 1]) * dx2 +
485  (6 * y[k] - h2 * sol[k]) * dx1) / (6 * h);
486 
487  return ret_val;
488 
489  } catch(exception &N) {
490  std::cout << "ERROR: EvaluaSpline en el compresor: " << FNumeroCompresor << std::endl;
491  std::cout << "Tipo de error: " << N.what() << std::endl;
492  throw Exception("ERROR: EvaluaSpline en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
493  }
494 }
495 
496 //---------------------------------------------------------------------------
497 // ---------- EvaluaRCSplines ----------- //
498 // Evalua la relacion de compresion de la curva interpolada para un //
499 // determinado 'Massflow' //
500 //---------------------------------------------------------------------------
501 
502 double TMapaComp::EvaluaRCSplines(double Massflow) {
503  try {
504  double ret_val = 0.;
505  if(Massflow < FGastoInt[0]) {
506  ret_val = FRelCompInt[0];
507  } else {
508  int k = 0;
509  int n = FNumPuntos;
510  //Determinacion del indice para evaluacion spline
511  while(FGastoInt[k] < Massflow && k < n) {
512  ++k;
513  }
514  k = k - 1;
515  if(k < 0)
516  k = 0;
517  if(k >= n)
518  k = n - 1;
519  double h = FGastoInt[k + 1] - FGastoInt[k];
520  double h2 = h * h;
521  double dx1 = FGastoInt[k + 1] - Massflow;
522  double dx2 = Massflow - FGastoInt[k];
523  double dx13 = pow3(dx1);
524  double dx23 = pow3(dx2);
525 
526  ret_val = (FCoefSplRC[k] * dx13 + FCoefSplRC[k + 1] * dx23 + (6 * FRelCompInt[k + 1] - h2 * FCoefSplRC[k + 1]) * dx2 +
527  (6 * FRelCompInt[k] - h2 * FCoefSplRC[k]) * dx1) / (6 * h);
528 
529  if(ret_val < 1.)
530  ret_val = 1.;
531  }
532 
533  return ret_val;
534 
535  } catch(exception &N) {
536  std::cout << "ERROR: EvaluaRCSpline en el compresor: " << FNumeroCompresor << std::endl;
537  std::cout << "Tipo de error: " << N.what() << std::endl;
538  throw Exception("ERROR: EvaluaRCSpline en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
539  }
540 }
541 
542 double TMapaComp::EvaluaRCHermite(double Massflow) {
543  try {
544  double ret_val = 0., t2 = 0., t3 = 0., t = 0., h00 = 0., h10 = 0., h01 = 0., h11 = 0., h = 0.;
545 
546  if(Massflow < FGastoInt[0]) {
547  ret_val = FRelCompInt[0];
548  } else if(Massflow >= FGastoInt[FNumPuntos]) {
549  ret_val = 1.;
550  } else {
551  int k = 0;
552  int n = FNumPuntos;
553  //Determinacion del indice para evaluacion spline
554  while(FGastoInt[k] < Massflow && k < n) {
555  ++k;
556  }
557  h = (FGastoInt[k] - FGastoInt[k - 1]);
558  t = (Massflow - FGastoInt[k - 1]) / h;
559  t2 = t * t;
560  t3 = t2 * t;
561  h00 = 2 * t3 - 3 * t2 + 1;
562  h10 = t3 - 2 * t2 + t;
563  h01 = -2 * t3 + 3 * t2;
564  h11 = t3 - t2;
565  ret_val = h00 * FRelCompInt[k - 1] + h * h10 * FCoefSplRC[k - 1] + h01 * FRelCompInt[k] + h * h11 * FCoefSplRC[k];
566  }
567  return ret_val;
568 
569  } catch(exception &N) {
570  std::cout << "ERROR: EvaluaRCSpline en el compresor: " << FNumeroCompresor << std::endl;
571  std::cout << "Tipo de error: " << N.what() << std::endl;
572  throw Exception("ERROR: EvaluaRCSpline en el compresor: " + std::to_string(FNumeroCompresor) + N.what());
573  }
574 }
575 
576 //---------------------------------------------------------------------------
577 //-------- PolOrtogonal -----------------------------------------------------
578 // Obtiene los coeficiente 'b', 'c' y 'd' del polinomio ortogonal utilizado//
579 // para ajustar el rendimiento de compresor donde 'nterms' es el numero de //
580 // terminos que tienen los coeficientes, 'npoint' es el numero de puntos de//
581 // rendimiento que tiene una curva determinada 'ma' representa el massflow //
582 //'rd' el rendimiento para dicho massflow y 'w' es un vector cuyas componentes//
583 // valen 1 Ajuste por minimos cuadrados //
584 //---------------------------------------------------------------------------
585 
586 void TMapaComp::PolOrtogonal(int nterms, int npoint, double *ma, double *rd, double *w, double *b, double *c,
587  double *d) {
588  double p = 0;
589  try {
590  for(int j = 0; j < nterms; ++j) {
591  b[j] = 0.;
592  d[j] = 0.;
593  FOrtp[j].s = 0.;
594  }
595  c[0] = 0.;
596  for(int i = 0; i < npoint; ++i) {
597  d[0] += rd[i] * w[i];
598  b[0] += ma[i] * w[i];
599  FOrtp[0].s += w[i];
600  }
601  d[0] /= FOrtp[0].s;
602  for(int i = 0; i < npoint; ++i) {
603  FOrtp[i].error = rd[i] - d[0];
604  }
605  b[0] /= FOrtp[0].s;
606  for(int i = 0; i < npoint; ++i) {
607  FOrtp[i].pjm1 = 1.;
608  FOrtp[i].pj = ma[i] - b[0];
609  }
610  int j = 0;
611  while(j < nterms - 1) {
612  ++j;
613  for(int i = 0; i < npoint; ++i) {
614  p = FOrtp[i].pj * w[i];
615  d[j] += FOrtp[i].error * p;
616  p *= FOrtp[i].pj;
617  b[j] += ma[i] * p;
618  FOrtp[j].s += p;
619  }
620  d[j] /= FOrtp[j].s;
621  for(int i = 0; i < npoint; ++i) {
622  FOrtp[i].error -= d[j] * FOrtp[i].pj;
623  }
624  if(j < nterms) {
625  b[j] /= FOrtp[j].s;
626  c[j] = FOrtp[j].s / FOrtp[j - 1].s;
627  for(int i = 0; i < npoint; ++i) {
628  p = FOrtp[i].pj;
629  FOrtp[i].pj = (ma[i] - b[j]) * FOrtp[i].pj - c[j] * FOrtp[i].pjm1;
630  FOrtp[i].pjm1 = p;
631  }
632  }
633  }
634  } catch(exception &N) {
635  std::cout << "ERROR: PolOrtogonal en el compresor: " << FNumeroCompresor << std::endl;
636  std::cout << "Tipo de error: " << N.what() << std::endl;
637  throw Exception("ERROR: PolOrtogonal en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
638  }
639 }
640 
641 //---------------------------------------------------------------------------
642 //--------InterpolaMapa -----------------------------------------------------
643 // Ajusta los coeficientes de una spline que representa la curva de regimen//
644 // 'rtc' a temperatura ambiente 'AmbientTemperature' e interpola los coeficiente del //
645 // polinomio ortogonal que ajusta el rendimiento //
646 //---------------------------------------------------------------------------
647 
648 void TMapaComp::InterpolaMapa(double rtc, double AmbientTemperature) {
649 
650  double DeltaN = 0.;
651  int Cent = 0;
652 
653  try {
654 
655 // Se adapta el regimen de giro a interpolar al mapa compresor
656  FRegComp = rtc * sqrt(FTempRef / AmbientTemperature);
657 
658 // Se halla el regimen de giro inferior al que se quiere interpolar
659 
660  FCurvInf = 0;
661  for(int i = 0; i < FNumCurvasReg - 1; ++i) {
662  if(FRegComp >= FRegimenCurva[i]) {
663  FCurvInf = i;
664  }
665  }
666  DeltaN = (FRegComp - FRegimenCurva[FCurvInf]) / FIncReg;
667 
668 // Se interpola linealmente el massflow de bombeo y el de RC = 1
669 
670  if(FRegComp > FRegMax) {
671  FGastoBombeoX = FGastoBombeo[FCurvInf + 1] + (FGastoBombeo[FCurvInf + 1] - FGastoBombeo[FCurvInf]) / FIncReg *
672  (FRegComp - FRegMax);
673  } else if(FRegComp < FRegMin) {
674  FGastoBombeoX = FGastoBombeo[FCurvInf + 1] - (FGastoBombeo[FCurvInf + 2] - FGastoBombeo[FCurvInf + 1]) / FIncReg *
675  (FRegMin - FRegComp);
676  } else {
677  FGastoBombeoX = FGastoBombeo[FCurvInf + 1] + (FGastoBombeo[FCurvInf + 2] - FGastoBombeo[FCurvInf + 1]) * DeltaN;
678  }
679  if(FGastoBombeoX < 0)
680  FGastoBombeoX = 0.;
681 
682  FRelCompBombeoX = EvaluaHermite(FGastoBombeoX, FNumCurvasReg + 1, FGastoBombeo, FRelCompBombeo, FCoefSplBombeo);
683 
684  if(FRegComp > FRegMax) {
685  FGastoRelComp1X = FGastoRelComp1[FCurvInf] + (FGastoRelComp1[FCurvInf] - FGastoRelComp1[FCurvInf - 1]) / FIncReg *
686  (FRegComp - FRegMax);
687  } else if(FRegComp < FRegMin) {
688  FGastoRelComp1X = FGastoRelComp1[FCurvInf] - (FGastoRelComp1[FCurvInf + 1] - FGastoRelComp1[FCurvInf]) / FIncReg *
689  (FRegMin - FRegComp);
690  } else {
691  FGastoRelComp1X = FGastoRelComp1[FCurvInf] + (FGastoRelComp1[FCurvInf + 1] - FGastoRelComp1[FCurvInf]) * DeltaN;
692  }
693 
694 // Calculo de la curva interpolada Relacion de compresion/Masa corregida
695 
696 // Se haya el ultimo massflow del que se tiene informacion para el regimen de giro inferior
697 
698  if(FRegComp < FRegMin) {
699  FNumPuntos = floor(FGastoRelComp1X / FIncGasto) + 1;
700  } else {
701  int j = 0;
702  while(FRelComp[FCurvInf][j] > 0. && j < FNumPuntosGasto - 1) {
703  ++j;
704  }
705  FNumPuntos = j;
706  }
707 
708 // Se interpola linealmente la relacion de compresion entre el regimen de giro inferior
709 // y superior excepto para el ultimo punto.
710 
711  if(FRegComp < FRegMin) {
712  for(int j = 0; j < FNumPuntos; ++j) {
713  FGastoInt[j] = FGastoMin + FIncGasto * j;
714  FRelCompInt[j] = FRelComp[FCurvInf][j] - (FRelComp[FCurvInf + 1][j] - FRelComp[FCurvInf][j]) / FIncReg *
715  (FRegMin - FRegComp);
716  if(FRelCompInt[j] < 1.)
717  FRelCompInt[j] = 1.;
718  }
719  } else {
720  for(int j = 0; j < FNumPuntos; ++j) {
721  FGastoInt[j] = FGastoMin + FIncGasto * j;
722  FRelCompInt[j] = FRelComp[FCurvInf][j] + (FRelComp[FCurvInf + 1][j] - FRelComp[FCurvInf][j]) * DeltaN;
723  }
724  }
725  FGastoInt[FNumPuntos] = FGastoRelComp1X;
726  FRelCompInt[FNumPuntos] = 1.;
727 
728  FCurrentPresMAX = FRelCompInt[0];
729  for(int j = 1; j < FNumPuntos; ++j) {
730  if(FRelCompInt[j] > FCurrentPresMAX)
731  FCurrentPresMAX = FRelCompInt[j];
732  }
733 
734 // Se obtienen los coeficientes que ajusta la spline de la curva interpolada.
735 
736  Hermite(FNumPuntos, FGastoInt, FRelCompInt, FCoefSplRC);
737 
738 // Se interpolan linealmente los coeficientes de cada polinomio ortogonal para el
739 // regimen de giro buscado
740  /*
741  FNumTerms=3;
742  Cent=2;
743  if(FNumCurvasRen[FCurvInf]<FNumTerms){
744  FNumTerms=FNumCurvasRen[FCurvInf];
745  --Cent;
746  }
747  if(FNumCurvasRen[FCurvInf+1]<FNumTerms){
748  FNumTerms=FNumCurvasRen[FCurvInf+1];
749  --Cent;
750  }
751  if(Cent>FNumTerms) FNumTerms=3;
752 
753  for(int i=0;i<FNumTerms;++i){
754  if(FRegComp>FRegMax){
755  FCoefbX[i]=FCoefbSup[FCurvInf+1][i]+(FCoefbSup[FCurvInf+1][i]-FCoefbInf[FCurvInf][i])/FIncReg*(FRegComp-FRegMax);
756  FCoefcX[i]=FCoefcSup[FCurvInf+1][i]+(FCoefcSup[FCurvInf+1][i]-FCoefcInf[FCurvInf][i])/FIncReg*(FRegComp-FRegMax);
757  FCoefdX[i]=FCoefdSup[FCurvInf+1][i]+(FCoefdSup[FCurvInf+1][i]-FCoefdInf[FCurvInf][i])/FIncReg*(FRegComp-FRegMax);
758  }else if(FRegComp<FRegMin){
759  FCoefbX[i]=FCoefbInf[FCurvInf][i]-(FCoefbSup[FCurvInf+1][i]-FCoefbInf[FCurvInf][i])/FIncReg*(FRegMin-FRegComp);
760  FCoefcX[i]=FCoefbInf[FCurvInf][i]-(FCoefcSup[FCurvInf+1][i]-FCoefcInf[FCurvInf][i])/FIncReg*(FRegMin-FRegComp);
761  FCoefdX[i]=FCoefbInf[FCurvInf][i]-(FCoefdSup[FCurvInf+1][i]-FCoefdInf[FCurvInf][i])/FIncReg*(FRegMin-FRegComp);
762  }else{
763  FCoefbX[i]=FCoefbInf[FCurvInf][i]+(FCoefbSup[FCurvInf+1][i]-FCoefbInf[FCurvInf][i])*DeltaN;
764  FCoefcX[i]=FCoefbInf[FCurvInf][i]+(FCoefcSup[FCurvInf+1][i]-FCoefcInf[FCurvInf][i])*DeltaN;
765  FCoefdX[i]=FCoefbInf[FCurvInf][i]+(FCoefdSup[FCurvInf+1][i]-FCoefdInf[FCurvInf][i])*DeltaN;
766  }
767  } */
768  } catch(exception &N) {
769  std::cout << "ERROR: InterpolaMapa en el compresor: " << FNumeroCompresor << std::endl;
770  std::cout << "Tipo de error: " << N.what() << std::endl;
771  throw Exception("ERROR: InterpolaMapa en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
772  }
773 }
774 
775 //---------------------------------------------------------------------------
776 // ---------EvaluaRendimiento -----------------------------------------------
777 // Evalua el rendimiento a partir de los coeficientes del polinomio //
778 // ortogonal de la curva de funcionamiento para un valor de //
779 // masa de aire 'MasaAire' //
780 //---------------------------------------------------------------------------
781 
782 double TMapaComp::EvaluaRendimiento(double MasaAire) {
783  int k = 0;
784  double prev2 = 0.;
785  double prev = 0.;
786  double ret_val = 0.;
787  try {
788  k = FNumTerms - 1;
789  ret_val = FCoefdX[k];
790  prev = 0.;
791  --k;
792  while(k >= 0) {
793  prev2 = prev;
794  prev = ret_val;
795  ret_val = FCoefdX[k] + (MasaAire - FCoefbX[k]) * prev - FCoefcX[k + 1] * prev2;
796  --k;
797  }
798  return ret_val;
799  } catch(exception &N) {
800  std::cout << "ERROR: EvaluaRendimiento en el compresor: " << FNumeroCompresor << std::endl;
801  std::cout << "Tipo de error: " << N.what() << std::endl;
802  throw Exception("ERROR: EvaluaRendimiento en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
803  }
804 }
805 
806 //---------------------------------------------------------------------------
807 //---------------------------------------------------------------------------
808 
809 double TMapaComp::EvaluaRendSplines(double MasaAire) {
810  double AireAdim = 0.;
811  double ret_val = 0.;
812  double RendInf = 0.;
813  double RendSup = 0.;
814  double DeltaN = 0.;
815  try {
816  if(MasaAire < FGastoBombeoX) {
817  //En Bombeo se toma como rendimiento el del ultimo punto.
818  DeltaN = (FRegComp - FRegimenCurva[FCurvInf]) / FIncReg;
819 
820  RendInf = FRendAdim[FCurvInf][0];
821  RendSup = FRendAdim[FCurvInf + 1][0];
822 
823  if(FRegComp > FRegMax) {
824  ret_val = RendSup + (RendSup - RendInf) / FIncReg * (FRegComp - FRegMax);
825  } else if(FRegComp < FRegMin) {
826  ret_val = RendInf - (RendSup - RendInf) / FIncReg * (FRegMin - FRegComp);
827  } else {
828  ret_val = RendInf + (RendSup - RendInf) * DeltaN;
829  }
830  } else if(MasaAire > FGastoRelComp1X) {
831  ret_val = FRendGastoMaximo;
832  } else {
833  AireAdim = (MasaAire - FGastoBombeoX) / (FGastoRelComp1X - FGastoBombeoX);
834 
835  RendInf = EvaluaHermite(AireAdim, FNumCurvasRenAd[FCurvInf], FGastoAdim[FCurvInf], FRendAdim[FCurvInf],
836  FCoefSplRend[FCurvInf]);
837 
838  RendSup = EvaluaHermite(AireAdim, FNumCurvasRenAd[FCurvInf + 1], FGastoAdim[FCurvInf + 1], FRendAdim[FCurvInf + 1],
839  FCoefSplRend[FCurvInf + 1]);
840 
841  DeltaN = (FRegComp - FRegimenCurva[FCurvInf]) / FIncReg;
842 
843  if(FRegComp > FRegMax) {
844  ret_val = RendSup + (RendSup - RendInf) / FIncReg * (FRegComp - FRegMax);
845  } else if(FRegComp < FRegMin) {
846  ret_val = RendInf - (RendSup - RendInf) / FIncReg * (FRegMin - FRegComp);
847  } else {
848  ret_val = RendInf + (RendSup - RendInf) * DeltaN;
849  }
850  }
851  return ret_val;
852  } catch(exception &N) {
853  std::cout << "ERROR: EvaluaRendSplines en el compresor: " << FNumeroCompresor << std::endl;
854  std::cout << "Tipo de error: " << N.what() << std::endl;
855  throw Exception("ERROR: EvaluaRendSplines en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
856  }
857 }
858 
859 //---------------------------------------------------------------------------
860 //---------------------------------------------------------------------------
861 
862 double TMapaComp::GetRelCompInt(int i) {
863  try {
864  return FRelCompInt[i];
865  } catch(exception &N) {
866  std::cout << "ERROR: GetRelCompInt en el compresor: " << FNumeroCompresor << std::endl;
867  std::cout << "Tipo de error: " << N.what() << std::endl;
868  throw Exception("ERROR: GetRelCompInt en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
869  }
870 }
871 
872 //---------------------------------------------------------------------------
873 //---------------------------------------------------------------------------
874 
875 double TMapaComp::GetGastoInt(int i) {
876  try {
877  return FGastoInt[i];
878  } catch(exception &N) {
879  std::cout << "ERROR: GetGastoInt en el compresor: " << FNumeroCompresor << std::endl;
880  std::cout << "Tipo de error: " << N.what() << std::endl;
881  throw Exception("ERROR: GetGastoInt en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
882  }
883 }
884 
885 //---------------------------------------------------------------------------
886 //---------------------------------------------------------------------------
887 
888 double TMapaComp::BuscaRegimen(double RC, double Massflow, double AmbientTemperature) {
889  try {
890  double RCSup = 0.;
891  double RCInf = 0.;
892  double ret_val = 0., val1 = 0., val2 = 0.;
893  double reg = 0.;
894  int i = 0, sup, inf;
895  while(RC > RCSup && i < FNumCurvasReg) {
896  reg = FRegimenCurva[i] * AmbientTemperature / FTempRef;
897  InterpolaMapa(reg, AmbientTemperature);
898 
899  if(Massflow < FGastoRelComp1[i]) {
900  RCSup = EvaluaRCSplines(Massflow);
901  } else {
902  RCSup = 1.;
903  }
904  sup = i;
905  ++i;
906  }
907  if(RC > RCSup) {
908  std::cout << "WARNING: El Punto se sale por arriba del mapa " << std::endl;
909  //sup=i-1;
910  }
911  if(sup > 0) {
912  inf = sup - 1;
913  reg = FRegimenCurva[inf] * AmbientTemperature / FTempRef;
914  InterpolaMapa(reg, AmbientTemperature);
915 
916  if(Massflow < FGastoRelComp1[inf]) {
917  RCInf = EvaluaRCSplines(Massflow);
918  } else {
919  std::cout << "WARNING: El siguiente punto puede estar mal calculado:" << std::endl;
920  std::cout << " Relacionde compresion: " << RC << std::endl;
921  std::cout << " Massflow masico: " << Massflow << std::endl;
922  std::cout << " Deberias aumentar el numero de curvas del mapa" << std::endl << std::endl;
923  RCInf = 1.;
924  }
925  ret_val = (RC - RCInf) / (RCSup - RCInf) * (FRegimenCurva[sup] - FRegimenCurva[inf]) + FRegimenCurva[inf];
926  } else {
927  inf = sup;
928  sup = inf + 1;
929  RCInf = RCSup;
930  InterpolaMapa(FRegimenCurva[sup], __units::degCToK(AmbientTemperature));
931  RCSup = EvaluaRCSplines(Massflow);
932  val1 = (RC - RCInf) / (RCSup - RCInf) * (FRegimenCurva[sup] - FRegimenCurva[inf]) + FRegimenCurva[inf];
933 
934  val2 = (RC - 1) / (RCInf - 1) * (FRegimenCurva[inf]);
935  ret_val = (val1 + val2) / 2;
936  std::cout << "WARNING: El siguiente punto puede estar mal calculado:" << std::endl;
937  std::cout << " Relacionde compresion: " << RC << std::endl;
938  std::cout << " Massflow masico: " << Massflow << std::endl;
939  std::cout << " Deberias incluir curvas con menor regimen de giro" << std::endl << std::endl;
940  }
941  return ret_val;
942  } catch(exception &N) {
943  std::cout << "ERROR: BuscaRegimen en el compresor: " << FNumeroCompresor << std::endl;
944  std::cout << "Tipo de error: " << N.what() << std::endl;
945  throw Exception("ERROR: BuscaRegimen en el compresor: "/* +std::to_string(FNumeroCompresor)+N.what()*/);
946  }
947 }
948 
949 //---------------------------------------------------------------------------
950 //---------------------------------------------------------------------------
951 
952 void TMapaComp::ImprimeMapa() {
953  std::cout << "Printing compresso map .";
954  FILE *fich;
955  fich = fopen("MapaSimple.txt", "w");
956  double inc = 0.;
957  double *massflow;
958  double **rc;
959  massflow = new double[200];
960  rc = new double*[200];
961  for(int j = 0; j < 200; j++) {
962  rc[j] = new double[FNumCurvasReg];
963  }
964  inc = (FGastoMax - FGastoMin) / 200;
965  for(int j = 0; j < 200; j++) {
966  massflow[j] = FGastoMin + inc * (double) j;
967  }
968  for(int i = 0; i < FNumCurvasReg; i++) {
969  std::cout << ".";
970  InterpolaMapa(FRegMin + FIncReg * (double) i, FTempRef);
971  for(int j = 0; j < 200; j++) {
972  rc[j][i] = EvaluaRCSplines(massflow[j]);
973  }
974  }
975 
976  for(int j = 0; j < 200; j++) {
977  if(j < FNumCurvasReg) {
978  fprintf(fich, "%lf\t%lf\t", FGastoBombeo[j], FRelCompBombeo[j]);
979  } else {
980  fprintf(fich, "\t\t");
981  }
982  fprintf(fich, "%lf\t", massflow[j]);
983  for(int i = 0; i < FNumCurvasReg; i++) {
984  fprintf(fich, "%lf\t", rc[j][i]);
985  }
986  fprintf(fich, "\n");
987  }
988  delete[] massflow;
989  for(int j = 0; j < 200; j++) {
990  delete[] rc[j];
991  }
992  fclose(fich);
993  std::cout << "Finished" << std::endl;
994 }
995 
996 //---------------------------------------------------------------------------
997 
998 #pragma package(smart_init)
stSpline
Definition: Globales.h:911
pow3
T pow3(T x)
Returns x to the power of 3.
Definition: Math_wam.h:101
stOrtoPol
Definition: Globales.h:925
TCompressorMap
Definition: TCompressorMap.h:12
Exception
Custom exception class.
Definition: Exception.hpp:39