Add heterogeneous reacting surface boundary conditions

Author: Thomas Jackson

docs/documentation/case.md

@@ -356,6 +356,9 @@ 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 |
+| `thermal_bc`         | Integer | Thermal boundary-condition selector: 0 = zero-normal-gradient temperature, 1 = prescribed wall temperature, 2 = reacting surface energy balance. |
+| `Twall`              | Real    | Prescribed wall temperature used when `thermal_bc = 1`. |
+| `surface_reaction`   | Integer | Heterogeneous surface-reaction flag: 0 = disabled, 1 = enabled. |
 | `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. |
@@ -394,6 +397,12 @@ 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.
 
+- `thermal_bc` selects the thermal immersed-boundary condition. A value of 0 applies a zero-normal-gradient temperature condition, 1 prescribes the wall temperature using `Twall`, and 2 solves the reacting-surface energy balance for the surface temperature. The `thermal_bc = 2` option requires `surface_reaction = 1`.
+
+- `Twall` specifies the prescribed surface temperature when `thermal_bc = 1` and must be positive in that case.
+
+- `surface_reaction` enables heterogeneous surface chemistry when set to 1. Surface reactions require `chemistry = T` and cannot be combined with `inj_species > 0`.
+
 - 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.
@@ -1184,13 +1193,17 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 | `chem_params%%adap_substeps`   | Logical | Per-rank adaptive sub-step count driven by local stiffness  |
 | `chem_params%%reaction_substeps_max` | Integer | Sub-step ceiling when `adap_substeps` is enabled       |
 | `cantera_file`                | String  | Cantera-format mechanism file (e.g., .yaml)              |
+| `surface_cantera_file`        | String  | Cantera-format mechanism file for heterogeneous surface chemistry |
+| `surface_phase`               | String  | Cantera interface phase name for heterogeneous surface chemistry |
 
 - `chem_params%%transport_model` specifies the methodology for calculating diffusion coefficients and other transport properties, `1` for mixture-average, `2` for Unity-Lewis
 - `chem_params%%reaction_substeps` controls how the reaction source is integrated. With `0` (default) the net production rates are added to the flow right-hand side and advanced by the flow time stepper (fine for hydrogen). With a value `> 0`, the reaction is instead integrated by operator splitting after each flow update: every cell's constant-density, constant-internal-energy reactor is advanced over the timestep with that many sub-steps of an **α-QSS** (quasi-steady-state) integrator — a matrix-free, Jacobian-free predictor–corrector (Mott/CHEMEQ2) that splits the net rate into creation/destruction parts and applies a Padé α-weighting, so it stays stable on stiff mechanisms where an explicit source diverges. This decouples the (often much faster) chemical timescale from the flow timestep and is required for stiff mechanisms — e.g. hydrocarbons such as GRI-Mech methane, which otherwise diverge on the first step
 - `chem_params%%adap_substeps` (default `F`) makes each rank choose its α-QSS sub-step count per flow step from the largest chemical stiffness among its own cells: the count sits at `reaction_substeps` (the floor) in inert or burned gas and rises toward `reaction_substeps_max` (the ceiling) only across the reaction front. It uses no MPI collectives. When enabled, `reaction_substeps >= 1` and `reaction_substeps_max >= reaction_substeps` are required
 
 - `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
 
+- `surface_cantera_file` and `surface_phase` specify the Cantera mechanism file and interface phase used for heterogeneous surface chemistry. These parameters must be specified together when a surface mechanism is used.
+
 ### 18. Chemistry-Specific Boundary Conditions
 
 | Parameter          | Type    | Description                                                                 |

src/common/m_derived_types.fpp

@@ -349,11 +349,21 @@ module m_derived_types
         real(wp), dimension(1:3,1:3) :: rotation_matrix
         !> matrix that converts from fluid reference frame to IB reference frame
         real(wp), dimension(1:3,1:3) :: rotation_matrix_inverse
-        integer :: airfoil_id  !< index into ib_airfoil(:) for airfoil geometry patches
-        integer :: model_id  !< index into stl_models(:) for STL/OBJ geometry patches
-        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
+        integer                      :: airfoil_id                    !< index into ib_airfoil(:) for airfoil geometry patches
+        integer                      :: model_id                      !< index into stl_models(:) for STL/OBJ geometry patches
+        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
+
+        ! Thermal immersed-boundary condition
+        !    0 = zero-normal-gradient temperature
+        !    1 = prescribed wall temperature (Twall)
+        !    2 = reacting surface energy balance
+        integer  :: thermal_bc
+        real(wp) :: Twall
+
+        ! Heterogeneous surface reaction    0 = none    1 = enabled
+        integer :: surface_reaction
         integer :: moving_ibm  !< 0 for no moving, 1 for moving, 2 for moving on forced path
         real(wp) :: v_blow  !< Wall-normal surface blowing speed (burning/transpiring IB surface); 0 = impermeable
         integer :: inj_species  !< Injected species index at a blowing surface (chemistry); 0 = mirror ambient

