1 #include"ACT_Sub_DLL.h"
6 double min(
double a,
double b) {
17 void FUNCTION_FOR_INTERPOLATION(
double *interpolated,
double *time_interpolated,
double *CAD_to_interpolate,
18 double *vector_to_interpolate,
int size_interpolated,
int size_to_interpolate,
double speed) {
20 int auxiliar = 0, counter = 0;
22 for(counter = 0; counter < size_interpolated; counter++) {
23 while((CAD_to_interpolate[auxiliar] * 60. / (360. * speed) < time_interpolated[counter])
24 && (counter < size_interpolated) && (auxiliar < size_to_interpolate - 2)) {
25 auxiliar = auxiliar + 1;
28 interpolated[counter] = vector_to_interpolate[auxiliar];
31 interpolated[counter] = vector_to_interpolate[auxiliar - 1] + (vector_to_interpolate[auxiliar] -
32 vector_to_interpolate[auxiliar - 1]) * (time_interpolated[counter] - CAD_to_interpolate[auxiliar - 1] * 60. /
33 (360. * speed)) / (CAD_to_interpolate[auxiliar] * 60. / (360. * speed) - CAD_to_interpolate[auxiliar - 1] * 60. /
42 void CALCULUS_OF_VIRTUAL_VELOCITY(
double *inj_velocity,
double *virtual_velocity,
double *dmf,
double *time_vector,
43 double rofuel,
double dc,
double n_holes,
double nozzle_d,
double D,
int size,
double PI,
double speed,
double *EOI_IM,
44 double inj_num,
double *SOI_IM,
double Piston_D,
double DBowl,
double CTM,
double *CAD,
double Kswirl) {
45 int counter = 0, auxiliar = 0;
49 virt = (
double*) malloc(size *
sizeof(
double));
50 w = (
double*) malloc(size *
sizeof(
double));
52 for(counter = 0; counter < size; counter++) {
53 virtual_velocity[counter] = -9999.;
54 w[counter] = (pow((Piston_D / DBowl),
55 2) * 0.75 * (DBowl / 2) * 2 * PI * speed / 60 * pow((1 / cosh(CAD[counter] / 100)), 3) + 1) * CTM * Kswirl;
58 for(counter = 0; counter < size; counter++) {
60 inj_velocity[counter] = dmf[counter] * 4. / (rofuel * n_holes * dc * PI * pow(nozzle_d, 2));
62 for(i = 0; i < inj_num; i++) {
63 if((time_vector[counter] > (SOI_IM[i] / (speed * 6))) && ((time_vector[counter] - (SOI_IM[i] / (speed * 6))) < 0.001)) {
66 virt[counter] = inj_velocity[counter] * (0.5274 * exp(-(time_vector[counter] - (SOI_IM[i] /
67 (speed * 6))) / 0.00023) + 1.);
68 if(virt[counter] > inj_velocity[counter]) {
69 inj_velocity[counter] = virt[counter];
74 for(auxiliar = counter; auxiliar < size; auxiliar++) {
75 virt[auxiliar] = inj_velocity[counter] * exp(-(time_vector[auxiliar] - time_vector[counter]) / (D * exp((
76 time_vector[auxiliar] - time_vector[counter]) / 0.006)));
77 if(virt[auxiliar] > virtual_velocity[auxiliar]) {
78 virtual_velocity[auxiliar] = virt[auxiliar];
96 void CALCULUS_OF_ACCUMULATED_INJ_RATE(
double *acu_dmf,
double *dmf,
double *time_vector,
int size) {
100 for(counter = 0; counter < size; counter++) {
102 acu_dmf[counter] = 0.;
104 acu_dmf[counter] = acu_dmf[counter - 1] + ((time_vector[counter] - time_vector[counter - 1]) * (dmf[counter] -
105 (dmf[counter] - dmf[counter - 1]) / 2.));
109 if(acu_dmf[size - 1] > 0) {
110 for(counter = 0; counter < size; counter++) {
111 acu_dmf[counter] = acu_dmf[counter] / acu_dmf[size - 1];
119 void STOICHIOMETRY_CONSTANTS(
double HC,
double *Kst1,
double *Kst2,
double *Kst3,
double *Kst4,
double *Kst5,
122 *Kst1 = (1 + HC / 4) * 32 / (12 + HC);
123 *Kst2 = pow(*Kst1, -1);
124 *Kst3 = 11 / (8 * (1 + HC / 4));
125 *Kst4 = 9 * HC / (32 * (1 + HC / 4));
126 *Kst5 = 44 / (12 + HC);
127 *Kst6 = 9 * HC / (12 + HC);
133 void CALCULUS_OF_NUMBER_ELEMENTS(
int *num_i_IM,
double *time_vector,
int size,
double speed,
double *SOI_IM,
134 double *EOI_IM,
int inj_num) {
136 int time_counter = 0, inj_counter = 0;
138 for(inj_counter = 0; inj_counter < inj_num; inj_counter++) {
139 num_i_IM[inj_counter] = 0;
140 for(time_counter = 0; time_counter < size; time_counter++) {
141 if((time_vector[time_counter] >= SOI_IM[inj_counter] * 60. / (360. * speed))
142 && (time_vector[time_counter] <= EOI_IM[inj_counter] * 60. / (360. * speed))) {
143 num_i_IM[inj_counter] = num_i_IM[inj_counter] + 1;
151 void CALCULUS_OF_POI(
double **POI_IM,
double **mfuel_i_IM,
double **mfuel_ij_IM,
double *acu_dmf,
double *time_vector,
152 int size,
double speed,
int *num_i_IM,
int num_j,
double *SOI_IM,
double* EOI_IM,
int inj_num,
154 int counter = 0, inj_counter = 0, aux = 0;
157 auxiliar_numi = (
int*) malloc(inj_num *
sizeof(
int*));
159 for(inj_counter = 0; inj_counter < inj_num; inj_counter++) {
160 auxiliar_numi[inj_counter] = -1;
163 for(inj_counter = 0; inj_counter < inj_num; inj_counter++) {
164 for(counter = 0; counter < size; counter++) {
165 if(time_vector[counter] >= SOI_IM[inj_counter] * 60. / (360. * speed)
166 && time_vector[counter] <= EOI_IM[inj_counter] * 60. / (360. * speed)) {
168 auxiliar_numi[inj_counter] = auxiliar_numi[inj_counter] + 1;
169 aux = auxiliar_numi[inj_counter];
170 POI_IM[inj_counter][aux] = time_vector[counter];
173 if((acu_dmf[counter - 1] < 0.1) && (acu_dmf[counter] >= 0.1)) {
174 elementcontrol[0][0].inj_number = inj_counter;
175 elementcontrol[0][0].num_i = aux;
176 elementcontrol[0][0].inj_number = inj_counter;
177 elementcontrol[1][0].num_i = aux;
178 elementcontrol[1][0].inj_number = inj_counter;
179 elementcontrol[1][0].num_i = aux;
180 elementcontrol[2][0].num_i = aux;
181 elementcontrol[2][0].inj_number = inj_counter;
182 elementcontrol[2][0].num_i = aux;
184 if((acu_dmf[counter - 1] < 0.3) && (acu_dmf[counter] >= 0.3)) {
185 elementcontrol[3][0].inj_number = inj_counter;
186 elementcontrol[3][0].num_i = aux;
187 elementcontrol[3][0].inj_number = inj_counter;
188 elementcontrol[4][0].num_i = aux;
189 elementcontrol[4][0].inj_number = inj_counter;
190 elementcontrol[4][0].num_i = aux;
191 elementcontrol[5][0].num_i = aux;
192 elementcontrol[5][0].inj_number = inj_counter;
193 elementcontrol[5][0].num_i = aux;
195 if((acu_dmf[counter - 1] < 0.5) && (acu_dmf[counter] >= 0.5)) {
196 elementcontrol[6][0].inj_number = inj_counter;
197 elementcontrol[6][0].num_i = aux;
198 elementcontrol[6][0].inj_number = inj_counter;
199 elementcontrol[7][0].num_i = aux;
200 elementcontrol[7][0].inj_number = inj_counter;
201 elementcontrol[7][0].num_i = aux;
202 elementcontrol[8][0].num_i = aux;
203 elementcontrol[8][0].inj_number = inj_counter;
204 elementcontrol[8][0].num_i = aux;
206 if((acu_dmf[counter - 1] < 0.7) && (acu_dmf[counter] >= 0.7)) {
207 elementcontrol[9][0].inj_number = inj_counter;
208 elementcontrol[9][0].num_i = aux;
209 elementcontrol[9][0].inj_number = inj_counter;
210 elementcontrol[10][0].num_i = aux;
211 elementcontrol[10][0].inj_number = inj_counter;
212 elementcontrol[10][0].num_i = aux;
213 elementcontrol[11][0].num_i = aux;
214 elementcontrol[11][0].inj_number = inj_counter;
215 elementcontrol[11][0].num_i = aux;
217 if((acu_dmf[counter - 1] < 0.9) && (acu_dmf[counter] >= 0.9)) {
218 elementcontrol[12][0].inj_number = inj_counter;
219 elementcontrol[12][0].num_i = aux;
220 elementcontrol[12][0].inj_number = inj_counter;
221 elementcontrol[13][0].num_i = aux;
222 elementcontrol[13][0].inj_number = inj_counter;
223 elementcontrol[13][0].num_i = aux;
224 elementcontrol[14][0].num_i = aux;
225 elementcontrol[14][0].inj_number = inj_counter;
226 elementcontrol[14][0].num_i = aux;
230 mfuel_i_IM[inj_counter][aux] = ((acu_dmf[counter] + acu_dmf[counter + 1]) / 2.) - acu_dmf[counter - 1];
231 }
else if(aux == num_i_IM[inj_counter] - 1) {
232 mfuel_i_IM[inj_counter][aux] = acu_dmf[counter + 1] - ((acu_dmf[counter - 1] + acu_dmf[counter]) / 2.);
234 mfuel_i_IM[inj_counter][aux] = ((acu_dmf[counter] + acu_dmf[counter + 1]) / 2.) - ((
235 acu_dmf[counter - 1] + acu_dmf[counter]) / 2.);
242 for(inj_counter = 0; inj_counter < inj_num; inj_counter++) {
243 for(counter = 0; counter < num_i_IM[inj_counter]; counter++) {
244 mfuel_ij_IM[inj_counter][counter] = mfuel_i_IM[inj_counter][counter] / num_j;
254 void CALCULATE_CYCLE(
double *roair,
double *CAD,
double delta_t,
double *V_cyl,
double VTDC,
int counter,
double speed,
255 double *p_cyl,
double *HRF,
double *acu_mf,
double *Mbb,
double *acu_Mbb,
double AFe,
double f,
double mfuel,
256 double mEGR,
double mairIVC,
double *T_cyl,
double HP,
double *Yair,
double *Yfuel,
double *Yburned,
double *U,
257 double *CV,
double *H_cooler,
double *H,
double T_Evaporation_fuel,
double inj_fuel_temp,
double PRECISION_ITERATION,
258 double *defor,
double *Rmixture,
double Atmosphere_press,
double *Gamma,
double PI,
double Runiv,
double Piston_D,
259 double S,
double Crank_L,
double Connecting_Rod_L,
double E,
double Piston_Axis_D,
double Piston_Crown_H,
double DBowl,
260 double VBowl,
double M_Connecting_Rod,
double M_P_R_PA,
double MW_air,
double MW_fuel,
double MW_burned,
261 double C_ESteel,
double C_Mech_Defor,
double CTM,
double WC1A,
double WC1B,
double C2,
double C_MBLBY,
double Cbb,
262 double TPIS,
double TCYL_HEAD,
double TCYL,
double *Qcylhead,
double *Qcyl,
double *Qpis) {
266 double average_Volume = 0.;
268 double average_Temperature = 0.;
270 double average_Pressure = 0.;
282 double delta_CAD = 0.;
289 acu_mf[counter] = acu_mf[counter] * mfuel;
292 acu_mf[counter] = acu_mf[counter] * mfuel;
293 acu_mf[counter - 1] = acu_mf[counter - 1] * mfuel;
296 acu_mf[counter] = acu_mf[counter] * mfuel;
297 acu_mf[counter - 1] = acu_mf[counter - 1] * mfuel;
300 delta_CAD = delta_t * 360. * speed / 60.;
301 V_cyl[counter] = VOLUME(CAD[counter], VTDC, PI, Piston_D, Crank_L, Connecting_Rod_L, E);
305 DEFORMATIONS(&(V_cyl[counter]), &(defor[counter]), p_cyl[counter], CAD[counter], delta_CAD, speed, PI, Piston_D, S,
306 Connecting_Rod_L, E, Piston_Axis_D, Piston_Crown_H, M_Connecting_Rod,
307 M_P_R_PA, C_ESteel, C_Mech_Defor);
309 MASIC_RATIO(&Yair[counter], &Yfuel[counter], &Yburned[counter], &Rmixture[counter], HRF[counter], acu_mf[counter],
310 acu_Mbb[counter], AFe, f, mfuel, mEGR, mairIVC, Runiv, MW_air, MW_fuel,
313 PROPERTIES(&U[counter], &CV[counter], T_cyl[counter], T_cyl[counter], Yair[counter], Yfuel[counter], Yburned[counter]);
315 MCYL = mairIVC + acu_mf[counter] - acu_Mbb[counter];
316 roair[counter] = MCYL / V_cyl[counter];
317 p_cyl[counter] = MCYL * Rmixture[counter] * T_cyl[counter] / V_cyl[counter];
321 *Gamma = (CV[counter - 1] + Rmixture[counter - 1]) / CV[counter - 1];
323 Mbb[counter] = BLOW_BY(p_cyl[counter - 1], T_cyl[counter - 1], Rmixture[counter - 1], delta_CAD, speed, *Gamma,
324 Atmosphere_press, Piston_D, C_MBLBY, Cbb);
326 acu_Mbb[counter] = acu_Mbb[counter - 1] + Mbb[counter];
328 DEFORMATIONS(&(V_cyl[counter]), &(defor[counter]), p_cyl[counter - 1], CAD[counter], delta_CAD, speed, PI, Piston_D, S,
329 Connecting_Rod_L, E, Piston_Axis_D, Piston_Crown_H, M_Connecting_Rod,
330 M_P_R_PA, C_ESteel, C_Mech_Defor);
332 MASIC_RATIO(&Yair[counter], &Yfuel[counter], &Yburned[counter], &Rmixture[counter], HRF[counter - 1], acu_mf[counter],
333 acu_Mbb[counter], AFe, f, mfuel, mEGR, mairIVC, Runiv, MW_air, MW_fuel,
336 PROPERTIES(&UANT, &CV[counter], T_cyl[counter - 1], T_cyl[counter - 1], Yair[counter], Yfuel[counter],
339 MCYL = mairIVC + acu_mf[counter] - acu_Mbb[counter];
341 roair[counter] = MCYL / V_cyl[counter];
343 average_Volume = 0.5 * (V_cyl[counter] + V_cyl[counter - 1]);
345 H_cooler[counter] = HEAT_COOLER(p_cyl[counter - 1], p_cyl[0], T_cyl[counter - 1], T_cyl[0], average_Volume, V_cyl[0],
346 delta_CAD, speed, VTDC, &H[counter], PI, Piston_D, S, DBowl, VBowl, CTM,
347 WC1A, WC1B, C2, TPIS, TCYL_HEAD, TCYL, CAD[counter], Qcylhead, Qcyl, Qpis, counter);
350 QC = (HRF[counter - 1] - HRF[counter - 2]) * QT;
351 DU = -(p_cyl[counter - 1] * (V_cyl[counter] - V_cyl[counter - 1])) + QC - H_cooler[counter];
353 T_cyl[counter] = T_cyl[counter - 1] + DU / MCYL / CV[counter];
355 ecg = -21776.6 * Runiv / MW_fuel + Runiv / MW_fuel * (-4.5826 * T_Evaporation_fuel + 0.12428 / 2. * pow(
356 T_Evaporation_fuel, 2.) - 0.00007233 / 3. * pow(T_Evaporation_fuel, 3.) + 0.000000016269 / 4. * pow(
357 T_Evaporation_fuel, 4.) - 26067.28 / T_Evaporation_fuel);
359 ent_ref = (Runiv / MW_fuel) * ((-4.5826 * 298.) + (0.12428 / 2. * pow(298., 2.)) - (0.00007233 / 3. * pow(298.,
360 3.)) + (0.000000016269 / 4. * pow(298., 4.)) - (26067.28 * pow(298., -1.)));
363 HFTiny = -1852564 + 2195 * (inj_fuel_temp);
369 T_cyl[counter] = T_cyl[counter] - AERR / CV[counter] / MCYL;
370 p_cyl[counter] = MCYL * Rmixture[counter] * T_cyl[counter] / V_cyl[counter];
371 PROPERTIES(&U[counter], &CV[counter], T_cyl[counter], T_cyl[counter - 1], Yair[counter], Yfuel[counter],
374 average_Temperature = 0.5 * (T_cyl[counter] + T_cyl[counter - 1]);
375 average_Pressure = 0.5 * (p_cyl[counter] + p_cyl[counter - 1]);
376 H_cooler[counter] = HEAT_COOLER(average_Pressure, p_cyl[0], average_Temperature, T_cyl[0], average_Volume, V_cyl[0],
377 delta_CAD, speed, VTDC, &H[counter], PI, Piston_D, S, DBowl, VBowl,
378 CTM, WC1A, WC1B, C2, TPIS, TCYL_HEAD, TCYL, CAD[counter], Qcylhead, Qcyl, Qpis, counter);
380 Wi = -average_Pressure * (V_cyl[counter] - V_cyl[counter - 1]);
382 PROPERTIES_FUEL(&uf, T_cyl[counter]);
384 AERR = (MCYL * (U[counter] - UANT)) - Wi - QC + H_cooler[counter] - ((HFTiny - uf) *
385 (acu_mf[counter] - acu_mf[counter - 1])) + (Rmixture[counter] * T_cyl[counter] * Mbb[counter]);
388 a = fabs(AERR / U[counter]);
392 }
while(a > PRECISION_ITERATION);
396 acu_mf[counter] = acu_mf[counter] / mfuel;
399 acu_mf[counter] = acu_mf[counter] / mfuel;
400 acu_mf[counter - 1] = acu_mf[counter - 1] / mfuel;
403 acu_mf[counter] = acu_mf[counter] / mfuel;
404 acu_mf[counter - 1] = acu_mf[counter - 1] / mfuel;
409 void CALCULUS_OF_MEAN_VARIABLES(
double *p_cyl,
double *T_cyl,
double *dp_da_cyl,
double *CAD,
double *pmax,
410 double *Tmax,
double *dp_da_max,
double *p_exit,
double *T_exit,
int size) {
412 int counter = 0, iter = 0, n = 0;
415 dp_da_aux = (
double*) malloc(size *
sizeof(
double));
418 for(counter = 0; counter < size; counter++) {
420 if(p_cyl[counter] > *pmax) {
421 *pmax = p_cyl[counter];
423 if(T_cyl[counter] > *Tmax) {
424 *Tmax = T_cyl[counter];
427 dp_da_aux[counter] = (p_cyl[counter] - p_cyl[counter - 1]) / (CAD[counter] - CAD[counter - 1]);
430 if(counter == size - 1) {
431 *p_exit = p_cyl[counter];
432 *T_exit = T_cyl[counter];
437 for(counter = 0; counter < size; counter++)
438 dp_da_cyl[counter] = dp_da_aux[counter];
448 for(counter = 0; counter < size - 1; counter++) {
449 if(dp_da_cyl[counter] > *dp_da_max) {
450 *dp_da_max = dp_da_cyl[counter];
460 double VOLUME(
double CAD,
double VTDC,
double PI,
double Piston_D,
double Crank_L,
double Connecting_Rod_L,
double E) {
462 double V_cyl = 0., AREA = 0., AUX = 0., A = 0.;
465 AREA = PI * Piston_D * Piston_D / 4.;
466 AUX = (Crank_L * sqrt(pow(1. + 1. / (Crank_L / Connecting_Rod_L), 2) - pow(E / Crank_L,
467 2.)) - Crank_L * (cos(A) + sqrt(pow(1. / (Crank_L / Connecting_Rod_L), 2.) - pow(sin(A) - E / Crank_L, 2.))));
468 V_cyl = (VTDC + AREA * AUX);
475 void DEFORMATIONS(
double *V_cyl,
double *DEFOR,
double p_cyl,
double CAD,
double delta_CAD,
double speed,
double PI,
476 double Piston_D,
double S,
double Connecting_Rod_L,
double E,
double Piston_Axis_D,
double Piston_Crown_H,
477 double M_Connecting_Rod,
double M_P_R_PA,
double C_ESteel,
double C_Mech_Defor) {
484 double Lanterior = 0.;
486 double Lposterior = 0.;
492 AVp = (PI * pow(Piston_D, 2.) / 4.) * C_Mech_Defor * (p_cyl / C_ESteel) * pow(Piston_D / Piston_Axis_D,
493 2.) * (Piston_Crown_H + Connecting_Rod_L + S / 2.);
495 auxalfa = (CAD - delta_CAD) * PI / 180.;
496 Lanterior = (sqrt(pow((S / 2.) + Connecting_Rod_L, 2.) - pow(E,
497 2.))) - ((S / 2.) * cos(auxalfa) + sqrt(pow(Connecting_Rod_L, 2.) - pow((S / 2.) * sin(auxalfa) - E, 2.)));
499 auxalfa = (CAD + delta_CAD) * PI / 180.;
500 Lposterior = (sqrt(pow((S / 2.) + Connecting_Rod_L, 2.) - pow(E,
501 2.))) - ((S / 2.) * cos(auxalfa) + sqrt(pow(Connecting_Rod_L, 2.) - pow((S / 2.) * sin(auxalfa) - E, 2.)));
503 auxalfa = (CAD) * PI / 180.;
504 Lactual = (sqrt(pow((S / 2.) + Connecting_Rod_L, 2.) - pow(E,
505 2.))) - ((S / 2.) * cos(auxalfa) + sqrt(pow(Connecting_Rod_L, 2.) - pow((S / 2.) * sin(auxalfa) - E, 2.)));
507 Acel = (Lanterior + Lposterior - Lactual * 2.) / pow((delta_CAD / (6. * speed)), 2.);
509 Msist = 0.33 * M_Connecting_Rod + M_P_R_PA;
511 AVi = Msist * Acel * C_Mech_Defor * 1. / C_ESteel * pow(Piston_D / Piston_Axis_D,
512 2.) * (Piston_Crown_H + Connecting_Rod_L + S / 2.);
514 *DEFOR = (AVp - AVi) * 1000000;
515 *V_cyl = *V_cyl + (AVp - AVi);
521 void MASIC_RATIO(
double *Yair,
double *Yfuel,
double *Yburned,
double *Rmixture,
double HRF,
double acu_mf,
522 double acu_Mbb,
double AFe,
double f,
double mfuel,
double mEGR,
double mairIVC,
double Runiv,
double MW_air,
523 double MW_fuel,
double MW_burned) {
525 double Yq2 = 0., Yq0 = 0.;
527 Yq2 = (1. + (AFe)) / (1. + (f) - ((acu_Mbb / 2.) / mfuel) * (1. / (1. + (mEGR) / (mairIVC))));
528 Yq0 = Yq2 * (mEGR) / (mairIVC);
530 *Yburned = (((mfuel + mfuel * AFe) * HRF) + ((mEGR) * Yq2) - ((acu_Mbb / 2.) * Yq0)) /
531 (acu_mf + mairIVC - acu_Mbb / 2.);
533 if(acu_mf >= mfuel * HRF) {
534 *Yfuel = (acu_mf - mfuel * HRF) / (acu_mf + mairIVC - acu_Mbb);
539 *Yair = 1. - *Yfuel - *Yburned;
540 *Rmixture = *Yair * Runiv / MW_air + *Yfuel * Runiv / MW_fuel + *Yburned * Runiv / MW_burned;
546 void PROPERTIES(
double *u,
double *CV,
double T_cyl,
double T_cyl_pre,
double Yair,
double Yfuel,
double Yburned) {
548 double cva = 0., cvf = 0., cvq = 0., average_Temperature = 0., ua = 0., uf = 0., uq = 0.;
550 average_Temperature = 0.5 * (T_cyl + T_cyl_pre);
552 cva = (-10.4199 * pow(average_Temperature, 0.5)) + 2522.88 - (67227.1 * pow(average_Temperature,
553 -0.5)) + (917124.4 * pow(average_Temperature, -1.)) - (4174853.6 * pow(average_Temperature, -1.5));
554 cvf = -256.4 + (6.95372 * average_Temperature) - (0.00404715 * pow(average_Temperature,
555 2.)) + (0.000000910259 * pow(average_Temperature, 3.)) + (1458487. * pow(average_Temperature, -2.));
556 cvq = 641.154 + (0.43045 * average_Temperature) - (0.0001125 * pow(average_Temperature,
557 2.)) + (0.000000008979 * pow(average_Temperature, 3.));
559 *CV = cva * Yair + cvf * Yfuel + cvq * Yburned;
561 ua = -4193697.9 - (6.9466 * pow(T_cyl, 1.5)) + (2522.88 * T_cyl) - (134454.16 * pow(T_cyl,
562 0.5)) + (917124.39 * log(T_cyl)) + (8349707.14 * pow(T_cyl, -0.5));
563 uf = -1445686.1 - (256.4 * T_cyl) + (3.47686 * pow(T_cyl, 2.)) - (0.00134905 * pow(T_cyl,
564 3.)) + (0.000000227565 * pow(T_cyl, 4.)) - (1458487. * pow(T_cyl, -1.));
565 uq = -3251495. + (1028.75 * T_cyl) - (0.15377 * pow(T_cyl, 2.)) + (0.000067895 * pow(T_cyl, 3.));
567 *u = ua * Yair + uf * Yfuel + uq * Yburned;
572 void PROPERTIES_FUEL(
double *uf,
double T_cyl) {
574 *uf = -1445686.1 - (256.4 * T_cyl) + (3.47686 * pow(T_cyl, 2.)) - (0.00134905 * pow(T_cyl,
575 3.)) + (0.000000227565 * pow(T_cyl, 4.)) - (1458487. * pow(T_cyl, -1.));
580 double HEAT_COOLER(
double p_cyl,
double pressureIVC,
double T_cyl,
double temperatureIVC,
double average_Volume,
581 double volumeIVC,
double delta_CAD,
double speed,
double VTDC,
double *H,
double PI,
double Piston_D,
double S,
582 double DBowl,
double VBowl,
double CTM,
double WC1A,
double WC1B,
double C2,
double TPIS,
double TCYL_HEAD,
double TCYL,
583 double CAD,
double *Qcylhead,
double *Qcyl,
double *Qpis,
int counter)
591 double Cylinder_capacity = 0.;
592 double Piston_area = 0.;
593 double Cylinder_head_area = 0.;
594 double H_cooler = 0.;
596 cm = 2. * S * speed / 60.;
597 C1 = CALCULATE_C1(cm, CTM, WC1A, WC1B, Piston_D, DBowl, speed, CAD, PI);
599 PA = pressureIVC * pow(volumeIVC / average_Volume, 1.36);
601 Cylinder_capacity = PI * Piston_D * Piston_D * S / 4.;
603 comb = ((Cylinder_capacity * temperatureIVC) / (pressureIVC * volumeIVC)) * (p_cyl - PA);
605 *H = 0.012 * pow(Piston_D, -0.2) * pow(T_cyl, -0.53) * pow(p_cyl, 0.8) * pow((C1 + C2 * comb), 0.8);
607 CALCULATE_AREAS(&Piston_area, &Cylinder_head_area, PI, Piston_D, DBowl, VBowl);
609 H_cooler = *H * (Piston_area * (T_cyl - TPIS) + Cylinder_head_area * (T_cyl - TCYL_HEAD) + 4. * ((
610 average_Volume - VTDC) / Piston_D * (T_cyl - TCYL))) * (delta_CAD) / 6. / speed;
612 Qcylhead[counter] = *H * Cylinder_head_area * (T_cyl - TCYL_HEAD) * (delta_CAD) / 6. / speed;
613 Qcyl[counter] = *H * 4. * ((average_Volume - VTDC) / Piston_D * (T_cyl - TCYL)) * (delta_CAD) / 6. / speed;
614 Qpis[counter] = *H * Piston_area * (T_cyl - TPIS) * (delta_CAD) / 6. / speed;
621 double CALCULATE_C1(
double cm,
double CTM,
double WC1A,
double WC1B,
double Piston_D,
double DBowl,
double speed,
622 double CAD,
double PI) {
627 double ratio_CTM = 0.;
630 KCTM = exp(-0.200679 * pow(CTM, 0.431262));
632 ratio_CTM = pow(DBowl / Piston_D, 2) * (1 / KCTM);
633 x_alfa = ratio_CTM + 1 / (pow(cosh(CAD / 100), 40) + ratio_CTM / (1 - ratio_CTM));
635 cu = x_alfa * (2 * PI * speed / 60) * CTM * pow(Piston_D / DBowl, 2) * KCTM * (DBowl / 2);
637 C1 = WC1A * cm + WC1B * cu;
644 void CALCULATE_AREAS(
double *Piston_area,
double *Cylinder_head_area,
double PI,
double Piston_D,
double DBowl,
647 *Piston_area = PI * pow(Piston_D, 2.) / 4. + VBowl * 4. / DBowl;
648 *Cylinder_head_area = PI * pow(Piston_D, 2.) / 4.;
653 double BLOW_BY(
double p_cyl,
double T_cyl,
double Rmixture,
double delta_CAD,
double speed,
double Gamma,
654 double Atmosphere_press,
double Piston_D,
double C_MBLBY,
double Cbb) {
657 double Pressure_up = 0.;
658 double Pressure_down = 0.;
659 double Pressure_critic = 0.;
662 if(p_cyl > Atmosphere_press) {
664 Pressure_down = Atmosphere_press;
666 Pressure_up = Atmosphere_press;
667 Pressure_down = p_cyl;
670 Pressure_critic = Pressure_up * pow(2. / (Gamma + 1.), (Gamma / (Gamma - 1.)));
672 if(Pressure_down < Pressure_critic) {
673 Pressure_down = Pressure_critic;
676 C_Z = ((2. * Gamma) / (Gamma - 1.)) * (pow(Pressure_down / Pressure_up, 2. / Gamma) - pow(Pressure_down / Pressure_up,
677 ((Gamma + 1.) / Gamma)));
678 BBy = Cbb * C_MBLBY * Piston_D * Pressure_up * pow(C_Z / (Rmixture * T_cyl), 0.5);
680 if(Atmosphere_press > p_cyl) {
684 BBy = BBy * delta_CAD / 6. / speed;
691 void CALCULUS_OF_IMP_HP(
double *complete_p_cyl,
double *complete_CAD,
double *p_cyl,
double *V_cyl,
692 double *complete_V_cyl,
double *complete_deform,
double *WI_HP,
double *IMP_HP,
int complete_size,
693 int complete_prev_size,
double delta_t,
double speed,
int size,
double IVC,
double EVO,
double VTDC,
694 double Cylinder_capacity,
double PI,
double Piston_D,
double S,
double Crank_L,
double Connecting_Rod_L,
double E,
695 double Piston_Axis_D,
double Piston_Crown_H,
double M_Connecting_Rod,
double M_P_R_PA,
double C_ESteel,
696 double C_Mech_Defor,
double inlet_pres,
double exhaust_pres) {
698 double delta_CAD = 0.;
703 delta_CAD = delta_t * 360. * speed / 60.;
704 complete_CAD[0] = IVC - complete_prev_size * delta_t * 360. * speed / 60.;
706 for(counter = 1; counter < complete_size; counter++) {
707 complete_CAD[counter] = complete_CAD[counter - 1] + delta_t * 360. * speed / 60.;
710 for(counter = 0; counter < complete_size; counter++) {
711 if(counter < complete_prev_size) {
712 complete_p_cyl[counter] = ((inlet_pres - p_cyl[0]) / (-180 - IVC)) * (complete_CAD[counter] - (-180)) + inlet_pres;
714 if((counter >= complete_prev_size) && (counter < complete_prev_size + size)) {
715 complete_p_cyl[counter] = p_cyl[counter - complete_prev_size];
717 if(counter >= complete_prev_size + size) {
718 complete_p_cyl[counter] = ((exhaust_pres - p_cyl[size - 1]) / (180 - EVO)) * (complete_CAD[counter] -
719 (180)) + exhaust_pres;
723 for(counter = 0; counter < complete_size; counter++) {
724 if(counter < complete_prev_size) {
725 complete_V_cyl[counter] = VOLUME(complete_CAD[counter], VTDC, PI, Piston_D, Crank_L, Connecting_Rod_L, E);
726 DEFORMATIONS(&complete_V_cyl[counter], &complete_deform[counter], complete_p_cyl[counter], complete_CAD[counter],
727 delta_CAD, speed, PI, Piston_D, S, Connecting_Rod_L, E, Piston_Axis_D,
728 Piston_Crown_H, M_Connecting_Rod, M_P_R_PA, C_ESteel, C_Mech_Defor);
730 if((counter >= complete_prev_size) && (counter < complete_prev_size + size)) {
731 complete_V_cyl[counter] = V_cyl[counter - complete_prev_size];
733 if(counter >= complete_prev_size + size) {
734 complete_V_cyl[counter] = VOLUME(complete_CAD[counter], VTDC, PI, Piston_D, Crank_L, Connecting_Rod_L, E);
735 DEFORMATIONS(&complete_V_cyl[counter], &complete_deform[counter], complete_p_cyl[counter], complete_CAD[counter],
736 delta_CAD, speed, PI, Piston_D, S, Connecting_Rod_L, E, Piston_Axis_D,
737 Piston_Crown_H, M_Connecting_Rod, M_P_R_PA, C_ESteel, C_Mech_Defor);
741 *WI_HP = *WI_HP + (complete_p_cyl[counter] + complete_p_cyl[counter - 1]) * 0.5 * (complete_V_cyl[counter] -
742 complete_V_cyl[counter - 1]);
745 *IMP_HP = *WI_HP / Cylinder_capacity;
751 void FUNCTION_NOX(
double *YNOeq_value,
double *KdYNO_value,
double **YNOeq,
double **KdYNO,
double temperature,
752 double mO2,
double mtotal) {
754 double di = 0., dj = 0.;
756 YNOeq[0][0] = 0.000000;
757 KdYNO[0][0] = 0.000000;
758 YNOeq[0][1] = 0.000034;
759 KdYNO[0][1] = 0.000000;
760 YNOeq[0][2] = 0.000054;
761 KdYNO[0][2] = 0.000000;
762 YNOeq[0][3] = 0.000074;
763 KdYNO[0][3] = 0.000000;
764 YNOeq[0][4] = 0.000000;
765 KdYNO[0][4] = 0.000000;
766 YNOeq[0][5] = 0.000000;
767 KdYNO[0][5] = 0.000000;
768 YNOeq[0][6] = 0.000000;
769 KdYNO[0][6] = 0.000000;
770 YNOeq[0][7] = 0.000000;
771 KdYNO[0][7] = 0.000000;
772 YNOeq[0][8] = 0.000000;
773 KdYNO[0][8] = 0.000000;
774 YNOeq[0][9] = 0.000000;
775 KdYNO[0][9] = 0.000000;
776 YNOeq[0][10] = 0.000000;
777 KdYNO[0][10] = 0.000000;
778 YNOeq[0][11] = 0.000000;
779 KdYNO[0][11] = 0.000000;
780 YNOeq[1][0] = 0.000000;
781 KdYNO[1][0] = 0.000000;
782 YNOeq[1][1] = 0.000068;
783 KdYNO[1][1] = 0.000000;
784 YNOeq[1][2] = 0.000090;
785 KdYNO[1][2] = 0.000000;
786 YNOeq[1][3] = 0.000113;
787 KdYNO[1][3] = 0.000000;
788 YNOeq[1][4] = 0.000135;
789 KdYNO[1][4] = 0.000000;
790 YNOeq[1][5] = 0.000160;
791 KdYNO[1][5] = 0.000000;
792 YNOeq[1][6] = 0.000186;
793 KdYNO[1][6] = 0.000000;
794 YNOeq[1][7] = 0.000210;
795 KdYNO[1][7] = 0.000000;
796 YNOeq[1][8] = 0.000000;
797 KdYNO[1][8] = 0.000000;
798 YNOeq[1][9] = 0.000000;
799 KdYNO[1][9] = 0.000000;
800 YNOeq[1][10] = 0.000000;
801 KdYNO[1][10] = 0.000000;
802 YNOeq[1][11] = 0.000000;
803 KdYNO[1][11] = 0.000000;
804 YNOeq[2][0] = 0.000000;
805 KdYNO[2][0] = 0.000000;
806 YNOeq[2][1] = 0.000080;
807 KdYNO[2][1] = 0.000000;
808 YNOeq[2][2] = 0.000131;
809 KdYNO[2][2] = 0.000000;
810 YNOeq[2][3] = 0.000168;
811 KdYNO[2][3] = 0.000000;
812 YNOeq[2][4] = 0.000197;
813 KdYNO[2][4] = 0.000000;
814 YNOeq[2][5] = 0.000224;
815 KdYNO[2][5] = 0.000000;
816 YNOeq[2][6] = 0.000251;
817 KdYNO[2][6] = 0.000000;
818 YNOeq[2][7] = 0.000277;
819 KdYNO[2][7] = 0.000000;
820 YNOeq[2][8] = 0.000304;
821 KdYNO[2][8] = 0.000000;
822 YNOeq[2][9] = 0.000000;
823 KdYNO[2][9] = 0.000000;
824 YNOeq[2][10] = 0.000000;
825 KdYNO[2][10] = 0.000000;
826 YNOeq[2][11] = 0.000000;
827 KdYNO[2][11] = 0.000000;
828 YNOeq[3][0] = 0.000000;
829 KdYNO[3][0] = 0.000000;
830 YNOeq[3][1] = 0.000056;
831 KdYNO[3][1] = 0.000000;
832 YNOeq[3][2] = 0.000111;
833 KdYNO[3][2] = 0.000000;
834 YNOeq[3][3] = 0.000165;
835 KdYNO[3][3] = 0.000000;
836 YNOeq[3][4] = 0.000220;
837 KdYNO[3][4] = 0.000000;
838 YNOeq[3][5] = 0.000274;
839 KdYNO[3][5] = 0.000000;
840 YNOeq[3][6] = 0.000329;
841 KdYNO[3][6] = 0.000000;
842 YNOeq[3][7] = 0.000366;
843 KdYNO[3][7] = 0.000000;
844 YNOeq[3][8] = 0.000396;
845 KdYNO[3][8] = 0.000000;
846 YNOeq[3][9] = 0.000000;
847 KdYNO[3][9] = 0.000000;
848 YNOeq[3][10] = 0.000000;
849 KdYNO[3][10] = 0.000000;
850 YNOeq[3][11] = 0.000000;
851 KdYNO[3][11] = 0.000000;
852 YNOeq[4][0] = 0.000001;
853 KdYNO[4][0] = 0.000000;
854 YNOeq[4][1] = 0.000150;
855 KdYNO[4][1] = 0.000000;
856 YNOeq[4][2] = 0.000226;
857 KdYNO[4][2] = 0.000000;
858 YNOeq[4][3] = 0.000280;
859 KdYNO[4][3] = 0.000000;
860 YNOeq[4][4] = 0.000324;
861 KdYNO[4][4] = 0.000000;
862 YNOeq[4][5] = 0.000368;
863 KdYNO[4][5] = 0.000000;
864 YNOeq[4][6] = 0.000411;
865 KdYNO[4][6] = 0.000000;
866 YNOeq[4][7] = 0.000455;
867 KdYNO[4][7] = 0.000000;
868 YNOeq[4][8] = 0.000500;
869 KdYNO[4][8] = 0.000000;
870 YNOeq[4][9] = 0.000551;
871 KdYNO[4][9] = 0.000000;
872 YNOeq[4][10] = 0.000000;
873 KdYNO[4][10] = 0.000000;
874 YNOeq[4][11] = 0.000000;
875 KdYNO[4][11] = 0.000000;
876 YNOeq[5][0] = 0.000001;
877 KdYNO[5][0] = 0.000000;
878 YNOeq[5][1] = 0.000089;
879 KdYNO[5][1] = 0.000000;
880 YNOeq[5][2] = 0.000173;
881 KdYNO[5][2] = 0.000000;
882 YNOeq[5][3] = 0.000258;
883 KdYNO[5][3] = 0.000000;
884 YNOeq[5][4] = 0.000342;
885 KdYNO[5][4] = 0.000000;
886 YNOeq[5][5] = 0.000426;
887 KdYNO[5][5] = 0.000000;
888 YNOeq[5][6] = 0.000510;
889 KdYNO[5][6] = 0.000000;
890 YNOeq[5][7] = 0.000590;
891 KdYNO[5][7] = 0.000000;
892 YNOeq[5][8] = 0.000644;
893 KdYNO[5][8] = 0.000000;
894 YNOeq[5][9] = 0.000696;
895 KdYNO[5][9] = 0.000000;
896 YNOeq[5][10] = 0.000000;
897 KdYNO[5][10] = 0.000000;
898 YNOeq[5][11] = 0.000000;
899 KdYNO[5][11] = 0.000000;
900 YNOeq[6][0] = 0.000002;
901 KdYNO[6][0] = 0.000000;
902 YNOeq[6][1] = 0.000235;
903 KdYNO[6][1] = 0.000000;
904 YNOeq[6][2] = 0.000351;
905 KdYNO[6][2] = 0.000000;
906 YNOeq[6][3] = 0.000438;
907 KdYNO[6][3] = 0.000000;
908 YNOeq[6][4] = 0.000511;
909 KdYNO[6][4] = 0.000000;
910 YNOeq[6][5] = 0.000585;
911 KdYNO[6][5] = 0.000000;
912 YNOeq[6][6] = 0.000658;
913 KdYNO[6][6] = 0.000000;
914 YNOeq[6][7] = 0.000727;
915 KdYNO[6][7] = 0.000000;
916 YNOeq[6][8] = 0.000790;
917 KdYNO[6][8] = 0.000000;
918 YNOeq[6][9] = 0.000858;
919 KdYNO[6][9] = 0.000000;
920 YNOeq[6][10] = 0.000931;
921 KdYNO[6][10] = 0.000000;
922 YNOeq[6][11] = 0.000000;
923 KdYNO[6][11] = 0.000000;
924 YNOeq[7][0] = 0.000005;
925 KdYNO[7][0] = 0.000000;
926 YNOeq[7][1] = 0.000136;
927 KdYNO[7][1] = 0.000000;
928 YNOeq[7][2] = 0.000254;
929 KdYNO[7][2] = 0.000000;
930 YNOeq[7][3] = 0.000372;
931 KdYNO[7][3] = 0.000000;
932 YNOeq[7][4] = 0.000489;
933 KdYNO[7][4] = 0.000000;
934 YNOeq[7][5] = 0.000607;
935 KdYNO[7][5] = 0.000000;
936 YNOeq[7][6] = 0.000725;
937 KdYNO[7][6] = 0.000000;
938 YNOeq[7][7] = 0.000843;
939 KdYNO[7][7] = 0.000001;
940 YNOeq[7][8] = 0.000961;
941 KdYNO[7][8] = 0.000001;
942 YNOeq[7][9] = 0.001066;
943 KdYNO[7][9] = 0.000001;
944 YNOeq[7][10] = 0.001140;
945 KdYNO[7][10] = 0.000001;
946 YNOeq[7][11] = 0.000000;
947 KdYNO[7][11] = 0.000000;
948 YNOeq[8][0] = 0.000010;
949 KdYNO[8][0] = 0.000000;
950 YNOeq[8][1] = 0.000352;
951 KdYNO[8][1] = 0.000000;
952 YNOeq[8][2] = 0.000525;
953 KdYNO[8][2] = 0.000000;
954 YNOeq[8][3] = 0.000653;
955 KdYNO[8][3] = 0.000001;
956 YNOeq[8][4] = 0.000762;
957 KdYNO[8][4] = 0.000001;
958 YNOeq[8][5] = 0.000871;
959 KdYNO[8][5] = 0.000001;
960 YNOeq[8][6] = 0.000981;
961 KdYNO[8][6] = 0.000001;
962 YNOeq[8][7] = 0.001081;
963 KdYNO[8][7] = 0.000002;
964 YNOeq[8][8] = 0.001177;
965 KdYNO[8][8] = 0.000002;
966 YNOeq[8][9] = 0.001270;
967 KdYNO[8][9] = 0.000002;
968 YNOeq[8][10] = 0.001362;
969 KdYNO[8][10] = 0.000002;
970 YNOeq[8][11] = 0.001474;
971 KdYNO[8][11] = 0.000003;
972 YNOeq[9][0] = 0.000017;
973 KdYNO[9][0] = 0.000000;
974 YNOeq[9][1] = 0.000210;
975 KdYNO[9][1] = 0.000001;
976 YNOeq[9][2] = 0.000376;
977 KdYNO[9][2] = 0.000001;
978 YNOeq[9][3] = 0.000543;
979 KdYNO[9][3] = 0.000002;
980 YNOeq[9][4] = 0.000709;
981 KdYNO[9][4] = 0.000003;
982 YNOeq[9][5] = 0.000875;
983 KdYNO[9][5] = 0.000003;
984 YNOeq[9][6] = 0.001041;
985 KdYNO[9][6] = 0.000004;
986 YNOeq[9][7] = 0.001208;
987 KdYNO[9][7] = 0.000005;
988 YNOeq[9][8] = 0.001374;
989 KdYNO[9][8] = 0.000005;
990 YNOeq[9][9] = 0.001523;
991 KdYNO[9][9] = 0.000006;
992 YNOeq[9][10] = 0.001629;
993 KdYNO[9][10] = 0.000007;
994 YNOeq[9][11] = 0.001748;
995 KdYNO[9][11] = 0.000009;
996 YNOeq[10][0] = 0.000029;
997 KdYNO[10][0] = 0.000000;
998 YNOeq[10][1] = 0.000508;
999 KdYNO[10][1] = 0.000003;
1000 YNOeq[10][2] = 0.000755;
1001 KdYNO[10][2] = 0.000005;
1002 YNOeq[10][3] = 0.000936;
1003 KdYNO[10][3] = 0.000006;
1004 YNOeq[10][4] = 0.001085;
1005 KdYNO[10][4] = 0.000008;
1006 YNOeq[10][5] = 0.001234;
1007 KdYNO[10][5] = 0.000010;
1008 YNOeq[10][6] = 0.001382;
1009 KdYNO[10][6] = 0.000012;
1010 YNOeq[10][7] = 0.001530;
1011 KdYNO[10][7] = 0.000014;
1012 YNOeq[10][8] = 0.001676;
1013 KdYNO[10][8] = 0.000016;
1014 YNOeq[10][9] = 0.001806;
1015 KdYNO[10][9] = 0.000018;
1016 YNOeq[10][10] = 0.001926;
1017 KdYNO[10][10] = 0.000020;
1018 YNOeq[10][11] = 0.002061;
1019 KdYNO[10][11] = 0.000025;
1020 YNOeq[11][0] = 0.000045;
1021 KdYNO[11][0] = 0.000001;
1022 YNOeq[11][1] = 0.000307;
1023 KdYNO[11][1] = 0.000005;
1024 YNOeq[11][2] = 0.000535;
1025 KdYNO[11][2] = 0.000010;
1026 YNOeq[11][3] = 0.000763;
1027 KdYNO[11][3] = 0.000015;
1028 YNOeq[11][4] = 0.000992;
1029 KdYNO[11][4] = 0.000020;
1030 YNOeq[11][5] = 0.001220;
1031 KdYNO[11][5] = 0.000026;
1032 YNOeq[11][6] = 0.001448;
1033 KdYNO[11][6] = 0.000031;
1034 YNOeq[11][7] = 0.001676;
1035 KdYNO[11][7] = 0.000036;
1036 YNOeq[11][8] = 0.001905;
1037 KdYNO[11][8] = 0.000041;
1038 YNOeq[11][9] = 0.002108;
1039 KdYNO[11][9] = 0.000046;
1040 YNOeq[11][10] = 0.002269;
1041 KdYNO[11][10] = 0.000056;
1042 YNOeq[11][11] = 0.002411;
1043 KdYNO[11][11] = 0.000067;
1044 YNOeq[12][0] = 0.000067;
1045 KdYNO[12][0] = 0.000002;
1046 YNOeq[12][1] = 0.000709;
1047 KdYNO[12][1] = 0.000024;
1048 YNOeq[12][2] = 0.001051;
1049 KdYNO[12][2] = 0.000036;
1050 YNOeq[12][3] = 0.001300;
1051 KdYNO[12][3] = 0.000047;
1052 YNOeq[12][4] = 0.001498;
1053 KdYNO[12][4] = 0.000059;
1054 YNOeq[12][5] = 0.001696;
1055 KdYNO[12][5] = 0.000070;
1056 YNOeq[12][6] = 0.001894;
1057 KdYNO[12][6] = 0.000082;
1058 YNOeq[12][7] = 0.002101;
1059 KdYNO[12][7] = 0.000098;
1060 YNOeq[12][8] = 0.002312;
1061 KdYNO[12][8] = 0.000116;
1062 YNOeq[12][9] = 0.002488;
1063 KdYNO[12][9] = 0.000131;
1064 YNOeq[12][10] = 0.002642;
1065 KdYNO[12][10] = 0.000143;
1066 YNOeq[12][11] = 0.002805;
1067 KdYNO[12][11] = 0.000166;
1068 YNOeq[13][0] = 0.000093;
1069 KdYNO[13][0] = 0.000007;
1070 YNOeq[13][1] = 0.000441;
1071 KdYNO[13][1] = 0.000035;
1072 YNOeq[13][2] = 0.000745;
1073 KdYNO[13][2] = 0.000068;
1074 YNOeq[13][3] = 0.001049;
1075 KdYNO[13][3] = 0.000100;
1076 YNOeq[13][4] = 0.001354;
1077 KdYNO[13][4] = 0.000133;
1078 YNOeq[13][5] = 0.001658;
1079 KdYNO[13][5] = 0.000165;
1080 YNOeq[13][6] = 0.001962;
1081 KdYNO[13][6] = 0.000198;
1082 YNOeq[13][7] = 0.002267;
1083 KdYNO[13][7] = 0.000230;
1084 YNOeq[13][8] = 0.002571;
1085 KdYNO[13][8] = 0.000263;
1086 YNOeq[13][9] = 0.002840;
1087 KdYNO[13][9] = 0.000297;
1088 YNOeq[13][10] = 0.003048;
1089 KdYNO[13][10] = 0.000349;
1090 YNOeq[13][11] = 0.003239;
1091 KdYNO[13][11] = 0.000415;
1092 YNOeq[14][0] = 0.000132;
1093 KdYNO[14][0] = 0.000021;
1094 YNOeq[14][1] = 0.000958;
1095 KdYNO[14][1] = 0.000156;
1096 YNOeq[14][2] = 0.001419;
1097 KdYNO[14][2] = 0.000238;
1098 YNOeq[14][3] = 0.001753;
1099 KdYNO[14][3] = 0.000305;
1100 YNOeq[14][4] = 0.002011;
1101 KdYNO[14][4] = 0.000368;
1102 YNOeq[14][5] = 0.002269;
1103 KdYNO[14][5] = 0.000430;
1104 YNOeq[14][6] = 0.002527;
1105 KdYNO[14][6] = 0.000492;
1106 YNOeq[14][7] = 0.002807;
1107 KdYNO[14][7] = 0.000597;
1108 YNOeq[14][8] = 0.003099;
1109 KdYNO[14][8] = 0.000724;
1110 YNOeq[14][9] = 0.003337;
1111 KdYNO[14][9] = 0.000815;
1112 YNOeq[14][10] = 0.003527;
1113 KdYNO[14][10] = 0.000865;
1114 YNOeq[14][11] = 0.003721;
1115 KdYNO[14][11] = 0.000958;
1116 YNOeq[15][0] = 0.000175;
1117 KdYNO[15][0] = 0.000057;
1118 YNOeq[15][1] = 0.000626;
1119 KdYNO[15][1] = 0.000212;
1120 YNOeq[15][2] = 0.001158;
1121 KdYNO[15][2] = 0.000372;
1122 YNOeq[15][3] = 0.001555;
1123 KdYNO[15][3] = 0.000550;
1124 YNOeq[15][4] = 0.001934;
1125 KdYNO[15][4] = 0.000729;
1126 YNOeq[15][5] = 0.002314;
1127 KdYNO[15][5] = 0.000908;
1128 YNOeq[15][6] = 0.002693;
1129 KdYNO[15][6] = 0.001088;
1130 YNOeq[15][7] = 0.003072;
1131 KdYNO[15][7] = 0.001267;
1132 YNOeq[15][8] = 0.003452;
1133 KdYNO[15][8] = 0.001447;
1134 YNOeq[15][9] = 0.003733;
1135 KdYNO[15][9] = 0.001636;
1136 YNOeq[15][10] = 0.003997;
1137 KdYNO[15][10] = 0.001895;
1138 YNOeq[15][11] = 0.004246;
1139 KdYNO[15][11] = 0.002222;
1140 YNOeq[16][0] = 0.000207;
1141 KdYNO[16][0] = 0.000109;
1142 YNOeq[16][1] = 0.001144;
1143 KdYNO[16][1] = 0.000754;
1144 YNOeq[16][2] = 0.001864;
1145 KdYNO[16][2] = 0.001323;
1146 YNOeq[16][3] = 0.002299;
1147 KdYNO[16][3] = 0.001679;
1148 YNOeq[16][4] = 0.002630;
1149 KdYNO[16][4] = 0.001974;
1150 YNOeq[16][5] = 0.002961;
1151 KdYNO[16][5] = 0.002269;
1152 YNOeq[16][6] = 0.003291;
1153 KdYNO[16][6] = 0.002565;
1154 YNOeq[16][7] = 0.003658;
1155 KdYNO[16][7] = 0.003131;
1156 YNOeq[16][8] = 0.004045;
1157 KdYNO[16][8] = 0.003843;
1158 YNOeq[16][9] = 0.004328;
1159 KdYNO[16][9] = 0.004167;
1160 YNOeq[16][10] = 0.004590;
1161 KdYNO[16][10] = 0.004491;
1162 YNOeq[16][11] = 0.004820;
1163 KdYNO[16][11] = 0.004800;
1164 YNOeq[17][0] = 0.000305;
1165 KdYNO[17][0] = 0.000389;
1166 YNOeq[17][1] = 0.000871;
1167 KdYNO[17][1] = 0.001117;
1168 YNOeq[17][2] = 0.001512;
1169 KdYNO[17][2] = 0.001698;
1170 YNOeq[17][3] = 0.002012;
1171 KdYNO[17][3] = 0.002534;
1172 YNOeq[17][4] = 0.002495;
1173 KdYNO[17][4] = 0.003402;
1174 YNOeq[17][5] = 0.002978;
1175 KdYNO[17][5] = 0.004271;
1176 YNOeq[17][6] = 0.003460;
1177 KdYNO[17][6] = 0.005140;
1178 YNOeq[17][7] = 0.003943;
1179 KdYNO[17][7] = 0.006009;
1180 YNOeq[17][8] = 0.004425;
1181 KdYNO[17][8] = 0.006877;
1182 YNOeq[17][9] = 0.004795;
1183 KdYNO[17][9] = 0.007841;
1184 YNOeq[17][10] = 0.005123;
1185 KdYNO[17][10] = 0.008969;
1186 YNOeq[17][11] = 0.005441;
1187 KdYNO[17][11] = 0.010388;
1188 YNOeq[18][0] = 0.000340;
1189 KdYNO[18][0] = 0.000618;
1190 YNOeq[18][1] = 0.001487;
1191 KdYNO[18][1] = 0.003637;
1192 YNOeq[18][2] = 0.002386;
1193 KdYNO[18][2] = 0.006320;
1194 YNOeq[18][3] = 0.002942;
1195 KdYNO[18][3] = 0.007955;
1196 YNOeq[18][4] = 0.003358;
1197 KdYNO[18][4] = 0.009197;
1198 YNOeq[18][5] = 0.003774;
1199 KdYNO[18][5] = 0.010439;
1200 YNOeq[18][6] = 0.004190;
1201 KdYNO[18][6] = 0.011681;
1202 YNOeq[18][7] = 0.004658;
1203 KdYNO[18][7] = 0.014287;
1204 YNOeq[18][8] = 0.005154;
1205 KdYNO[18][8] = 0.017620;
1206 YNOeq[18][9] = 0.005490;
1207 KdYNO[18][9] = 0.018637;
1208 YNOeq[18][10] = 0.005836;
1209 KdYNO[18][10] = 0.020254;
1210 YNOeq[18][11] = 0.006110;
1211 KdYNO[18][11] = 0.021109;
1212 YNOeq[19][0] = 0.000498;
1213 KdYNO[19][0] = 0.002188;
1214 YNOeq[19][1] = 0.001192;
1215 KdYNO[19][1] = 0.005203;
1216 YNOeq[19][2] = 0.001950;
1217 KdYNO[19][2] = 0.007094;
1218 YNOeq[19][3] = 0.002567;
1219 KdYNO[19][3] = 0.010552;
1220 YNOeq[19][4] = 0.003165;
1221 KdYNO[19][4] = 0.014212;
1222 YNOeq[19][5] = 0.003763;
1223 KdYNO[19][5] = 0.017873;
1224 YNOeq[19][6] = 0.004361;
1225 KdYNO[19][6] = 0.021533;
1226 YNOeq[19][7] = 0.004959;
1227 KdYNO[19][7] = 0.025193;
1228 YNOeq[19][8] = 0.005557;
1229 KdYNO[19][8] = 0.028854;
1230 YNOeq[19][9] = 0.006029;
1231 KdYNO[19][9] = 0.033036;
1232 YNOeq[19][10] = 0.006430;
1233 KdYNO[19][10] = 0.037389;
1234 YNOeq[19][11] = 0.006826;
1235 KdYNO[19][11] = 0.042859;
1236 YNOeq[20][0] = 0.000532;
1237 KdYNO[20][0] = 0.003029;
1238 YNOeq[20][1] = 0.001892;
1239 KdYNO[20][1] = 0.015222;
1240 YNOeq[20][2] = 0.002987;
1241 KdYNO[20][2] = 0.026207;
1242 YNOeq[20][3] = 0.003679;
1243 KdYNO[20][3] = 0.032808;
1244 YNOeq[20][4] = 0.004194;
1245 KdYNO[20][4] = 0.037504;
1246 YNOeq[20][5] = 0.004708;
1247 KdYNO[20][5] = 0.042201;
1248 YNOeq[20][6] = 0.005223;
1249 KdYNO[20][6] = 0.046897;
1250 YNOeq[20][7] = 0.005806;
1251 KdYNO[20][7] = 0.057298;
1252 YNOeq[20][8] = 0.006425;
1253 KdYNO[20][8] = 0.070743;
1254 YNOeq[20][9] = 0.006822;
1255 KdYNO[20][9] = 0.073583;
1256 YNOeq[20][10] = 0.007264;
1257 KdYNO[20][10] = 0.080335;
1258 YNOeq[20][11] = 0.007589;
1259 KdYNO[20][11] = 0.082249;
1260 YNOeq[21][0] = 0.000769;
1261 KdYNO[21][0] = 0.010316;
1262 YNOeq[21][1] = 0.001603;
1263 KdYNO[21][1] = 0.021430;
1264 YNOeq[21][2] = 0.002483;
1265 KdYNO[21][2] = 0.027015;
1266 YNOeq[21][3] = 0.003225;
1267 KdYNO[21][3] = 0.039702;
1268 YNOeq[21][4] = 0.003949;
1269 KdYNO[21][4] = 0.053304;
1270 YNOeq[21][5] = 0.004673;
1271 KdYNO[21][5] = 0.066906;
1272 YNOeq[21][6] = 0.005397;
1273 KdYNO[21][6] = 0.080508;
1274 YNOeq[21][7] = 0.006121;
1275 KdYNO[21][7] = 0.094110;
1276 YNOeq[21][8] = 0.006845;
1277 KdYNO[21][8] = 0.107712;
1278 YNOeq[21][9] = 0.007431;
1279 KdYNO[21][9] = 0.123569;
1280 YNOeq[21][10] = 0.007914;
1281 KdYNO[21][10] = 0.138603;
1282 YNOeq[21][11] = 0.008398;
1283 KdYNO[21][11] = 0.157431;
1284 YNOeq[22][0] = 0.000733;
1285 KdYNO[22][0] = 0.008818;
1286 YNOeq[22][1] = 0.002364;
1287 KdYNO[22][1] = 0.055774;
1288 YNOeq[22][2] = 0.003664;
1289 KdYNO[22][2] = 0.095179;
1290 YNOeq[22][3] = 0.004507;
1291 KdYNO[22][3] = 0.118901;
1292 YNOeq[22][4] = 0.005133;
1293 KdYNO[22][4] = 0.135043;
1294 YNOeq[22][5] = 0.005759;
1295 KdYNO[22][5] = 0.151184;
1296 YNOeq[22][6] = 0.006385;
1297 KdYNO[22][6] = 0.167326;
1298 YNOeq[22][7] = 0.007095;
1299 KdYNO[22][7] = 0.203906;
1300 YNOeq[22][8] = 0.007850;
1301 KdYNO[22][8] = 0.251396;
1302 YNOeq[22][9] = 0.008317;
1303 KdYNO[22][9] = 0.258692;
1304 YNOeq[22][10] = 0.008866;
1305 KdYNO[22][10] = 0.283291;
1306 YNOeq[22][11] = 0.009251;
1307 KdYNO[22][11] = 0.286651;
1308 YNOeq[23][0] = 0.001135;
1309 KdYNO[23][0] = 0.041490;
1310 YNOeq[23][1] = 0.002114;
1311 KdYNO[23][1] = 0.078217;
1312 YNOeq[23][2] = 0.003119;
1313 KdYNO[23][2] = 0.093379;
1314 YNOeq[23][3] = 0.003993;
1315 KdYNO[23][3] = 0.135042;
1316 YNOeq[23][4] = 0.004850;
1317 KdYNO[23][4] = 0.180119;
1318 YNOeq[23][5] = 0.005708;
1319 KdYNO[23][5] = 0.225196;
1320 YNOeq[23][6] = 0.006566;
1321 KdYNO[23][6] = 0.270273;
1322 YNOeq[23][7] = 0.007423;
1323 KdYNO[23][7] = 0.315350;
1324 YNOeq[23][8] = 0.008281;
1325 KdYNO[23][8] = 0.360426;
1326 YNOeq[23][9] = 0.008991;
1327 KdYNO[23][9] = 0.413773;
1328 YNOeq[23][10] = 0.009566;
1329 KdYNO[23][10] = 0.460907;
1330 YNOeq[23][11] = 0.010147;
1331 KdYNO[23][11] = 0.519839;
1332 YNOeq[24][0] = 0.001060;
1333 KdYNO[24][0] = 0.032900;
1334 YNOeq[24][1] = 0.002782;
1335 KdYNO[24][1] = 0.147795;
1336 YNOeq[24][2] = 0.003941;
1337 KdYNO[24][2] = 0.167231;
1338 YNOeq[24][3] = 0.005066;
1339 KdYNO[24][3] = 0.258906;
1340 YNOeq[24][4] = 0.005704;
1341 KdYNO[24][4] = 0.264210;
1342 YNOeq[24][5] = 0.006342;
1343 KdYNO[24][5] = 0.269515;
1344 YNOeq[24][6] = 0.007100;
1345 KdYNO[24][6] = 0.328396;
1346 YNOeq[24][7] = 0.007866;
1347 KdYNO[24][7] = 0.390753;
1348 YNOeq[24][8] = 0.008632;
1349 KdYNO[24][8] = 0.453110;
1350 YNOeq[24][9] = 0.009542;
1351 KdYNO[24][9] = 0.635565;
1352 YNOeq[24][10] = 0.010467;
1353 KdYNO[24][10] = 0.829667;
1354 YNOeq[24][11] = 0.011082;
1355 KdYNO[24][11] = 0.901448;
1356 YNOeq[25][0] = 0.001260;
1357 KdYNO[25][0] = 0.060880;
1358 YNOeq[25][1] = 0.003050;
1359 KdYNO[25][1] = 0.229015;
1360 YNOeq[25][2] = 0.004197;
1361 KdYNO[25][2] = 0.310817;
1362 YNOeq[25][3] = 0.005218;
1363 KdYNO[25][3] = 0.375920;
1364 YNOeq[25][4] = 0.006240;
1365 KdYNO[25][4] = 0.441022;
1366 YNOeq[25][5] = 0.007151;
1367 KdYNO[25][5] = 0.533430;
1368 YNOeq[25][6] = 0.008052;
1369 KdYNO[25][6] = 0.717659;
1370 YNOeq[25][7] = 0.008954;
1371 KdYNO[25][7] = 0.901887;
1372 YNOeq[25][8] = 0.009856;
1373 KdYNO[25][8] = 1.086116;
1374 YNOeq[25][9] = 0.010694;
1375 KdYNO[25][9] = 1.250362;
1376 YNOeq[25][10] = 0.011370;
1377 KdYNO[25][10] = 1.385262;
1378 YNOeq[25][11] = 0.012055;
1379 KdYNO[25][11] = 1.553527;
1380 YNOeq[26][0] = 0.001407;
1381 KdYNO[26][0] = 0.092724;
1382 YNOeq[26][1] = 0.002875;
1383 KdYNO[26][1] = 0.212964;
1384 YNOeq[26][2] = 0.004343;
1385 KdYNO[26][2] = 0.333203;
1386 YNOeq[26][3] = 0.005701;
1387 KdYNO[26][3] = 0.555891;
1388 YNOeq[26][4] = 0.006614;
1389 KdYNO[26][4] = 0.700366;
1390 YNOeq[26][5] = 0.007426;
1391 KdYNO[26][5] = 0.738973;
1392 YNOeq[26][6] = 0.008332;
1393 KdYNO[26][6] = 0.913561;
1394 YNOeq[26][7] = 0.009244;
1395 KdYNO[26][7] = 1.096972;
1396 YNOeq[26][8] = 0.010156;
1397 KdYNO[26][8] = 1.280383;
1398 YNOeq[26][9] = 0.011238;
1399 KdYNO[26][9] = 1.810014;
1400 YNOeq[26][10] = 0.012336;
1401 KdYNO[26][10] = 2.373222;
1402 YNOeq[26][11] = 0.013061;
1403 KdYNO[26][11] = 2.573348;
1404 YNOeq[27][0] = 0.001732;
1405 KdYNO[27][0] = 0.191921;
1406 YNOeq[27][1] = 0.003699;
1407 KdYNO[27][1] = 0.638817;
1408 YNOeq[27][2] = 0.004995;
1409 KdYNO[27][2] = 0.858891;
1410 YNOeq[27][3] = 0.006161;
1411 KdYNO[27][3] = 1.035088;
1412 YNOeq[27][4] = 0.007328;
1413 KdYNO[27][4] = 1.211285;
1414 YNOeq[27][5] = 0.008372;
1415 KdYNO[27][5] = 1.458364;
1416 YNOeq[27][6] = 0.009428;
1417 KdYNO[27][6] = 1.965043;
1418 YNOeq[27][7] = 0.010484;
1419 KdYNO[27][7] = 2.471722;
1420 YNOeq[27][8] = 0.011540;
1421 KdYNO[27][8] = 2.978402;
1422 YNOeq[27][9] = 0.012519;
1423 KdYNO[27][9] = 3.428585;
1424 YNOeq[27][10] = 0.013015;
1425 KdYNO[27][10] = 3.273183;
1426 YNOeq[27][11] = 0.014101;
1427 KdYNO[27][11] = 4.229005;
1428 YNOeq[28][0] = 0.001892;
1429 KdYNO[28][0] = 0.271459;
1430 YNOeq[28][1] = 0.003479;
1431 KdYNO[28][1] = 0.548492;
1432 YNOeq[28][2] = 0.005065;
1433 KdYNO[28][2] = 0.825524;
1434 YNOeq[28][3] = 0.006468;
1435 KdYNO[28][3] = 1.347835;
1436 YNOeq[28][4] = 0.007661;
1437 KdYNO[28][4] = 1.792081;
1438 YNOeq[28][5] = 0.008605;
1439 KdYNO[28][5] = 1.891306;
1440 YNOeq[28][6] = 0.009663;
1441 KdYNO[28][6] = 2.357629;
1442 YNOeq[28][7] = 0.010729;
1443 KdYNO[28][7] = 2.847770;
1444 YNOeq[28][8] = 0.011795;
1445 KdYNO[28][8] = 3.337912;
1446 YNOeq[28][9] = 0.012961;
1447 KdYNO[28][9] = 4.311444;
1448 YNOeq[28][10] = 0.014040;
1449 KdYNO[28][10] = 5.344972;
1450 YNOeq[28][11] = 0.015074;
1451 KdYNO[28][11] = 6.473152;
1452 YNOeq[29][0] = 0.002311;
1453 KdYNO[29][0] = 0.547036;
1454 YNOeq[29][1] = 0.004436;
1455 KdYNO[29][1] = 1.623547;
1456 YNOeq[29][2] = 0.005878;
1457 KdYNO[29][2] = 2.167642;
1458 YNOeq[29][3] = 0.007187;
1459 KdYNO[29][3] = 2.608745;
1460 YNOeq[29][4] = 0.008496;
1461 KdYNO[29][4] = 3.049848;
1462 YNOeq[29][5] = 0.009585;
1463 KdYNO[29][5] = 3.545103;
1464 YNOeq[29][6] = 0.010636;
1465 KdYNO[29][6] = 4.567760;
1466 YNOeq[29][7] = 0.011811;
1467 KdYNO[29][7] = 5.626905;
1468 YNOeq[29][8] = 0.012986;
1469 KdYNO[29][8] = 6.686050;
1470 YNOeq[29][9] = 0.014160;
1471 KdYNO[29][9] = 7.745196;
1472 YNOeq[29][10] = 0.014952;
1473 KdYNO[29][10] = 7.970170;
1474 YNOeq[29][11] = 0.015965;
1475 KdYNO[29][11] = 9.430666;
1476 YNOeq[30][0] = 0.002479;
1477 KdYNO[30][0] = 0.729234;
1478 YNOeq[30][1] = 0.004172;
1479 KdYNO[30][1] = 1.328519;
1480 YNOeq[30][2] = 0.005864;
1481 KdYNO[30][2] = 1.927803;
1482 YNOeq[30][3] = 0.007296;
1483 KdYNO[30][3] = 2.736098;
1484 YNOeq[30][4] = 0.008593;
1485 KdYNO[30][4] = 3.615589;
1486 YNOeq[30][5] = 0.009854;
1487 KdYNO[30][5] = 4.471673;
1488 YNOeq[30][6] = 0.011080;
1489 KdYNO[30][6] = 5.647474;
1490 YNOeq[30][7] = 0.012303;
1491 KdYNO[30][7] = 6.844018;
1492 YNOeq[30][8] = 0.013526;
1493 KdYNO[30][8] = 8.040562;
1494 YNOeq[30][9] = 0.014740;
1495 KdYNO[30][9] = 9.826371;
1496 YNOeq[30][10] = 0.016016;
1497 KdYNO[30][10] = 12.394166;
1498 YNOeq[30][11] = 0.015965;
1499 KdYNO[30][11] = 9.430666;
1500 YNOeq[31][0] = 0.002998;
1501 KdYNO[31][0] = 1.423809;
1502 YNOeq[31][1] = 0.004828;
1503 KdYNO[31][1] = 2.519829;
1504 YNOeq[31][2] = 0.006421;
1505 KdYNO[31][2] = 3.574074;
1506 YNOeq[31][3] = 0.007964;
1507 KdYNO[31][3] = 4.619481;
1508 YNOeq[31][4] = 0.009551;
1509 KdYNO[31][4] = 6.218774;
1510 YNOeq[31][5] = 0.011081;
1511 KdYNO[31][5] = 8.172902;
1512 YNOeq[31][6] = 0.012270;
1513 KdYNO[31][6] = 10.093628;
1514 YNOeq[31][7] = 0.013460;
1515 KdYNO[31][7] = 12.014354;
1516 YNOeq[31][8] = 0.014649;
1517 KdYNO[31][8] = 13.935080;
1518 YNOeq[31][9] = 0.015837;
1519 KdYNO[31][9] = 15.873938;
1520 YNOeq[31][10] = 0.016990;
1521 KdYNO[31][10] = 18.181945;
1522 YNOeq[31][11] = 0.015965;
1523 KdYNO[31][11] = 9.430666;
1524 YNOeq[32][0] = 0.002987;
1525 KdYNO[32][0] = 1.478557;
1526 YNOeq[32][1] = 0.004827;
1527 KdYNO[32][1] = 2.797413;
1528 YNOeq[32][2] = 0.006595;
1529 KdYNO[32][2] = 4.249854;
1530 YNOeq[32][3] = 0.008122;
1531 KdYNO[32][3] = 6.151483;
1532 YNOeq[32][4] = 0.009650;
1533 KdYNO[32][4] = 8.053111;
1534 YNOeq[32][5] = 0.011178;
1535 KdYNO[32][5] = 9.954740;
1536 YNOeq[32][6] = 0.012568;
1537 KdYNO[32][6] = 12.597257;
1538 YNOeq[32][7] = 0.013949;
1539 KdYNO[32][7] = 15.287844;
1540 YNOeq[32][8] = 0.015331;
1541 KdYNO[32][8] = 17.978430;
1542 YNOeq[32][9] = 0.016690;
1543 KdYNO[32][9] = 21.877635;
1544 YNOeq[32][10] = 0.016990;
1545 KdYNO[32][10] = 18.181945;
1546 YNOeq[32][11] = 0.015965;
1547 KdYNO[32][11] = 9.430666;
1548 YNOeq[33][0] = 0.003794;
1549 KdYNO[33][0] = 3.423591;
1550 YNOeq[33][1] = 0.005713;
1551 KdYNO[33][1] = 5.641213;
1552 YNOeq[33][2] = 0.007425;
1553 KdYNO[33][2] = 7.845004;
1554 YNOeq[33][3] = 0.009093;
1555 KdYNO[33][3] = 10.045869;
1556 YNOeq[33][4] = 0.010826;
1557 KdYNO[33][4] = 13.462808;
1558 YNOeq[33][5] = 0.012517;
1559 KdYNO[33][5] = 17.651895;
1560 YNOeq[33][6] = 0.013877;
1561 KdYNO[33][6] = 21.407047;
1562 YNOeq[33][7] = 0.015237;
1563 KdYNO[33][7] = 25.162198;
1564 YNOeq[33][8] = 0.016553;
1565 KdYNO[33][8] = 28.970749;
1566 YNOeq[33][9] = 0.017734;
1567 KdYNO[33][9] = 32.850350;
1568 YNOeq[33][10] = 0.016990;
1569 KdYNO[33][10] = 18.181945;
1570 YNOeq[33][11] = 0.015965;
1571 KdYNO[33][11] = 9.430666;
1572 YNOeq[34][0] = 0.003697;
1573 KdYNO[34][0] = 3.302310;
1574 YNOeq[34][1] = 0.005641;
1575 KdYNO[34][1] = 5.913553;
1576 YNOeq[34][2] = 0.007524;
1577 KdYNO[34][2] = 8.840037;
1578 YNOeq[34][3] = 0.009206;
1579 KdYNO[34][3] = 12.826540;
1580 YNOeq[34][4] = 0.010888;
1581 KdYNO[34][4] = 16.813044;
1582 YNOeq[34][5] = 0.012570;
1583 KdYNO[34][5] = 20.799547;
1584 YNOeq[34][6] = 0.014116;
1585 KdYNO[34][6] = 26.307950;
1586 YNOeq[34][7] = 0.015652;
1587 KdYNO[34][7] = 31.915095;
1588 YNOeq[34][8] = 0.017189;
1589 KdYNO[34][8] = 37.522241;
1590 YNOeq[34][9] = 0.017734;
1591 KdYNO[34][9] = 32.850350;
1592 YNOeq[34][10] = 0.016990;
1593 KdYNO[34][10] = 18.181945;
1594 YNOeq[34][11] = 0.015965;
1595 KdYNO[34][11] = 9.430666;
1596 YNOeq[35][0] = 0.004045;
1597 KdYNO[35][0] = 4.577814;
1598 YNOeq[35][1] = 0.006250;
1599 KdYNO[35][1] = 9.798461;
1600 YNOeq[35][2] = 0.008273;
1601 KdYNO[35][2] = 14.855960;
1602 YNOeq[35][3] = 0.010156;
1603 KdYNO[35][3] = 19.786978;
1604 YNOeq[35][4] = 0.012162;
1605 KdYNO[35][4] = 26.201035;
1606 YNOeq[35][5] = 0.013852;
1607 KdYNO[35][5] = 33.581057;
1608 YNOeq[35][6] = 0.015309;
1609 KdYNO[35][6] = 39.701992;
1610 YNOeq[35][7] = 0.016871;
1611 KdYNO[35][7] = 48.099433;
1612 YNOeq[35][8] = 0.018213;
1613 KdYNO[35][8] = 54.505888;
1614 YNOeq[35][9] = 0.017734;
1615 KdYNO[35][9] = 32.850350;
1616 YNOeq[35][10] = 0.016990;
1617 KdYNO[35][10] = 18.181945;
1618 YNOeq[35][11] = 0.015965;
1619 KdYNO[35][11] = 9.430666;
1620 YNOeq[36][0] = 0.004503;
1621 KdYNO[36][0] = 7.015211;
1622 YNOeq[36][1] = 0.006540;
1623 KdYNO[36][1] = 11.921544;
1624 YNOeq[36][2] = 0.008530;
1625 KdYNO[36][2] = 17.499392;
1626 YNOeq[36][3] = 0.010358;
1627 KdYNO[36][3] = 25.335258;
1628 YNOeq[36][4] = 0.012187;
1629 KdYNO[36][4] = 33.171123;
1630 YNOeq[36][5] = 0.014016;
1631 KdYNO[36][5] = 41.006989;
1632 YNOeq[36][6] = 0.015687;
1633 KdYNO[36][6] = 50.802898;
1634 YNOeq[36][7] = 0.017470;
1635 KdYNO[36][7] = 59.199654;
1636 YNOeq[36][8] = 0.018213;
1637 KdYNO[36][8] = 54.505888;
1638 YNOeq[36][9] = 0.017734;
1639 KdYNO[36][9] = 32.850350;
1640 YNOeq[36][10] = 0.016990;
1641 KdYNO[36][10] = 18.181945;
1642 YNOeq[36][11] = 0.015965;
1643 KdYNO[36][11] = 9.430666;
1644 YNOeq[37][0] = 0.004910;
1645 KdYNO[37][0] = 9.582819;
1646 YNOeq[37][1] = 0.007222;
1647 KdYNO[37][1] = 19.290391;
1648 YNOeq[37][2] = 0.009365;
1649 KdYNO[37][2] = 28.791551;
1650 YNOeq[37][3] = 0.011378;
1651 KdYNO[37][3] = 38.132688;
1652 YNOeq[37][4] = 0.013425;
1653 KdYNO[37][4] = 47.851436;
1654 YNOeq[37][5] = 0.015123;
1655 KdYNO[37][5] = 58.733176;
1656 YNOeq[37][6] = 0.016886;
1657 KdYNO[37][6] = 73.864675;
1658 YNOeq[37][7] = 0.018403;
1659 KdYNO[37][7] = 84.540857;
1660 YNOeq[37][8] = 0.018213;
1661 KdYNO[37][8] = 54.505888;
1662 YNOeq[37][9] = 0.017734;
1663 KdYNO[37][9] = 32.850350;
1664 YNOeq[37][10] = 0.016990;
1665 KdYNO[37][10] = 18.181945;
1666 YNOeq[37][11] = 0.015965;
1667 KdYNO[37][11] = 9.430666;
1668 YNOeq[38][0] = 0.005081;
1669 KdYNO[38][0] = 10.810942;
1670 YNOeq[38][1] = 0.007358;
1671 KdYNO[38][1] = 21.064593;
1672 YNOeq[38][2] = 0.009609;
1673 KdYNO[38][2] = 33.021948;
1674 YNOeq[38][3] = 0.011576;
1675 KdYNO[38][3] = 47.544379;
1676 YNOeq[38][4] = 0.013542;
1677 KdYNO[38][4] = 62.066810;
1678 YNOeq[38][5] = 0.015508;
1679 KdYNO[38][5] = 76.589241;
1680 YNOeq[38][6] = 0.017375;
1681 KdYNO[38][6] = 92.178864;
1682 YNOeq[38][7] = 0.018403;
1683 KdYNO[38][7] = 84.540857;
1684 YNOeq[38][8] = 0.018213;
1685 KdYNO[38][8] = 54.505888;
1686 YNOeq[38][9] = 0.017734;
1687 KdYNO[38][9] = 32.850350;
1688 YNOeq[38][10] = 0.016990;
1689 KdYNO[38][10] = 18.181945;
1690 YNOeq[38][11] = 0.015965;
1691 KdYNO[38][11] = 9.430666;
1692 YNOeq[39][0] = 0.005869;
1693 KdYNO[39][0] = 18.928014;
1694 YNOeq[39][1] = 0.008278;
1695 KdYNO[39][1] = 36.111013;
1696 YNOeq[39][2] = 0.010530;
1697 KdYNO[39][2] = 53.043002;
1698 YNOeq[39][3] = 0.012659;
1699 KdYNO[39][3] = 69.780392;
1700 YNOeq[39][4] = 0.014779;
1701 KdYNO[39][4] = 85.511999;
1702 YNOeq[39][5] = 0.016443;
1703 KdYNO[39][5] = 98.841862;
1704 YNOeq[39][6] = 0.018292;
1705 KdYNO[39][6] = 122.897344;
1706 YNOeq[39][7] = 0.018403;
1707 KdYNO[39][7] = 84.540857;
1708 YNOeq[39][8] = 0.018213;
1709 KdYNO[39][8] = 54.505888;
1710 YNOeq[39][9] = 0.017734;
1711 KdYNO[39][9] = 32.850350;
1712 YNOeq[39][10] = 0.016990;
1713 KdYNO[39][10] = 18.181945;
1714 YNOeq[39][11] = 0.015965;
1715 KdYNO[39][11] = 9.430666;
1716 YNOeq[40][0] = 0.006016;
1717 KdYNO[40][0] = 20.651069;
1718 YNOeq[40][1] = 0.008391;
1719 KdYNO[40][1] = 38.537591;
1720 YNOeq[40][2] = 0.010860;
1721 KdYNO[40][2] = 58.400257;
1722 YNOeq[40][3] = 0.013135;
1723 KdYNO[40][3] = 82.521261;
1724 YNOeq[40][4] = 0.015215;
1725 KdYNO[40][4] = 111.901382;
1726 YNOeq[40][5] = 0.017126;
1727 KdYNO[40][5] = 131.273590;
1728 YNOeq[40][6] = 0.018292;
1729 KdYNO[40][6] = 122.897344;
1730 YNOeq[40][7] = 0.018403;
1731 KdYNO[40][7] = 84.540857;
1732 YNOeq[40][8] = 0.018213;
1733 KdYNO[40][8] = 54.505888;
1734 YNOeq[40][9] = 0.017734;
1735 KdYNO[40][9] = 32.850350;
1736 YNOeq[40][10] = 0.016990;
1737 KdYNO[40][10] = 18.181945;
1738 YNOeq[40][11] = 0.015965;
1739 KdYNO[40][11] = 9.430666;
1740 YNOeq[41][0] = 0.006678;
1741 KdYNO[41][0] = 30.590106;
1742 YNOeq[41][1] = 0.009003;
1743 KdYNO[41][1] = 56.851468;
1744 YNOeq[41][2] = 0.011337;
1745 KdYNO[41][2] = 85.428380;
1746 YNOeq[41][3] = 0.013641;
1747 KdYNO[41][3] = 112.821628;
1748 YNOeq[41][4] = 0.015914;
1749 KdYNO[41][4] = 139.021404;
1750 YNOeq[41][5] = 0.017822;
1751 KdYNO[41][5] = 162.375287;
1752 YNOeq[41][6] = 0.018292;
1753 KdYNO[41][6] = 122.897344;
1754 YNOeq[41][7] = 0.018403;
1755 KdYNO[41][7] = 84.540857;
1756 YNOeq[41][8] = 0.018213;
1757 KdYNO[41][8] = 54.505888;
1758 YNOeq[41][9] = 0.017734;
1759 KdYNO[41][9] = 32.850350;
1760 YNOeq[41][10] = 0.016990;
1761 KdYNO[41][10] = 18.181945;
1762 YNOeq[41][11] = 0.015965;
1763 KdYNO[41][11] = 9.430666;
1764 YNOeq[42][0] = 0.007042;
1765 KdYNO[42][0] = 37.697174;
1766 YNOeq[42][1] = 0.009506;
1767 KdYNO[42][1] = 67.480598;
1768 YNOeq[42][2] = 0.012090;
1769 KdYNO[42][2] = 100.891368;
1770 YNOeq[42][3] = 0.014310;
1771 KdYNO[42][3] = 126.744041;
1772 YNOeq[42][4] = 0.016455;
1773 KdYNO[42][4] = 166.893523;
1774 YNOeq[42][5] = 0.017822;
1775 KdYNO[42][5] = 162.375287;
1776 YNOeq[42][6] = 0.018292;
1777 KdYNO[42][6] = 122.897344;
1778 YNOeq[42][7] = 0.018403;
1779 KdYNO[42][7] = 84.540857;
1780 YNOeq[42][8] = 0.018213;
1781 KdYNO[42][8] = 54.505888;
1782 YNOeq[42][9] = 0.017734;
1783 KdYNO[42][9] = 32.850350;
1784 YNOeq[42][10] = 0.016990;
1785 KdYNO[42][10] = 18.181945;
1786 YNOeq[42][11] = 0.015965;
1787 KdYNO[42][11] = 9.430666;
1788 YNOeq[43][0] = 0.007797;
1789 KdYNO[43][0] = 54.697991;
1790 YNOeq[43][1] = 0.010132;
1791 KdYNO[43][1] = 94.655766;
1792 YNOeq[43][2] = 0.012471;
1793 KdYNO[43][2] = 137.909658;
1794 YNOeq[43][3] = 0.014842;
1795 KdYNO[43][3] = 175.933649;
1796 YNOeq[43][4] = 0.017201;
1797 KdYNO[43][4] = 214.977433;
1798 YNOeq[43][5] = 0.017822;
1799 KdYNO[43][5] = 162.375287;
1800 YNOeq[43][6] = 0.018292;
1801 KdYNO[43][6] = 122.897344;
1802 YNOeq[43][7] = 0.018403;
1803 KdYNO[43][7] = 84.540857;
1804 YNOeq[43][8] = 0.018213;
1805 KdYNO[43][8] = 54.505888;
1806 YNOeq[43][9] = 0.017734;
1807 KdYNO[43][9] = 32.850350;
1808 YNOeq[43][10] = 0.016990;
1809 KdYNO[43][10] = 18.181945;
1810 YNOeq[43][11] = 0.015965;
1811 KdYNO[43][11] = 9.430666;
1812 YNOeq[44][0] = 0.008082;
1813 KdYNO[44][0] = 62.810606;
1814 YNOeq[44][1] = 0.010710;
1815 KdYNO[44][1] = 114.217438;
1816 YNOeq[44][2] = 0.013385;
1817 KdYNO[44][2] = 168.602766;
1818 YNOeq[44][3] = 0.015685;
1819 KdYNO[44][3] = 210.612953;
1820 YNOeq[44][4] = 0.017201;
1821 KdYNO[44][4] = 214.977433;
1822 YNOeq[44][5] = 0.017822;
1823 KdYNO[44][5] = 162.375287;
1824 YNOeq[44][6] = 0.018292;
1825 KdYNO[44][6] = 122.897344;
1826 YNOeq[44][7] = 0.018403;
1827 KdYNO[44][7] = 84.540857;
1828 YNOeq[44][8] = 0.018213;
1829 KdYNO[44][8] = 54.505888;
1830 YNOeq[44][9] = 0.017734;
1831 KdYNO[44][9] = 32.850350;
1832 YNOeq[44][10] = 0.016990;
1833 KdYNO[44][10] = 18.181945;
1834 YNOeq[44][11] = 0.015965;
1835 KdYNO[44][11] = 9.430666;
1836 YNOeq[45][0] = 0.008989;
1837 KdYNO[45][0] = 94.851480;
1838 YNOeq[45][1] = 0.011436;
1839 KdYNO[45][1] = 142.583015;
1840 YNOeq[45][2] = 0.013888;
1841 KdYNO[45][2] = 197.204589;
1842 YNOeq[45][3] = 0.016267;
1843 KdYNO[45][3] = 262.823904;
1844 YNOeq[45][4] = 0.017201;
1845 KdYNO[45][4] = 214.977433;
1846 YNOeq[45][5] = 0.017822;
1847 KdYNO[45][5] = 162.375287;
1848 YNOeq[45][6] = 0.018292;
1849 KdYNO[45][6] = 122.897344;
1850 YNOeq[45][7] = 0.018403;
1851 KdYNO[45][7] = 84.540857;
1852 YNOeq[45][8] = 0.018213;
1853 KdYNO[45][8] = 54.505888;
1854 YNOeq[45][9] = 0.017734;
1855 KdYNO[45][9] = 32.850350;
1856 YNOeq[45][10] = 0.016990;
1857 KdYNO[45][10] = 18.181945;
1858 YNOeq[45][11] = 0.015965;
1859 KdYNO[45][11] = 9.430666;
1860 YNOeq[46][0] = 0.009261;
1861 KdYNO[46][0] = 107.006637;
1862 YNOeq[46][1] = 0.011927;
1863 KdYNO[46][1] = 185.033724;
1864 YNOeq[46][2] = 0.014583;
1865 KdYNO[46][2] = 260.495804;
1866 YNOeq[46][3] = 0.016267;
1867 KdYNO[46][3] = 262.823904;
1868 YNOeq[46][4] = 0.017201;
1869 KdYNO[46][4] = 214.977433;
1870 YNOeq[46][5] = 0.017822;
1871 KdYNO[46][5] = 162.375287;
1872 YNOeq[46][6] = 0.018292;
1873 KdYNO[46][6] = 122.897344;
1874 YNOeq[46][7] = 0.018403;
1875 KdYNO[46][7] = 84.540857;
1876 YNOeq[46][8] = 0.018213;
1877 KdYNO[46][8] = 54.505888;
1878 YNOeq[46][9] = 0.017734;
1879 KdYNO[46][9] = 32.850350;
1880 YNOeq[46][10] = 0.016990;
1881 KdYNO[46][10] = 18.181945;
1882 YNOeq[46][11] = 0.015965;
1883 KdYNO[46][11] = 9.430666;
1884 YNOeq[47][0] = 0.010112;
1885 KdYNO[47][0] = 150.669136;
1886 YNOeq[47][1] = 0.012741;
1887 KdYNO[47][1] = 233.387507;
1888 YNOeq[47][2] = 0.015179;
1889 KdYNO[47][2] = 305.763808;
1890 YNOeq[47][3] = 0.016267;
1891 KdYNO[47][3] = 262.823904;
1892 YNOeq[47][4] = 0.017201;
1893 KdYNO[47][4] = 214.977433;
1894 YNOeq[47][5] = 0.017822;
1895 KdYNO[47][5] = 162.375287;
1896 YNOeq[47][6] = 0.018292;
1897 KdYNO[47][6] = 122.897344;
1898 YNOeq[47][7] = 0.018403;
1899 KdYNO[47][7] = 84.540857;
1900 YNOeq[47][8] = 0.018213;
1901 KdYNO[47][8] = 54.505888;
1902 YNOeq[47][9] = 0.017734;
1903 KdYNO[47][9] = 32.850350;
1904 YNOeq[47][10] = 0.016990;
1905 KdYNO[47][10] = 18.181945;
1906 YNOeq[47][11] = 0.015965;
1907 KdYNO[47][11] = 9.430666;
1908 YNOeq[48][0] = 0.010506;
1909 KdYNO[48][0] = 177.128602;
1910 YNOeq[48][1] = 0.013232;
1911 KdYNO[48][1] = 270.988925;
1912 YNOeq[48][2] = 0.015179;
1913 KdYNO[48][2] = 305.763808;
1914 YNOeq[48][3] = 0.016267;
1915 KdYNO[48][3] = 262.823904;
1916 YNOeq[48][4] = 0.017201;
1917 KdYNO[48][4] = 214.977433;
1918 YNOeq[48][5] = 0.017822;
1919 KdYNO[48][5] = 162.375287;
1920 YNOeq[48][6] = 0.018292;
1921 KdYNO[48][6] = 122.897344;
1922 YNOeq[48][7] = 0.018403;
1923 KdYNO[48][7] = 84.540857;
1924 YNOeq[48][8] = 0.018213;
1925 KdYNO[48][8] = 54.505888;
1926 YNOeq[48][9] = 0.017734;
1927 KdYNO[48][9] = 32.850350;
1928 YNOeq[48][10] = 0.016990;
1929 KdYNO[48][10] = 18.181945;
1930 YNOeq[48][11] = 0.015965;
1931 KdYNO[48][11] = 9.430666;
1932 YNOeq[49][0] = 0.011414;
1933 KdYNO[49][0] = 247.150173;
1934 YNOeq[49][1] = 0.013859;
1935 KdYNO[49][1] = 270.988925;
1936 YNOeq[49][2] = 0.015179;
1937 KdYNO[49][2] = 305.763808;
1938 YNOeq[49][3] = 0.016267;
1939 KdYNO[49][3] = 262.823904;
1940 YNOeq[49][4] = 0.017201;
1941 KdYNO[49][4] = 214.977433;
1942 YNOeq[49][5] = 0.017822;
1943 KdYNO[49][5] = 162.375287;
1944 YNOeq[49][6] = 0.018292;
1945 KdYNO[49][6] = 122.897344;
1946 YNOeq[49][7] = 0.018403;
1947 KdYNO[49][7] = 84.540857;
1948 YNOeq[49][8] = 0.018213;
1949 KdYNO[49][8] = 54.505888;
1950 YNOeq[49][9] = 0.017734;
1951 KdYNO[49][9] = 32.850350;
1952 YNOeq[49][10] = 0.016990;
1953 KdYNO[49][10] = 18.181945;
1954 YNOeq[49][11] = 0.015965;
1955 KdYNO[49][11] = 9.430666;
1956 YNOeq[50][0] = 0.011781;
1957 KdYNO[50][0] = 247.150173;
1958 YNOeq[50][1] = 0.013859;
1959 KdYNO[50][1] = 270.988925;
1960 YNOeq[50][2] = 0.015179;
1961 KdYNO[50][2] = 305.763808;
1962 YNOeq[50][3] = 0.016267;
1963 KdYNO[50][3] = 262.823904;
1964 YNOeq[50][4] = 0.017201;
1965 KdYNO[50][4] = 214.977433;
1966 YNOeq[50][5] = 0.017822;
1967 KdYNO[50][5] = 162.375287;
1968 YNOeq[50][6] = 0.018292;
1969 KdYNO[50][6] = 122.897344;
1970 YNOeq[50][7] = 0.018403;
1971 KdYNO[50][7] = 84.540857;
1972 YNOeq[50][8] = 0.018213;
1973 KdYNO[50][8] = 54.505888;
1974 YNOeq[50][9] = 0.017734;
1975 KdYNO[50][9] = 32.850350;
1976 YNOeq[50][10] = 0.016990;
1977 KdYNO[50][10] = 18.181945;
1978 YNOeq[50][11] = 0.015965;
1979 KdYNO[50][11] = 9.430666;
1980 YNOeq[51][0] = 0.012459;
1981 KdYNO[51][0] = 247.150173;
1982 YNOeq[51][1] = 0.013859;
1983 KdYNO[51][1] = 270.988925;
1984 YNOeq[51][2] = 0.015179;
1985 KdYNO[51][2] = 305.763808;
1986 YNOeq[51][3] = 0.016267;
1987 KdYNO[51][3] = 262.823904;
1988 YNOeq[51][4] = 0.017201;
1989 KdYNO[51][4] = 214.977433;
1990 YNOeq[51][5] = 0.017822;
1991 KdYNO[51][5] = 162.375287;
1992 YNOeq[51][6] = 0.018292;
1993 KdYNO[51][6] = 122.897344;
1994 YNOeq[51][7] = 0.018403;
1995 KdYNO[51][7] = 84.540857;
1996 YNOeq[51][8] = 0.018213;
1997 KdYNO[51][8] = 54.505888;
1998 YNOeq[51][9] = 0.017734;
1999 KdYNO[51][9] = 32.850350;
2000 YNOeq[51][10] = 0.016990;
2001 KdYNO[51][10] = 18.181945;
2002 YNOeq[51][11] = 0.015965;
2003 KdYNO[51][11] = 9.430666;
2005 if(temperature < 1250) {
2007 }
else if(temperature > 3750) {
2010 i = (temperature - 1200) / 50;
2012 if(mO2 / mtotal <= 0.11) {
2013 j = (mO2 / mtotal + 0.01) / 0.01;
2014 di = (temperature - 1200) / 50 - i;
2015 dj = (mO2 / mtotal + 0.01) / 0.01 - j;
2019 j = (0.10 + 0.01) / 0.01;
2020 di = (temperature - 1200) / 50 - i;
2021 dj = (0.10 + 0.01) / 0.01 - j;
2040 *YNOeq_value = (1 / dj * (1 / di * YNOeq[i][j] + 1 / (1 - di) * YNOeq[i + 1][j]) / (1 / di + 1 / (1 - di)) + 1 /
2041 (1 - dj) * (1 / di * YNOeq[i][j + 1] + 1 / (1 - di) * YNOeq[i + 1][j + 1]) / (1 / di + 1 / (1 - di))) / (1 / dj + 1 /
2044 *KdYNO_value = (1 / dj * (1 / di * KdYNO[i][j] + 1 / (1 - di) * KdYNO[i + 1][j]) / (1 / di + 1 / (1 - di)) + 1 /
2045 (1 - dj) * (1 / di * KdYNO[i][j + 1] + 1 / (1 - di) * KdYNO[i + 1][j + 1]) / (1 / di + 1 / (1 - di))) / (1 / dj + 1 /
2050 void FUNCTION_SOOT_C(
double *soot_pre,
double element_FI) {
2054 double *soot_pre_vector;
2057 soot_pre_vector = (
double*) malloc(11 *
sizeof(
double));
2058 FI_vector = (
double*) malloc(11 *
sizeof(
double));
2060 soot_pre_vector[0] = 0;
2061 soot_pre_vector[1] = 0.01;
2062 soot_pre_vector[2] = 0.02;
2063 soot_pre_vector[3] = 0.05;
2064 soot_pre_vector[4] = 0.1;
2065 soot_pre_vector[5] = 0.2;
2066 soot_pre_vector[6] = 0.6;
2067 soot_pre_vector[7] = 0.7;
2068 soot_pre_vector[8] = 0.8;
2069 soot_pre_vector[9] = 0.9;
2070 soot_pre_vector[10] = 1;
2084 if(element_FI >= 10) {
2088 for(i = 0; i < 10; i++) {
2089 if((dk == 0) && (element_FI >= FI_vector[i])) {
2094 *soot_pre = (element_FI - FI_vector[di]) * (soot_pre_vector[di + 1] - soot_pre_vector[di]) /
2095 (FI_vector[di + 1] - FI_vector[di]) + soot_pre_vector[di];
2099 free(soot_pre_vector);
2113 void ACT(
double *engine_parameters,
double *engine_model_constants,
double *test_variables,
double *injection_rate,
2114 double *CAD_injection_rate,
int size_inlet_inj,
int NIN,
double *SOI,
double *EOI,
int CAI,
double *CAD_exit,
2115 double *HRF_exit,
double *ROHR_exit,
double *p_cyl_exit,
double *dp_da_cyl_exit,
double *T_cyl_exit,
2116 double *H_cooler_exit,
double *mean_var_exit,
double *heat_transfer,
double *injection_rate_exit,
2117 double *accum_injection_rate_exit,
sINtype dataIN,
sOUTtype *dataOUT) {
2121 stLiquid = 0, stEvaporated = 1, stOvermixed = 2, stBurned_poor_premix = 3, stBurned_rich_premix = 4, stBurned_by_diffusion = 5, stBurned_by_second_diffusion = 6, stInto_diffusion_flame = 7
2124 int elem_i = 0, pulso_i = 0;
2126 double realelementmfreac = 0.;
2129 double PI = __cons::Pi;
2130 double Runiv = 8314.;
2131 FILE *fich, *foculto, *foculto2;
2134 int CalcRad[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
2137 double YO2aire = 0.23019;
2139 double Zst = 0., Z1 = 0., Z0 = 0., TCC = 0., TSC = 0.;
2147 double Piston_D = engine_parameters[0];
2148 double S = engine_parameters[1];
2149 double Crank_L = S / 2.;
2150 double Connecting_Rod_L = engine_parameters[3];
2152 double E = engine_parameters[4];
2153 double Piston_Axis_D = engine_parameters[5];
2154 double Piston_Crown_H = engine_parameters[6];
2155 double DBowl = engine_parameters[7];
2156 double VBowl = engine_parameters[8];
2157 double M_Connecting_Rod = engine_parameters[9];
2159 double M_P_R_PA = engine_parameters[10];
2161 double C_ESteel = engine_parameters[11];
2163 double C_Mech_Defor = engine_parameters[12];
2165 double C_MBLBY = engine_parameters[13];
2167 double GCRatio = engine_parameters[14];
2168 double n_holes = engine_parameters[15];
2169 double nozzle_d = engine_parameters[16];
2170 double dc = engine_parameters[17];
2171 double CTM = engine_parameters[18];
2172 double WC1A = engine_parameters[19];
2173 double WC1B = engine_parameters[20];
2174 double C2 = engine_parameters[21];
2176 double IVC = engine_parameters[22];
2177 double EVO = engine_parameters[23];
2179 double Kmixture1 = engine_parameters[24];
2183 if(Kmixture1 < 0.) {
2185 Kmixture1 = -Kmixture1;
2191 double Kmixture2 = engine_model_constants[0];
2193 double Kpmx1 = engine_model_constants[1];
2195 double Kpmx2 = engine_model_constants[2];
2196 double Kpmx3 = engine_model_constants[3];
2197 double Kpmx4 = engine_model_constants[4];
2198 double Kpmx5 = engine_model_constants[5];
2200 double KID1 = engine_model_constants[6];
2202 double KID2 = engine_model_constants[7];
2203 double KID3 = engine_model_constants[8];
2204 double KID4 = engine_model_constants[9];
2205 double KID5 = engine_model_constants[10];
2207 double KNOx1 = engine_model_constants[11];
2209 double KNOx2 = engine_model_constants[12];
2211 double EC = engine_model_constants[13];
2213 double KSOOTA1 = engine_model_constants[14];
2215 double KSOOTA2 = engine_model_constants[15];
2216 double KSOOTA3 = engine_model_constants[16];
2218 double KSOOTA4 = engine_model_constants[17];
2219 double KSOOTA5 = engine_model_constants[18];
2220 double KSOOTA6 = engine_model_constants[19];
2221 double KSOOTA7 = engine_model_constants[20];
2223 double KSOOTC1 = dataIN.KSOOTC1;
2230 double SOOT_EVO_A = 0.;
2231 double SOOT_EVO_C = 0.;
2237 double speed = test_variables[0];
2238 double mairadm = test_variables[1];
2240 double mairIVC = test_variables[2];
2242 double temperatureIVC = test_variables[3];
2244 double mfuel = test_variables[4];
2245 double PCR = test_variables[5];
2246 double inlet_pres = test_variables[6];
2247 double exhaust_pres = test_variables[7];
2248 double Cbb = test_variables[8];
2249 double Atmosphere_press = test_variables[9];
2250 double inj_fuel_temp = test_variables[10];
2252 double TCYL_HEAD = test_variables[11];
2253 double TCYL = test_variables[12];
2254 double TPIS = test_variables[13];
2256 double rateBG = 1. - mairadm / mairIVC;
2260 double rofuel = 830.;
2263 double mmfuel = 12. + HC;
2264 double Vc_factor = 0.99;
2266 double MW_air = 28.97;
2267 double MW_fuel = 152.2;
2268 double MW_burned = 29.13;
2269 double HP = 42920000.;
2270 double T_Evaporation_fuel = 489.35;
2275 D = D * speed / 2000.;
2280 double PRECISION_ITERATION = 0.000005;
2281 double delta_t = 0.;
2283 double mf_burned = 0.;
2284 double mf_burned_pmx = 0.;
2285 double Cylinder_capacity = 0.;
2296 double dp_da_max = 0;
2298 double YO2_bowl_i = 0.;
2303 double *time_vector;
2304 double *time_vector_exit;
2306 double *complete_CAD;
2309 double *virtual_velocity;
2310 double *inj_velocity;
2313 double *complete_p_cyl;
2334 double *complete_V_cyl;
2336 double *complete_deform;
2344 double *evol_Radiacion;
2350 double *mixture_correction;
2352 int inj_counter = 0;
2353 double rate_area = 0.;
2355 int i_aux = 0, j_aux = 0;
2361 double **mfuel_i_IM;
2362 double **mfuel_ij_IM;
2367 double YNOeq_value = 0.;
2368 double KdYNO_value = 0.;
2375 double soot_pre = 0.;
2376 double RES_FSN = 0.;
2380 struct stPropertiesElement {
2383 double dmtotal = 0.;
2390 double mf_reac = 0.;
2391 double mf_evap = 0.;
2397 double Rpmx_value = 0.;
2402 double mSOOT_A = 0.;
2403 double dSOOT_A = 0.;
2404 double mSOOT_B = 0.;
2405 double dSOOT_B = 0.;
2406 double mSOOT_C = 0.;
2407 double dSOOT_C = 0.;
2410 double Pcyl_POC = 0.;
2413 double soot_precursor = 0.;
2423 struct stRealElementComposition {
2431 double mf_reac = 0.;
2433 double dmf_reac = 0.;
2434 double mf_reac_pmx = 0.;
2443 double mtotal_IVC = 0.;
2444 double mO2_IVC = 0.;
2446 double mN2_IVC = 0.;
2448 double mCO2_IVC = 0.;
2450 double mH2O_IVC = 0.;
2452 double NOx_IVC = 0.;
2454 double mSOOT_IVC_B = 0.;
2456 double mSOOT_IVC_A = 0.;
2457 double mSOOT_IVC_C = 0.;
2462 double YCO2IVC = 0.;
2463 double YH2OIVC = 0.;
2467 double *mtotal_bowl;
2478 double *mSOOT_bowl_A;
2479 double *mSOOT_bowl_B;
2480 double *mSOOT_bowl_C;
2482 double mSOOT_bowl_A_i_burned = 0.;
2483 double mSOOT_bowl_B_i_burned = 0.;
2484 double mSOOT_bowl_C_i_burned = 0.;
2498 double *dmtotal_Gfactor;
2511 int complete_size = 0;
2513 int complete_prev_size = 0;
2515 int complete_post_size = 0;
2517 double auxiliar = 0.;
2518 double *vector_to_interpolate;
2520 double counter_CAD_1 = 0.;
2521 double counter_CAD_2 = 0.;
2523 int m = 0, i = 0, j = 0, aux = 0;
2524 double Y1 = 0., T1 = 0., RID1 = 0.;
2526 double maximum = 0.;
2527 double minimum = 0.;
2528 double aux_mfuel = 0., Tot = 0.;
2530 int element_value = 0;
2546 double Ang_Grab = 0.;
2549 if(RadCalc == 1 && Ang_Grab > -179.) {
2550 foculto = fopen(
"elementos.csv",
"w");
2554 foculto2 = fopen(
"data_eje.csv",
"w");
2555 if(foculto2 == NULL)
2561 delta_t = 2. * 1e-5;
2563 auxiliar = IVC * 60. / (360. * speed);
2566 auxiliar = auxiliar + delta_t;
2567 }
while(auxiliar <= EVO * 60. / (360. * speed));
2569 CAD = (
double*) malloc(size *
sizeof(
double));
2570 time_vector = (
double*) malloc(size *
sizeof(
double));
2572 for(counter = 0; counter < size; counter++) {
2574 time_vector[counter] = IVC * 60. / (360. * speed);
2577 time_vector[counter] = time_vector[counter - 1] + delta_t;
2578 CAD[counter] = CAD[counter - 1] + delta_t * 360. * speed / 60.;
2584 SOI_IM = (
double*) malloc(8 *
sizeof(
double));
2585 EOI_IM = (
double*) malloc(8 *
sizeof(
double));
2586 SOC_IM = (
double*) malloc(8 *
sizeof(
double));
2588 mtotal_bowl = (
double*) malloc(size *
sizeof(
double));
2589 mO2_bowl = (
double*) malloc(size *
sizeof(
double));
2590 mN2_bowl = (
double*) malloc(size *
sizeof(
double));
2591 mCO2_bowl = (
double*) malloc(size *
sizeof(
double));
2592 mH2O_bowl = (
double*) malloc(size *
sizeof(
double));
2593 mHC_bowl = (
double*) malloc(size *
sizeof(
double));
2594 Tsq_cyl = (
double*) malloc(size *
sizeof(
double));
2595 Tadib_cyl = (
double*) malloc(size *
sizeof(
double));
2596 mNOx_bowl = (
double*) malloc(size *
sizeof(
double));
2597 mCO_bowl = (
double*) malloc(size *
sizeof(
double));
2598 HC_bowl = (
double*) malloc(size *
sizeof(
double));
2599 mSOOT_bowl_A = (
double*) malloc(size *
sizeof(
double));
2600 mSOOT_bowl_B = (
double*) malloc(size *
sizeof(
double));
2601 mSOOT_bowl_C = (
double*) malloc(size *
sizeof(
double));
2603 mtotal_Vc = (
double*) malloc(size *
sizeof(
double));
2604 mO2_Vc = (
double*) malloc(size *
sizeof(
double));
2605 mN2_Vc = (
double*) malloc(size *
sizeof(
double));
2606 mCO2_Vc = (
double*) malloc(size *
sizeof(
double));
2607 mH2O_Vc = (
double*) malloc(size *
sizeof(
double));
2608 mNOx_Vc = (
double*) malloc(size *
sizeof(
double));
2609 mSOOT_Vc_A = (
double*) malloc(size *
sizeof(
double));
2610 mSOOT_Vc_B = (
double*) malloc(size *
sizeof(
double));
2611 mSOOT_Vc_C = (
double*) malloc(size *
sizeof(
double));
2612 mHC_Vc = (
double*) malloc(size *
sizeof(
double));
2613 dmtotal_Gfactor = (
double*) malloc(size *
sizeof(
double));
2615 dmf = (
double*) malloc(size *
sizeof(
double));
2616 acu_dmf = (
double*) malloc(size *
sizeof(
double));
2617 inj_velocity = (
double*) malloc(size *
sizeof(
double));
2618 virtual_velocity = (
double*) malloc(size *
sizeof(
double));
2620 HRF = (
double*) malloc(size *
sizeof(
double));
2621 HRF_PMX = (
double*) malloc(size *
sizeof(
double));
2622 ROHR = (
double*) malloc(size *
sizeof(
double));
2623 T_cyl = (
double*) malloc(size *
sizeof(
double));
2624 acu_Mbb = (
double*) malloc(size *
sizeof(
double));
2625 V_cyl = (
double*) malloc(size *
sizeof(
double));
2626 Mbb = (
double*) malloc(size *
sizeof(
double));
2627 Yair = (
double*) malloc(size *
sizeof(
double));
2628 Yburned = (
double*) malloc(size *
sizeof(
double));
2629 Yfuel = (
double*) malloc(size *
sizeof(
double));
2630 U = (
double*) malloc(size *
sizeof(
double));
2631 CV = (
double*) malloc(size *
sizeof(
double));
2632 H_cooler = (
double*) malloc(size *
sizeof(
double));
2633 Qcylhead = (
double*) malloc(size *
sizeof(
double));
2634 Qcyl = (
double*) malloc(size *
sizeof(
double));
2635 Qpis = (
double*) malloc(size *
sizeof(
double));
2636 H = (
double*) malloc(size *
sizeof(
double));
2637 defor = (
double*) malloc(size *
sizeof(
double));
2638 Rmixture = (
double*) malloc(size *
sizeof(
double));
2639 roair = (
double*) malloc(size *
sizeof(
double));
2640 p_cyl = (
double*) malloc(size *
sizeof(
double));
2641 dp_da_cyl = (
double*) malloc(size *
sizeof(
double));
2642 XLO = (
double*) malloc(size *
sizeof(
double));
2643 PEN_min = (
double*) malloc(size *
sizeof(
double));
2644 PEN_max = (
double*) malloc(size *
sizeof(
double));
2645 evol_mSoot = (
double*) malloc(size *
sizeof(
double));
2646 mSootCil = (
double*) malloc(size *
sizeof(
double));
2647 evol_Radiacion = (
double*) malloc(size *
sizeof(
double));
2649 aux_vector = (
double*) malloc(CAI *
sizeof(
double));
2653 YNOeq = (
double**) malloc(52 *
sizeof(
double*));
2654 for(counter = 0; counter < 52; counter++) {
2655 YNOeq[counter] = (
double*) malloc(12 *
sizeof(
double));
2658 KdYNO = (
double**) malloc(52 *
sizeof(
double*));
2659 for(counter = 0; counter < 52; counter++) {
2660 KdYNO[counter] = (
double*) malloc(12 *
sizeof(
double));
2665 for(counter = 0; counter < NIN; counter++) {
2666 SOI_IM[counter] = SOI[counter];
2667 EOI_IM[counter] = EOI[counter];
2668 SOC_IM[counter] = EVO;
2670 vector_to_interpolate = (
double*) malloc(size_inlet_inj *
sizeof(
double));
2671 for(counter = 0; counter < size_inlet_inj; counter++) {
2672 vector_to_interpolate[counter] = injection_rate[counter];
2675 FUNCTION_FOR_INTERPOLATION(dmf, time_vector, CAD_injection_rate, vector_to_interpolate, size, size_inlet_inj, speed);
2676 free(vector_to_interpolate);
2678 for(counter = 0; counter < size; counter++) {
2679 if(time_vector[counter] <= SOI_IM[0] * 60. / (360. * speed)) {
2682 if(time_vector[counter] >= EOI_IM[NIN - 1] * 60. / (360. * speed)) {
2685 for(inj_counter = 1; inj_counter < NIN; inj_counter++) {
2686 if((time_vector[counter] >= EOI_IM[inj_counter - 1] * 60. / (360. * speed))
2687 && (time_vector[counter] <= SOI_IM[inj_counter] * 60. / (360. * speed))) {
2692 if(dmf[counter] < 0) {
2702 for(counter = 0; counter < size - 1; counter++) {
2703 rate_area = rate_area + (time_vector[counter + 1] - time_vector[counter]) * min(dmf[counter],
2704 dmf[counter + 1]) + (time_vector[counter + 1] - time_vector[counter]) * pow(
2705 pow((dmf[counter + 1] - dmf[counter]) / 2, 2), 0.5);
2709 for(counter = 0; counter < size; counter++) {
2710 dmf[counter] = dmf[counter] * mfuel * 1e-6 / rate_area;
2718 CALCULUS_OF_ACCUMULATED_INJ_RATE(acu_dmf, dmf, time_vector, size);
2722 CALCULUS_OF_VIRTUAL_VELOCITY(inj_velocity, virtual_velocity, dmf, time_vector, rofuel, dc, n_holes, nozzle_d, D, size,
2723 PI, speed, EOI_IM, inj_num, SOI_IM, Piston_D, DBowl, CTM, CAD, Kswirl);
2727 STOICHIOMETRY_CONSTANTS(HC, &Kst1, &Kst2, &Kst3, &Kst4, &Kst5, &Kst6);
2731 num_i_IM = (
int*) malloc(inj_num *
sizeof(
int));
2736 CALCULUS_OF_NUMBER_ELEMENTS(num_i_IM, time_vector, size, speed, SOI_IM, EOI_IM, inj_num);
2739 num_i = num_i_IM[0];
2741 POI_IM = (
double**) malloc(inj_num *
sizeof(
double*));
2742 for(counter = 0; counter < inj_num; counter++) {
2743 num_i = num_i_IM[counter];
2744 POI_IM[counter] = (
double*) malloc(num_i *
sizeof(
double));
2747 POC_IM = (
double**) malloc(inj_num *
sizeof(
double*));
2748 for(counter = 0; counter < inj_num; counter++) {
2749 element_value = num_i_IM[counter] * num_j;
2750 POC_IM[counter] = (
double*) malloc(element_value *
sizeof(
double));
2753 mfuel_i_IM = (
double**) malloc(inj_num *
sizeof(
double*));
2754 for(counter = 0; counter < inj_num; counter++) {
2755 num_i = num_i_IM[counter];
2756 mfuel_i_IM[counter] = (
double*) malloc(num_i *
sizeof(
double));
2759 mfuel_ij_IM = (
double**) malloc(inj_num *
sizeof(
double*));
2760 for(counter = 0; counter < inj_num; counter++) {
2761 num_i = num_i_IM[counter];
2762 mfuel_ij_IM[counter] = (
double*) malloc(num_i *
sizeof(
double));
2766 for(counter = 0; counter < 15; counter++) {
2771 CALCULUS_OF_POI(POI_IM, mfuel_i_IM, mfuel_ij_IM, acu_dmf, time_vector, size, speed, num_i_IM, num_j, SOI_IM, EOI_IM,
2772 inj_num, elementcontrol);
2774 mixture_correction = (
double*) malloc(num_j *
sizeof(
double));
2776 element = (
struct stPropertiesElement * *) malloc(inj_num *
sizeof(
struct stPropertiesElement*));
2777 for(counter = 0; counter < inj_num; counter++) {
2778 element_value = num_i_IM[counter] * num_j;
2779 element[counter] = (
struct stPropertiesElement*) malloc(element_value *
sizeof(
struct stPropertiesElement));
2782 realelement = (
struct stRealElementComposition * *) malloc(inj_num *
sizeof(
struct stRealElementComposition*));
2783 for(counter = 0; counter < inj_num; counter++) {
2784 element_value = num_i_IM[counter] * num_j;
2785 realelement[counter] = (
struct stRealElementComposition*) malloc(element_value *
sizeof(
2786 struct stRealElementComposition));
2789 mixture_correction[0] = 0.319;
2790 mixture_correction[1] = 0.668;
2791 mixture_correction[2] = 1.;
2792 mixture_correction[3] = 1.397;
2793 mixture_correction[4] = 2.07;
2800 mtotal_IVC = mairIVC;
2801 mEGR = mairIVC * rateBG;
2803 YO2IVC = YO2aire * (1 - rateBG * (mfuel * 1.e-3 + mfuel * 1e-3 * AFe) / (mairIVC * (1. - rateBG) + mfuel * 0.001));
2804 YCO2IVC = mmCO2 / mmfuel * mfuel * 0.001 * rateBG / (mfuel * 0.001 + mairIVC * (1 - rateBG));
2805 YH2OIVC = HC / 2. * mmH2O / mmfuel * mfuel * 0.001 * rateBG / (mfuel * 0.001 + mairIVC * (1 - rateBG));
2806 YN2IVC = 1 - YO2IVC - YCO2IVC - YH2OIVC;
2808 mO2_IVC = mairIVC * YO2IVC;
2809 mN2_IVC = mairIVC * YN2IVC;
2810 mCO2_IVC = mairIVC * YCO2IVC;
2811 mH2O_IVC = mairIVC * YH2OIVC;
2813 NOx_IVC = test_variables[14] * 1e-6 * 1.587 * mairIVC;
2814 mSOOT_IVC_A = test_variables[15] * 1e-6 * mairIVC;
2815 mSOOT_IVC_B = test_variables[15] * 1e-6 * mairIVC;
2816 mSOOT_IVC_C = test_variables[15] * 1e-6 * mairIVC;
2820 Cylinder_capacity = PI * Piston_D * Piston_D / 4. * S;
2821 VTDC = Cylinder_capacity / (GCRatio - 1.);
2822 f = mairIVC * 1000. / mfuel;
2826 for(counter = 0; counter < size; counter++) {
2828 mtotal_bowl[counter] = mtotal_IVC;
2829 mO2_bowl[counter] = mO2_IVC;
2830 mN2_bowl[counter] = mN2_IVC;
2831 mCO2_bowl[counter] = mCO2_IVC;
2832 mH2O_bowl[counter] = mH2O_IVC;
2833 mHC_bowl[counter] = 0;
2834 Tsq_cyl[counter] = 0;
2835 Tadib_cyl[counter] = 0;
2836 mNOx_bowl[counter] = NOx_IVC;
2837 mCO_bowl[counter] = 0;
2838 HC_bowl[counter] = 0;
2839 mSOOT_bowl_A[counter] = mSOOT_IVC_A;
2840 mSOOT_bowl_B[counter] = mSOOT_IVC_B;
2841 mSOOT_bowl_C[counter] = mSOOT_IVC_C;
2845 acu_Mbb[counter] = 0;
2849 Yburned[counter] = 0;
2853 H_cooler[counter] = 0;
2856 Rmixture[counter] = 0;
2859 mtotal_Vc[counter] = 0;
2860 mO2_Vc[counter] = 0;
2861 mN2_Vc[counter] = 0;
2862 mCO2_Vc[counter] = 0;
2863 mH2O_Vc[counter] = 0;
2864 mNOx_Vc[counter] = 0;
2865 mSOOT_Vc_A[counter] = 0;
2866 mSOOT_Vc_B[counter] = 0;
2867 mHC_Vc[counter] = 0;
2868 dmtotal_Gfactor[counter] = 0;
2870 Qcylhead[counter] = 0;
2872 evol_mSoot[counter] = 0.;
2873 evol_Radiacion[counter] = 0.;
2878 T_cyl[0] = temperatureIVC;
2879 T_cyl[1] = temperatureIVC;
2885 for(m = 0; m < inj_num; m++) {
2886 for(i = 0; i < num_i_IM[m]; i++) {
2887 for(j = 0; j < num_j; j++) {
2889 aux = i * num_j + j;
2890 element[m][aux].state = stEvaporated;
2891 element[m][aux].mtotal = mfuel_ij_IM[m][i] * (mfuel * 1e-3);
2892 element[m][aux].mf_jet = mfuel_ij_IM[m][i] * (mfuel * 1e-3);
2893 element[m][aux].mf_reac = mfuel_ij_IM[m][i] * (mfuel * 1e-3);
2894 element[m][aux].mf_evap = mfuel_ij_IM[m][i] * (mfuel * 1e-3);
2896 element[m][aux].mO2 = 0;
2897 element[m][aux].mN2 = 0;
2898 element[m][aux].mCO2 = 0.;
2899 element[m][aux].mH2O = 0.;
2901 element[m][aux].dmtotal = 0;
2902 POC_IM[m][aux] = 8888;
2903 element[m][aux].RID = 0;
2904 element[m][aux].Rpmx = 0;
2905 element[m][aux].Rpmx_value = 0;
2906 element[m][aux].mNOx = 0;
2907 element[m][aux].dNOx = 0;
2908 element[m][aux].CO = 0;
2909 element[m][aux].HC = 0;
2910 element[m][aux].mSOOT_A = 0;
2911 element[m][aux].dSOOT_A = 0;
2912 element[m][aux].mSOOT_B = 0;
2913 element[m][aux].dSOOT_B = 0;
2914 element[m][aux].mSOOT_C = 0;
2915 element[m][aux].dSOOT_C = 0;
2916 element[m][aux].TSD = 0;
2917 element[m][aux].Tadib = 0.;
2918 element[m][aux].Pcyl_POC = 1;
2919 element[m][aux].TNOx = 0;
2920 element[m][aux].X = 0;
2921 element[m][aux].soot_precursor = 0;
2922 element[m][aux].C = 1;
2923 element[m][aux].FI = 0;
2924 element[m][aux].FRLOL = 0;
2925 element[m][aux].tLOL =
2928 realelement[m][aux].mtotal = 0;
2929 realelement[m][aux].mO2 = 0;
2930 realelement[m][aux].mO = 0;
2931 realelement[m][aux].mCH = 0;
2932 realelement[m][aux].mN2 = 0;
2933 realelement[m][aux].mCO2 = 0;
2934 realelement[m][aux].mH2O = 0;
2935 realelement[m][aux].mf_reac = 0;
2936 realelement[m][aux].dmf_reac = 0;
2937 realelement[m][aux].mf_reac_pmx = 0;
2942 for(i = 0; i < 15; i++) {
2943 for(j = 0; j < size; j++) {
2944 elementcontrol[i][j].mfuel = 0;
2945 elementcontrol[i][j].mfuel_real = 0;
2946 elementcontrol[i][j].mtotal = 0;
2947 elementcontrol[i][j].mO2 = 0;
2948 elementcontrol[i][j].mO2_real = 0;
2949 elementcontrol[i][j].mO = 0;
2950 elementcontrol[i][j].mCO2 = 0;
2951 elementcontrol[i][j].mH2O = 0;
2952 elementcontrol[i][j].mCH = 0;
2953 elementcontrol[i][j].TSD = 0;
2954 elementcontrol[i][j].Tadib = 0.;
2955 elementcontrol[i][j].dNOx = 0;
2956 elementcontrol[i][j].mNOx = 0;
2957 elementcontrol[i][j].mSOOT_A = 0;
2958 elementcontrol[i][j].mSOOT_B = 0;
2959 elementcontrol[i][j].mSOOT_C = 0;
2960 elementcontrol[i][j].X = 0;
2961 elementcontrol[i][j].FI = 0;
2968 if(RadCalc == 1 && Ang_Grab > -179.) {
2969 finterno = fopen(
"paquete.csv",
"w");
2970 if(finterno == NULL) {
2971 printf(
"Error abriendo fichero interno.");
2975 fprintf(finterno,
"Ang,m_Soot,dm_Soot,mf,dmf \n");
2982 for(counter = 0; counter < size; counter++) {
2988 CALCULATE_CYCLE(roair, CAD, delta_t, V_cyl, VTDC, counter, speed, p_cyl, HRF, acu_dmf, Mbb, acu_Mbb, AFe, f,
2989 mfuel * 1e-6, mEGR * 1e-3, mairIVC * 1e-3, T_cyl, HP, Yair, Yfuel, Yburned, U, CV,
2990 H_cooler, H, T_Evaporation_fuel, inj_fuel_temp, PRECISION_ITERATION, defor, Rmixture, Atmosphere_press, &Gamma, PI,
2991 Runiv, Piston_D, S, Crank_L, Connecting_Rod_L, E, Piston_Axis_D,
2992 Piston_Crown_H, DBowl, VBowl, M_Connecting_Rod, M_P_R_PA, MW_air, MW_fuel, MW_burned, C_ESteel, C_Mech_Defor, CTM, WC1A,
2993 WC1B, C2, C_MBLBY, Cbb, TPIS, TCYL_HEAD, TCYL, Qcylhead, Qcyl,
3000 dmtotal_Gfactor[counter] = (Vc_factor - 1) * (exp(-6.9 * pow(CAD[counter] / 180,
3001 2)) - exp(-6.9 * pow(CAD[counter - 1] / 180, 2))) * (mtotal_bowl[counter - 1]);
3003 mtotal_bowl[counter] = mtotal_bowl[counter - 1] + dmtotal_Gfactor[counter];
3004 mtotal_Vc[counter] = mtotal_Vc[counter - 1] - dmtotal_Gfactor[counter];
3006 if(dmtotal_Gfactor[counter] < 0) {
3008 mO2_bowl[counter] = mO2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mO2_bowl[counter - 1] / mtotal_bowl[counter -
3010 mN2_bowl[counter] = mN2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mN2_bowl[counter - 1] / mtotal_bowl[counter -
3012 mCO2_bowl[counter] = mCO2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mCO2_bowl[counter - 1] / mtotal_bowl[counter -
3014 mH2O_bowl[counter] = mH2O_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mH2O_bowl[counter - 1] / mtotal_bowl[counter -
3016 mNOx_bowl[counter] = mNOx_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mNOx_bowl[counter - 1] / mtotal_bowl[counter -
3018 mSOOT_bowl_A[counter] = mSOOT_bowl_A[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_bowl_A[counter - 1] /
3019 mtotal_bowl[counter - 1]);
3020 mSOOT_bowl_B[counter] = mSOOT_bowl_B[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_bowl_B[counter - 1] /
3021 mtotal_bowl[counter - 1]);
3022 mSOOT_bowl_C[counter] = mSOOT_bowl_C[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_bowl_C[counter - 1] /
3023 mtotal_bowl[counter - 1]);
3025 mHC_bowl[counter] = mHC_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mHC_bowl[counter - 1] / mtotal_bowl[counter -
3028 mO2_Vc[counter] = mO2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mO2_bowl[counter - 1] / mtotal_bowl[counter - 1]);
3029 mN2_Vc[counter] = mN2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mN2_bowl[counter - 1] / mtotal_bowl[counter - 1]);
3030 mCO2_Vc[counter] = mCO2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mCO2_bowl[counter - 1] / mtotal_bowl[counter -
3032 mH2O_Vc[counter] = mH2O_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mH2O_bowl[counter - 1] / mtotal_bowl[counter -
3034 mNOx_Vc[counter] = mNOx_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mNOx_bowl[counter - 1] / mtotal_bowl[counter -
3036 mSOOT_Vc_A[counter] = mSOOT_Vc_A[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_bowl_A[counter - 1] /
3037 mtotal_bowl[counter - 1]);
3038 mSOOT_Vc_B[counter] = mSOOT_Vc_B[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_bowl_B[counter - 1] /
3039 mtotal_bowl[counter - 1]);
3040 mSOOT_Vc_C[counter] = mSOOT_Vc_C[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_bowl_C[counter - 1] /
3041 mtotal_bowl[counter - 1]);
3043 mHC_Vc[counter] = mHC_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mHC_bowl[counter - 1] / mtotal_bowl[counter - 1]);
3047 mO2_bowl[counter] = mO2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mO2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3048 mN2_bowl[counter] = mN2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mN2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3049 mCO2_bowl[counter] = mCO2_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mCO2_Vc[counter - 1] / mtotal_Vc[counter -
3051 mH2O_bowl[counter] = mH2O_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mH2O_Vc[counter - 1] / mtotal_Vc[counter -
3053 mNOx_bowl[counter] = mNOx_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mNOx_Vc[counter - 1] / mtotal_Vc[counter -
3055 mSOOT_bowl_A[counter] = mSOOT_bowl_A[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_Vc_A[counter - 1] /
3056 mtotal_Vc[counter - 1]);
3057 mSOOT_bowl_B[counter] = mSOOT_bowl_B[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_Vc_B[counter - 1] /
3058 mtotal_Vc[counter - 1]);
3059 mSOOT_bowl_C[counter] = mSOOT_bowl_C[counter - 1] + dmtotal_Gfactor[counter] * (mSOOT_Vc_C[counter - 1] /
3060 mtotal_Vc[counter - 1]);
3062 mHC_bowl[counter] = mHC_bowl[counter - 1] + dmtotal_Gfactor[counter] * (mHC_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3064 mO2_Vc[counter] = mO2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mO2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3065 mN2_Vc[counter] = mN2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mN2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3066 mCO2_Vc[counter] = mCO2_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mCO2_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3067 mH2O_Vc[counter] = mH2O_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mH2O_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3068 mNOx_Vc[counter] = mNOx_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mNOx_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3069 mSOOT_Vc_A[counter] = mSOOT_Vc_A[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_Vc_A[counter - 1] / mtotal_Vc[counter
3071 mSOOT_Vc_B[counter] = mSOOT_Vc_B[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_Vc_B[counter - 1] / mtotal_Vc[counter
3073 mSOOT_Vc_C[counter] = mSOOT_Vc_C[counter - 1] - dmtotal_Gfactor[counter] * (mSOOT_Vc_C[counter - 1] / mtotal_Vc[counter
3076 mHC_Vc[counter] = mHC_Vc[counter - 1] - dmtotal_Gfactor[counter] * (mHC_Vc[counter - 1] / mtotal_Vc[counter - 1]);
3083 if(time_vector[counter] < (SOC_IM[0] / (6. * speed))) {
3084 Tsq_cyl[counter] = T_cyl[counter];
3086 Tsq_cyl[counter] = Tsq_cyl[counter - 1] * pow((p_cyl[counter] / p_cyl[counter - 1]), ((Gamma - 1.) / Gamma));
3092 YO2i = mO2_bowl[counter] / mtotal_bowl[counter];
3098 Tadib_cyl[counter] = Tsq_cyl[counter] + 37630.5 * YO2IVC / (YO2IVC + AFe * YO2aire);
3100 if(Tadib_cyl[counter] < 2600) {
3101 Tadib_cyl[counter] = Tadib_cyl[counter] - 1.5538 * 1e-7 * pow(Tadib_cyl[counter], 2.6774);
3103 Tadib_cyl[counter] = Tadib_cyl[counter] - 7.136 * 1e-10 * pow(Tadib_cyl[counter], 3.3596);
3107 Tadib_cyl[counter] *= 0.9;
3110 Zst = YO2i / (YO2i + AFe * YO2aire);
3111 Z0 = -T_cyl[counter] * Zst / (Tadib_cyl[counter] - T_cyl[counter]);
3112 Z1 = 1. + inj_fuel_temp * (1. - Zst) / (Tadib_cyl[counter] - inj_fuel_temp);
3115 double Rhoa = 0., Ta = 0., d0 = 0., EfectoO2 = 0., Pres = 0.,
3117 double x1 = 0., x2 = 0., t1 = 0., t2 = 0., t0 = 0.,
3120 YO2_bowl_i = YO2IVC;
3121 Uo_i = virtual_velocity[counter];
3122 Pres = p_cyl[counter];
3123 Ta = Tsq_cyl[counter];
3124 Rhoa = Pres / (Runiv / MW_air * Ta);
3126 EfectoO2 = (1. + 3.33 / YO2_bowl_i) / (1. + 3.33 / YO2aire);
3127 XLO[counter] = 1.15e11 * Uo_i * pow(Ta, -3.74) * pow(Rhoa, -0.85) * pow(d0 * 1000,
3128 0.34) * EfectoO2 / 1000. * 0.62 / Kmixture1;
3130 PEN_min[counter] = XLO[counter];
3131 PEN_max[counter] = 0.;
3135 if(time_vector[counter] >= SOI_IM[0] * 60. / (360. * speed)) {
3139 for(m = 0; m < inj_num; m++) {
3140 for(i = 0; i < num_i_IM[m]; i++) {
3141 for(j = 0; j < num_j; j++) {
3145 aux = i * num_j + j;
3147 if((POI_IM[m][i] <= time_vector[counter]) && (element[m][aux].mf_jet > 0) && (mO2_bowl[counter] > 0)
3148 && (element[m][aux].state != stOvermixed)) {
3152 element[m][aux].dmtotal = Kmixture1 * (pow(element[m][aux].mf_jet,
3153 2) / element[m][aux].mtotal) * virtual_velocity[counter] * pow(roair[counter], 0.5) * pow(
3154 nozzle_d * 1000, -1) * pow((mO2_bowl[counter] / mtotal_bowl[counter]),
3155 Kmixture2) * (1 / mixture_correction[j]) * delta_t;
3157 element[m][aux].mtotal = element[m][aux].mtotal + element[m][aux].dmtotal;
3158 element[m][aux].mO2 = element[m][aux].mO2 + element[m][aux].dmtotal * (mO2_bowl[counter] / mtotal_bowl[counter]);
3159 element[m][aux].mCO2 = element[m][aux].mCO2 + element[m][aux].dmtotal * (mCO2_bowl[counter] / mtotal_bowl[counter]);
3160 element[m][aux].mH2O = element[m][aux].mH2O + element[m][aux].dmtotal * (mH2O_bowl[counter] / mtotal_bowl[counter]);
3161 element[m][aux].mN2 = element[m][aux].mN2 + element[m][aux].dmtotal * (mN2_bowl[counter] / mtotal_bowl[counter]);
3162 element[m][aux].mf_reac = element[m][aux].mf_reac + element[m][aux].dmtotal * (mHC_bowl[counter] /
3163 mtotal_bowl[counter]);
3164 element[m][aux].mf_evap = element[m][aux].mf_evap + element[m][aux].dmtotal * (mHC_bowl[counter] /
3165 mtotal_bowl[counter]);
3169 element[m][aux].FI = element[m][aux].mf_reac * Kst1 / element[m][aux].mO2;
3170 element[m][aux].C = element[m][aux].mf_reac / element[m][aux].mtotal;
3174 element[m][aux].X += element[m][aux].C * virtual_velocity[counter] * delta_t;
3177 if(element[m][aux].tLOL < -10) {
3179 if(element[m][aux].X >= XLO[counter]) {
3181 x2 = element[m][aux].X;
3182 x1 = x2 - element[m][aux].C * virtual_velocity[counter] * delta_t;
3184 t1 = time_vector[counter - 1];
3185 t2 = time_vector[counter];
3186 t0 = (x2 * x2 * t1 - x1 * x1 * t2) / (x2 * x2 - x1 * x1);
3187 K = x2 / sqrt(t2 - t0);
3189 element[m][aux].tLOL = t0 + pow(XLO[counter] / K, 2);
3192 K = element[m][aux].C * x2 / XLO[counter];
3195 K = 1. / (1. / K - 1.);
3198 element[m][aux].FRLOL = K * Kst1 / (mO2_bowl[counter] / mtotal_bowl[counter]);
3199 FRLOL = element[m][aux].FRLOL;
3201 element[m][aux].tLOL = time_vector[counter];
3202 element[m][aux].FRLOL = element[m][aux].FI;
3203 FRLOL = element[m][aux].FRLOL;
3210 element[m][aux].mNOx = element[m][aux].mNOx + element[m][aux].dmtotal * (mNOx_bowl[counter] / mtotal_bowl[counter]);
3214 element[m][aux].dSOOT_A = element[m][aux].dmtotal * (mSOOT_bowl_A[counter] / mtotal_bowl[counter]);
3215 element[m][aux].dSOOT_B = element[m][aux].dmtotal * (mSOOT_bowl_B[counter] / mtotal_bowl[counter]);
3216 element[m][aux].dSOOT_C = element[m][aux].dmtotal * (mSOOT_bowl_C[counter] / mtotal_bowl[counter]);
3219 if((time_vector[counter] < (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3220 element[m][aux].mSOOT_A = element[m][aux].mSOOT_A + element[m][aux].dSOOT_A;
3221 mSOOT_bowl_A[counter] = mSOOT_bowl_A[counter] - element[m][aux].dSOOT_A;
3223 if((time_vector[counter] > (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3224 mSOOT_bowl_A[counter] = mSOOT_bowl_A[counter] - element[m][aux].dSOOT_A;
3229 if((time_vector[counter] < (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3230 element[m][aux].mSOOT_B = element[m][aux].mSOOT_B + element[m][aux].dSOOT_B;
3231 mSOOT_bowl_B[counter] = mSOOT_bowl_B[counter] - element[m][aux].dSOOT_B;
3233 if((time_vector[counter] > (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3234 mSOOT_bowl_B[counter] = mSOOT_bowl_B[counter] - element[m][aux].dSOOT_B;
3239 if((time_vector[counter] < (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3240 element[m][aux].mSOOT_C = element[m][aux].mSOOT_C + element[m][aux].dSOOT_C;
3241 mSOOT_bowl_C[counter] = mSOOT_bowl_C[counter] - element[m][aux].dSOOT_C;
3243 if((time_vector[counter] > (SOC_IM[m] / (6 * speed))) && (element[m][aux].FI > 1)) {
3244 mSOOT_bowl_C[counter] = mSOOT_bowl_C[counter] - element[m][aux].dSOOT_C;
3249 if((element[m][aux].state == stEvaporated) && (SOC_IM[m] == EVO)) {
3251 Y1 = -KID3 * pow(element[m][aux].mf_reac / element[m][aux].mtotal,
3252 KID4) * pow(element[m][aux].mO2 / element[m][aux].mtotal, KID5);
3254 if(element[m][aux].mf_reac / element[m][aux].mtotal == 1) {
3257 T1 = 1000. / (T_cyl[counter] * (1 - element[m][aux].mf_reac / element[m][aux].mtotal));
3261 RID1 = Y1 + 0.155 * T1 - 0.2175;
3262 }
else if(((T1 - 0.35) > 0.8) && (T1 < 2)) {
3263 RID1 = (Y1 + pow(exp((T1 - 0.35) / 5), 7) - 3.1);
3264 element[m][aux].RID = element[m][aux].RID + KID1 * pow(p_cyl[counter] * 1e-5, KID2) * pow(10, -RID1) * delta_t;
3266 RID1 = (Y1 + pow(log(T1 - 0.35), 3));
3267 element[m][aux].RID = element[m][aux].RID + KID1 * pow(p_cyl[counter] * 1e-5, KID2) * pow(10, -RID1) * delta_t;
3271 if((element[m][aux].RID >= 1) && (SOC_IM[m] == EVO)) {
3272 SOC_IM[m] = time_vector[counter] * 6 * speed;
3279 if((element[m][aux].state == stEvaporated) && (time_vector[counter] > (SOC_IM[m] / (6 * speed)))) {
3281 element[m][aux].Rpmx = element[m][aux].Rpmx + Kpmx1 * pow((mO2_bowl[counter] / mtotal_bowl[counter]),
3282 0.855) * ((1. / 4. - 1.) / pow(log(.1), 2) * pow(
3283 log(element[m][aux].FI), 2) + 1.) * exp(-Kpmx2 / T_cyl[counter]) * pow(roair[counter], Kpmx5) * pow(speed / 1000.,
3284 Kpmx3) * pow(PCR / 500., Kpmx4) * delta_t;
3288 if(element[m][aux].Rpmx > 1.)
3289 element[m][aux].Rpmx = 1.;
3292 if(element[m][aux].Pcyl_POC < 2.)
3294 element[m][aux].Pcyl_POC = p_cyl[counter];
3300 if((element[m][aux].state == stEvaporated) && (POI_IM[m][i] < SOC_IM[m] / (6 * speed))
3301 && (time_vector[counter] > (SOC_IM[m] / (6 * speed)))) {
3302 if(POI_IM[m][i] < SOC_IM[m] / (6 * speed)) {
3304 POC_IM[m][aux] = time_vector[counter];
3306 if(element[m][aux].FI > 1.) {
3308 realelement[m][aux].mf_reac_pmx = element[m][aux].mf_reac * (1. - 1. / element[m][aux].FI * element[m][aux].Rpmx);
3309 if(element[m][aux].Rpmx == 1.)
3311 element[m][aux].state = stBurned_rich_premix;
3314 element[m][aux].Pcyl_POC = p_cyl[counter];
3315 realelement[m][aux].mf_reac_pmx = element[m][aux].mf_reac * (1. - element[m][aux].Rpmx);
3316 if(element[m][aux].Rpmx == 1.)
3318 element[m][aux].state = stBurned_poor_premix;
3325 if((element[m][aux].state == stEvaporated) && (POI_IM[m][i] >= SOC_IM[m] / (6 * speed))) {
3326 element[m][aux].state = stInto_diffusion_flame;
3332 if((element[m][aux].state == stInto_diffusion_flame) && (POI_IM[m][i] >= SOC_IM[m] / (6 * speed))
3333 && (time_vector[counter] > (SOC_IM[m] / (6 * speed)))) {
3335 if(element[m][aux].FI <= 1) {
3336 POC_IM[m][aux] = time_vector[counter];
3337 element[m][aux].Pcyl_POC = p_cyl[counter];
3338 element[m][aux].state = stBurned_by_diffusion;
3345 if((element[m][aux].state == stBurned_rich_premix)) {
3347 if(element[m][aux].FI <= 1) {
3348 element[m][aux].state = stBurned_by_second_diffusion;
3349 POC_IM[m][aux] = time_vector[counter];
3350 element[m][aux].Pcyl_POC = p_cyl[counter];
3357 if((element[m][aux].state == stEvaporated) && (time_vector[counter] > (SOC_IM[m] / (6 * speed)))) {
3358 if(element[m][aux].FI > 1.) {
3360 realelement[m][aux].mO2 = element[m][aux].mO2 * (1. - element[m][aux].Rpmx);
3361 mf_old = realelement[m][aux].mf_reac;
3362 realelement[m][aux].mf_reac = element[m][aux].mf_reac * (1. - 1. / element[m][aux].FI * element[m][aux].Rpmx);
3363 realelement[m][aux].dmf_reac = realelement[m][aux].mf_reac - mf_old;
3366 realelement[m][aux].mO2 = element[m][aux].mO2 - element[m][aux].mf_reac * Kst1 * element[m][aux].Rpmx;
3367 mf_old = realelement[m][aux].mf_reac;
3368 realelement[m][aux].mf_reac = element[m][aux].mf_reac * (1. - element[m][aux].Rpmx);
3369 realelement[m][aux].dmf_reac = realelement[m][aux].mf_reac - mf_old;
3373 if((element[m][aux].state == stBurned_rich_premix) || (element[m][aux].state == stInto_diffusion_flame)) {
3374 realelement[m][aux].mO2 = 0;
3375 mf_old = realelement[m][aux].mf_reac;
3376 realelement[m][aux].mf_reac = element[m][aux].mf_reac * (1. - 1. / element[m][aux].FI);
3377 realelement[m][aux].dmf_reac = realelement[m][aux].mf_reac - mf_old;
3380 if((element[m][aux].state == stBurned_poor_premix) || (element[m][aux].state == stBurned_by_diffusion)
3381 || (element[m][aux].state == stBurned_by_second_diffusion)) {
3382 realelement[m][aux].mO2 = element[m][aux].mO2 - element[m][aux].mf_reac * Kst1;
3383 mf_old = realelement[m][aux].mf_reac;
3384 realelement[m][aux].mf_reac = 0;
3385 realelement[m][aux].dmf_reac = -mf_old;
3395 if(element[m][aux].FI >= 1.) {
3397 TCC = Tadib_cyl[counter] * (Z1 - element[m][aux].C) / (Z1 - Zst);
3400 TCC = Tadib_cyl[counter] * (element[m][aux].C - Z0) / (Zst - Z0);
3405 TSC = (T_cyl[counter] * (1. - element[m][aux].C) + inj_fuel_temp * element[m][aux].C);
3409 if(time_vector[counter] > (SOC_IM[m] / (6 * speed)))
3411 element[m][aux].Tadib = TCC;
3414 element[m][aux].Tadib = TSC;
3417 element[m][aux].TNOx = element[m][aux].Tadib;
3423 if((element[m][aux].state == stBurned_by_diffusion) || (element[m][aux].state == stBurned_by_second_diffusion)
3424 || (element[m][aux].state == stBurned_poor_premix)) {
3426 if((time_vector[counter] == POC_IM[m][aux]) && (element[m][aux].state != stBurned_poor_premix)) {
3427 mNOx_bowl[counter] = mNOx_bowl[counter] - element[m][aux].mNOx * EC;
3428 element[m][aux].mNOx = element[m][aux].mNOx - element[m][aux].mNOx * EC;
3431 FUNCTION_NOX(&YNOeq_value, &KdYNO_value, YNOeq, KdYNO, element[m][aux].TNOx, realelement[m][aux].mO2,
3432 element[m][aux].mtotal);
3434 if(YNOeq_value == 0) {
3438 if(YNOeq_value - element[m][aux].mNOx / element[m][aux].mtotal > 0) {
3439 element[m][aux].dNOx = (KNOx2) * element[m][aux].mtotal * KdYNO_value * (YNOeq_value - element[m][aux].mNOx /
3440 element[m][aux].mtotal) / YNOeq_value * delta_t;
3442 element[m][aux].dNOx = 0;
3445 element[m][aux].mNOx = element[m][aux].mNOx + element[m][aux].dNOx;
3446 mNOx_bowl[counter] = mNOx_bowl[counter] + element[m][aux].dNOx;
3454 element[m][aux].dSOOT_A = KSOOTA1 * realelement[m][aux].mf_reac * pow(p_cyl[counter],
3455 KSOOTA2) * exp(-KSOOTA3 / element[m][aux].Tadib) * delta_t;
3456 element[m][aux].mSOOT_A = element[m][aux].mSOOT_A + element[m][aux].dSOOT_A;
3460 element[m][aux].dSOOT_A = KSOOTA4 * element[m][aux].mSOOT_A * pow((realelement[m][aux].mO2 / element[m][aux].mtotal),
3461 KSOOTA5) * pow(p_cyl[counter], KSOOTA6) * exp(
3462 -KSOOTA7 / element[m][aux].Tadib) * delta_t;
3463 if(element[m][aux].dSOOT_A <= element[m][aux].mSOOT_A) {
3464 element[m][aux].mSOOT_A = element[m][aux].mSOOT_A - element[m][aux].dSOOT_A;
3466 element[m][aux].mSOOT_A = 0;
3474 if(realelement[m][aux].dmf_reac < 0.)
3475 element[m][aux].mSOOT_C *= realelement[m][aux].mf_reac / (realelement[m][aux].mf_reac - realelement[m][aux].dmf_reac);
3478 if(element[m][aux].tLOL > -10.) {
3481 if((time_vector[counter] - element[m][aux].tLOL) > 0.) {
3483 valor = KSOOTC1 * realelement[m][aux].mf_reac * pow(Pres / (14.8 * Runiv / MW_air * 1000.),
3484 2.2) * (-1. / element[m][aux].FRLOL + 0.5) * 1.e6 * exp(
3485 -2399. / Tadib_cyl[counter]) * delta_t;
3497 element[m][aux].mSOOT_C += valor;
3500 if(element[m][aux].mSOOT_C > realelement[m][aux].mf_reac)
3501 element[m][aux].mSOOT_C = realelement[m][aux].mf_reac;
3503 if(element[m][aux].mSOOT_C < 0.)
3505 element[m][aux].mSOOT_C = 0.;
3508 evol_mSoot[counter] += element[m][aux].mSOOT_C;
3513 if((SOC_IM[m] * 60. / (360. * speed) < time_vector[counter]) && (POI_IM[m][i] < SOC_IM[m] * 60. / (360. * speed))
3514 && (element[m][aux].state == stEvaporated)) {
3515 if(element[m][aux].FI < 0.5) {
3517 mHC_bowl[counter] = mHC_bowl[counter] + element[m][aux].mf_reac;
3518 element[m][aux].state = stOvermixed;
3519 realelement[m][aux].mf_reac = 0;
3520 element[m][aux].mf_reac = 0;
3521 element[m][aux].mf_jet = 0;
3529 for(i_aux = 0; i_aux < 5; i_aux++) {
3532 if((elementcontrol[j_aux][0].inj_number == m) && (elementcontrol[j_aux][0].num_i == i)) {
3533 aux = i * num_j + 0;
3535 elementcontrol[j_aux][counter].mtotal = element[m][aux].mtotal;
3536 elementcontrol[j_aux][counter].mfuel = element[m][aux].mf_reac;
3537 elementcontrol[j_aux][counter].mfuel_real = realelement[m][aux].mf_reac;
3538 elementcontrol[j_aux][counter].mO2 = element[m][aux].mO2;
3539 elementcontrol[j_aux][counter].mO2_real = realelement[m][aux].mO2;
3540 elementcontrol[j_aux][counter].TSD = element[m][aux].TSD;
3541 elementcontrol[j_aux][counter].Tadib = element[m][aux].Tadib;
3542 elementcontrol[j_aux][counter].X = element[m][aux].X;
3543 elementcontrol[j_aux][counter].FI = element[m][aux].FI;
3545 aux = i * num_j + 2;
3546 j_aux = i_aux * 3 + 1;
3547 elementcontrol[j_aux][counter].mtotal = element[m][aux].mtotal;
3548 elementcontrol[j_aux][counter].mfuel = element[m][aux].mf_reac;
3549 elementcontrol[j_aux][counter].mfuel_real = realelement[m][aux].mf_reac;
3550 elementcontrol[j_aux][counter].mO2 = element[m][aux].mO2;
3551 elementcontrol[j_aux][counter].mO2_real = realelement[m][aux].mO2;
3552 elementcontrol[j_aux][counter].TSD = element[m][aux].TSD;
3553 elementcontrol[j_aux][counter].Tadib = element[m][aux].Tadib;
3554 elementcontrol[j_aux][counter].X = element[m][aux].X;
3555 elementcontrol[j_aux][counter].FI = element[m][aux].FI;
3557 aux = i * num_j + 4;
3558 j_aux = i_aux * 3 + 2;
3559 elementcontrol[j_aux][counter].mtotal = element[m][aux].mtotal;
3560 elementcontrol[j_aux][counter].mfuel = element[m][aux].mf_reac;
3561 elementcontrol[j_aux][counter].mfuel_real = realelement[m][aux].mf_reac;
3562 elementcontrol[j_aux][counter].mO2 = element[m][aux].mO2;
3563 elementcontrol[j_aux][counter].mO2_real = realelement[m][aux].mO2;
3564 elementcontrol[j_aux][counter].TSD = element[m][aux].TSD;
3565 elementcontrol[j_aux][counter].Tadib = element[m][aux].Tadib;
3566 elementcontrol[j_aux][counter].X = element[m][aux].X;
3567 elementcontrol[j_aux][counter].FI = element[m][aux].FI;
3575 for(m = 0; m < inj_num; m++) {
3576 for(i = 0; i < num_i_IM[m]; i++) {
3577 for(j = 0; j < num_j; j++) {
3578 if(POI_IM[m][i] <= time_vector[counter]) {
3579 aux = i * num_j + j;
3581 if(((element[m][aux].state == stEvaporated) && (time_vector[counter] > (SOC_IM[m] / (6 * speed))))
3582 || (element[m][aux].state == stBurned_rich_premix) || (element[m][aux].state == stBurned_poor_premix)
3583 || (element[m][aux].state == stBurned_by_diffusion) || (element[m][aux].state == stBurned_by_second_diffusion)
3584 || (element[m][aux].state == stInto_diffusion_flame)) {
3585 mf_burned = mf_burned + (element[m][aux].mf_reac - realelement[m][aux].mf_reac);
3588 if(((element[m][aux].state == stEvaporated) && (time_vector[counter] > (SOC_IM[m] / (6 * speed))))
3589 || (element[m][aux].state == stBurned_rich_premix) || (element[m][aux].state == stBurned_poor_premix)
3590 || (element[m][aux].state == stBurned_by_second_diffusion)) {
3591 mf_burned_pmx = mf_burned_pmx + (element[m][aux].mf_reac - realelement[m][aux].mf_reac_pmx);
3600 HRF[counter] = mf_burned / (mfuel * 1.e-3);
3601 HRF_PMX[counter] = mf_burned_pmx / (mfuel * 1.e-3);
3604 mO2_bowl[counter] = mO2_bowl[counter] - (HRF[counter] - HRF[counter - 1]) * mfuel * 1.e-3 * Kst1;
3605 mCO2_bowl[counter] = mCO2_bowl[counter] + (HRF[counter] - HRF[counter - 1]) * mfuel * 1.e-3 * Kst5;
3606 mH2O_bowl[counter] = mH2O_bowl[counter] + (HRF[counter] - HRF[counter - 1]) * mfuel * 1.e-3 * Kst6;
3609 mSootCil[counter] = mSOOT_IVC_C * (mO2_bowl[counter] + mO2_Vc[counter]) / mO2_IVC;
3615 double lambda = 0., dlambda = 0., tgAng = 0., Teje = 0., Tllama = 0., Xeje = 0., radio = 0., aux1 = 0., Xmax = 0.;
3616 int k = 0, iini = 0, ifin = 0;
3618 if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3619 && (RadCalc == 1)) {
3622 fprintf(foculto,
"Time,%le,[s],Ang,%le,[cig],LOL,%le,[mm],FrLOL,%le\n#iny,Elem,r,x,msoot,mtot,rho,T,FI,mf_reac\n",
3623 time_vector[counter], Ang_Grab, XLO[counter] * 1000., FRLOL);
3626 for(m = 0; m < inj_num; m++) {
3627 for(i = 0; i < num_i_IM[m]; i++) {
3628 for(j = 0; j < num_j; j++) {
3632 aux = i * num_j + j;
3633 fprintf(foculto,
"%d,%d,%d,%le,%le,%le,%le,%le,%le,%le\n", m + 1, i + 1, j + 1, element[m][aux].X,
3634 element[m][aux].mSOOT_C, element[m][aux].mtotal, roair[counter],
3635 element[m][aux].Tadib, element[m][aux].FI, realelement[m][aux].mf_reac);
3644 Itot = (
double*) malloc(inj_num *
sizeof(
double));
3649 for(m = 0; m < inj_num; m++) {
3654 if(CalcRad[m] == 1 && RadCalc == 1) {
3657 aux = i * num_j + j;
3659 while((element[m][aux].mSOOT_C < 1.e-20) && (i < (num_i_IM[m] - 1))) {
3661 aux = i * num_j + j;
3664 while((element[m][aux].mSOOT_C > 0.) && (i < (num_i_IM[m] - 1))) {
3666 aux = i * num_j + j;
3668 if(i == (num_i_IM[m] - 1)) {
3669 if(element[m][aux].mSOOT_C > 0.)
3682 tgAng = tan(ANG_CHORRO / 2 * PI / 180.) * Kmixture1 / 0.75;
3687 Radiation[m][0].x = (element[m][iini * num_j + 4].X * (element[m][(iini - 1) * num_j + 4].FI - 1.) + element[m][(iini -
3688 1) * num_j + 4].X * (1. - element[m][iini * num_j + 4].FI)) / (element[m][(iini - 1) * num_j + 4].FI - element[m][iini *
3690 if(PEN_max[counter] < Radiation[m][0].x)
3691 PEN_max[counter] = Radiation[m][0].x;
3694 Radiation[m][0].x = element[m][4].X * Kst1 * element[m][4].mf_jet / element[m][4].mO2;
3695 if(PEN_max[counter] < element[m][4].X)
3696 PEN_max[counter] = element[m][4].X;
3701 Radiation[m][0].RFlame = 0.;
3702 Radiation[m][0].R = Radiation[m][0].x * tgAng;
3703 Radiation[m][0].np = 1;
3704 Radiation[m][0].dr = 1.e-3;
3705 Radiation[m][0].T[0] = Tadib_cyl[counter];
3706 Radiation[m][0].Xsoot[0] = 0.;
3707 Radiation[m][0].Tau[0] = 1.;
3708 Radiation[m][0].PTau[0] = 1.;
3709 Radiation[m][0].Ilambda = 0.;
3710 Radiation[m][0].Itot = 0.;
3712 for(i = (iini + 1); i <= ifin; i++) {
3714 aux = i * num_j + j;
3715 Radiation[m][i - iini].x = element[m][aux].X;
3716 aux1 = Radiation[m][i - iini].x / Radiation[m][0].x;
3719 aux1 = -log(aux1) / 4.6;
3722 Radiation[m][i - iini].RFlame = Radiation[m][i - iini].x * tgAng * sqrt(aux1);
3723 Radiation[m][i - iini].R = Radiation[m][i - iini].x * tgAng;
3725 if(Radiation[m][i - iini].RFlame > cero) {
3726 Radiation[m][i - iini].np = NR;
3728 Radiation[m][i - iini].dr = Radiation[m][i - iini].RFlame * 2. / (float) Radiation[m][i - iini].np;
3730 Teje = element[m][aux].Tadib;
3731 Tllama = Tadib_cyl[counter];
3732 Xeje = element[m][aux].mSOOT_C / element[m][aux].mtotal * p_cyl[counter] / (287. * Teje) /
3736 for(k = 0; k < Radiation[m][i - iini].np; k++) {
3737 radio = -Radiation[m][i - iini].RFlame + (2. * k + 1.) / 2. * Radiation[m][i - iini].dr;
3738 Radiation[m][i - iini].T[k] = Teje + (Radiation[m][0].x / Radiation[m][i - iini].x * (exp(-4.6 * pow(
3739 radio / Radiation[m][i - iini].R, 2)) - 1.) / (1. - Radiation[m][0].x / Radiation[m][i - iini].x)) * (Tllama - Teje);
3741 Radiation[m][i - iini].Xsoot[k] = Xeje * exp(-4.6 * pow(radio / Radiation[m][i - iini].RFlame, 2));
3744 Radiation[m][i - iini].np = 1.;
3745 Radiation[m][i - iini].dr = 0.001;
3746 Teje = element[m][aux].Tadib;
3748 Xeje = element[m][aux].mSOOT_C / element[m][aux].mtotal * p_cyl[counter] / (287. * Teje) / 1800.;
3752 Radiation[m][i - iini].T[0] = Teje;
3753 Radiation[m][i - iini].Xsoot[0] = Xeje;
3755 Radiation[m][i - iini].Itot = 0.;
3761 if(ifin < (num_i_IM[m] - 1)) {
3763 if(XLO[counter] < Radiation[m][ifin - iini].x) {
3764 Radiation[m][ifin + 1 - iini].x = XLO[counter];
3765 aux = (ifin + 1) * num_j + 4;
3767 Radiation[m][ifin + 1 - iini].T[0] = (element[m][aux].Tadib * (element[m][aux - num_j].X - Radiation[m][ifin + 1 -
3768 iini].x) + element[m][aux - num_j].Tadib * (Radiation[m][ifin + 1 - iini].x - element[m][aux].X)) /
3769 (element[m][aux - num_j].X - element[m][aux].X);
3772 Radiation[m][ifin + 1 - iini].x = Radiation[m][ifin - iini].x;
3773 Radiation[m][ifin + 1 - iini].T[0] = Radiation[m][ifin - iini].T[(Radiation[m][ifin - iini].np - 1) / 2];
3777 Radiation[m][ifin + 1 - iini].x = Radiation[m][ifin - iini].x;
3778 PEN_min[counter] = Radiation[m][ifin - iini].x;
3779 Radiation[m][ifin + 1 - iini].T[0] = Radiation[m][ifin - iini].T[(Radiation[m][ifin - iini].np - 1) / 2];
3781 aux1 = Radiation[m][ifin + 1 - iini].x / Radiation[m][0].x;
3784 aux1 = -log(aux1) / 4.6;
3787 Radiation[m][ifin + 1 - iini].RFlame = Radiation[m][ifin + 1 - iini].x * tgAng * sqrt(aux1);
3788 Radiation[m][ifin + 1 - iini].R = Radiation[m][ifin + 1 - iini].x * tgAng;
3789 Radiation[m][ifin + 1 - iini].np = 1;
3790 Radiation[m][ifin + 1 - iini].dr = 0.001;
3791 Radiation[m][ifin + 1 - iini].Xsoot[0] = 0.;
3792 Radiation[m][ifin + 1 - iini].Tau[0] = 1.;
3793 Radiation[m][ifin + 1 - iini].PTau[0] = 1.;
3794 Radiation[m][ifin + 1 - iini].Ilambda = 0.;
3795 Radiation[m][ifin + 1 - iini].Itot = 0.;
3797 if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3798 && (RadCalc == 1)) {
3800 fprintf(foculto2,
"Time,%le,[s],Ang,%le,[cig],LOL,%le,[mm]\n", time_vector[counter], Ang_Grab, XLO[counter] * 1000.);
3802 fprintf(foculto2,
"x");
3803 for(i = iini; i <= (ifin + 1); i++)
3804 fprintf(foculto2,
",%le", Radiation[m][i - iini].x);
3805 fprintf(foculto2,
"\n");
3807 fprintf(foculto2,
"Rllama");
3808 for(i = iini; i <= (ifin + 1); i++)
3809 fprintf(foculto2,
",%le", Radiation[m][i - iini].RFlame);
3810 fprintf(foculto2,
"\n");
3812 fprintf(foculto2,
"R");
3813 for(i = iini; i <= (ifin + 1); i++)
3814 fprintf(foculto2,
",%le", Radiation[m][i - iini].R);
3815 fprintf(foculto2,
"\n");
3817 fprintf(foculto2,
"Xsoot");
3818 fprintf(foculto2,
",0.0");
3819 for(i = iini + 1; i <= ifin; i++)
3820 fprintf(foculto2,
",%le", Radiation[m][i - iini].Xsoot[(NR - 1) / 2]);
3821 fprintf(foculto2,
",0.0\n");
3823 fprintf(foculto2,
"T");
3824 for(i = iini; i <= (ifin + 1); i++) {
3825 if(Radiation[m][i - iini].np == NR)
3826 aux1 = Radiation[m][i - iini].T[(NR - 1) / 2];
3828 aux1 = Radiation[m][i - iini].T[0];
3829 fprintf(foculto2,
",%le", aux1);
3831 fprintf(foculto2,
"\n");
3836 for(lambda = dlambda; lambda < 5.; lambda = lambda + dlambda) {
3837 for(i = (iini + 1); i <= ifin; i++) {
3839 for(k = 0; k < Radiation[m][i - iini].np; k++) {
3840 Radiation[m][i - iini].Tau[k] = exp(-6.3 * Radiation[m][i - iini].Xsoot[k] * 1.e6 * Radiation[m][i - iini].dr / pow(
3841 lambda, 1.22 - 0.245 * log(lambda)));
3843 Radiation[m][i - iini].PTau[k] = Radiation[m][i - iini].Tau[k];
3845 Radiation[m][i - iini].PTau[k] = Radiation[m][i - iini].Tau[k] * Radiation[m][i - iini].PTau[k - 1];
3850 Radiation[m][i - iini].Ilambda = 0.;
3851 for(k = 0; k < Radiation[m][i - iini].np; k++) {
3852 if(Radiation[m][i - iini].T[k] > cero) {
3854 aux1 = 3.743e8 / (PI * pow(lambda, 5) * (exp(1.4387e4 / (lambda * Radiation[m][i - iini].T[k])) - 1.));
3855 aux1 = aux1 * (1. - Radiation[m][i - iini].Tau[k]) * Radiation[m][i - iini].PTau[k];
3856 Radiation[m][i - iini].Ilambda += aux1 * 1.e6;
3860 Radiation[m][i - iini].Itot += Radiation[m][i - iini].Ilambda;
3863 if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3864 && (RadCalc == 1)) {
3867 fprintf(foculto2,
"%le", lambda);
3868 for(i = iini; i <= (ifin + 1); i++)
3869 fprintf(foculto2,
",%le", Radiation[m][i - iini].Ilambda);
3870 fprintf(foculto2,
"\n");
3873 fprintf(foculto2,
"%le", lambda);
3874 fprintf(foculto2,
",0.0");
3875 for(i = iini + 1; i <= ifin; i++)
3876 fprintf(foculto2,
",%le", -log(Radiation[m][i - iini].PTau[NR - 1]));
3877 fprintf(foculto2,
",0.0\n");
3883 for(i = (iini + 1); i <= ifin; i++) {
3885 Radiation[m][i - iini].Itot -= Radiation[m][i - iini].Ilambda / 2.;
3886 Radiation[m][i - iini].Itot *= dlambda * 1.e-6;
3889 Radiation[m][i - iini].Itot += Radiation[m][i - iini].Ilambda / 2. * (10. - lambda) * 1.e-6;
3893 if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3894 && (RadCalc == 1)) {
3896 fprintf(foculto2,
"Itot");
3898 for(i = iini; i <= (ifin + 1); i++)
3899 fprintf(foculto2,
",%le", Radiation[m][i - iini].Itot);
3900 fprintf(foculto2,
"\n");
3905 for(i = (iini + 2); i <= (ifin + 1); i++) {
3906 Itot[m] += (Radiation[m][i - iini - 1].Itot + Radiation[m][i - iini].Itot) / 2. * (Radiation[m][i - iini - 1].x -
3907 Radiation[m][i - iini].x) * (Radiation[m][i - iini - 1].RFlame + Radiation[m][i - iini].RFlame) * PI;
3910 for(i = (iini + 1); i <= (ifin + 1); i++) {
3911 Itot[m] += (Radiation[m][i - iini - 1].Itot + Radiation[m][i - iini].Itot) / 2. * (Radiation[m][i - iini - 1].x -
3912 Radiation[m][i - iini].x) * (Radiation[m][i - iini - 1].RFlame + Radiation[m][i - iini].RFlame) * PI;
3918 aux1 = 1. / (1. - pow(pow(Radiation[m][1].x / Radiation[m][0].x,
3919 2) * mO2_bowl[counter] / mtotal_bowl[counter] * Xmax / (0.23 * 0.004298537), 1. / 2.5));
3927 if((time_vector[counter] <= (Ang_Grab / (6 * speed))) && (time_vector[counter + 1] > (Ang_Grab / (6 * speed)))
3929 fprintf(foculto2,
"ITOT,%le\n", Itot[m]);
3937 evol_Radiacion[counter] = 0.;
3938 for(m = 0; m < inj_num; m++)
3939 evol_Radiacion[counter] += Itot[m];
3958 if(RadCalc == 1 && Ang_Grab > -179.)
3959 fprintf(finterno,
"%le,%le,%le,%le,%le\n", time_vector[counter] * 6. * speed, element[pulso_i][elem_i].mSOOT_C,
3960 element[pulso_i][elem_i].dSOOT_C, realelement[pulso_i][elem_i].mf_reac,
3961 realelement[pulso_i][elem_i].dmf_reac);
3964 if(PEN_min[counter] > PEN_max[counter])
3965 PEN_min[counter] = 0.;
3969 if(RadCalc == 1 && Ang_Grab > -179.) {
3977 if(RadCalc == 1 && Ang_Grab > -179.)
3980 mSOOT_bowl_A_i_burned = mSOOT_bowl_A[size - 1];
3981 mSOOT_bowl_B_i_burned = mSOOT_bowl_B[size - 1];
3982 mSOOT_bowl_C_i_burned = mSOOT_bowl_C[size - 1];
3984 for(m = 0; m < inj_num; m++) {
3985 for(i = 0; i < num_i_IM[m]; i++) {
3986 for(j = 0; j < num_j; j++) {
3987 aux = i * num_j + j;
3988 mSOOT_bowl_A[size - 1] = mSOOT_bowl_A[size - 1] + element[m][aux].mSOOT_A;
3989 mSOOT_bowl_B[size - 1] = mSOOT_bowl_B[size - 1] + element[m][aux].mSOOT_B;
3990 mSOOT_bowl_C[size - 1] = mSOOT_bowl_C[size - 1] + element[m][aux].mSOOT_C;
3992 if(element[m][aux].FI <= 1) {
3993 mSOOT_bowl_A_i_burned = mSOOT_bowl_A_i_burned + element[m][aux].mSOOT_A;
3994 mSOOT_bowl_B_i_burned = mSOOT_bowl_B_i_burned + element[m][aux].mSOOT_B;
3995 mSOOT_bowl_C_i_burned = mSOOT_bowl_C_i_burned + element[m][aux].mSOOT_C;
4001 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4002 vector_to_interpolate[0] = 0;
4003 for(counter = 1; counter < size; counter++) {
4004 vector_to_interpolate[counter] = (HRF[counter] - HRF[counter - 1]) * (mfuel * 1e-6) * HP /
4005 (time_vector[counter] - time_vector[counter - 1]);
4013 for(counter = 1; counter < size - 1; counter++) {
4014 ROHR[counter] = vector_to_interpolate[counter];
4016 free(vector_to_interpolate);
4019 complete_prev_size = 0;
4020 complete_post_size = 0;
4021 counter_CAD_1 = IVC;
4022 counter_CAD_2 = EVO;
4025 complete_prev_size = complete_prev_size + 1;
4026 counter_CAD_1 = counter_CAD_1 - delta_t * 360. * speed / 60.;
4027 }
while(counter_CAD_1 >= -180);
4028 complete_prev_size--;
4030 complete_post_size = complete_post_size + 1;
4031 counter_CAD_2 = counter_CAD_2 + delta_t * 360. * speed / 60.;
4032 }
while(counter_CAD_2 <= 180);
4033 complete_post_size--;
4035 complete_size = size + complete_prev_size + complete_post_size;
4037 complete_p_cyl = (
double*) malloc(complete_size *
sizeof(
double));
4038 complete_V_cyl = (
double*) malloc(complete_size *
sizeof(
double));
4039 complete_deform = (
double*) malloc(complete_size *
sizeof(
double));
4040 complete_CAD = (
double*) malloc(complete_size *
sizeof(
double));
4042 CALCULUS_OF_IMP_HP(complete_p_cyl, complete_CAD, p_cyl, V_cyl, complete_V_cyl, complete_deform, &WI_HP, &IMP_HP,
4043 complete_size, complete_prev_size, delta_t, speed, size, IVC, EVO, VTDC,
4044 Cylinder_capacity, PI, Piston_D, S, Crank_L, Connecting_Rod_L, E, Piston_Axis_D, Piston_Crown_H, M_Connecting_Rod,
4045 M_P_R_PA, C_ESteel, C_Mech_Defor, inlet_pres, exhaust_pres);
4047 CALCULUS_OF_MEAN_VARIABLES(p_cyl, T_cyl, dp_da_cyl, CAD, &pmax, &Tmax, &dp_da_max, &p_exit, &T_exit, size);
4050 mean_var_exit[0] = IMP_HP * 1e-5;
4051 mean_var_exit[1] = IMP_HP * Cylinder_capacity / (mfuel * 1e-6 * HP);
4054 mean_var_exit[2] = pmax * 1e-5;
4055 mean_var_exit[3] = Tmax;
4056 mean_var_exit[4] = p_exit * 1e-5;
4057 mean_var_exit[5] = T_exit;
4058 mean_var_exit[6] = dp_da_max * 1e-5;
4061 mean_var_exit[7] = mNOx_bowl[size - 3] * 1e6 / (1.587 * (mairIVC + mfuel * 1.e-3));
4063 mean_var_exit[8] = mO2_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4066 (*dataOUT).species_EVO[0] = mN2_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4067 (*dataOUT).species_EVO[1] = mO2_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4068 (*dataOUT).species_EVO[2] = mCO2_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4069 (*dataOUT).species_EVO[3] = mH2O_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4070 (*dataOUT).species_EVO[4] = mNOx_bowl[size - 3] / (mairIVC + mfuel * 1.e-3) / 1.587;
4071 (*dataOUT).species_EVO[5] = mCO_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4072 (*dataOUT).species_EVO[6] = mSOOT_bowl_A[size - 3] / (mairIVC + mfuel * 1.e-3);
4073 (*dataOUT).species_EVO[7] = mHC_bowl[size - 3] / (mairIVC + mfuel * 1.e-3);
4076 for(i = 0; i < 8; i++)
4077 Tot += (*dataOUT).species_EVO[i];
4078 for(i = 0; i < 8; i++)
4079 (*dataOUT).species_EVO[i] /= Tot;
4081 SOOT_EVO_A = mSOOT_bowl_A[size - 1] * 1.e6 / (mairIVC + mfuel * 1.e-3);
4084 mean_var_exit[10] = YSoot_to_FSN(SOOT_EVO_A);
4086 SOOT_EVO_C = (evol_mSoot[size - 1] + mSootCil[size - 1]) * 1.e6 / (mairIVC + mfuel * 1.e-3);
4089 mean_var_exit[9] = YSoot_to_FSN(SOOT_EVO_C);
4091 mean_var_exit[9] = 0.0930816 * pow(mean_var_exit[9], 2.68);
4093 if(mean_var_exit[9] > 10.)
4094 mean_var_exit[9] = 10.;
4096 heat_transfer[0] = 0;
4097 heat_transfer[1] = 0;
4098 heat_transfer[2] = 0;
4099 heat_transfer[3] = 0;
4101 for(counter = 0; counter < size - 1; counter++) {
4102 heat_transfer[0] = heat_transfer[0] + min(H_cooler[counter],
4103 H_cooler[counter + 1]) + pow(pow((H_cooler[counter + 1] - H_cooler[counter]) / 2, 2), 0.5);
4104 heat_transfer[1] = heat_transfer[1] + min(Qcylhead[counter],
4105 Qcylhead[counter + 1]) + pow(pow((Qcylhead[counter + 1] - Qcylhead[counter]) / 2, 2), 0.5);
4106 heat_transfer[2] = heat_transfer[2] + min(Qcyl[counter],
4107 Qcyl[counter + 1]) + pow(pow((Qcyl[counter + 1] - Qcyl[counter]) / 2, 2), 0.5);
4108 heat_transfer[3] = heat_transfer[3] + min(Qpis[counter],
4109 Qpis[counter + 1]) + pow(pow((Qpis[counter + 1] - Qpis[counter]) / 2, 2), 0.5);
4286 time_vector_exit = (
double*) malloc(CAI *
sizeof(
double));
4287 time_vector_exit[0] = -180 * 60 / (360. * speed);
4290 for(counter = 1; counter < CAI; counter++) {
4291 time_vector_exit[counter] = time_vector_exit[counter - 1] + 60 / ((CAI - 1) * speed);
4292 CAD_exit[counter] = CAD_exit[counter - 1] + 360. / (CAI - 1);
4295 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4296 for(counter = 0; counter < size; counter++) {
4297 vector_to_interpolate[counter] = p_cyl[counter];
4299 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4300 free(vector_to_interpolate);
4301 for(counter = 0; counter < CAI; counter++) {
4302 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4303 aux_vector[counter] = 0;
4305 p_cyl_exit[counter] = aux_vector[counter];
4309 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4310 for(counter = 0; counter < size; counter++) {
4311 vector_to_interpolate[counter] = dp_da_cyl[counter];
4314 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4315 free(vector_to_interpolate);
4316 for(counter = 0; counter < CAI; counter++) {
4317 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4318 aux_vector[counter] = 0;
4320 dp_da_cyl_exit[counter] = aux_vector[counter];
4324 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4325 for(counter = 0; counter < size; counter++) {
4326 vector_to_interpolate[counter] = T_cyl[counter];
4328 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4329 free(vector_to_interpolate);
4330 for(counter = 0; counter < CAI; counter++) {
4331 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4332 aux_vector[counter] = 0;
4334 T_cyl_exit[counter] = aux_vector[counter];
4338 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4339 for(counter = 0; counter < size; counter++) {
4340 vector_to_interpolate[counter] = HRF[counter];
4342 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4343 free(vector_to_interpolate);
4344 for(counter = 0; counter < CAI; counter++) {
4345 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4346 aux_vector[counter] = 0;
4348 HRF_exit[counter] = aux_vector[counter];
4352 if(aux_mfuel == 0.) {
4353 for(counter = 0; counter < CAI; counter++) {
4354 HRF_exit[counter] = 0.;
4359 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4360 for(counter = 0; counter < size; counter++) {
4361 vector_to_interpolate[counter] = HRF_PMX[counter];
4363 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4364 free(vector_to_interpolate);
4365 for(counter = 0; counter < CAI; counter++) {
4366 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4367 aux_vector[counter] = 0;
4369 (*dataOUT).HRF_PMX[counter] = aux_vector[counter];
4373 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4374 for(counter = 0; counter < size; counter++) {
4375 vector_to_interpolate[counter] = ROHR[counter];
4377 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4378 free(vector_to_interpolate);
4379 for(counter = 0; counter < CAI; counter++) {
4380 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4381 aux_vector[counter] = 0;
4383 ROHR_exit[counter] = aux_vector[counter];
4387 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4388 for(counter = 0; counter < size; counter++) {
4389 vector_to_interpolate[counter] = H_cooler[counter];
4391 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4392 free(vector_to_interpolate);
4393 for(counter = 0; counter < CAI; counter++) {
4394 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4395 aux_vector[counter] = 0;
4397 H_cooler_exit[counter] = aux_vector[counter] / (delta_t * speed * 360. / 60);
4401 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4402 for(counter = 0; counter < size; counter++) {
4403 vector_to_interpolate[counter] = dmf[counter];
4405 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4406 free(vector_to_interpolate);
4407 for(counter = 0; counter < CAI; counter++) {
4408 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4409 aux_vector[counter] = 0;
4411 injection_rate_exit[counter] = aux_vector[counter];
4415 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4416 for(counter = 0; counter < size; counter++) {
4417 vector_to_interpolate[counter] = acu_dmf[counter];
4420 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4421 free(vector_to_interpolate);
4422 for(counter = 0; counter < CAI; counter++) {
4423 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4424 aux_vector[counter] = 0;
4426 accum_injection_rate_exit[counter] = aux_vector[counter];
4430 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4431 for(counter = 0; counter < size; counter++) {
4432 vector_to_interpolate[counter] = evol_mSoot[counter];
4434 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4435 free(vector_to_interpolate);
4436 for(counter = 0; counter < CAI; counter++) {
4437 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4438 aux_vector[counter] = 0;
4440 (*dataOUT).evol_Soot[counter] = aux_vector[counter];
4444 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4445 for(counter = 0; counter < size; counter++) {
4446 vector_to_interpolate[counter] = mSootCil[counter];
4448 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4449 free(vector_to_interpolate);
4450 for(counter = 0; counter < CAI; counter++) {
4451 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4452 aux_vector[counter] = 0;
4454 (*dataOUT).evol_Soot_CIL[counter] = aux_vector[counter];
4458 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4459 for(counter = 0; counter < size; counter++) {
4460 vector_to_interpolate[counter] = evol_Radiacion[counter];
4462 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4463 free(vector_to_interpolate);
4464 for(counter = 0; counter < CAI; counter++) {
4465 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4466 aux_vector[counter] = 0;
4468 (*dataOUT).evol_Radiacion[counter] = aux_vector[counter];
4472 vector_to_interpolate = (
double*) malloc(size *
sizeof(
double));
4473 for(counter = 0; counter < size; counter++) {
4474 vector_to_interpolate[counter] = XLO[counter];
4476 FUNCTION_FOR_INTERPOLATION(aux_vector, time_vector_exit, CAD, vector_to_interpolate, CAI, size, speed);
4477 free(vector_to_interpolate);
4478 for(counter = 0; counter < CAI; counter++) {
4479 if((CAD_exit[counter] < IVC) || (CAD_exit[counter] > EVO)) {
4480 aux_vector[counter] = 0;
4482 (*dataOUT).evol_LOL[counter] = aux_vector[counter];
4518 free(dmtotal_Gfactor);
4522 free(evol_Radiacion);
4527 free(virtual_velocity);
4558 for(counter = 0; counter < 52; counter++)
4559 free(YNOeq[counter]);
4561 for(counter = 0; counter < 52; counter++)
4562 free(KdYNO[counter]);
4567 for(counter = 0; counter < inj_num; counter++)
4568 free(POI_IM[counter]);
4570 for(counter = 0; counter < inj_num; counter++)
4571 free(POC_IM[counter]);
4573 for(counter = 0; counter < inj_num; counter++)
4574 free(mfuel_ij_IM[counter]);
4576 for(counter = 0; counter < inj_num; counter++)
4577 free(mfuel_i_IM[counter]);
4580 for(counter = 0; counter < 15; counter++)
4581 free(elementcontrol[counter]);
4582 free(elementcontrol);
4584 free(mixture_correction);
4586 for(counter = 0; counter < inj_num; counter++)
4587 free(element[counter]);
4589 for(counter = 0; counter < inj_num; counter++)
4590 free(realelement[counter]);
4593 free(complete_p_cyl);
4594 free(complete_V_cyl);
4595 free(complete_deform);
4598 free(time_vector_exit);
4603 double YSoot_to_FSN(
double YSoot) {
4604 double x1 = 0., x2 = 0., x = 0., y = 0.;
4614 while((x2 - x1) > 0.01) {
4616 y = 1 / 1.2 / 0.405 * 4.95 * x * exp(0.38 * x);