1.频率下线化欧拉方程推导

有量纲形式的线化欧拉方程可以表示为

$\rho_{0}^{*}\left(\frac{\partial}{\partial t^{*}}+U^{*} \frac{\partial}{\partial x^{*}}\right) u^{*}+\rho_{0}^{*} \frac{\mathrm{d} U^{*}}{\mathrm{d} y^{*}} v^{*}=-\frac{\partial p^{*}}{\partial x^{*}}$ $\rho_{0}^{*}\left(\frac{\partial}{\partial t^{*}}+U^{*} \frac{\partial}{\partial x^{*}}\right) v^{*}=-\frac{\partial p^{*}}{\partial y^{*}}$ $\frac{1}{\rho_{0}^{*} c_{0}^{* 2}}\left(\frac{\partial}{\partial t^{*}}+U^{*} \frac{\partial}{\partial x^{*}}\right) p^{*}=-\left(\frac{\partial u^{*}}{\partial x^{*}}+\frac{\partial v^{*}}{\partial y^{*}}\right)$

用来无量纲的变量

$p=\frac{p^{*}}{\rho_{0}^{*} c_{0}^{* 2}}, \quad(x, y)=\left(\frac{x^{*}}{H^{*}}, \frac{y^{*}}{H^{*}}\right), \quad(u, v)=\left(\frac{u^{*}}{c_{0}^{*}}, \frac{v^{*}}{c_{0}^{*}}\right)$ $\omega=\frac{\omega^{*} H^{*}}{c_{0}^{*}}, \quad M(y)=M_{0} f(y)=\frac{U^{*}(y)}{c_{0}^{*}}, \quad t=\frac{t^{*} c_{0}^{*}}{H^{*}}$

无量纲线化欧拉方程推导：
动量方程x：

$\begin{array}{l} \frac{H^{*}}{\rho_{0}^{*} c_{0}^{* 2}}\left\{\rho_{0}^{*}\left(\frac{\partial}{\partial t^{*}}+U^{*} \frac{\partial}{\partial x^{*}}\right) u^{*}+\rho_{0}^{*} \frac{\mathrm{d} U^{*}}{\mathrm{d} y^{*}} v^{*}\right\}=-\frac{\partial p^{*}}{\partial x^{*}} \frac{H^{*}}{\rho_{0}^{*} c_{0}^{* 2}}=-\frac{\partial p}{\partial x} \\ \Leftrightarrow\left(\frac{\partial}{\partial t}+M(y) \frac{\partial}{\partial x}\right) u+v \frac{\mathrm{d} M(y)}{\mathrm{d} y}=-\frac{\partial p}{\partial x} \end{array}$

动量方程y：

$\left( \frac{\partial }{\partial t}\text{+}M\left( y \right)\frac{\partial }{\partial x} \right)v\text{=}-\frac{\partial p}{\partial y}$

能量方程：

$\begin{array}{l} \frac{1}{\rho_{0} c_{0}^{2}}\left(\frac{\partial}{\partial t^{*}}+U^{*} \frac{\partial}{\partial x^{*}}\right) p^{*}=-\left(\frac{\partial u^{*}}{\partial x^{*}}+\frac{\partial v^{*}}{\partial y^{*}}\right) \\ \Leftrightarrow \frac{H^{*}}{c_{0}^{*}}\left\{\frac{1}{\rho_{0}^{*} c_{0}^{* 2}}\left(\frac{\partial}{\partial t^{*}}+U^{*} \frac{\partial}{\partial x^{*}}\right) p^{*}\right\}=-\frac{H^{*}}{c_{0}^{*}}\left(\frac{\partial u^{*}}{\partial x^{*}}+\frac{\partial v^{*}}{\partial y^{*}}\right) \\ \Leftrightarrow\left(\frac{\partial}{\partial t}+M(y) \frac{\partial}{\partial x}\right) p=-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right) \end{array}$

考虑均匀流动 $M(y)=M_0 f(y)$ ，则有：
\begin{equation}
\left(\frac{\partial}{\partial t}+M_{0} f(y) \frac{\partial}{\partial x}\right) u+M_{0} \frac{\mathrm{d} f(y)}{\mathrm{d} y} v=-\frac{\partial p}{\partial x} \label{u}
\end{equation}

\begin{equation}
\left(\frac{\partial}{\partial t}+M_{0} f(y) \frac{\partial}{\partial x}\right) v=-\frac{\partial p}{\partial y} \label{v}
\end{equation}

\begin{equation}
\left(\frac{\partial}{\partial t}+M_{0} f(y) \frac{\partial}{\partial x}\right) p=-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\label{p}
\end{equation}
现在需要消去无量纲形式线化欧拉方程中的轴向速度项（消去轴向速度扰动项，可能是因为声衬影响的主要是法向的速度，将u消去后还可以减少计算成本），为消去关于u的时间和空间导数项,考虑对\eqref{u}整体做空间偏导数:
\begin{equation}
\begin{aligned}
\frac{\partial}{\partial x}\left\{\left(\frac{\partial}{\partial t}+M_{0} f(y) \frac{\partial}{\partial x}\right) u+M_{0} \frac{\mathrm{d} f(y)}{\mathrm{d} y} v\right\}=-\frac{\partial^{2} p}{\partial x^{2}} \\
\Leftrightarrow \frac{\partial^{2}}{\partial x \partial t} u+M_{0} f(y) \frac{\partial^{2} u}{\partial x^{2}}+M_{0} \frac{\mathrm{d} f(y)}{\mathrm{d} y} \frac{\partial v}{\partial x}=-\frac{\partial^{2} p}{\partial x^{2}}
\end{aligned} \label{ux}
\end{equation}
对式\eqref{p}整体做时间和空间上的偏导数：
\begin{equation}
\begin{aligned}
&\left(\frac{\partial}{\partial t}+M_{0} f(y) \frac{\partial}{\partial x}\right) p=-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\\
&\Leftrightarrow \frac{\partial}{\partial t}\left\{-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right\}+M_{0} f(y) \frac{\partial}{\partial x}\left\{-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right\}\\
&=\frac{\partial^{2} p}{\partial t^{2}}+2 M_{0} f(y) \frac{\partial^{2} p}{\partial x \partial t}+M_{0}^{2} f(y)^{2} \frac{\partial^{2} p}{\partial x^{2}}\\
&\Leftrightarrow-\frac{\partial^{2}}{\partial x \partial t} u-\frac{\partial^{2}}{\partial y \partial t} v-M_{0} f(y) \frac{\partial^{2} u}{\partial x^{2}}-M_{0} f(y) \frac{\partial^{2}}{\partial x \partial y} v\\
&=\frac{\partial^{2} p}{a t^{2}}+2 M_{0} f(y) \frac{\partial^{2} p}{\partial x \partial t}+M_{0}^{2} f(y)^{2} \frac{\partial^{2} p}{\partial x^{2}}
\end{aligned} \label{pxt}
\end{equation}
将式\eqref{ux}和式\eqref{pxt}相加:
\begin{equation}
\begin{aligned}
&\frac{\partial^{2}}{\partial x \partial t} u+M_{0} f(y) \frac{\partial^{2} u}{\partial x^{2}}+M_{0} \frac{\mathrm{d} f(y)}{\mathrm{d} y} \frac{\partial v}{\partial x}\\
&-\frac{\partial^{2}}{\partial x \partial t} u-\frac{\partial^{2}}{\partial y \partial t} v-M_{0} f(y) \frac{\partial^{2} u}{\partial x^{2}}-M_{0} f(y) \frac{\partial^{2}}{\partial x \partial y} v\\
&=-\frac{\partial^{2} p}{\partial x^{2}}+\frac{\partial^{2} p}{\partial t^{2}}+2 M_{0} f(y) \frac{\partial^{2} p}{\partial x \partial t}+M_{0}^{2} f(y)^{2} \frac{\partial^{2} p}{\partial x^{2}}\\
&\Leftrightarrow M_{0} \frac{\mathrm{d} f(y)}{\mathrm{d} y} \frac{\partial v}{\partial x}-\frac{\partial^{2} v}{\partial y \partial t}-M_{0} f(y) \frac{\partial^{2} v}{\partial x \partial y}\\
&=-\frac{\partial^{2} p}{\partial x^{2}}+\frac{\partial^{2} p}{\partial t^{2}}+2 M_{0} f(y) \frac{\partial^{2} p}{\partial x \partial t}+M_{0}^{2} f(y)^{2} \frac{\partial^{2} p}{\partial x^{2}}
\end{aligned} \label{tmp1-del-v}
\end{equation}

