0 引　言

基于热电材料泽贝克（Seebeck）效应的温差热电动力系统，可将热能直接转换成电能，热电转换受材料性能限制，效率有限^{[1 - 2]}，但热电转换具有无运动部件、高可靠、长寿命等优点，在许多应用领域具有竞争力，逐步被工业、民用领域广泛研究利用^{[3 - 4]}，如同位素温差发电装置、生物质燃烧利用、火电厂废热利用，汽车废热利用等^{[5 - 6]}。热电转换技术涉及热电材料、热电器件、模块和系统等一系列问题，仍存在许多亟待解决的问题，热电系统还有较大的优化空间。目前分析计算多以集总参数方法为主，或者直接给定热电器件的温度或热流边界条件，计算模型考虑因素不够全面^{[7 - 8]}。

本文基于多物理场耦合模拟方法，针对典型的三明治型热电结构，建立温差转换结构三维多物理场模型，以热电器件、液态金属加热板和水冷板等为对象，分析了液态金属加热板、水冷板材质和流动条件对温差发电输出功率影响。

1 温差热电系统模型 1.1 温差热电系统结构

参考典型三明治型温差热电系统，以其中最小单元为研究对象，基于COMSOL Multiphysics软件进行建模分析，结构如图1所示。下部为液态金属加热板，中间为热电器件，上部为水冷板。

温差热电系统主要包括液态金属加热板、热界面材料、热电器件、水冷板等。其中液态金属加热板和水冷板尺寸为100 mm×100 mm×20 mm，圆形流道截面直径为16 mm。

温差热电系统中使用的热电器件为碲化铋-方钴矿级联热电器件，该热电器件共32对热电臂，一对热电臂包含P型热电臂与N型热电臂，热电臂上部为方钴矿，下部为碲化铋材料；不同的热电臂对通过钼铜电极串联，其中N型与P型交替均匀排布；P型与N型热电臂截面均为4 mm×4 mm，高度均为14 mm，热电器件总体尺寸为45 mm×45 mm×14.5 mm，具体结构见图2。

热电器件冷端与水冷板之间的热界面材料为导热硅脂，导热系数约14 W/(m·K)，通过外压紧密连接。液态金属加热板材质为304不锈钢，热电器件热端与液态金属加热板之间使用0.5 mm厚度的石墨纸作为热界面材料。

1.2 多场耦合模型

由于热电的过程是一个持续过程，而且热电材料的热容量小，冷端热端热通量大，所以器件各部分温度动态变化可忽略，只需考虑稳态下器件的输出特性并进行优化，为进行器件的三维数值仿真，根据稳态下的能量守恒方程和电荷连续方程：

$\nabla \overrightarrow{q}=\stackrel{.}{q}，$

(1)

$\nabla \overrightarrow j = 0 。$

(2)

以及热电耦合本构方程：

$\overrightarrow q = \left[ \alpha \right]T\overrightarrow j - \left[ \kappa \right]\nabla T，$

(3)

$\overrightarrow j = - \left[ \sigma \right]\nabla \varphi - \left[ \sigma \right]\left[ \alpha \right]\nabla T 。$

(4)

可以得到描述温度和电势分布以及热电耦合效应的控制方程：

$\nabla \left( {\left[{\boldsymbol{ \kappa }}\right]\nabla T} \right) + \frac{{\overrightarrow j }}{\sigma } - T\overrightarrow j \left( {\frac{{\partial {\boldsymbol{\alpha}} }}{{\partial T}}} \right)\nabla T - T\overrightarrow j {\left( {\nabla{\boldsymbol{ \alpha }}} \right)_T} = 0 ，$

(5)

$\nabla \left( {\left[ {\boldsymbol{\sigma}} \right]\nabla \varphi {\text{ + }}\left[ {\boldsymbol{\sigma }}\right]\left[ {\boldsymbol{\alpha}} \right]\nabla T} \right) = 0。$

(6)

式中：$\overrightarrow q$为热通量矢量；$\mathop q\limits^.$为单位体积的热生成率；$\overrightarrow j$为电流密度矢量；[α]为Seebeck系数矩阵；[κ]为热导率系数矩阵；[σ]为电导率系数矩阵。式(5)左边的4项分别为：傅里叶热传导、电流流经热电臂时产生的焦耳热、汤姆逊热和帕尔贴热。

当给定电流值时，通过求解上述方程，可以得到电势分布和温度分布，再利用数学运算得到输出电压V、输出功率P、热端吸热量Q_h。进而通过下述公式计算热电转换效率η。通过改变回路中电流值，即可获得I-P、I-η曲线，通过Nelder-Mead法优化电流即可获得P_max和η_max。

