数值研究及解析建模：变磁场作用下矩形管道内铅-锂流动的MHD压降（用于核聚变反应堆） Silvia Iannoni Gianluca Camera Marcello Iasiello Nicola Bianco Giuseppe Di Gironimo

《Journal of Nuclear Engineering》：Numerical Investigation and Analytical Modeling of MHD Pressure Drop in Lead–Lithium Flows Within Rectangular Ducts Under Variable Magnetic Field for Nuclear Fusion Reactors Silvia Iannoni, Gianluca Camera, Marcello Iasiello, Nicola Bianco and Giuseppe Di Gironimo

　　摘要
繁殖毯是托卡马克的关键组成部分，主要负责从聚变反应中提取热量以及进行氚的繁殖，这是确保聚变反应堆燃料自给自足的关键。最近的技术进步导致了双冷却铅-锂（DCLL）繁殖毯的发展，这种毯子使用液态金属（具体来说是铅-锂共晶合金）作为传热介质和氚繁殖器，而氦气则用于冷却反应堆的结构部件。在托卡马克环境中，流动的电导流体与强磁场之间的相互作用会产生磁流体动力学（MHD）效应。这些效应表现为流体内部涡流的产生以及与磁场相互作用产生的洛伦兹力，这会导致额外的压力损失并降低传热效率。本研究探讨了在不同强度的外部均匀磁场作用下，铅-锂流体在矩形截面管道中经历的压力降。开发了一个分析模型来估算不同外部磁场强度下的总MHD诱导压力损失，并将其与使用COMSOL Multiphysics进行的相对计算流体动力学（CFD）模拟进行了对比。这种比较有助于验证分析预测的准确性，并更好地理解施加的磁场强度对整体压力降的影响。因此，该分析模型的目的是为未来初步设计提供合理准确的MHD压力损失估算工具。

1. 引言
追求可控核聚变作为一种可行的能源激发了大量研究，旨在开发能够在确保长期运行可靠性的同时维持高功率密度的技术。在磁约束聚变装置（如托卡马克）中，聚变反应释放出大量的能量，这些能量必须被有效提取，同时通过氚的繁殖来维持燃料的自给自足。这些要求对暴露在极端热、中子和电磁环境中的反应堆部件提出了严格的要求。在这个框架下，繁殖毯起着核心作用，因为它既负责热量的提取也负责氚的生产，直接影响未来聚变电厂的整体性能和可行性。当前的聚变路线图将繁殖毯概念的验证确定为实现示范性（DEMO）级反应堆及更高级别的关键步骤[1]。在各种正在研究的毯子设计中，双冷却铅-锂（DCLL）概念作为一种有前景的解决方案脱颖而出，它结合了液态金属繁殖器和冷却剂，并为结构部件提供了独立的氦冷却回路[2]。使用铅-锂共晶合金使得在高温下运行成为可能，提高了热效率，同时提供了有效的氚繁殖能力。然而，在强磁场中采用电导液态金属引入了额外的物理复杂性。在聚变反应堆环境中，施加的磁场与液态金属的运动之间的相互作用会产生磁流体动力学（MHD）现象，其特点是电流的诱导和由此产生的洛伦兹力对抗流动。这些效应显著改变了速度分布，抑制了湍流，从工程角度来看，导致压力损失大幅增加[3,4,5,6]。由此产生的MHD诱导压力降直接影响泵送功率需求，是液态金属基毯子概念的主要设计约束[7,8]。
矩形导管在这一背景下尤为重要，因为它们代表了繁殖毯内的内部冷却通道和流动路径。因此，导管几何形状、壁电属性和磁场强度对MHD流动行为的影响已经得到了广泛研究[9,10]。虽然高保真数值模拟可以提供三维电流路径和流动结构的详细见解，但它们通常计算要求过高，不适合早期设计活动。因此，仍然迫切需要紧凑的分析模型，能够在与聚变相关的条件下提供可靠的MHD压力损失估算，以支持初步设计和优化研究。在这项工作中，研究了在均匀外部磁场作用下，通过矩形截面导管的铅-锂流动中的MHD诱导压力降。研究重点量化了压力损失对磁场强度的依赖性，这些参数范围代表了DCLL繁殖毯的应用。开发了一个分析模型来估算沿导管的总MHD压力降，明确考虑了流体核心、近壁层和导电壁之间的电流重新分布。然后使用COMSOL Multiphysics进行的CFD模拟对分析预测进行了验证，以便评估模型的准确性及潜在的物理假设。
本研究建立在作者之前关于圆形导管中铅-锂流动的MHD压力降的数值和分析研究基础上[11]，将提出的建模方法扩展到矩形几何形状，并重点分析磁场强度的影响。所提出的方法旨在为估算矩形导管中的MHD压力降提供实用且计算效率高的工具，适用于聚变反应堆毯中液态金属部件的初步设计评估。

