Mean gas temperature source: mean_temp_wall (thermally developed flow) and immersed-boundary support

docs/documentation/case.md

@@ -1342,7 +1342,7 @@ This forcing is implemented for **2D cases only**; enabling it in 1D or 3D is re
 These drive a flow through a periodic direction with a mean pressure gradient, as in particle-resolved random-array studies (e.g. Tenneti, Garg & Subramaniam, IJMF 2011).
 Unlike `bf_x[y,z]`, the source \f$f_i = -\partial\bar p/\partial x_i\f$ is a force per unit volume, independent of density.
 It is added to the momentum equation, with work \f$\mathbf{u}\cdot\mathbf{f}\f$ added to the energy equation.
-Immersed boundaries are not supported with the mean pressure gradient.
+With immersed boundaries the force is not applied inside the bodies, and the fluid volume, mass flux and work in the controller sum over fluid cells only.
 
 When `mean_mflux(i)` is set, \f$f_i\f$ is updated once per time step to hold the fluid-averaged bulk mass flux \f$Q_i = \sum \rho u_i\,dV / \sum dV\f$ at \f$Q_i^\ast\f$:
 
@@ -1358,11 +1358,12 @@ Each forced direction must be periodic (`bc_x[y,z]%%beg = bc_x[y,z]%%end = -1`)
 | Parameter        | Type | Description                                                              |
 | ---:             | :---: | :---                                                                    |
 | `mean_qdot`      | Real | Uniform heat source per unit volume \f$\dot q\f$ (initial value if controlled; default 0) |
-| `mean_temp`      | Real | Target mass-averaged fluid temperature \f$T_b^\ast\f$; enables the controller |
+| `mean_temp`      | Real | Target mass-averaged fluid temperature \f$T_b^\ast\f$ (the mixing-cup \f$T_m\f$ with `mean_temp_wall`); enables the controller |
 | `mean_temp_gain` | Real | Controller gain \f$G\f$, required with `mean_temp`, \f$0 < G < 4/3\f$; needs a constant `dt` (no `cfl_adap_dt`) |
+| `mean_temp_wall` | Real | Wall temperature \f$T_s\f$ of a thermally developed state; makes the source local (requires `mean_temp`) |
 
 The thermal analogue of the mean pressure gradient: a uniform volumetric heat source that holds the bulk gas temperature of a periodic domain (statistically stationary heat transfer, cf. the thermally fully developed periodic formulation of Tenneti, Sun, Garg & Subramaniam, IJHMT 58, 2013, and Sun, Tenneti & Subramaniam, IJHMT 86, 2015).
-\f$\dot q\f$ is added to the energy equation. Immersed boundaries are not supported with the thermostat.
+\f$\dot q\f$ is added to the energy equation. With immersed boundaries the source is not applied inside the bodies, and the controller sums run over fluid cells only.
 It requires `chemistry`; the controlled quantity is the temperature, so gas-phase heat release lowers the source needed and the controller holds \f$T_b^\ast\f$ net of chemistry (\f$\dot q\f$ may become negative).
 
 When `mean_temp` is set, \f$\dot q\f$ is updated once per time step from the mass-averaged fluid temperature \f$T_b = \sum \rho T\,dV / \sum \rho\,dV\f$, the quantity the energy budget acts on:
@@ -1372,7 +1373,20 @@ When `mean_temp` is set, \f$\dot q\f$ is updated once per time step from the mas
 with \f$\overline{\rho c_v} = \sum \rho c_v\,dV / \sum dV\f$, the same controller as `mean_mflux` (deadbeat for \f$G = 1\f$ and a constant heat loss).
 \f$t\f$, \f$T_b\f$, the mixing-cup temperature \f$T_m = \sum \rho u_1 T\,dV / \sum \rho u_1\,dV\f$, \f$\dot q\f$, and the fluid volume \f$V_f\f$ (source power \f$\dot q V_f\f$) are written every step to `mean_temp.dat` in the case directory.
 As for `mean_dpdx`, the controller state is not saved in restart files: set `mean_qdot` on restart to the last value in `mean_temp.dat`.
