---
title: Harmonic-Balanced Navier–Stokes (HBNS)
url: https://www.emergentmind.com/topics/harmonic-balanced-navier-stokes-hbns
type: topic
---

# Harmonic-Balanced Navier–Stokes (HBNS)

Searching arXiv for the specified HBNS paper and closely related harmonic-balance Navier–Stokes work to ground the encyclopedia entry.
Harmonic-Balanced Navier–Stokes (HBNS) is a frequency-domain, stabilized finite-element formulation for incompressible Navier–Stokes flow that exploits strict time-periodicity to replace conventional time marching with a spectral discretization over one cycle. In the cardiovascular setting considered in "Introducing a Harmonic Balance Navier-Stokes Finite Element Solver to Accelerate Cardiovascular Simulations" [2411.14315], HBNS is designed for physically-stable time-periodic flows and targets patient-specific simulation workloads in diagnosis and surgical planning, where the computational cost of conventional transient solvers is a practical bottleneck. The method solves for all retained temporal harmonics simultaneously, uses a Galerkin/least-squares stabilization framework, and is reported to reduce wall-clock time by up to two orders of magnitude while maintaining excellent agreement with conventional solvers when the number of modes is sufficiently large to represent the imposed boundary conditions [2411.14315].

## 1. Governing setting and periodicity assumptions

The formulation begins from the incompressible Navier–Stokes equations in a spatial domain $\Omega$ over one cardiac cycle of period $T$, with fundamental frequency $\omega = 2\pi/T$:

$$
\rho \frac{\partial u}{\partial t} + \rho (u \cdot \nabla)u + \nabla p - \mu \Delta u = 0
\quad \text{in } \Omega \times (0,T],
$$

$$
\nabla \cdot u = 0
\quad \text{in } \Omega \times (0,T],
$$

with Dirichlet data $u=g$ on $\Gamma_g \times (0,T]$ and Neumann data $(-pn+\mu \partial u/\partial n)=h\,n$ on $\Gamma_h \times (0,T]$ [2411.14315]. The central assumption is strict periodicity,

$$
u(x,t+T)=u(x,t),
$$

together with an analogous periodic representation of the pressure field. This assumption is structural rather than incidental: the method is explicitly formulated for flows whose time dependence can be represented over a single cycle and repeated thereafter.

In the paper’s framing, the usefulness of this assumption is tied to cardiovascular flows with intrinsic periodicity. A plausible implication is that HBNS is most naturally deployed when the dominant unsteadiness is driven by a known physiological forcing frequency rather than by broadband, non-periodic, or instability-dominated dynamics. That interpretation is consistent with the stated scope of the method, which is limited to strict time-periodicity at known $\omega$ and Reynolds numbers below the laminar–turbulent transition [2411.14315].

## 2. Harmonic decomposition and frequency-domain reformulation

HBNS truncates the temporal field to $M$ modes, with $N \coloneqq 2M-1$ time samples, and expands the solution in Fourier modes:

$$
u(x,t) = \sum_{|m|<M} \hat{u}_m(x)e^{i m \omega t},
\qquad
p(x,t) = \sum_{|m|<M} \hat{p}_m(x)e^{i m \omega t}.
$$

The boundary data $g$ and $h$ are expanded in the same manner so as to interpolate their spectral content exactly up to mode $M-1$ [2411.14315]. Substituting these expansions into the governing equations and matching terms in $e^{i m \omega t}$ yields, for each harmonic index $m$,

$$
\rho (i m \omega)\hat{u}_m
+
\rho \sum_k (\hat{u}_k \cdot \nabla)\hat{u}_{m-k}
+
\nabla \hat{p}_m
-
\mu \Delta \hat{u}_m
=0,
$$

$$
\nabla \cdot \hat{u}_m = 0.
$$

In stacked form, over the coefficient vectors $U^*$ and $P^*$, this becomes

$$
\rho \Omega U^* + \rho A_j \frac{\partial U_i^*}{\partial x_j} + \nabla P^* - \mu \Delta U^* = 0,
\qquad
\nabla \cdot U^* = 0,
$$

where $\Omega = \operatorname{diag}(i(m-M)\omega)$ and $A_j$ is the convolution matrix $A_j(m,k)=\hat{u}_{m-k,j}$ [2411.14315].

A key computational observation is that direct assembly in harmonic space is $O(M^2)$. HBNS therefore applies the inverse discrete Fourier transform, $U = E^{-1}U^*$, to obtain a harmonic-balance form at $N$ real time points $t_n$:

$$
\rho H U_i + \rho U_j \frac{\partial U_i}{\partial x_j} + \nabla P - \mu \Delta U_i = 0,
\qquad
\nabla \cdot U = 0,
$$

