# An extension of the compound flow theory with shear and friction

Jan Van den Berghe1,2<EMAIL_ADDRESS>Miguel A. Mendez1  Yann
Bartosiewicz2  1 von Karman Institute for Fluid Dynamics, Waterloosesteenweg
72, 1640 Sint-Genesius-Rode, Belgium 2 Institute of Mechanics, Materials, and
Civil Engineering (iMMC), Université catholique de Louvain (UCLouvain), 1348
Louvain-la-Neuve, Belgium

###### Abstract

Compound flows consist of two or more parallel compressible streams inside a
duct. Its theoretical treatment has recently been applied to model ejectors,
but also to post-process and to interpret CFD results. More precisely, it has
been observed that compound flows can choke upstream of the geometrical
throat. While it is widely known that friction can push the sonic section
downstream, no mechanism has been identified to explain the movement in the
opposite direction.

In this work, we extend the existing compound flow theory with shear and
friction. The resulting model captures the up- and downstream movement of the
sonic section. We show analytically that shear between the streams is the
driving force behind the upstream movement by analysing the singularity of the
governing equations in the sonic section. Additionally, we provide a quasi-
third order approximation of the governing equations in the vicinity of the
sonic section to alleviate potential numerical issues. Finally, we compare the
predictions of the model with post-processed CFD results of a compound nozzle.
We study the individual effects of shear and friction through an inviscid
simulation for the isentropic case and a viscid simulation with a slip or no-
slip condition at the wall to include shear with or without friction.

The proposed extension offers new possibilities in modelling and analysis of
compound flows.

###### keywords:

Gas dynamics, choking, compressible parallel streams.

##  1 Introduction

Parallel compressible streams are at the heart of many applications, the most
classic example being jet engines with a high-speed core and a slower
surrounding jet. We will refer to such heterogeneous flows as a single
compound flow. In this paper, we are interested in internal compound flows,
which appear for example in ejectors or (sc)ramjets.

Recently, Lamberts et al. (2018) observed that compound flows can choke
_upstream_ of the geometrical throat, in contrast to the classic _downstream_
movement of the sonic section due to friction. Kracik & Dvorak (2023) showed
that friction has this effect also on compound flows, by numerically
increasing the wall roughness. Lamberts et al. (2018) attributed the shift
upstream to non-isentropic effects, but the inherent cause remained unknown.
Both authors above have tackled this question through 2/3D CFD, which has
allowed to observe the phenomenon, but which lacks the interpretability of
theoretical models.

The initial theoretical treatment of compound flows was proposed by Pearson et
al. (1958). The main idea consists of solving the 1D conservation equations
across each stream of the compound flow, complemented with the assumption that
the static pressure is uniform in each cross-section. This assumption is
required to obtain a unique solution. A priori, the conservation equations can
be fulfilled with any dividing line between the streams, as if we study the
flow in independent channels. The axial pressure profile in the streams will
depend on the evolution of the respective cross-sections, but mass, momentum
and energy will always be conserved. However, only one configuration of the
cross-sections will keep the static pressure equal in each cross-section.
Pearson et al. (1958) state that the assumption of equal pressure is
reasonable if the curvature of the streamlines remains limited. The major
result of their study lies in the fact that a compound flow can be choked with
a combination of subsonic and supersonic streams. They show that information
can not travel upstream in the subsonic components due to the assumption of
uniform static pressure.

The theory of Pearson et al. (1958) has been adopted and derived in different
ways by Hoge & Segars (1965), Bernstein et al. (1967), Agnone (1971) and Fage
(1976), through the singularity of the equations or through minimising the
impulse function (as defined in chapter 4 of Shapiro (1953)). Bernstein et al.
(1967) gave the name ‘compound flow theory’, which we will also adopt. These
works did not consider exchanges between the streams, hence the flow remained
isentropic.

Later work investigated the effect of forces on the compound flow. In
particular, Clark (1995) and Papamoschou (1996) used the correlation proposed
by Papamoschou & Roshko (1988) to calculate the shear stress between the
streams. The authors used the resulting compound flow equations _with_
momentum exchange in a 1D ejector model. Surprisingly, neither of the authors
investigated choking, which is nonetheless a crucial element in ejectors. An
earlier work by Otis (1976) treats mixing and choking of a compound flow, but
is limited to constant-area ducts. An important conclusion of that work is
that additional effects such as forces or heat exchange do not affect the
choking condition as derived in the original compound theory (i.e., the
equivalent of a unitary Mach number, cf. sections 2 and 2.3), because they
appear on the right hand side of the conservation equations. To the best of
our knowledge, a general treatment of compound choking with momentum exchange
in variable area ducts has not yet been proposed in literature.

In this work, we extend the compound flow theory with shear and friction and
investigate their effect on choking. We show analytically that the choked
section can indeed move upstream or downstream depending on the relative
magnitude of the shear and wall friction forces. Our analysis is largely
similar to so-called ‘generalised 1D flow’, which describes quasi-1D
compressible flow with any combination of area change, friction, heat exchange
and/or mass injection. Equations for generalised 1D flow of a _single_ stream
can be found in classic textbooks, see (Shapiro, 1953, chapter 8). We propose
a similar generalised treatment of _compound_ flows by extending the original
theory with momentum exchange and by formally analysing choking of such a
‘generalised compound flow’. Our contribution is both theoretical and
numerical: (1) we propose a simple theoretical model that captures the
position of the sonic section and (2) we propose a quasi-third order
approximation of the pressure gradient in the vicinity of the sonic section to
evade numerical issues near the singular sonic point.

The governing equations are derived in section 2 and converted into an
explicit system of Ordinary Differential Equations (ODEs), which are more
interpretable and more practical for numerical integration. Section 3 goes in
more detail on the singularity of the equation for the static pressure at the
sonic section. We present numerical test cases of increasing complexity in
section 4 and compare the results with the 1D predictions in section 5.
Section 6 closes the paper with conclusions and perspectives.

##  2 Model definition

In this section, we first derive the governing equations in their general form
without specifying the forces at play. This facilitates future developments.
Next, we define the shear and friction forces and simplify the resulting
expression for the static pressure gradient, which is singular in the sonic
section. The corresponding choking mechanism is analysed in the final
subsection.

###  2.1 Governing equations

We consider an internal flow consisting of two parallel streams, which we will
denote with indices $p$ and $s$ (cf. figure 1). The theory can be derived for
any number of streams, but this technicality does not alter the final
conclusions of this work. Moreover, most applications consist of two streams.
Therefore, and for the sake of simplicity, only two streams are considered in
the remainder of this work. We assume the fluid is identical in all streams
and behaves like an ideal gas. The main assumption for compound-compressible
flows states that the static pressure is uniform in all cross-sections, hence:

$p_{p}(x)=p_{s}(x)=p(x).$ (1)

where $p$ denotes the static pressure and $x$ is the axial coordinate. The
evolution of the individual cross-sections $A_{i}$ is unknown a priori and
should be calculated. They are subject to a geometrical constraint imposed by
the known contours of the channel $A(x)$:

$\sum_{i\in[p,s]}A_{i}(x)=A(x).$ (2)

The governing equations consist of conservation of mass, momentum and energy
in each stream $i\in[p,s]$. For a steady, quasi-1D flow, these read as
follows:

$\displaystyle\frac{d}{dx}(\rho_{i}A_{i}V_{i})$ $\displaystyle=0,$ (3)
$\displaystyle\frac{d}{dx}(\rho_{i}A_{i}V_{i}^{2})$
$\displaystyle=-A_{i}\frac{dp_{i}}{dx}+F_{i},$ (4)
$\displaystyle\frac{d}{dx}(\rho_{i}A_{i}V_{i}h_{t,i})$ $\displaystyle=0,$ (5)

where $\rho$ denotes the density, $V$ denotes the axial velocity, $F$ is the
force per unit length exerted on stream $i$ (other than the pressure) and
$h_{t}$ denotes the total enthalpy. This means that the streams do not
exchange mass and that the separating lines between the streams hold the
meaning of dividing _streamlines_. In this work, we do not account for heat
exchange, but focus on momentum exchange instead.
This system of equations can be converted to an explicit system of ODEs to
facilitate numerical integration. The derivation can be found in classic text
books, for example in chapter 8 of Shapiro (1953). We cast the equations in
terms of the static pressure $p$, the total pressure $p_{t}$ and the total
temperature $T_{t}$:

$\displaystyle\frac{1}{p_{i}}\frac{dp_{i}}{dx}$
$\displaystyle=\left[\frac{1+\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}}{1-M_{i}^{2}}\right]\frac{F_{i}}{A_{i}p_{i}}+\left[\frac{\gamma\mbox{{Ma}}_{i}^{2}}{1-\mbox{{Ma}}_{i}^{2}}\right]\frac{1}{A_{i}}\frac{dA_{i}}{dx},$
(6) $\displaystyle\frac{1}{p_{t,i}}\frac{dp_{t,i}}{dx}$
$\displaystyle=\frac{F_{i}}{A_{i}p_{i}},$ (7)
$\displaystyle\frac{1}{T_{t,i}}\frac{dT_{t,i}}{dx}$ $\displaystyle=0,$ (8)

where we introduced the Mach number Ma and the ratio of specific heat
capacities $\gamma$. These equations are practical because the static pressure
appears in the main assumption (cf. equation (1)) and will play a crucial role
in the analysis of choked conditions (cf. section 2.3). The equations for the
total pressure and total temperature are relatively simple compared to other
variables and, moreover, are never singular.
The equations (6)-(8) do not form a closed system of equations since the
evolution of the cross-sections $A_{i}$ is unknown (as opposed to a single
stream where $A_{i}=A$). The info we _do_ have is the profile of the channel
$A(x)$, e.g., a nozzle, and hence also its derivative in space. Therefore, we
can eliminate $\frac{dA_{i}}{dx}$ from the system by summing over the streams
$i$. We first invert equation (6) to isolate the gradient of the cross-section
$A_{i}$:

$\frac{dA_{i}}{dx}=\left[A_{i}\frac{1-\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}\right]\frac{1}{p}\frac{dp}{dx}-\left[\frac{1+\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}\right]\frac{F_{i}}{p},$
(9)

where we have omitted the index for the pressure, cf. equation (1). Note that
this equation can be integrated to compute the cross-section of a single
stream if the pressure gradient is known. Next, we sum over the streams $i$:

$\sum_{i\in[p,s]}\frac{dA_{i}}{dx}=\frac{dA}{dx}=\sum_{i\in[p,s]}\left[A_{i}\frac{1-\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}\right]\frac{1}{p}\frac{dp}{dx}-\sum_{i\in[p,s]}\left[\frac{1+\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}\right]\frac{F_{i}}{p}.$
(10)

The pressure gradient is the only remaining unknown variable and we obtain an
explicit equation by inverting again:

$\frac{1}{p}\frac{dp}{dx}=\dfrac{1}{\beta}\left(\dfrac{dA}{dx}+\sum\limits_{i\in[p,s]}\left[\dfrac{1+\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}\right]\dfrac{F_{i}}{p}\right)$
(11)

where $\beta$ denotes the compound choking indicator:

$\beta=\sum_{i\in[p,s]}A_{i}\dfrac{1-\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}$
(12)

