Add three 0-1D CheMFC example cases

+ nD_perfect_reactor + 1D_reactive_shocktube + 1D_inert_shocktube
MFlowCode · Oct 16, 2024 · 37046b5 · 37046b5
1 parent 96bb43d
commit 37046b5
Show file tree

Hide file tree

Showing 14 changed files with 621 additions and 0 deletions.
diff --git a/examples/1D_inert_shocktube/README.md b/examples/1D_inert_shocktube/README.md
@@ -0,0 +1,12 @@
+# 1D Multi-Component Inert Shock Tube
+
+References:
+> 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
+
+## Initial Condition
+
+![Initial Condition](initial.png)
+
+## Results
+
+![Results](result.png)
diff --git a/examples/1D_inert_shocktube/case.py b/examples/1D_inert_shocktube/case.py
@@ -0,0 +1,113 @@
+#!/usr/bin/env python3
+
+# References:
+# + https://doi.org/10.1016/j.compfluid.2013.10.014: 4.3. Multi-component inert shock tube
+
+import json
+import cantera as ct
+
+ctfile    = 'h2o2.yaml'
+sol_L     = ct.Solution(ctfile)
+sol_L.TPX =  400,  8000, 'H2:2,O2:1,AR:7'
+sol_R     = ct.Solution(ctfile)
+sol_R.TPX = 1200, 80000, 'H2:2,O2:1,AR:7'
+
+L    = 0.10
+Nx   = 400
+dx   = L / Nx
+dt   = 5e-9
+Tend = 40e-6
+
+chemistry  = True
+NT         = int(Tend / dt)
+SAVE_COUNT = 200
+NS         = NT // SAVE_COUNT
+
+case = {
+    # Logistics ================================================================
+    'run_time_info'                : 'T',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    'x_domain%beg'                 : -L/2,
+    'x_domain%end'                 : +L/2,
+    'm'                            : Nx,
+    'n'                            : 0,
+    'p'                            : 0,
+    'dt'                           : float(dt),
+    't_step_start'                 : 0,
+    't_step_stop'                  : NT,
+    't_step_save'                  : NS,
+    't_step_print'                 : NS,
+    'parallel_io'                  : 'F',
+    # ==========================================================================
+
+    # Simulation Algorithm Parameters ==========================================
+    'model_eqns'                   : 2,
+    'num_fluids'                   : 1,
+    'num_patches'                  : 2,
+    'mpp_lim'                      : 'F',
+    'mixture_err'                  : 'F',
+    'time_stepper'                 : 3,
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1E-16,
+    'weno_avg'                     : 'F',
+    'mapped_weno'                  : 'T',
+    'mp_weno'                      : 'T',
+    'riemann_solver'               : 1,
+    'wave_speeds'                  : 1,
+    'avg_state'                    : 2,
+    'bc_x%beg'                     :-2,
+    'bc_x%end'                     :-3,
+    # ==========================================================================
+
+    # Chemistry ================================================================
+    'chemistry'                    : 'F' if not chemistry else 'T',
+    'chem_params%advection'        : 'T',
+    'chem_params%diffusion'        : 'F',
+    'chem_params%reactions'        : 'T',
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    'format'                       : 1,
+    'precision'                    : 2,
+    'prim_vars_wrt'                : 'T',
+    # ==========================================================================
+
+    # ==========================================================================
+    'patch_icpp(1)%geometry'       : 1,
+    'patch_icpp(1)%x_centroid'     : -L/4,
+    'patch_icpp(1)%length_x'       : L/2,
+    'patch_icpp(1)%vel(1)'         : 0,
+    'patch_icpp(1)%pres'           : sol_L.P,
+    'patch_icpp(1)%alpha(1)'       : 1,
+    'patch_icpp(1)%alpha_rho(1)'   : sol_L.density,
+    # ==========================================================================
+
+    # ==========================================================================
+    'patch_icpp(2)%geometry'       : 1,
+    'patch_icpp(2)%x_centroid'     : L/4,
+    'patch_icpp(2)%length_x'       : L/2,
+    'patch_icpp(2)%vel(1)'         : 0,
+    'patch_icpp(2)%pres'           : sol_R.P,
+    'patch_icpp(2)%alpha(1)'       : 1,
+    'patch_icpp(2)%alpha_rho(1)'   : sol_R.density,
+    # ==========================================================================
+
+    # Fluids Physical Parameters ===============================================
+    'fluid_pp(1)%gamma'            : 1.0E+00/(4.4E+00-1.0E+00),
+    'fluid_pp(1)%pi_inf'           : 0,
+    # ==========================================================================
+
+    # Chemistry ================================================================
+    'cantera_file'                 : ctfile,
+    # ==========================================================================
+}
+
+if chemistry:
+    for i in range(len(sol_L.Y)):
+        case[f'patch_icpp(1)%Y({i+1})'] = sol_L.Y[i]
+        case[f'patch_icpp(2)%Y({i+1})'] = sol_R.Y[i]
+
+if __name__ == '__main__':
+    print(json.dumps(case))
diff --git a/examples/1D_inert_shocktube/initial.png b/examples/1D_inert_shocktube/initial.png
diff --git a/examples/1D_inert_shocktube/result.png b/examples/1D_inert_shocktube/result.png
diff --git a/examples/1D_inert_shocktube/viz.py b/examples/1D_inert_shocktube/viz.py
@@ -0,0 +1,65 @@
+import mfc.viz
+
+import os
+
+import subprocess
+import seaborn as sns
+import matplotlib.pyplot as plt
+from tqdm import tqdm
+
+from case import sol_L as sol
+
+case = mfc.viz.Case(".")
+
+os.makedirs("viz", exist_ok=True)
+
+#sns.set_theme(style=mfc.viz.generate_cpg_style())
+
+Y_VARS = ["H2", "O2", "H2O", "N2"]
+
+variables = [
+    ("rho", "prim.1"),
+    ("u_x", "prim.2"),
+    ("p",   "prim.3"),
+    ("E",   "cons.3"),
+    *[(f"Y_{name}", f"prim.{5 + sol.species_index(name)}") for name in Y_VARS],
+    ("T",   "prim.15"),
+]
+
+for variable in tqdm(variables, desc="Loading Variables"):
+    case.load_variable(*variable)
+
+for step in tqdm(case.get_timesteps(), desc="Rendering Frames"):
+    fig, axes = plt.subplots(2, 3, figsize=(16, 9))
+
+    def pad_ylim(ylim, pad=0.1):
+        return (ylim[0] - pad*(ylim[1] - ylim[0]), ylim[1] + pad*(ylim[1] - ylim[0]))
+
+    case.plot_step(step, "rho", ax=axes[0, 0])
+    axes[0, 0].set_ylim(*pad_ylim(case.get_minmax_time("rho")))
+    axes[0, 0].set_ylabel("$\\rho$")
+    case.plot_step(step, "u_x", ax=axes[0, 1])
+    axes[0, 1].set_ylim(*pad_ylim(case.get_minmax_time("u_x")))
+    axes[0, 1].set_ylabel("$u_x$")
+    case.plot_step(step, "p",   ax=axes[1, 0])
+    axes[1, 0].set_ylim(*pad_ylim(case.get_minmax_time("p")))
+    axes[1, 0].set_ylabel("$p$")
+    for y in Y_VARS:
+        case.plot_step(step, f"Y_{y}", ax=axes[1, 1], label=y)
+    axes[1, 1].set_ylim(0, 1.1*max(case.get_minmax_time(f"Y_{y}")[1] for y in Y_VARS))
+    axes[1, 1].set_ylabel("$Y_k$")
+    case.plot_step(step, "T",   ax=axes[1, 2])
+    axes[1, 2].set_ylim(*pad_ylim(case.get_minmax_time("T")))
+    axes[1, 2].set_ylabel("$T$")
+    case.plot_step(step, "E",   ax=axes[0, 2])
+    axes[0, 2].set_ylim(*pad_ylim(case.get_minmax_time("E")))
+    axes[0, 2].set_ylabel("$E$")
+
+    plt.tight_layout()
+    plt.savefig(f"viz/{step:06d}.png")
+    plt.close()
+
+subprocess.run([
+    "ffmpeg", "-y", "-framerate", "60", "-pattern_type", "glob", "-i",
+    "viz/*.png", "-c:v", "libx264", "-pix_fmt", "yuv420p", "viz.mp4"
+])
diff --git a/examples/1D_reactive_shocktube/README.md b/examples/1D_reactive_shocktube/README.md
@@ -0,0 +1,14 @@
+# 1D Multi-Component Reactive Shock Tube
+
+References:
+> 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
+
+> 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
+
+## Initial Condition
+
+![Initial Condition](initial.png)
+
+## Results
+
+![Results](result.png)
diff --git a/examples/1D_reactive_shocktube/case.py b/examples/1D_reactive_shocktube/case.py
@@ -0,0 +1,114 @@
+#!/usr/bin/env python3
+
+# References:
+# + https://doi.org/10.1016/j.ijhydene.2023.03.190:  Verification of numerical method
+# + https://doi.org/10.1016/j.compfluid.2013.10.014: 4.7. Multi-species reactive shock tube
+
+import json
+import cantera as ct
+
+ctfile    = 'h2o2.yaml'
+sol_L     = ct.Solution(ctfile)
+sol_L.DPX =   0.072,  7173, 'H2:2,O2:1,AR:7'
+
+sol_R     = ct.Solution(ctfile)
+sol_R.DPX = 0.18075, 35594, 'H2:2,O2:1,AR:7'
+
+u_l = 0
+u_r = -487.34
+
+L  = 0.12
+Nx = 800
+dx = L/Nx
+dt = dx/abs(u_r)*0.01
+Tend=230e-6
+
+NT=int(Tend/dt)
+SAVE_COUNT=100
+NS=NT//SAVE_COUNT
+
+chemistry = True
+
+case = {
+    # Logistics ================================================================
+    'run_time_info'                : 'T',
+    # ==========================================================================
+
+    # Computational Domain Parameters ==========================================
+    'x_domain%beg'                 : 0,
+    'x_domain%end'                 : L,
+    'm'                            : Nx,
+    'n'                            : 0,
+    'p'                            : 0,
+    'dt'                           : float(dt),
+    't_step_start'                 : 0,
+    't_step_stop'                  : NT,
+    't_step_save'                  : NS,
+    't_step_print'                 : NS,
+    'parallel_io'                  : 'F',
+
+    # Simulation Algorithm Parameters ==========================================
+    'model_eqns'                   : 2,
+    'num_fluids'                   : 1,
+    'num_patches'                  : 2,
+    'mpp_lim'                      : 'F',
+    'mixture_err'                  : 'F',
+    'time_stepper'                 : 3,
+    'weno_order'                   : 5,
+    'weno_eps'                     : 1E-16,
+    'weno_avg'                     : 'F',
+    'mapped_weno'                  : 'T',
+    'mp_weno'                      : 'T',
+    'riemann_solver'               : 1,
+    'wave_speeds'                  : 1,
+    'avg_state'                    : 2,
+    'bc_x%beg'                     :-2,
+    'bc_x%end'                     :-3,
+
+    # Chemistry ================================================================
+    'chemistry'                    : 'F' if not chemistry else 'T',
+    'chem_params%advection'        : 'T',
+    'chem_params%diffusion'        : 'F',
+    'chem_params%reactions'        : 'T',
+    'cantera_file'                 : ctfile,
+    # ==========================================================================
+
+    # Formatted Database Files Structure Parameters ============================
+    'format'                       : 1,
+    'precision'                    : 2,
+    'prim_vars_wrt'                : 'T',
+    # ==========================================================================
+
+    # ==========================================================================
+    'patch_icpp(1)%geometry'       : 1,
+    'patch_icpp(1)%x_centroid'     : L/4,
+    'patch_icpp(1)%length_x'       : L/2,
+    'patch_icpp(1)%vel(1)'         : u_l,
+    'patch_icpp(1)%pres'           : sol_L.P,
+    'patch_icpp(1)%alpha(1)'       : 1,
+    'patch_icpp(1)%alpha_rho(1)'   : sol_L.density,
+    # ==========================================================================
+
+    # ==========================================================================
+    'patch_icpp(2)%geometry'       : 1,
+    'patch_icpp(2)%x_centroid'     : 3*L/4,
+    'patch_icpp(2)%length_x'       : L/2,
+    'patch_icpp(2)%vel(1)'         : u_r,
+    'patch_icpp(2)%pres'           : sol_R.P,
+    'patch_icpp(2)%alpha(1)'       : 1,
+    'patch_icpp(2)%alpha_rho(1)'   : sol_R.density,
+    # ==========================================================================
+
+    # Fluids Physical Parameters ===============================================
+    'fluid_pp(1)%gamma'            : 1.0E+00/(4.4E+00-1.0E+00),
+    'fluid_pp(1)%pi_inf'           : 0,
+    # ==========================================================================
+}
+
+if chemistry:
+    for i in range(len(sol_L.Y)):
+        case[f'patch_icpp(1)%Y({i+1})'] = sol_L.Y[i]
+        case[f'patch_icpp(2)%Y({i+1})'] = sol_R.Y[i]
+
+if __name__ == '__main__':
+    print(json.dumps(case))
diff --git a/examples/1D_reactive_shocktube/initial.png b/examples/1D_reactive_shocktube/initial.png
diff --git a/examples/1D_reactive_shocktube/result.png b/examples/1D_reactive_shocktube/result.png
diff --git a/examples/1D_reactive_shocktube/viz.py b/examples/1D_reactive_shocktube/viz.py
@@ -0,0 +1,65 @@
+import mfc.viz
+
+import os
+
+import subprocess
+import seaborn as sns
+import matplotlib.pyplot as plt
+from tqdm import tqdm
+
+from case import sol_L as sol
+
+case = mfc.viz.Case(".")
+
+os.makedirs("viz", exist_ok=True)
+
+sns.set_theme(style=mfc.viz.generate_cpg_style())
+
+Y_VARS = ["H2", "O2", "H2O", "N2"]
+
+variables = [
+    ("rho", "prim.1"),
+    ("u_x", "prim.2"),
+    ("p",   "prim.3"),
+    ("E",   "cons.3"),
+    *[(f"Y_{name}", f"prim.{5 + sol.species_index(name)}") for name in Y_VARS],
+    ("T",   "prim.15"),
+]
+
+for variable in tqdm(variables, desc="Loading Variables"):
+    case.load_variable(*variable)
+
+for step in tqdm(case.get_timesteps(), desc="Rendering Frames"):
+    fig, axes = plt.subplots(2, 3, figsize=(16, 9))
+
+    def pad_ylim(ylim, pad=0.1):
+        return (ylim[0] - pad*(ylim[1] - ylim[0]), ylim[1] + pad*(ylim[1] - ylim[0]))
+
+    case.plot_step(step, "rho", ax=axes[0, 0])
+    axes[0, 0].set_ylim(*pad_ylim(case.get_minmax_time("rho")))
+    axes[0, 0].set_ylabel("$\\rho$")
+    case.plot_step(step, "u_x", ax=axes[0, 1])
+    axes[0, 1].set_ylim(*pad_ylim(case.get_minmax_time("u_x")))
+    axes[0, 1].set_ylabel("$u_x$")
+    case.plot_step(step, "p",   ax=axes[1, 0])
+    axes[1, 0].set_ylim(*pad_ylim(case.get_minmax_time("p")))
+    axes[1, 0].set_ylabel("$p$")
+    for y in Y_VARS:
+        case.plot_step(step, f"Y_{y}", ax=axes[1, 1], label=y)
+    axes[1, 1].set_ylim(0, 1.1*max(case.get_minmax_time(f"Y_{y}")[1] for y in Y_VARS))
+    axes[1, 1].set_ylabel("$Y_k$")
+    case.plot_step(step, "T",   ax=axes[1, 2])
+    axes[1, 2].set_ylim(*pad_ylim(case.get_minmax_time("T")))
+    axes[1, 2].set_ylabel("$T$")
+    case.plot_step(step, "E",   ax=axes[0, 2])
+    axes[0, 2].set_ylim(*pad_ylim(case.get_minmax_time("E")))
+    axes[0, 2].set_ylabel("$E$")
+
+    plt.tight_layout()
+    plt.savefig(f"viz/{step:06d}.png")
+    plt.close()
+
+subprocess.run([
+    "ffmpeg", "-y", "-framerate", "60", "-pattern_type", "glob", "-i",
+    "viz/*.png", "-c:v", "libx264", "-pix_fmt", "yuv420p", "viz.mp4"
+])