2016, Vol. 33| 磁流体方腔槽道流整体线性稳定性数值方法研究 |  [PDF全文] |
2. 中国科学院大学物理学院, 北京 100049
2. University of Chinese Academy of Sciences, Beijing 100049, China
磁场作用下的金属导电流体(磁流体)方腔槽道流在钢铁工业的连续铸钢工艺[1, 2]和磁约束聚变反应堆的液态包层系统[3, 4]中有重要的工程应用背景.尤其是对于液态锂铅包层系统，磁流体方腔槽道流是其中的重要流动类型.而同时液态锂铅包层是最具现实可行性和发展潜力的包层概念，其中的液态金属流动的磁流体动力学效应是值得关注的重要问题，由磁流体动力学效应带来的压降问题、分流影响以及流动传热效率问题对液态锂铅包层系统的实现起着至为关键的作用.
导电流体方腔槽道流在磁场作用下有不同于普通绝缘流体的动力学效应.首先，在磁场的作用下磁流体方腔槽道流的速度流型发生变化，会产生拐点[5]、剪切层[6]甚至射流[7]，这些速度流型结构会使得流动更加不稳定；其次，磁场的作用有抑制三维扰动的能力，使流场沿磁场方向更加均匀分布，由此恒定磁场作用下的磁流体流场往往是一个准二维结构[8, 9]；换句话说，恒定磁场的作用会抑制磁流体方腔槽道流中湍流的形成，使层流至湍流的转捩发生在更高的雷诺数下[10].这两种效应使得磁场作用下的磁流体方腔槽道流层流至湍流的转捩在某些情况下可以从线性失稳开始，Kinet和Knaepen[11]对磁流体Hunt射流的数值模拟证实了存在这种从线性失稳开始的转捩过程.
磁流体方腔槽道流的整体线性稳定性研究对于层流至湍流转捩机制的认识具有重要的理论指导价值.Priede等人[12]首先采用Chebyshev谱配置法对磁流体Hunt槽道流进行了整体线性稳定性分析，他们利用了方腔截面的对称性，结果表明，大Hartmann数下的临界雷诺数近似为Rec=112(以平均速度作为特征速度)，可见Hunt槽道流是在较小的雷诺数下发生失稳的.Priede等人[13]进一步对完全导电壁面的磁流体方腔槽道流进行了整体线性稳定性分析，同样发现在中等雷诺数下就可以发生线性失稳.对于壁面完全绝缘的磁流体槽道流，在磁场强度较弱的情况下对任意雷诺数流动是线性稳定的，而在强磁场作用下流动会发生线性失稳，但发生在非常高的雷诺数下(104量级).Krasnov[14]采用伴随优化方法考察了壁面完全绝缘的磁流体方腔槽道流的线性最优瞬态增长，发现磁场对最优瞬态增长有抑制作用，并且最优时空结构出现在侧层(与磁场方向相平行的壁面)附近.
磁流体方腔槽道流的整体线性稳定性特征模态通常采用谱配置法来计算，但还未见用于计算非正则模线性最优瞬态增长.这主要是由于谱配置法为了计算大量叠加所需的线性模态，需要构造非常庞大的满矩阵，存储量和计算开销非常巨大.为此，Paredes等人[15]提出可以将基于非均匀网格的高阶FD-q差分格式运用于方腔槽道流的整体线性稳定性特征模态的计算，他们的计算分析表明高阶FD-q差分格式相对于谱方法计算效率可以达到1万倍的提升.本文重点比较Chebyshev谱配置法和高阶FD-q差分格式的计算特点，在验证高阶FD-q差分格式高效性的基础上，直接通过线性模态叠加的方式来计算磁流体方腔槽道流的非正则模线性最优瞬态增长，计算结果与Krasnov[14]采用伴随优化方法得到的结果一致.
1 磁流体方腔槽道流控制方程 1.1 准静态近似下的磁流体控制方程在常压力梯度的作用下，不可压粘性导电流体在垂直均匀磁场的作用下沿方腔槽道流动(简称磁流体方腔槽道流).如图 1所示，流体沿x方向流动，外加磁场方向沿z方向，坐标原点位于方形截面中心处.磁流体密度为ρ，运动粘性为ν，电导率为σ；正方形截面半高度(也是半宽度)为L；外加均匀磁场BeB沿z方向.垂直磁场的上下壁面是理想完全导电的，而平行的左右侧壁面则是绝缘的，这样的流动称为磁流体Hunt流动[7].适用于低磁雷诺数下的准静态近似不可压缩磁流体无量纲控制方程:
$$\begin{array}{c} \nabla \cdot \boldsymbol{v} = 0,\\ \frac{{\partial \boldsymbol{v}}}{{\partial t}} + \left( {\boldsymbol{v} \cdot \nabla } \right)\boldsymbol{v} = - \nabla p + \frac{1}{{Re}}\Delta \boldsymbol{v} + \frac{{H{a^2}}}{{Re}}\left( {\boldsymbol{J} \times {\boldsymbol{e}_{\rm{B}}}} \right),\\ \nabla \cdot \boldsymbol{J} = 0,\\ \boldsymbol{J} = - \nabla \phi + \boldsymbol{v} \times {\boldsymbol{e}_{\rm{B}}}, \end{array}$$ 其中，v=(u,v,w)，J=(Jx,Jy,Jz)，p，ϕ，分别表示速度矢量、电流密度矢量、压力及电势.控制方程中的无量纲参数有Reynolds数和Hartmann数，它们的定义为 $$Re = {U_0}L/v,\;\;Ha = BL\sqrt {\sigma /\rho v} .$$ 这里， U0为基本流速度的最大值.Hartmann数为零时就退化成无磁场下的方腔槽道流动. |
| 图1 磁流体Hunt槽道流流动模型 Fig. 1 Magnetohydrodynamic(MHD)Hunt duct flow model |
壁面无量纲边界条件为
无滑移边界：v×n=0，
无穿透边界：v·n=0，
导电壁面：ϕ=0，
绝缘壁面：n·Δϕ=0.
1.2 磁流体基本流方程磁流体方腔槽道基本流的控制方程是一个二维偏微分方程，形式为
$$\begin{array}{c} \nabla _{2{\rm{d}}}^2\bar u - H{a^2}\left( {\bar u + \frac{{\partial \bar \phi }}{{\partial y}}} \right) = Re\frac{{\partial \bar p}}{{\partial x}},\\ \nabla _{2{\rm{d}}}^2\bar \phi = - \frac{{\partial \bar u}}{{\partial y}}. \end{array}$$磁流体方腔槽道基本流沿截面可以主要分成三个区域：核心区、Hartmann层和侧层.中心大部分区域构成核心区，在此区域内粘性力与惯性力可以忽略，压力梯度平衡电磁力，流向速度均匀分布；Hartmann层为垂直于磁场的壁面边界层，粘性力平衡压力梯度与电磁力之和；侧层(也称为Shercliff层)为平行于磁场的壁面边界层，电磁力很小，粘性力平衡压力梯度.
磁流体Hunt基本流及绝缘方腔流基本流均存在级数展开形式的精确解，也可以采用谱配置方法进行高精度的计算求解.图 2给出了Ha=10和Ha=100时的磁流体Hunt基本流的流向速度分布，可以看出大Hartmann数下，槽道中心区域的流动速度很小，而流动通量集中在两侧壁面附近，形成M型的Hunt射流.图 3给出了Hunt基本流在z=0截面上的流向速度分布，最大速度就发生在这个截面上，并且随着Hartmann数的增大，最大速度的位置越靠近侧壁面.当槽道四个壁面均为绝缘壁面时，可以计算Ha=0(无外加磁场)，Ha=10和Ha=100时的基本流速度分布，如图 4所示.无外加磁场影响时，流向速度在y和z方向上呈现近似抛物线型的分布.而随着Ha数的增大，磁场对流动的影响增强，其边界层明显变薄，同时核心区内速度趋于均匀分布，形成近似“柱塞流”(如图 4(c)所示).
 |