As discussed in literature, e.g., by Bernstein et al. (1967), $\beta>0$
corresponds to a compound-subsonic flow (in absence of forces, the pressure
gradient and the area gradient share the same sign, e.g., leading to an
expansion in a convergent channel), $\beta=0$ corresponds to compound-sonic
flow and $\beta<0$ to compound-supersonic flow. A convenient variable is the
so-called equivalent Mach number $\mbox{{Ma}}_{eq}$, introduced by Hedges &
Hill (1974):

$\mbox{{Ma}}_{eq}=\left(\gamma\frac{\beta}{A}+1\right)^{-\frac{1}{2}},$ (13)

which equals one for a compound choked flow. Section 2.3 analyses choking in
more details. It is important to note that the forces do not affect the
denominator of equation (11), as pointed out by Otis (1976). The compound
choking criterion of $\beta=0$ thus remains unchanged with respect to the
isentropic case.

Equations (7), (8), (9) and (11) form a closed system for the (uniform) static
pressure $p$ and the total pressure $p_{t,i}$, the total temperature $T_{t,i}$
and the cross-section $A_{i}$ in each stream. These variables fully define the
state of the flow and other variables can easily be calculated from them using
classic gas dynamic relations. The explicit nature of the equations
facilitates numerical integration with a shooting method by iterating on the
static pressure at the inlet to match an imposed static pressure at the
outlet.

Figure 1: We consider an axisymmetric nozzle with a contour $A(x)$ and two
streams, which have the same total temperature $T_{t}$, but can have a
different total pressure $p_{t}$. They exchange momentum through a shear force
$F_{ps}$ and a wall friction force $F_{w}$. The static pressure $p$ is assumed
to be uniform in each cross-section. We solve for the evolution of the flow
variables $p$, $p_{t,p}$, $p_{t,s}$, $T_{t,p}$, $T_{t,s}$ and the individual
cross-sections $A_{p}$ and $A_{s}$.

###  2.2 Model closure

We further assume that the compound flow is axisymmetric, with the primary
stream surrounded by the secondary stream. Hence, wall friction acts only on
the secondary stream:

$F_{w}=-\frac{1}{2}f_{w}\rho_{s}V_{s}^{2}l_{w}=-\frac{1}{2}f_{w}\gamma
p\mbox{{Ma}}_{s}^{2}l_{w},$ (14)

where the friction per unit length is denoted by $F_{w}$ and depends on the
friction coefficient $f_{w}$, the dynamic pressure of the secondary stream and
the wetted perimeter of the cross-section $l_{w}$. In an axisymmetric
geometry, $l_{w}$ can be related to the cross-section $A$ as follows:

$l_{w}=2\pi r=2\sqrt{\pi A},$ (15)

where $r$ denotes the radius of the total cross-section $A$. The shear is
defined in a similar way and depends on the difference in dynamic pressure.
This particular modelling choice facilitates its definition through the
uniform static pressure $p$:

$F_{ps}=\frac{1}{2}f_{ps}\gamma
p\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)l_{ps},$ (16)

where the shear force $F_{ps}$ is positive if the primary stream has a higher
Mach number (which is the most common case), $f_{ps}$ denotes the shear
coefficient and $l_{ps}$ denotes the length of the interface between both
streams. The latter coincides with the perimeter of the primary stream, hence:

$l_{ps}=2\pi r_{p}=2\sqrt{\pi A_{p}}.$ (17)

The shear forces counteracts the primary stream if it is positive, so the
primary force is defined as follows:

$F_{p}=-F_{ps}.$ (18)

The shear force is equal and opposite for the secondary stream, on which also
the (negative) friction force is applied:

$F_{p}=F_{ps}+F_{w}.$ (19)

Substituting the forces above in equation (11) eventually leads to the
following equations for the static pressure gradient:

$\frac{1}{p}\frac{dp}{dx}=\dfrac{1}{\beta}\left(\dfrac{dA}{dx}+\frac{1}{2}f_{ps}l_{ps}\frac{\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)^{2}}{\mbox{{Ma}}_{p}^{2}\mbox{{Ma}}_{s}^{2}}-\frac{1}{2}f_{w}l_{w}\left(1+\left(\gamma-1\right)\mbox{{Ma}}_{s}^{2}\right)\right).$
(20)

We recall that a compound-subsonic flow is characterised by a positive
compound choking indicator ($\beta>0$). As a result, the flow tends to
accelerate in a convergent section and due to friction and tends to recover
pressure in a divergent section and due to shear. Since the shear force
opposes the stream with the highest Mach number, it tends to equalise the Mach
numbers in both streams. The momentum exchange between both streams thus makes
the flow more uniform and leads to a net increase of static pressure. The
inverse is true for a compound-supersonic flow. The equation is singular in
the sonic section ($\beta=0$), which is analysed below.

###  2.3 Choking mechanism

The pressure gradient in a sonic section is bounded, which is only possible
according to equation (20) if both the numerator and the denominator equal
zero. In absence of forces, this leads to the classic result that a nozzle is
choked in the throat where the gradient of the cross-section vanishes (see for
example Bernstein et al. (1967)). In the more general case, the relative
magnitudes of the force terms in equation (20) determine the position of the
choked section.

Conveniently, the effect of shear and friction appear in two distinct terms.
Moreover, the shear term is always positive and the friction term is always
negative. This leads to some fundamental implications.

Friction has a negative contribution in the numerator, so in absence of shear,
the gradient of the cross-section must be positive. This implies that the
sonic section is situated in a divergent channel, hence downstream of the
throat. Likewise, equation (20) implies that shear pushes the sonic section
upstream towards a convergent section of the nozzle. The relative magnitude of
both terms will determine in which direction the sonic section will move. This
analysis offers a plausible explanation of the observed movement of the sonic
section by Lamberts et al. (2018) and Kracik & Dvorak (2023). We investigate
this effect in further detail by comparing with axisymmetric CFD simulations
in section 5.

The compound choking condition has also been derived by minimising the impulse
function (cf. Hoge & Segars (1965)), which is given by the following equation:

$J=\sum_{i\in[p,s]}\rho_{i}A_{i}V_{i}^{2}+pA_{i}=\sum_{i\in[p,s]}pA_{i}(1+\gamma\mbox{{Ma}}_{i}^{2}).$
(21)

Its derivative in space reads as follows (substituting the derivatives of the
cross-sections and the Mach numbers with the governing equations):

$\frac{dJ}{dx}=\left[\sum_{i\in[p,s]}A_{i}\frac{1-\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}\right]\frac{dp}{dx}-\left[\sum_{i\in[p,s]}F_{i}\frac{1-\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}\right],$
(22)

where we recognise the compound choking indicator $\beta$ in the first term.
Hence, the impulse function reaches an extreme value in the sonic section in
absence of forces, as derived by Hoge & Segars (1965). However, the second
term is not necessarily zero in case of shear and friction. We fill in the
forces defined in section 2.2:

$\frac{dJ}{dx}=\left[\sum_{i\in[p,s]}A_{i}\frac{1-\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}\right]\frac{dp}{dx}-\frac{1}{2}f_{ps}l_{ps}p\frac{\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)^{2}}{\mbox{{Ma}}_{p}^{2}\mbox{{Ma}}_{s}^{2}}+\frac{1}{2}f_{w}l_{w}p(1-\mbox{{Ma}}_{s}^{2}).$
(23)

Evidently, the terms on the right hand side generally do not add up to zero,
so the section with the minimal impulse function does not coincide with the
sonic section as defined through $\beta$. As shown by Bernstein et al. (1967),
the latter corresponds to the condition that compound waves can no longer
travel upstream, which constitutes the fundamental principle of choking.

##  3 Approximations to circumvent the singularity at the sonic section

So far, we have analysed choking by requiring the numerator to equal zero if
$\beta=0$. However, this leaves the pressure gradient undetermined $(0/0)$. In
this section, we develop an analytical expression for the pressure gradient in
the sonic section and a quasi-third order approximation in its vicinity. We
denote the numerator of equation (20) by $N$ to ease the notation of the
following derivation. We expand the numerator and the denominator in a Taylor
series around the sonic section located at $x=x_{t}$ up to the third order:

$\frac{1}{p}\frac{dp}{dx}=\frac{N}{\beta}=\frac{N\rvert_{t}+(x-x_{t})\left.\frac{dN}{dx}\right\rvert_{t}+\frac{1}{2}(x-x_{t})^{2}\left.\frac{d^{2}N}{dx^{2}}\right\rvert_{t}+\frac{1}{6}(x-x_{t})^{3}\left.\frac{d^{3}N}{dx^{3}}\right\rvert_{t}+...}{\beta\rvert_{t}+(x-x_{t})\left.\frac{d\beta}{dx}\right\rvert_{t}+\frac{1}{2}(x-x_{t})^{2}\left.\frac{d^{2}\beta}{dx^{2}}\right\rvert_{t}+\frac{1}{6}(x-x_{t})^{3}\left.\frac{d^{3}\beta}{dx^{3}}\right\rvert_{t}+...}.$
(24)

Both the numerator $N$ and the compound choking indicator $\beta$ equal zero
in the sonic section, so the equation above simplifies to the following
expression:

$\frac{1}{p}\frac{dp}{dx}=\frac{N}{\beta}=\frac{\left.\frac{dN}{dx}\right\rvert_{t}+\frac{1}{2}(x-x_{t})\left.\frac{d^{2}N}{dx^{2}}\right\rvert_{t}+\frac{1}{6}(x-x_{t})^{2}\left.\frac{d^{3}N}{dx^{3}}\right\rvert_{t}+...}{\left.\frac{d\beta}{dx}\right\rvert_{t}+\frac{1}{2}(x-x_{t})\left.\frac{d^{2}\beta}{dx^{2}}\right\rvert_{t}+\frac{1}{6}(x-x_{t})^{2}\left.\frac{d^{3}\beta}{dx^{3}}\right\rvert_{t}+...}.$
(25)

The first derivatives suffice to evaluate the pressure gradient in the sonic
section ($x=x_{t}$). However, numerical issues may arise in its vicinity
($x=x_{t}+\epsilon$), since the pressure gradient will stem from the division
of a small number by a small number in equation (20). The first order
approximation is constant and strictly only valid in the sonic section, so its
accuracy diminishes relatively quickly further away from the sonic section. A
second order approximation explicitly depends on the axial coordinate $x$ and
is more accurate in the vicinity of the sonic section. Therefore, we also
develop an analytical expression for the second derivatives. These expressions
are too long to include in the main body of this work, but can be found in
appendices A and B. The pressure gradient can thus be evaluated at different
orders of approximation: a first and a second order approximation can be
calculated analytically. Practically, a quasi-third order approximation can be
obtained by computing the third derivatives numerically, using finite
differences between the computed second derivatives.

To compute the required derivatives, we expand the derivatives of the
numerator and the denominator. This operation brings several derivatives
forwards, including those of the Mach numbers $\mbox{{Ma}}_{p}$ and
$\mbox{{Ma}}_{s}$, the cross-sections $A_{p}$ and $A_{s}$ (through $\beta$,
cf. equation (12)), the perimeter $l_{ps}$ of the interface between both
streams and the wetted perimeter of the wall $l_{w}$. The latter two can be
linked back to the derivatives of the cross-sections (cf. equations (17) and
(15)). The derivatives of the cross-sections can be substituted by equation
(9), which re-introduces the pressure gradient in the equation above. The Mach
numbers can be retrieved from the static and the total pressure as follows:

