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