Unity Lewis Implementation

docs/documentation/case.md

@@ -993,6 +993,22 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 
 - `bc_y%beg = -2` to enable reflective boundary conditions
 
+### 17. Chemistry
+
+| Parameter                     | Type    | Description                                              |
+| ---:                          | :---:   | :---                                                     |
+| `chemistry`                   | Logical | Enable chemistry simulation                              |
+| `chem_params%diffusion`       | Logical | Enable multispecies diffusion                            |
+| `chem_params%reactions`       | Logical | Enable chemical reactions                                |
+| `chem_params%gamma_method`    | Integer | Methodology for calculating the heat capacity ratio      |
+| `chem_params%transport_model` | Integer | Methodology for calculating the diffusion coefficients   |
+| `cantera_file`                | String  | Cantera-format mechanism file (e.g., .yaml)              |
+
+- `chem_params%transport_model` specifies the methodology for calculating diffusion coefficients and other transport properties,`1` for mixture-average, `2` for Unity-Lewis
+
+- `cantera_file` specifies the chemical mechanism file. If the file is part of the standard Cantera library, only the filename is required. Otherwise, the file must be located in the same directory as your `case.py` file
+
+
 ## Enumerations
 
 ### Boundary conditions

examples/1D_multispecies_diffusion/case.py

@@ -0,0 +1,78 @@
+#!/usr/bin/env python3
+# References:
+# + https://doi.org/10.1016/j.compfluid.2013.10.014: 4.4. Multicomponent diffusion test case
+
+import json
+import argparse
+import math
+import cantera as ct
+
+ctfile = "gri30.yaml"
+sol_L = ct.Solution(ctfile)
+sol_L.TPX = 300, 8000, "O2:2,N2:2,H2O:5"
+
+L = 0.05
+Nx = 100
+dx = L / Nx
+dt = 0.3e-6
+Tend = 0.05
+
+NT = int(Tend / dt)
+SAVE_COUNT = 2000
+NS = 2000
+case = {
+    "run_time_info": "T",
+    "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",
+    "model_eqns": 2,
+    "num_fluids": 1,
+    "num_patches": 1,
+    "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": 2,
+    "wave_speeds": 2,
+    "avg_state": 1,
+    "bc_x%beg": -1,
+    "bc_x%end": -1,
+    "viscous": "F",
+    "chemistry": "T",
+    "chem_params%diffusion": "T",
+    "chem_params%reactions": "F",
+    "chem_params%transport_model": 2,  # Unity-Lewis
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "chem_wrt_T": "T",
+    "patch_icpp(1)%geometry": 1,
+    "patch_icpp(1)%hcid": 182,
+    "patch_icpp(1)%x_centroid": L / 2,
+    "patch_icpp(1)%length_x": L,
+    "patch_icpp(1)%vel(1)": "0",
+    "patch_icpp(1)%pres": 1.01325e5,
+    "patch_icpp(1)%alpha(1)": 1,
+    "patch_icpp(1)%alpha_rho(1)": 1,
+    "fluid_pp(1)%gamma": 1.0e00 / (1.9326e00 - 1.0e00),
+    "fluid_pp(1)%pi_inf": 0,
+    "cantera_file": ctfile,
+}
+
+for i in range(len(sol_L.Y)):
+    case[f"chem_wrt_Y({i + 1})"] = "T"
+    case[f"patch_icpp(1)%Y({i+1})"] = 0.0
+if __name__ == "__main__":
+    print(json.dumps(case))

src/common/m_derived_types.fpp

@@ -456,6 +456,7 @@ module m_derived_types
         !> gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
         !> gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
         integer :: gamma_method
+        integer :: transport_model
     end type chemistry_parameters
 
     !> Lagrangian bubble parameters

src/post_process/m_global_parameters.fpp

@@ -303,6 +303,7 @@ module m_global_parameters
     real(wp) :: rhoref, pref
     !> @}
 
+    type(chemistry_parameters) :: chem_params
     !> @name Bubble modeling variables and parameters
     !> @{
     integer :: nb
@@ -420,6 +421,9 @@ contains
             #:endfor
         #:endfor
 
+        chem_params%gamma_method = 1
+        chem_params%transport_model = 1
+
         ! Fluids physical parameters
         do i = 1, num_fluids_max
             fluid_pp(i)%gamma = dflt_real

src/pre_process/m_global_parameters.fpp

@@ -220,6 +220,7 @@ module m_global_parameters
 
     real(wp) :: rhoref, pref !< Reference parameters for Tait EOS
 
+    type(chemistry_parameters) :: chem_params
     !> @name Bubble modeling
     !> @{
     integer :: nb
@@ -580,6 +581,9 @@ contains
             patch_ib(i)%rotation_matrix_inverse = patch_ib(i)%rotation_matrix
         end do
 
+        chem_params%gamma_method = 1
+        chem_params%transport_model = 1
+
         ! Fluids physical parameters
         do i = 1, num_fluids_max
             fluid_pp(i)%gamma = dflt_real

src/simulation/m_global_parameters.fpp

@@ -649,6 +649,7 @@ contains
         chem_params%diffusion = .false.
         chem_params%reactions = .false.
         chem_params%gamma_method = 1
+        chem_params%transport_model = 1
 
         num_bc_patches = 0
         bc_io = .false.

src/simulation/m_mpi_proxy.fpp

@@ -125,7 +125,7 @@ contains
                 call MPI_BCAST(chem_params%${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
             #:endfor
 
-            #:for VAR in [ 'gamma_method' ]
+            #:for VAR in [ 'gamma_method', 'transport_model' ]
                 call MPI_BCAST(chem_params%${VAR}$, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
             #:endfor
         end if

toolchain/mfc/run/case_dicts.py

@@ -355,7 +355,7 @@ def analytic(self):
 for var in [ 'diffusion', 'reactions' ]:
     SIMULATION[f'chem_params%{var}'] = ParamType.LOG
 
-for var in [ 'gamma_method' ]:
+for var in [ 'gamma_method', 'transport_model']:
     SIMULATION[f'chem_params%{var}'] = ParamType.INT
 
 for var in ["R0ref", "p0ref", "rho0ref", "T0ref", "ss", "pv", "vd",

toolchain/mfc/test/cases.py

@@ -1030,7 +1030,7 @@ def foreach_example():
                            "2D_backward_facing_step",
                            "2D_forward_facing_step",
                            "1D_convergence",
-                           "3D_IGR_33jet"]
+                           "3D_IGR_33jet","1D_multispecies_diffusion"]
             if path in casesToSkip:
                 continue
             name = f"{path.split('_')[0]} -> Example -> {'_'.join(path.split('_')[1:])}"