$\mbox{{Ma}}_{i}^{2}=\frac{2}{\gamma-1}\left(\left(\frac{p_{t,i}}{p}\right)^{\frac{\gamma-1}{\gamma}}-1\right).$
(26)

Hence, the gradients of the Mach number can be recast in terms of static and
total pressure gradients; the latter being substituted by equation (7). The
resulting expression for the pressure gradient consists of a combination of
explicit terms and implicit terms with the pressure gradient also appearing on
the right hand side. More specifically, the first derivatives of the numerator
and the denominator are proportional to the pressure gradient. Therefore, a
quadratic equation appears for the pressure gradient in the sonic section,
which is specified in appendix A. The two solutions correspond to a subsonic
or a supersonic flow downstream, depending on the sign of the pressure
gradient.

Numerically, one has to decide when to switch to the approximation. This can
be achieved by imposing a tolerance on $\beta$, e.g., using the approximation
when $-\epsilon<\beta<\epsilon$, with $\epsilon$ a small positive number. A
less crude condition is based on an estimate of the truncation error in the
Taylor series. We compute the first derivative after the truncation in the
numerator and the denominator to estimate their corresponding uncertainties
$\delta N$ and $\delta\beta$. These two uncertainties are propagated to the
pressure gradient as follows:

$\delta\left(\frac{1}{p}\frac{dp}{dx}\right)=\sqrt{\left(\frac{1}{\beta}\right)^{2}\left(\delta
N\right)^{2}+\left(-\frac{N}{\beta^{2}}\right)^{2}\left(\delta\beta\right)^{2}}.$
(27)

The equation above is analytically available for estimating the uncertainty on
the first order approximation, cf. the second derivatives in appendix B, but
the third derivatives need to be estimated through finite differences for the
second order approximation. The switch from the direct equation (20) to the
approximation can then be based on the relative uncertainty on the static
pressure gradient. A threshold of 1 % is used to obtain the results in section
5.

The state at the sonic section is required to compute the approximations, so
practically, the equations should first be integrated up to the sonic section
without any approximation of the pressure gradient. This is an iterative
process with a shooting method, where the static pressure at the inlet is
updated until the sonic condition is reached _and_ the numerator of equation
(20) also equal zero. The resulting sonic state can then be used to compute
the desired approximation (and the states in its vicinity in case finite
differences are required for the third derivatives). The shooting method
should then be repeated since the governing equations change slightly compared
to the direct method. Generally, it has been found that circumventing the
singularity with the approximations not only renders the integration more
stable, but is also beneficial for conservation of mass in the final solution
(note that the system of equations (7), (8), (9) and (11) is not
conservative).

##  4 Validation test cases with axisymmetric simulations

We compare the predictions of the 1D theoretical model above with post-
processed axisymmetic RANS simulations of a compound nozzle. Below, we
describe the numerical set-up and the test cases.

###  4.1 Numerical set-up

The computational domain is shown in figure 2. It consists of an axisymmetric
convergent-divergent nozzle with two distinct inlets. The simulation is two-
dimensional, which can be carried out in OpenFOAM by generating a 3D geometry
in the shape of a wedge with a single cell in the tangential direction. The
angle of the wedge equals 5 degrees. The contour of the upper wall is defined
as a cosine between $\pi/4$ and $\pi$ with appropriate scaling to reach the
dimensions shown in figure 2. The primary inlet extends from the centerline to
halfway the inlet. The numerical domain is intentionally asymmetrical with the
inlet closer to the throat to limit the extent of the shear layer at the sonic
section and thus to obtain more distinct states. The mesh has been generated
using Gmsh coupled to Python, which allows full control of the number of
elements and of the refinement near the walls.

The simulations are performed with an adapted version of the compressible
transient solver rhoCentralFoam in OpenFOAM v9
111https://openfoam.org/version/9/. The original solver is three-dimensional
and solves for swirling motion in an axisymmetric setting. The 1D compound
flow theory (cf. section 2) does not account for this effect, so it is not
desired in the 2D simulations. Therefore, the tangential component of the
velocity is initialised as zero and is not updated by the modified solver at
later times. A steady solution is obtained by letting the solver advance for a
sufficiently long time.

The boundary conditions at the inlets consist of the total pressure and total
temperature. The primary inlet is chosen to be supersonic, so also the
velocity is imposed. We arbitrarily choose a turbulence intensity of 5 %, with
a turbulent mixing length of 0.1 mm, which corresponds to one hundredth of the
radius at the inlet. We refer to table 1 in the next section for the numerical
values of the test cases. The flow is supersonic at the exit, so a zero
gradient condition is imposed. The wall is treated with a slip or no-slip
condition depending on the test case. It is always assumed to be adiabatic.
The sides of the wedge receive the ‘wedge’ boundary condition to inform
OpenFOAM of the axisymmetry. Sutherland’s law is used for the viscosity.

The resulting 2D flow fields need to be post-processed to compare them with
the predictions of the 1D compound flow theory. The first step consists of
identifying the dividing streamline between the streams. Its radial position
$r_{div}$ can be found at any point $x$ by integrating the axial mass flux as
follows:

$\dot{m}_{p}(x)=\int_{0}^{r_{div}}\rho(x,r)u(x,r)\ 2\pi r\
dr,\quad\dot{m}_{s}(x)=\int_{r_{div}}^{R}\rho(x,r)u(x,r)\ 2\pi r\ dr,$ (28)

where $u$ denotes the axial component of the velocity vector. The radius of
the dividing streamline is known at the primary inlet as it is part of the
definition of the numerical domain. This allows to calculate the primary mass
flow rate, and consequently, the evolution of the dividing streamline in $x$.
The second step consists of averaging the relevant flow variables across the
ejector’s section, i.e.

$\displaystyle\hat{\rho}_{p,s}(x)$
$\displaystyle=\frac{1}{A_{p,s}}\int_{A_{p,s}}\rho\ dA_{p,s},$ (29)
$\displaystyle\hat{u}_{p,s}(x)$
$\displaystyle=\frac{1}{\hat{\rho}_{p,s}A_{p,s}}\int_{A_{p,s}}\rho u\
dA_{p,s},$ (30) $\displaystyle\hat{e_{t}}_{p,s}(x)$
$\displaystyle=\frac{1}{\hat{\rho}_{p,s}A_{p,s}}\int_{A_{p,s}}\rho e_{t}\
dA_{p,s},$ (31)

where $e_{t}=c_{V}T+\frac{1}{2}V^{2}$ denotes the total specific internal
energy and $c_{V}$ is the specific heat capacity at constant volume. The
density-averaged velocity and internal energy ensure that the 1D mass flow
rate and internal energy match their integrated 2D counterpart. The averaged
Mach number follows from the temperature and the velocity:

$\displaystyle\hat{T}_{p,s}=$
$\displaystyle\frac{1}{c_{V}}\left(\hat{e_{t}}_{p,s}-\frac{1}{2}\hat{u}_{p,s}^{2}\right),$
(32) $\displaystyle\hat{\mbox{{Ma}}}_{p,s}=$
$\displaystyle\frac{\hat{u}_{p,s}}{\sqrt{\gamma R\hat{T}_{p,s}}},$ (33)

where $R$ denotes the specific gas constant.

Figure 2: The computational domain consists of an axisymmetric convergent-
divergent nozzle with two distinct inlets. The primary stream (index $p$) is
surrounded by the secondary stream (index $s$). The primary inlet extends
halfway of the complete inlet. All dimensions are in millimeter. The coarse
mesh is displayed for visibility.

###  4.2 Selected Test Cases

Three simulations are carried out to cover various conditions. The boundary
conditions are listed in table 1. The first simulation is inviscid and thus
serves as reference for the original compound flow theory without momentum
exchange. Secondly, shear is added numerically by switching to RANS
simulations with a slip condition at the wall. Finally, friction is added
through a no-slip condition.

A mesh convergence study is carried out on the viscid simulation with the no-
slip condition. The number of cells is doubled in each direction, leading to
meshes with $400\times 25$, $800\times 52$ and $1600\times 106$ cells. The
resulting distributions of $y^{+}$ near the wall and of the Mach number at the
centerline are shown for the three meshes on figure 3. The medium and fine
mesh show good agreement and both are wall-resolved with $y^{+}<1$. The fine
mesh is computationally not too expensive, so the results in the next section
are based on this mesh to profit from its higher resolution.

Simulation |  $p_{t,p}$ [bar] |  $p_{p}$ [bar] |  $T_{t,p}$ [K] |  $p_{t,s}$ [bar] |  $T_{t,s}$ [K] |  $TI$ [%] | Wall | Comments
---|---|---|---|---|---|---|---|---
1 | 3.0 | 1.2 | 300 | 1.5 | 300 | n/a | slip | inviscid
2 | 3.0 | 1.2 | 300 | 1.5 | 300 | 5 | slip | viscid
3 | 3.0 | 1.2 | 300 | 1.5 | 300 | 5 | no-slip | viscid
Table 1: Boundary conditions of the simulations presented in section 5. They
respectively cover an isentropic case, a case with shear but without friction
and a case with both forces.

(a) $y^{+}$ of the first cell near the wall.

(b) Mach number at the centerline.

Figure 3: The coarse, medium and fine mesh contain $400\times 25$, $800\times
52$ and $1600\times 106$ elements respectively. The medium and fine meshes
lead to a wall-resolved simulation with $y^{+}<1$ and to good agreement on the
Mach number at the centerline. The fine mesh is retained for the other
simulations and the post-processing.

##  5 Results

###  5.1 The effect of shear and friction on compound choking

Figure 4 shows the Mach number fields of the three simulations in table 1,
along with the dividing streamline in white and the sonic line in black. No
shear layer develops in the inviscid case, so the two streams remain clearly
distinct across the dividing streamline. The sonic line ends on the wall
downstream of the geometrical throat (indicated by the arrows). The compound
flow thus has a subsonic and a supersonic stream at the throat. The secondary
stream is accelerated through shear in the viscid simulations, causing the
sonic line to move radially outwards. As a result, the supersonic portion of
the compound flow in the throat is larger compared to the inviscid case. The
sonic line no longer reaches the wall when a no-slip condition is imposed. A
boundary layer develops, which tends to slow down the flow, as evidenced by
the lower maximal Mach number.

Figure 4: Mach number contours of the three simulations (cf. table 1). The
black line indicates the sonic line and tends towards the wall as the compound
flow becomes supersonic. The dividing streamline is indicated in white. The
arrows indicate the geometrical throat. The numerical domain (cf. figure 2)
has been reflected for visibility.

The averaged Mach number is computed with equations (28) - (33) and is
displayed with diamonds in figure 5 as the reference. Not all points are
displayed for visibility. We recall that the equivalent Mach number reflects
the complete compound flow and equals unity if $\beta=0$ (cf. equation (13)).
The procedure laid out in section 3 is applied with the same boundary
conditions as in table 1. Notably, this means that the Mach numbers at the
inlet are not imposed as strong boundary conditions, but are an indirect
result from the iterative procedure to find the static pressure at the inlet
that leads to a correctly choked flow (i.e., a numerator and denominator which
both equal 0). The only difference in the three cases are the shear and wall
friction coefficient which are used to obtain the 1D predictions. In the
inviscid case, both naturally equal 0. The shear coefficient equal 0.005 in
both viscid cases and the wall friction coefficient also equals 0.005 in the
simulation with the no-slip condition. These constant values have been found
through trial and error.

