isothermal and adiabatic wall BCs for reacting flows

docs/documentation/case.md

@@ -358,6 +358,8 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi
 | `airfoil_id`         | Integer | Index into `ib_airfoil` array for NACA airfoil geometry patches. |
 | `model_id`           | Integer | Index into `stl_models` array for STL/OBJ geometry patches. |
 | `slip`               | Logical | Apply a slip boundary |
+| `isothermal`         | Logical | Impose an isothermal wall (fixed temperature) at the patch surface (reacting flows). |
+| `twall`              | Real    | Wall temperature [K] for an isothermal immersed boundary. |
 | `moving_ibm`         | Integer | Sets the method used for IB movement. |
 | `vel(i)`             | Real    | Initial velocity of the moving IB in the i-th direction. |
 | `angular_vel(i)`     | Real    | Initial angular velocity of the moving IB in the i-th direction. |
@@ -396,6 +398,8 @@ Additional details on this specification can be found in [NACA airfoil](https://
 
 - `slip` applies a slip boundary to the surface of the patch if true and a no-slip boundary condition to the surface if false.
 
+- `isothermal` imposes an isothermal wall at the immersed surface for reacting flows: the ghost-point temperature is set so that the surface is held at the fixed wall temperature `twall` (in Kelvin). When false (the default), the wall is adiabatic (zero normal temperature gradient). Isothermal immersed boundaries model heat conduction and therefore require ``chemistry = 'T'`` and ``chem_params%%diffusion = 'T'``, and `twall` must be strictly positive. To keep the ghost-point image-point interpolation well-posed, avoid placing the patch flush with a domain boundary — extend it slightly past the domain so its faces lie outside the computational domain.
+
 - For STL/OBJ geometry (geometry 5 or 12), set `model_id` to index into the `stl_models` array and specify `model_filepath`, `model_scale`, `model_translate`, and `model_threshold` on that entry.
 
 - `moving_ibm` sets the method by which movement will be applied to the immersed boundary. Using 0 will result in no movement. Using 1 will result 1-way coupling where the boundary moves at a constant rate and applied forces to the fluid based upon its own motion. In 1-way coupling, the fluid does not apply forces back onto the IB. Using 2 will result in 2-way coupling, where the boundary pushes on the fluid and the fluid pushes back on the boundary via pressure and viscous forces. If external forces are applied, the boundary will also experience those forces.

examples/2D_ibm_Thermal_Flatplate/case.py

@@ -0,0 +1,89 @@
+#!/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",
+    "fd_order": 2,
+    "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,
+    # Overhang the domain so the IB's vertical faces sit outside [0, Lx]. Inside the
+    # domain the only IB surface is the horizontal top (normal +y), so ghost-point
+    # image points reflect straight into fluid — avoids ill-conditioned corner
+    # stencils at the domain edge that diverge between CPU and GPU.
+    "patch_ib(1)%length_x": 1.2 * 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))

src/common/include/2dHardcodedIC.fpp

@@ -542,7 +542,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

@@ -355,6 +355,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,14 @@ module m_variables_conversion
     $:GPU_DECLARE(create='[gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps]')
 #endif
 
+    ! Ghost-point bookkeeping for immersed boundaries with chemistry: a flag
+    ! marking cells that are IB ghost points and the pressure imposed there,
+    ! so the conservative-to-primitive conversion can re-impose it (see
+    ! s_convert_conservative_to_primitive_variables and m_ibm).
+    type(integer_field), public :: ghost_points_index
+    type(scalar_field), public  :: pressure_ghost_point
+    $:GPU_DECLARE(create='[ghost_points_index, pressure_ghost_point]')
+
     real(wp), allocatable, dimension(:)   :: Gs_vc
     integer, allocatable, dimension(:)    :: bubrs_vc
     real(wp), allocatable, dimension(:,:) :: Res_vc
@@ -361,6 +369,30 @@ contains
 
             $:GPU_UPDATE(device='[Res_vc, Re_idx, Re_size]')
         end if
