1.   引　言

作为超材料研究设计的一大理论基础, 变换理论提升了我们对于自然现象的操控能力. 变换光学^[1-5]的提出为我们开启了通往变换理论的大门. 自基于麦克斯韦方程形式经坐标变换后的不变性产生的变换光学^[1-5], 到声学领域中的变换声学^[6-10], 再到热学领域的变换热学^[11-20], 变换理论不断被丰富.

为了将变换理论拓展到流动控制领域, 有学者通过将坐标变换理论与Brinkman-Stokes方程、达西定律相结合的方式设计出具有各向异性渗透率的流动隐身斗篷^[21,22], 实现了物体在多孔介质流场中的隐身. 为推动变换理论在非多孔介质中的应用, 研究人员^[23,24]证明得出斯托克斯方程满足坐标变换后的形式不变性, 进而建立了变换理论一个新的分支—变换流体动力学. 目前已有学者依据变换流体动力学, 通过调控流体动力黏度实现了蠕动流下的物体隐身^[23]和复杂形状物体的流动隐身^[25]; 通过张量化流体动力黏度, 得到流体动力黏度的空间分布, 在引导流体流动的同时避免流体与物体之间的接触, 达到物体隐身与减小阻力的效果^[23]. 为克服上述各向异性流动超材料在实验实现上的困难, 部分学者利用散射抵消理论^[26]、对流扩散抵消理论^[27]、外加电场^[28]等方式实现了各向同性流动超材料下的物体隐身效果. 而随着神经网络学科的发展, 也有学者通过深度强化学习^[29]的方法实现流动隐身的效果. 当然, 变换流体动力学不仅可以应用于流动隐身超材料的设计, 还可为流动旋转^[30]与流动聚集^[31]装置的设计提供理论支持.

与电磁聚集^[32,33]、热流聚集^[12,34]及水波聚集^[35]异曲同工, 流动聚集^[31,36]可将流体汇聚于指定区域, 实现放大流速的效果. 现阶段关于流动聚集的研究局限于汇集单一方向的流体, 这一点极大地限制了流动聚集器的实用功能. 倘若可以在聚集流体的同时, 改变流体在聚集域内的流向, 即任意调控汇聚流体的速度方向, 那么流动聚集器的适用范围将会大大提高. 基于上述想法, 本文耦合流动旋转与聚集功能, 依托变换流体动力学, 推导流体动力黏度在指定流域的空间分布, 完成文丘里效应旋聚器(下文简称旋聚器)的理论设计. 通过数值模拟验证得出: 本文所设计文丘里效应旋聚器可在蠕动流状态下同时实现对流体的聚集和旋转功能, 打破传统流动聚集器仅可汇集单一方向流体的限制. 本研究可为流动旋聚器相关实验奠定理论基础；可为磁滞现象的相关研究拓展至旋转滞后现象带来启发；此外, 也可为超材料的非互易坐标变换设计提供新的思路.

2.   文丘里效应旋聚器的理论设计