The key finding in these results is the movement of the sonic section, which
is highlighted by the vertical dashed line and commented below. The sonic
section is located in the geometrical throat in the inviscid case, as
predicted by the original compound flow theory (see for example Bernstein et
al. (1967)). The simulation with shear but without wall friction shows a
modest shift upstream towards a convergent section, as expected in the
analysis of section 2.3. Conversely, friction overcomes the effect of shear in
the final simulation and pushes the sonic section in a divergent part of the
nozzle.

(a) Inviscid, sonic at $x/L=0.333$.

(b) Slip, sonic at $x/L=0.317$.

(c) No slip, sonic at $x/L=0.372$.

Figure 5: Comparison of the Mach number computed with the 1D compound flow
theory and from the post-processed CFD results. Shear pushes the choked
section upstream in absence of friction (slip condition), but is overcome by
the friction term in equation (20), so the flow chokes in a diverging section.
The results were obtained with a constant shear coefficient $f_{ps}=0.005$ and
a constant friction coefficient $f_{w}=0.005$.

###  5.2 Sensitivity analysis

The solutions in the previous section allow to compute each term in the
equation for the static pressure gradient (cf. equation (20)) separately. The
evolution of the components in the numerator is shown on figure 6 for the slip
and no-slip case. The inviscid case is trivial since the gradient of the
cross-section is the only term that remains. We have normalised the terms with
the maximal gradient of the cross-section, as this is a fixed quantity for the
three simulations. The evolution of the shear and friction terms supports the
theoretical considerations in section 2.3, the first always being positive and
the second being negative. As a result, the zero-crossing of the numerator
moves up- or downstream. The gradient of the cross-section balances the effect
of the forces in either direction by moving the sonic section in a convergent
or divergent section. The fact that the sonic section matches with the zero-
crossing of the numerator stems from the iterative procedure to find the
choked solution, as detailed in section 3. Note that shear tends to decrease
from left to right, as the difference in Mach number decreases, while friction
tends to increase as the flow accelerates through the nozzle.

(a) Slip, sonic at $x/L=0.317$.

(b) No slip, sonic at $x/L=0.372$.

Figure 6: Components of the numerator in equation (20), normalised by the
maximal area gradient. The positive shear term pushes the zero-crossing of the
sum upstream, while the negative friction term pushes it downstream. The
inviscid case is not shown since it consists only of the area gradient.

##  6 Conclusion

We have presented an extension of the compound flow theory with shear and wall
friction. The resulting 1D model captures the up- and downstream movement of
the sonic section observed by Lamberts et al. (2018) and Kracik & Dvorak
(2023). Moreover, the 1D predictions have been validated numerically against
2D axisymmetric CFD simulations with appropriate boundary conditions to
include or eliminate the effect of shear and friction. To the knowledge of the
authors, our work is the first to propose an explanation for choking of a
compound flow _upstream_ of the geometrical throat.

The Taylor expansion of the equation for the static pressure gradient allows
to circumvent its singularity, which avoids numerical instabilities and
increases the accuracy of the numerical predictions. We have provided
analytical expressions for all required derivatives for a second order
approximation and a quasi-third order approximation can be obtained using
finite differences. These expressions, and also the governing equations, have
first been derived in a general way in terms of a force $F_{i}$. This explicit
choice facilitates future developments in case new forces should be included
or if the current forces should be defined differently. This would be the case
for a compound flow with more than 2 streams or for a planar geometry instead
of the axisymmetric one considered in this work. Another open path includes
heat exchange and different total temperatures at the inlet.

Finally, we believe that the developed model does not only bring more insight
in the choking of compound flows, but also opens perspectives in modelling,
for example in the field of ejectors.

[Acknowledgements]J. Van den Berghe is supported by a F.R.S.-FNRS FRIA grant.

##  Appendix A First order approximation of the static pressure gradient near
the sonic section

We first derive the expressions for the first derivative of the required
quantities.

### Shear

The first derivative of the shear force is given by the following equation:

$\dfrac{dF_{ps}}{dx}=C_{Fps0}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)+C_{Fps1},$
(34)

where

$\displaystyle C_{Fps0}=$ $\displaystyle
f_{ps}p\sqrt{\pi}\sqrt{A_{p}}\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)$
(35)
$\displaystyle+\dfrac{f_{ps}p\sqrt{\pi}\sqrt{A_{p}}}{2}\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)\left(\dfrac{1}{\mbox{{Ma}}_{p}^{2}}-1\right),$
(36) $\displaystyle C_{Fps1}=$ $\displaystyle
f_{ps}p\sqrt{\pi}\sqrt{A_{p}}\Bigg{[}\dfrac{2F_{p}}{pA_{p}}\left(\dfrac{\left(\gamma-1\right)\mbox{{Ma}}_{p}^{2}}{2}+1\right)-\dfrac{2F_{s}}{pA_{s}}\left(\dfrac{\left(\gamma-1\right)\mbox{{Ma}}_{s}^{2}}{2}+1\right)\Bigg{]}$
(37)
$\displaystyle-\dfrac{f_{ps}\gamma\sqrt{\pi}}{2\sqrt{A_{p}}}\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)F_{p}\left(\dfrac{1}{\gamma\mbox{{Ma}}_{p}^{2}}+\dfrac{\gamma-1}{\gamma}\right).$
(38)

### Wall friction

The first derivative of the wall friction force is given by the following
equation:

$\dfrac{dF_{w}}{dx}=c_{Fw0}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}+c_{Fw1},$
(39)

where

$\displaystyle c_{Fw0}=$ $\displaystyle-
f_{w}\sqrt{\pi}\sqrt{A}p\left(\mbox{{Ma}}_{s}^{2}-2\right),$ (40)
$\displaystyle c_{Fw1}=$ $\displaystyle-\dfrac{\sqrt{\pi}\left(f_{w}\gamma
A_{s}\mbox{{Ma}}_{s}^{2}p\dfrac{dA}{dx}+\left(2f_{w}\gamma-2f_{w}\right)AF_{s}\mbox{{Ma}}_{s}^{2}+4f_{w}AF_{s}\right)}{2\sqrt{A}A_{s}}.$
(41)

### The numerator of the equation for the static pressure (cf. equation (11))

We introduce the following notation for the numerator (cf. equations (11) and
(25)):

$N=\dfrac{dA}{dx}+\sum\limits_{i\in[p,s]}\left[\dfrac{1+\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}\right]\dfrac{F_{i}}{p}=\dfrac{dA}{dx}+\sum\limits_{i\in[p,s]}Z_{i}.$
(42)

We intentionally use the more general equation (11) instead of the equation
(20) for the two streams to facilitate potential modifications to the forces
in future work. The following equations will still hold for different forces:

$\dfrac{dZ_{i}}{dx}=C_{Z0,i}\left(\dfrac{dF_{i}}{dx}\right)+C_{Z1,i}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}+C_{Z2,i},$
(43)

where

$\displaystyle C_{Z0,i}=$
$\displaystyle\dfrac{1}{p}\left(\dfrac{1}{\gamma\mbox{{Ma}}_{i}^{2}}+\dfrac{\gamma-1}{\gamma}\right),$
(44) $\displaystyle C_{Z1,i}=$
$\displaystyle\dfrac{\left(\left(\gamma-\gamma^{2}\right)\mbox{{Ma}}_{i}^{4}-\mbox{{Ma}}_{i}^{2}+2\right)F_{i}}{\gamma^{2}p\mbox{{Ma}}_{i}^{4}},$
(45) $\displaystyle C_{Z2,i}=$
$\displaystyle-\dfrac{2F_{i}^{2}\left(\dfrac{\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}}{2}+1\right)}{\gamma^{2}p^{2}A_{i}\mbox{{Ma}}_{i}^{4}}.$
(46)

### The compound choking indicator $\beta$

We introduce the following notation for the compound choking indicator (cf.
equation (12)):

$\beta=\sum\limits_{i}\beta_{i}=\sum\limits_{i}A_{i}\dfrac{1-\mbox{{Ma}}_{i}^{2}}{\gamma\mbox{{Ma}}_{i}^{2}}.$
(47)

The first derivative of $\beta$ is given by the following equation:

$\dfrac{d\beta_{i}}{dx}=C_{\beta
0,i}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}+C_{\beta 1,i},$ (48)

where

$\displaystyle C_{\beta 0,i}=$
$\displaystyle\dfrac{A_{i}\left({\mbox{{Ma}}_{i}^{4}}+\left(\gamma-3\right)\mbox{{Ma}}_{i}^{2}+3\right)}{\gamma^{2}{\mbox{{Ma}}_{i}^{4}}},$
(49) $\displaystyle C_{\beta 1,i}=$
$\displaystyle\dfrac{\left(\gamma-1\right){\mbox{{Ma}}_{i}^{4}}+\left(3-2\gamma\right)\mbox{{Ma}}_{i}^{2}-3}{\gamma^{2}{\mbox{{Ma}}_{i}^{4}}}\left(\dfrac{F_{i}}{p}\right).$
(50)

### Static pressure gradient

The first derivative of the static pressure in the sonic section is given by
the following quadratic equation:

$\left[\sum\limits_{i}C_{\beta
0,i}\right]\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)^{2}+\left[\sum\limits_{i}C_{\beta
1,i}-\sum\limits_{i}C_{Z1,i}\right]\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)-\dfrac{d^{2}A}{dx^{2}}-\sum\limits_{i}C_{Z0,i}\left(\dfrac{dF_{i}}{dx}\right)-\sum\limits_{i}C_{Z2,i}=0$
(51)

Note that the first derivatives of the force $F_{i}$ also give rise to the
pressure gradient (see above), so they should be included in the linear term
of the quadratic equation. The first order approximation in the vicinity of
the sonic section consists of the same constant gradient in the throat given
by the equation above.

##  Appendix B Second order approximation of the static pressure gradient near
the sonic section

The second order approximation requires the first and the second derivatives
of the quantities from the section A. In this section, we do not expand the
first derivatives for conciseness, but instead refer to section A for their
computation.

### Cross-section

$\displaystyle\dfrac{d^{2}A_{i}}{dx^{2}}=$
$\displaystyle\dfrac{A_{i}}{\gamma}\left(\dfrac{1}{\mbox{{Ma}}_{i}^{2}}-1\right)\dfrac{d}{dx}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}+C_{A^{2}0,i},$
(52) $\displaystyle C_{A^{2}0,i}=$
$\displaystyle-\dfrac{A_{i}}{\gamma{\mbox{{Ma}}_{i}^{4}}}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}\left(\dfrac{d\mbox{{Ma}}_{i}^{2}}{dx}\right)+\dfrac{F_{i}}{p}\left(\dfrac{1}{\gamma\mbox{{Ma}}_{i}^{2}}+\dfrac{\gamma-1}{\gamma}\right)\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)$
$\displaystyle+\dfrac{1}{\gamma}\left(\dfrac{1}{\mbox{{Ma}}_{i}^{2}}-1\right)\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}\left(\dfrac{dA_{i}}{dx}\right)+\dfrac{F_{i}}{\gamma
p{\mbox{{Ma}}_{i}^{4}}}\left(\dfrac{d\mbox{{Ma}}_{i}^{2}}{dx}\right)-\dfrac{1}{p}\left(\dfrac{1}{\gamma\mbox{{Ma}}_{i}^{2}}+\dfrac{\gamma-1}{\gamma}\right)\left(\dfrac{dF_{i}}{dx}\right).$
(53)

