Arrhenius temperature-dependent shear viscosity mu(T) = exp(C + D/T)

docs/documentation/case.md

@@ -489,6 +489,9 @@ A particle cloud is a compact specification of a bed of identical circular (2D)
 | `Re(1)` * | Real   | Shear viscosity of fluid.                      |
 | `Re(2)` * | Real   | Volume viscosity of fluid.                     |
 | `k_therm` | Real   | Thermal conductivity of fluid (Fourier heat conduction). |
+| `visc_model` | Integer | Shear viscosity model of fluid (0 constant, 1 Arrhenius \f$\mu(T)\f$). |
+| `visc_c`  | Real   | Arrhenius viscosity coefficient \f$C\f$ of fluid. |
+| `visc_d`  | Real   | Arrhenius viscosity coefficient \f$D\f$ of fluid. |
 | `cv`   ** | Real   | Sffened-gas parameter $c_v$ of fluid.          |
 | `qv`   ** | Real   | Stiffened-gas parameter $q$ of fluid.          |
 | `qvp`  ** | Real   | Stiffened-gas parameter $q'$ of fluid.         |
@@ -513,6 +516,8 @@ Details of implementation of viscosity in MFC can be found in \cite Coralic15.
 
 - `fluid_pp(i)%%k_therm` sets the thermal conductivity of the $i$-th fluid, in units consistent with the rest of the (non-dimensional) case. A positive value on any fluid activates Fourier heat conduction, which adds \f$\nabla\cdot(k\nabla T)\f$ to the energy equation using the thermal-equilibrium mixture temperature and \f$k = \sum_i \alpha_i k_i\f$ (see @ref equations "Equations"). It requires `fluid_pp(i)%%cv` to be positive on every fluid that sets it (the mixture temperature is undefined without \f$c_v\f$), `model_eqns = 2` or `model_eqns = 3` (the mixture conductivity is weighted by volume fractions that `model_eqns = 1` does not carry), and `fluid_pp(i)%%eos` to be the stiffened-gas or ideal-gas equation of state. Heat conduction is independent of `viscous`: it can be enabled in an otherwise inviscid run. It is not supported with `igr`, nor with `chemistry` (which already carries its own mixture-averaged conduction through `chem_params%%diffusion`).
 
+- `fluid_pp(i)%%visc_model = 1` makes the shear viscosity of the $i$-th fluid temperature dependent, \f$\mu_i(T) = \exp(C_i + D_i/T)\f$ with \f$C_i\f$ = `fluid_pp(i)%%visc_c` and \f$D_i\f$ = `fluid_pp(i)%%visc_d` (nonzero), evaluated at the thermal-equilibrium mixture temperature that heat conduction uses. It replaces \f$1/\mathrm{Re}_i(1)\f$ in the mixture \f$1/\mathrm{Re} = \sum_i \alpha_i \mu_i\f$ on the left and right Riemann states, so the face viscosity follows the reconstructed temperature, and the viscous time-step limit sees the cell value. `fluid_pp(i)%%Re(1)` must still be set (`viscous` requires it, as for `non_newtonian`) but its value is unused for this fluid; `Re(2)` (bulk viscosity) stays constant. It requires `viscous = T`, `model_eqns = 2` or `3`, `riemann_solver = 1` or `2`, `fluid_pp(i)%%cv > 0` on every fluid, and the stiffened-gas or ideal-gas equation of state. It is not supported with `non_newtonian`, `cyl_coord`, `bubbles_euler`, `ib`, `igr`, or `chemistry`. The exponent \f$C_i + D_i/T\f$ is capped at 50 so that \f$\mu\f$ stays finite as \f$T \to 0\f$. At an isothermal wall (`isothermal_in`/`isothermal_out`) the wall-face viscosity is evaluated from the reflected ghost state, i.e. the adjacent interior temperature, not from `Twall`.
+
 - `fluid_pp(i)%%cv`, `fluid_pp(i)%%qv`, and `fluid_pp(i)%%qvp` define $c_v$, $q$, and $q'$ as parameters of $i$-th fluid that are used in stiffened gas equation of state.
 
 - `fluid_pp(i)%%G` is required for `hypoelasticity`.

docs/module_categories.json

@@ -23,6 +23,7 @@
             "m_eos",
             "m_viscous",
             "m_hb_function",
