OpenWAM
Tfql.cpp
1 /*--------------------------------------------------------------------------------*\
2 ==========================|
3  \\ /\ /\ // O pen | OpenWAM: The Open Source 1D Gas-Dynamic Code
4  \\ | X | // W ave |
5  \\ \/_\/ // A ction | CMT-Motores Termicos / Universidad Politecnica Valencia
6  \\/ \// M odel |
7  ----------------------------------------------------------------------------------
8  License
9 
10  This file is part of OpenWAM.
11 
12  OpenWAM is free software: you can redistribute it and/or modify
13  it under the terms of the GNU General Public License as published by
14  the Free Software Foundation, either version 3 of the License, or
15  (at your option) any later version.
16 
17  OpenWAM is distributed in the hope that it will be useful,
18  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20  GNU General Public License for more details.
21 
22  You should have received a copy of the GNU General Public License
23  along with OpenWAM. If not, see <http://www.gnu.org/licenses/>.
24 
25 
26  \*--------------------------------------------------------------------------------*/
27 
28 //---------------------------------------------------------------------------
29 #pragma hdrstop
30 
31 #include "Tfql.h"
32 
33 //---------------------------------------------------------------------------
34 //---------------------------------------------------------------------------
35 
36 Tfql::Tfql(int ncilin) {
37 
38  Fncilin = ncilin;
39 
40 }
41 
42 //---------------------------------------------------------------------------
43 //---------------------------------------------------------------------------
44 
45 Tfql::~Tfql() {
46  if(Fmax != NULL)
47  delete[] Fmax;
48 
49  if(Fpar_dist != NULL) {
50  for(int i = 0; i < Fnparametros; ++i) {
51  delete[] Fpar_dist[i];
52  }
53  delete[] Fpar_dist;
54  }
55 
56  if(Flm != NULL) {
57  for(int i = 0; i < Fnwiebe; ++i) {
58  delete[] Flm[i];
59  }
60  delete[] Flm;
61  }
62 
63  if(Flc != NULL) {
64  for(int i = 0; i < Fnwiebe; ++i) {
65  delete[] Flc[i];
66  }
67  delete[] Flc;
68  }
69 
70  if(Flb != NULL) {
71  for(int i = 0; i < Fnwiebe; ++i) {
72  delete[] Flb[i];
73  }
74  delete[] Flb;
75  }
76 
77  if(Fla != NULL) {
78  for(int i = 0; i < Fnwiebe; ++i) {
79  delete[] Fla[i];
80  }
81  delete[] Fla;
82  }
83 
84  if(Fli != NULL) {
85  for(int i = 0; i < Fnwiebe; ++i) {
86  delete[] Fli[i];
87  }
88  delete[] Fli;
89  }
90 
91  if(Flab != NULL) {
92  for(int i = 0; i < Fnwiebe; ++i) {
93  for(int j = 0; j < Fnley; ++j) {
94  delete[] Flab[i][j];
95  }
96  delete[] Flab[i];
97  }
98  delete[] Flab;
99  }
100 
101  if(Fcombustion != NULL)
102  delete[] Fcombustion;
103 
104 }
105 
106 //---------------------------------------------------------------------------
107 //---------------------------------------------------------------------------
108 
109 void Tfql::lee_leylib(char *Ruta, FILE *fich) {
110  char *Filefql;
111  char *Datosfql;
112 
113  try {
114 
115  for(int i = 0; i <= (int) strlen(Ruta); i++) {
116  Datosfql[i] = Ruta[i];
117  }
118 
119  fscanf(fich, "%s ", &Filefql);
120  strcat(Datosfql, Filefql);
121 
122  Fichfql = fopen(Datosfql, "r");
123  if((Fichfql = fopen(Datosfql, "r")) == NULL) {
124  std::cout << "ERROR: Fichero de leyes de liberacion de calor no cargado";
125  } else {
126  fscanf(Fichfql, "%d ", &Fnley);
127  fscanf(Fichfql, "%d ", &Fnwiebe);
128  fscanf(Fichfql, "%d ", &Fnparametros);
129  Fpar_dist = new double *[Fnparametros];
130  for(int i = 0; i <= Fnparametros - 1; ++i) {
131  Fpar_dist[i] = new double[Fnley];
132  }
133  Fmax = new double[Fnparametros];
134  for(int k = 0; k <= Fnparametros - 1; k++) {
135  Fmax[k] = 0;
136  }
137  Flm = new double*[Fnwiebe];
138  for(int i = 0; i <= Fnwiebe - 1; ++i) {
139  Flm[i] = new double[Fnley];
140  }
141  Flc = new double*[Fnwiebe];
142  for(int i = 0; i <= Fnwiebe - 1; ++i) {
143  Flc[i] = new double[Fnley];
144  }
145  Flb = new double*[Fnwiebe];
146  for(int i = 0; i <= Fnwiebe - 1; ++i) {
147  Flb[i] = new double[Fnley];
148  }
149  Fli = new double*[Fnwiebe];
150  for(int i = 0; i <= Fnwiebe - 1; ++i) {
151  Fli[i] = new double[Fnley];
152  }
153  Fla = new double*[Fnwiebe];
154  for(int i = 0; i <= Fnwiebe - 1; ++i) {
155  Fla[i] = new double[Fnley];
156  }
157 
158  Flab = new double **[Fnwiebe];
159  for(int i = 0; i <= Fnwiebe - 1; ++i) {
160  Flab[i] = new double *[Fnley];
161  for(int j = 0; j <= Fnley - 1; ++j) {
162  Flab[i][j] = new double[Fncilin];
163  }
164  }
165 
166  Fcombustion = new bool[Fncilin];
167  for(int i = 0; i <= Fncilin - 1; ++i) {
168  Fcombustion[i] = false;
169  }
170 
171  printf("INFO: N. de leyes: %d\n", Fnley);
172  printf("INFO: N. de Wiebes: %d\n", Fnwiebe);
173 
174  for(int i = 0; i < Fnley; ++i) {
175  for(int j = 0; j <= Fnparametros - 1; ++j) {
176  fscanf(Fichfql, "%lf ", &Fpar_dist[j][i]);
177  if(Fmax[j] < fabs(Fpar_dist[j][i])) {
178  Fmax[j] = fabs(Fpar_dist[j][i]);
179  }
180  }
181  for(int j = 0; j <= Fnwiebe - 1; ++j) {
182  fscanf(Fichfql, "%lf %lf %lf %lf %lf ", &Flm[j][i], &Flc[j][i], &Flb[j][i], &Fli[j][i], &Fla[j][i]);
183  Fla[j][i] = 720. + Fla[j][i];
184  }
185 
186  }
187  }
188  fclose(Fichfql);
189  } catch(exception &N) {
190  std::cout << "ERROR: lee_leylib (DLL)" << std::endl;
191  std::cout << "Tipo de error: " << N.what() << std::endl;
192  throw Exception(N.what());
193  }
194 }
195 
196 //---------------------------------------------------------------------------
197 //---------------------------------------------------------------------------
198 
199 double Tfql::fql(double x, int j, int i) {
200 
201  try {
202 
203  double ret_val, wiebe = 0., c;
204 
205  ret_val = 0.;
206 
207  if(x > 720) {
208  c = x / 720;
209  x = x - 720 * floor(c);
210  }
211 
212  /* Depende del orden de encendido */
213  /* Para un cuatro cilindros dW10b*/
214  //if(i==0) x=x;
215  if(i == 1)
216  x = x - 540;
217  if(i == 2)
218  x = x - 180;
219  if(i == 3)
220  x = x - 360;
221 
222  /* Para un seis cilindros VOLVO de Vicente
223  if(i==0) x=x;
224  if(i==1) x=x-480;
225  if(i==2) x=x-240;
226  if(i==3) x=x-600;
227  if(i==4) x=x-120;
228  if(i==5) x=x-360;
229  */
230  while(x < 180)
231  x = x + 720;
232  while(x > 900)
233  x = x - 720;
234 
235  if(x > Fang0 + 719.75 && x < Fang0 + 720.25)
236  Fcombustion[i] = true;
237 
238  if(Fcombustion[i]) {
239  for(int k = 0; k <= Fnwiebe - 1; ++k) {
240  wiebe = fun_wiebe(x, Flm[k][j], Flc[k][j], Fli[k][j], Flab[k][j][i]);
241  ret_val += wiebe * Flb[k][j];
242  }
243  } else {
244  ret_val = 0.;
245  }
246  return ret_val;
247 
248  }
249 
250  catch(exception &N) {
251  std::cout << "ERROR: fql (DLL)" << std::endl;
252  std::cout << "Tipo de error: " << N.what() << std::endl;
253  throw Exception(N.what());
254  }
255 }
256 
257 //---------------------------------------------------------------------------
258 //---------------------------------------------------------------------------
259 
260 double Tfql::fun_wiebe(double x, double m, double c, double ia, double a0)
261 
262 {
263 
264  try {
265 
266  double ret_val = 0., xx = 0., xxx = 0.;
267 
268  if(x <= a0) {
269  ret_val = 0;
270  } else {
271 
272  xx = (x - a0) / ia;
273  xxx = pow(xx, m + 1);
274  if((xxx * c) > 10.) { /* Josevi */
275  ret_val = 1.;
276  } else {
277  ret_val = 1. - 1. / exp(xxx * c);
278  }
279 
280  }
281  return ret_val;
282  }
283 
284  catch(exception &N) {
285  std::cout << "ERROR: fun_wiebe (DLL)" << std::endl;
286  std::cout << "Tipo de error: " << N.what() << std::endl;
287  throw Exception(N.what());
288  }
289 }
290 
291 //---------------------------------------------------------------------------
292 //---------------------------------------------------------------------------
293 
294 void Tfql::calcula_angulos_combustion(double *parametros, int i) {
295 
296  try {
297  int k = 0;
298  double up, down, ang01 = 0., tras = 0., a = 0., b = 0., c = 0., r = 0., s = 0., d = 0., densidad = 0., egr = 0., e;
299  double dist = 0.;
300 
301  up = 0.;
302  down = 0.;
303 
304  if(Fnwiebe == 4) {
305  k = 1;
306  } else {
307  k = 0;
308  }
309 
310  for(int j = 0; j < Fnley; ++j) {
311  b = 0;
312  /* for(int i=0;i<=Fnparametros-1;++i){
313  a=pow((parametros[i]-Fpar_dist[i][j])/Fmax[i],2.);
314  b+=a;
315  } */
316 
317  //c=pow((parametros[1]-Fpar_dist[1][j])/Fmax[1],2.);
318  r = pow2((parametros[2] - Fpar_dist[2][j]) / Fmax[2]); // Regimen
319  s = pow2((parametros[5] - Fpar_dist[5][j]) / Fmax[5]); // SOI
320  d = pow2((parametros[6] - Fpar_dist[6][j]) / Fmax[6]); // Dosado absoluto
321  egr = pow2((parametros[7] - Fpar_dist[7][j]) / Fmax[7]); // Massflow de EGR
322  densidad = pow2((parametros[8] - Fpar_dist[8][j]) / Fmax[8]); // Density en 2
323  b = c + r + s + d + egr + densidad;
324 
325  dist = sqrt(b);
326  if(dist < 1e-15)
327  dist = 1e-15;
328  e = 2 + 1 / dist;
329  if(b < 0.008) {
330  dist = 1e-15;
331  } else {
332  dist = pow(b, e);
333  }
334  if(dist < 1e-15)
335  dist = 1e-15;
336  up += (Fla[k][j] / dist);
337  down += (1 / dist);
338  }
339  ang01 = up / down;
340 
341  Fang0 = ang01 - 720.;
342  Ffinc = ang01 - 720.;
343 
344  for(int j = 0; j < Fnley; ++j) {
345  tras = ang01 - Fla[k][j];
346  for(int n = 0; n <= Fnwiebe - 1; ++n) {
347  Flab[n][j][i] = Fla[n][j] + tras;
348  }
349  if(Flab[0][j][i] - 720. < Fang0) {
350  Fang0 = Flab[0][j][i] - 720.; // ?????????? anado -720
351  }
352  if((Flab[Fnwiebe - 1][j][i] + Fli[Fnwiebe - 1][j] - 720.) > Ffinc) {
353  Ffinc = Flab[Fnwiebe - 1][j][i] + Fli[Fnwiebe - 1][j] - 720.;
354  }
355  }
356  }
357 
358  catch(exception &N) {
359  std::cout << "ERROR: fql (DLL)" << std::endl;
360  std::cout << "Tipo de error: " << N.what() << std::endl;
361  throw Exception(N.what());
362  }
363 }
364 
365 //---------------------------------------------------------------------------
366 //---------------------------------------------------------------------------
367 
368 double Tfql::calcula_fql(double *parametros, double x, int i) {
369 
370  try {
371  double up, down, dist = 0., ret_val = 0., a = 0., b = 0., ley = 0., e;
372  int j = 0;
373 
374  up = 0.;
375  down = 0.;
376 
377  for(j = 0; j < Fnley; ++j) {
378  b = 0;
379  for(int k = 0; k < Fnparametros - 4; ++k) {
380  a = pow2((parametros[k] - Fpar_dist[k][j]) / Fmax[k]);
381  b += a;
382  }
383 
384  dist = sqrt(b);
385  e = 2 + 1 / dist;
386  dist = pow(b, e);
387  ley = fql(x, j, i);
388  if(ley > 1.)
389  printf("WARNING: Ley n. %d > 1. en %lf\n", j, x);
390  up += ley / dist;
391  down += 1. / dist;
392 
393  }
394 
395  ret_val = up / down;
396 
397  return ret_val;
398 
399  }
400 
401  catch(exception &N) {
402  std::cout << "ERROR: calcula_fql (DLL)" << std::endl;
403  std::cout << "Tipo de error: " << N.what() << std::endl;
404  throw Exception(N.what());
405  }
406 }
407 
408 //-----------------------------------------------------------------------------
409 //-----------------------------------------------------------------------------
410 
411 void Tfql::lee_leylib2(FILE *BaseDatos) {
412 
413  fscanf(BaseDatos, "%d ", &Fnley);
414  fscanf(BaseDatos, "%d ", &Fnwiebe);
415  fscanf(BaseDatos, "%d ", &Fnparametros);
416  Fpar_dist = new double *[Fnparametros];
417  for(int i = 0; i <= Fnparametros - 1; ++i) {
418  Fpar_dist[i] = new double[Fnley];
419  }
420  Fmax = new double[Fnparametros];
421  for(int k = 0; k <= Fnparametros - 1; k++) {
422  Fmax[k] = 0;
423  }
424  Flm = new double*[Fnwiebe];
425  for(int i = 0; i <= Fnwiebe - 1; ++i) {
426  Flm[i] = new double[Fnley];
427  }
428  Flc = new double*[Fnwiebe];
429  for(int i = 0; i <= Fnwiebe - 1; ++i) {
430  Flc[i] = new double[Fnley];
431  }
432  Flb = new double*[Fnwiebe];
433  for(int i = 0; i <= Fnwiebe - 1; ++i) {
434  Flb[i] = new double[Fnley];
435  }
436  Fli = new double*[Fnwiebe];
437  for(int i = 0; i <= Fnwiebe - 1; ++i) {
438  Fli[i] = new double[Fnley];
439  }
440  Fla = new double*[Fnwiebe];
441  for(int i = 0; i <= Fnwiebe - 1; ++i) {
442  Fla[i] = new double[Fnley];
443  }
444  Flab = new double **[Fnwiebe];
445  for(int i = 0; i <= Fnwiebe - 1; ++i) {
446  Flab[i] = new double *[Fnley];
447  for(int j = 0; j <= Fnley - 1; ++j) {
448  Flab[i][j] = new double[Fncilin];
449  }
450  }
451 
452  printf("INFO: N. de leyes: %d\n", Fnley);
453  printf("INFO: N. de Wiebes: %d\n", Fnwiebe);
454 
455  for(int i = 0; i < Fnley; ++i) {
456  for(int j = 0; j <= Fnparametros - 1; ++j) {
457  fscanf(BaseDatos, "%lf ", &Fpar_dist[j][i]);
458  if(Fmax[j] < fabs(Fpar_dist[j][i])) {
459  Fmax[j] = fabs(Fpar_dist[j][i]);
460  }
461  }
462  for(int j = 0; j <= Fnwiebe - 1; ++j) {
463  fscanf(BaseDatos, "%lf %lf %lf %lf %lf ", &Flm[j][i], &Flc[j][i], &Flb[j][i], &Fli[j][i], &Fla[j][i]);
464  Fla[j][i] = 720. + Fla[j][i];
465  }
466 
467  }
468 
469 }
470 
471 #pragma package(smart_init)
472 
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