diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 57a51d62c..bda07e4d6 100644 --- a/docs/documentation/case.md +++ b/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,6 +1193,8 @@ 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 @@ -1191,6 +1202,8 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met: - `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 | diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index d80c885b4..9d9564dd3 100644 --- a/src/common/m_derived_types.fpp +++ b/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 diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index d4a5c559b..7b01bd5c5 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/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 diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 45c1a3b50..a505505d5 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/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 diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index 0c1bc8f2d..d51f9112b 100644 --- a/src/simulation/m_checker.fpp +++ b/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 diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index e1a258766..d477b6801 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/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 diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index a853e82cc..42c626d2c 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -21,7 +21,10 @@ module m_ibm use m_model use m_patch_geometries use m_collisions - use m_thermochem, only: num_species, gas_constant, get_mixture_molecular_weight, get_mixture_energy_mass + use m_thermochem, only: num_species, gas_constant, molecular_weights, get_mixture_molecular_weight, get_mixture_energy_mass, & + & get_mixture_thermal_conductivity_mixavg, get_species_mass_diffusivities_mixavg + + use m_surface_thermochem, only: get_surface_net_production_rates, get_surface_reaction_heat_flux implicit none @@ -162,17 +165,21 @@ contains real(wp), dimension(3) :: r_IP, v_IP, pb_IP, mv_IP real(wp), dimension(18) :: nmom_IP real(wp), dimension(12) :: presb_IP, massv_IP - real(wp), dimension(10) :: Ys_IP + real(wp), dimension(10) :: Ys_IP, Ys_g #:else real(wp), dimension(num_fluids) :: Gs real(wp), dimension(num_fluids) :: alpha_rho_IP, alpha_IP real(wp), dimension(nb) :: r_IP, v_IP, pb_IP, mv_IP real(wp), dimension(nb*nmom) :: nmom_IP real(wp), dimension(nb*nnode) :: presb_IP, massv_IP - real(wp), dimension(num_species) :: Ys_IP + real(wp), dimension(num_species) :: Ys_IP, Ys_g + real(wp) :: W_species(num_species) #:endif real(wp) :: T_IP, mw_IP, e_IP !< Image-point temperature, mixture MW, and mass-specific internal energy (chemistry) real(wp) :: v_blow_eff !< Effective surface blowing speed (after any pressure-coupled burn-rate scaling) + real(wp), dimension(num_species) :: Ys_s + real(wp) :: T_s, T_g, mw_s, mw_g, rho_s, mdot_s, v_stefan, d + logical :: surface_converged ! Primitive variables at the image point associated with a ghost point, interpolated from surrounding fluid cells. real(wp), dimension(3) :: norm !< Normal vector from GP to IP @@ -223,7 +230,8 @@ contains $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, & & r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, & & innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, patch_id, & - & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff, vel_sum_g, E_ghost]') + & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff, Ys_g, Ys_s, T_s, T_g, mw_s, mw_g, rho_s, mdot_s, & + & v_stefan, d, surface_converged, vel_sum_g, E_ghost]') do i = 1, num_gps gp = ghost_points(i) j = gp%loc(1) @@ -267,6 +275,52 @@ contains alpha_rho_IP(1) = pres_IP*mw_IP/(T_IP*gas_constant) end if + ! Thermal and heterogeneous reacting-surface boundary conditions. + v_stefan = 0._wp + surface_converged = .false. + + if (chemistry .and. patch_ib(patch_id)%inj_species == 0) then + ! Intrinsic gas state at the image point: rho_IP = (alpha*rho)_IP / alpha_IP + call get_mixture_molecular_weight(Ys_IP, mw_IP) + T_IP = pres_IP*mw_IP*alpha_IP(1)/(alpha_rho_IP(1)*gas_constant) + + if (patch_ib(patch_id)%surface_reaction == 0) then + ! Inert surface: zero species flux. + Ys_g(:) = Ys_IP(:) + + ! thermal_bc = 0: zero normal temperature gradient thermal_bc = 1: prescribed surface temperature Twall + T_g = T_IP + 2._wp*real(patch_ib(patch_id)%thermal_bc, kind=wp)*(patch_ib(patch_id)%Twall - T_IP) + + call get_mixture_molecular_weight(Ys_g, mw_g) + alpha_rho_IP(1) = alpha_IP(1)*pres_IP*mw_g/(gas_constant*T_g) + else + ! Heterogeneous reacting surface. + d = abs(real(gp%levelset, kind=wp)) + + W_species(:) = molecular_weights(:) + + call s_solve_surface(pres_IP, T_IP, patch_ib(patch_id)%Twall, d, Ys_IP, W_species, & + & patch_ib(patch_id)%thermal_bc, Ys_s, T_s, mdot_s, surface_converged) + + if (surface_converged) then + call get_mixture_molecular_weight(Ys_s, mw_s) + + ! Intrinsic gas density at the reacting surface. + rho_s = pres_IP*mw_s/(gas_constant*T_s) + if (rho_s > 0._wp) v_stefan = mdot_s/rho_s + + Ys_g(:) = 2._wp*Ys_s(:) - Ys_IP(:) + T_g = 2._wp*T_s - T_IP + + call get_mixture_molecular_weight(Ys_g, mw_g) + alpha_rho_IP(1) = alpha_IP(1)*pres_IP*mw_g/(gas_constant*T_g) + else + Ys_g(:) = Ys_IP(:) + T_g = T_IP + end if + end if + end if + dyn_pres = 0._wp ! Set q_prim_vf params at GP so that mixture vars calculated properly @@ -361,6 +415,13 @@ contains if (buf > 0._wp) vel_g = vel_g + v_blow_eff*norm/buf end if + if (chemistry .and. patch_ib(patch_id)%inj_species == 0 .and. patch_ib(patch_id)%surface_reaction == 1 & + & .and. surface_converged) then + norm(1:3) = gp%levelset_norm + buf = sqrt(sum(norm**2)) + if (buf > 0._wp) vel_g = vel_g + v_stefan*norm/buf + end if + ! Set momentum vel_sum_g = 0._wp $:GPU_LOOP(parallelism='[seq]') @@ -384,18 +445,29 @@ contains ! Set Energy if (chemistry) then - ! Mirror the reacting-mixture state at the ghost point: interpolated species, - ! plus a thermodynamically consistent conserved energy from the mixture EOS. - ! (The gamma*pres_IP closure below is only valid for a calorically perfect gas - ! and yields an out-of-range temperature when inverted against the Cantera model.) - mw_IP = 0._wp - call get_mixture_molecular_weight(Ys_IP, mw_IP) - T_IP = pres_IP*mw_IP/(rho*gas_constant) - call get_mixture_energy_mass(T_IP, Ys_IP, e_IP) - $:GPU_LOOP(parallelism='[seq]') - do q = 1, num_species - q_cons_vf(eqn_idx%species%beg + q - 1)%sf(j, k, l) = rho*Ys_IP(q) - end do + ! Use the reconstructed thermal/species ghost state for an inert + ! thermal surface or a converged heterogeneous reacting surface. + if (patch_ib(patch_id)%inj_species == 0 .and. (patch_ib(patch_id)%surface_reaction == 0 & + & .or. surface_converged)) then + + call get_mixture_energy_mass(T_g, Ys_g, e_IP) + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_species + q_cons_vf(eqn_idx%species%beg + q - 1)%sf(j, k, l) = rho*Ys_g(q) + end do + else + ! Ordinary chemistry/injection, or fallback after a failed + ! heterogeneous surface solve: retain the image-point state. + mw_IP = 0._wp + call get_mixture_molecular_weight(Ys_IP, mw_IP) + T_IP = pres_IP*mw_IP/(rho*gas_constant) + call get_mixture_energy_mass(T_IP, Ys_IP, e_IP) + $:GPU_LOOP(parallelism='[seq]') + do q = 1, num_species + q_cons_vf(eqn_idx%species%beg + q - 1)%sf(j, k, l) = rho*Ys_IP(q) + end do + end if + q_cons_vf(eqn_idx%E)%sf(j, k, l) = rho*e_IP + dyn_pres else call s_compute_energy(pres_IP, alpha_rho_IP, alpha_IP, vel_sum_g, E_ghost) @@ -1603,4 +1675,301 @@ contains end subroutine s_finalize_ibm_module + !> Species flux residual at a heterogeneous reacting surface. + subroutine s_surface_species_residual(pres, T_s, d, Ys_IP, Ys_s, W_species, R_species, omega_s, mdot_s) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: pres, T_s, d + real(wp), intent(in) :: Ys_IP(num_species), Ys_s(num_species), W_species(num_species) + real(wp), intent(out) :: R_species(num_species), omega_s(num_species), mdot_s + real(wp) :: mw_IP, mw_s, rho_s, sum_BG + real(wp) :: Xs_IP(num_species), Xs_s(num_species) + real(wp) :: D_s(num_species), B_s(num_species), G_s(num_species) + integer :: k + + call get_mixture_molecular_weight(Ys_IP, mw_IP) + call get_mixture_molecular_weight(Ys_s, mw_s) + rho_s = pres*mw_s/(gas_constant*T_s) + + do k = 1, num_species + Xs_IP(k) = Ys_IP(k)*mw_IP/W_species(k) + Xs_s(k) = Ys_s(k)*mw_s/W_species(k) + end do + + call get_species_mass_diffusivities_mixavg(pres, T_s, Ys_s, D_s) + call get_surface_net_production_rates(rho_s, T_s, Ys_s, omega_s) + + mdot_s = 0._wp + do k = 1, num_species + mdot_s = mdot_s + W_species(k)*omega_s(k) + end do + + sum_BG = 0._wp + do k = 1, num_species + B_s(k) = rho_s*D_s(k)*W_species(k)/mw_s + G_s(k) = (Xs_IP(k) - Xs_s(k))/d + sum_BG = sum_BG + B_s(k)*G_s(k) + end do + + do k = 1, num_species + R_species(k) = -B_s(k)*G_s(k) + Ys_s(k)*(sum_BG + mdot_s) - W_species(k)*omega_s(k) + end do + + end subroutine s_surface_species_residual + + !> Surface energy residual: gas-side conduction balances heterogeneous reaction heat. Radiation and solid-side conduction are + !! omitted. + subroutine s_surface_energy_residual(pres, T_IP, T_s, d, Ys_s, R_energy) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: pres, T_IP, T_s, d + real(wp), intent(in) :: Ys_s(num_species) + real(wp), intent(out) :: R_energy + real(wp) :: mw_s, rho_s, k_s, q_rxn + + call get_mixture_molecular_weight(Ys_s, mw_s) + rho_s = pres*mw_s/(gas_constant*T_s) + + call get_mixture_thermal_conductivity_mixavg(T_s, Ys_s, k_s) + call get_surface_reaction_heat_flux(rho_s, T_s, Ys_s, q_rxn) + + R_energy = k_s*(T_s - T_IP)/d - q_rxn + + end subroutine s_surface_energy_residual + + !> Assemble the Newton residual for Ns species, with temperature appended only when it is solved. + subroutine s_surface_residual(pres, T_IP, T_s, d, Ys_IP, Ys_s, W_species, solve_temperature, flux_scale, energy_scale, R, & + & R_species, omega_s, mdot_s) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: pres, T_IP, T_s, d, flux_scale, energy_scale + real(wp), intent(in) :: Ys_IP(num_species), Ys_s(num_species), W_species(num_species) + logical, intent(in) :: solve_temperature + real(wp), intent(out) :: R(num_species + 1), R_species(num_species), omega_s(num_species), mdot_s + real(wp) :: R_energy + integer :: k + + call s_surface_species_residual(pres, T_s, d, Ys_IP, Ys_s, W_species, R_species, omega_s, mdot_s) + + R = 0._wp + do k = 1, num_species - 1 + R(k) = R_species(k)/flux_scale + end do + R(num_species) = sum(Ys_s) - 1._wp + + if (solve_temperature) then + call s_surface_energy_residual(pres, T_IP, T_s, d, Ys_s, R_energy) + R(num_species + 1) = R_energy/energy_scale + end if + + end subroutine s_surface_residual + + !> Newton solve for a reacting surface. thermal_bc=0: zero-normal-gradient T; 1: prescribed T; 2: energy balance. + subroutine s_solve_surface(pres, T_IP, T_wall, d, Ys_IP, W_species, thermal_bc, Ys_s, T_s, mdot_s, converged) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), intent(in) :: pres, T_IP, T_wall, d + real(wp), intent(in) :: Ys_IP(num_species), W_species(num_species) + integer, intent(in) :: thermal_bc + real(wp), intent(out) :: Ys_s(num_species), T_s, mdot_s + logical, intent(out) :: converged + integer, parameter :: max_iter = 30, max_backtrack = 20 + real(wp), parameter :: fd_eps_Y = 1.e-7_wp, fd_eps_T = 1.e-6_wp + real(wp), parameter :: tol = 1.e-8_wp, T_min = 200._wp, T_max = 5000._wp + real(wp), parameter :: Y_tol = 100._wp*epsilon(1._wp) + real(wp) :: A(num_species + 1, num_species + 1), rhs(num_species + 1), delta(num_species + 1) + real(wp) :: R(num_species + 1), R_pert(num_species + 1), R_trial(num_species + 1) + real(wp) :: R_species(num_species), R_species_pert(num_species), R_species_trial(num_species) + real(wp) :: omega_s(num_species), omega_pert(num_species), omega_trial(num_species) + real(wp) :: Ys_pert(num_species), Ys_trial(num_species) + real(wp) :: mdot_pert, mdot_trial, T_pert, T_trial + real(wp) :: flux_scale, energy_scale, R_energy, dx, lambda, norm_R, norm_trial + logical :: solve_temperature, linear_success, accepted + integer :: nsolve, iter, j, iback + + converged = .false. + Ys_s = Ys_IP + + select case (thermal_bc) + case (0) + T_s = T_IP + solve_temperature = .false. + case (1) + T_s = T_wall + solve_temperature = .false. + case (2) + T_s = min(max(T_IP, T_min), T_max) + solve_temperature = .true. + case default + T_s = T_IP + omega_s = 0._wp + mdot_s = 0._wp + return + end select + + nsolve = num_species + merge(1, 0, solve_temperature) + + call s_surface_species_residual(pres, T_s, d, Ys_IP, Ys_s, W_species, R_species, omega_s, mdot_s) + flux_scale = max(maxval(abs(R_species)), 1.e-12_wp) + energy_scale = 1._wp + if (solve_temperature) then + call s_surface_energy_residual(pres, T_IP, T_s, d, Ys_s, R_energy) + energy_scale = max(abs(R_energy), 1._wp) + end if + + call s_surface_residual(pres, T_IP, T_s, d, Ys_IP, Ys_s, W_species, solve_temperature, flux_scale, energy_scale, R, & + & R_species, omega_s, mdot_s) + norm_R = maxval(abs(R(1:nsolve))) + if (norm_R < tol) then + converged = .true. + return + end if + + do iter = 1, max_iter + do j = 1, num_species + Ys_pert = Ys_s + T_pert = T_s + dx = fd_eps_Y*max(abs(Ys_s(j)), 1._wp) + if (Ys_s(j) + dx > 1._wp) dx = -dx + Ys_pert(j) = Ys_pert(j) + dx + + call s_surface_residual(pres, T_IP, T_pert, d, Ys_IP, Ys_pert, W_species, solve_temperature, flux_scale, & + & energy_scale, R_pert, R_species_pert, omega_pert, mdot_pert) + A(1:nsolve,j) = (R_pert(1:nsolve) - R(1:nsolve))/dx + end do + + if (solve_temperature) then + Ys_pert = Ys_s + dx = fd_eps_T*max(abs(T_s), 1._wp) + T_pert = T_s + dx + if (T_pert > T_max) then + dx = -dx + T_pert = T_s + dx + end if + + call s_surface_residual(pres, T_IP, T_pert, d, Ys_IP, Ys_pert, W_species, solve_temperature, flux_scale, & + & energy_scale, R_pert, R_species_pert, omega_pert, mdot_pert) + A(1:nsolve,nsolve) = (R_pert(1:nsolve) - R(1:nsolve))/dx + end if + + rhs(1:nsolve) = -R(1:nsolve) + call s_solve_surface_linear_system(A, rhs, delta, nsolve, linear_success) + if (.not. linear_success) return + + lambda = 1._wp + accepted = .false. + do iback = 1, max_backtrack + Ys_trial = Ys_s + lambda*delta(1:num_species) + T_trial = T_s + if (solve_temperature) T_trial = T_s + lambda*delta(nsolve) + + if (minval(Ys_trial) < -Y_tol .or. maxval(Ys_trial) > 1._wp + Y_tol .or. T_trial < T_min .or. T_trial > T_max) then + lambda = 0.5_wp*lambda + cycle + end if + + where (Ys_trial < 0._wp) Ys_trial = 0._wp + where (Ys_trial > 1._wp) Ys_trial = 1._wp + + call s_surface_residual(pres, T_IP, T_trial, d, Ys_IP, Ys_trial, W_species, solve_temperature, flux_scale, & + & energy_scale, R_trial, R_species_trial, omega_trial, mdot_trial) + norm_trial = maxval(abs(R_trial(1:nsolve))) + if (norm_trial < norm_R) then + accepted = .true. + exit + end if + lambda = 0.5_wp*lambda + end do + + if (.not. accepted) return + + Ys_s = Ys_trial + T_s = T_trial + R = R_trial + R_species = R_species_trial + omega_s = omega_trial + mdot_s = mdot_trial + norm_R = norm_trial + + if (norm_R < tol) then + converged = .true. + return + end if + end do + + end subroutine s_solve_surface + + !> Small dense linear solve with partial pivoting for the local surface Newton system. + subroutine s_solve_surface_linear_system(A, b, x, nsolve, success) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: nsolve + real(wp), intent(inout) :: A(num_species + 1, num_species + 1) + real(wp), intent(inout) :: b(num_species + 1) + real(wp), intent(out) :: x(num_species + 1) + logical, intent(out) :: success + real(wp) :: factor, pivot_value, tmp, row_tmp(num_species + 1) + integer :: i, j, k, pivot + + success = .true. + x = 0._wp + + do k = 1, nsolve - 1 + pivot = k + pivot_value = abs(A(k, k)) + do i = k + 1, nsolve + if (abs(A(i, k)) > pivot_value) then + pivot = i + pivot_value = abs(A(i, k)) + end if + end do + if (pivot_value <= epsilon(1._wp)) then + success = .false. + return + end if + + if (pivot /= k) then + row_tmp(1:nsolve) = A(k,1:nsolve) + A(k,1:nsolve) = A(pivot,1:nsolve) + A(pivot,1:nsolve) = row_tmp(1:nsolve) + tmp = b(k) + b(k) = b(pivot) + b(pivot) = tmp + end if + + do i = k + 1, nsolve + factor = A(i, k)/A(k, k) + A(i, k) = 0._wp + do j = k + 1, nsolve + A(i, j) = A(i, j) - factor*A(k, j) + end do + b(i) = b(i) - factor*b(k) + end do + end do + + if (abs(A(nsolve, nsolve)) <= epsilon(1._wp)) then + success = .false. + return + end if + + x(nsolve) = b(nsolve)/A(nsolve, nsolve) + do i = nsolve - 1, 1, -1 + tmp = b(i) + do j = i + 1, nsolve + tmp = tmp - A(i, j)*x(j) + end do + if (abs(A(i, i)) <= epsilon(1._wp)) then + success = .false. + return + end if + x(i) = tmp/A(i, i) + end do + + end subroutine s_solve_surface_linear_system + end module m_ibm diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index f8fe85a9b..1040d02e6 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -168,7 +168,7 @@ contains do i = 1, num_ibs #:for VAR in [ 'radius', 'length_x', 'length_y', 'length_z', & & 'x_centroid', 'y_centroid', 'z_centroid', 'slip', 'mass', 'v_blow', & - & 'burn_rate_exp', 'burn_rate_pref'] + & 'Twall', 'burn_rate_exp', 'burn_rate_pref'] call MPI_BCAST(patch_ib(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) #:endfor #:for VAR in ['vel', 'angular_vel', 'angles'] @@ -179,6 +179,8 @@ contains 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) end do ! manual: ib_airfoil (kept manual alongside patch_ib) diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 6e0843352..699569602 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -392,8 +392,11 @@ contains ! allocated (not default-initialized) and s_reduce_ib_patch_array copies the whole ! struct into patch_ib, overwriting the defaults from ! s_assign_default_values_to_user_inputs -- so anything left unset here reaches the - ! solver as uninitialized memory (a nonzero v_blow injects a garbage wall-normal - ! velocity and NaNs the field). + ! solver as uninitialized memory. Therefore all thermal, species, and blowing + ! surface-condition fields must be initialized explicitly. + particle_cloud_ibs(ib_idx)%thermal_bc = 0 + particle_cloud_ibs(ib_idx)%Twall = 0._wp + particle_cloud_ibs(ib_idx)%surface_reaction = 0 particle_cloud_ibs(ib_idx)%v_blow = 0._wp particle_cloud_ibs(ib_idx)%inj_species = 0 particle_cloud_ibs(ib_idx)%burn_rate_exp = 0._wp diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 15f496211..8ea816d9f 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -646,6 +646,8 @@ def _load(): # Chemistry _r("cantera_file", STR, {"chemistry"}) + _r("surface_cantera_file", STR, {"chemistry"}) + _r("surface_phase", STR, {"chemistry"}) _r("chemistry", LOG, {"chemistry"}) # Condensed-phase reactive burn (programmed pressure burn on the multi-fluid model) @@ -951,9 +953,9 @@ def _load(): # grow patch_ib beyond this at runtime, but those entries are never in the namelist. _ib_tags = {"ib"} _ib_attrs: Dict[str, tuple] = {} - for a in ["geometry", "moving_ibm", "airfoil_id", "model_id", "inj_species"]: + for a in ["geometry", "moving_ibm", "airfoil_id", "model_id", "inj_species", "thermal_bc", "surface_reaction"]: _ib_attrs[a] = (INT, _ib_tags) - for a, pt in [("radius", REAL), ("slip", LOG), ("mass", REAL), ("v_blow", REAL), ("burn_rate_exp", REAL), ("burn_rate_pref", REAL)]: + for a, pt in [("radius", REAL), ("slip", LOG), ("mass", REAL), ("Twall", REAL), ("v_blow", REAL), ("burn_rate_exp", REAL), ("burn_rate_pref", REAL)]: _ib_attrs[a] = (pt, _ib_tags) for j in range(1, 4): _ib_attrs[f"angles({j})"] = (REAL, _ib_tags) diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index 91c8f453e..8ae9c15b8 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -152,6 +152,8 @@ "files_dir": "The relative path to the directory containing the extrusion IC files", "file_extension": "The last 6 digits of the extrusion files prim.XX.YY.123456.dat", "cantera_file": "Cantera mechanism file for chemistry", + "surface_cantera_file": "Cantera mechanism file for heterogeneous surface chemistry", + "surface_phase": "Cantera interface phase name for heterogeneous surface chemistry", "old_grid": "Use grid from previous simulation", "old_ic": "Use initial conditions from previous simulation", "t_step_old": "Time step to restart from", diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index ac8438833..61705ffc1 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -49,7 +49,12 @@ def _registry(): return REGISTRY -IGNORE = ["cantera_file", "chemistry"] +IGNORE = [ + "cantera_file", + "surface_cantera_file", + "surface_phase", + "chemistry", +] ALL = _ParamTypeMapping() CASE_OPTIMIZATION = [n for n, p in _registry().all_params.items() if p.case_optimization] SCHEMA = _registry().get_json_schema() diff --git a/toolchain/mfc/run/input.py b/toolchain/mfc/run/input.py index dd3a122cf..5affcda12 100644 --- a/toolchain/mfc/run/input.py +++ b/toolchain/mfc/run/input.py @@ -66,6 +66,390 @@ def get_cantera_solution(self): raise common.MFCException(f"Cantera file '{cantera_file}' not found. Searched: {', '.join(candidates)}.") + def get_cantera_surface(self): + # Lazy import to avoid slow startup for commands that don't need chemistry + import cantera as ct + import yaml + + surface_file = self.params.get("surface_cantera_file") + surface_phase = self.params.get("surface_phase") + + if surface_file is None and surface_phase is None: + return None + + if surface_file is None or surface_phase is None: + raise common.MFCException("surface_cantera_file and surface_phase must be specified together.") + + candidates = [ + surface_file, + os.path.join(self.dirpath, surface_file), + os.path.join(common.MFC_MECHANISMS_DIR, surface_file), + ] + + gas = self.get_cantera_solution() + + for candidate in candidates: + if not os.path.isfile(candidate): + continue + + try: + with open(candidate, "r", encoding="utf-8") as stream: + mechanism = yaml.safe_load(stream) + + phases = mechanism.get("phases", []) + + interface_data = None + for phase in phases: + if phase.get("name") == surface_phase: + interface_data = phase + break + + if interface_data is None: + raise common.MFCException(f"Surface phase '{surface_phase}' was not found in '{candidate}'.") + + adjacent_names = interface_data.get("adjacent-phases", []) + + adjacent = [] + + for phase_name in adjacent_names: + if phase_name == gas.name: + adjacent.append(gas) + else: + adjacent.append(ct.Solution(candidate, phase_name)) + + return ct.Interface( + candidate, + surface_phase, + adjacent=adjacent, + ) + + except common.MFCException: + raise + except Exception as e: + cons.print(f"[dim] Cantera: skipping surface mechanism " f"'{candidate}': {e}[/dim]") + continue + + raise common.MFCException(f"Cantera surface file '{surface_file}' with phase " f"'{surface_phase}' could not be loaded. " f"Searched: {', '.join(candidates)}.") + + def generate_surface_thermochem(self, sol, surface, directive_str=None) -> str: + """Generate the MFC heterogeneous surface-chemistry Fortran module.""" + + if directive_str == "mp": + gpu_routine_define = "#define GPU_ROUTINE(name) !$omp declare target device_type(any)" + elif directive_str == "acc": + gpu_routine_define = "#define GPU_ROUTINE(name) !$acc routine seq" + else: + gpu_routine_define = "#define GPU_ROUTINE(name) ! name" + + if surface is None: + num_species = len(sol.species_names) + + lines = [ + "! This file is automatically generated by the MFC toolchain.", + "! Do not edit manually.", + "", + gpu_routine_define, + "", + "module m_surface_thermochem", + "", + " use m_precision_select, only: wp", + "", + " implicit none", + "", + " private", + " public :: get_surface_net_production_rates", + " public :: get_surface_reaction_heat_flux", + "", + "contains", + "", + " subroutine get_surface_net_production_rates( &", + " density, temperature, mass_fractions, omega_s)", + "", + " GPU_ROUTINE(get_surface_net_production_rates)", + "", + " real(wp), intent(in) :: density", + " real(wp), intent(in) :: temperature", + f" real(wp), intent(in) :: mass_fractions({num_species})", + f" real(wp), intent(out) :: omega_s({num_species})", + "", + " omega_s = 0._wp", + "", + " end subroutine get_surface_net_production_rates", + "", + " subroutine get_surface_reaction_heat_flux( &", + " density, temperature, mass_fractions, q_rxn)", + "", + " GPU_ROUTINE(get_surface_reaction_heat_flux)", + "", + " real(wp), intent(in) :: density", + " real(wp), intent(in) :: temperature", + f" real(wp), intent(in) :: mass_fractions({num_species})", + " real(wp), intent(out) :: q_rxn", + "", + " q_rxn = 0._wp", + "", + " end subroutine get_surface_reaction_heat_flux", + "", + "end module m_surface_thermochem", + "", + ] + + return "\n".join(lines) + + gas_species = sol.species_names + gas_index = {name: i + 1 for i, name in enumerate(gas_species)} + + # Collect non-gas species participating in heterogeneous reactions. + # Their thermodynamic data are needed for reaction enthalpies. + nongas_species = {} + + for phase in surface.adjacent.values(): + for species in phase.species(): + if species.name not in gas_index: + nongas_species[species.name] = species + + # Verify that all species appearing in the surface reactions belong to + # either the gas mechanism, the surface phase, or an adjacent phase. + for reaction in surface.reactions(): + for name in set(reaction.reactants) | set(reaction.products): + if name in gas_index: + continue + if name in surface.species_names: + raise common.MFCException(f"Surface-site species '{name}' in heterogeneous reaction " "stoichiometry are not currently supported.") + if name in nongas_species: + continue + + raise common.MFCException(f"Surface reaction species '{name}' is not present in any " "phase associated with the surface mechanism.") + + # Return a Fortran expression for h/(R*T) for a non-gas species. + # Cantera's NasaPoly2 stores: + # [Tmid, a1_high, ..., a7_high, a1_low, ..., a7_low] + # and + # h/(R*T) = a1 + a2*T/2 + a3*T^2/3 + a4*T^3/4 + # + a5*T^4/5 + a6/T. + def nongas_h_rt_expressions(species): + thermo = species.thermo + + if thermo is None or thermo.__class__.__name__ != "NasaPoly2": + raise common.MFCException(f"Surface thermochemistry for non-gas species " f"'{species.name}' requires Cantera NasaPoly2 thermo data.") + + coeffs = list(thermo.coeffs) + if len(coeffs) != 15: + raise common.MFCException(f"Unexpected NasaPoly2 coefficient count for " f"surface species '{species.name}'.") + + tmid = coeffs[0] + high = coeffs[1:8] + low = coeffs[8:15] + + def h_rt(a): + return ( + f"({a[0]:.16e}_wp" + f" + ({a[1]:.16e}_wp)*temperature/2._wp" + f" + ({a[2]:.16e}_wp)*temperature**2/3._wp" + f" + ({a[3]:.16e}_wp)*temperature**3/4._wp" + f" + ({a[4]:.16e}_wp)*temperature**4/5._wp" + f" + ({a[5]:.16e}_wp)/temperature)" + ) + + return tmid, h_rt(low), h_rt(high) + + # Precompute generated h/(R*T) expressions for participating + # non-gas species from adjacent phases. + nongas_h_rt = {} + reaction_species = set() + for reaction in surface.reactions(): + reaction_species.update(reaction.reactants) + reaction_species.update(reaction.products) + + for name in reaction_species: + if name in gas_index: + continue + if name in nongas_species: + nongas_h_rt[name] = nongas_h_rt_expressions(nongas_species[name]) + + lines = [ + "! This file is automatically generated by the MFC toolchain.", + "! Do not edit manually.", + "", + gpu_routine_define, + "", + "module m_surface_thermochem", + "", + " use m_precision_select, only: wp", + " use m_thermochem, only: gas_constant, get_species_enthalpies_rt", + "", + " implicit none", + "", + " private", + " public :: get_surface_net_production_rates", + " public :: get_surface_reaction_heat_flux", + "", + "contains", + "", + " subroutine get_surface_net_production_rates( &", + " density, temperature, mass_fractions, omega_s)", + "", + " GPU_ROUTINE(get_surface_net_production_rates)", + "", + " real(wp), intent(in) :: density", + " real(wp), intent(in) :: temperature", + f" real(wp), intent(in) :: mass_fractions({len(gas_species)})", + f" real(wp), intent(out) :: omega_s({len(gas_species)})", + "", + f" real(wp) :: concentrations({len(gas_species)})", + " real(wp) :: rate", + "", + " omega_s = 0._wp", + " concentrations = 0._wp", + ] + + for i, species in enumerate(sol.species()): + lines.append(f" concentrations({i + 1}) = " f"density*mass_fractions({i + 1})/" f"{species.molecular_weight:.16e}_wp") + + lines.append("") + + def append_reaction_rate(lines, reaction_number, reaction): + rate = reaction.rate + + if not hasattr(rate, "pre_exponential_factor"): + raise common.MFCException(f"Surface reaction {reaction_number} does not use a " "supported Arrhenius rate expression.") + + A = rate.pre_exponential_factor + b = rate.temperature_exponent + Ea = rate.activation_energy + + lines.extend( + [ + f" ! Surface reaction {reaction_number}", + f" rate = {A:.16e}_wp", + ] + ) + + if b != 0.0: + lines.append(f" rate = rate*temperature**({b:.16e}_wp)") + + if Ea != 0.0: + lines.append(f" rate = rate*exp(-{Ea:.16e}_wp/" "(gas_constant*temperature))") + + # Cantera reaction orders override stoichiometric reactant orders. + orders = dict(reaction.orders) + for name in orders: + if name in surface.species_names: + raise common.MFCException(f"Surface reaction {reaction_number} specifies an explicit " f"reaction order for surface-site species '{name}', which is " "not currently supported.") + + for name, nu in reaction.reactants.items(): + if name not in gas_index: + continue + + order = orders.get(name, nu) + lines.append(f" rate = rate*concentrations({gas_index[name]})**" f"({order:.16e}_wp)") + + # Explicit orders may include gas species not appearing as reactants. + for name, order in orders.items(): + if name not in gas_index or name in reaction.reactants: + continue + + lines.append(f" rate = rate*concentrations({gas_index[name]})**" f"({order:.16e}_wp)") + + for reaction_number, reaction in enumerate(surface.reactions(), start=1): + append_reaction_rate(lines, reaction_number, reaction) + + for name in gas_species: + nu = reaction.products.get(name, 0.0) - reaction.reactants.get(name, 0.0) + if nu == 0.0: + continue + + lines.append(f" omega_s({gas_index[name]}) = " f"omega_s({gas_index[name]}) " f"+ ({nu:.16e}_wp)*rate") + + lines.append("") + + lines.extend( + [ + " end subroutine get_surface_net_production_rates", + "", + " subroutine get_surface_reaction_heat_flux( &", + " density, temperature, mass_fractions, q_rxn)", + "", + " GPU_ROUTINE(get_surface_reaction_heat_flux)", + "", + " real(wp), intent(in) :: density", + " real(wp), intent(in) :: temperature", + f" real(wp), intent(in) :: mass_fractions({len(gas_species)})", + " real(wp), intent(out) :: q_rxn", + "", + f" real(wp) :: h0_rt({len(gas_species)})", + f" real(wp) :: concentrations({len(gas_species)})", + " real(wp) :: rate", + " real(wp) :: delta_h", + " real(wp) :: h_rt_nongas", + "", + " call get_species_enthalpies_rt(temperature, h0_rt)", + "", + " concentrations = 0._wp", + ] + ) + + for i, species in enumerate(sol.species()): + lines.append(f" concentrations({i + 1}) = " f"density*mass_fractions({i + 1})/" f"{species.molecular_weight:.16e}_wp") + + lines.extend( + [ + "", + " q_rxn = 0._wp", + "", + ] + ) + + for reaction_number, reaction in enumerate(surface.reactions(), start=1): + append_reaction_rate(lines, reaction_number, reaction) + lines.append(" delta_h = 0._wp") + + # Gas species contribution to reaction enthalpy: + # h_k = (h_k/RT) * R * T, with R in J/(kmol K). + for name in gas_species: + nu = reaction.products.get(name, 0.0) - reaction.reactants.get(name, 0.0) + if nu == 0.0: + continue + + lines.append(f" delta_h = delta_h + ({nu:.16e}_wp)" f"*h0_rt({gas_index[name]})" "*gas_constant*temperature") + + # Non-gas species contribution. For Kestel this is graphite. + for name in nongas_species: + nu = reaction.products.get(name, 0.0) - reaction.reactants.get(name, 0.0) + if nu == 0.0: + continue + + if name not in nongas_h_rt: + raise common.MFCException(f"No supported thermochemistry is available for " f"non-gas surface-reaction species '{name}'.") + + tmid, low_expr, high_expr = nongas_h_rt[name] + + lines.extend( + [ + f" if (temperature <= {tmid:.16e}_wp) then", + f" h_rt_nongas = {low_expr}", + " else", + f" h_rt_nongas = {high_expr}", + " end if", + f" delta_h = delta_h + ({nu:.16e}_wp)" "*h_rt_nongas*gas_constant*temperature", + ] + ) + + # Positive q_rxn denotes heat released by an exothermic reaction. + lines.append(" q_rxn = q_rxn - delta_h*rate") + lines.append("") + + lines.extend( + [ + " end subroutine get_surface_reaction_heat_flux", + "", + "end module m_surface_thermochem", + "", + ] + ) + + return "\n".join(lines) + def generate_fpp(self, target) -> None: # Lazy import to avoid slow startup for commands that don't need chemistry import pyrometheus as pyro @@ -95,11 +479,35 @@ def generate_fpp(self, target) -> None: # Write the generated Fortran code to the m_thermochem.f90 file with the chosen precision sol = self.get_cantera_solution() + surface = self.get_cantera_surface() if target.name == "simulation" else None + if surface is not None: + cons.print(f"Loaded Cantera surface phase '{surface.name}' " f"with {surface.n_reactions} reaction(s).") thermochem_code = pyro.FortranCodeGenerator().generate("m_thermochem", sol, pyro.CodeGenerationOptions(scalar_type=real_type, directive_offload=directive_str)) common.file_write(os.path.join(modules_dir, "m_thermochem.f90"), thermochem_code, True) + if target.name == "simulation": + surface_thermochem_code = self.generate_surface_thermochem( + sol, + surface, + directive_str, + ) + + surface_thermochem_path = os.path.join( + modules_dir, + "m_surface_thermochem.f90", + ) + + common.file_write( + surface_thermochem_path, + surface_thermochem_code, + True, + ) + + if surface is not None: + cons.print(f"Generated m_surface_thermochem.f90 with " f"{surface.n_reactions} surface reaction(s).") + cons.unindent() def validate_constraints(self, target) -> None: