Harmonic-Balanced Navier–Stokes (HBNS)

• Harmonic-Balanced Navier–Stokes (HBNS) is a frequency-domain finite-element method that leverages periodicity to transform transient flow problems into spectral representations.
• It employs a harmonic decomposition with FFT/IFFT operations to reduce the computational cost by up to two orders of magnitude compared to conventional solvers.
• The method uses a stabilized Galerkin/least-squares formulation, maintaining high accuracy in patient-specific cardiovascular simulations for laminar flows.

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" (Jia et al., 2024), 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 (Jia et al., 2024).

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]$ (Jia et al., 2024). 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 $T$0 and Reynolds numbers below the laminar–turbulent transition (Jia et al., 2024).

2. Harmonic decomposition and frequency-domain reformulation

HBNS truncates the temporal field to $T$1 modes, with $T$2 time samples, and expands the solution in Fourier modes:

$T$3

The boundary data $T$4 and $T$5 are expanded in the same manner so as to interpolate their spectral content exactly up to mode $T$6 (Jia et al., 2024). Substituting these expansions into the governing equations and matching terms in $T$7 yields, for each harmonic index $T$8,

$T$9

$\omega = 2\pi/T$0

In stacked form, over the coefficient vectors $\omega = 2\pi/T$1 and $\omega = 2\pi/T$2, this becomes

$\omega = 2\pi/T$3

where $\omega = 2\pi/T$4 and $\omega = 2\pi/T$5 is the convolution matrix $\omega = 2\pi/T$6 (Jia et al., 2024).

A key computational observation is that direct assembly in harmonic space is $\omega = 2\pi/T$7. HBNS therefore applies the inverse discrete Fourier transform, $\omega = 2\pi/T$8, to obtain a harmonic-balance form at $\omega = 2\pi/T$9 real time points $$\rho \frac{\partial u}{\partial t} + \rho (u \cdot \nabla)u + \nabla p - \mu \Delta u = 0 \quad \text{in } \Omega \times (0,T],$$0:

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

where $$\rho \frac{\partial u}{\partial t} + \rho (u \cdot \nabla)u + \nabla p - \mu \Delta u = 0 \quad \text{in } \Omega \times (0,T],$$2 is a skew-symmetric time-coupling matrix and $$\rho \frac{\partial u}{\partial t} + \rho (u \cdot \nabla)u + \nabla p - \mu \Delta u = 0 \quad \text{in } \Omega \times (0,T],$$3 (Jia et al., 2024). In this representation, convection and diffusion act pointwise in time, whereas the temporal coupling is encoded by $$\rho \frac{\partial u}{\partial t} + \rho (u \cdot \nabla)u + \nabla p - \mu \Delta u = 0 \quad \text{in } \Omega \times (0,T],$$4. Because $$\rho \frac{\partial u}{\partial t} + \rho (u \cdot \nabla)u + \nabla p - \mu \Delta u = 0 \quad \text{in } \Omega \times (0,T],$$5 is applied through FFT/IFFT operations, the time-coupling cost is reduced to $$\rho \frac{\partial u}{\partial t} + \rho (u \cdot \nabla)u + \nabla p - \mu \Delta u = 0 \quad \text{in } \Omega \times (0,T],$$6.

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

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

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

The Galerkin/least-squares form seeks $$\rho \frac{\partial u}{\partial t} + \rho (u \cdot \nabla)u + \nabla p - \mu \Delta u = 0 \quad \text{in } \Omega \times (0,T],$$9 such that, for all $$\nabla \cdot u = 0 \quad \text{in } \Omega \times (0,T],$$0,

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

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

where the momentum-operator residual is

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

and $$\nabla \cdot u = 0 \quad \text{in } \Omega \times (0,T],$$4 is a diagonal stabilization parameter, one value per $$\nabla \cdot u = 0 \quad \text{in } \Omega \times (0,T],$$5, defined by

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

with $$\nabla \cdot u = 0 \quad \text{in } \Omega \times (0,T],$$7 the contravariant metric, $$\nabla \cdot u = 0 \quad \text{in } \Omega \times (0,T],$$8, and $$\nabla \cdot u = 0 \quad \text{in } \Omega \times (0,T],$$9 a shape constant, approximately $u=g$0 for tetrahedra (Jia et al., 2024).

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 $u=g$1 (Jia et al., 2024).

Spatial discretization uses equal-order linear tetrahedral interpolation $u=g$2 for all components:

$u=g$3

Dirichlet conditions are enforced by directly prescribing nodal values at the $u=g$4 time points. The resulting global system has $u=g$5 unknowns per node (Jia et al., 2024). 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 $u=g$6 is split as

$u=g$7

where $u=g$8 has the same sparsity pattern as a steady-state Navier–Stokes Jacobian and is block-diagonal in time, with cost $u=g$9, whereas $\Gamma_g \times (0,T]$0 encodes the $\Gamma_g \times (0,T]$1-coupling through