### Mach number

$\displaystyle\dfrac{d^{2}\mbox{{Ma}}_{i}^{2}}{dx^{2}}=$
$\displaystyle-\dfrac{2}{\gamma}\left(\dfrac{\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}}{2}+1\right)\dfrac{d}{dx}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}+C_{M^{2}0,i},$
(54) $\displaystyle C_{M^{2}0,i}=$
$\displaystyle-\dfrac{\left(\gamma-1\right)}{\gamma}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}\left(\dfrac{d\mbox{{Ma}}_{i}^{2}}{dx}\right)+\dfrac{\left(\gamma-1\right)}{\gamma}\dfrac{F_{i}}{pA_{i}}\left(\dfrac{d\mbox{{Ma}}_{i}^{2}}{dx}\right)$
$\displaystyle+\dfrac{2}{\gamma
pA_{i}}\left(\dfrac{\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}}{2}+1\right)\left(\dfrac{dF_{i}}{dx}\right)-\dfrac{2F_{i}}{\gamma
pA_{i}}\left(\dfrac{\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}}{2}+1\right)\left(\dfrac{1}{p}\dfrac{dp}{dx}\right).$
(55)

### Shear

The second derivative of the shear force is given by the following equation:

$\dfrac{d^{2}F_{ps}}{dx^{2}}=C_{Fps^{2}0}\left(\dfrac{d}{dx}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}\right)+C_{Fps^{2}1},$
(56)

where

$\displaystyle C_{Fps^{2}0}=$ $\displaystyle
f_{ps}p\sqrt{\pi}\sqrt{A_{p}}\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)+\dfrac{f_{ps}p\sqrt{\pi}\sqrt{A_{p}}\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)}{2}\left(\dfrac{1}{\mbox{{Ma}}_{p}^{2}}-1\right),$
(57) $\displaystyle C_{Fps^{2}1}=$ $\displaystyle f_{ps}\gamma
p\sqrt{\pi}\sqrt{A_{p}}\left(C_{M^{2}0,p}-C_{M^{2}0,s}\right)+\dfrac{f_{ps}\gamma
p\sqrt{\pi}\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)}{2\sqrt{A_{p}}}C_{A^{2}0,p}$
$\displaystyle+\dfrac{f_{ps}\gamma
p\sqrt{\pi}\left(\dfrac{dA_{p}}{dx}\right)}{\sqrt{A_{p}}}\left(\dfrac{d\mbox{{Ma}}_{p}^{2}}{dx}-\dfrac{d\mbox{{Ma}}_{s}^{2}}{dx}\right)+2f_{ps}\gamma\sqrt{\pi}\sqrt{A_{p}}p\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}\left(\dfrac{d\mbox{{Ma}}_{p}^{2}}{dx}-\dfrac{d\mbox{{Ma}}_{s}^{2}}{dx}\right)$
$\displaystyle-\dfrac{f_{ps}\gamma
p\sqrt{\pi}\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)}{4{{A_{p}}^{\frac{3}{2}}}}\left(\dfrac{dA_{p}}{dx}\right)^{2}+\dfrac{f_{ps}\gamma\sqrt{\pi}\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right)p}{\sqrt{A_{p}}}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}\left(\dfrac{dA_{p}}{dx}\right)$
$\displaystyle+f_{ps}\gamma\sqrt{\pi}\sqrt{A_{p}}\left(\mbox{{Ma}}_{p}^{2}-\mbox{{Ma}}_{s}^{2}\right){{\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}}^{2}}p,$
(58)

### Wall friction

The second derivative of the wall friction force is given by the following
equation:

$\dfrac{d^{2}F_{w}}{dx^{2}}=c_{Fw^{2}0}\left(\dfrac{d}{dx}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}\right)+c_{Fw^{2}1}$
(60)

where

$\displaystyle c_{Fw^{2}0}$
$\displaystyle=f_{w}\sqrt{\pi}\sqrt{A}p\left(2-\mbox{{Ma}}_{s}^{2}\right),$
(61) $\displaystyle c_{Fw^{2}1}=$ $\displaystyle-
f_{w}\gamma\sqrt{\pi}\sqrt{A}pC_{M^{2}0,s}-\dfrac{f_{w}\gamma\sqrt{\pi}p}{\sqrt{A}}\left(\dfrac{dA}{dx}\right)\left(\dfrac{d\mbox{{Ma}}_{s}^{2}}{dx}\right)-2f_{w}\gamma\sqrt{\pi}\sqrt{A}p\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}\left(\dfrac{d\mbox{{Ma}}_{s}^{2}}{dx}\right)$
$\displaystyle-\dfrac{f_{w}\gamma\sqrt{\pi}\mbox{{Ma}}_{s}^{2}p}{2\sqrt{A}}\left(\dfrac{d^{2}A}{dx^{2}}\right)+\dfrac{f_{w}\gamma\sqrt{\pi}\mbox{{Ma}}_{s}^{2}p}{4{{A}^{\frac{3}{2}}}}\left(\dfrac{dA}{dx}\right)^{2}$
$\displaystyle-\dfrac{f_{w}\gamma\sqrt{\pi}\mbox{{Ma}}_{s}^{2}p}{\sqrt{A}}\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}\left(\dfrac{dA}{dx}\right)-f_{w}\gamma\sqrt{\pi}\sqrt{A}\mbox{{Ma}}_{s}^{2}{{\bigg{(}\dfrac{1}{p}\dfrac{dp}{dx}\bigg{)}}^{2}}p$
(62)

### The numerator of the equation for the static pressure (cf. equation (11))

The second derivative of the force term $Z_{i}$ (cf. equation (42)) is given
by the following equation:

$\displaystyle\dfrac{d^{2}Z_{i}}{dx^{2}}=$ $\displaystyle
C_{Z^{2}0,i}\left(\dfrac{d}{dx}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)\right)+C_{Z^{2}1,i}\left(\dfrac{{d^{2}F_{i}}}{dx^{2}}\right)+C_{Z^{2}2,i}\left(\dfrac{dF_{i}}{dx}\right)$
$\displaystyle+C_{Z^{2}3,i}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)^{2}+C_{Z^{2}4,i}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)+C_{Z^{2}5,i},$
(63)

where

$\displaystyle C_{Z^{2}0,i}=$
$\displaystyle-\dfrac{\left(\gamma^{2}-\gamma\right)F_{i}\mbox{{Ma}}_{i}^{4}+F_{i}\mbox{{Ma}}_{i}^{2}-2F_{i}}{\gamma^{2}p\mbox{{Ma}}_{i}^{4}},$
(64) $\displaystyle C_{Z^{2}1,i}=$
$\displaystyle\dfrac{\left(\gamma-1\right)\mbox{{Ma}}_{i}^{2}+1}{\gamma
p\mbox{{Ma}}_{i}^{2}},$ (65) $\displaystyle C_{Z^{2}2,i}=$
$\displaystyle-\dfrac{\left(2\gamma^{2}-2\gamma\right)\mbox{{Ma}}_{i}^{4}+2\mbox{{Ma}}_{i}^{2}-4}{\gamma^{2}p\mbox{{Ma}}_{i}^{4}}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)-\dfrac{\left(3\gamma-3\right)F_{i}\mbox{{Ma}}_{i}^{2}+6F_{i}}{\gamma^{2}p^{2}A_{i}\mbox{{Ma}}_{i}^{4}},$
(66) $\displaystyle C_{Z^{2}3,i}=$
$\displaystyle-\dfrac{\left(\gamma^{2}-\gamma^{3}\right)F_{i}\mbox{{Ma}}_{i}^{6}-F_{i}\mbox{{Ma}}_{i}^{4}+\left(6-2\gamma\right)F_{i}\mbox{{Ma}}_{i}^{2}-8F_{i}}{\gamma^{3}p\mbox{{Ma}}_{i}^{6}},$
(67) $\displaystyle C_{Z^{2}4,i}=$
$\displaystyle-\dfrac{\left(1-\gamma^{2}\right)F_{i}^{2}\mbox{{Ma}}_{i}^{4}+\left(5\gamma-9\right)F_{i}^{2}\mbox{{Ma}}_{i}^{2}+14F_{i}^{2}}{\gamma^{3}p^{2}A_{i}\mbox{{Ma}}_{i}^{6}},$
(68) $\displaystyle C_{Z^{2}5,i}=$
$\displaystyle\dfrac{\left(3\gamma-3\right){{F_{i}}^{3}}\mbox{{Ma}}_{i}^{2}\operatorname{+}6{{F_{i}}^{3}}}{\gamma^{3}p^{3}A_{i}^{2}\mbox{{Ma}}_{i}^{6}}.$
(69)

We group the terms excluding the second derivatives as follows:

$C_{Z^{2}2-5,i}=C_{Z^{2}2,i}\left(\dfrac{dF_{i}}{dx}\right)+C_{Z^{2}3,i}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)^{2}+C_{Z^{2}4,i}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)+C_{Z^{2}5,i}$
(70)

### The compound choking indicator $\beta$

The second derivative of $\beta_{i}$ is given by the following equation:

$\dfrac{d^{2}\beta_{i}}{dx^{2}}=C_{\beta^{2}0,i}\left(\dfrac{d}{dx}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)\right)+C_{\beta^{2}1,i}\left(\dfrac{dF_{i}}{dx}\right)+C_{\beta^{2}2,i}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)^{2}+C_{\beta^{2}3,i}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)+C_{\beta^{2}4,i}$
(71)

where

$\displaystyle C_{\beta^{2}0,i}=$
$\displaystyle\dfrac{A_{i}{\mbox{{Ma}}_{i}^{4}}+\left(\gamma-3\right)A_{i}\mbox{{Ma}}_{i}^{2}+3A_{i}}{\gamma^{2}\mbox{{Ma}}_{i}^{4}},$
(72) $\displaystyle C_{\beta^{2}1,i}=$
$\displaystyle\dfrac{\left(\gamma-1\right){\mbox{{Ma}}_{i}^{4}}+\left(3-2\gamma\right){\mbox{{Ma}}_{i}^{2}}-3}{\gamma^{2}p\mbox{{Ma}}_{i}^{4}},$
(73) $\displaystyle C_{\beta^{2}2,i}=$
$\displaystyle\dfrac{-A_{i}\mbox{{Ma}}_{i}^{6}+\left(\gamma^{2}-5\gamma+7\right)A_{i}{\mbox{{Ma}}_{i}^{4}}+\left(9\gamma-18\right)A_{i}\mbox{{Ma}}_{i}^{2}+15A_{i}}{\gamma^{3}\mbox{{Ma}}_{i}^{6}},$
(74) $\displaystyle C_{\beta^{2}3,i}=$
$\displaystyle\dfrac{\left(1-\gamma^{2}\right)F_{i}\mbox{{Ma}}_{i}^{6}-\left(2\gamma^{2}-10\gamma+10\right)F_{i}\mbox{{Ma}}_{i}^{4}+\left(30-19\gamma\right)F_{i}\mbox{{Ma}}_{i}^{2}-27F_{i}}{\gamma^{3}p\mbox{{Ma}}_{i}^{6}},$
(75) $\displaystyle C_{\beta^{2}4,i}=$
$\displaystyle\dfrac{\left(2{\gamma^{2}}-5\gamma+3\right)F_{i}^{2}\mbox{{Ma}}_{i}^{4}+\left(10\gamma-12\right)F_{i}^{2}\mbox{{Ma}}_{i}^{2}+12F_{i}^{2}}{\gamma^{3}p^{2}A_{i}\mbox{{Ma}}_{i}^{6}}.$
(76)