+
+        ! These fields back the IB ghost-point pressure re-imposition, which is only
+        ! read (and written) when chemistry is enabled - skip them for non-reacting IB.
+        if (ib .and. chemistry) then
+            if (p > 0) then
+                @:ALLOCATE(ghost_points_index%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))
+                @:ALLOCATE(pressure_ghost_point%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))
+            else if (n > 0) then
+                @:ALLOCATE(ghost_points_index%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))
+                @:ALLOCATE(pressure_ghost_point%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))
+            else  ! 1D
+                @:ALLOCATE(ghost_points_index%sf(-buff_size:m+buff_size, 0:0, 0:0))
+                @:ALLOCATE(pressure_ghost_point%sf(-buff_size:m+buff_size, 0:0, 0:0))
+            end if
+            ghost_points_index%sf = 0
+            pressure_ghost_point%sf = 0._wp
+
+            @:ACC_SETUP_SFs(ghost_points_index)
+            @:ACC_SETUP_SFs(pressure_ghost_point)
+            ! @:ALLOCATE only creates (uninitialized) device memory, so push the
+            ! host zeros to the device; otherwise the ghost-point flag is garbage
+            ! on the device and triggers spurious pressure overrides below.
+            $:GPU_UPDATE(device='[ghost_points_index%sf, pressure_ghost_point%sf]')
+        end if
 #endif
 
         if (bubbles_euler) then
@@ -468,12 +500,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
@@ -501,6 +534,21 @@ contains
         real(wp)               :: E, D                     !< Prim/Cons variables within Newton-Raphson iteration
         real(wp)               :: f, dGa_dW, dp_dW, df_dW  !< Functions within Newton-Raphson iteration
         integer                :: iter                     !< Newton-Raphson iteration counter
+        logical                :: apply_ib_ghost_pres      !< Re-impose stored IB ghost-point pressure (chemistry IBM)
+
+        ! For chemistry immersed boundaries the ghost-point energy encodes the wall
+        ! temperature, so the pressure recovered from the EOS is not the interpolated
+        ! image-point pressure. Re-impose the stored ghost-point pressure, except on
+        ! the very first stage of the first step when it has not been populated yet.
+        apply_ib_ghost_pres = .false.
+#ifdef MFC_SIMULATION
+        if (chemistry .and. ib .and. n > 0) then
+            apply_ib_ghost_pres = .true.
+            if (present(t_step) .and. present(stage)) then
+                if (t_step == 0 .and. stage == 1) apply_ib_ghost_pres = .false.
+            end if
+        end if
+#endif
 
         $:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, &
                             & rhoYks, B, pres, vftmp, nbub_sc, G_K, T, pres_mag, Ga, B2, m2, S, W, dW, E, D, f, dGa_dW, dp_dW, &
@@ -675,6 +723,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 (apply_ib_ghost_pres) 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 +1290,10 @@ contains
 
 #ifdef MFC_SIMULATION
         @:DEALLOCATE(gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps, Gs_vc)
+        if (ib .and. chemistry) then
+            @:DEALLOCATE(ghost_points_index%sf)
+            @:DEALLOCATE(pressure_ghost_point%sf)
+        end if
         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

@@ -344,6 +344,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_checker.fpp

@@ -53,6 +53,9 @@ contains
         @:PROHIBIT(ib_state_wrt .and. .not. ib, "ib_state_wrt requires ib to be enabled")
         @:PROHIBIT(many_ib_patch_parallelism .and. .not. ib, "many_ib_patch_parallelism requires ib to be enabled")
 
+        @:PROHIBIT(chemistry .and. ib .and. (bubbles_euler .or. qbmm), &
+                   & "immersed boundaries with chemistry are not supported together with bubbles_euler/qbmm: the ghost-point image interpolation uses a bubble/QBMM call path that does not carry the chemistry state")
+
         @:PROHIBIT(bf_spatial_support .and. (n == 0 .or. p /= 0), &
                    & "bf_spatial_support is implemented for 2D only (it forces mom%beg and mom%beg+1)")
 

src/simulation/m_global_parameters.fpp

@@ -665,6 +665,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