$\eta {\text{ = }}\frac{P}{{{Q_h}}} = \frac{{IV}}{{{Q_h}}} \times 100\%。$

(7)

由于热电材料性能参数具有温度依赖性，上述控制方程具有强非线性，且对于一些特殊的应用环境其边界条件是关于温度和位置的复杂函数。

基于有限元法（Finite Element Method, FEM）对上述控制方程进行数值求解，利用变分原理，将温度和电势离散为：

$T = \overrightarrow N \overrightarrow {{T_e}}，$

(8)

$\varphi {\text{ = }}\overrightarrow N \overrightarrow {{\varphi _e}}。$

(9)

式中：$\overrightarrow N$为单元形状函数矢量；$\overrightarrow {{T_e}}$为节点温度矢量；$\overrightarrow {{\varphi _e}}$为节点电势矢量；相应的热电耦合有限元方程可离散为：

$\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}^{{{TT}}}}}&0 \\ {{{\boldsymbol{K}}^{\varphi {\rm T}}}}&{{{\boldsymbol{K}}^{\varphi \varphi }}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {{T_e}} \\ {{\varphi _e}} \end{array}} \right\} = \left\{ {\begin{array}{*{20}{c}} {{Q^L}} \\ {{I^L}} \end{array}} \right\}。$

(10)

式中：${{\boldsymbol{K}}^{\rm{TT}}}$为热刚度矩阵；${K^{\varphi \varphi }}$为电刚度矩阵；${{\boldsymbol{K}}^{\varphi {\rm T}}}$为Seebeck刚度矩阵；Q^L为热载荷矢量；I^L为电流载荷矢量。其中，

${{\boldsymbol{K}}^{TT}} = \int {_V\nabla \overrightarrow N \left[ \kappa \right]} \nabla \overrightarrow N {\mathrm{d}}V ，$

(11)

${{\boldsymbol{K}}^{\varphi \varphi }} = \int {_V\nabla \overrightarrow N \left[ \sigma \right]} \nabla \overrightarrow N {\rm d}V，$

(12)

${{\boldsymbol{K}}^{\varphi T}} = \int {_V\nabla \overrightarrow N \left[ \sigma \right]} \left[ \alpha \right]\nabla \overrightarrow N {\rm d}V。$

(13)

基于COMSOL Multiphysics 6.0的稳态建模、热电效应多物理场对温差热电系统的三维多物理场耦合进行仿真，求解不同工况下热电器件的输出特性，具体计算原理流程如图3所示。

1.3 计算网格

基于以上几何拓扑及不同的物理域特点，在COMSOL Multiphysics软件中，采用结构网格和非结构网格相结合的策略对计算域进行离散。如图4所示，流体区域使用结构网格，加热工质和冷却工质使用六面体网格，通过控制分布，在流体进出口处形成结构化的面网格，进一步使用扫略，完成流体区域的网格初步划分，最后添加边界层，以更好地适应流体区域的计算。液态金属加热板与水冷板，则使用四面体网格进行划分，控制网格的大小及角细化范围，得到用于计算的网格。

热电器件中，热电材料域使用控制分布的四边形面网格沿着高度方向生成结构化的六面体网格，与热电材料接触的电极区域，则先生成面网格，进一步根据面网格扫略生成棱柱形体网格。考虑热电器件与液态金属加热板和水冷板的部件尺寸相差较大，热电器件与液态金属加热板、水冷板使用装配形式网格，避免器件电极、极板网格造成液态金属加热板、水冷板网格局部尺寸过小，降低网格质量。最终计算使用的网格约包含100万单元，最小网格质量0.21，满足计算需求。

1.4 模型设置

在温差热电系统中，主要考虑流固传热耦合及热电材料热电耦合2种耦合过程，其中热电器件主要考虑三维热电耦合，温差热电系统其他部分考虑流固耦合传热。针对流体域，在给定进出口参数后，考虑弯管形式管道，进一步选择RANS k-ε湍流模型。在热电器件对应的固体域部分，考虑固体传热、电磁、电场等物理特性。

多物理场仿真计算中的具体模型及边界条件如下：

1）液态金属和冷却介质区域考虑流场变化影响传热，使用单相RNAS k-ε湍流模型，设置流体域的进出口参数，选择速度入口条件，设定初始速度为0；

2）液态金属加热板、液态金属、水冷板、冷却介质及热电器件使用固体和流体传热物理模型，设定液态金属入口温度和冷却介质入口温度，即在耦合过程忽略热电器件与液态金属加热板和水冷板之间界面热阻；