以原始空间${\boldsymbol{x}}(x, y, z)$和物理空间${\boldsymbol{x}} ^{\prime}(x', y', z')$分别表示坐标变换前后的空间, 则对于稳态、不可压缩、无体积力的蠕动流, 控制牛顿流体流动的连续性方程和动量守恒方程可表示为

$$\nabla \cdot {\boldsymbol{u}} = 0 ,$$

(1)

$$\mu \nabla^{2} {\boldsymbol{u}} = \nabla {p} ,$$

(2)

式中, ${\boldsymbol{u}}$, $\mu$, $p$分别表示速度、动力黏度与压强.

根据变换流体动力学^[23,24], 方程(1)和(2)在坐标变换前后方程形式保持一致, 进而可推导出坐标变换后控制牛顿流体流动的连续性方程和动量守恒方程为

$$\nabla' \cdot {\boldsymbol{u}}' = 0 ,$$

(3)

$$\begin{array}{l} {\boldsymbol{\mu}}' \nabla^{\prime 2} {\boldsymbol{u}}' = \nabla' p' , \end{array}$$

(4)

式中, ${\boldsymbol{u}}'$, ${\boldsymbol{\mu}}'$, $p'$分别表示流体在物理空间下的速度、动力黏度张量与压强. 其中, ${\boldsymbol{u}}^{\prime}= {\boldsymbol{J}}\cdot{\boldsymbol{u}} / {\rm{det}}({\boldsymbol{J}})$.

$$\begin{array}{l} {\boldsymbol{\mu}}^{\prime} = \operatorname {det}({\boldsymbol{J}})\cdot ({\boldsymbol{J}})^{ - 1}\mu {\boldsymbol{J}}^{ - {\rm{T}}}, \end{array}$$

(5)

$$\begin{split} {\boldsymbol{J}} = \;&\frac{\partial\left(x^{\prime}, y^{\prime}, z^{\prime}\right)}{\partial(x, y, z)} = \ {\boldsymbol{J}}_{x^{\prime} r^{\prime}} {\boldsymbol{J}}_{r^{\prime} r} {\boldsymbol{J}}_{r x} \\ =\;& \frac{\partial\left(x^{\prime}, y^{\prime}, z^{\prime}\right)}{\partial\left(r^{\prime}, \theta^{\prime}, z^{\prime}\right)} \frac{\partial\left(r^{\prime}, \theta^{\prime}, z^{\prime}\right)}{\partial(r, \theta, z)} \frac{\partial(r, \theta, z)}{\partial(x, y, z)}. \end{split}$$

(6)

流体动力黏度张量表达式如(5)式所示, 式中${\boldsymbol{J}}$为雅可比空间变换矩阵, 其表达式如(6)式所示. ${\boldsymbol{J}}^{-1}$和${\boldsymbol{J}}^{-{\rm{T}}}$分别为雅可比空间变换矩阵的逆矩阵及其逆矩阵的转置. 变换后的动力黏度张量连通了原始空间与物理空间. 接下来我们将根据旋聚器的功能设计要求, 具体化空间之间的联系.

图 1(a)所示为旋聚器的坐标变换原理. 为实现旋转聚集功能, 将原始空间($0 < r \leqslant R_2$)区域处坐标空间压缩至($0 < r \leqslant R_1$), 将($R_2 < r \leqslant$$R_3$)区域空间拉伸至($R_1 < r \leqslant R_3$); 接着我们将原始空间($0 < r \leqslant R_1$)区域处坐标空间旋转$\theta_0$度, 将($R_1 <$$r \leqslant R_3$)区域空间自$R_3$处的0度逐渐旋转至$R_1$处的$\theta_0$度. 相应地, 图 1(a)中坐标变换所反映的空间映射关系如(7)式所示. 旋转角度$\theta_0$可根据实际需要预先给定, 半径$R_1, R_2, R_3$如图 1(b)所示. 结合(6)式中雅可比空间变换矩阵的一般形式与(7)式中的映射关系, 可进一步得出反映这种映射关系的雅可比空间变换矩阵, 如(8)式所示.

$$\left\{\begin{aligned} &r^{\prime} = \frac{R_3-R_1}{R_3-R_2} r- \frac{R_2-R_1}{R_3-R_2} R_3, \\ &\theta^{\prime} = \theta-\theta_{0} \frac{R_3-r}{R_3-R_1}, \\ &z^{\prime} = z. \end{aligned}\right.$$

(7)

$$\boldsymbol{J} = \left[\begin{array}{ccc} \dfrac{\partial x^{\prime}}{\partial r^{\prime}} & \dfrac{\partial x^{\prime}}{r^{\prime} \partial \theta^{\prime}} & \dfrac{\partial x^{\prime}}{\partial z^{\prime}}\\ \dfrac{\partial y^{\prime}}{\partial r^{\prime}} & \dfrac{\partial y^{\prime}}{r^{\prime} \partial \theta^{\prime}} & \dfrac{\partial y^{\prime}}{\partial z^{\prime}}\\ \dfrac{\partial z^{\prime}}{\partial r^{\prime}} & \dfrac{\partial z^{\prime}}{r^{\prime} \partial \theta^{\prime}} & \dfrac{\partial z^{\prime}}{\partial z^{\prime}} \end{array}\right]\left[\begin{array}{ccc} \dfrac{R_3-R_1}{R_3-R_2} & 0 & 0\\ \dfrac{r^{\prime} \theta_{0}}{R_3-R_1} & \dfrac{r^{\prime}}{r} & 0\\ 0 & 0 & 1 \end{array}\right]\left[\begin{array}{ccc} \dfrac{\partial r}{\partial x} & \dfrac{\partial r}{\partial y} & \dfrac{\partial r}{\partial z}\\ \dfrac{r \partial \theta}{\partial x} & \dfrac{r \partial \theta}{\partial y} & \dfrac{r \partial \theta}{\partial z}\\ \dfrac{\partial z}{\partial x} & \dfrac{\partial z}{\partial y} & \dfrac{\partial z}{\partial z} \end{array}\right].$$

(8)

通过分析(8)式发现, 雅可比空间变换矩阵中包含旋转矩阵 $\boldsymbol{R}\left(\theta^{\prime}\right)$ 和$\boldsymbol{R}^{\rm{T}}\left(\theta \right)$. 前者为物理空间下的旋转矩阵, 后者为原始空间下的旋转矩阵. 从而可将(8)式进一步简化为(9)式。

$$\boldsymbol{J} = \boldsymbol{R}\left(\theta^{\prime}\right)\left[\begin{array}{ccc} \dfrac{R_3-R_1}{R_3-R_2} & 0 & 0\\ \dfrac{r^{\prime} \theta_{0}}{R_3-R_1} & \dfrac{r^{\prime}}{r} & 0\\ 0 & 0 & 1 \end{array}\right]\boldsymbol{R}^{\mathrm{T}}(\theta),$$

(9)

其中旋转矩阵的具体形式为

$$\begin{array}{l} \begin{aligned} \boldsymbol{R}\left({\theta}^{\prime}\right) = \left[\begin{array}{ccc} \cos\theta^{\prime} & -\sin\theta^{\prime} & 0\\ \sin\theta^{\prime} & \cos\theta^{\prime} & 0\\ 0 & 0 & 1\\ \end{array}\right] \end{aligned} \end{array}.$$

(10)

结合(5)式, (7)式和(9)式, 可得设计旋聚器所需的流体动力黏度在柱坐标下的分布规律为

$$\begin{array}{l} \begin{aligned} {\boldsymbol{\mu}}^{\prime} & = \left[\begin{array}{ccc} A & B & 0 \\ C & D & 0 \\ 0 & 0 & I \end{array}\right] \mu \end{aligned} \end{array},$$

(11)

其中

$$\begin{split} A = \;&\frac{{\theta_{0}^2}{(R_3 - R_2)^2}r^{\prime}} {(R_2 - R_1)^2(R_3 - R_1)^2}\big [r^{\prime}+\frac{R_3(R_2 - R_1)}{R_3 - R_2}\big ] \\& + \frac{(R_3 - R_2){r^{\prime}}}{(R_3 - R_2){r}+R_3(R_2 - R_1)},\\ B = \;&C = \frac{\theta_{0} \big[(R_3 - R_2)r^{\prime} + (R_2 - R_1)R_3\big]} {{(R_3-R_1)}^{2}},\\ D =\;& 1 + \frac{(R_2 - R_1)R_3} {(R_3 - R_2)r^{\prime}} ,\\ I =\;& \frac{(R_3 - R_1)r^{\prime}} {(R_3 - R_2)r^{\prime} + (R_2 - R_1)R_3} \frac{R_3 - R_1} {R_3 - R_2}. \end{split}$$

观察方程(11)可以发现, 流体动力黏度在物理空间下为一张量, 呈现各向异性的特点. 为实现该种情况下流体的流动效果, 参考等效介质理论^[37]和流动超材料文献^[23]设计的微米柱方法, 在流域中布置具有一定规律的微米柱阵列, 获得流体在各向异性动力黏度下的等效流场. 至此, 基于变换流体动力学设计的旋聚器所需参数均已获取. 接下来将在仿真软件COMSOL Multiphysics 中验证旋聚器工作效果.

3.   数值模拟及结果分析

设定流域整体尺寸$L \times H \times D$分别为长10 cm, 高10 cm, 深50 µm, 其中深度方向为流域在$z$轴方向上的高度. 考虑到流域深度$D \ll L$, $H$, 来流为Hele-Shaw 流动^[38,39], 即可视为二维流动. 如图 1(b)所示, 流体沿$x$方向自左向右流动, 进出口压差$\Delta p = 2$ kPa $(p_1 > p_2)$. $y$方向上下两边界为无滑移壁面. 背景流域流体为293.15 K下的液态水, 该状态下水的动力黏度$\mu = 10^{ - 3}$ Pa·s, 密度$\rho =$$997.1$ kg/m³.

下面根据图(2)所示的速度场分布图来比较传统流动聚集^[30]、旋转^[29]与本文设计的旋聚器之间的异同. 图 2(a)为无超材料的背景流场, 流线与等压线分别平行和垂直于$x$方向. 图 2(b)—图2(d)为不同功能下的速度场分布, 可以看出三者周围流域速度场均未受到干扰. 图 2(b)所示为传统流动旋转功能下的速度场分布, 在中心区域($0 < r \leqslant$$R_1$), 流体旋转了$\theta_0 = \pi/2$ ${\rm{rad}}$. 图 2(c)为传统流动聚集器的速度场分布图, 在环形区域($R_1 < r \leqslant R_3$), 流体逐渐朝着中心区域流动并最终汇聚于中心域. 本文设计的旋聚器的速度场分布如图 2(d)所示. 对比图 2(d)和图 2(c)可以发现传统流动聚集与流动旋转聚集的异同之处在于: 二者均可放大中心域($0 < r \leqslant R_1$)的流体流速, 说明二者均具备文丘里效应这一特性; 但是传统的流动聚集器只能将流体汇集于中心域, 无法改变此处流体的方向, 而文丘里效应旋聚器将流线和等压线汇集于中心区域的同时, 可改变流体在中心域的流向, 表明它对流体流向的控制程度更高. 为具体表现文丘里效应旋聚器在控制流体流向方面的灵活性, 下文将展示不同旋转预设角$\theta_0$对流体流向的影响.

图3所示为文丘里效应旋聚器在旋转角度$\theta_0$为$0 —$$\pi$ rad下的速度场分布. 在环形区域($R_1 <$$r \leqslant R_3$), 流体的流线与等压线呈逆时针旋转, 旋转角$\theta_0$越大, 流场的旋转程度越强. 可以看出, 不论环形区域内的流场如何旋转, 流体均会汇聚于中心区域, 说明文丘里效应旋聚器对流体施加的旋转和聚集效果相互兼容. 同时, 旋转角度$\theta_0$的变化不会影响文丘里效应旋聚器周围流域的速度场分布. 以上结果表明文丘里旋聚器在汇集流体的同时, 可灵活控制中心区域流体流向.

需要特别注意的是, 不论$\theta_0$为何值, 文丘里效应旋聚器中心域处流速方向均不能旋转至预定的角度, 出现旋转滞后现象. 如当$\theta_0 = \pi/2$ ${\rm{rad}}$ 时, 流体在旋聚器中的实际旋转角度约为$\pi/5$ ${\rm{rad}}$. 当$\theta_0 = \pi$ ${\rm{rad}}$ 时, 流体在旋聚器中的实际旋转角度约为$2\pi/5$ ${\rm{rad}}$. 经研究发现, 上述滞后现象为空间坐标变换顺序不同所导致. 当先对空间进行聚集再进行旋转时(即(7)式), 聚集空间与旋转空间会出现部分重叠区域, 致使旋转滞后现象的出现. 但当先对空间进行旋转再进行聚集时(见(12)式), 这种滞后现象便可被消除. 根据(12)式中新的空间映射关系以及(5)式和(6)式, 重新推导出对应的动力黏度张量表达式(13). 图4所示为坐标变换顺序为先旋转再聚集时的速度场分布. 从 图4可以看出, 旋聚器的旋转角度滞后现象被消除. 由此说明, 在使用坐标变换设计不同功能的超材料时, 需注意坐标变换的空间顺序关系.

$$\left\{\begin{aligned} &r^{\prime} = \frac{R_3 - R_1}{R_3 - R_2} r - \frac{R_2 - R_1}{R_3 - R_2} R_3, \\ &\theta^{\prime} = \theta - \theta_{0} \frac{R_3 - r}{R_3 - R_2} ,\\ &z^{\prime} = z, \end{aligned}\right.$$

(12)

$$\begin{array}{l} \begin{aligned} {\boldsymbol{\mu}}^{\prime} & = \left[\begin{array}{ccc} E & F & 0 \\ G & H & 0 \\ 0 & 0 & I\\ \end{array}\right] \mu \end{aligned} \end{array},$$

(13)

其中

$$\begin{split}&\\ E =\;& \frac{{\theta_0}^{2} r^{\prime} \left[(R_3 - R_2)r^{\prime} + (R_2 - R_1)R_3 \right]}{{(R_3 - R_1)}^{2} (R_3 - R_2)}\\ &+ \frac{r^{\prime}(R_3 - R_2)} {(R_3 - R_2)r^{\prime} + (R_2 - R_1)R_3},\\ G =\;& F = \frac{\theta_{0} \left[(R_3 - R_2)r^{\prime} + (R_2 - R_1)R_3\right]} {(R_3 - R_1)(R_3 - R_2)},\\ H = \;&1 + \frac{(R_2 - R_1)R_3} {(R_3 - R_2)r^{\prime}}, \\ I =\;& \frac{(R_3 - R_1)r^{\prime}} {(R_3 - R_2)r^{\prime} + (R_2 - R_1)R_3} \frac{R_3 - R_1} {R_3 - R_2}. \end{split}$$

为定量分析文丘里效应旋聚器对流速的放大效果, 同时考虑到$\theta_0$仅对流体流向产生影响, 所以我们取$\theta_0 = 0$ 时、不同$R_2/R_1$比值下的流场(图4(a))与背景流场(图 2(a))进行对比研究, 结果如图5所示. 可以看出, 流速在流域$|x/R_3 | \leqslant1$逐渐增加并在中心区域达到最大值, 该值为背景流域速度的$R_2/R_1$倍. 而流速在流域$|x/R_3 | > 1$均保持一致. 上述现象表明: 文丘里效应旋聚器具备放大流体流速的功能; 同时它的存在不会影响其周围流域的速度场分布, 可实现对附近流域的零干扰与自我隐身的效果.

由于在实际工作环境中, 来流往往会出现非均匀流动, 所以文丘里效应旋聚器应对非均匀流时的效果同样备受关注. 为此, 将图1中的进出口边界条件改成流体自左上入口流至右下出口, 其余条件保持不变, 此时的流场分布如图6所示. 观察图6可以发现, 文丘里效应旋聚器在来流为非均匀流时, 仍可稳定地发挥对于流体的旋转和放大流速的效果, 同时不影响其外流场的原有状态. 由此可知, 本文设计的文丘里效应旋聚器可以出色地应对不同状态的来流, 进一步提升了它的应用价值.

综上所述, 本文关于旋聚器的设计工作可以为变换光学^[1,2]、变换声学^[6,7]、变换热学^[12,13,20] 等变换理论所指导的旋转聚集器设计提供新的思路, 也可为不同功能超材料的设计奠定理论基础. 除此之外, 本文中旋转滞后现象与磁滞现象^[40,41]存在相似之处, 这可为磁滞现象的相关研究拓展至旋转滞后现象带来启发.文中所发现的空间坐标变换存在非互易性这一结论, 不仅适用于变换流体动力学, 同样也可拓展至不同领域下的变换理论.

4.   结　论

基于变换流体动力学, 通过耦合流动聚集和流动旋转功能, 设计出具备文丘里效应的流动旋转聚集器, 实现对任意方向上流体流速放大的功能. 通过推导文丘里效应旋聚器设计所需要的空间映射关系, 结合雅可比空间变换矩阵, 得出携带空间坐标变换信息的流体动力黏度张量, 实现了文丘里效应旋聚器的理论设计. 此外, 从数值模拟层面上验证并定性分析了文丘里效应旋聚器可在蠕动流状态下同时实现对流动流体的聚集和旋转功能, 通过改变预设旋转角$\theta_0$展现了文丘里效应旋聚器在调控流体流向上优越的灵活性; 同时, 我们发现并解释了流动旋转和聚集功能耦合时因作用空间相互干涉引起的旋转滞后现象的本质, 进一步得出空间坐标变换存在非互易性这一结论. 通过调整文丘里效应旋聚器的几何尺寸, 定量分析了文丘里效应旋聚器放大流速的功能, 以及它对背景流域零干扰和完美隐身的性能. 最后, 本文验证得出文丘里效应旋集器同样胜任于非均匀流场. 我们相信, 本文关于文丘里效应旋聚器的理论设计和数值模拟工作: 1)可为流动流体聚集和旋转功能在实验层面上的实现提供理论基础; 2)可为磁滞现象的相关研究拓展至旋转滞后现象带来启发; 3)可为超材料的非互易坐标变换设计提供新的思路.