OpenWAM
TCilindro4T.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 "TCilindro4T.h"
32 #include "TBloqueMotor.h"
33 #include "TCCCilindro.h"
34 #include "TTubo.h"
35 #include "Constantes.h"
36 
37 //#include <cmath>
38 // ---------------------------------------------------------------------------
39 // ---------------------------------------------------------------------------
40 
41 TCilindro4T::TCilindro4T(TBloqueMotor *Engine, int nc, bool ThereIsEGR) :
42  TCilindro(Engine, ThereIsEGR) {
43 
44  FNumeroCilindro = nc;
45  FAcumMasaPorAdm = 0;
46  FAcumMasaPorEsc = 0;
47  FAcumMasaEGR = 0;
48 
49  FCalor.Liberado = 0;
50 
51  FAnguloRetrasoCombustion = 2.;
52  FPrimerInstanteCicloCerrado = false;
53  FNumIny = 0;
54 }
55 
56 // ---------------------------------------------------------------------------
57 // ---------------------------------------------------------------------------
58 
59 TCilindro4T::~TCilindro4T() {
60 
61 }
62 
63 // ---------------------------------------------------------------------------
64 // ---------------------------------------------------------------------------
65 
66 void TCilindro4T::ActualizaPropiedades(double TiempoActual) {
67  try {
68  double masavalesc, MasaAdmInstante = 0., MasaEscInstante = 0.;
69  double FraccionCC = 0., MasaCortocircuitoAdm = 0., MasaCortocircuitoEsc = 0.;
70  int NumeroUnionesEntrante = 0;
71  double FraccionMasicaAcum = 0., mfquefin = 0.;
72  double z = 0.;
73  double ctorbadmp = 0.;
74 
75  if(FCicloCerrado & FPrimerInstanteCicloCerrado) {
76 
77  FPrimerInstanteCicloCerrado = false;
78 
79  if(FMfControlled)
80  FMasaFuel = FMfController->Output(TiempoActual);
81  else if(FMotor->getCombustible() == nmMEC)
82  FMasaFuel = FMotor->getMasaFuel();
83 
84  // Calculo Inicio y Fin de la Combustion.
85  FMfint = FMasaFuel; // kg/cc
86  FMaint = FMasaPorAdmision; // kg/cc
87  FRegInt = FMotor->getRegimen();
88 
89  if(FMotor->getACT()) {
90 
91  FMotor->NewInjectionData(TiempoActual);
92 
93  VariableInicialesCicloACT();
94 
95  for(int i = 0; i < FIN; i++) {
96  FSOP[i] = FMotor->getInjecPulse(i).Angulo;
97  FMFI[i] = FMotor->getInjecPulse(i).Masa * FMotor->getMasaFuel() * 1e6;
98 
99  }
100 
101  CALCULUS_OF_INJECTION_RATE(FIN, FSOP, FMFI, FSOI, FEOI, Ftest_variables[0], FCAI, Ftest_variables[5], FA_TASA, FB_TASA,
102  FC_TASA, FD_TASA, Finjection_rate, FCAD_injection_rate);
103 
104  ACT(Fengine_parameters, Fengine_model_constants, Ftest_variables, Finjection_rate, FCAD_injection_rate, Fsize_inlet_inj,
105  FIN, FSOI, FEOI, FCAI, FCAD_exit, FHRF_exit, FROHR_exit,
106  Fp_cyl_exit, Fdp_da_cyl_exit, FT_cyl_exit, FH_cooler_exit, Fmean_var_exit, Fheat_transfer, Finjection_rate_exit,
107  Faccum_injection_rate_exit, dataIN, &dataOUT);
108 
109  // Normalize total mass fraction.
110  double ACT_Correction = 0;
111  for(int i = 0; i < 8; i++) {
112  ACT_Correction += dataOUT.species_EVO[i];
113  }
114  for(int i = 0; i < 8; i++) {
115  dataOUT.species_EVO[i] = dataOUT.species_EVO[i] / ACT_Correction;
116  }
117 
118  FPresionMedAdm = 0.;
119  FPresionMedEsc = 0.;
120  }
121 
122  if(FCalcComb == nmFQL) {
123  InicioFinCombustion();
124  if(FMotor->getFTipoDatosIny() == 1) {
125  FNumIny = FMotor->getFNumeroInyecciones();
126  FTInyeccion = FMotor->getFTIny();
127  FPercentInyeccion = FMotor->getFPercentIny();
128  FAnguloInjeccion = FMotor->getFAngIny();
129  for(int i = 0; i < FNumIny; i++) {
130  FTInyeccion[i] = FTInyeccion[i] * 1e-3;
131  FPercentInyeccion[i] = FPercentInyeccion[i] / 100;
132  //Si la inyección empieza en angulo negativo, hay que sumar 720º
133  if(FAnguloInjeccion[i] < 0) {
134  FAnguloInjeccion[i] = FAnguloInjeccion[i] + FMotor->getAngTotalCiclo();
135  }
136 // else {
137 // FAnguloInjeccion[i] = FAnguloInjeccion[i];
138 // }
139  }
140  } else if(FMotor->getFTipoDatosIny() == 2) {
141  FNumIny = 1.;
142  FAnguloInjeccion.resize(FNumIny);
143  FTInyeccion.resize(FNumIny);
144  FPercentInyeccion.resize(FNumIny);
145  if(FMotor->getFAngIniIny() < 0) {
146  FAnguloInjeccion[0] = FMotor->getFAngIniIny() + FMotor->getAngTotalCiclo();
147  } else {
148  FAnguloInjeccion[0] = FMotor->getFAngIniIny();
149  }
150  FPercentInyeccion[0] = 1.;
151  } else {
152  // Si no hay datos del angulo de la inyección, se estiman de la FQL
153  //Solo hay inyección piloto si hay 4 Wiebes
154  if(FMotor->getLeyQuemadoBD()[0].Wiebes.size() == 4) {
155  FNumIny = 2.;
156  FAnguloInjeccion.resize(FNumIny);
157  FTInyeccion.resize(FNumIny);
158  FPercentInyeccion.resize(FNumIny);
159  FTInyeccion[0] = 5 / (FMotor->getRegimen() * 6.);
160  FTInyeccion[1] = (FFinComb - FIniComb) / (FMotor->getRegimen() * 6.) / 4.;
161  FPercentInyeccion[0] = FMotor->getLeyQuemadoBD()[0].Wiebes[0].Beta;
162  FPercentInyeccion[1] = 1. - FPercentInyeccion[0];
163  //Si la combustión principal empieza en angulo negativo, hay que sumar 720º
164  FAnguloInjeccion[1] = FIniComb - FAnguloRetrasoCombustion;
165  if(FAnguloInjeccion[1] < 0)
166  FAnguloInjeccion[1] += FMotor->getAngTotalCiclo();
167  FAnguloInjeccion[0] = FIniComb - 20;
168  if(FAnguloInjeccion[0] < 0)
169  FAnguloInjeccion[0] += FMotor->getAngTotalCiclo();
170 
171  } else {
172  FNumIny = 1;
173  FAnguloInjeccion.resize(FNumIny);
174  FTInyeccion.resize(FNumIny);
175  FPercentInyeccion.resize(FNumIny);
176  if(FMotor->getCombustible() == nmMEC) {
177  FTInyeccion[0] = (FFinComb - FIniComb) / (FMotor->getRegimen() * 6.) / 4.;
178  FAnguloInjeccion[0] = FIniComb - FAnguloRetrasoCombustion;
179  } else {
180  FTInyeccion[0] = (720 - FIniComb - FDistribucion.CA) / (FMotor->getRegimen() * 6.) / 2.;
181  FAnguloInjeccion[0] = FDistribucion.CA + 2.;
182  }
183  if(FAnguloInjeccion[0] < 0)
184  FAnguloInjeccion[0] += FMotor->getAngTotalCiclo();
185  FPercentInyeccion[0] = 1;
186  }
187  }
188  } else if(FCalcComb == nmACT) {
189  FNumIny = 1;
190  FAnguloInjeccion[0] = FMotor->getInjecPulse(0).Angulo;
191  FIniComb = FMotor->getInjecPulse(0).Angulo;
192  FFinComb = FDistribucion.AE;
193  }
194  // Queda una ultima opcion, que se realice el calculo con la dll. En ese caso se el valor a
195  // FIniComb y FFinComb ya se le habra dado a partir de la DLL con una property write.
196 
197  std::cout << "INFO: Fuel mass: " << FMasaFuel * 1e6 << " (mg)" << std::endl;
198  std::cout << "____________________________________________________" << std::endl;
199  }
200 
201  if(FMotor->getGammaCalculation() == nmGammaConstante) {
202  if(FMotor->getSpeciesModel() == nmCalculoCompleto) {
203  if(FMotor->getSpeciesNumber() == 9) {
204  FFraccionMasicaEspecieFuel = 0; // No se tiene en cuenta el combustible
205  } else if(FMotor->getSpeciesNumber() == 10) {
206  FFraccionMasicaEspecieFuel = FFraccionMasicaEspecie[7];
207  }
208  FRMezcla = CalculoCompletoRMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2],
209  FFraccionMasicaEspecieFuel, nmComposicionTemperatura,
210  FMotor->getCombustible());
211  FCpMezcla = CalculoCompletoCpMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2],
212  FFraccionMasicaEspecieFuel, __units::degCToK(FTemperature),
213  nmComposicionTemperatura, FMotor->getCombustible());
214  FGamma = CalculoCompletoGamma(FRMezcla, FCpMezcla, nmComposicionTemperatura);
215 
216  } else if(FMotor->getSpeciesModel() == nmCalculoSimple) {
217  FRMezcla = CalculoSimpleRMezcla(FComposicionCicloCerrado[0], FComposicionCicloCerrado[1], FMotor->getGammaCalculation(),
218  FMotor->getCombustible());
219  //FRMezcla = 287*FComposicionCicloCerrado[2] + 55.95*FComposicionCicloCerrado[1] + 285.4*FComposicionCicloCerrado[0];
220  FCvMezcla = CalculoSimpleCvMezcla(__units::degCToK(FTemperature), FComposicionCicloCerrado[0],
221  FComposicionCicloCerrado[1], nmComposicionTemperatura, FMotor->getCombustible());
222  FGamma = CalculoSimpleGamma(FRMezcla, FCvMezcla, nmComposicionTemperatura);
223  }
224 
225  } else if(FMotor->getGammaCalculation() == nmComposicion || FMotor->getGammaCalculation() == nmComposicionTemperatura) {
226  if(FMotor->getSpeciesModel() == nmCalculoCompleto) {
227  if(FMotor->getSpeciesNumber() == 9) {
228  FFraccionMasicaEspecieFuel = 0; // No se tiene en cuenta el combustible
229  } else if(FMotor->getSpeciesNumber() == 10) {
230  FFraccionMasicaEspecieFuel = FFraccionMasicaEspecie[7];
231  }
232  FRMezcla = CalculoCompletoRMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2],
233  FFraccionMasicaEspecieFuel, FMotor->getGammaCalculation(),
234  FMotor->getCombustible());
235  FCpMezcla = CalculoCompletoCpMezcla(FFraccionMasicaEspecie[0], FFraccionMasicaEspecie[1], FFraccionMasicaEspecie[2],
236  FFraccionMasicaEspecieFuel, __units::degCToK(FTemperature),
237  FMotor->getGammaCalculation(), FMotor->getCombustible());
238  FGamma = CalculoCompletoGamma(FRMezcla, FCpMezcla, FMotor->getGammaCalculation());
239 
240  } else if(FMotor->getSpeciesModel() == nmCalculoSimple) {
241 
242  FRMezcla = CalculoSimpleRMezcla(FComposicionCicloCerrado[0], FComposicionCicloCerrado[1], FMotor->getGammaCalculation(),
243  FMotor->getCombustible());
244  //FRMezcla = 287*FComposicionCicloCerrado[2] + 55.95*FComposicionCicloCerrado[1] + 285.4*FComposicionCicloCerrado[0];
245 
246  FCvMezcla = CalculoSimpleCvMezcla(__units::degCToK(FTemperature), FComposicionCicloCerrado[0],
247  FComposicionCicloCerrado[1], FMotor->getGammaCalculation(), FMotor->getCombustible());
248  FGamma = CalculoSimpleGamma(FRMezcla, FCvMezcla, FMotor->getGammaCalculation());
249 
250  }
251  }
252  FGamma1 = __Gamma::G1(FGamma);
253  FGamma2 = __Gamma::G2(FGamma);
254  FGamma4 = __Gamma::G4(FGamma);
255  FGamma6 = __Gamma::G6(FGamma);
256 
257  /* if(FNumeroCilindro==4){
258  printf(" %lf %lf\n", FAnguloActual,FGamma);
259  }
260 
261  if(FNumeroCilindro==4 & FAnguloActual>355.35){
262  printf(" \n");
263  } */
264 
265  // ACTUALIZAMOS LOS VALORES DE TIEMPO Y ANGULO AL TIEMPO ACTUAL
266  FTime0 = FTime1;
267  FTime1 = TiempoActual;
268  FDeltaT = FTime1 - FTime0;
269  FDeltaAngulo = 360. * FMotor->getRegimen() / 60. * FDeltaT;
270  FAnguloAnterior = FAnguloActual;
271  FAnguloActual = FAnguloAnterior + FDeltaAngulo;
272  if(FAnguloActual > FMotor->getAngTotalCiclo()) {
273  FAnguloActual -= FMotor->getAngTotalCiclo();
274  FNumeroCiclo++;
275  // std::cout << "INFO: El cilindro " << FNumeroCilindro << " comienza el ciclo " << FNumeroCiclo << std::endl;
276  }
277  // SE CENTRA LA COMBUSTION EN 0
278  if(FAnguloActual > 360.) {
279  FAnguloComb0 = FAnguloActual - FMotor->getAngTotalCiclo() - FDeltaAngulo;
280  FAnguloComb = FAnguloActual - FMotor->getAngTotalCiclo();
281  } else {
282  FAnguloComb0 = FAnguloComb;
283  FAnguloComb = FAnguloActual;
284  }
285 
286  if(FAnguloActual >= FDistribucion.AE && FAnguloAnterior < FDistribucion.AE) {
287  FCicloCerrado = false;
288  FResMediosCilindro.CalorCombustionSUM = FCalor.LiberadoTotal;
289  FCalor.Liberado = 0.;
290  FCalor.LiberadoTotal = 0.;
291  FCalor.FQL = 0.;
292  FCalor.FQL0 = 0.;
293  FFuelTotal = 0;
294  }
295 
296  // CALCULO DE LA PRESION MEDIA DE ADMISION Y DE ESCAPE PARA ACT
297  if(FMotor->getACT()) {
298  for(int i = 0; i < FNumeroUnionesAdm; i++) {
299  FPresionMedAdm += FCCValvulaAdm[i]->GetTuboExtremo(0).Pipe->GetPresion(0) * FDeltaT / FNumeroUnionesAdm;
300  }
301  for(int i = 0; i < FNumeroUnionesEsc; i++) {
302  FPresionMedEsc += FCCValvulaEsc[i]->GetTuboExtremo(0).Pipe->GetPresion(0) * FDeltaT / FNumeroUnionesEsc;
303  }
304  FTimeAcumAct += FDeltaT;
305  }
306 
307  /* =============================== */
308  /* BALANCE DE MASA EN EL CILINDRO */
309  /* =============================== */
310  FMasa0 = FMasa;
311 
312  if(!FCicloCerrado) {
313  // GASTO POR ADMISION
314  for(int i = 0; i < FNumeroUnionesAdm; i++) {
315  FMasaValvAdm = -dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getMassflow() *
316  FDeltaT; // signo - xq el massflow entrante se calcula negativo en la BC. Asi en el cilindro sera positivo si es entrante y viceversa.
317  if(FAnguloActual - FDeltaAngulo < FDistribucion.CA
318  && FAnguloActual > FDistribucion.CA) { // FDistribucion.CA es el CrankAngle de cierre de la admision.
319  FMasaValvAdm = FMasaValvAdm * (FDistribucion.CA - FAnguloActual + FDeltaAngulo) / FDeltaAngulo;
320  }
321  MasaAdmInstante += FMasaValvAdm;
322  FMasa += FMasaValvAdm;
323  /* SUMATORIO DE LA MASA QUE ATRAVIESA CADA UNA DE LAS VALVULAS DE ADMISION */
324  FAcumMasaPorAdm += FMasaValvAdm;
325 
326  /* Transporte de especies quimicas */
327  for(int j = 0; j < FMotor->getSpeciesNumber() - FIntEGR; j++) {
328  FMasaEspecie[j] += FCCValvulaAdm[i]->GetFraccionMasicaEspecie(j) * FMasaValvAdm;
329  }
330  if(FHayEGR)
331  FAcumMasaEGR += FCCValvulaAdm[i]->GetFraccionMasicaEspecie(FMotor->getSpeciesNumber() - 1) * FMasaValvAdm;
332  }
333  // GASTO POR ESCAPE
334  for(int i = 0; i < FNumeroUnionesEsc; i++) {
335  masavalesc = -dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getMassflow() *
336  FDeltaT; // signo - xq el massflow entrante se calcula negativo en la BC. Asi en el cilindro sera positivo si es entrante y viceversa.
337  if(FAnguloActual - FDeltaAngulo < FDistribucion.CE && FAnguloActual > FDistribucion.CE) {
338  masavalesc = masavalesc * (FDistribucion.CE - FAnguloActual + FDeltaAngulo) /
339  FDeltaAngulo; // FDistribucion.CE es el CrankAngle de cierre del escape.
340  }
341  MasaEscInstante += masavalesc;
342  FMasa += masavalesc;
343  FAcumMasaPorEsc += masavalesc;
344 
345  /* Transporte de especies quimicas */
346  for(int j = 0; j < FMotor->getSpeciesNumber() - FIntEGR; j++) {
347  FMasaEspecie[j] += FCCValvulaEsc[i]->GetFraccionMasicaEspecie(j) * masavalesc;
348  }
349  }
350  }
351 
352  /* ================================= */
353  /* MASA DE BLOW-BY */
354  /* ================================= */
355  if(FPressure > FMotor->getPresionAmb() && FMotor->getGeometria().CDBlowBy > 0.) {
356  if(FMotor->getPresionAmb() / FPressure < pow(2. / FGamma2, FGamma6)) {
357  FPresionCarter = FPressure * pow(2. / FGamma2, FGamma4 / 2.);
358 
359  } else {
360  FPresionCarter = FMotor->getPresionAmb();
361  }
362  z = FGamma4 * (pow(FPresionCarter / FPressure, 2 / FGamma) - pow(FPresionCarter / FPressure, FGamma2 / FGamma));
363  FGastoBlowBy = FMotor->getGeometria().CDBlowBy * 3.5e-5 * FMotor->getGeometria().Diametro * __units::BarToPa(
364  FPressure) * sqrt(z / FRMezcla / __units::degCToK(FTemperature));
365  } else {
366  FGastoBlowBy = 0.;
367  }
368  FMasaBlowBy = FGastoBlowBy * FDeltaT;
369 
370  /* ================================= */
371  /* INYECCION DE COMBUSTIBLE (MEC) */
372  /* ================================= */
373 
374  //Inyeccion con datos de Angulo y duracion
375  if(!FInyeccion) {
376  for(int i = 0; i < FNumIny; i++) {
377  if(FAnguloActual > FAnguloInjeccion[i] && FAnguloAnterior <= FAnguloInjeccion[i]) {
378  // En el angulo de begining de la combusti�n se empieza a introducir el combustible
379  // Se pasa a estado de injeccion verdadero.
380  // FMasaFuel = FMotor->getMasaFuel();
381  FInyeccion = true;
382  FTasaFuel = 0.;
383  ind = i;
384  }
385  }
386  }
387 
388  if(FInyeccion) {
389  if(!FCicloCerrado) {
390  FInyeccion = false;
391  FFuelTotal = 0.;
392  FFuelInstant = 0.;
393  FFuelAcum = 0.;
394  }
395  if(FMotor->getACT()) {
396  FTasaFuel = Interp1(FAnguloComb, FCAD_injection_rate, Finjection_rate, FCAI) / 1000.;
397  FFuelInstant = FTasaFuel * FDeltaT;
398 
399  if(FAnguloComb > FCAD_injection_rate[FCAI - 1]) {
400  FInyeccion = false;
401  }
402  } else {
403  // Se reparte la inyeccion de combustible durante 1.3 ms
404  // FFuelInstant representa la cantidad de masa de fuel introducida en un instante de calculo.
405  if(FMotor->getFTipoDatosIny() == 2) {
406  FTasaFuel = FMotor->TasaInyInterp(FAnguloComb);
407  } else {
408  FTasaFuel = FMasaFuel * FPercentInyeccion[ind] / FTInyeccion[ind];
409  }
410  FFuelInstant = FTasaFuel * FDeltaT;
411  // Se va acumulando la masa de fuel para comprobar que no super el valor de combustible fijado por cilindro y ciclo.
412  }
413  if((FFuelAcum + FFuelInstant) > FMasaFuel * FPercentInyeccion[ind]) {
414  // La inyeccion se corta cuando se alcanza el valor de combustible total.
415  FFuelInstant = FMasaFuel * FPercentInyeccion[ind] - FFuelAcum;
416  FInyeccion = false;
417  FFuelAcum = 0.;
418  } else {
419  FFuelAcum += FFuelInstant;
420  FFuelTotal += FFuelInstant;
421  }
422  } else {
423  FFuelInstant = 0.;
424  }
425 
426  FMasa = FMasa + FFuelInstant - FMasaBlowBy;
427 
428  // AL CIERRE DE LA VALVULA DE ADMISION SE OBTIENE LA MASA TOTAL POR ADMISION Y ESCAPE,
429  // LA MASA ATRAPADA EN EL CILINDRO, EL COMBUSTIBLE (PARA LOS MEP EN FUNCION DE
430  // LA MASA POR ADMISION) LA PRESION Y LOS ANGULOS DE INICIO Y FIN DE LA COMBUSTION
431  if(FAnguloActual > FDistribucion.CA && FAnguloAnterior <= FDistribucion.CA && !FCicloCerrado) {
432 
433  FCicloCerrado = true;
434  FPrimerInstanteCicloCerrado = true;
435  FMasaPorAdmision = FAcumMasaPorAdm;
436 
437  FAcumMasaPorAdm = 0;
438  FMasaPorEscape = FAcumMasaPorEsc;
439  FAcumMasaPorEsc = 0;
440  FMasaEGR = FAcumMasaEGR;
441  FAcumMasaEGR = 0;
442  FMasaAtrapada = FMasa;
443 
444  FPresionRCA = FPressure;
445 
446  if(FMotor->getACT()) {
447  if(FMotor->getSpeciesNumber() == 9) {
448  FSpecies_IVC[0] = FFraccionMasicaEspecie[7]; // N2
449  } else if(FMotor->getSpeciesNumber() == 10) {
450  FSpecies_IVC[0] = FFraccionMasicaEspecie[8]; // N2
451  }
452  FSpecies_IVC[1] = FFraccionMasicaEspecie[0]; // O2
453  FSpecies_IVC[2] = FFraccionMasicaEspecie[1]; // CO2
454  FSpecies_IVC[3] = FFraccionMasicaEspecie[2]; // H2O
455  FSpecies_IVC[4] = FFraccionMasicaEspecie[5]; // NOx
456  FSpecies_IVC[5] = FFraccionMasicaEspecie[6]; // CO
457  FSpecies_IVC[6] = FFraccionMasicaEspecie[4]; // SOOT
458  FSpecies_IVC[7] = FFraccionMasicaEspecie[3]; // HC
459  }
460 
461  // Transporte de Especies Quimicas. Se inician las masas de cada especie en el ciclo cerrado.
462  if(FMotor->getSpeciesModel() == nmCalculoCompleto) {
463  if(FMotor->getSpeciesNumber() == 9) {
464  for(int i = 0; i < FMotor->getSpeciesNumber() - FIntEGR; i++) {
465  FFraccionComienzoCicloCerrado[i] = FFraccionMasicaEspecie[i];
466  }
467  FComposicionCicloCerrado[1] = 0.;
468  // No llega combustible de los tubos
469  // FComposicionCicloCerrado[2]=FFraccionMasicaEspecie[0]+FFraccionMasicaEspecie[1]+FFraccionMasicaEspecie[2]+FFraccionMasicaEspecie[7]; // Aire fresco
470  FComposicionCicloCerrado[2] = FFraccionMasicaEspecie[0] / FMotor->GetComposicionAtmosfera(0);
471  FComposicionCicloCerrado[0] = 1 - FComposicionCicloCerrado[2]; // Gases Quemados
472 
473  } else if(FMotor->getSpeciesNumber() == 10) {
474  for(int i = 0; i < FMotor->getSpeciesNumber() - FIntEGR; i++) {
475  FFraccionComienzoCicloCerrado[i] = FFraccionMasicaEspecie[i];
476  }
477  FComposicionCicloCerrado[1] = FFraccionMasicaEspecie[7];
478  // Combustible
479  // FComposicionCicloCerrado[2]=FFraccionMasicaEspecie[0]+FFraccionMasicaEspecie[1]+FFraccionMasicaEspecie[2]+FFraccionMasicaEspecie[8]; // Aire fresco
480  FComposicionCicloCerrado[2] = FFraccionMasicaEspecie[0] / FMotor->GetComposicionAtmosfera(0);
481  FComposicionCicloCerrado[0] = 1 - FComposicionCicloCerrado[1] - FComposicionCicloCerrado[2]; // Gases Quemados
482  }
483  } else if(FMotor->getSpeciesModel() == nmCalculoSimple) {
484  if(FMotor->getSpeciesNumber() == 3) {
485  FComposicionCicloCerrado[1] = 0.;
486  // No llega combustible de los tubos
487  FComposicionCicloCerrado[2] = FFraccionMasicaEspecie[1];
488  // Aire fresco
489  FComposicionCicloCerrado[0] = FFraccionMasicaEspecie[0];
490  // Gases Quemados
491 
492  } else if(FMotor->getSpeciesNumber() == 4) {
493  FComposicionCicloCerrado[1] = FFraccionMasicaEspecie[1];
494  // Combustible
495  FComposicionCicloCerrado[2] = FFraccionMasicaEspecie[2];
496  // Aire fresco
497  FComposicionCicloCerrado[0] = FFraccionMasicaEspecie[0];
498  // Gases Quemados
499  }
500  }
501  for(int j = 0; j < 3; j++) {
502  FMasaEspecieCicloCerrado[j] = FMasa * FComposicionCicloCerrado[j];
503  }
504 
505  if(FHayEGR)
506  FFraccionMasicaEspecie[FMotor->getSpeciesNumber() - 1] = 1.;
507 
508  if(FMotor->getCombustible() == nmMEP) {
509  CalculaFuelMEP(FMasaEspecieCicloCerrado[2]);
510  }
511 
512  std::cout << "INFO: End of Gas-exchange process in cylinder " << FNumeroCilindro << std::endl;
513  std::cout << "____________________________________________________" << std::endl;
514  std::cout << "INFO: Intake mass: " << FMasaPorAdmision * 1e3 << " (g)" << std::endl;
515  std::cout << "INFO: Exhaust mass: " << FMasaPorEscape * 1e3 << " (g)" << std::endl;
516  std::cout << "INFO: Trapped mass: " << FMasaAtrapada * 1e3 << " (g)" << std::endl;
517  // std::cout << "INFO: Masa de combustible: " << FMasaFuel*1e6 << " (mg)" << std::endl;
518  std::cout << "INFO: Pressure at I.C.: " << FPresionRCA << " (bar)" << std::endl;
519  // std::cout << "____________________________________________________" << std::endl;
520  }
521 
522  /* ================================= */
523  /* BALANCE DE ENERGIA EN EL CILINDRO */
524  /* ================================= */
525 
526  FVolumen0 = FVolumen;
527  FVolumen = CalculaVolumen(FAnguloActual);
528 
529  // CALCULO DE LA TRANSMISION DE CALOR
530 
531  if(FAnguloActual >= 0. && FAnguloActual <= 180.) {
532  Fequis = Fratioctm + 1 / (pow(cosh(FAnguloActual / 100.), 40.) + Fratioctm / (1 - Fratioctm));
533  } else if(FAnguloActual >= 540. && FAnguloActual <= 720.) {
534  Fequis = Fratioctm + 1 / (pow(cosh((720. - FAnguloActual) / 100.), 40.) + Fratioctm / (1 - Fratioctm));
535  } else {
536  Fequis = 0.;
537  }
538  FCu = __units::RPMToRad_s(FMotor->getRegimen()) * Fequis * Fctorbadmp * FKctm;
539  if(FMotor->getGeometria().DiametroBowl > 0.0) {
540  FCu *= pow2(FMotor->getGeometria().Diametro) / (2. * FMotor->getGeometria().DiametroBowl);
541  }
542 
543  // MOTOR SIN COMBUSTION O EN ARRASTRE
544  if(FAnguloComb < FIniComb || FAnguloComb > FFinComb || FMasaFuel == 0.) {
545  Fc2 = 0;
546 
547  // MOTOR EN COMBUSTION
548  } else {
549  Fc2 = 0.001; // antes 3.24e-3 (ahora 0.001) MOTOR DE CAMARA ABIERTA ; 6.22e-3 MOTOR DE C�MARA DIVIDIDA
550  }
551 
552  // COEFICIENTE DE PELICULA DE WOSCHNI
553  double deltaP = __units::BarToPa(FPressure - FPresionRCA * pow(FVolumenCA / FVolumen, 1.36));
554  if(deltaP < 0.)
555  deltaP = 0.;
556 
557  FCm = CalculaCm();
558 
559  if(FAnguloActual >= 180. && FAnguloActual <= 360.) {
560  // CARRERA DE ESCAPE
561  Fh = 1.2e-2 * pow(__units::BarToPa(FPressure), 0.8) * pow(__units::degCToK(FTemperature),
562  -0.53) * pow(FMotor->getGeometria().Diametro, -0.2) * pow(
563  (6.18 * FCm + 0.417 * FCu) * FMotor->getAjusteTranCalEsc() + Fc2 * deltaP * FMotor->getGeometria().CilindradaUnitaria /
564  FMasaAtrapada / FRMezcla, 0.8);
565 
566  } else if(FAnguloActual >= 360. && FAnguloActual <= 540.) {
567  // CARRERA DE ADMISION
568  Fh = 1.2e-2 * pow(__units::BarToPa(FPressure), 0.8) * pow(__units::degCToK(FTemperature),
569  -0.53) * pow(FMotor->getGeometria().Diametro, -0.2) * pow(
570  (6.18 * FCm + 0.417 * FCu) * FMotor->getAjusteTranCalAdm() + Fc2 * deltaP * FMotor->getGeometria().CilindradaUnitaria /
571  FMasaAtrapada / 287, 0.8);
572  } else { // CARRERAS DE COMPRESION Y EXPANSION
573  Fh = 1.2e-2 * pow(__units::BarToPa(FPressure), 0.8) * pow(__units::degCToK(FTemperature),
574  -0.53) * pow(FMotor->getGeometria().Diametro, -0.2) * pow(
575  (FMotor->getWoschni().cw1 * FCm + FMotor->getWoschni().cw2 * FCu) + Fc2 * deltaP *
576  FMotor->getGeometria().CilindradaUnitaria / FMasaAtrapada / 287, 0.8);
577  }
578 
579  if(FMotor->getCalculoPared() != nmTempFija)
580  CalculaTemperaturasPared();
581 
582  // CALOR TRANSMITIDO A LAS PAREDES
583  FCalor.TransCilindro = Fh * (FTempPared[0].Cylinder - FTemperature) * 4 * (FVolumen / 2. + FVolumen0 / 2. -
584  FMotor->getGeometria().VCC) / FMotor->getGeometria().Diametro;
585  FCalor.TransPiston = Fh * (FTempPared[0].Piston - FTemperature) * FMotor->getGeometria().AreaPiston;
586  FCalor.TransCulata = Fh * (FTempPared[0].Culata - FTemperature) * FMotor->getGeometria().AreaCulata;
587  FCalor.TransTotal = (FCalor.TransCilindro + FCalor.TransPiston + FCalor.TransCulata) * FDeltaT;
588 
589  // CALCULO DEL CALOR LIBERADO
590  if(FCicloCerrado) {
591  if(FAnguloComb < FIniComb || FAnguloComb > FFinComb || FMasaFuel == 0.) {
592  FCalor.Liberado = 0.;
593  FCalor.LiberadoTotal = 0.;
594  FCalor.FQL = 0.;
595  FCalor.FQL0 = 0.;
596 
597  FMasaEspecieCicloCerrado[0] = FMasaEspecieCicloCerrado[0] - FComposicionCicloCerrado[0] * FMasaBlowBy;
598  // Gases Quemados
599  FMasaEspecieCicloCerrado[1] = FMasaEspecieCicloCerrado[1] + FFuelInstant;
600  // Combustible
601  FMasaEspecieCicloCerrado[2] = FMasa - FMasaEspecieCicloCerrado[0] - FMasaEspecieCicloCerrado[1];
602  // Aire Fresco
603 
604  FComposicionCicloCerrado[0] = FMasaEspecieCicloCerrado[0] / FMasa;
605  FComposicionCicloCerrado[1] = FMasaEspecieCicloCerrado[1] / FMasa;
606  FComposicionCicloCerrado[2] = FMasaEspecieCicloCerrado[2] / FMasa;
607 
608  } else {
609 
610  if(FPrimeraCombustion) {
611  if(FCalcComb == nmFQL) {
612  Fecg0 = 0.;
613  FecgTotal = 0.;
614  //(FHcl - FUfgasoil) * (FFuelInstant + FFuelInstantPil) / FDeltaAngulo;
615  //energía del combustible gaseoso (evaporación), hay que integrala
616  FCalor.FQL0 = CalculaCalorLiberado(FAnguloComb - FDeltaAngulo);
617  } else {
618  FCalor.FQL0 = FCalor.FQL;
619  }
620  FPrimeraCombustion = false;
621  }
622  if(FCalcComb == nmFQL) {
623  FUfgasoil = CalculoUfgasoil(FTemperature);
624  Fecg = (FHcl - FUfgasoil) * FFuelInstant;
625  //energía del combustible gaseoso (evaporación), hay que integrarla
626  FCalor.FQL = CalculaCalorLiberado(FAnguloComb);
627  FecgInt = (Fecg + Fecg0) / 2 * FDeltaT;
628  FecgTotal += (Fecg + Fecg0) / 2 * FDeltaT;
629  Fecg0 = Fecg;
630  } else if(FCalcComb == nmACT) {
631  FCalor.FQL = Interp1(FAnguloComb, FCAD_exit, FHRF_exit, FCAI);
632  // printf("%lf %lf\n",FAnguloComb,FCalor.FQL);
633  }
634  FCalor.Liberado = (FCalor.FQL - FCalor.FQL0) * FMfint * FMotor->getPoderCalorifico() *
635  FMotor->getRendimientoCombustion();
636  FCalor.LiberadoTotal += FCalor.Liberado;
637  if(FCalor.Liberado < 0)
638  FCalor.Liberado = 0.;
639 
640  // Transporte de Especies Quimicas
641  FMfquem = FCalor.Liberado / FMotor->getPoderCalorifico();
642  FMairequem = FMfquem / FDosadoEstequiometrico;
643  if(FMairequem > (FMasaEspecieCicloCerrado[2])) {
644  // Se necesita mas aire del que tenemos para completar la combustion.
645  mfquefin = FMasaEspecieCicloCerrado[2] * FDosadoEstequiometrico;
646  if(!FSaturado) {
647  printf("WARNING: Cylinder %d does not have enough air at CrankAngle %lf\n", FNumeroCilindro, FAnguloActual);
648  FCalor.Liberado = (FCalor.FQL - FCalor.FQL0) * mfquefin * FMotor->getPoderCalorifico() *
649  FMotor->getRendimientoCombustion();
650 
651  FSaturado = true;
652 
653  FMasaEspecieCicloCerrado[0] = FMasaEspecieCicloCerrado[0] + FMairequem + mfquefin - FComposicionCicloCerrado[0] *
654  FMasaBlowBy;
655  // Gases Quemados
656  FMasaEspecieCicloCerrado[1] = FMasaEspecieCicloCerrado[1] - mfquefin - FComposicionCicloCerrado[1] * FMasaBlowBy +
657  FFuelInstant;
658  // Combustible
659  FMasaEspecieCicloCerrado[2] = 0.; // Aire Fresco
660 
661  } else {
662  FCalor.Liberado = 0.;
663  }
664 
665  } else {
666  FSaturado = false;
667 
668  FMasaEspecieCicloCerrado[0] = FMasaEspecieCicloCerrado[0] + FMairequem + FMfquem - FComposicionCicloCerrado[0] *
669  FMasaBlowBy;
670  // Gases Quemados
671  FMasaEspecieCicloCerrado[1] = FMasaEspecieCicloCerrado[1] - FMfquem - FComposicionCicloCerrado[1] * FMasaBlowBy +
672  FFuelInstant;
673 // if (FMasaEspecieCicloCerrado[1] < 0) {
674 // FMasaEspecieCicloCerrado[1] = 0;
675 // }
676  // Combustible
677  FMasaEspecieCicloCerrado[2] = FMasaEspecieCicloCerrado[2] - FMairequem - FComposicionCicloCerrado[2] * FMasaBlowBy;
678  // Aire Fresco
679  }
680 
681  FComposicionCicloCerrado[0] = FMasaEspecieCicloCerrado[0] / FMasa;
682  FComposicionCicloCerrado[1] = FMasaEspecieCicloCerrado[1] / FMasa;
683  FComposicionCicloCerrado[2] = FMasaEspecieCicloCerrado[2] / FMasa;
684 
685  FCalor.FQL0 = FCalor.FQL;
686  }
687  }
688 
689  /* Transporte de Especies Quimicas */
690  if(FCicloCerrado) {
691  if(FMotor->getSpeciesModel() == nmCalculoCompleto) {
692  if(FMotor->getACT()) {
693  double Yfuel = (FFuelTotal - FMfint * FCalor.FQL) / FMasa;
694  if(Yfuel < 0)
695  Yfuel = 0;
696 
697  FFraccionMasicaEspecie[0] = (FSpecies_IVC[1] * (1 - FCalor.FQL) + dataOUT.species_EVO[1] * FCalor.FQL) * (1 - Yfuel);
698  FFraccionMasicaEspecie[1] = (FSpecies_IVC[2] * (1 - FCalor.FQL) + dataOUT.species_EVO[2] * FCalor.FQL) * (1 - Yfuel);
699  FFraccionMasicaEspecie[2] = (FSpecies_IVC[3] * (1 - FCalor.FQL) + dataOUT.species_EVO[3] * FCalor.FQL) * (1 - Yfuel);
700  FFraccionMasicaEspecie[3] = (FSpecies_IVC[7] * (1 - FCalor.FQL) + dataOUT.species_EVO[7] * FCalor.FQL) * (1 - Yfuel);
701  FFraccionMasicaEspecie[4] = (FSpecies_IVC[6] * (1 - FCalor.FQL) + dataOUT.species_EVO[6] * FCalor.FQL) * (1 - Yfuel);
702  FFraccionMasicaEspecie[5] = (FSpecies_IVC[4] * (1 - FCalor.FQL) + dataOUT.species_EVO[4] * FCalor.FQL) * (1 - Yfuel);
703  FFraccionMasicaEspecie[6] = (FSpecies_IVC[5] * (1 - FCalor.FQL) + dataOUT.species_EVO[5] * FCalor.FQL) * (1 - Yfuel);
704 
705  if(FMotor->getSpeciesNumber() == 9) {
706  FFraccionMasicaEspecie[3] += Yfuel;
707  FFraccionMasicaEspecie[7] = (FSpecies_IVC[0] * (1 - FCalor.FQL) + dataOUT.species_EVO[0] * FCalor.FQL) * (1 - Yfuel);
708 
709  } else if(FMotor->getSpeciesNumber() == 10) {
710  FFraccionMasicaEspecie[7] = Yfuel;
711  FFraccionMasicaEspecie[8] = (FSpecies_IVC[0] * (1 - FCalor.FQL) + dataOUT.species_EVO[0] * FCalor.FQL) * (1 - Yfuel);
712  }
713  for(int j = 0; j < FMotor->getSpeciesNumber() - FIntEGR; j++) {
714  FMasaEspecie[j] = FFraccionMasicaEspecie[j] * FMasa;
715  }
716  } else {
717  if(FMotor->getSpeciesNumber() == 9) {
718  if(!FSaturado) {
719  FMolesCombQuemado = FMfquem / (FXComb * 12.01 + FYComb * 1.01 + FZComb * 16);
720  // Moles de combustible quemado en el incremento temporal calculado
721  } else {
722  FMolesCombQuemado = mfquefin / (FXComb * 12.01 + FYComb * 1.01 + FZComb * 16);
723  // Moles de combustible quemado en el incremento temporal calculado
724  }
725  FMasaO2Reactivos = FMolesCombQuemado * (FYComb / 4 + FXComb - FZComb / 2) * __PM::O2;
726  FMasaN2Reactivos = FMolesCombQuemado * (FYComb / 4 + FXComb - FZComb / 2) * FRelacionMolarN2_O2 * __PM::N2;
727  // FMasaH2OReactivos=FMolesCombQuemado*(FYComb/4+FXComb-FZComb/2)*FRelacionMolarH2O_O2*__PM::H2O;
728  FMasaCO2Productos = FMolesCombQuemado * FXComb * __PM::CO2;
729  FMasaH2OProductos = FMolesCombQuemado * (FYComb / 2
730  /* +(FYComb/4+FXComb-FZComb/2)*FRelacionMolarH2O_O2 */) * __PM::H2O;
731  FMasaN2Productos = FMolesCombQuemado * (FYComb / 4 + FXComb - FZComb / 2) * FRelacionMolarN2_O2 * __PM::N2;
732 
733  FFraccionMasicaEspecie[0] = (FMasaEspecie[0] - FMasaO2Reactivos - FFraccionMasicaEspecie[0] * FMasaBlowBy) / FMasa;
734  // Y O2
735  FFraccionMasicaEspecie[1] = (FMasaEspecie[1] + FMasaCO2Productos - FFraccionMasicaEspecie[1] * FMasaBlowBy) / FMasa;
736  // Y CO2
737  FFraccionMasicaEspecie[2] = (FMasaEspecie[2] + FMasaH2OProductos - FFraccionMasicaEspecie[2] *
738  FMasaBlowBy /* -FMasaH2OReactivos */) / FMasa;
739  // Y H2O
740  FFraccionMasicaEspecie[7] = (FMasaEspecie[7] - FFraccionMasicaEspecie[7] * FMasaBlowBy) / FMasa;
741  // Y N2
742  FFraccionMasicaEspecie[3] = FComposicionCicloCerrado[1];
743  // Combustible que queda sin quemar
744  FFraccionMasicaEspecie[4] = 0;
745  FFraccionMasicaEspecie[5] = 0;
746  FFraccionMasicaEspecie[6] = 0;
747  for(int j = 0; j < FMotor->getSpeciesNumber() - FIntEGR; j++) {
748  FMasaEspecie[j] = FFraccionMasicaEspecie[j] * FMasa;
749  }
750 
751  } else if(FMotor->getSpeciesNumber() == 10) {
752  if(!FSaturado) {
753  FMolesCombQuemado = FMfquem / (FXComb * 12.01 + FYComb * 1.01 + FZComb * 16);
754  // Moles de combustible quemado en el incremento temporal calculado
755  } else {
756  FMolesCombQuemado = mfquefin / (FXComb * 12.01 + FYComb * 1.01 + FZComb * 16);
757  // Moles de combustible quemado en el incremento temporal calculado
758  }
759  FMasaO2Reactivos = FMolesCombQuemado * (FYComb / 4 + FXComb - FZComb / 2) * __PM::O2;
760  FMasaN2Reactivos = FMolesCombQuemado * (FYComb / 4 + FXComb - FZComb / 2) * FRelacionMolarN2_O2 * __PM::N2;
761  // FMasaH2OReactivos=FMolesCombQuemado*(FYComb/4+FXComb-FZComb/2)*FRelacionMolarH2O_O2*__PM::H2O;
762  FMasaCO2Productos = FMolesCombQuemado * FXComb * __PM::CO2;
763  FMasaH2OProductos = FMolesCombQuemado * (FYComb / 2
764  /* +(FYComb/4+FXComb-FZComb/2)*FRelacionMolarH2O_O2 */) * __PM::H2O;
765  FMasaN2Productos = FMolesCombQuemado * (FYComb / 4 + FXComb - FZComb / 2) * FRelacionMolarN2_O2 * __PM::N2;
766 
767  FFraccionMasicaEspecie[0] = (FMasaEspecie[0] - FMasaO2Reactivos - FFraccionMasicaEspecie[0] * FMasaBlowBy) / FMasa;
768  // Y O2
769  FFraccionMasicaEspecie[1] = (FMasaEspecie[1] + FMasaCO2Productos - FFraccionMasicaEspecie[1] * FMasaBlowBy) / FMasa;
770  // Y CO2
771  FFraccionMasicaEspecie[2] = (FMasaEspecie[2] + FMasaH2OProductos - FFraccionMasicaEspecie[2] *
772  FMasaBlowBy /* -FMasaH2OReactivos */) / FMasa;
773  // Y H2O
774  FFraccionMasicaEspecie[8] = (FMasaEspecie[8] - FFraccionMasicaEspecie[8] * FMasaBlowBy) / FMasa;
775  // Y N2
776  FFraccionMasicaEspecie[3] = 0;
777  FFraccionMasicaEspecie[4] = 0;
778  FFraccionMasicaEspecie[5] = 0;
779  FFraccionMasicaEspecie[6] = 0;
780  FFraccionMasicaEspecie[7] = FComposicionCicloCerrado[1];
781  for(int j = 0; j < FMotor->getSpeciesNumber() - FIntEGR; j++) {
782  FMasaEspecie[j] = FFraccionMasicaEspecie[j] * FMasa;
783  }
784 
785  }
786  }
787  } else if(FMotor->getSpeciesModel() == nmCalculoSimple) {
788  if(FMotor->getSpeciesNumber() == 3) {
789  FFraccionMasicaEspecie[0] = FComposicionCicloCerrado[0] + FComposicionCicloCerrado[1];
790  FFraccionMasicaEspecie[1] = FComposicionCicloCerrado[2];
791 
792  for(int j = 0; j < FMotor->getSpeciesNumber() - FIntEGR; j++) {
793  FMasaEspecie[j] = FFraccionMasicaEspecie[j] * FMasa;
794  }
795  } else if(FMotor->getSpeciesNumber() == 4) {
796  FFraccionMasicaEspecie[0] = FComposicionCicloCerrado[0];
797  FFraccionMasicaEspecie[1] = FComposicionCicloCerrado[1];
798  FFraccionMasicaEspecie[2] = FComposicionCicloCerrado[2];
799 
800  for(int j = 0; j < FMotor->getSpeciesNumber() - FIntEGR; j++) {
801  FMasaEspecie[j] = FFraccionMasicaEspecie[j] * FMasa;
802  }
803  }
804  }
805  } else if(!FCicloCerrado) {
806  for(int j = 0; j < FMotor->getSpeciesNumber() - 2; j++) {
807  FFraccionMasicaEspecie[j] = (FMasaEspecie[j] - FFraccionMasicaEspecie[j] * FMasaBlowBy) / FMasa;
808  FraccionMasicaAcum += FFraccionMasicaEspecie[j];
809  FMasaEspecie[j] = FFraccionMasicaEspecie[j] * FMasa;
810  }
811  FFraccionMasicaEspecie[FMotor->getSpeciesNumber() - 2] = 1. - FraccionMasicaAcum;
812  FMasaEspecie[FMotor->getSpeciesNumber() - 2] = FFraccionMasicaEspecie[FMotor->getSpeciesNumber() - 2] * FMasa;
813  if(FHayEGR) {
814  FFraccionMasicaEspecie[FMotor->getSpeciesNumber() - 1] = (FMasaEspecie[FMotor->getSpeciesNumber() - 1] -
815  FFraccionMasicaEspecie[FMotor->getSpeciesNumber() - 1] * FMasaBlowBy) / FMasa;
816  FMasaEspecie[FMotor->getSpeciesNumber() - 1] = FFraccionMasicaEspecie[FMotor->getSpeciesNumber() - 1] * FMasa;
817  }
818  }
819  /* FIN Transporte de Especies Quimicas */
820 
821  if(FAnguloActual > FDistribucion.CA && FAnguloAnterior <= FDistribucion.CA) {
822  if(FMasaFuel == 0.) {
823  FAFR = 0.;
824  } else {
825  FAFR = FMasaEspecieCicloCerrado[2] / FMasaFuel;
826  // Masa Aire Fresco Inicial al cierre/Masa Fuel
827  }
828  }
829 
830 // if (FMotor->getGammaCalculation() == nmGammaConstante) {
831 // if (FMotor->getSpeciesModel() == nmCalculoCompleto) {
832 // if (!FCicloCerrado) {
834 // FRMezcla = CalculoSimpleRMezcla(0.1, nmComposicionTemperatura)
835 // FCvMezcla = CalculoSimpleCvMezcla(FTemperature + 273., 0.1, nmComposicionTemperatura);
836 //
837 // }
838 // else if (FCicloCerrado) {
841 // FRMezcla = CalculoSimpleRMezcla(FComposicionCicloCerrado[0], nmComposicionTemperatura);
842 // FCvMezcla = CalculoSimpleCvMezcla(FTemperature + 273., FComposicionCicloCerrado[0], nmComposicionTemperatura);
843 // }
844 // }
845 // else if (FMotor->getSpeciesModel() == nmCalculoSimple) {
848 // FRMezcla = CalculoSimpleRMezcla(FComposicionCicloCerrado[0], nmComposicionTemperatura);
849 // FCvMezcla = CalculoSimpleCvMezcla(FTemperature + 273., FComposicionCicloCerrado[0], nmComposicionTemperatura);
850 //
851 // }
852 // FGamma = CalculoSimpleGamma(FRMezcla, FCvMezcla, nmComposicionTemperatura);
853 // FGamma1 = __Gamma::G1(FGamma);
854 // FGamma2 = __Gamma::G2(FGamma);
855 // FGamma4 = __Gamma::G4(FGamma);
856 // FGamma6 = __Gamma::G6(FGamma);
857 // }
858  bool PrimerPaso = true;
859  double ASon0 = FAsonido;
860  double ASon1 = FAsonido;
861  double Temp0 = FTemperature;
862  double Temp1 = FTemperature;
863  double MasTemp0 = 1 / __units::degCToK(FTemperature) / FMasa0;
864  double MasTemp1 = 1 / __units::degCToK(FTemperature) / FMasa0;
865  double MasTempMed = 0.;
866  bool CotaError = false;
867  double H1 = 0.;
868  double H0 = 0.;
869  double Energia = 0.;
870  double Error = 0.;
871  double Diff = 0.;
872 
873  // ITERACION PARA BUSCAR EL ESTADO TERMODINAMICO FINAL
874  while(!CotaError) {
875  // ENTALPIA POR ADMISION;
876  H1 = 0.;
877  for(int i = 0; i < FNumeroUnionesAdm; i++) {
878  if(dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getMassflow() != 0
879  && dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getSentidoFlujo() == nmEntrante) {
880  if(PrimerPaso) {
881  H1 += EntalpiaEntrada(dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getSpeedsound(),
882  dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getVelocity(),
883  -dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getMassflow() * FDeltaT, ASon1 / __cons::ARef, FMasa0);
884  } else {
885  H1 += EntalpiaEntrada(dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getSpeedsound(),
886  dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getVelocity(),
887  -dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getMassflow() * FDeltaT, ASon1 / __cons::ARef, FMasa);
888  }
889  }
890  }
891  // ENTALPIA POR ESCAPE;
892  for(int i = 0; i < FNumeroUnionesEsc; i++) {
893  if(dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getMassflow() != 0
894  && dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getSentidoFlujo() == nmEntrante) {
895  if(PrimerPaso) {
896  H1 += EntalpiaEntrada(dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getSpeedsound(),
897  dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getVelocity(),
898  -dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getMassflow() * FDeltaT, ASon1 / __cons::ARef, FMasa0);
899  } else {
900  H1 += EntalpiaEntrada(dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getSpeedsound(),
901  dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getVelocity(),
902  -dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getMassflow() * FDeltaT, ASon1 / __cons::ARef, FMasa);
903  }
904  }
905  }
906  if(PrimerPaso) {
907  H0 = H1;
908  PrimerPaso = false;
909  }
910  MasTempMed = (MasTemp0 + MasTemp1) / 2;
911  Energia = FVolumen0 * FMasa / FVolumen / FMasa0 * exp((H1 + H0) / 2. + (FCalor.TransTotal + FCalor.Liberado + FecgInt)
912  * (MasTempMed / FRMezcla));
913  Temp1 = __units::KTodegC(__units::degCToK(FTemperature) * pow(Energia, FGamma1));
914  Error = (Diff = Temp1 - Temp0, fabs(Diff)) / Temp1;
915  if(Error > 1e-6) {
916  MasTemp1 = 1. / (__units::degCToK(Temp1) * FMasa);
917  Temp0 = Temp1;
918  } else {
919  CotaError = true;
920  }
921  }
922  // ACTUALIZA NUEVAS PROPIEDADES CILINDRO
923 
924  // if (FNumeroCilindro == 1) {
925  // std::cout << "INFO: Imposed pressure at E.O.: " << FPressure << " (bar)" << std::endl;
926  // std::cout << "INFO: Imposed temperature at E.O.: " << FTemperature << " (\260C)" << std::endl;
927  // }
928 
929  FTemperatura0 = FTemperature;
930  FTemperature = Temp1;
931 
932  FAsonido0 = FAsonido;
933  FAsonido = sqrt(FGamma * FRMezcla * __units::degCToK(FTemperature));
934 
935  FPresion0 = FPressure;
936  FPressure = __units::PaToBar(__units::degCToK(FTemperature) * FMasa * FRMezcla / FVolumen);
937  if(FAnguloActual > FDistribucion.AE && FAnguloAnterior <= FDistribucion.AE) {
938  std::cout << "INFO: Begin gas-exchange process in cylinder " << FNumeroCilindro << std::endl;
939  std::cout << "____________________________________________________" << std::endl << std::endl;
940  std::cout << "INFO: Calculated pressure at E.O.: " << FPressure << " (bar)" << std::endl;
941  std::cout << "INFO: Calculated temperature at E.O.: " << FTemperature << " (\260C)" << std::endl;
942  if(FMotor->getCalculoDePAAE() == nmPAAEImpuesta) {
943  FPressure = FMotor->getPresionAAE();
944  FTemperature = __units::KTodegC(__units::BarToPa(FPressure) * FVolumen / FMasa / FRMezcla);
945  FAsonido = sqrt(FRMezcla * FGamma * __units::degCToK(FTemperature));
946  std::cout << "INFO: Imposed pressure at E.O.: " << FPressure << " (bar)" << std::endl;
947  std::cout << "INFO: Imposed temperature at E.O.: " << FTemperature << " (\260C)" << std::endl;
948  }
949  if(FMotor->getImponerComposicionAE()) {
950  for(int i = 0; i < FMotor->getSpeciesNumber() - 1; i++) {
951  FFraccionMasicaEspecie[i] = FMotor->GetComposicionInicial(i);
952  }
953  }
954  std::cout << "____________________________________________________" << std::endl << std::endl;
955  }
956  /* if(FNumeroCilindro==3){
957  printf(" %lf %lf %lf \n ", FGamma,FRMezcla,FFraccionMasicaEspecie[0]);
958  } */
959 
960  // CALCULA EL PAR INSTANTANEO PRODUCIDO POR EL CILINDRO
961  double a = __units::DegToRad(FAnguloActual);
962  double L = FMotor->getGeometria().Biela;
963  double R = FMotor->getGeometria().Carrera / 2.;
964  double e = FMotor->getGeometria().Excentricidad / 1000;
965  double area = __geom::Circle_area(FMotor->getGeometria().Diametro);
966 
967  double b = asin((e - R * sin(a)) / L);
968 
969  // CALCULO TRABAJO Y PAR
970  FParInstantaneo = __units::BarToPa(FPressure - FMotor->getPresionAmb()) * area * R * sin(a - b) / cos(b);
971 
972  FPreMed = (FPressure + FPresion0) / 2.;
973 
974  FTrabajoNetoACUM += __units::BarToPa(FPreMed) * (FVolumen - FVolumen0);
975  if(FAnguloActual > 180 && FAnguloActual <= 540) {
976  FTrabajoBombeoACUM += __units::BarToPa(FPreMed) * (FVolumen - FVolumen0);
977  }
978  if(FAnguloActual < FAnguloAnterior) {
979  FTrabajoNeto = FTrabajoNetoACUM;
980  FTrabajoNetoACUM = 0.;
981  FPMN = __units::PaToBar(FTrabajoNeto / FMotor->getGeometria().CilindradaUnitaria);
982  FTrabajoBombeo = FTrabajoBombeoACUM;
983  FTrabajoBombeoACUM = 0.;
984  FPMB = __units::PaToBar(FTrabajoBombeo / FMotor->getGeometria().CilindradaUnitaria);
985  FPMI = FPMN - FPMB;
986  }
987 
988  FDensidadReferencia = FMotor->getTuboRendVol()->GetGamma(FMotor->getNodoMedio()) * __units::BarToPa(
989  FMotor->getTuboRendVol()->GetPresion(FMotor->getNodoMedio())) / pow2(
990  FMotor->getTuboRendVol()->GetAsonido(FMotor->getNodoMedio()) * __cons::ARef);
991 
992  /* ================================= */
993  /* MODELO DE CORTOCIRCUITO */
994  /* ================================= */
995 
996  if(FAnguloActual > FDistribucion.AA && FAnguloActual < FDistribucion.CE) {
997  if(MasaAdmInstante > 0. && MasaEscInstante < 0.) {
998  // Cortocircuito con sentido Admision hacia Escape. (massflow positivo)
999  MasaCortocircuitoAdm = MasaAdmInstante * FAlphaEscape / __cons::Pi;
1000  MasaCortocircuitoEsc = MasaEscInstante * (__cons::Pi_2 - FAlphaEscape) / __cons::Pi;
1001  if(fabs(MasaCortocircuitoAdm) < fabs(MasaCortocircuitoEsc)) {
1002  FMasaCortocircuito = fabs(MasaCortocircuitoAdm);
1003  } else {
1004  FMasaCortocircuito = fabs(MasaCortocircuitoEsc);
1005  }
1006  FGastoCortocircuito = FMasaCortocircuito / FDeltaT;
1007  for(int j = 0; j < FMotor->getSpeciesNumber() - 2; j++) {
1008  for(int i = 0; i < FNumeroUnionesAdm; i++) {
1009  if(dynamic_cast<TCCCilindro*>(FCCValvulaAdm[i])->getSentidoFlujo() == nmEntrante) {
1010  FraccionCC += FCCValvulaAdm[i]->GetFraccionMasicaEspecie(j);
1011  NumeroUnionesEntrante++;
1012  }
1013  }
1014  FraccionCC = FraccionCC / NumeroUnionesEntrante;
1015  // MasaEscInstante tiene signo - cuando el flujo sale por el escape.
1016  FComposicionSaliente[j] = FFraccionMasicaEspecie[j] * (MasaEscInstante + FMasaCortocircuito) / MasaEscInstante -
1017  FraccionCC * FMasaCortocircuito / MasaEscInstante;
1018  FraccionCC = 0.;
1019  NumeroUnionesEntrante = 0;
1020  }
1021  if(FHayEGR)
1022  FComposicionSaliente[FMotor->getSpeciesNumber() - 1] = FFraccionMasicaEspecie[FMotor->getSpeciesNumber() - 1];
1023  } else if(MasaAdmInstante < 0. && MasaEscInstante > 0.) {
1024  // Cortocircuito con sentido Escape hacia Admision. (massflow negativo)
1025  MasaCortocircuitoAdm = MasaAdmInstante * (__cons::Pi_2 - FAlphaAdmision) / __cons::Pi;
1026  MasaCortocircuitoEsc = MasaEscInstante * FAlphaAdmision / __cons::Pi;
1027  if(fabs(MasaCortocircuitoAdm) < fabs(MasaCortocircuitoEsc)) {
1028  FMasaCortocircuito = -fabs(MasaCortocircuitoAdm);
1029  } else {
1030  FMasaCortocircuito = -fabs(MasaCortocircuitoEsc);
1031  }
1032  FGastoCortocircuito = FMasaCortocircuito / FDeltaT;
1033  for(int j = 0; j < FMotor->getSpeciesNumber() - 2; j++) {
1034  for(int i = 0; i < FNumeroUnionesEsc; i++) {
1035  if(dynamic_cast<TCCCilindro*>(FCCValvulaEsc[i])->getSentidoFlujo() == nmEntrante) {
1036  FraccionCC += FCCValvulaEsc[i]->GetFraccionMasicaEspecie(j);
1037  NumeroUnionesEntrante++;
1038  }
1039  }
1040  FraccionCC = FraccionCC / NumeroUnionesEntrante;
1041  // La masa por la valvula de admision sera negativa, por ser saliente del cilindro.
1042  // La masa de cortocircuito de escape a admision es negativa.
1043  FComposicionSaliente[j] = FFraccionMasicaEspecie[j] * (MasaAdmInstante - FMasaCortocircuito) / MasaAdmInstante +
1044  FraccionCC * FMasaCortocircuito / MasaAdmInstante;
1045  FraccionCC = 0.;
1046  NumeroUnionesEntrante = 0;
1047  }
1048  if(FHayEGR)
1049  FComposicionSaliente[FMotor->getSpeciesNumber() - 1] = FFraccionMasicaEspecie[FMotor->getSpeciesNumber() - 1];
1050  } else {
1051  // No hay cortocircuito, pues en este instante entra o sale massflow por todas las valvulas.
1052  FMasaCortocircuito = 0.;
1053  FGastoCortocircuito = 0.;
1054  for(int j = 0; j < FMotor->getSpeciesNumber() - FIntEGR; j++) {
1055  FComposicionSaliente[j] = FFraccionMasicaEspecie[j];
1056  }
1057  }
1058  } else {
1059  // No hay cortocircuito, pues alguna de las valvulas no esta abierta.
1060  FMasaCortocircuito = 0.;
1061  FGastoCortocircuito = 0.;
1062  for(int j = 0; j < FMotor->getSpeciesNumber() - FIntEGR; j++) {
1063  FComposicionSaliente[j] = FFraccionMasicaEspecie[j];
1064  }
1065  }
1066 
1067  } catch(exception & N) {
1068  std::cout << "ERROR: TCilindro4T::ActualizaPropiedades en cilindro: " << FNumeroCilindro << std::endl;
1069  std::cout << "Tipo de error: " << N.what() << std::endl;
1070  throw Exception(N.what());
1071  }
1072 }
1073 
1074 // ---------------------------------------------------------------------------
1075 // ---------------------------------------------------------------------------
1076 
1077 double TCilindro4T::EntalpiaEntrada(double ASonEnt, double VelEnt, double MasEnt, double ASonCil, double MasCil) {
1078  try {
1079  double xx = 0., yy = 0., ret_val = 0.;
1080 
1081  if(fabs(MasEnt) != 0.) {
1082  xx = (ASonEnt * ASonEnt / ASonCil / ASonCil - 1.) / FGamma1;
1083  yy = VelEnt * VelEnt / ASonCil / ASonCil / 2.;
1084  ret_val = FGamma * MasEnt * (xx + yy) / MasCil;
1085  } else {
1086  ret_val = 0.;
1087  }
1088  return ret_val;
1089  } catch(exception & N) {
1090  std::cout << "ERROR: TCilindro4T::EntalpiaEntrada en cilindro: " << FNumeroCilindro << std::endl;
1091  std::cout << "Tipo de error: " << N.what() << std::endl;
1092  throw Exception(N.what());
1093  }
1094 }
1095 
1096 // ---------------------------------------------------------------------------
1097 // ---------------------------------------------------------------------------
1098 
1099 void TCilindro4T::VariableInicialesCicloACT() {
1100  try {
1101  FPresionMedAdm /= FTimeAcumAct;
1102  FPresionMedEsc /= FTimeAcumAct;
1103  FTimeAcumAct = 0.;
1104 
1105  Ftest_variables[0] = FRegInt; // Speed
1106  if(FMasaPorAdmision != 0)
1107  Ftest_variables[1] = (FMasaPorAdmision - FMasaEGR) * 1000;
1108  // Measured air mass
1109  Ftest_variables[2] = FMasaAtrapada * 1000; // In-cylinder air mass at inlet valve closing
1110  Ftest_variables[3] = __units::degCToK(FTemperature); // In-cylinder temperature at inlet valve closing
1111  Ftest_variables[4] = FMasaFuel * 1e6; // Fuel injected mass
1112  Ftest_variables[5] = FMotor->getInjectionSys().InjectPressure;
1113  // Injection pressure
1114  Ftest_variables[6] = __units::BarToPa(FPresionMedAdm); // Inlet pressure
1115  Ftest_variables[7] = __units::BarToPa(FPresionMedEsc); // Exhaust pressure
1116  // Ftest_variables[8]=FMotor->getGeometria().CDBlowBy; //Blow-by coefficient
1117  // Ftest_variables[9]=FMotor->getPresionAmb()*1e5; //Atmosphere pressure
1118  Ftest_variables[10] = 318.15; // Fuel injection temperature
1119  Ftest_variables[11] = __units::degCToK(FTempPared[0].Culata);
1120  // Cylinder head temperature
1121  Ftest_variables[12] = __units::degCToK(FTempPared[0].Cylinder);
1122  // Cylinder temperature
1123  Ftest_variables[13] = __units::degCToK(FTempPared[0].Piston); // Piston temperature
1124  if(FMotor->getSpeciesModel() == nmCalculoCompleto) {
1125  Ftest_variables[14] = FFraccionMasicaEspecie[5] * 1e6 * 1.587;
1126  // NOx concentration at IVC
1127  Ftest_variables[15] = FFraccionMasicaEspecie[4] * 1e6;
1128  // SOOT concentration at IVC
1129  Ftest_variables[16] = FFraccionMasicaEspecie[6] * 1e6;
1130  // CO concentration at IVC
1131  Ftest_variables[17] = FFraccionMasicaEspecie[3] * 1e6;
1132  // HC concentration at IVC
1133  }
1134  } catch(exception & N) {
1135  std::cout << "ERROR: TCilindro4T::VariableInicialesCicloACT en cilindro: " << FNumeroCilindro << std::endl;
1136  std::cout << "Tipo de error: " << N.what() << std::endl;
1137  throw Exception(N.what());
1138  }
1139 }
1140 
1141 #pragma package(smart_init)
1142 
Constantes.h
TTubo::GetGamma
double GetGamma(int i) const
Gets the specific heat capacities ratio at a given cell.
Definition: TTubo.cpp:5444
TCilindro
Definition: TCilindro.h:59
TTubo::GetPresion
double GetPresion(int i) const
Gets the fluid pressure.
Definition: TTubo.cpp:5468
Exception
Custom exception class.
Definition: Exception.hpp:39
TCCCilindro
Foo.
Definition: TCCCilindro.h:46
TTubo.h
TBloqueMotor
Definition: TBloqueMotor.h:43
pow2
T pow2(T x)
Returns x to the power of 2.
Definition: Math_wam.h:88
TController::Output
virtual double Output(double Time)=0
TTubo::GetAsonido
double GetAsonido(int i) const
Gets the speed of sound.
Definition: TTubo.cpp:5412