消去上式\eqref{tmp1-del-v}中的关于v的二阶时间偏导数，对式\eqref{v}求空间偏导数并与式\eqref{tmp1-del-v}相加：
\begin{equation}
\begin{aligned}
&\frac{\partial}{\partial y}\left\{\left(\frac{\partial}{\partial t}+M_{0} f(y) \frac{\partial}{\partial x}\right) v\right\}=-\frac{\partial^{2} p}{\partial y^{2}}\\
&\Leftrightarrow \frac{\partial^{2} v}{\partial y \partial t}+M_{0} f(y) \frac{\partial^{2} v}{\partial x \partial y}+M_{0} \frac{\partial v}{\partial x} \frac{d f(y)}{d y}=-\frac{\partial^{2} p}{\partial y^{2}}\\
&\Leftrightarrow 2 M_{0} \frac{\partial v}{\partial x} \frac{d f(y)}{d y}=-\frac{\partial^{2} p}{\partial y^{2}}-\frac{\partial^{2} p}{\partial x^{2}}+\frac{\partial^{2} p}{\partial t^{2}}+2 M_{0} f(y) \frac{\partial^{2} p}{\partial x \partial t}+M_{0}^{2} f(y)^{2} \frac{\partial^{2} p}{\partial x^{2}}
\end{aligned} \label{时间final-p}
\end{equation}
假设变量可以用以下的形式表示：
\begin{equation}
\begin{array}{l}
p=P(y) \exp (-\mathrm{i} k x) \exp (\mathrm{i} \omega t) \\
v=V(y) \exp (-\mathrm{i} k x) \exp (\mathrm{i} \omega t) \\
q=Q(y) \exp (-\mathrm{i} k x) \exp (\mathrm{i} \omega t)
\end{array} \label{频域变量}
\end{equation}

\eqref{频域变量}中 $q=i {\partial p}/{\partial x}$ ，将\eqref{频域变量}代入\eqref{v}和\eqref{时间final-p}可得：
\begin{equation}
\begin{array}{l}
\left[\mathrm{i} \omega-\mathrm{i} k M_{0} f(y)\right] V(y) \exp (-\mathrm{i} k x) \exp (\mathrm{i} \omega t)=-\frac{\mathrm{d} P(y)}{\mathrm{d} y} \exp (-\mathrm{i} k x) \exp (\mathrm{i} \omega t) \\
\Leftrightarrow \mathrm{i}\left[\omega-k M_{0} f(y)\right] V(y)=-\frac{\mathrm{d} P(y)}{\mathrm{d} y}
\end{array}
\end{equation}
以及
\begin{equation}
\begin{aligned}
&2 M_{0} \frac{\partial v}{\partial x} \frac{\mathrm{d} f(y)}{\mathrm{d} y}=-\frac{\partial^{2} p}{\partial y^{2}}-\frac{\partial^{2} p}{\partial x^{2}}+\frac{\partial^{2} p}{\partial t^{2}}+2 M_{0} f(y) \frac{\partial^{2} p}{\partial x \partial t}+M_{0}^{2} f(y)^{2} \frac{\partial^{2} p}{\partial x^{2}}\\
&\Leftrightarrow-2 i k M_{0} V(y) \frac{\mathrm{d} f(y)}{\mathrm{d} y}=-\frac{\mathrm{d}^{2} P(y)}{\mathrm{d} y^{2}}+k^{2} P(y)-\omega^{2} P(y)+2 M_{0} f(y) k \omega P(y)-M_{0}^{2} f(y)^{2} k^{2} P(y)\\
&\Leftrightarrow\left[1-M_{0}^{2} f(y)^{2}\right] k^{2} P(y)+2 \omega M_{0} f(y) k P(y)-\omega^{2} P(y)-\frac{\mathrm{d}^{2} P(y)}{\mathrm{d} y^{2}}=-2 i M_{0} k V(y) \frac{\mathrm{d} f(y)}{\mathrm{d} y}
\end{aligned}
\end{equation}
将频域下方程写成矩阵形式：
\begin{equation}
k\left(\begin{array}{ccc}
\mathbf{I}-M_{0}^{2} \mathbf{f}^{2} & 2 \mathbf{i} M_{0} \mathbf{f}_{\mathbf{a}} & 0 \\
0 & \mathbf{i} M_{0} \mathbf{f} & 0 \\
0 & 0 & \mathbf{I}
\end{array}\right)\left(\begin{array}{l}
\mathbf{Q} \\
\mathbf{V} \\
\mathbf{P}
\end{array}\right)=\left(\begin{array}{ccc}
-2 \omega M_{0} \mathbf{f} & 0 & \omega^{2} \mathbf{I}+\mathbf{D}_{2} \\
0 & \mathbf{i} \omega \mathbf{I} & \mathbf{D}_{1} \\
\mathbf{I} & 0 & 0
\end{array}\right)\left(\begin{array}{l}
\mathbf{Q} \\
\mathbf{V} \\
\mathbf{P}
\end{array}\right) \label{final-矩阵}
\end{equation}
其中，式\eqref{final-矩阵}的变量的具体表述为：