where $H = E^{-1}\Omega E$ is a skew-symmetric time-coupling matrix and $U_j=\operatorname{diag}(U_j(t_n))$ [2411.14315]. In this representation, convection and diffusion act pointwise in time, whereas the temporal coupling is encoded by $H$. Because $H$ is applied through FFT/IFFT operations, the time-coupling cost is reduced to $O(N\log N)$.

This reformulation is one of the distinguishing features of HBNS. It preserves the nonlinear coupling of harmonics while avoiding direct quadratic-in-mode assembly, and it converts the temporal problem into a coupled system over a finite set of phase points within one period.

## 3. Stabilized weak form and finite-element discretization

The weak formulation is posed on the trial and test spaces

$$
U \in [H^1(\Omega)]^{3N}, \qquad P \in [L^2(\Omega)]^{N},
$$

$$
W \in [H^1_0(\Omega)]^{3N}, \qquad Q \in [L^2(\Omega)]^{N}.
$$

The Galerkin/least-squares form seeks $(U,P)$ such that, for all $(W,Q)$,

$$
\int_\Omega W^T[\rho H U + \rho U_j \partial U/\partial x_j]\,d\Omega
+
\int_\Omega \nabla W : (-PI+\mu \nabla U)\,d\Omega
+
\int_\Omega Q(\nabla \cdot U)\,d\Omega
$$

$$
+
\sum_{\text{elements } e}\int_e \tau\,L(W)^T L(U,P)\,d\Omega
=
\int_{\Gamma_h} W^T(h\,n)\,d\Gamma,
$$

where the momentum-operator residual is

$$
L(U,P)=\rho H U + \rho U_j \frac{\partial U}{\partial x_j} + \nabla P - \mu \Delta U
$$

and $\tau$ is a diagonal stabilization parameter, one value per $t_n$, defined by

$$
\tau_n =
\left[
\big(U(t_n)\cdot \xi\, U(t_n)\big)+C_I \nu^2 \xi:\xi
\right]^{-1/2},
$$

with $\xi$ the contravariant metric, $\nu=\mu/\rho$, and $C_I$ a shape constant, approximately $3$ for tetrahedra [2411.14315].

The stated role of the GLS term is twofold. First, it permits stable solution in convection-dominant flows. Second, it allows convenient use of the same interpolation functions for velocity and pressure. The paper further notes that the GLS contribution simultaneously provides SUPG and PSPG effects and formally recovers the steady-state GLS formulation in the limit $\omega \to 0$ [2411.14315].

Spatial discretization uses equal-order linear tetrahedral interpolation $N_A(x)$ for all components:

$$
U_i(x) \approx \sum_A N_A(x) U_{Ai},
\qquad
P(x) \approx \sum_A N_A(x) P_A.
$$

Dirichlet conditions are enforced by directly prescribing nodal values at the $N$ time points. The resulting global system has $4N$ unknowns per node [2411.14315]. This equal-order construction is significant because pressure stabilization is built into the formulation rather than enforced through mixed interpolation pairs.

## 4. Algebraic structure, nonlinear solution, and computational scaling

After assembly, the complex-valued tangent matrix $L$ is split as

$$
L = C + P,
$$

where $P$ has the same sparsity pattern as a steady-state Navier–Stokes Jacobian and is block-diagonal in time, with cost $O(N)$, whereas $C$ encodes the $H$-coupling through

$$
C_{AB} = \int_\Omega \rho N_A H N_B\,d\Omega.
$$

The $C$ blocks are dense in time but are applied via FFT/IFFT in $O(N\log N)$ [2411.14315]. The paper therefore states that overall assembly-and-solve cost scales nearly $O(N)$.

The nonlinear algebraic system is solved by Newton–Raphson. At iteration $k$,

$$
[C+P(U^{(k)})]\Delta y = -r(U^{(k)}),
\qquad
y^{(k+1)} = y^{(k)} + \Delta y,
$$

with residual vector $r$ and tangent $L$ [2411.14315]. To enhance convergence, a pseudo-time derivative

$$
\rho \frac{\partial U}{\partial \tilde{t}}
\approx
\rho \frac{U^{k+1}-U^k}{\Delta \tilde{t}}
$$

is added and treated with a generalized-$\alpha$ scheme. The pseudo-time step $\Delta \tilde{t}$ is chosen so that the convective CFL is $O(1)$. Each linear system is solved with GMRES using a Jacobi preconditioner on $P$, and iterations continue until $\|r\|$ drops by $10^{-3}$ [2411.14315].

