不同金属导体间的温差可以通过热电效应将热量转换成热电电流，电流与磁场作用产生的洛伦兹力会驱动金属流体流动^[1]，这一流动被称为热电磁流体力学(thermoelectric magnetohydrodynamic，TEMHD)流动，从而将导电流固耦合系统中的温差转换成金属流体的动能，进而影响流场中的换热效率。液态金属作为高效的热输运介质，在受控热核聚变、电磁冶金以及高热流密度冷却器件等领域有重要应用价值。在碳中和背景下，基于当前我国多煤、缺油、少气的基础能源结构，发展清洁的受控热核聚变能源是十分必要的。基于此背景，本文研究热电效应在导电管道内产生的流动现象及这一流动现象对换热效率的影响。

托科马克装置是当前最具应用前景的磁约束受控热核聚变反应堆，包层和流动液态金属壁中的液态金属扮演着能量输运和氚增殖等重要角色。强磁场下产生的磁流体力学(magneto-hydrodynamics，MHD)效应会阻碍液态金属的流动，而热电电流与磁场作用则会产生驱动金属流体运动的洛伦兹力，在一定程度上抵消MHD效应对于流动的阻碍作用。在经典的MHD流动中，磁场强度是抑制流动的主要因素，磁场越强，MHD效应越明显，流动受到的阻力越大；而在TEMHD流动中，磁场强度既会增强MHD效应，也会强化热电效应产生的洛伦兹力。由欧姆定律以及洛伦兹力的表达式可知：磁阻尼效应产生的洛伦兹力与磁场强度的平方正相关，而热电效应产生的洛伦兹力与磁场强度正相关。因此，随磁场强度的变化液态金属会呈现不同的流动特征。

Xu等^[2]的研究表明不同磁场强度下热电效应会使金属流体呈现不同的流动结构，实验中流固界面上的温度梯度是由于瑞利伯纳德对流引起的：在中等磁场强度下，实验观察到了MP(magneto-precessional)流动模态，这一流动模态使瑞利伯纳德对流在垂直于磁场的平面内发生偏转；而在强磁场强度下MP流动模态转变为多涡胞MCMC(multi-cellular magneto-convection)流动模态，热电效应的影响消失。Zhang等^[3]实验研究了恒定温差下热电效应产生的TEMHD流动，观察到梯度磁场作用下热电效应产生的多涡旋环流结构，实验中测得的特征速度与温差成近似线性关系。Jaworski等^[4]在实验中观察到热电效应产生的旋转流动，且当磁场反向时，液态锂的流动也反向。这些研究表明热电效应产生的TEMHD流动可以呈现复杂的流动结构，如何利用这一效应强化金属流体在磁场作用下的换热效率具有重要的研究意义。管道是包层结构中液态金属流动的重要部件，液态金属在管道内的流动因此被广泛研究，Chen等^[5]模拟真实尺寸大哈特曼数和雷诺数下液态金属在包层内复杂管道中的流动情况，研究了壁面导电性对于MHD压降、流动形态的影响。Hussam和Sheard^[6]模拟磁场作用下管道内圆柱的尺寸、圆柱于壁面的间距以及雷诺数对于金属流体的流动和换热的影响，通过比较传热强化效果以及压降的大小确定了强化传热的最优圆柱位置。Modestov等^[7]通过实验和数值计算研究外加电流作用下矩形槽道内液态金属的流动和换热特性，结果表明：竖直方向上的洛伦兹力使金属流体中形成了涡结构，强化金属流体间的混合效率，使其温度分布更均匀，换热效果增强。Sahu等^[8]通过实验研究不锈钢薄壁管道中，磁场作用下铅锂流体的流动传热特性，管道底壁具有均匀的表面热流。根据浮力、电磁力和惯性力的相对强弱，得到3种流动状态。模拟结果表明：在浮力占主导地位的流态中，准二维湍流占主导地位，决定了整个传热机制；当电磁力占主导地位时，浮力被磁场抑制。在惯性力占主导地位的区域，电磁力和浮力不会决定液态金属的传热机制。Akhmedagaev等^[9]采用直接数值模拟和线性稳定性分析，研究底部均匀加热和外加横向水平磁场作用下水平管道内的混合对流，结合二维近似和全三维分析研究与核聚变反应堆的液态金属包层的典型工况下的流动传热特性，湍流状态下涡结构会使得流场中产生很大的温度波动；并且证明了在高哈特曼数下对湍流流动的准二维近似会产生较大的误差。Mistrangelo等^[10]讨论了与核聚变相关的磁流体动力学效应与其他物理现象，如氚输运、腐蚀和热传递耦合的研究现状。雷天扬等^[11]研究强磁场作用下微槽道内的液态金属流动和换热问题。热电效应作为一种驱动力，对于管道内金属流体的流动和换热特性的影响仍需进一步的研究。基于此，本文对这一问题进行建模，并通过数值计算研究其流动和换热规律。

