Appendix B: Transient Solver Mathematical Foundation
This appendix provides a high-level overview of the physical and mathematical foundation of the R-THYM transient hydraulic engine.
The R-THYM web application is directly integrated with the open-source RTHYM-MOC engine. This ensures that the solver physics, numerical accuracy, and boundary condition logic remain mathematically equivalent across both the JavaScript web environment and the C++/Python engines.
For detailed mathematical derivations—including the governing partial differential equations, Courant grid discretization, unsteady friction (Brunone) filter formulations, and the specific quadratic proofs for all boundary conditions (valves, pumps, and hydropneumatic tanks)—please refer to the complete theoretical documentation hosted in the RTHYM-MOC GitHub repository.
1. Solver Architecture & Scope
The engine utilizes the Method of Characteristics (MOC), a widely established numerical technique for computing the propagation of pressure waves through pressurized pipe networks. The solver is designed to accurately model:
- 1-D flow (cross-sectional average velocities)
- Elastic pipe walls and slightly compressible liquids
- Fully turbulent, single-phase liquid conditions
- Relatively small velocity changes (the acoustic approximation)
Internal Unit System
All calculations are performed internally using US customary units. The API boundary automatically handles conversion from user-facing inputs (such as GPM or psi) to the internal system:
| Quantity | Internal Unit | API Boundary Equivalent |
|---|---|---|
| Length, diameter | ft | inches (for diameter) |
| Piezometric head | ft | - |
| Velocity | ft/s | - |
| Flow | ft³/s (CFS) | GPM |
| Pressure | ft HGL | psi |
| Gravitational acceleration g | 32.2 ft/s² | - |
2. Fundamental Benchmark: The Joukowsky Relation
While the complete system of MOC compatibility equations (C+ and C-) is handled by the open-source library, the theoretical maximum pressure surge for an instantaneous flow stoppage remains the primary physical benchmark used to validate the solver.
This instantaneous pressure change is governed by the Joukowsky equation:
$$\Delta H = -\frac{a}{g} \Delta V$$For a complete, instantaneous valve closure from an initial velocity $V_0$, the maximum head rise is:
$$\Delta H = \frac{a \cdot V_0}{g}$$Where: * $H$ = piezometric head (ft) * $V$ = cross-sectional average velocity (ft/s) * $a$ = pressure wave speed (ft/s) * $g$ = gravitational acceleration (32.2 ft/s²)
By default, when no specific pipe material elasticity ($E$) is defined, the engine applies a conservative rigid-pipe approximation with a wave speed of $a = 4{,}000 \text{ ft/s}$.
3. MOC Simulation Settings
The four controls on Simulation → MOC Simulation Settings... (MOC mode only) change how the open-source engine builds its grid and how it treats vapor cavities. They do not apply to EPS. You cannot change them while a run is locked — click Reset first.
3.1 Simulation Grid Timestep
The slider sets the MOC time step $\Delta t$. Available values:
| Slider | $\Delta t$ | Typical use |
|---|---|---|
| Coarse | 0.01 s | Default. Fast screening when cavitation is not the study. |
| 0.005 s, 0.002 s | Intermediate. | |
| Fine | 0.001 s | Minimum $\Delta t$ that can run DVCM. |
| 0.0005 s, 0.0002 s | Finer collapse tracking. | |
| High-Fi | 0.0001 s | Severe column separation / secondary collapse spikes. |
A smaller $\Delta t$ resolves faster valve and pump events and (with DVCM) cavity collapse. It also multiplies the number of time steps and, unless you cap segments, the number of spatial points.
The characteristic grid must satisfy the Courant–Friedrichs–Lewy condition. For a pipe of length $L$ and design wave speed $a$:
$$C_r = \frac{a \cdot \Delta t}{\Delta x}$$How $C_r$ is enforced depends on Grid Discretization Mode. As a rule of thumb, $\Delta t$ should be at least an order of magnitude smaller than the fastest boundary event (valve stroke or pump trip). For DVCM, start at $0.001\text{ s}$ or finer and halve $\Delta t$ until peak pressures change by less than about 1%.
3.2 Cavitation Model
When local head falls to vapor pressure, the two models differ:
| Standard Clamping (Legacy) | High-Fidelity Column Separation (DVCM) | |
|---|---|---|
| What happens at vapor pressure | Head is clamped to the vapor floor. No cavity volume is stored. | The node switches to a vapor-cavity regime. Pocket volume is integrated: $V_c = \int (Q_\text{out} - Q_\text{in})\,dt$. |
| Column collision | Not modeled. Secondary collapse overpressure is missing. | When $V_c$ returns to zero the liquid columns rejoin and a collapse spike is computed. |
| Timestep | Stable at coarse $\Delta t$ (including 0.01 s). | Requires $\Delta t \le 0.001\text{ s}$. Coarser values revert to legacy clamping. |
| When to use | Pressures stay above vapor, or you only need a floor for screening. | Pump trips, slam closures, and surge-protection sizing where collapse peaks matter. |
Legacy clamping is a safeguard against non-physical sub-vapor heads. It is not a substitute for DVCM when the design question is “how hard do the columns hit when the cavity collapses?”
[!IMPORTANT] If you select DVCM and then move the timestep slider coarser than 0.001 s, R-THYM disables DVCM and shows the amber notice on the dialog. Choose Fine or High-Fi before relying on collapse spikes.
3.3 Grid Discretization Mode
Each pipe must be split into an integer number of segments $N$. The two modes trade wave-travel accuracy against numerical damping.
Courant Interpolation (variable $C_r \le 1$) — R-THYM default.
- Design wave speed $a$ is kept (no systematic $a$ distortion).
- Characteristics that do not land on a grid node are interpolated along the reach ($C_r \le 1$).
- Interpolation adds a small amount of numerical damping, so peaks can be slightly lower than a $C_r = 1$ grid.
- Use for screening, large models, and when travel time (phase) matters more than the last psi of peak.
- Pipes shorter than $a \cdot \Delta t$ fall back to one segment with $a$ scaled so $C_r = 1$.
Wave-Speed Adjustment (fixed $C_r = 1$).
- Target spacing $\Delta x = a \cdot \Delta t$.
- $N = \mathrm{round}(L / \Delta x)$ (at least 1, or 2 if a segment cap is on).
- Adjusted wave speed $a' = L / (N \cdot \Delta t)$ so $C_r = 1$ exactly.
There is no interpolation damping, so peak magnitudes are better preserved — the usual choice for safety-critical surge sizing. The cost is a small change in travel time. Distortion $|a' - a|/a$ of a few percent is normal; tens of percent means the cap or $\Delta t$ is fighting the pipe length (see §3.4).
3.4 Maximum Segments per Pipe
This cap limits $N$ per pipe so a fine $\Delta t$ does not explode the grid on long reaches.
| Value | Meaning |
|---|---|
| 100 (default) | Interactive default. Enough for typical plant piping; long transmission mains are coarsened. |
| 10–2000 | Raise the cap for more spatial detail, or lower it to save time. |
| 0 | Uncapped. $N$ follows $\Delta t$ and $a$ only. |
With a cap $N_\max > 0$:
$$N = \max\!\bigl(2,\ \min(N_\text{uncapped},\, N_\max)\bigr)$$Capping reduces spatial resolution. It does not fix a $\Delta t$ that is too coarse for DVCM volume integration. In Wave-Speed Adjustment, a tight cap plus a tiny $\Delta t$ can force $a' \gg a$ (large travel-time error). In Courant Interpolation, $a$ stays at design and the cap mainly coarsens $\Delta x$.
[!TIP] For a short, junction-heavy network, leave the cap at 100 (or 0). For a long line at High-Fi $\Delta t$, keep the cap (or interpolation) so the browser stays interactive, then remove the cap or refine $\Delta t$ for the sign-off run.
The live System Overview sidebar repeats the active grid mode, segment cap, and cavitation model during an MOC run.
4. Other Numerical Parameters
These are used by the engine but are not on the MOC Simulation Settings dialog. Defaults are conservative:
| Parameter | Symbol | Default | Notes |
|---|---|---|---|
| Vapor pressure | $p_\text{vap}$ | −14.0 psi | Transient Acceptance Criteria field. Solver cavitation floor and the post-run acceptance check. |
| Boundary-layer time constant | $\tau_{BL}$ | 0.5 s | IIR filter constant for unsteady friction (USF). |
| Brunone USF coefficient | $k_\text{Bru}$ | 0.0 | 0 = steady friction only; typical calibrated range: 0.02–0.15 |
Further Reading: Governing PDEs, compatibility equations, USF filters, and boundary-condition algebra are in the RTHYM-MOC documentation.