3）热电器件域使用热电效应模型，由于固体传热过程已经在前述的固体和流体传热模型中实现，此处删除热电效应多物理场带入的固体传热物理场。

热电器件中存在诸多界面，如热电材料与电极的接触界面，热电材料与热电材料的接触界面，这些界面相对材料本身，会引入额外的接触电阻，因此需通过热阻模型进行补偿；在热电效应中还包含电流场，因热电器件陶瓷极板不导电，所以电流场区域相对固体传热物理场区域范围需要更新调整，以减少不必要的耦合计算。

与固体传热过程类似，器件内的不同材料的接触界面也存在接触电阻，需要通过设置电接触模型，考虑对应的界面接触电阻的影响。

新型热电器件中的热电材料物性随着温度会发生较大的变化，此处通过COMSOL Multiphysics 软件中函数插值功能考虑。

2 计算结果分析 2.1 热电器件最优电流

热电器件最优工作电流分析时，固定热电器件热端和冷端温度分别为T_heater=550℃，T_cooler=30℃，通过探针变量读取不同工作电流下热电器件的负载电压，结合电流与电压计算输出电功率，通过探针变量读取热电器件热端界面的输入热流积分值，得到固定温场下热电器件的输入热功率值；通过热电器件输出电功率与输入热功率得到转换效率。以转换效率为目标函数，以负载电流为控制变量和参数，开展最优效率负载电流优化分析，优化计算结果如图5和图6所示。在负载电流为5.05 A条件下，取得最优转换效率10.75%，在最优转换效率工况负载电压为2.82 V。

2.2 水冷板材料设置的影响

确定热电器件的最优效率负载电流后，分析不同温差热电系统水冷板材料对传热效果与热电器件的输出功率的影响，计算所得温度和电势分布结果如图7所示。在水冷板材质为铜时，热电器件热端为491.6℃，冷端为36.08℃，加热工质进出口温度为499.85℃/497.11℃，降低2.74℃，冷却水进出口温度分别为29.85℃/30.95℃，温升为1.10℃，热电器件输出电功率约为14.74 W，系统冷端传热功率约119 W。同样加热工质入口和冷却介质入口速度和温度条件下，将水冷板材料更换为铝合金后，热电器件热端为491.07℃，冷端为37.07℃，加热工质进出口温度为499.85℃/497.12℃，降低2.73℃，冷却水进出口温度分别为29.851℃/30.952℃，温升为1.101℃，热电器件输出电功率约为14.737 W，系统冷端传热功率约为119 W。

经过对比可以发现，水冷板材质对整个传热过程影响较小，根据温度分布可知，系统主要热阻在热电器件上，水冷板材质不同，在计算给定情况下仅导致系统温度分布产生轻微差别，并未对传热和热电转换产生明显影响。

2.3 不同入口条件下温差热电系统性能分析

在不锈钢液态金属加热板和铜水冷板材质设置条件下，分别改变液态金属入口速度和冷却工质入口速度，得到其热电器件输出电功率和转换效率随冷、热介质入口速度的变化。计算结果如图8所示，输出功率和效率随着入口速度增加先快速上升，后上升速度减缓，并趋于饱和。在变化平缓的范围内，冷水介质或者液态金属产生的热阻可以忽略，无需继续优化对应的部分。

3 结　语

本文针对温差热电系统，基于COMSOL Multiphysics 软件，使用多物理场场耦合的计算分析方法，重点开展碲化铋-方钴矿级联器件的最优工作电流分析，温差热电系统材料设置分析和冷热两端工质流动情况的分析。最终得到结论如下：

1）考虑热电材料变物性影响，不同热电材料、电极之间的电阻、热阻因素，使用热电效应多物理场耦合分析方法，设定热电器件冷热两端温度，分析计算最优热电转换效率对应的负载电流，得到该热电器件在550℃/25℃温差条件下最优转换效率可达10.75%，对应负载电流为5.05 A，负载电压为2.82 V。

2）进一步纳入液态金属加热板、加热工质、水冷板、冷却工质，使用固体和流体流动和传热耦合模型，分析不同水冷板材料对传热和热电过程的影响。经计算发现，水冷板材质对传热过程影响较弱，铝材质或铜材质对整体结构的传热功率、电功率基本无影响。

3）冷却工质、加热工质入口流速变化对热电器件输出电功率及转换效率在达到一定速度后，因流体介质热阻相对有限，入口速度增长引起的热电器件输出和转换效率变化基本可以忽略。