isothermal and adiabatic wall BCs for reacting flows

examples/2D_Thermal_Flatplate/case.py

@@ -15,7 +15,7 @@
     "run_time_info": "T",
     "x_domain%beg": 0.0,
     "x_domain%end": Lx,
-    "y_domain%beg": 0.0,
+    "y_domain%beg": 0,
     "y_domain%end": Ly,
     "m": 699,
     "n": 699,

examples/2D_ibm_Thermal_Flatplate/case.py

@@ -0,0 +1,84 @@
+#!/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": -Ly / 10,
+    "y_domain%end": Ly,
+    "m": 399,
+    "n": 399,
+    "p": 0,
+    "dt": 4e-08,
+    "t_step_start": 0,
+    "t_step_stop": 150000,
+    "t_step_save": 30000,
+    "num_patches": 1,
+    "model_eqns": "5eq",
+    "alt_soundspeed": "F",
+    "num_fluids": 1,
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": "rk3",
+    "mp_weno": "F",
+    "weno_order": 5,
+    "weno_eps": 1e-16,
+    "riemann_solver": "hllc",
+    "wave_speeds": "direct",
+    "avg_state": "arithmetic",
+    "bc_x%beg": -7,
+    "bc_x%end": -3,
+    "bc_y%beg": -3,
+    "bc_y%end": -3,
+    "format": "silo",
+    "precision": "double",
+    "prim_vars_wrt": "T",
+    "parallel_io": "T",
+    "chemistry": "T",
+    "chem_params%diffusion": "T",
+    "chem_params%reactions": "F",
+    "cantera_file": ctfile,
+    "chem_wrt_T": "T",
+    "ib": "T",
+    "num_ibs": 1,
+    "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": 10 * 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,
+    "patch_ib(1)%geometry": 3,
+    "patch_ib(1)%x_centroid": Lx / 2,
+    "patch_ib(1)%y_centroid": -Ly / 20,
+    "patch_ib(1)%length_x": Lx,
+    "patch_ib(1)%length_y": Ly / 10,
+    "patch_ib(1)%slip": "F",
+    "patch_ib(1)%isothermal": "T",
+    "patch_ib(1)%Twall": 600.0,
+    "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))

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,
+    "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": "5eq",
+    "alt_soundspeed": "F",
+    "num_fluids": 1,
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": "rk3",
+    "mp_weno": "F",
+    "weno_order": 5,
+    "weno_eps": 1e-16,
+    "riemann_solver": "hllc",
+    "wave_speeds": "direct",
+    "avg_state": "arithmetic",
+    "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": "silo",
+    "precision": "double",
+    "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/2dHardcodedIC.fpp

@@ -380,7 +380,11 @@
         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
+        if (y_cc(j) >= 0.0_wp) then
+            q_prim_vf(eqn_idx%mom%beg)%sf(i, j, 0) = u_mean
+        else
+            q_prim_vf(eqn_idx%mom%beg)%sf(i, j, 0) = 0.0_wp
+        end if
         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

src/common/m_derived_types.fpp

@@ -354,6 +354,8 @@ module m_derived_types
         real(wp) :: length_x, length_y, length_z  !< Dimensions of the patch. x,y,z Lengths.
         real(wp) :: radius  !< Dimensions of the patch. radius.
         logical :: slip
+        logical :: isothermal
+        real(wp) :: Twall
         integer :: moving_ibm  !< 0 for no moving, 1 for moving, 2 for moving on forced path
         real(wp) :: mass, moment  !< mass and moment of inertia of object used to compute forces in 2-way coupling
         real(wp), dimension(1:3) :: force, torque  !< vectors for the computed force and torque values applied to an IB

src/common/m_variables_conversion.fpp

@@ -45,6 +45,16 @@ module m_variables_conversion
     $:GPU_DECLARE(create='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps]')
 #endif
 
