Skip to content

Commit 07d8008

Browse files
henryleberrehenryleberre
committed
Add three 0-1D CheMFC example cases
+ nD_perfect_reactor + 1D_reactive_shocktube + 1D_inert_shocktube
1 parent 5feeb65 commit 07d8008

File tree

14 files changed

+621
-0
lines changed

14 files changed

+621
-0
lines changed

‎examples/1D_inert_shocktube/README.md‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
# 1D Multi-Component Inert Shock Tube
2+
3+
References:
4+
> P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Comput. & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
5+
6+
## Initial Condition
7+
8+
![Initial Condition](initial.png)
9+
10+
## Results
11+
12+
![Results](result.png)

‎examples/1D_inert_shocktube/case.py‎

Lines changed: 113 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,113 @@
1+
#!/usr/bin/env python3
2+
3+
# References:
4+
# + https://doi.org/10.1016/j.compfluid.2013.10.014: 4.3. Multi-component inert shock tube
5+
6+
import json
7+
import cantera as ct
8+
9+
ctfile = 'h2o2.yaml'
10+
sol_L = ct.Solution(ctfile)
11+
sol_L.TPX = 400, 8000, 'H2:2,O2:1,AR:7'
12+
sol_R = ct.Solution(ctfile)
13+
sol_R.TPX = 1200, 80000, 'H2:2,O2:1,AR:7'
14+
15+
L = 0.10
16+
Nx = 400
17+
dx = L / Nx
18+
dt = 5e-9
19+
Tend = 40e-6
20+
21+
chemistry = True
22+
NT = int(Tend / dt)
23+
SAVE_COUNT = 200
24+
NS = NT // SAVE_COUNT
25+
26+
case = {
27+
# Logistics ================================================================
28+
'run_time_info' : 'T',
29+
# ==========================================================================
30+
31+
# Computational Domain Parameters ==========================================
32+
'x_domain%beg' : -L/2,
33+
'x_domain%end' : +L/2,
34+
'm' : Nx,
35+
'n' : 0,
36+
'p' : 0,
37+
'dt' : float(dt),
38+
't_step_start' : 0,
39+
't_step_stop' : NT,
40+
't_step_save' : NS,
41+
't_step_print' : NS,
42+
'parallel_io' : 'F',
43+
# ==========================================================================
44+
45+
# Simulation Algorithm Parameters ==========================================
46+
'model_eqns' : 2,
47+
'num_fluids' : 1,
48+
'num_patches' : 2,
49+
'mpp_lim' : 'F',
50+
'mixture_err' : 'F',
51+
'time_stepper' : 3,
52+
'weno_order' : 5,
53+
'weno_eps' : 1E-16,
54+
'weno_avg' : 'F',
55+
'mapped_weno' : 'T',
56+
'mp_weno' : 'T',
57+
'riemann_solver' : 1,
58+
'wave_speeds' : 1,
59+
'avg_state' : 2,
60+
'bc_x%beg' :-2,
61+
'bc_x%end' :-3,
62+
# ==========================================================================
63+
64+
# Chemistry ================================================================
65+
'chemistry' : 'F' if not chemistry else 'T',
66+
'chem_params%advection' : 'T',
67+
'chem_params%diffusion' : 'F',
68+
'chem_params%reactions' : 'T',
69+
# ==========================================================================
70+
71+
# Formatted Database Files Structure Parameters ============================
72+
'format' : 1,
73+
'precision' : 2,
74+
'prim_vars_wrt' : 'T',
75+
# ==========================================================================
76+
77+
# ==========================================================================
78+
'patch_icpp(1)%geometry' : 1,
79+
'patch_icpp(1)%x_centroid' : -L/4,
80+
'patch_icpp(1)%length_x' : L/2,
81+
'patch_icpp(1)%vel(1)' : 0,
82+
'patch_icpp(1)%pres' : sol_L.P,
83+
'patch_icpp(1)%alpha(1)' : 1,
84+
'patch_icpp(1)%alpha_rho(1)' : sol_L.density,
85+
# ==========================================================================
86+
87+
# ==========================================================================
88+
'patch_icpp(2)%geometry' : 1,
89+
'patch_icpp(2)%x_centroid' : L/4,
90+
'patch_icpp(2)%length_x' : L/2,
91+
'patch_icpp(2)%vel(1)' : 0,
92+
'patch_icpp(2)%pres' : sol_R.P,
93+
'patch_icpp(2)%alpha(1)' : 1,
94+
'patch_icpp(2)%alpha_rho(1)' : sol_R.density,
95+
# ==========================================================================
96+
97+
# Fluids Physical Parameters ===============================================
98+
'fluid_pp(1)%gamma' : 1.0E+00/(4.4E+00-1.0E+00),
99+
'fluid_pp(1)%pi_inf' : 0,
100+
# ==========================================================================
101+
102+
# Chemistry ================================================================
103+
'cantera_file' : ctfile,
104+
# ==========================================================================
105+
}
106+
107+
if chemistry:
108+
for i in range(len(sol_L.Y)):
109+
case[f'patch_icpp(1)%Y({i+1})'] = sol_L.Y[i]
110+
case[f'patch_icpp(2)%Y({i+1})'] = sol_R.Y[i]
111+
112+
if __name__ == '__main__':
113+
print(json.dumps(case))