$\begin{aligned} &\mathbf{f}=\left[\begin{array}{cccc} f & 0 & \cdots & 0 \\ 0 & f & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & f \end{array}\right], \mathbf{f}^{2}=\left[\begin{array}{cccc} f^{2} & 0 & \cdots & 0 \\ 0 & f^{2} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & f^{2} \end{array}\right]\\ &\mathbf{f}_{a}=\left[\begin{array}{cccc} \mathrm{d} f / \mathrm{d} y & 0 & \cdots & 0 \\ 0 & \mathrm{d} f / \mathrm{d} y & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \mathrm{d} f / \mathrm{d} y \end{array}\right] \end{aligned}$

2.散射矩阵如何构造

本小节试着构造一个区域的散射矩阵，由多个区域组成的散射矩阵

2.1 单一区域中散射矩阵的构造

在给出了求解的广义特征值问题后，现在针对带空腔细节结构的管道进行分析，考虑如下图所示的管道
[cavity]
区域j中的第n个特征向量，可以表示为 $\left(\boldsymbol{Q}_{n}^{j}, \boldsymbol{V}_{n}^{j}, \boldsymbol{P}_{n}^{j}\right) ^T$
在每个区域内，模态解可以通过区域内各个特征向量的线性组合得到：
\begin{equation}
\left.\begin{array}{rl}
Q^{j}(x) & =\sum_{n=1}^{N} C_{n}^{j} Q_{n}^{j} \exp \left(-\mathrm{i} k_{n}^{j} x\right) \\
P^{j}(x) & =\sum_{n=1}^{N} C_{n}^{j} P_{n}^{j} \exp \left(-\mathrm{i} k_{n}^{j} x\right) \\
V^{j}(x) & =\sum_{n=1}^{N} C_{n}^{j} V_{n}^{j} \exp \left(-\mathrm{i} k_{n}^{j} x\right)
\end{array}\right\}
\end{equation}

管道中I区域采用 $N_1$ 个点离散，区域II采用 $N_2$ 个点离散。假设空腔中没有速度扰动，所以区域I中有 $N=3N_1$ 个模态，区域II中有 $N=2N_2+N_1$ 个模态。
The modes in each duct segment are then matched using the continuity of pressure $p$, velocity $v$ and $\partial {p}/\partial{x}$ at the interfaces between zones I and II, and $\partial {p}/\partial{x}=0$ on the vertical walls inside the cavity. The continuity and wall conditions can be put in the form of a large matrix that links all the incoming waves in the cell to outgoing waves and to all the internal variables. From this large matrix, the scattering matrix is written:
\begin{equation}
\left(\begin{array}{c}
\boldsymbol{C}_{2}^{+} \\
\boldsymbol{C}_{1}^{-}
\end{array}\right)=\boldsymbol{S}\left(\begin{array}{c}
\boldsymbol{C}_{1}^{+} \\
\boldsymbol{C}_{2}^{-}
\end{array}\right) \label{散射矩阵}
\end{equation}