1 物理模型与数值方法 1.1 计算物理模型

本文涉及的数学符号由表 1进行了汇总。

液态金属锂和不锈钢的物性参数由表 2给出。

计算物理模型以及几何尺寸如图 1所示，方形管道内一部分壁面由导电不锈钢组成(图 1阴影部分)，其余壁面由绝缘材料组成。管道上表面的温度为T_h，下表面温度为T_c，且T_h>T_c，在聚变堆中，偏滤器的上表面承受来自等离子体的高温热流，下表面一般会用冷却介质进行冷却，以保证液态金属不会发生蒸发现象。由于绝缘材料不会产生热电效应，本文的计算中忽略了绝缘固体材料对于温度场的影响。计算模型的初始条件以及边界条件将在下文给出，方形管道的几何尺寸为[-a, a]×[0, 16a]×[-a/2, a/2]，不锈钢导电壁面的厚度为d，且d=a/10。磁场方向沿x方向，流动方向沿y方向，在真实的聚变堆中，磁场方向也是沿着第一壁的展向方向。

1.2 控制方程

保证计算过程中电荷的守恒性是磁流体动力学模拟的一个关键问题，本文采用计算流体力学开源软件OpenFOAM^[12]中的PISO^[13]算法求解速度压力耦合关系，采用Ni等^[14]发展的有限体积框架下任意网格的相容守恒格式计算网格面心处的守恒型电流密度通量，并由这一电流密度通过守恒型插值得到网格体心处的洛伦兹力。该洛伦兹力作为显式的源项进入离散代数方程组的右端项，流固界面处电流密度和热通量的连续性条件采用分区迭代算法来处理，详细的求解过程已在Chen等^[15]的论文中介绍，此处不再赘述。

1.2.1 流体域的控制方程

对于大部分工程应用中的液态金属流动过程，其磁雷诺数很小$\left(R_{m}=\mu_{m} \sigma U L \ll 1\right)$，诱导磁场相对于外加磁场可以忽略不计。在Quasi-static Approximation^[16]这一假设下，得到磁流体的运动和传热方程如下

$\begin{gathered} \nabla \cdot \boldsymbol{u}=0 . \end{gathered}$

(1)

$\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \cdot \nabla \boldsymbol{u}=-\frac{1}{\rho_{\mathrm{f}}} \nabla p+\nabla \cdot\left(\nu_{\mathrm{f}} \nabla \boldsymbol{u}\right)+\frac{1}{\rho_{\mathrm{f}}} \boldsymbol{J} \times \boldsymbol{B} .$

(2)

$\rho_{\mathrm{f}} C_{p \mathrm{f}}\left(\frac{\partial T}{\partial t}+\boldsymbol{u} \cdot \nabla T\right)=\nabla\left(\lambda_{\mathrm{f}} \nabla T\right) .$

(3)

考虑到SEEBECK效应产生的热电电流，将欧姆定律和电荷守恒定律分别由以下2个公式表示