‎examples/1D_inert_shocktube/initial.png‎

58.2 KB
Loading

‎examples/1D_inert_shocktube/result.png‎

81.1 KB
Loading

‎examples/1D_inert_shocktube/viz.py‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
import mfc.viz
2+
3+
import os
4+
5+
import subprocess
6+
import seaborn as sns
7+
import matplotlib.pyplot as plt
8+
from tqdm import tqdm
9+
10+
from case import sol_L as sol
11+
12+
case = mfc.viz.Case(".")
13+
14+
os.makedirs("viz", exist_ok=True)
15+
16+
#sns.set_theme(style=mfc.viz.generate_cpg_style())
17+
18+
Y_VARS = ["H2", "O2", "H2O", "N2"]
19+
20+
variables = [
21+
("rho", "prim.1"),
22+
("u_x", "prim.2"),
23+
("p", "prim.3"),
24+
("E", "cons.3"),
25+
*[(f"Y_{name}", f"prim.{5 + sol.species_index(name)}") for name in Y_VARS],
26+
("T", "prim.15"),
27+
]
28+
29+
for variable in tqdm(variables, desc="Loading Variables"):
30+
case.load_variable(*variable)
31+
32+
for step in tqdm(case.get_timesteps(), desc="Rendering Frames"):
33+
fig, axes = plt.subplots(2, 3, figsize=(16, 9))
34+
35+
def pad_ylim(ylim, pad=0.1):
36+
return (ylim[0] - pad*(ylim[1] - ylim[0]), ylim[1] + pad*(ylim[1] - ylim[0]))
37+
38+
case.plot_step(step, "rho", ax=axes[0, 0])
39+
axes[0, 0].set_ylim(*pad_ylim(case.get_minmax_time("rho")))
40+
axes[0, 0].set_ylabel("$\\rho$")
41+
case.plot_step(step, "u_x", ax=axes[0, 1])
42+
axes[0, 1].set_ylim(*pad_ylim(case.get_minmax_time("u_x")))
43+
axes[0, 1].set_ylabel("$u_x$")
44+
case.plot_step(step, "p", ax=axes[1, 0])
45+
axes[1, 0].set_ylim(*pad_ylim(case.get_minmax_time("p")))
46+
axes[1, 0].set_ylabel("$p$")
47+
for y in Y_VARS:
48+
case.plot_step(step, f"Y_{y}", ax=axes[1, 1], label=y)
49+
axes[1, 1].set_ylim(0, 1.1*max(case.get_minmax_time(f"Y_{y}")[1] for y in Y_VARS))
50+
axes[1, 1].set_ylabel("$Y_k$")
51+
case.plot_step(step, "T", ax=axes[1, 2])
52+
axes[1, 2].set_ylim(*pad_ylim(case.get_minmax_time("T")))
53+
axes[1, 2].set_ylabel("$T$")
54+
case.plot_step(step, "E", ax=axes[0, 2])
55+
axes[0, 2].set_ylim(*pad_ylim(case.get_minmax_time("E")))
56+
axes[0, 2].set_ylabel("$E$")
57+
58+
plt.tight_layout()
59+
plt.savefig(f"viz/{step:06d}.png")
60+
plt.close()
61+
62+
subprocess.run([
63+
"ffmpeg", "-y", "-framerate", "60", "-pattern_type", "glob", "-i",
64+
"viz/*.png", "-c:v", "libx264", "-pix_fmt", "yuv420p", "viz.mp4"
65+
])