src/pre_process/m_global_parameters.fpp

@@ -339,6 +339,11 @@ contains
             patch_ib(i)%airfoil_id = 0
             patch_ib(i)%model_id = 0
             patch_ib(i)%slip = .false.
+
+            patch_ib(i)%thermal_bc = 0
+            patch_ib(i)%Twall = 0._wp
+            patch_ib(i)%surface_reaction = 0
+
             patch_ib(i)%v_blow = 0._wp
             patch_ib(i)%inj_species = 0
             patch_ib(i)%burn_rate_exp = 0._wp

src/pre_process/m_mpi_proxy.fpp

@@ -125,13 +125,15 @@ 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', 'v_blow', &
+                & 'length_x', 'length_y', 'length_z', 'radius', 'Twall', 'v_blow', &
                 & 'burn_rate_exp', 'burn_rate_pref']
                 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)%inj_species, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(patch_ib(i)%thermal_bc, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(patch_ib(i)%surface_reaction, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
             call MPI_BCAST(patch_ib(i)%slip, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         end do
 

src/simulation/m_checker.fpp

@@ -34,7 +34,7 @@ contains
             end if
         end if
 
-        if (ib .and. chemistry) then
+        if (ib) then
             call s_check_inputs_ib_injection
         end if
 
@@ -101,16 +101,36 @@ contains
 
     end subroutine s_check_inputs_nvidia_uvm
 
-    !> Validates that each burning immersed-boundary patch injects a species index within the mechanism. inj_species indexes the
-    !! image-point mass-fraction array Ys_IP(1:num_species) in m_ibm; an out-of-range value is an out-of-bounds write (silent
-    !! corruption). Only reachable with chemistry.
+    !> Validates immersed-boundary injection, thermal, and heterogeneous surface-reaction parameters.
     impure subroutine s_check_inputs_ib_injection
 
         integer :: i
 
         do i = 1, num_ibs
-            @:PROHIBIT(patch_ib(i)%inj_species > num_species, &
-                       & "patch_ib inj_species must be <= num_species (it indexes the image-point species mass fractions; an out-of-range value writes out of bounds)")
+            ! Basic parameter ranges
+            @:PROHIBIT(patch_ib(i)%inj_species < 0 .or. patch_ib(i)%inj_species > num_species, &
+                       & "patch_ib inj_species must be in [0,num_species]")
+            @:PROHIBIT(patch_ib(i)%thermal_bc < 0 .or. patch_ib(i)%thermal_bc > 2, "patch_ib thermal_bc must be 0, 1, or 2")
+            @:PROHIBIT(patch_ib(i)%surface_reaction < 0 .or. patch_ib(i)%surface_reaction > 1, &
+                       & "patch_ib surface_reaction must be 0 or 1")
+
+            ! Thermal immersed-boundary condition
+            !    0 = zero-normal-gradient temperature
+            !    1 = prescribed wall temperature (Twall)
+            !    2 = reacting surface energy balance
+            if (patch_ib(i)%thermal_bc == 1) then
+                @:PROHIBIT(patch_ib(i)%Twall <= 0._wp, "patch_ib Twall must be > 0 when thermal_bc = 1")
+            end if
+
+            if (patch_ib(i)%thermal_bc == 2) then
+                @:PROHIBIT(patch_ib(i)%surface_reaction /= 1, "patch_ib thermal_bc = 2 requires surface_reaction = 1")
+            end if
+
+            ! Heterogeneous surface reaction    0 = none    1 = enabled
+            if (patch_ib(i)%surface_reaction == 1) then
+                @:PROHIBIT(.not. chemistry, "patch_ib surface_reaction = 1 requires chemistry = T")
+                @:PROHIBIT(patch_ib(i)%inj_species > 0, "patch_ib surface_reaction = 1 cannot be combined with inj_species > 0")
+            end if
         end do
 
     end subroutine s_check_inputs_ib_injection

src/simulation/m_global_parameters.fpp

@@ -660,6 +660,11 @@ contains
             patch_ib(i)%airfoil_id = 0
             patch_ib(i)%model_id = 0
             patch_ib(i)%slip = .false.
+
+            patch_ib(i)%thermal_bc = 0
+            patch_ib(i)%Twall = 0._wp
+            patch_ib(i)%surface_reaction = 0
+
             patch_ib(i)%v_blow = 0._wp
             patch_ib(i)%inj_species = 0
             patch_ib(i)%burn_rate_exp = 0._wp