$\boldsymbol{J}=\sigma_{\mathrm{f}}\left(-\nabla \varphi+\boldsymbol{u} \times \boldsymbol{B}-S_{\mathrm{f}} \nabla T\right) .$

(4)

$\nabla \cdot \boldsymbol{J}=0 .$

(5)

由欧姆定律式(4)和电荷守恒定律式(5)得到的电势泊松方程来计算电势。对于纳维斯托克斯方程组中的相关物理量，采用a, υ_f/a, a²/υ_f, ρ_fυ_f²/a²对长度、速度、时间和压力进行无量纲化，采用B, υ_fB, σ_fυ_fB/a, ΔT对磁场、电势、电流密度和温度进行无量纲化，相关无量纲参数主要有：雷诺数(Re)，哈特曼数(Ha)以及热电数(Te)。

1.2.2 固体域的控制方程

固体域的热传导、欧姆定律和电荷守恒定律可写为

$\rho_{\mathrm{s}} C_{p \mathrm{~s}} \frac{\partial T}{\partial t}=\nabla \cdot\left(\lambda_{\mathrm{s}} \nabla T\right) .$

(9)

$\boldsymbol{J}=\sigma_{\mathrm{s}}\left(-\nabla \phi-S_{\mathrm{s}} \nabla T\right) .$

(10)

$\nabla \cdot \boldsymbol{J}=0 .$

(11)

固体区域的电势由欧姆定律式(10)和电荷守恒定律式(11)得到的电势泊松方程来计算。

1.2.3 物理边界条件

管道入口处的速度为均匀来流，大小为U_in；出口处速度为充分发展，壁面处的速度满足无滑移边界条件。