‎examples/1D_reactive_shocktube/README.md‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
# 1D Multi-Component Reactive Shock Tube
2+
3+
References:
4+
> P. J. Martínez Ferrer, R. Buttay, G. Lehnasch, and A. Mura, “A detailed verification procedure for compressible reactive multicomponent Navier–Stokes solvers”, Comput. & Fluids, vol. 89, pp. 88–110, Jan. 2014. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.compfluid.2013.10.014
5+
6+
> H. Chen, C. Si, Y. Wu, H. Hu, and Y. Zhu, “Numerical investigation of the effect of equivalence ratio on the propagation characteristics and performance of rotating detonation engine”, Int. J. Hydrogen Energy, Mar. 2023. Accessed: Oct. 13, 2024. [Online]. Available: https://doi.org/10.1016/j.ijhydene.2023.03.190
7+
8+
## Initial Condition
9+
10+
![Initial Condition](initial.png)
11+
12+
## Results
13+
14+
![Results](result.png)

‎examples/1D_reactive_shocktube/case.py‎

Lines changed: 114 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,114 @@
1+
#!/usr/bin/env python3
2+
3+
# References:
4+
# + https://doi.org/10.1016/j.ijhydene.2023.03.190: Verification of numerical method
5+
# + https://doi.org/10.1016/j.compfluid.2013.10.014: 4.7. Multi-species reactive shock tube
6+
7+
import json
8+
import cantera as ct
9+
10+
ctfile = 'h2o2.yaml'
11+
sol_L = ct.Solution(ctfile)
12+
sol_L.DPX = 0.072, 7173, 'H2:2,O2:1,AR:7'
13+
14+
sol_R = ct.Solution(ctfile)
15+
sol_R.DPX = 0.18075, 35594, 'H2:2,O2:1,AR:7'
16+
17+
u_l = 0
18+
u_r = -487.34
19+
20+
L = 0.12
21+
Nx = 800
22+
dx = L/Nx
23+
dt = dx/abs(u_r)*0.01
24+
Tend=230e-6
25+
26+
NT=int(Tend/dt)
27+
SAVE_COUNT=100
28+
NS=NT//SAVE_COUNT
29+
30+
chemistry = True
31+
32+
case = {
33+
# Logistics ================================================================
34+
'run_time_info' : 'T',
35+
# ==========================================================================
36+
37+
# Computational Domain Parameters ==========================================
38+
'x_domain%beg' : 0,
39+
'x_domain%end' : L,
40+
'm' : Nx,
41+
'n' : 0,
42+
'p' : 0,
43+
'dt' : float(dt),
44+
't_step_start' : 0,
45+
't_step_stop' : NT,
46+
't_step_save' : NS,
47+
't_step_print' : NS,
48+
'parallel_io' : 'F',
49+
50+
# Simulation Algorithm Parameters ==========================================
51+
'model_eqns' : 2,
52+
'num_fluids' : 1,
53+
'num_patches' : 2,
54+
'mpp_lim' : 'F',
55+
'mixture_err' : 'F',
56+
'time_stepper' : 3,
57+
'weno_order' : 5,
58+
'weno_eps' : 1E-16,
59+
'weno_avg' : 'F',
60+
'mapped_weno' : 'T',
61+
'mp_weno' : 'T',
62+
'riemann_solver' : 1,
63+
'wave_speeds' : 1,
64+
'avg_state' : 2,
65+
'bc_x%beg' :-2,
66+
'bc_x%end' :-3,
67+
68+
# Chemistry ================================================================
69+
'chemistry' : 'F' if not chemistry else 'T',
70+
'chem_params%advection' : 'T',
71+
'chem_params%diffusion' : 'F',
72+
'chem_params%reactions' : 'T',
73+
'cantera_file' : ctfile,
74+
# ==========================================================================
75+
76+
# Formatted Database Files Structure Parameters ============================
77+
'format' : 1,
78+
'precision' : 2,
79+
'prim_vars_wrt' : 'T',
80+
# ==========================================================================
81+
82+
# ==========================================================================
83+
'patch_icpp(1)%geometry' : 1,
84+
'patch_icpp(1)%x_centroid' : L/4,
85+
'patch_icpp(1)%length_x' : L/2,
86+
'patch_icpp(1)%vel(1)' : u_l,
87+
'patch_icpp(1)%pres' : sol_L.P,
88+
'patch_icpp(1)%alpha(1)' : 1,
89+
'patch_icpp(1)%alpha_rho(1)' : sol_L.density,
90+
# ==========================================================================
91+
92+
# ==========================================================================
93+
'patch_icpp(2)%geometry' : 1,
94+
'patch_icpp(2)%x_centroid' : 3*L/4,
95+
'patch_icpp(2)%length_x' : L/2,
96+
'patch_icpp(2)%vel(1)' : u_r,
97+
'patch_icpp(2)%pres' : sol_R.P,
98+
'patch_icpp(2)%alpha(1)' : 1,
99+
'patch_icpp(2)%alpha_rho(1)' : sol_R.density,
100+
# ==========================================================================
101+
102+
# Fluids Physical Parameters ===============================================
103+
'fluid_pp(1)%gamma' : 1.0E+00/(4.4E+00-1.0E+00),
104+
'fluid_pp(1)%pi_inf' : 0,
105+
# ==========================================================================
106+
}
107+
108+
if chemistry:
109+
for i in range(len(sol_L.Y)):
110+
case[f'patch_icpp(1)%Y({i+1})'] = sol_L.Y[i]
111+
case[f'patch_icpp(2)%Y({i+1})'] = sol_R.Y[i]
112+
113+
if __name__ == '__main__':
114+
print(json.dumps(case))