矢量 $C_{1}^{\pm}$ 包含了$x=0$处沿流动方向（与流动方向相反）传播的波的模态系数，矢量 $C_{2}^{\pm}$ 表征$x=W$的模态系数。
\begin{equation}
S=\left(\begin{array}{cc}
T^{+} & R^{-} \\
R^{+} & T^{-}
\end{array}\right)
\end{equation}
$T^{+}\left(2 N_{1} \times 2 N_{1}\right)$,$R^{+}\left(N_{1} \times 2 N_{1}\right)$,$T^{-}\left(N_{1} \times N_{1}\right)$和$R^{-}\left(2 N_{1} \times N_{1}\right)$为沿流动方向和与流动方向相反的传输和反射矩阵。在区域I中共存在$3N_{1}$个模态，$N_{1}$个声学模态沿x正方向传播，$N_{1}$个升邪模态沿x负方向传播，$N_{1}$个流动模态沿x方向传播。

$ C_{2}^{+}$表示的是$x=W$处往下游传播的波的系数，它包含了$N_{1}$个声学模态和$N_{1}$个流动模态，所以沿流动方向的传输矩阵$T^{+}\left(2 N_{1} \times 2 N_{1}\right)$和反射矩阵$R^{-}\left(2 N_{1} \times N_{1}\right)$的行维度为$2N_{1}$。
$ C_{1}^{-}$表示的是$x=0$处往上游传播的波的系数，它只有声学模态，由$x=0$处传输的波和$x=W$传输回来的声波决定，所以$R^{+}\left(N_{1} \times 2 N_{1}\right)$和$T^{-}\left(N_{1} \times N_{1}\right)$的行维度为$N_{1}$
T和R列维度的大小由$ C_{1}^{+}$和$ C_{2}^{-}$决定。

2.2 多个区域中散射矩阵的构造

经过多个区域的散射矩阵可以由各个区域里的散射矩阵$\mathbf{S}^{(1)}$、$\mathbf{S}^{(2)}$组成：
\begin{equation}
\mathbf{S}^{(1+2)}=\left(\begin{array}{cc}
\boldsymbol{T}^{+(2)} \boldsymbol{E} \boldsymbol{T}^{+(1)} & \boldsymbol{R}^{-(2)}+\boldsymbol{T}^{+(2)} \boldsymbol{R}^{-(1)} \boldsymbol{F} \boldsymbol{T}^{-(2)} \\
\boldsymbol{R}^{+(1)}+\boldsymbol{T}^{-(1)} \boldsymbol{R}^{+(2)} \boldsymbol{E} \boldsymbol{T}^{+(1)} & \boldsymbol{T}^{-(1)} \boldsymbol{F} \boldsymbol{T}^{-(2)}
\end{array}\right)
\end{equation}
其中$\boldsymbol{E}$，$\boldsymbol{F}$的具体形式如下：
\begin{equation}
\begin{aligned}
&\boldsymbol{E}=\left(\boldsymbol{I}-\boldsymbol{R}^{-(1)} \boldsymbol{R}^{+(2)}\right)^{-1}\\
&\boldsymbol{F}=\left(\boldsymbol{I}-\boldsymbol{R}^{+(2)} \boldsymbol{R}^{-(1)}\right)^{-1}
\end{aligned}
\end{equation}

2.3 经过一个周期后散射矩阵的构造（Bloch modes of the periodic system）