2. 文献综述
关于电导液体中MHD流动的研究可以追溯到Hartmann的开创性工作，他首次分析了均匀磁场存在下的层流，并强调了流动阻力的显著变化[3]。紧接着，Shercliff将理论描述扩展到横向磁场下的导管和管道配置，为壁限制MHD流动提供了关键参考[4]。Moreau后来建立了成熟的理论框架[5]，Davidson进一步通过使用无量纲参数（如Hartmann数和交互参数）对MHD流动的物理解释和缩放进行了阐述[6]。从更应用的角度来看，Bühler对通道和容器中的MHD流动进行了全面研究，特别关注了与导管流动相关的Hartmann数和边界层[12]。对中等Hartmann数和交互参数下的圆形导管中液态金属MHD流动的进一步分析有助于明确实际感兴趣的流动状态[13]。Smolentsev及其同事提供了关于聚变相关液态金属毯中MHD热流体动力学问题的综述，强调了电流闭合路径、壁电属性和三维效应对压力降的影响[7]。在聚变路线图和相关测试需求中，多次强调了繁殖毯对DEMO开发和验证的重要性[1]。在这个背景下，DCLL概念因其高温运行和氚繁殖的潜力而受到了广泛关注，同时也带来了显著的MHD驱动的液压难题[2]。
除了总体毯子考虑之外，还有几项数值研究针对代表毯子通道的导管配置进行了研究。Sterl首次模拟了矩形导管中的液态金属MHD流动[14]，为横向磁场对非圆形截面中的速度分布和压力降的影响提供了参考。最近，Smolentsev及其同事的三维模拟研究了具有导电壁的矩形导管中的MHD流动，证实了三维电流路径和壁效应的重要性[15]。最近的研究还解决了毯子接入导管和复杂通道的问题，在直接设计相关的几何形状中提供了压力降特性[16]。其他研究关注了毯子歧管和系统级设计约束，为液态金属毯组件提供了指导方针和数值分析[9,17]。还在非均匀磁场下研究了弯曲管道配置，以代表现实托卡马克布局的典型几何特征[10]。尽管文献提供了大量关于毯子相关通道中MHD诱导压力降的数值证据，但适用于初步设计的紧凑分析工具仍然相对有限，特别是在涉及几何复杂性和电磁电流重新分布的情况下[8]。最近的研究提出了基于相关性的方法和三维MHD压力损失的优化策略[18]，进一步证实了对适用预测工具的需求。在最近的一项数值和分析研究中，作者研究了在均匀且单位磁场下，通过圆形导管的铅-锂流动中的MHD压力降，开发了简化的分析相关性，并将其与COMSOL模拟进行了对比[11]。本研究基于此工作，采用了相同的方法，将分析模型扩展到不同的几何配置，并进行了一项参数研究，改变的是磁场强度而不是导管倾斜角度。

3. 问题实施
3.1. 几何形状和材料
所研究的导管具有外部矩形横截面，尺寸为???? ×????（???? =65 mm 和 ???? =60 mm），壁厚为???? =3 mm。它由沿x方向排列的水平部分（长度为6 ????）和沿y方向排列的垂直部分（长度为3.5 ????）组成。这两部分通过内外侧的弯头连接，弯头的半径为???? =????（图1）。图1. 导管几何形状。分析中考虑使用铅-锂作为流体，导管壁的材料为EUROFER-97 [19]。相应的物理属性在表1中报告。在整个研究中，几何配置保持不变。
表1. 铅-锂和EUROFER-97的物理属性 [10,19]。

3.2. 物理和边界条件
PbLi流动在无感应近似[20]下进行分析，该近似忽略了由感应电流产生的磁场对施加磁场的影响。因此，控制方程是通过将Navier-Stokes方程与Maxwell方程耦合得到的，形成了方程（1）-（6）中报告的系统。
动量：???(??·?)???=????+???2??+?? （1）
质量连续性：?·??=0 （2）
洛伦兹力：??=??×?? （3）
广义欧姆定律：??=???(??+??×??) （4）
静电场：??=???? （5）
电荷守恒：?·??=0 （6）

边界条件包括规定的平均入口速度（??in =0.02 ?m/s）、出口压力（??out =0，作为参考相对压力）以及导管壁的无滑移条件（??wall =0）。假设壁的外表面是电绝缘的（?????? ·?? =0）[11]。
流体域受到沿y方向（???? =???? =0）定向的均匀环形磁场的作用。磁场强度在所有模拟中从???? =0.5 T变化到???? =2.5 T，增量为0.1 T。