+#ifndef MFC_SIMULATION
+    type(integer_field), public :: ghost_points_index
+    type(scalar_field), public  :: pressure_ghost_point
+    $:GPU_DECLARE(create='[ghost_points_index, pressure_ghost_point]')
+#else
+    type(integer_field), public :: ghost_points_index
+    type(scalar_field), public  :: pressure_ghost_point
+    $:GPU_DECLARE(create='[ghost_points_index, pressure_ghost_point]')
+#endif
+
     real(wp), allocatable, dimension(:)   :: Gs_vc
     integer, allocatable, dimension(:)    :: bubrs_vc
     real(wp), allocatable, dimension(:,:) :: Res_vc
@@ -361,6 +371,20 @@ contains
 
             $:GPU_UPDATE(device='[Res_vc, Re_idx, Re_size]')
         end if
+
+        if (p > 0) then
+            @:ALLOCATE(ghost_points_index%sf(0:m, 0:n, 0:p))
+            @:ALLOCATE(pressure_ghost_point%sf(0:m, 0:n, 0:p))
+        else
+            @:ALLOCATE(ghost_points_index%sf(0:m, 0:n, 0:0))
+            @:ALLOCATE(pressure_ghost_point%sf(0:m, 0:n, 0:0))
+        end if
+
+        ghost_points_index%sf = 0
+        pressure_ghost_point%sf = 0.0_wp
+
+        @:ACC_SETUP_SFs(ghost_points_index)
+        @:ACC_SETUP_SFs(pressure_ghost_point)
 #endif
 
         if (bubbles_euler) then
@@ -468,12 +492,13 @@ contains
 
     !> Convert conserved variables (rho*alpha, rho*u, E, alpha) to primitives (rho, u, p, alpha). Conversion depends on model_eqns:
     !! each model has different variable sets and EOS.
-    subroutine s_convert_conservative_to_primitive_variables(qK_cons_vf, q_T_sf, qK_prim_vf, ibounds)
+    subroutine s_convert_conservative_to_primitive_variables(qK_cons_vf, q_T_sf, qK_prim_vf, ibounds, t_step, stage)
 
         type(scalar_field), dimension(sys_size), intent(in)    :: qK_cons_vf
         type(scalar_field), intent(inout)                      :: q_T_sf
         type(scalar_field), dimension(sys_size), intent(inout) :: qK_prim_vf
         type(int_bounds_info), dimension(1:3), intent(in)      :: ibounds
+        integer, optional, intent(in)                          :: t_step, stage
 
         #:if USING_AMD and not MFC_CASE_OPTIMIZATION
             real(wp), dimension(3) :: alpha_K, alpha_rho_K
@@ -675,6 +700,14 @@ contains
                     call s_compute_pressure(qK_cons_vf(eqn_idx%E)%sf(j, k, l), qK_cons_vf(eqn_idx%alf)%sf(j, k, l), dyn_pres_K, &
                                             & pi_inf_K, gamma_K, rho_K, qv_K, rhoYks, pres, T, pres_mag=pres_mag)
 
+#ifdef MFC_SIMULATION
+                    if (.not. (t_step == 0 .and. stage == 1)) then
+                        if (ghost_points_index%sf(j, k, l) == 1) then
+                            pres = pressure_ghost_point%sf(j, k, l)
+                        end if
+                    end if
+#endif
+
                     qK_prim_vf(eqn_idx%E)%sf(j, k, l) = pres
 
                     if (chemistry) then
@@ -1234,6 +1267,8 @@ contains
 
 #ifdef MFC_SIMULATION
         @:DEALLOCATE(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps, Gs_vc)
+        @:DEALLOCATE(ghost_points_index%sf)
+        @:DEALLOCATE(pressure_ghost_point%sf)
         if (bubbles_euler) then
             @:DEALLOCATE(bubrs_vc)
         end if

src/pre_process/m_check_ib_patches.fpp

