Add thermal_bc = 3: a lumped solid whose wall temperature follows its heat balance (after #1956)

docs/documentation/case.md

@@ -363,8 +363,11 @@ 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`. |
+| `thermal_bc`         | Integer | Thermal boundary-condition selector: 0 = zero-normal-gradient temperature, 1 = prescribed wall temperature, 2 = reacting surface energy balance, 3 = lumped solid whose temperature evolves. |
+| `Twall`              | Real    | Wall temperature: prescribed (`thermal_bc = 1`) or initial (`thermal_bc = 3`). |
+| `rho_solid`, `cp_solid` | Real | Solid density [kg/m³] and heat capacity [J/kg/K] (`thermal_bc = 3`). |
+| `emissivity`, `T_rad` | Real   | Surface emissivity and radiative surroundings temperature [K] (`thermal_bc = 3`; default 0, no radiation). |
+| `heat_power`         | Real    | Heat supplied to the body [W; W/m per unit depth in 2D], e.g. Joule heating (`thermal_bc = 3`). |
 | `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. |
@@ -416,7 +419,7 @@ 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`. A non-zero `thermal_bc` requires `chemistry = T` and cannot be combined with `inj_species > 0`, since the thermal condition is applied by the chemistry ghost-state reconstruction, which an injecting surface bypasses.
+- `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`, 2 solves the reacting-surface energy balance for the surface temperature, and 3 treats the body as one lumped solid whose temperature, starting at `Twall`, evolves as m c_s dT/dt = Q_in + c_s (T − 298.15 K) Ṁ_out + `heat_power` − ε σ A (T⁴ − `T_rad`⁴). Q_in and Ṁ_out are the energy into and mass out of the body summed over every face between its cells and fluid cells, from the same total fluxes the flow update applies, weighted by the Runge–Kutta stages; over a step the body gains exactly the energy the fluid loses, reaction heat and the enthalpy of the gasified carbon included. m = `rho_solid`·V, with V from the geometry (circle per unit depth, sphere, or cylinder), and A is the surface area from the `ib_surface_wrt` surface points. One temperature per body is valid while the Biot number h·R/k_solid is small (graphite particles and mm rods). The temperature is clamped to the thermodynamic window [200, 5000] K and is carried across restarts in `restart_data/ib_state`. The `thermal_bc = 2` option requires `surface_reaction = 1`. A non-zero `thermal_bc` requires `chemistry = T` and cannot be combined with `inj_species > 0`, since the thermal condition is applied by the chemistry ghost-state reconstruction, which an injecting surface bypasses.
 
 - `Twall` specifies the prescribed surface temperature when `thermal_bc = 1` and must be positive in that case.
 
@@ -773,6 +776,7 @@ To restart the simulation from $k$-th time step, see @ref running "Restarting Ca
 | `heat_ratio_wrt`        | Logical | Add the specific heat ratio to the database	|
 | `ib_force_wrt`          | Logical | Record the immersed-boundary force history to `D/ib_forces.dat` (default off) |
 | `ib_force_stride`       | Integer | Stride, in time steps, of the per-step immersed-boundary force record (default 1) |
+| `ib_surface_wrt`        | Logical | Write the wall temperature and gasified mass flux at each thermal/reacting IB surface point at every save (default off) |
 | `ib_state_wrt`          | Logical | Parameter to handle writing IB state on saves and outputting the state as a point mesh to SILO files. |
 | `pi_inf_wrt`            | Logical | Add the liquid stiffness function to the database |
 | `pres_inf_wrt`          | Logical | Add the liquid stiffness to the formatted database	 |
@@ -844,6 +848,8 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 
 - `ib_force_wrt` records the force, torque and kinematics of every immersed boundary in a single shared text file, `D/ib_forces.dat`, described below. It is off by default: the history is written every step, which at large rank counts is a cost a run should opt into rather than inherit. `ib_force_stride` writes only every N-th step, for runs long enough that the history itself becomes large.
 
+- `ib_surface_wrt` writes, at every save, one text file per rank, `D/ib_surface_<rank>_<save>.dat`, with a line per surface point of each chemistry IB (`thermal_bc` /= 0 or `surface_reaction` = 1) and no header, columns `x y z ib nx ny nz area T_wall mdot`. `ib` is the global IB index, `(nx, ny, nz)` the level-set normal, `T_wall` the surface temperature [K] and `mdot` the gasified mass flux [kg/m²/s] from the surface solve (0 on an inert surface). `area` [m², or m per unit depth in 2D] is the surface each point stands for: summing `area`·`mdot` over an IB's lines gives its mass loss rate [kg/s], and `mdot`/ρ_solid is the local surface regression rate. The points are the ghost points within two cell sizes h = (cell volume)^(1/d) of the surface, each weighted by cell volume/(2h); for circles, spheres and cylinder sides the weight is also scaled by (R/(R − depth))^(d−1), the area ratio between the surface and the layer the point sits in. The areas then sum to the surface area without bias (a sphere at R = 12h: within 0.3%; a circle at R = 20h: within about 1%); other shapes keep the uncorrected band, which underestimates a convex surface by about (d−1)·h/R.
+
 #### Immersed-boundary force history {#sec-ib-force-history}
 
 `D/ib_forces.dat` holds one fixed-width record per body per written step. Its twenty columns are

examples/2D_ibm_poiseuille_nn/compare_analytic.py

@@ -95,7 +95,7 @@ def main():
     t_last = int(os.path.basename(dats[-1]).split("_")[1].split(".")[0])
     ib_state = os.path.join(RESTART_DIR, f"ib_state_{t_last}.dat")
     if os.path.exists(ib_state):
-        rec = np.fromfile(ib_state, dtype=np.float64).reshape(2, 20)
+        rec = np.fromfile(ib_state, dtype=np.float64).reshape(2, -1)
         f_ana = RHO * G_X * 0.5 * (Y_HI - Y_LO) * L_X  # tau_w*L_x, nominal H
         print(f"IBM x-force per wall: {rec[0, 1]:.4e} (bottom), {rec[1, 1]:.4e} (top); " f"analytic tau_w*L_x = {f_ana:.4e}")
         print(f"  force ratios vs analytic: {rec[0, 1] / f_ana:.3f}, {rec[1, 1] / f_ana:.3f}")

src/common/m_constants.fpp

@@ -21,6 +21,10 @@ module m_constants
     !! wall temperature outside it.
     real(wp), parameter :: T_surface_min = 200._wp
     real(wp), parameter :: T_surface_max = 5000._wp
+
+    !> Reals per IB in restart_data/ib_state: time, force(3), torque(3), vel(3), angular_vel(3), angles(3), centroid(3), radius,
+    !! Twall. Twall carries a thermal_bc = 3 body's evolving temperature across a restart.
+    integer, parameter :: ib_state_nfields = 21
     !> Radius cutoff to avoid division by zero for 3D spherical harmonic patch (geometry 14)
     real(wp), parameter :: small_radius = 1.e-32_wp
     integer, parameter  :: num_stcls_min = 5        !< Minimum # of stencils

src/common/m_derived_types.fpp

@@ -363,8 +363,12 @@ module m_derived_types
         !    0 = zero-normal-gradient temperature
         !    1 = prescribed wall temperature (Twall)
         !    2 = reacting surface energy balance
+        !    3 = lumped solid: Twall evolves with the body's heat balance
         integer  :: thermal_bc
         real(wp) :: Twall
+        real(wp) :: rho_solid, cp_solid  !< Solid density [kg/m^3] and heat capacity [J/kg/K] (thermal_bc = 3)
+        real(wp) :: emissivity, T_rad    !< Surface emissivity and radiative surroundings temperature [K] (thermal_bc = 3)
+        real(wp) :: heat_power           !< Heat supplied to the body, e.g. Joule heating [W; W/m per unit depth in 2D]
 
         ! Heterogeneous surface reaction    0 = none    1 = enabled
         integer :: surface_reaction

src/common/m_helper.fpp

@@ -10,17 +10,17 @@ module m_helper
 
     use m_derived_types
     use m_global_parameters
-    use m_constants, only: BC_PERIODIC
+    use m_constants, only: BC_PERIODIC, ib_state_nfields
     use ieee_arithmetic  !< For checking NaN
 
     implicit none
 
     private
-    public :: s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, s_simpson, s_transcoeff, &
-        & s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, f_create_transform_matrix, &
-        & f_create_bbox, s_print_2D_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, double_factorial, factorial, &
-        & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul, &
-        & s_prng_splitmix32, f_local_rank_owns_location
+    public :: s_pack_ib_state, s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, &
+        & s_simpson, s_transcoeff, s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, &
+        & f_create_transform_matrix, f_create_bbox, s_print_2D_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, &
+        & double_factorial, factorial, f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, &
+        & s_prng, modmul, s_prng_splitmix32, f_local_rank_owns_location
 
 contains
 
@@ -771,4 +771,25 @@ contains
 
     end function f_local_rank_owns_location
 
+    !> One IB's restart_data/ib_state record; see ib_state_nfields for the layout.
+    pure subroutine s_pack_ib_state(ib_patch, time, buf)
+
+        type(ib_patch_parameters), intent(in)              :: ib_patch
+        real(wp), intent(in)                               :: time
+        real(wp), dimension(ib_state_nfields), intent(out) :: buf
+
+        buf(1) = time
+        buf(2:4) = ib_patch%force(1:3)
+        buf(5:7) = ib_patch%torque(1:3)
+        buf(8:10) = ib_patch%vel(1:3)
+        buf(11:13) = ib_patch%angular_vel(1:3)
+        buf(14:16) = ib_patch%angles(1:3)
+        buf(17) = ib_patch%x_centroid
+        buf(18) = ib_patch%y_centroid
+        buf(19) = ib_patch%z_centroid
+        buf(20) = ib_patch%radius
+        buf(21) = ib_patch%Twall
+
+    end subroutine s_pack_ib_state
+
 end module m_helper

src/post_process/m_data_output.fpp

@@ -13,7 +13,8 @@ module m_data_output
     use m_helper
     use m_variables_conversion
     use m_eos
-    use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, format_silo, format_binary, precision_single
+    use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, format_silo, format_binary, precision_single, &
+        & ib_state_nfields
 
     implicit none
 
@@ -1363,8 +1364,7 @@ contains
         character(len=len_trim(case_dir) + 3*name_len) :: file_loc
 
 #ifdef MFC_MPI
-        integer, parameter                              :: NFIELDS_PER_IB = 20
-        real(wp)                                        :: ib_buf(NFIELDS_PER_IB)
+        real(wp)                                        :: ib_buf(ib_state_nfields)
         real(wp), dimension(:,:), allocatable           :: ib_data
         logical                                         :: file_exist
         character(LEN=4*name_len), dimension(num_procs) :: meshnames
@@ -1389,7 +1389,7 @@ contains
             nBodies = num_ibs
 
             if (nBodies > 0) then
-                allocate (ib_data(nBodies, NFIELDS_PER_IB))
+                allocate (ib_data(nBodies, ib_state_nfields))
                 allocate (px(nBodies), py(nBodies), pz(nBodies))
                 allocate (force_x(nBodies), force_y(nBodies), force_z(nBodies))
                 allocate (torque_x(nBodies), torque_y(nBodies), torque_z(nBodies))

src/pre_process/m_data_output.fpp

@@ -22,7 +22,7 @@ module m_data_output
     use m_boundary_io
     use m_thermochem, only: species_names
     use m_helper
-    use m_constants, only: model_eqns_5eq, precision_single
+    use m_constants, only: model_eqns_5eq, precision_single, ib_state_nfields
 
     implicit none
 
@@ -786,16 +786,10 @@ contains
 
             type(ib_patch_parameters), intent(in) :: ib_patch
             integer, intent(in)                   :: gbl_id
-            real(wp), dimension(20)               :: ib_buf
-
-            ib_buf = 0._wp
-            ib_buf(8:10) = ib_patch%vel
-            ib_buf(11:13) = ib_patch%angular_vel
-            ib_buf(14:16) = ib_patch%angles
-            ib_buf(17) = ib_patch%x_centroid
-            ib_buf(18) = ib_patch%y_centroid
-            ib_buf(19) = ib_patch%z_centroid
-            ib_buf(20) = ib_patch%radius
+            real(wp), dimension(ib_state_nfields) :: ib_buf
+
+            call s_pack_ib_state(ib_patch, 0._wp, ib_buf)
+            ib_buf(2:7) = 0._wp  ! no force or torque before the first step
 
             if (file_per_process) write (file_unit) gbl_id
             write (file_unit) ib_buf

src/pre_process/m_global_parameters.fpp

@@ -343,6 +343,11 @@ contains
 
             patch_ib(i)%thermal_bc = 0
             patch_ib(i)%Twall = 0._wp
+            patch_ib(i)%rho_solid = 0._wp
+            patch_ib(i)%cp_solid = 0._wp
+            patch_ib(i)%emissivity = 0._wp
+            patch_ib(i)%T_rad = 0._wp
+            patch_ib(i)%heat_power = 0._wp
             patch_ib(i)%surface_reaction = 0
 
             patch_ib(i)%v_blow = 0._wp

src/simulation/m_data_output.fpp

@@ -19,7 +19,7 @@ module m_data_output
     use m_delay_file_access
     use m_ibm
     use m_boundary_common
-    use m_constants, only: model_eqns_5eq, precision_single
+    use m_constants, only: model_eqns_5eq, precision_single, ib_state_nfields
 
     implicit none
 
@@ -1142,8 +1142,7 @@ contains
         integer, dimension(MPI_STATUS_SIZE)  :: status
         logical                              :: file_exist, dir_check
         integer                              :: i, ib_idx
-        integer, parameter                   :: NFIELDS_PER_IB = 20
-        real(wp)                             :: ib_buf(NFIELDS_PER_IB)
+        real(wp)                             :: ib_buf(ib_state_nfields)
         integer                              :: file_unit
         character(len=10)                    :: t_step_string
 
@@ -1175,16 +1174,7 @@ contains
             write (file_unit) num_local_ibs
             do i = 1, num_local_ibs
                 ib_idx = local_ib_patch_ids(i)
-                ib_buf(1) = mytime
-                ib_buf(2:4) = patch_ib(ib_idx)%force(1:3)
-                ib_buf(5:7) = patch_ib(ib_idx)%torque(1:3)
-                ib_buf(8:10) = patch_ib(ib_idx)%vel(1:3)
-                ib_buf(11:13) = patch_ib(ib_idx)%angular_vel(1:3)
-                ib_buf(14:16) = patch_ib(ib_idx)%angles(1:3)
-                ib_buf(17) = patch_ib(ib_idx)%x_centroid
-                ib_buf(18) = patch_ib(ib_idx)%y_centroid
-                ib_buf(19) = patch_ib(ib_idx)%z_centroid
-                ib_buf(20) = patch_ib(ib_idx)%radius
+                call s_pack_ib_state(patch_ib(ib_idx), mytime, ib_buf)
 
                 write (file_unit) patch_ib(ib_idx)%gbl_patch_id
                 write (file_unit) ib_buf
@@ -1211,21 +1201,12 @@ contains
 
             do i = 1, num_local_ibs
                 ib_idx = local_ib_patch_ids(i)
-                ib_buf(1) = mytime
-                ib_buf(2:4) = patch_ib(ib_idx)%force(1:3)
-                ib_buf(5:7) = patch_ib(ib_idx)%torque(1:3)
-                ib_buf(8:10) = patch_ib(ib_idx)%vel(1:3)
-                ib_buf(11:13) = patch_ib(ib_idx)%angular_vel(1:3)
-                ib_buf(14:16) = patch_ib(ib_idx)%angles(1:3)
-                ib_buf(17) = patch_ib(ib_idx)%x_centroid
-                ib_buf(18) = patch_ib(ib_idx)%y_centroid
-                ib_buf(19) = patch_ib(ib_idx)%z_centroid
-                ib_buf(20) = patch_ib(ib_idx)%radius
+                call s_pack_ib_state(patch_ib(ib_idx), mytime, ib_buf)
 
                 ! Global IB index determines position in file
-                disp = int(patch_ib(ib_idx)%gbl_patch_id - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK
+                disp = int(patch_ib(ib_idx)%gbl_patch_id - 1, MPI_OFFSET_KIND)*int(ib_state_nfields, MPI_OFFSET_KIND)*WP_MOK
 
-                call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, NFIELDS_PER_IB, mpi_p, status, ierr)
+                call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, ib_state_nfields, mpi_p, status, ierr)
             end do
 
             call MPI_FILE_CLOSE(ifile, ierr)
@@ -1240,8 +1221,7 @@ contains
         integer, intent(in)                  :: t_step
         character(LEN=path_len + 2*name_len) :: file_loc
         integer                              :: i, ios, file_unit
-        integer, parameter                   :: NFIELDS_PER_IB = 20
-        real(wp)                             :: ib_buf(NFIELDS_PER_IB)
+        real(wp)                             :: ib_buf(ib_state_nfields)
 
         call s_create_directory(trim(case_dir) // '/restart_data')
 
@@ -1252,16 +1232,7 @@ contains
         if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc))
 
         do i = 1, num_ibs
-            ib_buf(1) = mytime
-            ib_buf(2:4) = patch_ib(i)%force(1:3)
-            ib_buf(5:7) = patch_ib(i)%torque(1:3)
-            ib_buf(8:10) = patch_ib(i)%vel(1:3)
-            ib_buf(11:13) = patch_ib(i)%angular_vel(1:3)
-            ib_buf(14:16) = patch_ib(i)%angles(1:3)
-            ib_buf(17) = patch_ib(i)%x_centroid
-            ib_buf(18) = patch_ib(i)%y_centroid
-            ib_buf(19) = patch_ib(i)%z_centroid
-            ib_buf(20) = patch_ib(i)%radius
+            call s_pack_ib_state(patch_ib(i), mytime, ib_buf)
 
             write (file_unit) ib_buf
         end do

src/simulation/m_global_parameters.fpp

@@ -518,6 +518,7 @@ contains
         ib_coefficient_of_friction = dflt_real
         ib_state_wrt = .false.
         ib_force_wrt = .false.
+        ib_surface_wrt = .false.
         ib_force_stride = 1
         many_ib_patch_parallelism = .false.
 
@@ -695,6 +696,11 @@ contains
 
             patch_ib(i)%thermal_bc = 0
             patch_ib(i)%Twall = 0._wp
+            patch_ib(i)%rho_solid = 0._wp
+            patch_ib(i)%cp_solid = 0._wp
+            patch_ib(i)%emissivity = 0._wp
+            patch_ib(i)%T_rad = 0._wp
+            patch_ib(i)%heat_power = 0._wp
             patch_ib(i)%surface_reaction = 0
 
             patch_ib(i)%v_blow = 0._wp