|
| 1 | +import csv |
| 2 | +import cantera as ct |
| 3 | +from tqdm import tqdm |
| 4 | + |
| 5 | +import mfc.viz |
| 6 | +from case import dt, NS, Tend, SAVE_COUNT, sol |
| 7 | + |
| 8 | + |
| 9 | +case = mfc.viz.Case(".", dt) |
| 10 | + |
| 11 | +for name in tqdm(sol.species_names, desc="Loading Variables"): |
| 12 | + case.load_variable(f"Y_{name}", f"prim.{5 + sol.species_index(name)}") |
| 13 | +case.load_variable("rho", "prim.1") |
| 14 | + |
| 15 | +time_save = Tend/SAVE_COUNT |
| 16 | + |
| 17 | +oh_idx = sol.species_index("OH") |
| 18 | +def generate_ct_saves() -> tuple: |
| 19 | + reactor = ct.IdealGasReactor(sol) |
| 20 | + reactor_network = ct.ReactorNet([reactor]) |
| 21 | + |
| 22 | + ct_time = 0.0 |
| 23 | + ct_ts, ct_Ys, ct_rhos = [0.0], [reactor.thermo.Y], [reactor.thermo.density] |
| 24 | + |
| 25 | + while ct_time < Tend: |
| 26 | + reactor_network.advance(ct_time + time_save) |
| 27 | + ct_time += time_save |
| 28 | + ct_ts.append(ct_time) |
| 29 | + ct_Ys.append(reactor.thermo.Y) |
| 30 | + ct_rhos.append(reactor.thermo.density) |
| 31 | + |
| 32 | + return ct_ts, ct_Ys, ct_rhos |
| 33 | + |
| 34 | +ct_ts, ct_Ys, ct_rhos = generate_ct_saves() |
| 35 | + |
| 36 | +with open("mfc.csv", "w") as f: |
| 37 | + writer = csv.writer(f) |
| 38 | + keys = ["t"] + list(set(case.get_data()[0].keys()) - set(["x"])) |
| 39 | + writer.writerow(keys) |
| 40 | + for i, t_step in enumerate(sorted(case.get_timesteps())): |
| 41 | + t = t_step * dt |
| 42 | + row = [t] + [case.get_data()[t_step][key][0] for key in keys[1:]] |
| 43 | + writer.writerow(row) |
| 44 | + |
| 45 | +with open("cantera.csv", "w") as f: |
| 46 | + writer = csv.writer(f) |
| 47 | + keys = ["t"] + [f"Y_{_}" for _ in list(sol.species_names)] + ["rho"] |
| 48 | + writer.writerow(keys) |
| 49 | + for step in range(len(ct_ts)): |
| 50 | + row = [ct_ts[step]] + [ ct_Ys[step][i] for i in range(len(sol.species_names)) ] + [ct_rhos[step]] |
| 51 | + print([ct_ts[step]], row) |
| 52 | + writer.writerow(row) |
| 53 | + |
| 54 | +def find_induction_time(ts: list, Ys: list, rhos: list) -> float: |
| 55 | + for t, y, rho in zip(ts, Ys, rhos): |
| 56 | + if (y * rho / sol.molecular_weights[oh_idx]) >= 1e-6: |
| 57 | + return t |
| 58 | + |
| 59 | + return None |
| 60 | + |
| 61 | +skinner_induction_time = 0.052e-3 |
| 62 | +ct_induction_time = find_induction_time( |
| 63 | + ct_ts, |
| 64 | + [y[oh_idx] for y in ct_Ys], |
| 65 | + [rho for rho in ct_rhos] |
| 66 | +) |
| 67 | +mfc_induction_time = find_induction_time( |
| 68 | + sorted(case.get_timestamps()), |
| 69 | + [ case.get_data()[step]["Y_OH"][0] for step in sorted(case.get_timesteps()) ], |
| 70 | + [ case.get_data()[step]["rho"][0] for step in sorted(case.get_timesteps()) ] |
| 71 | +) |
| 72 | + |
| 73 | +print("Induction Times ([OH] >= 1e-6):") |
| 74 | +print(f" + Skinner et al.: {skinner_induction_time} s") |
| 75 | +print(f" + Cantera: {ct_induction_time} s") |
| 76 | +print(f" + (Che)MFC: {mfc_induction_time} s") |
0 commit comments