相应的无量纲参数是雷诺数、Hartmann数和交互参数，分别为：
Hartmann数：Ha=??????√???? （7）
雷诺数：Re=????in?????? （8）
交互参数：N=????????2????in （9）
这里，??表示PbLi的电导率，??是其动态粘度，??是其密度，??in是流体的入口速度，????是特征长度。后者定义为 ????=???2，其中液压直径???由 ???=23????3????????+???? 给出，其中 ???? =???? ?2???? =59 mm 和 ???? =???? ?2???? =54 mm。所有模拟中的雷诺数保持恒定，等于Re =5237，与磁场强度无关。相反，Hartmann数Ha和交互参数N明确依赖于磁场强度。当???? =0.5 T时，Ha和N的最小值分别为379和31.6，而当?? =2.5 T时，最大值分别为Ha =1890和N =790。
尽管雷诺数值在形式上将流动置于过渡区域（2300
3.3. 计算网格和数值收敛性
计算域使用comsol生成的三维非结构化四面体网格进行网格划分。 <104）[21]，但相对较高的n值表明洛伦兹力显著主导了惯性效应。实际上，hartmann数的平方代表了洛伦兹力与粘性力的比率，而交互参数量化了电磁力与惯性力的比率[5]。在这些强磁场下的液态金属流动条件下，随着交互参数的增加，洛伦兹力倾向于抑制速度波动并稳定流动[22,23,24]。特别是，控制mhd导管流动转变的相关参数已被证明是比率 ??=Re/Ha，该比率基于Hartmann层厚度，其临界值约为?? ≈380 [24]。在本研究中，所研究的整个磁场范围内，re ha的相应值远低于这一阈值，使得湍流效应可以忽略不计，因此流动基本上是层流的。 3.3. 计算网格和数值收敛性>
3.3. 计算网格和数值收敛性
计算域使用comsol生成的三维非结构化四面体网格进行网格划分。>特别关注了具有强速度梯度和强烈电磁效应区域的求解，即靠近壁面的边界层（Hartmann层和Shercliff层[3,4]）。采用了非结构化的三维混合网格，主要由四面体元素和棱柱体元素组成，特别是在靠近壁面的区域。最终离散化包含了大约5·10^6个元素和10^6个节点。随着施加的磁场强度的增加，预期MHD效应会逐渐增强，导致Hartmann层和Shercliff层变得更薄，速度和电流密度的梯度也变得更陡峭。为了确保在这些日益严格的条件下数值稳定性和收敛性，随着磁场强度的增加，全局网格密度也随之系统地增加。采用这种细化策略是为了充分解析与更高Hartmann数相关的逐渐变薄的边界层，并防止在相互作用参数较大时收敛性下降。为了进一步提高需要的空间分辨率，在管道的入口截面处进行了局部网格细化。其余体积使用非结构化的四面体公式进行离散化，确保整个域内的网格一致性以及精细区域和整体区域之间的平滑过渡。所得到的网格在保持整个域内数值稳定的同时，提供了对Hartmann层和Shercliff层的足够分辨率。通过在对具有代表性的磁场强度的逐步细化网格上进行模拟并监测沿管道的压力降来评估网格的独立性。最终网格根据磁场强度的不同，包含了大约??2(10^6)个元素，作为数值精度和计算效率之间的折中方案。

4. 分析模型
开发了一个分析模型来估计导电PbLi流与外部施加的磁场相互作用所产生的MHD压力损失。该模型的构建遵循了[11]中采用的概念框架：它基于控制管道横截面内和穿过壁面的电流产生、重新分布和闭合的物理机制。根据无感应假设[11]对动量方程进行无量纲化处理，得到1??(???·?)???=?????+1Ha2?2???+???×??? (10)，不同力贡献的相对重要性可以通过相互作用参数N和Hartmann数Ha来评估。对于本研究相关的参数范围，其特点是N的值较大，顺流方向的压力梯度主要由洛伦兹力项平衡，而惯性贡献则相对次要。因此，MHD引起的压力梯度可以直接与电磁力密度?? ×??相关联。为了简化问题，我们设定???? =??和??in =??。忽略电势的轴向梯度，假设诱导电流位于管道横截面上，其闭合路径如图2所示。由于?? ×??项产生的电压，电流在流体中被诱导，并倾向于沿着由截面几何形状和流体及壁材的电学性质决定的优先路径闭合。特别是，我们假设一部分电流流经截面的核心区域，产生制动洛伦兹力，而另一部分电流则通过靠近壁面的层和管道壁闭合，产生局部加速的力效应。

