Harmonic Balance

Harmonic Balance equations describe periodic unsteady flows where the unsteady frequencies are known beforehand. They are suited to modeling periodically repeating flow fields that typically occur in turbomachinery such as compressors, turbines, and fans.

The Harmonic Balance method in Simcenter STAR-CCM+ is a full decomposition of the Navier-Stokes equations in the frequency domain. The unsteady, transient flow is represented in the frequency domain as a Fourier series in time. All transport equations for momentum, energy, and turbulence are decomposed into the frequency domain on the basis of fundamental driving modes, usually a blade-passing frequency or repeating wake modes. Steady-state equations representing the unsteady solution at discrete time levels in a single unsteady period are solved to obtain the Fourier coefficients.

The number of time levels required depends on the number of modes retained in the problem. The steady-state solution in every time level is implicitly coupled at the periodic boundaries by the physical time derivatives. The linear system is then subjected to approximate factorization to achieve implicit coupling between time levels.

Governing Equations

The governing equations for harmonic balance are the Navier-Stokes equations in integral form for a rigid, arbitrary control volume $V$ with differential surface area $d\mathbf{a}$ in a relative frame of reference rotating steadily with angular velocity $\mathbf{\Omega }$.

$$\underset{V}{\int }\frac{\partial \mathbf{W}}{\partial t}dV+\oint [\mathbf{F}-\mathbf{G}]\cdot d\mathbf{a}=\underset{V}{\int }\mathbf{H}dV$$

(5042)

where:

• $\mathbf{W}$ is the solution vector of conservation variables, $\mathbf{W}={[\rho ,\rho \mathbf{v},\rho E]}^T$.
• $\mathbf{F}$ is the inviscid flux, $\mathbf{F}={[\rho {\mathbf{v}}_r,\rho \mathbf{v}\otimes {\mathbf{v}}_r+p\mathbf{I},\rho E{\mathbf{v}}_r+p\mathbf{v}]}^T$.
• $\mathbf{G}$ is the viscous flux, $\mathbf{G}={[0,\mathbf{T},\mathbf{T}\cdot {\mathbf{v}}_r+\mathbf{q}]}^T$.
• $\mathbf{H}$ is the source vector, $\mathbf{H}={[S_u,\rho \mathbf{\Omega }\otimes \mathbf{v},S_u]}^T$.

and:

• $\rho$ is the density.
• $\mathbf{v}$ is the absolute velocity.
• ${\mathbf{v}}_r$ is the relative velocity, ${\mathbf{v}}_r=\mathbf{v}-\mathbf{r}\otimes \mathbf{\Omega }$.
• $E$ is the total enthalpy.
• $p$ is the pressure.
• $\mathbf{T}$ is the shear stress tensor.
• $\mathbf{q}$ is the heat flux.
• $S_u$ is a user-defined source term.

Harmonic Balance Equations

Since the solution $\mathbf{W}$ to the governing equations is periodic in time, it can be represented by a Fourier series:

$$\mathbf{W}(\mathbf{x},t)={\displaystyle \sum _{m=-M}^M{\hat{\mathbf{W}}}_m\left(\mathbf{x}\right)e^{i\omega mt}}$$

(5043)

where:

• $\omega$ is the fundamental frequency of the disturbance.
• $M$ is the number of harmonics retained in the solution.

${\hat{\mathbf{W}}}_m$ are the Fourier coefficients, uniquely determined from the discrete Fourier transform:

$${\hat{\mathbf{W}}}_m\left(\mathbf{x}\right)=\frac{1}{N}{\displaystyle \sum _{n=0}^{N-1}{\mathbf{W}}_n^{*}(\mathbf{x},t_n)e^{-i\omega m{t}_n}}$$

(5044)

where ${\mathbf{W}}_n^{*}$ are a set of $N=2M+1$ solutions at discrete time levels $t_n=nT/N$ distributed throughout one period of unsteadiness $T$.

At any location in the flowfield domain, the time-level solutions can be transformed into Fourier coefficients and vice versa using a discrete Fourier transform operator $E$ and its corresponding inverse $E^{-1}$ as follows:

$$\hat{\mathbf{W}}=E{\mathbf{W}}^{*}$$

(5045)

$${\mathbf{W}}^{*}=E^{-1}\hat{\mathbf{W}}$$

(5046)

where $E$ and $E^{-1}$ are square matrices of dimension $N\times N$ and the Fourier coefficients and time-level solutions have been assembled into the vectors $\hat{\mathbf{W}}$ and ${\mathbf{W}}^{*}$:

$$\begin{array}{c}\hat{\mathbf{W}}={[{\hat{\mathbf{W}}}_{-M,}{\hat{\mathbf{W}}}_{-M+1}\dots ,{\hat{\mathbf{W}}}_{M-1,}{\hat{\mathbf{W}}}_M]}^T\\ {\mathbf{W}}^{*}={[{\mathbf{W}}_0^{*},{\mathbf{W}}_1^{*},\dots ,{\mathbf{W}}_{N-2}^{*},{\mathbf{W}}_{N-1}^{*}]}^T\end{array}$$

The solutions at each discrete time level are obtained by applying the governing equations Eqn. (5042) to all the ${\mathbf{W}}^{*}$ simultaneously:

$$\underset{V}{\int }\frac{\partial {\mathbf{W}}^{*}}{\partial t}dV+\oint [{\mathbf{F}}^{*}-{\mathbf{G}}^{*}]\cdot d\mathbf{a}=\underset{V}{\int }{\mathbf{H}}^{*}dV$$

(5047)

where the flux and source vectors ${\mathbf{F}}^{*}$, ${\mathbf{G}}^{*}$, and ${\mathbf{H}}^{*}$ are evaluated using the corresponding time-level solution. Example:

$${\mathbf{F}}^{*}={\left[\mathbf{F}\right({\mathbf{W}}_0^{*}),\mathbf{F}({\mathbf{W}}_1^{*}),\mathrm{...},\mathbf{F}({\mathbf{W}}_{N-2}^{*}),\mathbf{F}({\mathbf{W}}_{N-1}^{*}\left)\right]}^T$$

The time derivative in Eqn. (5047) is evaluated by differentiating Eqn. (5047) with respect to time, then using Eqn. (5045) as follows:

$$\frac{\partial {\mathbf{W}}^{*}}{\partial t}=\frac{\partial E^{-1}}{\partial t}\hat{\mathbf{W}}=\frac{\partial E^{-1}}{\partial t}E{\mathbf{W}}^{*}$$

or

$$\frac{\partial {\mathbf{W}}^{*}}{\partial t}=D{\mathbf{W}}^{*}$$

(5048)

where $D$ is the pseudo-spectral $N\times N$ operator. Substituting Eqn. (5048) for the time derivative in Eqn. (5047) yields the final harmonic balance equation:

$$\underset{V}{\int }D{\mathbf{W}}^{*}dV+\oint [{\mathbf{F}}^{*}-{\mathbf{G}}^{*}]\cdot d\mathbf{a}=\underset{V}{\int }{\mathbf{H}}^{*}dV$$

(5049)