Add patch_ib%thermal_layer: start a hot IB in its conduction solution

docs/documentation/case.md

@@ -365,6 +365,7 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi
 | `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`. |
+| `thermal_layer`      | Real    | Thickness of the initial conductive layer around a `thermal_bc = 1` circle or sphere [m]; 0 = none. |
 | `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. |
@@ -420,6 +421,8 @@ Additional details on this specification can be found in [NACA airfoil](https://
 
 - `Twall` specifies the prescribed surface temperature when `thermal_bc = 1` and must be positive in that case.
 
+- `thermal_layer` = δ > 0 starts the gas around a `thermal_bc = 1` circle (geometry 2) or sphere (geometry 8) in the conduction solution for a wall brought to `Twall` a time t₀ = δ²/(4α) earlier: T = T∞ + (`Twall` − T∞)(R/r)^((d−1)/2) erfc((r − R)/δ), with pressure and composition unchanged. It is exact for a sphere and the leading-order term for a cylinder. Without it a hot wall meets ambient gas as a temperature step at t = 0, and the near-wall cells heat faster than they can expand; on coarse grids that transient can drive the state non-physical. Resolve δ with a few cells. Where layers overlap, the IB with the largest temperature change sets the cell, and periodic directions use the nearest image.
+
 - `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.

src/common/m_derived_types.fpp

@@ -366,6 +366,7 @@ module m_derived_types
         !    2 = reacting surface energy balance
         integer  :: thermal_bc
         real(wp) :: Twall
+        real(wp) :: thermal_layer  !< Initial conductive-layer thickness around a thermal_bc = 1 IB [m]; 0 = none
 
         ! Heterogeneous surface reaction    0 = none    1 = enabled
         integer :: surface_reaction

src/pre_process/m_global_parameters.fpp

@@ -343,6 +343,7 @@ contains
 
             patch_ib(i)%thermal_bc = 0
             patch_ib(i)%Twall = 0._wp
+            patch_ib(i)%thermal_layer = 0._wp
             patch_ib(i)%surface_reaction = 0
 
             patch_ib(i)%v_blow = 0._wp

src/pre_process/m_initial_condition.fpp

@@ -126,6 +126,9 @@ contains
         if (mixlayer_perturb) call s_perturb_mixlayer(ic%q_prim_vf)
         if (simplex_perturb) call s_perturb_simplex(ic%q_prim_vf)
         if (chemistry) call s_compute_T_from_primitives(ic%q_T_sf, ic%q_prim_vf, idwint)
+        if (chemistry .and. ib) then
+            if (any(patch_ib(1:num_ibs)%thermal_layer > 0._wp)) call s_apply_ib_thermal_layers(ic%q_prim_vf, ic%q_T_sf)
+        end if
 
         if (elliptic_smoothing .and. chemistry) then
             call s_elliptic_smoothing(ic%q_prim_vf, ic%bc_type, ic%q_T_sf)
@@ -143,6 +146,56 @@ contains
 
     end subroutine s_generate_initial_condition
 
+    !> Start the gas around each prescribed-temperature circle/sphere IB in the conduction solution for a wall impulsively brought
+    !! to Twall a time t0 ago, T = T_inf + (Twall - T_inf)*(R/r)^((d-1)/2)*erfc((r - R)/delta), delta = 2*sqrt(alpha*t0). Exact for
+    !! a sphere and the leading-order term for a cylinder (Carslaw & Jaeger). Pressure and composition are kept, so density scales
+    !! as T_inf/T. Where layers overlap, the IB with the largest temperature change sets the cell.
+    impure subroutine s_apply_ib_thermal_layers(q_prim_vf, q_T_sf)
+
+        type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf
+        type(scalar_field), intent(inout)                      :: q_T_sf
+        ! erfc(6) ~ 2e-17: beyond 6 delta a layer changes nothing
+        real(wp), parameter :: reach = 6._wp
+        real(wp)            :: x(3), c(3), span(3), r, R_ib, dT, dT_i, T_inf
+        logical             :: periodic(3)
+        integer             :: i, j, k, l, d
+
+        span = [x_domain%end - x_domain%beg, y_domain%end - y_domain%beg, z_domain%end - z_domain%beg]
+        periodic = [bc_x%beg == BC_PERIODIC, bc_y%beg == BC_PERIODIC, bc_z%beg == BC_PERIODIC]
+
+        do l = 0, p
+            do k = 0, n
+                do j = 0, m
+                    x = 0._wp
+                    x(1) = x_cc(j)
+                    if (n > 0) x(2) = y_cc(k)
+                    if (p > 0) x(3) = z_cc(l)
+                    T_inf = q_T_sf%sf(j, k, l)
+                    dT = 0._wp
+                    do i = 1, num_ibs
+                        if (.not. patch_ib(i)%thermal_layer > 0._wp) cycle
+                        c = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, patch_ib(i)%z_centroid]
+                        do d = 1, num_dims
+                            c(d) = x(d) - c(d)
+                            if (periodic(d)) c(d) = c(d) - span(d)*nint(c(d)/span(d))
+                        end do
+                        r = norm2(c(1:num_dims))
+                        R_ib = patch_ib(i)%radius
+                        if (r - R_ib > reach*patch_ib(i)%thermal_layer) cycle
+                        dT_i = (patch_ib(i)%Twall - T_inf)*(R_ib/max(r, R_ib))**(0.5_wp*real(num_dims - 1, &
+                                & wp))*erfc(max(r - R_ib, 0._wp)/patch_ib(i)%thermal_layer)
+                        if (abs(dT_i) > abs(dT)) dT = dT_i
+                    end do
+                    if (dT /= 0._wp) then
+                        q_prim_vf(eqn_idx%cont%beg)%sf(j, k, l) = q_prim_vf(eqn_idx%cont%beg)%sf(j, k, l)*T_inf/(T_inf + dT)
+                        q_T_sf%sf(j, k, l) = T_inf + dT
+                    end if
+                end do
+            end do
+        end do
+
+    end subroutine s_apply_ib_thermal_layers
+
     !> Deallocation procedures for the module
     impure subroutine s_finalize_initial_condition_module
 

src/pre_process/m_mpi_proxy.fpp

@@ -125,7 +125,7 @@ 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', 'Twall', 'v_blow', &
+                & 'length_x', 'length_y', 'length_z', 'radius', 'Twall', 'thermal_layer', 'v_blow', &
                 & 'burn_rate_exp', 'burn_rate_pref']
                 call MPI_BCAST(patch_ib(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
             #:endfor

src/simulation/m_global_parameters.fpp

@@ -695,6 +695,7 @@ contains
 
             patch_ib(i)%thermal_bc = 0
             patch_ib(i)%Twall = 0._wp
+            patch_ib(i)%thermal_layer = 0._wp
             patch_ib(i)%surface_reaction = 0
 
             patch_ib(i)%v_blow = 0._wp

src/simulation/m_start_up.fpp

@@ -1204,6 +1204,7 @@ contains
         ! reaction and blowing fields must be set here rather than left as whatever patch_ib held.
         ib_patch%thermal_bc = 0
         ib_patch%Twall = 0._wp
+        ib_patch%thermal_layer = 0._wp
         ib_patch%surface_reaction = 0
         ib_patch%v_blow = 0._wp
         ib_patch%inj_species = 0