5. 数值结果
5.1. 数值压力降和速度场评估
首先在COMSOL上进行了没有磁场效应的CFD模拟，以测试代码并为后续分析建立一个可靠的基准案例。图4显示了在没有磁场效应的情况下，对于?? =0的层流情况，????平面上的速度大小分布。流动模式与经典流体动力学行为[25]一致，表现出核心区域的加速和由于粘性效应导致的壁面附近速度的降低。在弯道内壁下游发生流动分离，而在垂直截面的外壁上流动被压缩并因此加速。对于层流情况计算得到的总压力降为Δ????????=5.9·10^-1Pa (38)。
随后，如前所述，进行了多次CFD模拟：在保持几何形状不变的情况下，磁场强度从???? =0.5 T变化到???? =2.5 T，每次增量为0.1 T。图5显示了作为施加磁场强度函数的CFD模拟得到的总MHD压力降。观察到Δ??MHD随B的增加而单调增加，相应的压力降比层流情况下的Δ??值高两个到四个数量级。这一趋势是非线性的，符合经典MHD缩放定律，即电磁效应与磁场强度的平方和相互作用参数成正比[5,6]。随着磁场的增加，惯性效应逐渐减弱，流动进入一个受MHD强烈支配的阶段，导致流动阻力显著增加[12]。
图6显示了?? =0.5 T时的速度大小分布。与纯流体动力学层流情况相比，由于洛伦兹力的作用，速度场已经出现了明显的改变。核心区域变得部分扁平化，动量在弯道处重新分布。特别是，靠近内壁的分离区域发生了位移，垂直截面开始出现由于磁场引起的分层现象。尽管电磁效应已经可见，惯性在流动动力学中仍然起着重要作用。
图6. 在?? =0、?? =0.5 T时，????平面上的速度大小分布。在中等磁场强度?? =1.6 T时（见图7），MHD效应变得占主导地位，强烈改变了速度场。整体流动受到显著制动，导致核心速度降低，水平截面上的速度分布更加平坦，同时在靠近壁面的区域形成了局部加速区。值得注意的是，在有限元（FEM）模拟中，所有固体壁面都施加了经典的无滑移边界条件，因此在边界处的速度严格为零。观察到的靠近壁面的明显加速现象对应于在横向磁场作用下的MHD管道流动中典型的M形速度分布，其中洛伦兹力抑制了核心流动，而在靠近壁面的薄边界层内形成了局部速度最大值，随后在壁面处速度降为零（见图8）。图7显示了在?? =0时，?? =1.6 T的?????平面上的速度大小分布。图8显示了yz平面上的轴向速度。在弯曲处，流动分离现象被大大抑制，速度场通过惯性力、粘性力和电磁力之间的平衡而变得更加有序。对于考虑的最高磁场强度?? =2.5 T（见图9），流动受到电磁效应的强烈控制。速度场在弯曲的内壁附近显示出明显的加速现象，并且存在一个速度减小的受限区域。虽然水平部分几乎保持均匀，但垂直部分形成了一个中心速度减快的狭窄区域，且靠近壁面的加速作用更强烈。图9显示了在?? =0时，?? =2.5 T的?????平面上的速度大小分布。总体而言，这三个速度图在所研究的磁场范围内显示出相似的流动特征。特别是，随着B的增加，弯曲内壁附近的加速现象变得更加明显，而垂直部分的速度梯度逐渐增强。在图10中，平均速度分量????、????和????作为垂直坐标y的函数在弯曲下游的直立垂直部分进行了绘制；为了获得????、????和????的值，将垂直部分划分为30个?????平面，并在每个平面上提取了速度分量u、v和w的空间平均值。值得注意的是，平均而言，垂直方向上的速度分量变为负值，这表明由于洛伦兹力的作用形成了再循环区域，导致流动方向发生了反转。这一现象也可以在图6、图7、图8和图9中观察到：在所有情况下，都观察到了弯曲下游垂直部分外壁附近的局部加速现象。这种现象可以这样解释：由于弯曲的作用产生了横向速度分量，这些横向速度分量所关联的洛伦兹力主要沿+x方向作用，位于?????平面内，因此与垂直部分的主要流动方向不对齐。因此，这种力不会直接导致沿流动方向的压降，而是促进了动量的侧向重新分布，使流体向外壁方向移动，从而导致局部流动加速和核心减速。实际上，图10清楚地显示了顺流方向的分量????相对于横向分量仍然是主导的。图10显示了通过连续截面平均速度场得到的垂直部分中沿垂直坐标y的平均速度分量????、????和????。图11和图12分别展示了在?? =0.5 T的情况下，水平和垂直部分上?????平面和?????平面上的速度大小分布。在两个平面上，管道核心区域都观察到了制动效应，而在靠近壁面的层内形成了局部加速。这种流动组织与预期的MHD行为[5,6,12]定性一致，其中洛伦兹力导致动量重新分布，表现为核心速度降低和壁面速度增加，这与分析电流分布模型背后的物理假设相符。图11显示了?? =0.5 T时xz平面上的速度大小分布。图12显示了?? =0.5 T时yz平面上的速度大小分布。