-All active directions must be periodic and Cartesian.
+The domain must be Cartesian and periodic in \f$x\f$; \f$y\f$ and \f$z\f$ may be periodic or walls (`-15`, `-16`).
+
+With `mean_temp_wall`, the source instead gives the thermally fully developed state of flow past isothermal walls at \f$T_s\f$ (Sun, Tenneti & Subramaniam, IJHMT 86, 2015): \f$T = T_s + \psi\, e^{-\beta x}\f$ with \f$\psi\f$ periodic.
+MFC solves for \f$\psi = T - T_s\f$ with the source of the \f$e^{-\beta x}\f$ factor,
+
+\f[ S = \beta\left[\rho c_p u_1 \psi - \partial_x(\lambda\psi) - \lambda\,\partial_x\psi + \beta\lambda\psi\right], \f]
+
+with local \f$c_p\f$ and conductivity \f$\lambda\f$ and central differences in \f$x\f$.
+The controller (on \f$T_m\f$, below) sets the source total \f$\dot q V_f\f$; \f$\beta\f$ is the root of \f$\dot q V_f = \beta(A_1 + \beta A_2)\f$, with \f$A_1 = \sum [\rho c_p u_1 \psi - \partial_x(\lambda\psi) - \lambda\,\partial_x\psi]\,dV\f$ (the same discrete terms the source applies, evaluated once per step at its first stage and held through the Runge-Kutta stages like \f$\dot q\f$ and \f$\beta\f$, so the balance holds with immersed boundaries and varying properties) and \f$A_2 = \sum \lambda\psi\,dV\f$, that tends to \f$\dot q V_f/A_1\f$ as \f$A_2 \to 0\f$. It requires `chem_params%%diffusion`, `mean_mflux(1) > 0` (a held streamwise flow for the mixing-cup temperature), and every wall (`bc_y,z` of `-15`, `-16`, and every immersed boundary) isothermal at `mean_temp_wall`.
+At steady state, \f$\beta\f$ is the streamwise decay rate of \f$T - T_s\f$ and the walls' heat equals \f$-\dot q V_f\f$.
+\f$\beta\f$, \f$A_1\f$ and \f$A_2\f$ are written as columns 7-9 of `mean_temp.dat`.
+In this mode the controller holds the mixing-cup temperature \f$T_m\f$ at `mean_temp`, not \f$T_b\f$: in flows with streamlines of very different speed (e.g. periodic particle arrays) the slowest-decaying mode can sit on a fast, unobstructed streak that carries little mass but much of the flux, and holding the mass mean let it grow without bound, while the flux-weighted \f$T_m\f$ feeds it back.
+Keep \f$|T_s - T_m|\f$ small (a few percent of \f$T_m\f$): the cold or hot core of the developed state is several times deeper than \f$|T_s - T_m|\f$, and large amplitudes bias \f$\beta\f$ through the property changes.
+Known issue: only the plane channel (`examples/2D_thermal_channel`, \f$\mathrm{Nu} = 7.541\f$ for constant properties) is verified. In flows with streamlines of very different speed the hot and cold wall cases relax at different rates and \f$\beta\f$ is not converged; treat such results as unverified.
 
 ### 14. Magnetohydrodynamics (MHD)
 

docs/documentation/equations.md

@@ -811,6 +811,7 @@ With a target bulk mass flux \f$Q_i^\ast\f$, \f$f_i^{n+1} = f_i^n + (G/\Delta t)
 \f[\frac{\partial E}{\partial t}\bigg|_\text{mts} = \dot q, \qquad \dot q^{n+1} = \dot q^n + \frac{G\,\overline{\rho c_v}}{\Delta t}\left(T_b^\ast - 2T_b^n + T_b^{n-1}\right),\f]
 
 with \f$T_b\f$ the mass-averaged fluid temperature.
+With `mean_temp_wall` \f$= T_s\f$, \f$\dot q\f$ becomes the local \f$\beta[\rho c_p u_1 \psi - \partial_x(\lambda\psi) - \lambda\partial_x\psi + \beta\lambda\psi]\f$, \f$\psi = T - T_s\f$, with \f$\beta\f$ from \f$\dot q V_f = \beta(\sum\rho c_p u_1\psi\,dV + \beta\sum\lambda\psi\,dV)\f$ and the controller on the mixing-cup \f$T_m\f$: the periodic part of the thermally developed \f$T = T_s + \psi e^{-\beta x}\f$.
 
 ---
 

examples/2D_thermal_channel/case.py