+            "m_arrhenius_viscosity",
             "m_surface_tension",
             "m_reactive_burn",
             "m_bubbles",

examples/2D_Couette_Variable_Viscosity/.gitignore

@@ -0,0 +1,6 @@
+# MFC run outputs and the per-grid run copies live here; regenerate with run_suite.sh.
+runs/
+*.log
+# validate.py output
+figures/
+summary.json

examples/2D_Couette_Variable_Viscosity/README.md

@@ -0,0 +1,155 @@
+# 2D Couette flow with temperature-dependent viscosity
+
+A unit-level validation of MFC's Arrhenius viscosity model
+`fluid_pp%visc_model = 1`, `mu(T) = exp(C + D/T)`, against an **exact** steady
+solution. This is the analytical gold standard for the variable-viscosity
+feature: it pins down the spatial order of accuracy on a problem with a known
+closed-form answer, so a sign error or discretization bug in the variable-`mu`
+viscous flux shows up directly.
+
+## The physics, and why it tests mu(T)
+
+A single fluid fills the gap `0 <= y <= H` between two no-slip walls:
+
+- bottom wall (`y = 0`): fixed, held at temperature `T0 = 300 K` (cold),
+- top wall (`y = H`): sliding at speed `U`, held at temperature `T1 = 400 K` (hot),
+- periodic in the streamwise (`x`) direction.
+
+Conduction sustains a temperature gradient across the gap, and because
+`mu(T) = exp(C + D/T)` with `D > 0`, the viscosity **falls as temperature rises**
+(the liquid/Arrhenius law). The hot fluid near the moving wall is therefore
+*thinner* than the cold fluid near the fixed wall.
+
+In steady state the shear stress `mu(T) du/dy` is uniform across the gap (no
+pressure gradient, no body force). Uniform stress with a non-uniform `mu` forces
+`du/dy` to vary inversely with `mu`: the velocity climbs slowly through the cold,
+viscous lower half and steeply through the hot, thin upper half. **The velocity
+profile is curved.** A constant-viscosity solver gives a straight line
+`u = U y/H` with `u(H/2) = U/2`; here the exact mid-gap velocity is
+`u(H/2) ≈ 0.352 U`. That ~30% deficit *is* the `mu(T)` signal, and it must match
+the exact profile to within the discretization error.
+
+## The exact solution (`reference.py`)
+
+The incompressible steady state solves the coupled two-point BVP
+
+```
+momentum:  d/dy [ mu(T) du/dy ] = 0
+energy:    d/dy [ k  dT/dy ] + mu(T) (du/dy)^2 = 0
+```
+
+with `u(0)=0, u(H)=U, T(0)=T0, T(H)=T1`. `reference.py` solves it with
+`scipy.solve_bvp` to ~1e-10. The viscous-heating term `mu(T)(du/dy)^2` couples
+the energy equation back to the momentum equation, so the reference also exercises
+the energy side of `mu(T)`. At this case's Brinkman number (`Br ≈ 0.014`) the
+heating bump is small (~0.17 K above the linear profile), and the closed-form
+`Br -> 0` limit (linear `T`, `u = U * integral_0^y dy'/mu / integral_0^H dy'/mu`)
+agrees with the coupled solver to ~0.02 % of `U` -- an independent check on the
+solver.
+
+## Compressible-solver setup
+
+MFC is a compressible code; temperature is not stored but follows from the
+stiffened-gas EOS `T = (p + p_inf)/((gam-1) rho cv)`. Two consequences shape the
+deck (`case.py`, constants in `couette_config.py`):
+
+- **Soft sound speed.** A deliberately small `cv` puts the sound speed near
+  `c ≈ 10` so the flow sits at `Ma ≈ 0.1`. Low Mach keeps the incompressible
+  reference faithful and the acoustic time step affordable.
+- **Temperature via a density profile.** The target linear `T(y)` is imposed by
+  initializing a `y`-varying density `rho(y) = (p0+p_inf)/((gam-1) cv T(y))` at
+  uniform pressure -- a valid mechanical equilibrium (no gravity).
+- **Modes.** The deck uses `model_eqns = 3` (6-equation) and `riemann_solver = 2`
+  (HLLC); `visc_model = 1` also runs with `model_eqns = 2` and HLL. MFC evaluates
+  `mu(T)` at the mixture temperature of the reconstructed left and right states at
+  each cell face, so the face viscosity follows the reconstructed temperature.
+  Fourier conduction (`fluid_pp(1)%k_therm > 0`) with isothermal no-slip walls
+  (`bc_y = -16`, `isothermal_in/out`) holds the walls at `T0` and `T1`.
+- **Compact viscous gradient (`weno_Re_flux = F`, `weno_avg = F`).** `weno_Re_flux`
+  is designed for *discontinuous* viscosity at material interfaces. For a single
+  fluid with smooth `mu(T)` it leaves the `2*dy` (odd-even) mode of the
+  wall-normal velocity undamped -- and because this flow has no advection
+  (`v ~ 0`), the HLLC upwinding has nothing to damp either -- so `u(y)` grows a
+  grid-scale checkerboard (~5 % of `U` at `Ny = 33`). The compact central
+  gradient damps it; the resolved profile then matches the exact solution to the
+  discretization error. (The two flags must be set together: `weno_avg = F` with
+  `weno_Re_flux = T` disables the viscous coupling entirely.)
+
+The flow starts from rest (`u = 0`) and develops toward the curved steady profile,
+so the match is *not* baked into the initial condition.
+
+Dimensionless parameters: `Re ≈ 20`, `Ma ≈ 0.1`, `Pr = 1`, `Br ≈ 0.014`,
+viscosity contrast `mu(T0)/mu(T1) ≈ 3.5`.
+
+## Running it
+
+```bash
+# 1. Run the grid-refinement sweep. The first grid is a normal building run
+#    (the analytic density IC compiles into pre_process, so the build is
+#    case-specific); the remaining grids then reuse that binary (--no-build)
+#    and run concurrently, each in its own runs/n<N>/ on one rank.
+#    Default grids: Ny = 33, 65, 97.
+cd examples/2D_Couette_Variable_Viscosity
+./run_suite.sh                 # or: ./run_suite.sh 32 64 96 128
+
+# 2. Compare to the exact solution; writes summary.json + figures/
+python3 validate.py
+```
+
+A single grid can be run directly, choosing the wall-normal resolution with an
+environment variable:
+
+```bash
+COUETTE_N=64 ./mfc.sh run examples/2D_Couette_Variable_Viscosity/case.py -n 1
+```
+
+## Files
+
+| file | role |
+|------|------|
+| `couette_config.py` | single source of truth for all physical constants |
+| `reference.py`      | exact coupled-BVP solution (`scipy.solve_bvp`) |
+| `case.py`           | MFC input deck, parametrized by `COUETTE_N` |
+| `run_suite.sh`      | runs the grid sweep into `runs/n<N>/` |
+| `validate.py`       | compares MFC to the exact solution; figures + `summary.json` |
+
+## Results
+
+The MFC steady velocity matches the exact `mu(T)` Couette profile and converges at
+second order; temperature is exact to ~1e-6 and the runs are fully steady
+(`unsteadiness ~ 1e-9`). The viscosity contrast bends `u(H/2)` to **0.352** versus
+the constant-viscosity **0.500** -- the signal a constant-`mu` code cannot produce.
+
+| Ny | L2(u)/U  | Linf(u)/U | L2(T)/dT |
+|----|----------|-----------|----------|
+| 33 | 1.33e-4  | 1.83e-4   | 1.58e-6  |
+| 65 | 3.41e-5  | 4.85e-5   | 3.93e-7  |
+| 97 | 1.53e-5  | 2.20e-5   | 1.75e-7  |
+
+Observed spatial order (L2 velocity) = **2.00**. `validate.py` writes to `figures/`: profile
+overlay (`couette_profiles.png`), L2 convergence (`couette_convergence.png`), and
+the spatial relative error (`couette_error.png`) -- the latter shows the
+pointwise `(u - u_exact)/U` shrinking with refinement and collapsing under an
+`Ny^2` rescaling, i.e. the second-order accuracy holds everywhere in the gap, not
+just in the L2 norm, and the full data to `summary.json`.
+
+## Scope and the companion benchmark
+
+This case validates the **exponential/Arrhenius (liquid) law** that MFC actually
+implements -- the right gold standard for temperature-dependent droplet
+viscosity. It is a unit test: `mu(T)` is isolated on an exact answer, free of the
+EOS/buoyancy coupling that entangles system-level benchmarks.
+
+The natural *system-level* companion is the low-Mach differentially-heated cavity
+(Le Quéré / Vierendeels), where boundary-layer asymmetry is driven by the
+viscosity contrast and validated against tabulated wall Nusselt numbers. That
+benchmark uses **Sutherland's law (a gas law)**, which MFC does not yet
+implement, so it is a separate follow-up (add the Sutherland model, then the
+cavity case).
+
+### Cost note
+
+The time step is explicit-stability limited: acoustic (`~dy`) on coarse grids and
+**viscous (`~dy^2`)** once `Ny` exceeds ~60. The default sweep stays in the cheap
+range; `Ny = 129` is viscous-CFL limited (~221k steps, hours on a single rank) and
+is optional.