$\Gamma_g \times (0,T]$2

The $\Gamma_g \times (0,T]$3 blocks are dense in time but are applied via FFT/IFFT in $\Gamma_g \times (0,T]$4 (Jia et al., 2024). The paper therefore states that overall assembly-and-solve cost scales nearly $\Gamma_g \times (0,T]$5.

The nonlinear algebraic system is solved by Newton–Raphson. At iteration $\Gamma_g \times (0,T]$6,

$\Gamma_g \times (0,T]$7

with residual vector $\Gamma_g \times (0,T]$8 and tangent $\Gamma_g \times (0,T]$9 (Jia et al., 2024). To enhance convergence, a pseudo-time derivative

$(-pn+\mu \partial u/\partial n)=h\,n$0

is added and treated with a generalized-$(-pn+\mu \partial u/\partial n)=h\,n$1 scheme. The pseudo-time step $(-pn+\mu \partial u/\partial n)=h\,n$2 is chosen so that the convective CFL is $(-pn+\mu \partial u/\partial n)=h\,n$3. Each linear system is solved with GMRES using a Jacobi preconditioner on $(-pn+\mu \partial u/\partial n)=h\,n$4, and iterations continue until $(-pn+\mu \partial u/\partial n)=h\,n$5 drops by $(-pn+\mu \partial u/\partial n)=h\,n$6 (Jia et al., 2024).

From an implementation perspective, the paper emphasizes that existing steady-solver codes can be adapted through the addition of a single $(-pn+\mu \partial u/\partial n)=h\,n$7-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 (Jia et al., 2024). In all cases, the comparison baseline is a conventional time-marching solver.

Case	Flow and discretization details	Reported HBNS result
Glenn pulmonary flow	$(-pn+\mu \partial u/\partial n)=h\,n$8; mesh $(-pn+\mu \partial u/\partial n)=h\,n$9 M elements; conventional time marching used $\Gamma_h \times (0,T]$0 steps over $\Gamma_h \times (0,T]$1 cycles	$\Gamma_h \times (0,T]$2 took $\Gamma_h \times (0,T]$3 h versus $\Gamma_h \times (0,T]$4 h; $\Gamma_h \times (0,T]$5 speed-up
Cerebral arteries	$\Gamma_h \times (0,T]$6; mesh $\Gamma_h \times (0,T]$7 M elements; inlet has ten significant Fourier components; conventional time marching used $\Gamma_h \times (0,T]$8 steps over $\Gamma_h \times (0,T]$9 cycles	$u(x,t+T)=u(x,t),$0 took $u(x,t+T)=u(x,t),$1 h versus $u(x,t+T)=u(x,t),$2 h; $u(x,t+T)=u(x,t),$3 faster
Left main coronary artery	$u(x,t+T)=u(x,t),$4; mesh $u(x,t+T)=u(x,t),$5 M elements; conventional time marching used $u(x,t+T)=u(x,t),$6 steps	$u(x,t+T)=u(x,t),$7 took $u(x,t+T)=u(x,t),$8 h versus $u(x,t+T)=u(x,t),$9 h; $T$00 speed-up

For the Glenn case, the inlet SVC flow is described as smooth, with energy cascading to higher harmonics. Using $T$01, the integrated velocity error $T$02 is approximately $T$03, $T$04, and $T$05, respectively, relative to time marching. The outlet LPA flow-rate error falls below $T$06 at $T$07 and below $T$08 at $T$09 (Jia et al., 2024).

For the cerebral artery case, the inlet contains ten significant Fourier components and is exact at $T$10. The domain-time-integral velocity error is less than or equal to $T$11 for $T$12 (Jia et al., 2024).

For the left main coronary case, the inflow profile contains sharp kinks and therefore demands $T$13 to achieve velocity error below $T$14; by contrast, the steady case $T$15 yields approximately $T$16 error. The LAD outlet flow-rate error falls below $T$17 at $T$18 despite an inlet truncation of approximately $T$19, which the paper attributes to viscous filtering. Memory overhead remains modest: HBNS at $T$20 uses approximately $T$21 the memory of the conventional solver because most memory is spent on mesh and sparsity structures rather than the unknown vector itself (Jia et al., 2024).

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 $T$22 minutes, with up to two orders-of-magnitude cost reduction when enough modes are retained to represent the imposed boundary conditions accurately (Jia et al., 2024). 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 $T$23-coupling matrix and FFT calls (Jia et al., 2024). 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 $T$24 and is intended for Reynolds numbers below the laminar–turbulent transition, given as approximately $T$25–$T$26. Flows with geometry-driven instabilities or turbulence lie beyond its current remit. Resistance (Windkessel) outlet conditions have not yet been implemented (Jia et al., 2024). 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 (Jia et al., 2024). 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.