From an implementation perspective, the paper emphasizes that existing steady-solver codes can be adapted through the addition of a single $H$-coupling matrix and FFT calls. This suggests that the method’s barrier to adoption lies less in re-engineering the finite-element infrastructure than in accommodating coupled temporal modes within an otherwise familiar incompressible-flow solver architecture.

## 5. Cardiovascular test cases and reported performance

The method is evaluated on three patient-specific physiological cases: Glenn pulmonary flow, cerebral arteries, and left main coronary arteries [2411.14315]. In all cases, the comparison baseline is a conventional time-marching solver.

| Case | Flow and discretization details | Reported HBNS result |
|---|---|---|
| Glenn pulmonary flow | $Re \approx 1{,}000$; mesh $\approx 1.95$ M elements; conventional time marching used $4{,}000$ steps over $4$ cycles | $N=19$ took $\approx 1.2$ h versus $\approx 50$ h; $\approx 40\times$ speed-up |
| Cerebral arteries | $Re \approx 120$; mesh $\approx 0.73$ M elements; inlet has ten significant Fourier components; conventional time marching used $4{,}000$ steps over $4$ cycles | $N=19$ took $\approx 0.4$ h versus $\approx 12$ h; $\approx 30\times$ faster |
| Left main coronary artery | $Re \approx 100$; mesh $\approx 1.34$ M elements; conventional time marching used $3{,}200$ steps | $N=25$ took $\approx 1.0$ h versus $\approx 30$ h; $\approx 30\times$ speed-up |

For the Glenn case, the inlet SVC flow is described as smooth, with energy cascading to higher harmonics. Using $N=7,13,19$, the integrated velocity error $E_{\Omega,T}^u$ is approximately $12\%$, $7\%$, and $3\%$, respectively, relative to time marching. The outlet LPA flow-rate error falls below $3\%$ at $N=7$ and below $1\%$ at $N=13$ [2411.14315].

For the cerebral artery case, the inlet contains ten significant Fourier components and is exact at $N=19$. The domain-time-integral velocity error is less than or equal to $2\%$ for $N \ge 19$ [2411.14315].

For the left main coronary case, the inflow profile contains sharp kinks and therefore demands $N \approx 25$ to achieve velocity error below $5\%$; by contrast, the steady case $N=1$ yields approximately $40\%$ error. The LAD outlet flow-rate error falls below $3\%$ at $N=13$ despite an inlet truncation of approximately $8\%$, which the paper attributes to viscous filtering. Memory overhead remains modest: HBNS at $N=25$ uses approximately $2.4\times$ the memory of the conventional solver because most memory is spent on mesh and sparsity structures rather than the unknown vector itself [2411.14315].

Across the three examples, the paper also states more generally that conventional time marching takes more than ten hours, whereas HBNS can produce a solution in approximately $30$ minutes, with up to two orders-of-magnitude cost reduction when enough modes are retained to represent the imposed boundary conditions accurately [2411.14315]. The test cases collectively indicate that the method’s accuracy is strongly conditioned by how faithfully the truncated harmonic basis represents the input waveform.

## 6. Scope, limitations, and prospective extensions

The paper identifies several advantages of HBNS: spectral accuracy in time, dramatic reduction in wall-clock time for periodic flows, built-in high-frequency noise filtering, and straightforward adaptation of existing steady-solver codes via a single $H$-coupling matrix and FFT calls [2411.14315]. These properties place HBNS within a class of methods that trade full transient resolution for a global-in-period representation of recurrent dynamics.

The same source also states the chief assumptions and current limits. HBNS requires strict time-periodicity at a known $\omega$ and is intended for Reynolds numbers below the laminar–turbulent transition, given as approximately $Re \lesssim 2{,}000$–$3{,}000$. Flows with geometry-driven instabilities or turbulence lie beyond its current remit. Resistance (Windkessel) outlet conditions have not yet been implemented [2411.14315]. Accordingly, HBNS should not be construed as a generic replacement for time marching in arbitrary unsteady flow; its applicability is tied to the existence of a stable periodic orbit that can be resolved with a finite harmonic truncation.

Potential extensions listed in the paper include fluid–structure interaction in the harmonic domain, application to respiratory or other physiologically periodic flows, coupling HBNS output to lumped-parameter models or machine-learning surrogates, and incorporation of richer stabilization, including two-parameter or variational-multiscale formulations, for high-Womersley-number regimes [2411.14315]. This suggests a broader methodological agenda in which harmonic balance is treated not only as an acceleration strategy for cardiovascular CFD, but also as a reusable space-time framework for periodic multiphysics and reduced-order coupling.

Source: https://www.emergentmind.com/topics/harmonic-balanced-navier-stokes-hbns