@@ -0,0 +1,114 @@
+#!/usr/bin/env python3
+"""
+2D plane channel between isothermal no-slip walls (bc_y = -16) at Tw, periodic in x, driven at a bulk mass flux
+(mean_mflux). The thermostat with mean_temp_wall = Tw holds the gas in the thermally fully developed state
+T = Tw + psi(y) exp(-beta x), with the mixing-cup temperature at T0. Constant-property reference (Shah & London 1978):
+Nu = q_w D_h / (k (Tw - T_m)) = 7.541, D_h = 2H, T_m the mixing-cup temperature; axial conduction changes it by
+under 0.02 % for Pe_Dh > 100. q_w = -mean_qdot V_f / (2 L_x) at steady state (mean_temp.dat). A hot/cold wall pair
+(--tw above and below --t0) cancels the viscous heating to first order.
+"""
+
+import argparse
+import json
+import math
+
+import cantera as ct
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--ny", type=int, default=32, help="cells across the gap H")
+parser.add_argument("--nx", type=int, default=25, help="cells along x (WENO5 needs at least 25 per direction)")
+parser.add_argument("--tw", type=float, default=330.0, help="wall temperature [K]")
+parser.add_argument("--t0", type=float, default=300.0, help="initial gas and target mixing-cup temperature [K]")
+parser.add_argument("--mach", type=float, default=0.1, help="bulk Mach number")
+parser.add_argument("--re", type=float, default=150.0, help="Re_Dh = rho U_b 2H / mu")
+parser.add_argument("--tend", type=float, default=0.5, help="run time [H^2 / alpha]; developed in about 0.1")
+parser.add_argument("--saves", type=int, default=10)
+args, _ = parser.parse_known_args()
+
+gas = ct.Solution("h2o2.yaml")
+gas.TPY = args.t0, ct.one_atm, "O2:0.233,N2:0.767"
+rho, mu, k = gas.density, gas.viscosity, gas.thermal_conductivity
+a = math.sqrt(gas.cp_mass / gas.cv_mass * ct.one_atm / rho)
+U = args.mach * a
+H = args.re * mu / (2.0 * rho * U)
+dy = H / args.ny
+Lx = args.nx * dy
+dt = 0.4 * dy / (U + a)
+nt = args.saves * math.ceil(args.tend * H**2 * rho * gas.cp_mass / k / dt / args.saves)
+
+case = {
+    "run_time_info": "F",
+    "x_domain%beg": 0.0,
+    "x_domain%end": Lx,
+    "y_domain%beg": 0.0,
+    "y_domain%end": H,
+    "m": args.nx - 1,
+    "n": args.ny - 1,
+    "p": 0,
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": nt,
+    "t_step_save": nt // args.saves,
+    "num_patches": 1,
+    "model_eqns": 2,
+    "num_fluids": 1,
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": 3,
+    "weno_order": 5,
+    "weno_eps": 1e-16,
+    "mapped_weno": "T",
+    "weno_Re_flux": "T",
+    "mp_weno": "F",
+    "weno_avg": "T",
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "fd_order": 4,
+    "bc_x%beg": -1,
+    "bc_x%end": -1,
+    "bc_y%beg": -16,
+    "bc_y%end": -16,
+    "bc_y%isothermal_in": "T",
+    "bc_y%isothermal_out": "T",
+    "bc_y%Twall_in": args.tw,
+    "bc_y%Twall_out": args.tw,
+    "viscous": "T",
+    "chemistry": "T",
+    "chem_params%diffusion": "T",
+    "chem_params%reactions": "F",
+    "chem_params%transport_model": 1,
+    "cantera_file": "h2o2.yaml",
+    "chem_wrt_T": "T",
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "F",
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%x_centroid": 0.5 * Lx,
+    "patch_icpp(1)%y_centroid": 0.5 * H,
+    "patch_icpp(1)%length_x": Lx,
+    "patch_icpp(1)%length_y": H,
+    "patch_icpp(1)%vel(1)": U,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%pres": ct.one_atm,
+    "patch_icpp(1)%alpha_rho(1)": rho,
+    "patch_icpp(1)%alpha(1)": 1.0,
+    "fluid_pp(1)%gamma": 1.0 / (1.4 - 1.0),
+    "fluid_pp(1)%eos": "stiffened_gas",
+    "fluid_pp(1)%pi_inf": 0.0,
+    "fluid_pp(1)%Re(1)": 1.0 / mu,
+    # Poiseuille drive and the developed wall loss (Nu 7.541) as initial controller states
+    "mean_dpdx(1)": -12.0 * mu * U / H**2,
+    "mean_mflux(1)": rho * U,
+    "mean_dpdx_gain": 0.05,
+    "mean_qdot": -7.541 * k * (args.tw - args.t0) / H**2,
+    "mean_temp": args.t0,
+    "mean_temp_gain": 0.05,
+    "mean_temp_wall": args.tw,
+}
+for i, y in enumerate(gas.Y):
+    case[f"patch_icpp(1)%Y({i + 1})"] = y
+    case[f"chem_wrt_Y({i + 1})"] = "T"
+
+print(json.dumps(case))