| 图2 磁流体Hunt基本流的流向速度分布 Fig. 2 Stream-wise velocity profiles of MHD Hunt duct flow |
 |
| 图3 磁流体Hunt基本流在z=0截面上的流向速度分布 Fig. 3 Stream-wise velocity profiles along the cross section z=0 of MHD Hunt duct flow |
 |
| 图4 不同磁场强度下绝缘方腔槽道流的基本流速度分布 Fig. 4 Stream-wise velocity profiles of MHD duct flow with perfectly insulating walls at different magnetic field strength |
磁流体方腔槽道流分解成基本流和扰动量的叠加，形式为
$$\begin{array}{*{20}{c}} {\boldsymbol{v} = \left( {\bar u,0,0} \right) + \boldsymbol{v}',}&{p = \bar p + p',}&{\phi = \bar \phi = \phi ',} \end{array}$$ 代入磁流体无量纲控制方程后，忽略二阶扰动小量后可以得到小扰动线化方程.对扰动量进行正则模形式展开 $$\left( {u',v',w',p',\phi '} \right) = \left( {\hat u,\hat v,\hat w,\hat p,\hat \phi } \right)\exp \left[{{\rm{i}}\left( {\alpha x - \omega t} \right)} \right].$$ 将正则模展开进一步代入小扰动线化方程，就可以得到整体线性稳定性方程 $$\begin{array}{c} {\rm{i}}\alpha \hat u + {\partial _y}\hat v + {\partial _z}\hat w = 0,\\ \frac{1}{{Re}}\left( {{{\rm{D}}^2} - {\alpha ^2}} \right)\hat u - {\rm{i}}\alpha \bar u\hat u - {\partial _y}\bar u\hat v - {\partial _z}\bar u\hat w - {\rm{i}}\alpha {\rm{\hat p - }}\frac{{H{a^2}}}{{Re}}\left( {\hat u + \frac{{\partial \hat \phi }}{{\partial y}}} \right) = - {\rm{i}}\omega \hat u,\\ \frac{1}{{Re}}\left( {{{\rm{D}}^2} - {k^2}} \right)\hat v - {\rm{i}}\alpha \bar u\hat v - {\partial _y}\hat p - \frac{{H{a^2}}}{{Re}}\left( {\hat v - {\rm{i}}\alpha \hat \phi } \right) = - {\rm{i}}\omega \hat v,\\ \frac{1}{{Re}}\left( {{{\rm{D}}^2} - {k^2}} \right)\hat w - {\rm{i}}\alpha \bar u\hat w - {\partial _z}\hat p = - {\rm{i}}\omega \hat w,\\ \left( {{{\rm{D}}^2} - {k^2}} \right)\hat \phi - {\rm{i}}\alpha \hat v + {\partial _y}\hat u = 0 \end{array}$$ 以及边界条件： $$\begin{array}{c} {\rm{Hunt流}}:\begin{array}{*{20}{c}} {\hat u = \hat v = \hat w = 0}&{\hat \phi = 0,}&{z = \pm 1,} \end{array}\\ \begin{array}{*{20}{c}} {\hat u = \hat v = \hat w = 0}&{\frac{{\partial \hat \phi }}{{\partial y}} = 0,}&{y = \pm 1;} \end{array} \end{array}$$ $$\begin{array}{c} {\rm{绝缘方腔流}}:\begin{array}{*{20}{c}} {\hat u = \hat v = \hat w = 0}&{\frac{{\partial \hat \phi }}{{\partial z}} = 0,}&{z = \pm 1,} \end{array}\\ \begin{array}{*{20}{c}} {\hat u = \hat v = \hat w = 0}&{\frac{{\partial \hat \phi }}{{\partial y}} = 0,}&{y = \pm 1;} \end{array} \end{array}$$ 这里，D2=∂yy+∂zz. 1.4 线性最优瞬态增长定义如果只对扰动量进行关于波数的正则模展开，形式如下
$$\left[{u',v',w',p',\phi '} \right] = \left[{\hat u,\hat v,\hat w,\hat p,\hat \phi } \right]\left( {y,z,t} \right){{\rm{e}}^{i\alpha x}},$$ 代入线化小扰动方程后，可以得到一个初值问题 $$\frac{{{\rm{d}}\boldsymbol{q}}}{{{\rm{d}}t}} = L\boldsymbol{q},$$ 这里，$\boldsymbol{q} = {\left( {\hat u,\hat v,\hat w,\hat p,\hat \phi } \right)^{\rm{T}}}$，L为线性稳定性算子.初值问题的解可以写成 $$\boldsymbol{q}\left( t \right) = {\boldsymbol{q}_0}\exp \left( {t{\rm{\boldsymbol{q}}}} \right),$$ 线性最优瞬态增长定义为 $$G\left( {t,\alpha ,Re,Ha} \right) = \mathop {\max }\limits_{{\boldsymbol{q}_0}} \left\{ {\left\| {q\left( t \right)} \right\|_{\rm{E}}^2/\left\| {{\boldsymbol{q}_0}} \right\|_{\rm{E}}^2} \right\} = \mathop {\max }\limits_{{\boldsymbol{q}_0}} \left\{ {\left\| {{\boldsymbol{q}_0}\exp \left( {t\boldsymbol{L}} \right)} \right\|_{\rm{E}}^2/\left\| {{\boldsymbol{q}_0}} \right\|_{\rm{E}}^2} \right\} = \left\| {\exp \left( {t{\rm{\boldsymbol{L}}}} \right)} \right\|_{\rm{E}}^2.$$ 通常，内积模可以定义成扰动动能形式 $$\left\| {q\left( t \right)} \right\|_{\rm{E}}^2 = \smallint \smallint \left( {\hat u\left( {y,z,t} \right){{\hat u}^ + }\left( {y,z,t} \right) + \hat v\left( {y,z,t} \right){{\hat v}^ + }\left( {y,z,t} \right) + \hat w\left( {y,z,t} \right){{\hat w}^ + }\left( {y,z,t} \right)} \right){\rm{d}}y{\rm{d}}z.$$ 这里，上标“+”表示复共轭，空间积分沿槽道截面进行.进一步地，将初值问题的解投影到N个增长率最大的特征谱上，即 $$\boldsymbol{q}\left( t \right) \approx \sum\limits_{n = 1}^N {{\boldsymbol{\kappa} _n}\left( t \right){{\widetilde {\bf{\boldsymbol{q}}}}_n},} $$ 其中，展开系数κn(t)满足 $$\frac{{\rm{d}}}{{{\rm{d}}t}}\left( {\begin{array}{*{20}{c}} {{\kappa _1}}\\ \vdots \\ {{\kappa _N}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} { - {\rm{i}}{\omega _1}}&{}&{}\\ {}& \ddots &{}\\ {}&{}&{ - {\rm{i}}{\omega _N}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{\kappa _1}}\\ \vdots \\ {{\kappa _N}} \end{array}} \right) = - {\rm{i}}\omega \kappa ,$$ 将特征谱展开代入内积模，可以得到 $$\left\| {\boldsymbol{q}\left( t \right)} \right\|_{\rm{E}}^2 = \smallint \smallint {\boldsymbol{q}^{\rm{H}}}M \boldsymbol{q}{\rm{d}}y{\rm{d}}z \approx \smallint \smallint {\left( {\sum\limits_{n = 1}^N {{\kappa _n}\left( t \right){{\tilde q}_n}} } \right)^{\rm{H}}}M\left( {\sum\limits_{m = 1}^N {{\boldsymbol{\kappa} _m}\left( t \right){{\tilde q}_m}} } \right){\rm{d}}y{\rm{d}}z = \sum\limits_{n,m = 1}^N {{\boldsymbol{\kappa} _n}{{\left( t \right)}^ + }{\boldsymbol{Q}_{mn}}{\boldsymbol{\kappa} _m}\left( t \right),} $$ 其中，正定矩阵${\boldsymbol{Q}_{mn}} = \smallint \smallint \tilde q_n^{\rm{H}}M{\tilde q_m}{\rm{d}}y{\rm{d}}z$可以进行Cholesky分解Q=FFH.由此，内积模可以近似成L2模的形式 $$\left\| {\boldsymbol{q}\left( t \right)} \right\|_{\rm{E}}^2 \approx {\boldsymbol{\kappa} ^{\rm{H}}}\boldsymbol{Q\kappa} = {\boldsymbol{\kappa} ^{\rm{H}}}{\boldsymbol{F}^{\rm{H}}}\boldsymbol{F\kappa} = {\left\langle {\boldsymbol{F\kappa} ,\boldsymbol{F\kappa} } \right\rangle _2} = \left\| {\boldsymbol{F\kappa} } \right\|_2^2,$$ 线性最优瞬态增长可以简化成 $$G\left( {t,\alpha ,Re,Ha} \right) = \mathop {\max }\limits_{{\boldsymbol{\kappa} _0}} \left\{ {\left\| {\boldsymbol{F\kappa} \left( t \right)} \right\|_2^2/\left\| {\boldsymbol{F}{\boldsymbol{\kappa} _0}} \right\|_2^2} \right\} = \mathop {\max }\limits_{{\boldsymbol{\kappa} _0}} \left\{ {\left\| {F\exp \left( { - i\omega t} \right)} \right\|_2^2/\left\| {\boldsymbol{F}{\boldsymbol{\kappa} _0}} \right\|_2^2} \right\}.$$ 这样，线性最优瞬态增长又可以转化成求解 $$G\left( {t,\alpha ,Re,Ha} \right) = \left\| {\boldsymbol{F}\exp \left( { - {\rm{i}}\boldsymbol{\omega} t} \right){\boldsymbol{F}^{ - 1}}} \right\|_2^2,$$ 其中，ω是由最不稳定的N个整体线性稳定性模态特征值构成的对角阵，而L2模的计算可以采用奇异值分解的方式进行[16, 17]. 2 磁流体方腔槽道流整体线性稳定性数值方法整体线性稳定性计算方法可以分为时间步进方法[18](time-stepping approach)和广义特征值方法[19](generalized eigenvalue problem)，如图 5左图所示.时间步进方法采用高精度数值方法模拟线化的Navier-Stokes方程，考察任意小扰动的时间发展，由此采用特征值迭代方法计算整体最不稳定的模态；该方法对流动区域形状没有限制，但是它一般只能考虑有限的几个最不稳定的特征模态.广义特征值方法是通过Laplace变换将时域变换为频域，构成了一个大型矩阵形式的广义特征值问题，然后计算其整体特征模态.广义特征值方法可以考虑大量的流场特征模态，但却通常应用于简单的流动几何区域.
 |
| 图5 流动整体线性稳定性分析的计算方法流程图 Fig. 5 Numerical algorithm flow charts for linear global stability analysis of fluid flows |
基于方腔槽道流的简单几何流动区域，它的整体线性稳定性研究一般采用广义特征值方法，其核心算法在于整体线性稳定性方程(二维偏微分方程组)的空间离散和大型广义矩阵特征值的计算.从现有空间离散方法来看(图 5右图)，主要分成低阶精度格式和高阶精度格式.低阶精度格式包括有限差分、有限元和有限体积，需要更密的网格划分以实现计算收敛，因而整体线性稳定性分析时具有很大的计算量.而高阶精度格式有谱配置法、谱元法和高精度有限差分法，这些方法是当前整体线性稳定性研究普遍采用的方法.这里将主要介绍Chebyshev谱配置法和基于非均匀网格的高阶FD-q差分方法.
2.1 Chebyshev谱配置法首先，定义在区间[-1, 1]上的n阶Chebyshev多项式为
$${T_n}\left( \xi \right) = \cos \left( {n\arccos \xi } \right).$$ 根据代数精度的需要，我们将配置点取在N阶Chebyshev多项式的极值点(通常称为Chebyshev-Gauss-Lobatto谱配置点)上： $${\xi _j} = \cos \left( {j{\rm{\pi }}/N} \right),j = 0,1,\cdots ,N.$$ 然后，根据磁流体方腔槽道流的整体线性稳定性方程是一个二维偏微分方程的特征值问题，可以对待求特征函数以Chebyshev多项式为基函数进行展开，以电势特征函数为例 $$\hat \phi \left( {y,z} \right) = \sum\limits_{n = 0}^{{N_y}} {\sum\limits_{m = 0}^{{N_z}} {{a_{nm}}{T_n}\left( \eta \right){T_m}\left( \zeta \right),} } $$ 可以对展开式进行各阶偏导数计算，如一阶导数 $$\begin{array}{l} \frac{{\partial \hat \phi \left( {y,z} \right)}}{{\partial y}} = \left( {\frac{{{\rm{d}}\eta }}{{{\rm{d}}y}}} \right)\sum\limits_{n = 0}^{{N_y}} {\sum\limits_{m = 0}^{{N_z}} {{a_{nm}}{{T'}_n}\left( \eta \right){T_m}\left( \zeta \right),} } \\ \frac{{\partial \hat \phi \left( {y,z} \right)}}{{\partial y}} = \left( {\frac{{{\rm{d}}\zeta }}{{{\rm{d}}z}}} \right)\sum\limits_{n = 0}^{{N_y}} {\sum\limits_{m = 0}^{{N_z}} {{a_{nm}}{T_n}\left( \eta \right){{T'}_m}\left( \zeta \right),} } \end{array}$$ 根据链式法则，还可以给出高阶导数的表达式.为了计算配置点上的这些偏导数，关键是要计算Chebyshev多项式的各阶导数，这可以通过Chebyshev多项式低阶与高阶的递推关系计算，即 $$\begin{array}{l} \begin{array}{*{20}{c}} {T_0^{\left( k \right)}\left( \xi \right) = 0,}&{T_1^{\left( k \right)}\left( \xi \right) = T_0^{\left( {k - 1} \right)}\left( \xi \right)} \end{array},\\ T_0^{\left( k \right)}\left( \xi \right) = 2nT_{n - 1}^{\left( {k - 1} \right)}\left( \xi \right) + \frac{n}{{n - 2}}T_{n - 2}^{\left( k \right)}\left( \xi \right). \end{array}$$ 最后，将展开函数及其高阶导数表达式代入整体线性稳定性方程，在内部配置点上满足整体线性稳定性方程，而在边界配置点上用边界条件来替代，这样就可以构造一个大型的广义特征值矩阵 $$\boldsymbol{Ax} = \boldsymbol{\omega Bx}$$对于原始变量的磁流体方腔槽道流的整体线性稳定性方程，大型矩阵A和B的存储规模为2×16×52×N4×10-9Gbytes.这里，2表示大型矩阵个数，16表示一个双精度复数所需字节数，5为整体线性稳定性方程个数，N为单个方向的配置点个数.磁流体方腔槽道流属于典型的较大雷诺数下的边界层问题，因此对于大Hartmann数时N>50才能达到足够的精度，这样广义特征值矩阵的存储规模将大于5Gbytes.为了减小大型矩阵的存储规模，由此减小计算量，可以通过消去某些原始变量来减小整体线性稳定性方程的个数.本文的计算将采用消去流向速度特征函数(利用连续性方程)，将磁流体方腔槽道流的整体线性稳定性方程由5个减小为4个.
2.2 高阶FD-q差分格式高阶FD-q差分格式最先由Hermanns和Hernandez[20]提出，Paredes[15]将其运用于整体线性稳定性特征模态的计算.FD-q差分格式与其它主要差分格式的不同之处在于采用了非均匀网格，Hermanns和Hernandez[20]建立的非均匀网格采用一种迭代算法构造.实际上，非均匀网格可以不必采用这种复杂算法构造，可以直接采用Chebyshev谱配置点来作为非均匀网格点，这样的构造便于采用Chebyshev-Gauss-Lobatto积分公式进行高精度的数值积分，从而有利于线性最优瞬态增长的计算.如果为了进一步提高Cheyshev谱配置点的非均匀网格在边界上的分辨率，还可以采用Kosloff-Tal-Ezer变换[21].这样，方腔截面[-1, 1]×[-1, 1]上的配置点可以表示成
$$\begin{array}{l} {y_i} = \arcsin \left( { - {\alpha _0}\cos \left( {{\rm{\pi }}i/N} \right)} \right)/\arcsin {\alpha _0},i = 0,1,\cdots ,N,\\ {z_j} = \arcsin \left( { - {\alpha _0}\cos \left( {{\rm{\pi }}j/N} \right)} \right)/\arcsin {\alpha _0},j = 0,1,\cdots ,N, \end{array}$$ 这里，参数α0为正数，用来控制网格点的伸缩.沿某个坐标方向(如y方向)，计算区域可以被分成N+1子区间Ωi(如图 6所示)，定义如下 $${\Omega _i} = \left[{\left( {{y_{i - 1}} + {y_i}} \right)/2,\left( {{y_i} + {y_{i + 1}}} \right)/2} \right],i = 1,\cdots ,N - 1.$$ 每个区间Ωi包含相应的配置点，然后在每个区间Ωi上进行q阶拉格朗日插值Ii(y)， $$\begin{array}{*{20}{c}} {{I_i}\left( y \right) = \sum\limits_{j = {s_i}}^{{s_i} + q} {{l_{ij}}\left( y \right){\phi _j},} }&{{l_{ij}}\left( y \right) = \prod\limits_{m = 0,j \ne {s_i} + m}^q {\left( {y - {y_{{s_i} + m}}} \right)/\left( {{y_j} - {y_{{s_i} + m}}} \right)} } \end{array}.$$ 进一步，q阶拉格朗日插值可以写成重心插值(barycentric interpolation)， $$\begin{array}{*{20}{c}} {{I_i}\left( y \right) = \sum\limits_{j = {s_i}}^{{s_i} + q} {\frac{{{w_j}}}{{y - {y_j}}}{\phi _j}} /\sum\limits_{j = {s_i}}^{{s_i} + q} {\frac{{{w_j}}}{{y - {y_j}}}} }&{{w_j} = \prod\limits_{m = 0,j \ne {s_i} + m}^q {{{\left( {y - {y_{{s_i} + m}}} \right)}^{ - 1}}} } \end{array}.$$ 需要强调的是，拉格朗日重心插值已经被证明具有更好的计算效率和数值稳定性[27, 28].当q为偶数阶时,插值模板 $$\left\{ {{s_i}} \right\} = {\rm{\{ }}\underbrace {0,\cdots ,0}_{q/2\;\;{\rm{times}}},\underbrace {0,1,\cdots ,N - q,}_{{\rm{centered}}\;{\rm{FD}}}\underbrace {N - q,\cdots ,N - q}_{q/2\;\;{\rm{times}}}{\rm{\} }},$$ 而当q为奇数阶时,插值模板 $$\left\{ {{s_i}} \right\} = {\rm{\{ }}\underbrace {0,\cdots ,0}_{\left( {q - 1} \right)/2\;\;{\rm{times}}},\underbrace {0,1,\cdots ,N - q,}_{{\rm{centered}}\;{\rm{FD}}}\underbrace {N - q,\cdots ,N - q}_{\left( {q + 1} \right)/2\;\;{\rm{times}}}{\rm{\} }},$$ 同样，FD-q差分格式还可以在z方向进行，并且还可以在每个方向上进行求导，以实现高阶导数的空间离散. |
| 图6 FD-q差分格式的非均匀网格分布(由Chebyshev谱配置点构成),q=6和N=10 Fig. 6 Non-uniform grid point distributions (consist of Chebyshev collocation points) of FD-q finite difference scheme with q=6 and N=10 |
采用Chebyshev谱配置法或高阶FD-q差分格式对磁流体方腔槽道流动整体线性稳定性方程进行空间离散后，最终转化为广义矩阵特征值问题(Ax=ωBx)的求解，而这类大型的广义特征值问题往往采用Arnoldi迭代算法.首先采用shift-and-invert策略(参见Theofilis[22]),将广义矩阵特征值问题化为一般矩阵特征值问题((A-μB)-1Bx=θx)，然后对(A-μB)采用LU分解求逆，可以采用UMFPACK[23]或MUMPS[24]软件包，最后使用Lehoucq等人开发的并行版本ARPACK软件包[25](隐式重启Arnoldi迭代方法)计算一般矩阵的特征值.
3 数值结果验证与分析 3.1 磁流体Hunt槽道流的整体线性稳定性在磁流体Hunt槽道流中，与磁场垂直的壁面是导电的，而与磁场平行的壁面是绝缘的.磁流体Hunt槽道流在绝缘壁面附近出现Hunt射流，磁场强度对速度分布具有显著影响，Hartmann数越大射流越靠近壁面，甚至出现速度拐点，出现具有KH不稳定的速度分布.
表 1给出了磁流体方腔Hunt槽道流在参数α=0.8,Ha=10,Re=2 000下的不同配置点个数计算得到整体线性最不稳定特征谱，而表 2给出的是在参数α=5,Ha=100,Re=1 000下线性最不稳定特征谱.比较表 1和表 2，可以看出大Hartmann数需要更多的配置点数来达到相同的收敛精度.这主要是边界层效应引起的，两侧壁附近的Shercliff边界层厚度正比于Ha-1/2，同时上下壁面附近的Hartmann边界层厚度更正比Ha-1.进一步，可以清楚地发现表 2的第二组特征谱对应于Priede等人[12]给出的最不稳定特征谱，而实际上真正的最不稳定模态为表 2的第一组特征谱.
| 表1 谱配置法计算得到的Hunt槽道流整体线性最不稳定特征谱 (α=0.8, Ha=10, Re=2 000) Table 1 Least stable spectrum of linear global stability for MHD Hunt duct flows by spectral collocation method (α=0.8, Ha=10, Re=2 000) |
 |
| 点击放大 |
| 表2 谱配置法计算得到的Hunt槽道流整体线性最不稳定特征谱 (α=5, Ha=100, Re=1 000) Table 2 Least stable spectrum of linear global stability for MHD Hunt duct flows by spectral collocation method (α=5, Ha=100, Re=1 000). |
 |
| 点击放大 |
表 3给出了采用FD-q差分格式计算得到的在参数α=5, Ha=100,Re=1 000下Hunt槽道流的整体线性最不稳定特征谱，与表 2对照表明：相比谱配置方法，达到相同的收敛精度FD-q差分格式需要更多的非均匀网格点，但是采用稀疏矩阵存储技术使得计算效率得到了显著的提高[15].
| 表3 FD-q差分格式计算得到的Hunt槽道流整体线性最不稳定特征谱 (α=5, Ha=100, Re=1 000) Table 3 Least stable spectrum of linear global stability for MHD Hunt duct flows by FD-q finite difference scheme (α=5, Ha=100, Re=1 000) |
 |
| 点击放大 |
图 7给出了低磁场作用下(Ha=10)磁流体Hunt流的整体线性稳定性特征.左图为稳定性中性边界曲线，从曲线中可以确定临界雷诺数为2 018.2和临界波数为0.83；右图为流向速度扰动的三维特征时空结构(Re=2 500和α=0.8).从图中可以看出，低磁场作用下(Ha=10)磁流体Hunt流在中等雷诺数下最先发生失稳，最不稳定模态发生在小波数区域，因此发生的是流体动力学长波不稳定；流向速度扰动的三维特征时空结构偏向侧壁面，沿方腔截面成左右和上下反对称特征结构.
 |
|
图7 低磁场作用下(Ha=10)磁流体Hunt流的整体线性稳定性 (左：中性边界曲线；右：流向速度扰动的三维特征时空结构，Re=2 500，α=0.8) Fig. 7 Linear global stability of MHD Hunt duct flows under action of weak magnetic field(Ha=10) (left: neutral boundary curve; right: three-dimensional spatial-temporal eigenstructure of streamwise velocity perturbation at Re=2 500 and α=0.8) |
图 8给出了中强磁场作用下(Ha=100)磁流体Hunt流的整体线性稳定性特征.左图为稳定性中性边界曲线，从曲线中可以确定临界雷诺数为1 170和临界波数为4.9；右图为流向速度扰动的三维特征时空结构(Re=2 500和α=4.0).从图中可以看出，中强磁场作用下磁流体Hunt流可以在更小的雷诺数下发生线性失稳，并且具有两个非常靠近的失稳模态；最不稳定模态发生在较大的波数区域，因此发生的是流体动力学短波不稳定；流向速度扰动的三维特征时空结构更加靠近侧壁面，沿方腔截面成上下对称结构，而两个失稳模态分别对应左右对称和反对称结构(图中给出的是反对称的时空结构).
 |
|
图8 中强磁场作用下(Ha=100)磁流体Hunt流的整体线性稳定性 (左：中性边界曲线；右：流向速度扰动的三维特征时空结构，Re=2 500，α=4.0) Fig. 8 Linear global stability of MHD Hunt duct flows under the action of moderately strong magnetic field(Ha=100) (left: neutral boundary curves; right: three-dimensional spatial-temporal eigenstructure of streamwise velocity perturbation at Re=2 500 and α=4.0) |
对于壁面完全绝缘的磁流体槽道流，磁场对速度分布改变相对导电壁面要小，一般不会产生射流速度分布，流动在非常高的雷诺数下才可能线性失稳.因此，对其主要应该开展线性最优瞬态增长研究.为了验证高阶FD-q差分格式用于计算线性最优瞬态增长的准确性，我们首先计算了无磁场作用下三种流向波数(α =0,0.5,1.0)下的线性最优瞬态增长，如图 9所示.采用的流动雷诺数与Biau等人[26]一致，计算得到的线性最优瞬态增长曲线与Biau等人的结果一致，波数为零时发生最大的最优瞬态增长，同时还计算得到了一致的初始最优截面速度扰动空间结构.
 |
|
图9 无磁场作用下槽道流的线性最优瞬态增长(Ha=0, Re=3 315.2) (左：不同时刻下的最优增长；右：波数为零时初始最优截面速度扰动的空间结构) Fig. 9 Linear optimal transient growths of duct flows without magnetic field with Ha=0 and Re=3 315.2 (left: optimal growths as functions of time at different wave numbers; right: spatial structure of initial optimal velocity perturbations along the cross section) |
图 10给出了不同波数和不同Hartmann数下壁面完全绝缘磁流体槽道流的线性最优瞬态增长情况.从图中可以看出最大最优瞬态增长并不发生在波数为零的情况下，同时磁场对最优瞬态增长有抑制作用.进一步地，图 11给出了波数为零时壁面完全绝缘磁流体槽道流的线性最优瞬态增长的初始最优截面速度扰动分布，从图中可以看出初始最优截面速度扰动成涡状结构，位于侧壁面附近，并且随着磁场增大(从Ha=10到Ha=50)，涡状结构尺度变小，也更加贴近侧壁面.这些结果与Krasnov[14]采用伴随优化方法得到的结果相一致.
 |
| 图10 不同波数和不同Hartmann数下壁面完全绝缘磁流体槽道流的线性最优瞬态增长 Fig. 10 Linear optimal transient growths of MHD duct flows with perfectly insulating walls at different wave numbers and Hartmann numbers |
 |
| 图11 壁面完全绝缘磁流体槽道流的线性最优瞬态增长的初始最优截面速度扰动空间结构(α=0) Fig. 11 Spatial structure of initial optimal velocity perturbations along the cross section for linear optimal transient growths of MHD duct flows with perfectly insulating walls |
考察了两种适用于磁流体方腔槽道流整体线性稳定性研究的数值计算方法：Chebyshev谱配置法和基于非均匀网格的高阶FD-q差分格式.Chebyshev谱配置法收敛快且精度高，广泛应用于方腔槽道流动的整体线性稳定性计算，但是其需要构造非常庞大的满矩阵，存储量和计算开销巨大；高阶FD-q差分格式采用的是非均匀网格，可以保证其网格精度收敛远好于均匀网格的高阶有限差分格式，同时由于差分格式可以采用稀疏矩阵进行存储，显著降低了存储量和计算开销，由此使得高阶FD-q差分格式也非常适合方腔槽道流动的整体线性稳定性研究.本文采用的非均匀网格由基于Kosloff-Tal-Ezer变换的Cheyshev谱配置点构成，既可以提高在边界上的网格分辨率，又便于进行高精度的数值积分，由此可以用于方腔槽道流动线性最优瞬态增长的计算.此外，采用拉格朗日重心插值技术提高了计算效率和数值稳定性.
为了验证两种数值计算方法的精确性和有效性，采用Chebyshev谱配置法实现了磁流体Hunt流整体线性稳定性特征谱的计算，给出了低磁场(Ha=10)和中强磁场(Ha=100)作用下的线性稳定性边界曲线，与Priede等人[12]的计算结果相一致，验证了磁流体Hunt方腔槽道流失稳可以发生在较低的雷诺数下，层流至湍流转捩可以是从线性失稳开始.同时，采用高阶FD-q差分格式实现了绝缘壁面磁流体槽道流的线性最优瞬态增长，揭示了不同波数和不同Hartmann数下的线性最优瞬态增长特征，计算结果与Krasnov等人[14]采用伴随优化方法得到的结果一致.
| [1] | DAVIDSON P A. Magnetohydrodynamics in materials processing[J]. Annu Rev Fluid Mech, 1999, 31:273-300. |
| [2] | THOMAS B G, ZHANG L. Mathematical modeling of fluid flow in continuous casting:a review[J]. ISIJ Intl, 2001,41:1181-1193. |
| [3] | BÜHLER L. Liquid metal magnetohydrodynamics for fusion blankets[M]//MOLOKOV S, MOREAU R, MOFFATT H K, ed. Magnetohydrodynamics:historical evolution and trends, Springer, 2007:171-194. |
| [4] | 倪明玖.磁约束核聚变关键能量转换部件的磁流体力学研究[J].力学与实践,2010, 32(5):1-9. |
| [5] | KAKUTANI T. The hydrodynamic stability of the modified plane Couette flow in the presence of a transverse magnetic field[J]. J Phys Soc Japan, 1964, 19:1041-1057. |
| [6] | LEHNERT B. On the behaviour of an electrically conductive liquid in a magnetic field[J]. Ark Fys,1952, 5:69-90. |
| [7] | HUNT J C R. Magnetohydrodynamic flow in rectangular ducts[J]. J Fluid Mech, 1965, 21:577-590. |
| [8] | MOFFATT H K.On the suppression of turbulence by a uniform magnetic field[J]. J Fluid Mech, 1967, 28:571-592. |
| [9] | DAVIDSON P A. Magnetic damping of jets and vortices[J]. J Fluid Mech, 1995, 299:153-186. |
| [10] | SHATROV V, GERBETH G. Marginal turbulent magnetohydrodynamic flow in a square duct[J]. Physics of Fluids, 2010, 22:084101. |
| [11] | KINET M, KNAEPEN B, MOLOKOV S. Instabilities and transition in magnetohydrodynamic flows in ducts with electrically conducting walls[J]. Physical Review Letters, 2009, 103(15):154501. |
| [12] | PRIEDE J, ALEKSANDROVA S, MOLOKOV S. Linear stability of Hunt's flow[J]. Journal of Fluid Mechanics, 2010, 649:115. |
| [13] | PRIEDE J, ALEKSANDROVA, MOLOKOV S. Linear stability of magnetohydrodynamic flow in a perfectly conducting rectangular duct[J]. Journal of Fluid Mechanics, 2012, 708:111-127. |
| [14] | KRASNOV D, ZIKANOV O, ROSSI M, et al. Optimal linear growth in magnetohydrodynamic duct flow[J]. Journal of Fluid Mechanics, 2010, 653:273-299. |
| [15] | PAREDES P, HERMANNS M, CLAINCHE S L, THEOFILIS V. Order 104 speedup in global linear instability analysis using matrix formation[J]. Comput Methods Appl Mech Engrg, 2013, 253:287-304. |
| [16] | SCHMID P J, HENNINGSON D S. Stability and transition in shear flows[M]. New York:Springer-Verlag,2001,142. |
| [17] | SCHMID P J. Nonmodal Stability Theory[J]. Ann Rev Fluid Mech, 2007, 39:129-162. |
| [18] | TUCKERMAN L S, BARKLEY D. Bifurcation analysis for time-steppers[M]//DOEDEL E,TUCKERMAN L S, ed. Numerical methods for bifurcation problems and large-scale dynamical systems, New York:Springer, 2000. |
| [19] | THEOFILIS V. Advances in global linear instability analysis of nonparallel and three-dimensional flows[J]. Prog Aero Sci, 2003, 39:249-315. |
| [20] | HERMANNS M, HERNÀNDEZ J A. Stable high-order finite-difference methods based on non-uniform grid point distributions[J]. Int J Numer Meth Fluids, 2008, 56:233-255. |
| [21] | KOSLOFF D, TAL-EZER H. A modified Chebyshev pseudospectral method with an O(N-1) time step restriction[J]. J Comput Phys, 1993, 104(2):457-469. |
| [22] | THEOFILIS V. Global linear instability[J]. Annu Rev Fluid Mech, 2011, 43:319-352. |
| [23] | DAVIS T A. Algorithm 832:UMFPACK, an unsymmetric-pattern multifrontal method[J]. ACM Transactions on Mathematical Software, 2004,30(2):196-199. |
| [24] | AMESTORY P, DUFF I, KOSTER J. A fully asynchronous multifrontal solver using distributed dynamic scheduling[J]. SIAM J Matrix Anal, 2001, A 1:15-41. |
| [25] | LEHOUCQ R B, SORENSEN D C, YANG C. ARPACK users guide:Solution of large-scale eigenvalue problems with implicitly restarted Arnoldi methods[M]. Philadelphia:SIAM,1998. |
| [26] | BIAU D, SOUEIDA H, BPTTARO A. Transition to turbulence in duct flow[J]. J Fluid Mech, 2008, 596:133-142. |
| [27] | BERRUT J P, TREFETHEN LN. Barycentric Lagrange interpolation[J]. SIAM Review, 2004, 46:501-517. |
| [28] | HIGHAM N. The numerical stability of barycentric Lagrange interpolation[J]. IMA Journal of Numerical Analysis, 2004, 24:547-556. |
引用本文 |