$\left\{\begin{array}{l} \varphi_{\mathrm{w}}=\varphi_{\mathrm{f}}=\varphi_{\mathrm{s}}, J_{\mathrm{s} n}=J_{\mathrm{f} n} \\ T_{\mathrm{w}}=T_{\mathrm{f}}=T_{\mathrm{s}}, q_{\mathrm{s} n}=q_{\mathrm{f} n} \end{array} .\right.$

(12)

$\frac{\partial \varphi}{\partial n}=(\boldsymbol{u} \times \boldsymbol{B}-S \nabla T) .$

(13)

在流固耦合系统中，耦合界面处的物理量要满足式(12)中描述的热和电的连续性和通量守恒关系。在流固界面处，两侧的法向电流密度和热流密度相等，由此可得到流固耦合壁面处的电势和温度表达式，推导过程已在Chen等^[15]介绍过，此处不再重复列出。为满足电荷守恒关系，管道入口和出口的法向电流通量为零(J_n=0)，由式(4)、式(5)可得该边界条件的表达式为式(13)；其余绝缘管壁也适用这一边界条件。管道上、下表面的温度分别为T_h, T_c，管道入口处温度为T_c，其余壁面为绝热边界条件。

2 模型验证

Hunt^[17]给出磁场作用下金属流体在壁面导电的方管中流动的解析解，通常用这一解析解来验证算法。对于Ha=1 000的Hunt-Ⅱ流动，采用Ni等^[14]发展的相容守恒格式以及分区迭代算法计算得到的数值解与Hunt解析解^[17]的对比，在划分网格时，哈特曼层内至少分布6个均匀网格，侧层内至少分布8个均匀网格，对比结果可参考Chen等^[15]的文章，此处不再赘述。

3 结果与讨论 3.1 对Hartmann流动的分析

Hartmann流动是两无限大平板间液态金属的理想流动，其示意图如图 2(a)所示，平板的厚度为$d$，无量纲壁面导电比为$C_{\mathrm{w}}=$ $\left(d \sigma_{\mathrm{w}}\right) /\left(a \sigma_{\mathrm{f}}\right)$，磁场方向沿$y$方向，在充分发展的假设下，流速仅在板间($y$方向)发生变化。本节对文献中给出的不同壁面导电比，由压力梯度驱动的Hartmann流动以及由热电效应驱动的Hartmann流动解析解进行分析。

Muller等^[18]给出了由恒定压力梯度驱动的Hartmann流动解析解, Shercliff^[1]基于薄壁假设以及速度温度单向耦合假设，得到板间平均流速表达式。无量纲流速与磁场及壁面导电性的关系如图 2(b)所示，图中虚线表示的是由压力梯度驱动的Hartmann流动，实线表示的是由热电效应驱动的Hartmann流动。在相同条件下，壁面导电性越强，流动强度越小。当压力梯度保持不变时，磁场越大，阻碍流动的洛仑兹力越强，因此流动速度越小；而对于热电效应驱动的流动，当磁场强度小于临界磁场时，速度随磁场增强而增大，当磁场强度大于临界磁场时，流动速度随磁场增大而减小。基于速度温度单向耦合假设的Hartmann解析解反映了热电效应的作用规律，实际上，速度和温度通过热电效应耦合在一起，是双向耦合的，而在推导解析解时，将这一关系简化为单向耦合关系，即速度不会影响温度分布。在某些参数范围内，这一假设可能是合理的，但是其适用范围无法通过理论给出，只能通过实验或数值计算得到的数据来确定。

3.2 不考虑热电效应的方管流动

本文建模时忽略了浮力效应的作用，因此在不考虑热电效应和浮力效应的管道流动中，速度与温度是解耦的，本节计算1.1节中的物理模型在Re=764.5, Ha=25、50、100、200等工况下的流动，研究速度分布、压降等规律，该雷诺数下，管道内为层流。对于该物理模型，在沿管道流向的非导电段，金属流体在磁场作用下会表现出Shercliff流动的主要特征，而在导电段，金属流体在磁场作用下会表现出Hunt流动的特征。

Re=764.5, Ha=50、100的速度分布云图、电流线和洛伦兹力云图如图 3所示，随着磁场增大，管道内的速度截面在侧壁出现了Hunt流动中典型的M型分布，这是由于侧层附近的电流线与磁场方向的夹角逐渐减小，洛伦兹力逐渐减小。

侧壁中心轴线处的速度分布如图 4(a)所示，图中的位置坐标和速度分量分别采用长度a和入口速度U_in进行无量纲化，在侧壁附近洛伦兹力很小，因此流动速度会增大，出现M型速度分布，在中心流动区域，电流线与磁场几乎是垂直的，因此会产生阻碍流动的洛伦兹力，且磁场越强，阻碍流动的洛伦兹力越强；哈特曼层附近的电流线方向与中心流动区域的电流线方向相反，使得粘性边界层内的洛伦兹力的方向与中心主流区域的洛伦兹力方向相反，如图 3所示，在黏性边界层内的洛伦兹力促进了液态金属的流动，所以哈特曼层内的速度梯度非常大，且随磁场增大，哈特曼层内的速度变化越剧烈。

管道中心轴线处的沿程压降由以下公式计算得到

$\Delta p=\left(p-p_{\text {out }}\right) /\left(\frac{1}{2} \rho_{\mathrm{f}} U_{\text {in }}^{2}\right) .$

(14)

不同磁场强度下，管道中心轴线处的沿程压降(Δp)如图 4(b)所示。当管道入口流速不变时，磁场越大沿程压降越大，说明管道中金属流体所受的洛伦兹力变大，且压力梯度与洛伦兹力相平衡。流经导电管道时，金属流体受到的洛伦兹力比非导电段更强，所以压降更大。

本节首先对不考虑Seebeck效应时，部分导电管道内金属流体的流动特征进行了分析，以下将对Seebeck效应作用下部分导电管道内金属流体的流动和换热规律进行讨论。

3.3 热电效应作用下的方管流动

对于这一流动问题，相关的无量纲控制参数为雷诺数、哈特曼数以及热电数, 对这一问题进行参数化研究才可以得到不同无量纲控制参数下的主要流动特征，本文作为初步的工作，计算了部分导电管道内，无量纲参数为Re=764.5, Ha=25、50、100、200等工况下的热电效应作用下的流动。图 5给出了不同磁场强度下y=8a处的温度分布云图、速度分布云图、电流线和洛伦兹力分布云图。

对于小雷诺数的情况，可以认为流动是由热电效应驱动的，热电效应产生的洛伦兹力所做的功转化成了流体的动能，流场中的主导电流是热电效应产生的热电电流；而对于大雷诺数的情况，流动是由惯性力所驱动的，热电效应产生的洛伦兹力不是主导作用，流场中的主导电流是流体运动产生的感应电流，而不是热电电流。热电效应产生的热电电流在磁场作用下将温度场与速度场耦合在一起，在本节的计算工况中，Re=764.5，当磁场强度较小时(Ha=25)，流动是非稳态的，温度、速度以及洛伦兹力等物理量的分布是不对称的，且随着时间发生变化；当磁场较大时(Ha≥50)，流动是稳态的，流动结构具有较好的对称性，这一点可以从图 5中的云图和电流线观察到。这说明在小磁场下，热电效应产生了非稳态的流动，而随着磁场增强，非稳态的流动逐渐转变为稳态流动，这也说明在小磁场时，金属流体的速度和温度是双向耦合的，不能简化成单向耦合关系，而在大磁场情况下，这一单向耦合假设才可能是正确的。非稳态流动到稳态流动的转变过程与雷诺数、哈特曼数以及热电数等无量纲参数相关，后续将会进行研究。

对图 5进行进一步分析可知，流固界面上的热电效应在方管的4个角点附近非常显著，且热电电流与磁场作用产生了与流动方向相反的洛伦兹力，哈特曼层附近的洛伦兹力分布与不考虑热电效应时的洛伦兹力分布明显不一样。在洛仑兹力的作用下，4个角点处形成与流动方向相反的回流涡结构，热电效应产生的回流涡形态和流线如图 6(a)和6(b)所示，这一回流涡结构使得管道截面处的温度呈现哑铃型分布。

对比图 4(a)和图 6(c)中侧壁中心线处的速度曲线可知：对于较小的磁场，热电效应使流场中的速度分布发生了明显变化。角点处流向回流涡结构的形成使中心区域的流速增加，进而促进了中心区域的对流换热；而随着磁场增强，流动由非稳态变成稳态，中心区域的流动被磁阻尼效应主导，热电效应的作用尺度逐渐由整个流动平面收缩至边界层附近，出现非对称的M型速度分布。

4 结论

本文以液态金属锂和不锈钢为研究对象，对不考虑热电效应时的流动现象以及热电效应作用下的流动现象进行模拟，磁场方向沿展向。

1) 不考虑热电效应时，研究了较大雷诺数(Re=745.6)下，不同磁场强度作用(Ha=25、50、100、200)下的管道内液态金属流动现象，研究表明流场中惯性力与洛伦兹力相平衡，导电段的侧壁会产生对称的M型速度分布；随磁场强度增大，管道沿程压降也增加，且流经导电段时，压降损失更大，流动现象以及相关的变化规律与前人的研究工作符合。

2) 热电效应作用下，发现热电效应与磁场作用产生的洛伦兹力在方管的4个角点处形成与流动方向相反的回流涡，且回流涡主要分布在角点的壁面附近，回流涡结构增大了中心区域的流速，促进了中心区域金属流体的对流换热。随着磁场强度增大，流动由非稳态流动变为稳态流动，速度与温度从双向耦合转变为单向耦合，中心区域的流动被磁阻尼效应主导，热电效应的作用尺度逐渐收缩至边界处附近，侧壁会出现非对称的M型速度分布。

3) 热电效应在管道截面上产生的不均匀洛伦兹力分布使金属流体产生了局部加速现象，生成涡结构来强化换热。这一过程可以充分利用聚变堆中等离子体反应生成的热量。相比于在管道中放置绕流物产生涡结构，这种方式不会产生额外的流动阻力，而且温度梯度越大的区域热电效应越明显，这意味着高温区域的热量可以及时输运至低温区域。