5.2 分析与数值压降的比较
将MHD引起的压降的分析预测与整个磁场强度范围内通过计算机流体力学（CFD）模拟得到的结果进行了比较。图13报告了从分析模型和数值模拟中获得的总MHD压降Δ???MHD作为磁场强度B的函数。总体而言，分析模型与计算预测的MHD引起的压降吻合良好，不仅在正确趋势上一致，而且在所研究的磁场范围内也保持了正确的数量级。在较高磁场强度B下观察到轻微的高估，因为更强的电磁效应促进了更为显著的三维流动结构，而这些结构仅部分被分析模型的简化假设所描述。在更高的磁场强度下，特别是在管道的垂直部分，观察到了明显的速度梯度和壁面局部加速效应。在这个区域，平均流动方向大致与施加的磁场对齐，同时在截面内局部产生了不可忽视的垂直于磁场方向的速度分量。因此，相关的洛伦兹力主要沿着横向方向作用，这与经典MHD理论[5,12]一致，并不会直接导致顺流方向的压降。这也体现在分析模型中，该模型仅考虑了事先已知的主要流动方向。根据分析公式，MHD引起的压降损失来源于管道的水平和弯曲部分，在这些部分垂直于磁场方向的速度分量显著，并产生了顺流方向的洛伦兹力贡献。因此，垂直部分没有明确包括在MHD引起的压降评估中。为了评估这部分的可能贡献，通过比较弯曲直下方区域的面积平均压力与垂直管道出口处的压力，从CFD模拟中评估了垂直部分沿线的压力变化。在所研究的磁场范围内，所得到的压力变化大约为10?1 Pa，这比管道全程的总压降小几个数量级。这些结果证实，忽略垂直部分的压降不会影响整体结果。为了进一步测试分析模型的鲁棒性，还对不同的壁电导率（即不同的????,??值）进行了敏感性分析。将管道材料改为CuCrZr（?? =4.64 ×105 S m?1）或F82H（?? =1.37 ×106 S m?1）时，分析预测仍然随施加的磁场变化而变化，见图14。图14显示了不同管道材料的计算与分析MHD压降比较。与EUROFER材料相比，观察到了更大的差异，尽管对于分析估计来说仍然是可以接受的。实际上，EUROFER、CuCrZr和F82H的最大百分比误差分别为11.3%、28.4%和13.7%，其中CuCrZr的最大百分比误差仅为42.8 Pa。这些结果证实，即使在参数（如壁面导电率）发生变化时，该模型也适用于初步设计阶段。

6. 结论
在这项工作中，提出了一种对在均匀磁场作用下滑动的铅-锂（Lead–Lithium）流通过矩形管道时MHD引起的压降进行的数值和分析研究。通过改变固定范围内的磁场强度进行了参数分析，从而表征了所得到的压降和流动结构。CFD模拟突显了压降对磁场强度的强烈依赖性，随着相互作用参数的增加，洛伦兹力逐渐主导了流动动力学。数值结果还揭示了在MHD条件下与弯曲几何形状相关的特征流动现象，包括管道壁附近的局部加速以及核心区域和靠近壁面区域之间的动量重新分布。为了估计MHD引起的压降，开发了一个基于管道横截面电流重新分布等效电路表示的分析模型。分析预测与计算结果吻合良好，再现了整个磁场强度范围内压降的整体趋势和正确数量级。在较高磁场强度B下观察到略微的高估，这可能与分析公式中采用的简化假设有关。总体而言，所提出的分析公式为估计矩形管道中的MHD引起的压降提供了一个实用且紧凑的工具，这些管道代表了DCLL增殖毯层通道。尽管对电流路径和流动组织进行了简化假设，该模型仍能在需要快速评估MHD效应的初步设计和参数研究中表现良好。未来的发展可能会着眼于将当前方法扩展到不同的管道几何形状和磁场配置，例如空间变化的磁场，以增强其在更现实毯层布局中的适用性。