当几何、边界条件和控制方程在x方向关于长度W周期对称时，Floquet–Bloch理论认为任何的物理量可以表示成：
\begin{equation}
\phi(x+W, y)=\phi_{B}(x+W, y) \mathrm{e}^{-\mathrm{i} k_{B}(x+W)}=\phi_{B}(x, y) \mathrm{e}^{-\mathrm{i} k_{B}(x+W)}=\phi(x, y) \mathrm{e}^{-\mathrm{i} k_{B} W}
\end{equation}
$\phi_{B}(x, y)$为在x方向上关于长度W对称的物理场。在研究的周期结构上下游边界处（$x=0$,$x=W$）,存在以下的关系式：
\begin{equation}
\phi(W, y)=\phi(0, y) \mathrm{e}^{-i k_{B} W}
\end{equation}
其矩阵形式为：
\begin{equation}
\left(\begin{array}{l}
\boldsymbol{C}_{2}^{+} \\
\boldsymbol{C}_{2}^{2}
\end{array}\right)=\mathrm{e}^{-\mathrm{i} \boldsymbol{k}_{\boldsymbol{B}} \mathrm{W}}\left(\begin{array}{l}
\boldsymbol{C}_{1}^{+} \\
\boldsymbol{C}_{1}^{1}
\end{array}\right) \label{周期关系}
\end{equation}
考虑式\eqref{散射矩阵}所示的散射矩阵，将\eqref{散射矩阵}重写成：
\begin{equation}
M_{1}\left(\begin{array}{c}
C_{1}^{+} \\
C_{1}^{-}
\end{array}\right)=M_{2}\left(\begin{array}{c}
C_{2}^{+} \\
C_{2}^{-}
\end{array}\right) \label{散射矩阵重写}
\end{equation}
其中
\begin{equation}
M_{1}=\left(\begin{array}{cc}
T^{+} & 0 \\
-R^{+} & I
\end{array}\right), \quad M_{2}=\left(\begin{array}{cc}
I & -R^{-} \\
0 & T^{-}
\end{array}\right)
\end{equation}
结合\eqref{散射矩阵重写}和\eqref{周期关系}，可以得到经过一个周期结构后的散射矩阵：
\begin{equation}
M_{1}\left(\begin{array}{c}
\boldsymbol{C}_{1}^{+} \\
\boldsymbol{C}_{1}^{-}
\end{array}\right)=\mathrm{e}^{-\mathrm{i} k_{B} W} \boldsymbol{M}_{2}\left(\begin{array}{c}
\boldsymbol{C}_{1}^{+} \\
\boldsymbol{C}_{1}^{-}
\end{array}\right)
\end{equation}

2.4 散射矩阵中$\boldsymbol{S}$的构造——Mode matching

区域交界面处满足连续性条件
空腔内部垂直壁面上满足$\partial {p}/\partial{x}=0$

\begin{equation}
\underbrace{\left(\begin{array}{cc}
-\mathbf{Q}_{1}^{-} & \mathbf{Q}_{2}^{+} \\
\mathbf{0} & \mathbf{V}_{2}^{+} \\
-\mathbf{V}_{1}^{-} & \mathbf{P}_{2}^{+}\left(1: N_{1},:\right)
\end{array}\right)}_{\mathbf{S}_{1}}\left(\begin{array}{c}
\mathbf{C}_{1}^{-} \\
\mathbf{C}_{2}^{+}
\end{array}\right)=\underbrace{\left(\begin{array}{cc}
\mathbf{Q}_{1}^{+} & -\mathbf{Q}_{2}^{-} \\
\mathbf{0} & -\mathbf{V}_{2}^{-} \\
\mathbf{V}_{1}^{+} & -\mathbf{P}_{2}^{-}\left(1: N_{1},:\right)
\end{array}\right)}_{\mathbf{S}_{2}}\left(\begin{array}{c}
\mathbf{C}_{1}^{+} \\
\mathbf{C}_{2}^{-}
\end{array}\right)
\end{equation}

\begin{equation}
\begin{array}{l}
R^{+}=\frac{\mathbf{C}_{1}^{-}(1)}{\mathbf{C}_{1}^{+}(1)}=\mathbf{S}(1,1), \quad T^{-}=\frac{\mathbf{C}_{1}^{-}(1)}{\mathbf{C}_{2}^{-(1)}}=\mathbf{S}\left(1,2 N_{1}+1\right) \\
T^{+}=\frac{\mathbf{C}_{2}^{+}(1)}{\mathbf{C}_{1}^{+}(1)}=\mathbf{S}\left(N_{1}+1,1\right), \quad R^{-}=\frac{\mathbf{C}_{2}^{+}(1)}{\mathbf{C}_{2}^{-}(1)}=\mathbf{S}\left(N_{1}+1,2 N_{1}+1\right)
\end{array}
\end{equation}

3.Multimodal method

By matching the pressure and axial velocity at the interface between different uniform segments, scattering matrices are obtained for each individual segment;
these are then combined to construct a global scattering matrix for multiple segments.