‎examples/1D_reactive_shocktube/initial.png‎

48.5 KB
Loading

‎examples/1D_reactive_shocktube/result.png‎

83.1 KB
Loading

‎examples/1D_reactive_shocktube/viz.py‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
import mfc.viz
2+
3+
import os
4+
5+
import subprocess
6+
import seaborn as sns
7+
import matplotlib.pyplot as plt
8+
from tqdm import tqdm
9+
10+
from case import sol_L as sol
11+
12+
case = mfc.viz.Case(".")
13+
14+
os.makedirs("viz", exist_ok=True)
15+
16+
sns.set_theme(style=mfc.viz.generate_cpg_style())
17+
18+
Y_VARS = ["H2", "O2", "H2O", "N2"]
19+
20+
variables = [
21+
("rho", "prim.1"),
22+
("u_x", "prim.2"),
23+
("p", "prim.3"),
24+
("E", "cons.3"),
25+
*[(f"Y_{name}", f"prim.{5 + sol.species_index(name)}") for name in Y_VARS],
26+
("T", "prim.15"),
27+
]
28+
29+
for variable in tqdm(variables, desc="Loading Variables"):
30+
case.load_variable(*variable)
31+
32+
for step in tqdm(case.get_timesteps(), desc="Rendering Frames"):
33+
fig, axes = plt.subplots(2, 3, figsize=(16, 9))
34+
35+
def pad_ylim(ylim, pad=0.1):
36+
return (ylim[0] - pad*(ylim[1] - ylim[0]), ylim[1] + pad*(ylim[1] - ylim[0]))
37+
38+
case.plot_step(step, "rho", ax=axes[0, 0])
39+
axes[0, 0].set_ylim(*pad_ylim(case.get_minmax_time("rho")))
40+
axes[0, 0].set_ylabel("$\\rho$")
41+
case.plot_step(step, "u_x", ax=axes[0, 1])
42+
axes[0, 1].set_ylim(*pad_ylim(case.get_minmax_time("u_x")))
43+
axes[0, 1].set_ylabel("$u_x$")
44+
case.plot_step(step, "p", ax=axes[1, 0])
45+
axes[1, 0].set_ylim(*pad_ylim(case.get_minmax_time("p")))
46+
axes[1, 0].set_ylabel("$p$")
47+
for y in Y_VARS:
48+
case.plot_step(step, f"Y_{y}", ax=axes[1, 1], label=y)
49+
axes[1, 1].set_ylim(0, 1.1*max(case.get_minmax_time(f"Y_{y}")[1] for y in Y_VARS))
50+
axes[1, 1].set_ylabel("$Y_k$")
51+
case.plot_step(step, "T", ax=axes[1, 2])
52+
axes[1, 2].set_ylim(*pad_ylim(case.get_minmax_time("T")))
53+
axes[1, 2].set_ylabel("$T$")
54+
case.plot_step(step, "E", ax=axes[0, 2])
55+
axes[0, 2].set_ylim(*pad_ylim(case.get_minmax_time("E")))
56+
axes[0, 2].set_ylabel("$E$")
57+
58+
plt.tight_layout()
59+
plt.savefig(f"viz/{step:06d}.png")
60+
plt.close()
61+
62+
subprocess.run([
63+
"ffmpeg", "-y", "-framerate", "60", "-pattern_type", "glob", "-i",
64+
"viz/*.png", "-c:v", "libx264", "-pix_fmt", "yuv420p", "viz.mp4"
65+
])

0 commit comments

Comments
 (0)