This post presents the equations for constant area mixing over a general meridonal surface, followed by a fixed-point iteration method for their solution. The control surface may be any surface of revolution, and mixed-out velocity may have components in the axial and radial directions. There are a lot of equations, but I think they might be useful for anyone else that needs to implement fully-general averaging of planes from a turbomachinery flow field.

Motivation

Computational fluid dynamics (CFD) simulations give the local flow field around a turbine or compressor geometry of interest. Sometimes, detail is useful for understanding the flow physics, but sometimes detail does not help with interpreting broader trends. In particular, we may want to compare a CFD solution to a lower-fidelity model or design tool that is faster to run but (we hope) captures the important features of the flow.

We can move forward with our comparison by averaging the flow from the CFD solution: converting a two-dimensional flow field over a cut plane into just one flow state. As noted by Cumpsty and Horlock (2006), we will neccesarily lose some information in this process, and have a choice to make regarding which quantities to conserve.

One choice of average is to ‘mix out’ the flow at constant area. (Sometimes, a constant static pressure assumption is used, but I don’t like this because it hides the lossyness of the mixing process as a change in geometry, not a change in the flow). We imagine the flow continuing down an infinite frictionless duct until it reaches a uniform state by natural mixing. The end state will conserve mass, moment of momentum, and energy but will have a reduced pressure and increased entropy because mixing is irreversible.

Cumpsty and Horlock note that mixed-out averaging is pessimistic, because it assumes that the mixing process is wholly irreversible. In theory, a non-uniformity passing through a downstream blade row could be brought to a uniform state by a reversible process. I don’t think this is a problem because I don’t believe that reversible mixing is significant in practice. Anyway, if our steady CFD uses mixing planes between blade rows we have already assumed that circumferential non-uniformities mix out entirely irreversibly.

Burdett and Povey (2022), in their utterly comprehensive paper, argue that mixing out to a single state is not physically rational in swirling flow. Even in an infinite annular duct, the flow will never reach uniformity because radial equilibrium requires a spanwise static pressure gradient. This is a fair point, but there is nothing we can do about it if we want an averaged value to compare to a low-order model.

My opinion is that constant-area mixing is the most sensible way to average a flow field. It preserves the quantities we care about: conserving angular momentum means that integrated blade loadings are matched; and conserving energy means that work is also matched. It does not preserve entropy change, but including mixing losses in efficiences is justifiable.

Governing equations

The governing equations that define the mixed-out flow are, in three-dimensional polar coordinates:

$$ \newcommand{\rovx}{\rho V_x} \newcommand{\mix}[1]{\tilde{#1}} \newcommand{\m}{\dot{m}} \newcommand{\Pr}{\dot{P}_r} \newcommand{\Px}{\dot{P}_x} \newcommand{\Pt}{\dot{P}_{r\theta}} \newcommand{\vx}{V_x} \newcommand{\vr}{V_r} \newcommand{\vt}{V_\theta} \newcommand{\vmagm}{\mix{V}} \newcommand{\vmm}{\mix{V_m}} \newcommand{\pm}{\mix{p}} \newcommand{\vrm}{\mix{V_r}} \newcommand{\vxm}{\mix{V_x}} \newcommand{\vtm}{\mix{V_\theta}} \newcommand{\rom}{\mix{\rho}} \newcommand{\ho}{h_0} \newcommand{\hstatm}{\mix{h}} \newcommand{\hom}{\mix{h_0}} \newcommand{\Ax}{A_x} \newcommand{\Ar}{A_r} \newcommand{\I}{{I}} \newcommand{\Itot}{\dot{\mathcal{I}}} \newcommand{\dee}[1]{\,\mathrm{d}{#1}} \newcommand{\Im}{\mix{\I}} \newcommand{\rovxsq}{\rho V_x^2} \newcommand{\rovrsq}{\rho V_r^2} \newcommand{\rovr}{\rho V_r} \newcommand{\Bm}{\mix{\beta}} \newcommand{\tanBm}{\tan\Bm} \newcommand{\tansqBm}{{\tan^2\Bm}} \newcommand{\iter}[2]{{{#1}^{(#2)}}} \newcommand{\Asx}{A_x^*} \newcommand{\Asr}{A_r^*} \begin{align} &\iint \rovx \dee{\Ax} + \iint \rovr \dee{\Ar} &&= \rom \vxm \Ax + \rom \vrm \Ar &&= \m\,,\\ & \iint \left(\rovxsq + p\right) \dee{\Ax} + \iint \rovr \vx \dee{\Ar} &&= \left(\rom \vxm^2 + \pm \right) \Ax + \rom \vrm \Ar &&= \Px\,,\\ & \iint \rovx \vr \dee{\Ax} + \iint \left(\rovrsq + p \right) \dee{\Ar} &&= \rom \vxm \vrm \Ax + \left( \rom \vrm^2 + \pm \right) \Ar &&= \Pr\,,\\ & \iint \rovx r \vt \dee{\Ax} + \iint \rovr r \vt \dee{\Ar} &&= \rom \vxm \mix{r} \vtm \Ax + \rom \vrm \mix{r} \vtm \Ar &&=\Pt\,, \\ & \iint \rovx \I \dee{\Ax} + \iint \rovr \I \dee{\Ar} &&= \rom \vxm \Im\Ax + \rom \vrm \Im \Ar &&= \Itot\,, \end{align} $$

where the rothalpy $$ \I = \ho - \Omega r \vt = h + \frac{1}{2}V^2 - \Omega r \vt \, , $$ and \(\mix{r}\) is the radius at which the mixed-out state is specified. Here, we shall set \(\mix{r}\) to the mid-span value. Each of these equations embodies conservation of a physical quantity. The first column is a total flow rate given by integration of a non-uniform flux over the surface. The second column is the corresponding flow rate for the mixed out state, denoted by \(\tilde{\Box}\). The third column is a scalar variable that we will use as a shorthand for the total flow rate of each quantity:

• mass flow, \(\m\,\);
• axial momentum flow, \(\Px\,\);
• radial momentum flow, \(\Pr\,\);
• moment of momentum flow, \(\Pt\) (in a slight abuse of notation);
• rothalpy flow, \(\Itot\,\).

Solving the system

Various numerical intergration methods are available for computing the flux integrals. In my case, I have a CFD solution at each node of a rectangular grid. Given that the solution is only second-order accurate, a simple midpoint rule is sufficient. The total flow rate is: the average of four nodal flux values for each face, dotted with the face area vector, summed over all faces.

We seek a solution for the mixed out state given total flow rates from the system of five non-linear equations. We also need an equation of state for the fluid, to relate changes in enthalpy to density and pressure. Fixed point iteration is robust and easy to code. Computational expense is not much of a concern here; taking the flux integrals dominates the runtime.

[Update 2023-07-29: The fixed point iteration I describe below will only converge to a subsonic mixed-out state, even when provided with a supersonic initial guess! This is a consequence of the existence of two physically plausible solutions satisfying the conservation laws. I expect that it would be possible to reformulate the fixed-point iteration in a way that would converge to the supersonic solution, but I did not go down that route. Instead, I use the iteration step below, but apply the Nelder-Mead optimiser to update \(\rom\) and \(\Bm\) to drive the residuals of mass, momentum, and energy conservation below a tolerance.]

Any of the \(A\) or \(V\) components can in general be zero, so we need to branch on the relative values of \(\Ax\) and \(\Ar\) to avoid numerical issues.

Mostly axial \(\Ax > \Ar\)

If \(\Ax > \Ar\) it is safe to assume \(\Ax \ne 0\), \(\vxm \ne 0\) and finite \(\tanBm\). The flow is predominantly in the axial direction, however, \(\Ar\) and \(\vrm\) may be zero. We can neatly eliminate \(\vrm\) from most of the conservation equations in favour of \(\vxm\) and \(\tanBm=\vrm/\vxm\):

$$ \begin{align} \m &= \rom \vxm \left( \Ax + \Ar \tanBm \right) &&= \rom \vxm \Asx \\ \Px &= \pm \Ax + \rom \vxm^2 \left( \Ax + \Ar \tanBm \right) &&= \pm \Ax + \rom \vxm^2 \Asx\\ \Pr &= \pm \Ar + \rom \vrm^2 \left( \frac{\Ax}{\tanBm} + \Ar \right) &&= \pm \Ar + \rom \vrm^2 \frac{\Asx}{\tanBm}\\ \Pt &= \rom \vxm \mix{r} \vtm \left( \Ax + \Ar\tanBm \right) &&= \rom \vxm \mix{r} \vtm \Asx\\ \Itot &= \rom \vxm \Im \left( \Ax + \Ar \tanBm \right) &&= \rom \vxm \Im \Asx \end{align} $$ where $$ \Asx = \Ax + \Ar \tanBm $$ Here, \(\Asx\) is a reference area defined such that the total flow due to both axial and radial components is given by a mixed-out axial flux times the reference area.

We can take an arithmetic average of nodal density and pitch angle to get good initial guesses for the mixed out flow, \(\iter{\rom}{0}\) and \(\iter{\Bm}{0}\). We need to make sure that \(\iter{\Bm}{0}\ne 0\), otherwise \(\vrm\) will get stuck at zero and the iterations will not converge. Proceed as follows:

1. Use conservation of mass to set the axial velocity $$ \iter{\vxm}{n} = \frac{\m}{\iter{\rom}{n}\Asx}\,. $$
2. Given an axial velocity, we can use conservation of axial momentum to get the mixed-out pressure, and conservation of moment of momentum to calculate tangential velocity: $$ \iter{\pm}{n} = \frac{\Px - \iter{\rom}{n} \iter{\vxm^2}{n} \Asx}{\Ax}\,, \quad \iter{\vtm}{n} = \frac{\Pt}{\iter{\mix{r} \rom}{n} \iter{\vxm}{n} \Asx}\,, \quad $$
3. Mixed-out pressure then yields radial velocity squared using conservation of radial momentum: $$ \iter{\vrm^2}{n} = \frac{\Pr - \iter{\pm}{n} \Ar}{\iter{\rom}{n} \Asx}\tanBm \,. $$
4. Now we have all velocity components, calculate meridonal velocity and velocity magnitude: $$ \iter{\vmm^2}{n} = \iter{\vxm^2}{n} + \iter{\vrm^2}{n} \,, \quad \iter{\vmagm^2}{n} = \iter{\vxm^2}{n} + \iter{\vrm^2}{n} + \iter{\vtm^2}{n} $$
5. Use conservation of rothalpy and velocity magnitude to get static enthalpy: $$ \iter{\Im}{n} = \frac{\Itot}{\iter{\rom}{n} \iter{\vxm}{n} \Asx}\,, \quad \iter{\hstatm}{n} = \iter{\Im}{n} + \mix{r} \iter{\vtm}{n} \Omega - \frac{1}{2}\iter{\vmagm^2}{n} $$
6. Finally, we are in a position to update our initial guesses. Evaluate the equation of state for density; use trigonometry for pitch angle: $$ \iter{\rom}{n+1} = \rho\left(\iter{\pm}{n}, \iter{\hstatm}{n}\right)\,, \quad \iter{\Bm}{n+1} = \arccos \frac{\iter{\vxm}{n}}{\iter{\vmm}{n}}\,. $$
7. Finish if the changes are smaller than a tolerance, otherwise return to 1.

Mostly radial \(\Ar > \Ax\)

In the case where the flow is predominantly axial, \(\Ar > \Ax\). Now, \(\Ar\ne 0\), \(\vrm\ne 0\) and \(\tanBm \ne 0\). It is possible that \(\Ax=0\) or \(\vxm=0\), so we cannot divide by them as in the algorithm above. Instead, we express the conservation equations using a radial reference area $$ \begin{align} \m &= \rom \vrm \left( \frac{\Ax}{\tanBm} + \Ar \right) &&= \rom \vrm \Asr \\ \Px &= \pm \Ax + \rom \vxm^2 \left( \Ax+ \Ar\tanBm \right) &&= \pm \Ax + \rom \vxm^2 \Asr\tanBm\\ \Pr &= \pm \Ar + \rom \vrm^2 \left( \frac{\Ax}{\tanBm} + \Ar \right) &&= \pm \Ar + \rom \vrm^2 \Asr\\ \Pt &= \rom \vxm \mix{r} \vtm \left( \frac{\Ax}{\tanBm} + \Ar\right) &&= \rom \vrm \mix{r} \vtm \Asr\\ \Itot &= \rom \vrm \Im \left( \frac{\Ax}{\tanBm} + \Ar \right) &&= \rom \vrm \Im \Asr \end{align} $$ where $$ \Asr = \frac{\Ax}{\tanBm} + \Ar = \frac{\Asx}{\tanBm}\,. $$ Analogous to the previous case, the total flow is given by a radial flux times the radial reference area. With the same initial guesses for \(\iter{\rom}{0}\) and \(\iter{\Bm}{0}\), the solution procedure is then:

1. Use conservation of mass to set the radial velocity $$ \iter{\vrm}{n} = \frac{\m}{\iter{\rom}{n}\Asr}\,. $$
2. Given a radial velocity, we can use conservation of radial momentum to get the mixed-out pressure, and conservation of moment of momentum to calculate tangential velocity: $$ \iter{\pm}{n} = \frac{\Pr - \iter{\rom}{n} \iter{\vrm^2}{n} \Asr}{\Ar}\,, \quad \iter{\vtm}{n} = \frac{\Pt}{\iter{\mix{r} \rom}{n} \iter{\vrm}{n} \Asr}\,, \quad $$
3. Mixed-out pressure then yields axial velocity squared using conservation of axial momentum: $$ \iter{\vxm^2}{n} = \frac{\Px - \iter{\pm}{n} \Ax}{\iter{\rom}{n} \Asr \tanBm} \,. $$
4. Now we have all velocity components, calculate meridonal velocity and velocity magnitude: $$ \iter{\vmm^2}{n} = \iter{\vxm^2}{n} + \iter{\vrm^2}{n} \,, \quad \iter{\vmagm^2}{n} = \iter{\vxm^2}{n} + \iter{\vrm^2}{n} + \iter{\vtm^2}{n} $$
5. Use conservation of rothalpy and velocity magnitude to get static enthalpy: $$ \iter{\Im}{n} = \frac{\Itot}{\iter{\rom}{n} \iter{\vrm}{n} \Asr}\,, \quad \iter{\hstatm}{n} = \iter{\Im}{n} + \mix{r} \iter{\vtm}{n} \Omega - \frac{1}{2}\iter{\vmagm^2}{n} $$
6. Update the initial guesses, using equation of state and velocity components: $$ \iter{\rom}{n+1} = \rho\left(\iter{\pm}{n}, \iter{\hstatm}{n}\right)\,, \quad \iter{\Bm}{n+1} = \arcsin \frac{\iter{\vrm}{n}}{\iter{\vmm}{n}}\,. $$
7. Finish if the changes are smaller than a tolerance, otherwise return to 1.

An aside on multi-row CFD

Suppose we have a completely axial flow with \(\vr=0\) sampled on a meridional plane. If the plane slopes at all, i.e. \(x\) is not constant, then \(\Ar\ne0\,\). Conservation of radial momentum becomes:

$$ \iint p \dee{\Ar} = \rom \vxm \vrm \Ax + \left( \rom \vrm^2 + \pm \right) \Ar $$

If the flow on this sloping plane is non-uniform, then in general \(p\ne\pm\) and we will need a non-zero mixed-out raidial velocity \(\vrm\ne 0\) to keep the equation balanced! Another way of thinking about this is that mixing is equivalent to applying an opposing drag force against the flow; drag is an irreversible process that generates entropy (Denton 1993). If the plane slopes in the meridional view, then this force will have a component in the radial direction.

As pointed out to me by Chris Clark, when using mixing planes to connect adjacent blade rows in a steady CFD calculation, the shape of the mixing plane is a modelling choice that will affect the results. You probably want your mixing plane to be orthogonal to the annulus. But if the blades are swept, then the mixing plane will not be located the same distance away up the span.

Conclusion

In this post, we have discussed motivation for mixed-out averaging, then looked at the governing equations and numerical algorithms for solving them. For a robust solution, two different fixed-point iteration methods are needed, depending on the flow direction. The equations also show that the shape of the mixing surface can introduce velocity components that were not present in the non-uniform flow.

References

• Cumpsty, N. A., and Horlock, J. H. (2005). “Averaging Nonuniform Flow for a Purpose.” J. Turbomach. 128(1). doi:10.1115/1.2098807.
• Burdett, D., and Povey, T. (2022). “Analysis of Averaging Methods for Nonuniform Total Pressure Fields.” J. Turbomach. 144(5). doi:10.1115/1.4053020.