We group the terms excluding the second derivative of the pressure as follows:

$C_{\beta^{2}1-4,i}=C_{\beta^{2}1,i}\left(\dfrac{dF_{i}}{dx}\right)+C_{\beta^{2}2,i}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)^{2}+C_{\beta^{2}3,i}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)+C_{\beta^{2}4,i}$
(77)

### Static pressure gradient

The second derivative of the static pressure in the sonic section is given by:

$\displaystyle\dfrac{d}{dx}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)=\frac{1}{2\left(\dfrac{d\beta}{dx}\right)^{2}}\Bigg{[}$
$\displaystyle\left(\dfrac{d\beta}{dx}\left(\sum\limits_{i}C_{Z^{2}0,i}\right)-\left(\sum\limits_{i}C_{\beta^{2}0,i}\right)\left(\dfrac{d^{2}A}{dx^{2}}+\sum\limits_{i}\dfrac{dZ_{i}}{dx}\right)\right)\left(\dfrac{d}{dx}\left(\dfrac{1}{p}\dfrac{dp}{dx}\right)\right)$
$\displaystyle+\dfrac{d\beta}{dx}\left(\sum\limits_{i}C_{Z^{2}1,i}\right)\left(\dfrac{{d^{2}F_{i}}}{dx^{2}}\right)+\dfrac{d\beta}{dx}\left(\dfrac{d^{3}A}{dx^{3}}+\sum\limits_{i}C_{Z^{2}2-5,i}\right)$
$\displaystyle-\left(\sum\limits_{i}C_{\beta^{2}1-4,i}\right)\left(\dfrac{d^{2}A}{dx^{2}}+\sum\limits_{i}\dfrac{dZ_{i}}{dx}\right)\Bigg{]}.$
(78)

Note that the second derivative of the forces $F_{i}$ still gives rise to a
second derivative of the static pressure, but it depends on the definition of
the forces. The second order approximation of the first derivative of the
static pressure follows from the Taylor expansion in equation (25), for which
all terms can now be calculated up to the second derivatives. A quasi-third
order approximation can be obtained by computing the second derivatives on
multiple grid points and applying finite differences to approximate the
required third derivatives.

## References

  * Agnone (1971) Agnone, Anthony M 1971 Comment on” choked flow: A generalization of the concept and some experimental data”. AIAA Journal 9 (4), 768.
  * Bernstein et al. (1967) Bernstein, Arthur, Heiser, WH & Hevenor, Charles 1967 Compound-compressible nozzle flow. Journal of Applied Mechanics 34, 548––554.
  * Clark (1995) Clark, Larry 1995 Application of compound flow analysis to supersonic ejector-mixer performance prediction. In 33rd Aerospace Sciences Meeting and Exhibit, p. 645.
  * Fage (1976) Fage, E 1976 Apparent subsonic choking in dual-stream nozzles. AIAA Journal 14 (5), 681–683.
  * Hedges & Hill (1974) Hedges, KR & Hill, PG 1974 Compressible flow ejectors: Part ii—flow field measurements and analysis. Journal of Fluids Engineering 96 (3), 272––281.
  * Hoge & Segars (1965) Hoge, Harold J & Segars, Ronald A 1965 Choked flow-a generalization of the concept and some experimental data. AIAA Journal 3 (12), 2177–2183.
  * Kracik & Dvorak (2023) Kracik, Jan & Dvorak, Vaclav 2023 Effect of wall roughness on secondary flow choking in supersonic air ejector with adjustable motive nozzle. International Journal of Heat and Fluid Flow 103, 109168.
  * Lamberts et al. (2018) Lamberts, Olivier, Chatelain, Philippe, Bourgeois, Nicolas & Bartosiewicz, Yann 2018 The compound-choking theory as an explanation of the entrainment limitation in supersonic ejectors. Energy 158, 524–536.
  * Otis (1976) Otis, DR 1976 Choking and mixing of two compressible fluid streams. Journal of Fluids Engineering 98 (2), 311–317.
  * Papamoschou (1996) Papamoschou, Dimitri 1996 Analysis of partially mixed supersonic ejector. Journal of propulsion and power 12 (4), 736–741.
  * Papamoschou & Roshko (1988) Papamoschou, Dimitri & Roshko, Anatol 1988 The compressible turbulent shear layer: an experimental study. Journal of fluid Mechanics 197, 453–477.
  * Pearson et al. (1958) Pearson, H, Holiday, JB & Smith, SF 1958 A theory of the cylindrical ejector propelling nozzle”, jr aeron. The Journal of the Royal Aeronautical Society 62, 746–751.
  * Shapiro (1953) Shapiro, Ascher H 1953 The dynamics and thermodynamics of compressible fluid flow. New York: Ronald Press .


[affiliation=1,2]XiaopengWang [affiliation=1,*]RuiboFu
[affiliation=1,]ZhengqiWen [affiliation=1,2]ZhiyongWang
[affiliation=4]YuankunXie [affiliation=2]YukunLiu [affiliation=3]JianhuaTao
[affiliation=1]XuefeiLiu [affiliation=1]YongweiLi [affiliation=1,2]XinQi
[affiliation=1,2]YiLu [affiliation=5]ShuchenShi

# Genuine-Focused Learning using Mask AutoEncoder for Generalized Fake Audio
Detection

###### Abstract

The generalization of Fake Audio Detection (FAD) is critical due to the
emergence of new spoofing techniques. Traditional FAD methods often focus
solely on distinguishing between genuine and known spoofed audio. We propose a
Genuine-Focused Learning (GFL) framework guided, aiming for highly generalized
FAD, called GFL-FAD. This method incorporates a Counterfactual Reasoning
Enhanced Representation (CRER) based on audio reconstruction using the Mask
AutoEncoder (MAE) architecture to accurately model genuine audio features. To
reduce the influence of spoofed audio during training, we introduce a genuine
audio reconstruction loss, maintaining the focus on learning genuine data
features. In addition, content-related bottleneck (BN) features are extracted
from the MAE to supplement the knowledge of the original audio. These BN
features are adaptively fused with CRER to further improve robustness. Our
method achieves state-of-the-art performance with an EER of 0.25% on
ASVspoof2019 LA.

###### keywords:

fake audio detection, genuine-focused learning, counterfactual reasoning
enhanced representation, mask autoencoder

##  1 Introduction

Fake Audio Detection (FAD) is vital for ensuring information security and
preventing fraud. With the emergence of new spoofing techniques[1], developing
robust FAD measures for Automatic Speaker Verification (ASV) [2] systems has
become increasingly urgent. FAD systems must have the ability to generalize to
adapt to the evolving and unpredictable nature of spoofing attacks in the real
world. To promote the development of generalized research on FAD, the ASVspoof
[3, 4] initiative and its series of challenges have amassed a substantial
database of genuine and spoofed utterances. Of particular note, the test sets
contain numerous spoofing methods not present in the training sets, providing
valuable resources for evaluating the performance and generalization
challenges of FAD systems.

Currently, the detection of spoofed audio is mainly based on a two-step
approach involving feature extraction and classification tasks. The choice of
front-end feature extraction is critical for the model's generalizability.
Traditional approaches involve manual extraction of features such as first-
order spectral features [5, 6, 7, 8], second-order spectral features [9, 10,
11], and methods that fuse spectral features [12, 13]. These approaches have
somewhat improved the robustness of FAD; however, their capability to detect
unknown spoofing methods is still limited, and their generalizability needs
further enhancement.

In recent years, pre-trained models based on first-order spectral features
have demonstrated remarkable noise robustness and generalization capabilities
in FAD tasks [14, 15, 16, 17]. They use an encoder-only architecture, trained
with contrastive loss for the purpose of self-supervised learning [18, 19].
Their outstanding performance is due to pre-training on a large corpus of
genuine data with background noise. In contrast, pre-trained models based on
second-order spectral features use a different architecture [20, 21, 22],
typically an encoder-decoder structure where the decoder supports the encoder
during self-supervised training. However, similar to first-order spectral
features, when applied to FAD tasks, these models also function merely as
feature extractors, without further modeling the genuine audio patterns. These
patterns are present in the FAD training data. In addition, research on better
exploiting second-order spectral features pre-training for transfer learning
to FAD tasks is still very limited.

In this study, we introduce a novel Genuine-Focused Learning (GFL) framework
guided to enhance the generalization capability of FAD, called GFL-FAD. Unlike
previous self-supervised models, we adopt an encoder-decoder architecture. To
magnify the fine differences between genuine and spoofed audio, a
Counterfactual Reasoning Enhanced Representation (CRER) method is proposed in
the light of the audio reconstruction mechanism, where Mask AutoEncoder (MAE)
[23] architecture is employed for accurately modeling genuine audio features.
To reduce the influence of spoofed audio during training, we propose the
Genuine Audio Reconstruction Loss (GAR Loss) to maintain the focus on learning
genuine audio features, thereby deepening the understanding of the essence of
genuine audio and improving the detection ability of the model against unknown
spoofing techniques. In addition, content-related bottleneck (BN) features are
extracted from MAE to supplement the knowledge of the original audio, and we
use an attention mechanism to effectively fuse the BN features of MAE with
CRER, which improves the robustness of the model. Evaluated in the
ASVspoof2019 Logical Access (LA) scenario, our approach achieves state-of-the-
art performance with an Equal Error Rate (EER) of 0.25%.

  * •

We propose the GFL-FAD architecture, focusing on genuine audio to amplify the
differences between genuine and spoofed audio, thereby enhancing detection
capabilities.

  * •

We propose the GAR Loss to further focus on genuine learning and better
capture the distribution of genuine features, thereby improving the
generalizability of the model.

  * •

We use an attention mechanism to effectively fuse the BN features of MAE with
CRER, thereby increasing the robustness of the model.

Figure 1: Illustration of the overall architecture for the proposed GFL-FAD.
We start by segmenting patches from the spectral representation. A portion of
these patches are encoded by the MAE encoder, resulting in BN features. The
remaining patches are masked and concatenated with the BN features before
being fed to the MAE decoder for spectrogram reconstruction. We refer to the
reconstructed spectral features as CRER. If the input audio is genuine, an
additional GAR Loss is computed to maintain focus on modeling genuine audio
knowledge. Through fusion, the CRER is combined with the BN features and
passed to the back-end classification network for final detection.

##  2 Proposed method

The overall framework of our approach is shown in Figure 1. Unlike previous
self-supervised models, which are primarily used as feature extractors, we
jointly train the decoder of the MAE to enhance the discrimination between
genuine and spoofed audio(Section 2.1). To improve the method's focus on
genuine, we introduce the GAR Loss (Section 2.2). In addition, to preserve the
original information, we use a fusion module to combine the BN features of MAE
with CRER (Section 2.3) for the final classification task.

###  2.1 Counterfactual Reasoning Enhanced Representation