examples/2D_Couette_Variable_Viscosity/case.py

@@ -0,0 +1,145 @@
+#!/usr/bin/env python3
+"""Variable-viscosity planar Couette flow -- MFC input deck.
+
+Validates the temperature-dependent (Arrhenius) viscosity feature
+fluid_pp%visc_model = 1, mu(T) = exp(C + D/T), against the exact steady solution
+in reference.py. See README.md for the full description.
+
+Setup (all constants from couette_config.py):
+  - 2D channel, periodic in x, walls top and bottom in y.
+  - Bottom wall (y=0): no-slip, fixed, isothermal at T0 (cold -> viscous).
+  - Top wall   (y=H): no-slip, sliding at U, isothermal at T1 (hot -> thin).
+  - Single fluid; model_eqns = 3 + riemann_solver = 2 (HLLC).
+  - viscous on, and Fourier conduction (fluid_pp(1)%k_therm > 0) pins the wall
+    temperatures and sustains the temperature gradient that drives mu(T).
+  - Initial condition: quiescent (u=0), uniform pressure, with a y-varying
+    density that encodes the linear conduction temperature profile via the EOS.
+    The flow then develops from rest toward the curved steady Couette profile.
+
+The wall-normal grid n is set by the COUETTE_N environment variable so the
+validation harness can run a spatial-refinement sweep on a single deck.
+"""
+
+import json
+import os
+
+import couette_config as cfg
+
+# Wall-normal resolution (refined by the convergence sweep); streamwise grid is
+# fixed and small because there is no x-gradient.
+# WENO5 needs >= 25 cells per dimension, so the streamwise (x) grid is pinned at
+# that minimum (there is no x-gradient) and the sweep keeps n >= 32.
+n_cells = int(os.environ.get("COUETTE_N", "64"))
+m_cells = int(os.environ.get("COUETTE_M", "24"))
+
+# Square cells: size the (periodic, gradient-free) x-extent so dx == dy on every
+# grid. This keeps the acoustic CFL set by dy alone and the time step ~ 1/n.
+dy = cfg.H / (n_cells + 1)
+Lx = dy * (m_cells + 1)
+dx = Lx / (m_cells + 1)
+
+# Time step: minimum of acoustic, viscous, and thermal explicit-stability limits.
+c_max = (cfg.gam * (cfg.gam - 1.0) * cfg.cv * cfg.T1) ** 0.5  # sound speed, hot wall
+dt_acoustic = 0.4 * min(dx, dy) / c_max
+dt_viscous = 0.4 * dy**2 / (2.0 * cfg.nu_max)
+alpha_max = cfg.k_therm / (cfg.rho_of_T(cfg.T0) * cfg.cp)  # thermal diffusivity, cold wall
+dt_thermal = 0.4 * dy**2 / (2.0 * alpha_max)
+dt = min(dt_acoustic, dt_viscous, dt_thermal)
+
+t_step_stop = int(round(cfg.t_end / dt))
+if "COUETTE_NSTEPS" in os.environ:  # short override for timing probes
+    t_step_stop = int(os.environ["COUETTE_NSTEPS"])
+t_step_save = max(1, t_step_stop // 40)  # ~40 snapshots to confirm steadiness
+
+# Linear conduction temperature profile T(y) = T0 + (T1-T0)*y/H imposed through
+# the density proxy rho(y) = (p0+p_inf)/((gam-1)*cv*T(y)). Written with plain
+# decimals (no 'e' literals -- MFC would expand a bare 'e' to Euler's number).
+gm1_cv = (cfg.gam - 1.0) * cfg.cv
+slope = (cfg.T1 - cfg.T0) / cfg.H
+rho_expr = f"({cfg.p0 + cfg.p_inf})/(({gm1_cv})*({cfg.T0} + ({slope})*y))"
+
+case = {
+    # -- Run info --
+    "run_time_info": "F",
+    # -- Domain (periodic x, walls in y) --
+    "x_domain%beg": 0.0,
+    "x_domain%end": Lx,
+    "y_domain%beg": 0.0,
+    "y_domain%end": cfg.H,
+    "m": m_cells,
+    "n": n_cells,
+    "p": 0,
+    "cyl_coord": "F",
+    # -- Time stepping --
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": t_step_stop,
+    "t_step_save": t_step_save,
+    # -- Model: 6-equation + HLLC (visc_model = 1 also runs with model_eqns = 2 or HLL) --
+    "model_eqns": 3,
+    "num_fluids": 1,
+    "alt_soundspeed": "F",
+    "mpp_lim": "F",
+    "mixture_err": "T",
+    "time_stepper": 3,
+    # -- Spatial scheme --
+    "weno_order": 5,
+    "weno_eps": 1e-16,
+    "mapped_weno": "T",
+    "null_weights": "F",
+    "mp_weno": "F",
+    # Compact (non-WENO) viscous gradient. weno_Re_flux is meant for discontinuous
+    # viscosity at material interfaces; for a single fluid with smooth mu(T) it
+    # leaves the 2*dy (odd-even) mode of the wall-normal velocity undamped -- in
+    # this advection-free flow there is no convective upwinding to remove it -- so
+    # u(y) develops a grid-scale checkerboard. The compact central gradient damps
+    # it and the profile matches the exact solution. (Both must be set together:
+    # weno_avg=F with weno_Re_flux=T disables the viscous coupling entirely.)
+    "weno_Re_flux": os.environ.get("COUETTE_WENO_RE", "F"),
+    "weno_avg": os.environ.get("COUETTE_WENO_AVG", "F"),
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    # -- Boundary conditions: periodic x; no-slip isothermal walls in y --
+    "bc_x%beg": -1,
+    "bc_x%end": -1,
+    "bc_y%beg": -16,  # no-slip wall, bottom (fixed)
+    "bc_y%end": -16,  # no-slip wall, top (sliding)
+    "bc_y%ve1": cfg.U,  # top-wall x-velocity (the Couette drive)
+    "bc_y%isothermal_in": "T",
+    "bc_y%Twall_in": cfg.T0,
+    "bc_y%isothermal_out": "T",
+    "bc_y%Twall_out": cfg.T1,
+    # -- Output --
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "T",
+    # -- Physics --
+    "viscous": "T",
+    "chemistry": "F",
+    # -- Fluid: ideal gas with Arrhenius mu(T) = exp(C + D/T) --
+    "fluid_pp(1)%gamma": 1.0 / (cfg.gam - 1.0),
+    "fluid_pp(1)%pi_inf": cfg.gam * cfg.p_inf / (cfg.gam - 1.0),
+    "fluid_pp(1)%cv": cfg.cv,
+    "fluid_pp(1)%Re(1)": 1.0 / cfg.mu_ref,  # nominal Re; visc_model=1 overrides with mu(T)
+    "fluid_pp(1)%k_therm": cfg.k_therm,
+    "fluid_pp(1)%visc_model": 1,
+    "fluid_pp(1)%visc_c": cfg.C,
+    "fluid_pp(1)%visc_d": cfg.D,
+    # -- Initial condition: quiescent, uniform p, linear-T density profile --
+    "num_patches": 1,
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%x_centroid": Lx / 2,
+    "patch_icpp(1)%y_centroid": cfg.H / 2,
+    "patch_icpp(1)%length_x": Lx,
+    "patch_icpp(1)%length_y": cfg.H,
+    "patch_icpp(1)%vel(1)": 0.0,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%pres": cfg.p0,
+    "patch_icpp(1)%alpha_rho(1)": rho_expr,
+    "patch_icpp(1)%alpha(1)": 1.0,
+}
+
+if __name__ == "__main__":
+    print(json.dumps(case))