Isothermal Boundary Conditions for Reacting Flows

CMakeLists.txt

@@ -1026,4 +1026,4 @@ site_name(SITE_NAME)
 
 configure_file(
     "${CMAKE_CURRENT_SOURCE_DIR}/toolchain/cmake/configuration.cmake.in"
-    "${CMAKE_CURRENT_BINARY_DIR}/configuration.txt")
+    "${CMAKE_CURRENT_BINARY_DIR}/configuration.txt")

docs/documentation/case.md

@@ -1068,8 +1068,20 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 
 - `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
 
+### 18. Chemistry-Specific Boundary Conditions
 
-### 18. GPU Performance (NVIDIA UVM)
+| Parameter          | Type    | Description                                                                 |
+| ---:               | :----:  | :---                                                                        |
+| `bc_[x,y,z]%%isothermal_in`    | Logical | Enable isothermal wall at the domain entrance (minimum coordinate).         |
+| `bc_[x,y,z]%%isothermal_out`   | Logical | Enable isothermal wall at the domain exit (maximum coordinate).             |
+| `bc_[x,y,z]%%Twall_in`         | Real    | Temperature [K] of the entrance isothermal wall.                            |
+| `bc_[x,y,z]%%Twall_out`        | Real    | Temperature [K] of the exit isothermal wall.                                |
+
+This boundary condition can be used for fixed-temperature (isothermal) walls at the domain extremities. It is exclusively available for reacting flows and requires chemistry to be enabled. It properly evaluates heat and species fluxes at the interface when ``chemistry = 'T'``, ``chem_params%%diffusion = 'T'``, and the corresponding domain boundary is set to a slip wall (`bc_[x,y,z]%%[beg,end]` = -15) or a no-slip wall (`bc_[x,y,z]%%[beg,end]` = -16).
+
+
+
+### 19. GPU Performance (NVIDIA UVM)
 
 | Parameter                  | Type    | Description                                              |
 | ---:                       | :---:   | :---                                                     |

examples/2D_Thermal_Flatplate/case.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+import json
+import math
+
+import cantera as ct
+
+Lx = 0.05
+Ly = 0.05
+
+ctfile = "h2o2.yaml"
+sol_L = ct.Solution(ctfile)
+sol_L.TPY = 1125, ct.one_atm, "O2:0.21,N2:0.79"
+# Configuring case dictionary
+case = {
+    "run_time_info": "T",
+    "x_domain%beg": 0.0,
+    "x_domain%end": Lx,
+    "y_domain%beg": 0.0,
+    "y_domain%end": Ly,
+    "m": 699,
+    "n": 699,
+    "p": 0,
+    "dt": 4.0e-08,
+    "t_step_start": 0,
+    "t_step_stop": 75000,
+    "t_step_save": 4500,
+    "num_patches": 1,
+    "model_eqns": 2,
+    "alt_soundspeed": "F",
+    "num_fluids": 1,
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": 3,
+    "mp_weno": "F",
+    "weno_order": 5,
+    "weno_eps": 1e-16,
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "bc_x%beg": -7,
+    "bc_x%end": -3,
+    "bc_y%beg": -16,
+    "bc_y%end": -3,
+    "bc_y%isothermal_in": "T",
+    "bc_y%Twall_in": 600.0,
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "T",
+    "chemistry": "T",
+    "chem_params%diffusion": "T",
+    "chem_params%reactions": "F",
+    "cantera_file": ctfile,
+    "chem_wrt_T": "T",
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%hcid": 291,
+    "patch_icpp(1)%x_centroid": Lx / 2,
+    "patch_icpp(1)%y_centroid": Ly / 2,
+    "patch_icpp(1)%length_x": Lx,
+    "patch_icpp(1)%length_y": Ly,
+    "patch_icpp(1)%vel(1)": 0,
+    "patch_icpp(1)%vel(2)": 0,
+    "patch_icpp(1)%pres": 101325,
+    "patch_icpp(1)%alpha_rho(1)": 1.00,
+    "patch_icpp(1)%alpha(1)": 1,
+    "fluid_pp(1)%gamma": 1.0e00 / (1.4e00 - 1.0e00),
+    "fluid_pp(1)%pi_inf": 0.0e00,
+    "viscous": "T",
+    "fluid_pp(1)%Re(1)": 100000,
+}
+for i in range(len(sol_L.Y)):
+    case[f"chem_wrt_Y({i + 1})"] = "T"
+    case[f"patch_icpp(1)%Y({i + 1})"] = sol_L.Y[i]
+
+if __name__ == "__main__":
+    print(json.dumps(case))

src/common/include/1dHardcodedIC.fpp

@@ -51,6 +51,21 @@
         molar_mass_inv = y1/31.998_wp + y2/18.01508_wp + y3/16.04256_wp + y4/28.0134_wp
 
         q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 1.01325_wp*(10.0_wp)**5/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)
+
+    case(191)  ! 1D Dual Isothermal case
+
+        q_prim_vf(eqn_idx%E)%sf(i, 0, 0) = 101325.0_wp
+        q_prim_vf(eqn_idx%mom%beg)%sf(i, 0, 0) = 0.0_wp
+        q_prim_vf(eqn_idx%species%beg)%sf(i, 0, 0) = 1.0_wp
+
+        if (x_cc(i) <= 0.025_wp) then
+            temp = 700.0_wp + ((1000.0_wp - 700.0_wp)/0.025_wp)*x_cc(i)
+        else
+            temp = 1200.0_wp + ((900.0_wp - 1000.0_wp)/0.025_wp)*(x_cc(i) - 0.025_wp)
+        end if
+
+        molar_mass_inv = 1.0_wp/2.01588_wp
+        q_prim_vf(eqn_idx%cont%beg)%sf(i, 0, 0) = 101325.0_wp/(temp*8.3144626_wp*1000.0_wp*molar_mass_inv)
     case default
         call s_int_to_str(patch_id, iStr)
         call s_mpi_abort("Invalid hcid specified for patch " // trim(iStr))

src/common/include/2dHardcodedIC.fpp

@@ -8,8 +8,16 @@
     real(wp) :: sinA, cosA
     real(wp) :: r_sq
 
+    ! # 291 - Shear/Thermal Layer Case
+    real(wp) :: delta_shear, u_max, u_mean
+    real(wp) :: T_wall, T_inf, P_atm, T_loc
+    real(wp) :: delta_th, R_mix
+    real(wp) :: Y_N2, Y_O2, MW_N2, MW_O2
+    real(wp) :: bottom_blend_u, bottom_blend_T
+
     ! # 207
     real(wp) :: sigma, gauss1, gauss2
+
     ! # 208
     real(wp) :: ei, d, fsm, alpha_air, alpha_sf6
 
@@ -305,6 +313,32 @@
             q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, &
                       & 0) = 112.99092883944267*((0.1/0.3))*x_cc(i)*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2)))
         end if
+    case (291)  ! Isothermal Flat Plate
+        T_inf = 1125.0_wp
+        T_wall = 600.0_wp
+        P_atm = 101325.0_wp
+
+        ! Boundary/Shear Layer thicknesses
+        delta_th = 0.0003_wp  ! Thermal BL thickness
+        delta_shear = 8e-3_wp  ! Velocity BL thickness
+
+        u_max = 50.0_wp  ! Freestream Velocity (m/s)
+
+        MW_N2 = 28.0134e-3_wp
+        MW_O2 = 31.999e-3_wp
+        Y_N2 = 0.767_wp
+        Y_O2 = 0.233_wp
+        R_mix = 8.314462618_wp*((Y_N2/MW_N2) + (Y_O2/MW_O2))
+        bottom_blend_u = tanh(y_cc(j)/delta_shear)
+        bottom_blend_T = tanh(y_cc(j)/delta_th)
+        u_mean = u_max*bottom_blend_u
+        T_loc = T_wall + (T_inf - T_wall)*bottom_blend_T
+        q_prim_vf(eqn_idx%cont%beg)%sf(i, j, 0) = P_atm/(R_mix*T_loc)
+        q_prim_vf(eqn_idx%mom%beg)%sf(i, j, 0) = u_mean
+        q_prim_vf(eqn_idx%mom%end)%sf(i, j, 0) = 0.0_wp
+        q_prim_vf(eqn_idx%E)%sf(i, j, 0) = P_atm
+        q_prim_vf(eqn_idx%species%beg)%sf(i, j, 0) = Y_O2
+        q_prim_vf(eqn_idx%species%end)%sf(i, j, 0) = Y_N2
     case default
         if (proc_rank == 0) then
             call s_int_to_str(patch_id, iStr)