The CRER is extracted through the decoder of the MAE. Our MAE architecture is
based on the Audio Masked AutoEncoder (audio-MAE)[21], with the encoder
comprising 12 layers of ViT-Base (ViT-B), and the decoder being a 16 layers
transformer that includes a local attention mechanism. The spectrogram input
to the MAE must first be preprocessed, including segmented into small patches,
each of which is mapped to a $C$-dimensional embedding vector. All embedding
vectors are represented as $X=[x_{1},x_{2},\ldots,x_{N}]\in\mathbb{R}^{C\times
N}$, where $N=F\times T$ is the number of patches, $F$ is the number of sub-
bands, and $T$ is the number of temporal segments. Position embedding is added
to each patch to capture its location and the relative positional
relationships between patches. The formula is as follows:

$X^{p}=X+E_{pos}$ (1)

The sequence
$X^{p}=[x^{p}_{1},x^{p}_{2},\ldots,x^{p}_{N}]\in\mathbb{R}^{C\times N}$
consists of patch embedding with integrated positional information, while
$E_{pos}\in\mathbb{R}^{C\times N}$ denotes the positional embedding for all
patches. After obtaining $X^{p}$, certain areas of time and frequency domains
are randomly masked, and the unmasked patches are fed into the encoder to
encode the audio representations directly as content-related BN features. The
masked patches, along with the encoded audio representations, are then input
into the decoder for spectrogram reconstruction. The decoded features serve as
the CRER.

###  2.2 Genuine Audio Reconstruction Loss

A key component of our training strategy is the inclusion of the GAR Loss,
specifically designed for GFL, to enhance the model's ability to discriminate
between genuine and spoofed audio. The formula is as follows:

$L_{\text{GAR}}=\frac{\sum_{i=1}^{N}\text{mask}_{i}\cdot
L_{\text{recon},i}}{\sum_{i=1}^{N}\text{mask}_{i}}$ (2)

where $\text{mask}_{i}$ denotes the authenticity of sample $i$ (label 1 is the
genuine sample), and $L_{\text{recon},i}$ denotes the reconstruction loss for
the $i$-th sample. The formula for the reconstruction loss is as follows:

$L_{\text{recon},i}=\text{MSE}(\text{prediction}_{i},\text{input}_{i})$ (3)

We use the Mean Squared Error (MSE) as the reconstruction loss, aiming to
minimize the MSE between the CRER and the original spectrogram. In the
reconstruction loss, we calculate the loss only for the masked patch blocks,
rather than the entire spectrogram. Additionally, to better optimize the
model, we standardize the reconstructed CRER by its mean and variance.

###  2.3 Fusion Strategy for Feature Characteristics

In our fusion strategy, we employ the encoder segment of the Transformer model
[24] to integrate the BN features from MAE with CRER. We utilize the multi-
head attention mechanism of the Transformer, which conducts several self-
attention operations concurrently, to highlight the unique characteristics and
focal points of each feature set. The attention computation is formulated as
follows:

$\operatorname{Attention}(Q,K,V)=\operatorname{softmax}\left(\frac{QK^{T}}{\sqrt{d_{k}}}\right)V$
(4)

Here, $d_{k}$ denotes the dimension of the keys. CRER and BN features are
first mapped to the same dimension through two separate linear layers. the
embeddings of BN features to serve as queries ($Q$), with CRER acting as both
keys ($K$) and values ($V$). This configuration allows for a targeted fusion
of features, emphasizing the significant aspects of genuine audio present in
the BN features while enhancing them with the nuanced distinctions provided by
CRER.

The multi-head attention mechanism generates $h$ distinct representations of
($Q,K,V$), each undergoing a scaled dot-product attention calculation. The
outcomes from each head are then concatenated and passed through a feed-
forward layer, as described by:

$\text{MultiHead}(Q,K,V)=\text{Concat}(head_{1},\ldots,head_{h})W^{O}$ (5)
$head_{i}=\text{Attention}(QW_{i}^{Q},KW_{i}^{K},VW_{i}^{V})$ (6)

This fusion strategy emphasizes the strengths of each feature set, ensuring
that the integrated representation benefits from the detailed resolution of
the BN features and the discriminative power of the CRER.

###  2.4 Loss Function

In our framework, the total loss function $L_{\text{total}}$ is a blend of
Cross-Entropy Loss $L_{\text{CE}}$ for classification and Genuine Audio
Reconstruction Loss $L_{\text{GAR}}$ for audio reconstruction, formulated as:

$L_{\text{total}}=L_{\text{CE}}+\alpha\cdot L_{\text{GAR}}$ (7)

Here, $L_{\text{CE}}$ evaluates the model's ability to classify audio as
genuine or spoofed, directly impacting the accuracy of fake audio detection.
$L_{\text{GAR}}$ focuses on the reconstruction of genuine audio, ensuring that
the model retains the ability to accurately reproduce the characteristics of
genuine audio signals. The coefficient $\alpha$, initially set to 0.01, is a
non-learnable parameter that we investigate in ablation experiments to
understand its impact on performance. This optimization enhances the model's
focus on classification for spoofed audio while maintaining the integrity of
genuine audio reconstruction.

##  3 Experimental setup

###  3.1 Dataset details

Our experiments were conducted using the ASVspoof2019 LA dataset, which is
recognized for its extensive collection of genuine and spoofed audio samples,
produced via advanced text-to-speech and voice conversion techniques. The
evaluation set is particularly noted for encompassing a broader range of
spoofing attacks, featuring 13 distinct scenarios (A07-A19). This dataset and
the distribution of attacks across the partitions are summarized in Table 1.
Additionally, acknowledging the influence of initialization on the performance
of spoofing detection systems, the results presented are averaged over
multiple runs with varied random seeds, ensuring a reliable and consistent
performance evaluation.

Table 1: Overview of the ASVspoof2019 LA dataset Partition | genuine | Spoof | Attacks
---|---|---|---
Training | 2,580 | 22,800 | A01 - A06
Development | 2,548 | 22,296 | A01 - A06
Evaluation | 7,355 | 63,882 | A07 - A19

###  3.2 Implementation details

In our experimental setup, a raw waveform comprising 64,600 frames
(approximately 4 seconds) was acquired for analysis. We use a window size of
25 ms with a hop length of 10 ms to transform the waveform into 128 mel-bank
features. We initialized our MAE model using the pre-trained parameters of
audio-MAE [21]. We used AASIST [5] as our back-end classifier. Subsequently,
our proposed GFL-FAD method was implemented and trained using the PyTorch
framework for 100 epochs. The training process was conducted on an NVIDIA RTX
4090 GPU, utilizing a batch size of 16 to efficiently process the data. To
optimize the training procedure, we adopted the AdamW [25] optimizer with a
learning rate of $5\times 10^{-6}$. Additionally, a cosine annealing learning
rate decay scheme was employed to dynamically adjust the learning rate
throughout the training process, ensuring optimal model convergence and
performance.

For the evaluation of our proposed system, we utilized two principal metrics:
the EER and the Minimum Tandem Detection Cost Function (min t-DCF) [26], with
a particular emphasis on the EER for its direct interpretability. The min
t-DCF is also reported as it represents the standard for the ASVspoof
challenges, providing a comprehensive assessment in conjunction with ASV
systems.

##  4 Results

Figure 2: Visualization of High-Dimensional Representations of Genuine and Spoofed Audio Samples using T-SNE  Table 2: Performance Comparison of First-Order and Second-Order Spectral Methods With and Without Pre-Training in ASVspoof 2019 LA Evaluation System | Front-end | min | EER
---|---|---|---
|  | t-DCF | (%)
No pre-training |  |  |
LCNN-LSTM-sum [11] | LFCC | 0.0524 | 1.92
Attention+Resnet [27] | FFT | 0.0510 | 1.87
MCG-Res2Net50 [10] | CQT | 0.0520 | 1.78
FFT-L-SENet [28] | FFT | 0.0368 | 1.14
Raw PC-DARTS [29] | Raw waveform | 0.0517 | 1.77
Res-TSSDNet [8] | Raw waveform | 0.0481 | 1.64
RawGAT-ST [7] | Raw waveform | 0.0335 | 1.06
S2pecNet [12] | Raw & LFCC | 0.0240 | 0.84
AASIST [5] | Raw waveform | 0.0275 | 0.83
Graph-ST [9] | LFB | 0.0166 | 0.58
DFSincNet [30] | Raw waveform | 0.0176 | 0.52
F0-SENet34 [31] | STFT | 0.0143 | 0.43
pre-training |  |  |
HuRawNet2 [32] | Raw waveform | 0.1393 | 1.96
Wav2vec2+LGF [33] | Raw waveform | - | 1.28
Wav2vec2+AASIST [34] | Raw waveform | - | 0.90
Wav2vec2+LLGF [33] | Raw waveform | - | 0.86
Wav2vec2+VIB [35] | Raw waveform | 0.0107 | 0.40
Wav2vec2+ASP [15] | Raw waveform | - | 0.31
Ours | LFB | 0.0071 | 0.25

###  4.1 Performance Comparison

We visualized the high-dimensional representations of the model at its best
performance using T-SNE, as shown in Figure 2. The visualization confirms the
accurate fitting of genuine audio by the model, highlighting the effectiveness
of our approach in learning the patterns of genuine data. However, the
visualization also shows that some spoofed samples were incorrectly classified
as genuine. This may be because the artifacts in the spectra of these spoofed
samples are not pronounced, making it difficult for the model to accurately
detect them even when the differences are amplified.

In our performance analysis, we listed the performance of SOTA methods based
on first-order and second-order spectral features in terms of EER and min
t-DCF metrics on the ASVspoof2019 dataset, as shown in Table 2. F0-SENet34
[31] achieved the best performance by fusing the fundamental frequency with
real and imaginary spectrograms in scenarios without pre-training. In
scenarios with pre-trained models, the studies were based on first-order
spectral features pre-training, with the Wav2vec2+ASP [15] method achieving
the best performance. Our GFL-FAD, based on second-order spectral features
pre-training, achieved the SOTA performance among all of the above methods. In
addition, we compared the performance of our GFL-FAD method with that of using
only the encoder model pre-trained with audio-MAE in scenarios using the same
classifier, which also demonstrated the effectiveness of our proposed
approach. The specific performance details are provided in the ablation
experiments.

###  4.2 Ablation experiments

The ablation experiments performed on the ASVspoof2019 LA dataset demonstrate
the significant contributions of each component within our GFL-FAD approach,
as detailed in Table 3. Removing the GAR Loss component slightly increases the
EER to 0.39% and the min t-DCF to 0.0115, indicating the importance of GAR
Loss within the framework. Furthermore, the absence of the BN features (EN)
leads to an increase in the EER to 0.41%, demonstrating the positive impact of
these features on system performance. It is noteworthy that not using the CRER
features (DE), which essentially treats the pre-trained model encoder as a
mere feature extractor, results in the worst performance. We also explored the
optimal setting of the hyperparameter alpha, which balances detection and GAR
Loss. The results in Table 4 show that the best performance is achieved when
alpha is set to 0.01. In addition, we investigated the effect of different
mask ratios (time + frequency) on FAD, as shown in Figure 3. The study
revealed that excessively high mask ratios prevent MAE from learning effective
audio representations due to a significant loss of spectral information.
Conversely, excessively low mask ratios oversimplify the reconstruction task
and prevent effective learning by the MAE. Ultimately, the best performance
was achieved with a mask ratio of 0.3.