@@ -38,6 +38,7 @@ contains
             if (i <= num_ibs) then
                 call s_int_to_str(i, iStr)
                 @:PROHIBIT(patch_ib(i)%geometry == dflt_int, "IB patch undefined. patch_ib("//trim(iStr)//")%geometry must be set.")
+                call s_check_ib_patch_thermal(i)
 
                 ! Constraints on the geometric initial condition patch parameters
                 if (patch_ib(i)%geometry == 2) then
@@ -72,6 +73,19 @@ contains
 
     end subroutine s_check_ib_patches
 
+    impure subroutine s_check_ib_patch_thermal(patch_id)
+
+        integer, intent(in) :: patch_id
+
+        call s_int_to_str(patch_id, iStr)
+
+        if (.not. patch_ib(patch_id)%isothermal) return
+
+        @:PROHIBIT(f_is_default(patch_ib(patch_id)%Twall) .or. patch_ib(patch_id)%Twall <= 0._wp, &
+                   & 'in immersed boundary patch ' // trim(iStr) // ': patch_ib('//trim(iStr)//')%Twall must be set and > 0')
+
+    end subroutine s_check_ib_patch_thermal
+
     !> Verify that the geometric parameters of the circle patch have been consistently inputted.
 
     impure subroutine s_check_circle_ib_patch_geometry(patch_id)
@@ -232,7 +246,8 @@ contains
         @:PROHIBIT((.not. f_is_default(patch_ib(patch_id)%x_centroid)) .or. (.not. f_is_default(patch_ib(patch_id)%y_centroid)) &
                    & .or. (.not. f_is_default(patch_ib(patch_id)%z_centroid)) &
                    & .or. (.not. f_is_default(patch_ib(patch_id)%length_x)) .or. (.not. f_is_default(patch_ib(patch_id)%length_y)) &
-                   & .or. (.not. f_is_default(patch_ib(patch_id)%length_z)) .or. (.not. f_is_default(patch_ib(patch_id)%radius)), &
+                   & .or. (.not. f_is_default(patch_ib(patch_id)%length_z)) .or. (.not. f_is_default(patch_ib(patch_id)%radius)) &
+                   & .or. patch_ib(patch_id)%isothermal .or. (.not. f_is_default(patch_ib(patch_id)%Twall)), &
                    & 'in inactive IB patch ' // trim(iStr))
 
     end subroutine s_check_inactive_ib_patch_geometry

src/pre_process/m_global_parameters.fpp

@@ -305,6 +305,8 @@ contains
             patch_ib(i)%airfoil_id = 0
             patch_ib(i)%model_id = 0
             patch_ib(i)%slip = .false.
+            patch_ib(i)%isothermal = .false.
+            patch_ib(i)%Twall = dflt_real
 
             ! Variables to handle moving immersed boundaries, defaulting to no movement
             patch_ib(i)%moving_ibm = 0

src/pre_process/m_mpi_proxy.fpp

@@ -125,14 +125,14 @@ contains
             call MPI_BCAST(patch_ib(i)%geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
 
             #:for VAR in [ 'x_centroid', 'y_centroid', 'z_centroid', &
-                & 'length_x', 'length_y', 'length_z', 'radius']
+                & 'length_x', 'length_y', 'length_z', 'radius', 'Twall']
                 call MPI_BCAST(patch_ib(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
             #:endfor
             call MPI_BCAST(patch_ib(i)%airfoil_id, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
             call MPI_BCAST(patch_ib(i)%model_id, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
             call MPI_BCAST(patch_ib(i)%slip, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(patch_ib(i)%isothermal, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         end do
-
         ! manual: ib_airfoil (kept manual alongside patch_ib)
         do i = 1, num_ib_airfoils_max
             #:for VAR in ['c', 'p', 't', 'm']

src/simulation/m_global_parameters.fpp

@@ -581,6 +581,8 @@ contains
             patch_ib(i)%airfoil_id = 0
             patch_ib(i)%model_id = 0
             patch_ib(i)%slip = .false.
+            patch_ib(i)%isothermal = .false.
+            patch_ib(i)%Twall = dflt_real
 
             ! Variables to handle moving immersed boundaries, defaulting to no movement
             patch_ib(i)%moving_ibm = 0