Figure 3: Performance of GFL-FAD at Different Mask Ratios Table 3: Ablation Experiments of Our Contributions Method | EER(% ) | min t-DCF
---|---|---
GFL-FAD | 0.25 | 0.0071
w/o GAR | 0.39 | 0.0115
w/o EN | 0.41 | 0.0119
w/o DE | 0.56 | 0.0184
Table 4: Selection of hyperparameter $\alpha$ $\alpha$ | EER(% ) | min t-DCF
---|---|---
1 | 0.54 | 0.0153
0.1 | 0.41 | 0.0122
0.01 | 0.25 | 0.0071

##  5 Conclusion

In this study, we present GFL-FAD, a novel framework that uses GFL to enhance
the generalization capabilities of FAD. This approach marks the first
application of an encoder-decoder structure, specifically the MAE, to FAD
tasks. Our method enhances the subtle differences between real and spoofed
audio through the CRER, which is informed by an audio reconstruction
mechanism. To minimize the impact of spoofed audio during training, we propose
the GAR loss, which ensures focused learning of genuine audio features. This
deepens our understanding of the nature of genuine audio and improves the
model's ability to detect novel spoofing methods. In addition, we extract BN
features from the MAE encoder and fuse them with CRER using an attention
mechanism to improve the robustness of the model. Our GFL-FAD demonstrates
state-of-the-art performance with an EER of 0.25% in the challenging
ASVspoof2019 LA scenario. In the future, we aim to explore more effective ways
to model genuine audio patterns and extend our methodology to other speech
anti-spoofing contexts, such as the detection of partially spoofed speech.

##  6 Acknowledgements

This work is supported by the National Natural Science Foundation of China
(NSFC) (No. 62101553, No. 62306316, No. U21B20210, No. 62201571).

## References

  * [1] R. Fu, J. Tao, Z. Wen, J. Yi, and T. Wang, ``Focusing on attention: prosody transfer and adaptative optimization strategy for multi-speaker end-to-end speech synthesis,'' in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 6709–6713.
  * [2] D. A. Reynolds, ``An overview of automatic speaker recognition technology,'' in _2002 IEEE international conference on acoustics, speech, and signal processing_ , vol. 4. IEEE, 2002, pp. IV–4072.
  * [3] M. Todisco, X. Wang, V. Vestman, M. Sahidullah, H. Delgado, A. Nautsch, J. Yamagishi, N. Evans, T. Kinnunen, and K. A. Lee, ``Asvspoof 2019: Future horizons in spoofed and fake audio detection,'' in _INTERSPEECH 2019-20th Annual Conference of the International Speech Communication Association_ , 2019.
  * [4] X. Wang, J. Yamagishi, M. Todisco, H. Delgado, A. Nautsch, N. Evans, M. Sahidullah, V. Vestman, T. Kinnunen, K. A. Lee _et al._ , ``Asvspoof 2019: A large-scale public database of synthesized, converted and replayed speech,'' _Computer Speech & Language_, vol. 64, p. 101114, 2020.
  * [5] J.-w. Jung, H.-S. Heo, H. Tak, H.-j. Shim, J. S. Chung, B.-J. Lee, H.-J. Yu, and N. Evans, ``Aasist: Audio anti-spoofing using integrated spectro-temporal graph attention networks,'' in _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2022, pp. 6367–6371.
  * [6] H. Tak, J.-W. Jung, J. Patino, M. Kamble, M. Todisco, and N. Evans, ``End-to-end spectro-temporal graph attention networks for speaker verification anti-spoofing and speech deepfake detection,'' in _ASVSPOOF 2021, Automatic Speaker Verification and Spoofing Countermeasures Challenge_. ISCA, 2021, pp. 1–8.
  * [7] H. Tak, J. weon Jung, J. Patino, M. Kamble, M. Todisco, and N. Evans, ``End-to-end spectro-temporal graph attention networks for speaker verification anti-spoofing and speech deepfake detection,'' in _Proc. 2021 Edition of the Automatic Speaker Verification and Spoofing Countermeasures Challenge_ , 2021, pp. 1–8.
  * [8] G. Hua, A. B. J. Teoh, and H. Zhang, ``Towards end-to-end synthetic speech detection,'' _IEEE Signal Processing Letters_ , vol. 28, pp. 1265–1269, 2021.
  * [9] F. Chen, S. Deng, T. Zheng, Y. He, and J. Han, ``Graph-based spectro-temporal dependency modeling for anti-spoofing,'' in _ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2023, pp. 1–5.
  * [10] X. Li, X. Wu, H. Lu, X. Liu, and H. Meng, ``Channel-Wise Gated Res2Net: Towards Robust Detection of Synthetic Speech Attacks,'' in _Proc. Interspeech 2021_ , 2021, pp. 4314–4318.
  * [11] X. Wang and J. Yamagishi, ``A comparative study on recent neural spoofing countermeasures for synthetic speech detection,'' _Interspeech 2021_ , 2021.
  * [12] P. Wen, K. Hu, W. Yue, S. Zhang, W. Zhou, and Z. Wang, ``Robust Audio Anti-Spoofing with Fusion-Reconstruction Learning on Multi-Order Spectrograms,'' in _Proc. INTERSPEECH 2023_ , 2023, pp. 271–275.
  * [13] C. Fan, J. Xue, S. Dong, M. Ding, J. Yi, J. Li, and Z. Lv, ``Subband fusion of complex spectrogram for fake speech detection,'' _Speech Communication_ , vol. 155, p. 102988, 2023.
  * [14] Y. Xie, H. Cheng, Y. Wang, and L. Ye, ``Domain generalization via aggregation and separation for audio deepfake detection,'' _IEEE Transactions on Information Forensics and Security_ , 2023.
  * [15] J. W. Lee, E. Kim, J. Koo, and K. Lee, ``Representation Selective Self-distillation and wav2vec 2.0 Feature Exploration for Spoof-aware Speaker Verification,'' in _Proc. Interspeech 2022_ , 2022, pp. 2898–2902.
  * [16] Y. Xie, H. Cheng, Y. Wang, and L. Ye, ``Learning a self-supervised domain-invariant feature representation for generalized audio deepfake detection,'' in _Proc. INTERSPEECH_ , vol. 2023, 2023, pp. 2808–2812.
  * [17] Z. Cai and M. Li, ``Integrating frame-level boundary detection and deepfake detection for locating manipulated regions in partially spoofed audio forgery attacks,'' _Computer Speech & Language_, vol. 85, p. 101597, 2024.
  * [18] A. Baevski, Y. Zhou, A. Mohamed, and M. Auli, ``wav2vec 2.0: A framework for self-supervised learning of speech representations,'' _Advances in neural information processing systems_ , vol. 33, pp. 12 449–12 460, 2020.
  * [19] S. Chen, C. Wang, Z. Chen, Y. Wu, S. Liu, Z. Chen, J. Li, N. Kanda, T. Yoshioka, X. Xiao _et al._ , ``Wavlm: Large-scale self-supervised pre-training for full stack speech processing,'' _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, no. 6, pp. 1505–1518, 2022.
  * [20] A. Baade, P. Peng, and D. Harwath, ``MAE-AST: Masked Autoencoding Audio Spectrogram Transformer,'' in _Proc. Interspeech 2022_ , 2022, pp. 2438–2442.
  * [21] P.-Y. Huang, H. Xu, J. Li, A. Baevski, M. Auli, W. Galuba, F. Metze, and C. Feichtenhofer, ``Masked autoencoders that listen,'' _Advances in Neural Information Processing Systems_ , vol. 35, pp. 28 708–28 720, 2022.
  * [22] Y. Gong, C.-I. Lai, Y.-A. Chung, and J. Glass, ``Ssast: Self-supervised audio spectrogram transformer,'' in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 36, no. 10, 2022, pp. 10 699–10 709.
  * [23] K. He, X. Chen, S. Xie, Y. Li, P. Dollár, and R. Girshick, ``Masked autoencoders are scalable vision learners,'' in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2022, pp. 16 000–16 009.
  * [24] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, ``Attention is all you need,'' _Advances in neural information processing systems_ , vol. 30, 2017.
  * [25] I. Loshchilov and F. Hutter, ``Decoupled weight decay regularization,'' in _International Conference on Learning Representations_ , 2018.
  * [26] T. Kinnunen, H. Delgado, N. Evans, K. A. Lee, V. Vestman, A. Nautsch, M. Todisco, X. Wang, M. Sahidullah, J. Yamagishi _et al._ , ``Tandem assessment of spoofing countermeasures and automatic speaker verification: Fundamentals,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 28, pp. 2195–2210, 2020.
  * [27] H. Ling, L. Huang, J. Huang, B. Zhang, and P. Li, ``Attention-based convolutional neural network for asv spoofing detection.'' in _Interspeech_ , 2021, pp. 4289–4293.
  * [28] J. Xue, C. Fan, Z. Lv, J. Tao, J. Yi, C. Zheng, Z. Wen, M. Yuan, and S. Shao, ``Audio deepfake detection based on a combination of f0 information and real plus imaginary spectrogram features,'' in _Proceedings of the 1st International Workshop on Deepfake Detection for Audio Multimedia_ , 2022, pp. 19–26.
  * [29] W. Ge, J. Patino, M. Todisco, and N. Evans, ``Raw differentiable architecture search for speech deepfake and spoofing detection,'' in _ASVSPOOF 2021, Automatic Speaker Verification and Spoofing Countermeasures Challenge_. ISCA, 2021, pp. 22–28.
  * [30] B. Huang, S. Cui, J. Huang, and X. Kang, ``Discriminative frequency information learning for end-to-end speech anti-spoofing,'' _IEEE Signal Processing Letters_ , vol. 30, pp. 185–189, 2023.
  * [31] J. Xue, C. Fan, Z. Lv, J. Tao, J. Yi, C. Zheng, Z. Wen, M. Yuan, and S. Shao, ``Audio deepfake detection based on a combination of f0 information and real plus imaginary spectrogram features,'' in _Proceedings of the 1st International Workshop on Deepfake Detection for Audio Multimedia_ , 2022, pp. 19–26.
  * [32] L. Li, T. Lu, X. Ma, M. Yuan, and D. Wan, ``Voice deepfake detection using the self-supervised pre-training model hubert,'' _Applied Sciences_ , vol. 13, no. 14, p. 8488, 2023.
  * [33] X. Wang and J. Yamagishi, ``Investigating Self-Supervised Front Ends for Speech Spoofing Countermeasures,'' in _Proc. The Speaker and Language Recognition Workshop (Odyssey 2022)_ , 2022, pp. 100–106.
  * [34] H. Tak, M. Todisco, X. Wang, J. weon Jung, J. Yamagishi, and N. Evans, ``Automatic speaker verification spoofing and deepfake detection using wav2vec 2.0 and data augmentation,'' in _Proc. The Speaker and Language Recognition Workshop (Odyssey 2022)_ , 2022, pp. 112–119.
  * [35] Y. Eom, Y. Lee, J. S. Um, and H. R. Kim, ``Anti-Spoofing Using Transfer Learning with Variational Information Bottleneck,'' in _Proc. Interspeech 2022_ , 2022, pp. 3568–3572.