我与传热学的50年！
来源：《APPL THERM ENG》
作者：陶文铨 教授
1957年秋，也就是中国交通大学从上海迁往西安的第二年，我考入了交通大学。1962年夏，我从西安交通大学动力机械制造工程系毕业。我很幸运地被录取为研究生，在杨士敏教授的指导下学习热传递(HT)，杨教授于1953年，在美国芝加哥伊利诺伊理工学院雅各布教授指导下获得博士学位。1966年底，我获得了研究生毕业证书，并加入了西安交大的教师队伍。1980年10月至1982年12月，我有幸作为访问学者在美国明尼苏达大学传热实验室与E. M. Sparrow教授一起工作。我于1983年初回国，此后一直在西安交大工作。作为HT和流体流动(FF)的教授和研究员，我在这一学科工作了50多年。我的主要研究方向是数值传热和增强对流传热。 
对流项的离散化方法是有限体积法中的一个重要问题。稳定是主要挑战之一。对于稳定格式，无论网格步长和速度有多大，数值解都是稳定的，而条件稳定格式可能导致数值解振荡。到1987年为止，已经提出了几种分析方案稳定性的方法。然而，没有人可以应用于分析所有现有的方案。基于物理考虑，提出了符号保存原则。其基本思想是:利用所研究的方法对一维瞬态扩散对流方程的对流项进行离散化，扩散项采用中心差分(CD)进行离散化。从物理角度考虑，下一时刻点i和时刻点n的原始扰动在相邻两个点处的诱导值，记为 
，记为，应满足以下条件: 
的符号应与原始扰动相同[1]。而使上述不等式成立的条件可以确定该格式的临界佩莱特数，超过该格式的数值解将振荡。这是一个非常简单且具有物理意义的要求。利用这一原理可以充分分析所有现有方案的稳定性，无一例外。 在上世纪末本世纪初，我们开始关注方案的构建，并提出了SCSD[2]和SGSD[3]，后者具有二阶精度且无条件稳定，伪扩散影响较小。在格式的开发过程中，提出了一个问题:现有的对流项格式都是单独提出的。在FVM框架下是否有统一的方法来构造方案?在文献[4,5]中成功地解决了这个问题，提出了FVM中构造方案的一般方法，现有的所有方案都可以由该方法导出。
当时方案构建的另一个问题是方案能否同时具有绝对稳定性和高精度。当时的主流观点是，一个方案有两个特点是矛盾的，低阶格式可以是绝对稳定的，而高阶格式只能是条件稳定的。但事实并非如此。 文献[4,5]中已经证明，一旦一种方案的插值系数能够满足某些条件，则一种方案可以同时是绝对稳定的，且精度大于2阶。 有界性是方案构造中的另一个问题。如果数值解不超过或低于物理过程固有的控制值，则所采用的对流格式称为具有对流有界性。较为著名的判据是对流有界条件(CBC)的Gaskell-Lau判据[6]。这个CBC被认为是充分和必要的。[7]研究首先指出这种CBC是不必要的，提出了一种新的有界方案，该方案不满足G-L CBC所需的条件。后来的例子[8]表明，满足G-L CBC的方案可能导致物理上不合理的解，并提出了一种改进判据。有趣的是，其他作者在没有意识到这种改进的CBC的情况下提出的一些新方案自动实现了这种改进的CBC。现有的所有保有界方案都满足这种改进的CBC。 对于不可压缩流体流动的数值求解，SIMPLE算法得到了广泛的应用。它是半隐式的，因为当确定一个速度的修正时，它的邻近点的影响完全被忽略。这种半隐式特征影响了收敛速度和鲁棒性，自1972年SIMPLE算法提出以来，出现了许多变体。然而，所有这些简单族算法都不能完全克服半隐式特征。关键是在SIMPLE-family算法中，速度是通过在中间修正中添加一个小修正来更新的。CLEAR算法在文献[9,10]中提出，该算法通过动量方程重新求解更新后的速度，而不是添加一个小的修正。大大提高了算法的收敛速度和鲁棒性。后来在[11,12]中，在CLEAR算法的基础上发明了IDEAL算法，其中执行多个双重内部迭代(一次针对压力，另一次针对速度)。算法的收敛速度和鲁棒性得到了显著提高。 对于开放系统中N-S方程的数值求解，在如何确定出口边界的问题上存在争议。一个重要的观点是，出口边界应该定位在没有再循环流动的地方，否则数值解将没有意义。这种方法不仅需要较大的计算空间，而且有时难以实现。在电视屏幕从炉膛取出后的冷却过程这一实际工程问题的刺 激下，我们提出了一种处理这种情况的实用方法[13]:出口边界可以定位在需要的任何位置，垂直于出口边界的速度分量由局部质量守恒决定，平行于出口边界的速度分量由垂直于出口边界的零阶导数决定。然后用解域中的总质量守恒条件来修正出口处的局部法向速度。我们将不同数值处理的结果与试验结果进行了比较，发现所提出方法的结果与试验数据吻合较好，而其他方法的数值结果与试验数据存在质的差异。 在二维流场和传热的模拟中，经常会遇到与边界没有任何联系的孤岛位于场中。共轭法是将流体区和固体区同时求解的一种有效方法，将卖区视为一种特殊的高粘度流体。如何用这种共轭方法处理这样的孤岛，困扰了我们团队好几年。最后计算出数值细节如下[14,15]:(1)在每次迭代时，将孤岛上的u0、v0初始值设置为零;(2)通过设置动量方程主对角线元素的很大系数，使岛中的u*、v*为零;(3)通过对系数d设置很小的值，使u '、v '接近于零;(4)将岛内接近于零的速度传递到边界，每次迭代时设置扩散系数很大，采用调和平均进行界面插值。 第一类边界条件是数值模拟的最有利条件。对于第二种和第三种边界条件，可以使用所谓的附加源项方法(ASTM)，其中在边界处传递的能量被视为与边界相邻的第一个控制体的内部源。添加此附加源项后，边界相邻控制体与边界的连接可以完全切断，求解代数方程时，该边界可视为第一类边界。这是一种非常有效的数值处理方法，可以节省可观的计算时间。[16]文献中介绍了ASTM，并为其合理性提供了一个推导过程。此外，该方法在许多实际问题中得到了应用。在[17]中，ASTM用于处理制冷室不同表面之间的表面辐射；[18]中采用ASTM数值方法确定了开槽翅片效率；[19]中采用该数值方法求解了共轭单、双壳体内的自然对流。 为了克服FVM最严重的缺点，不规则域网格生成方法的研究在20世纪80年代到90年代引起了很大的研究兴趣。我们团队的工作包括以下内容。首先在教材的第2版中，专门安排了一整章详细介绍了体拟坐标(body- fitting coordinates, BFC)的方法。结果表明，在众多HT和FF数值模拟教材中，[16]中给出BFC内容最为详细。文献[20]提出了一种二维Delaunay三角形网格自动生成方法。[21]采用多曲面法生成截面均匀、发散和收敛的扭曲方形风管网格。[22]提出了一种生成各向异性三角形的方法，不同方向的拉伸量差异很大。通过所提出的方法，生成的三角形的纵横比可达2000。这种网格对于紊流边界层的流动计算是非常有用的。 沸腾传热的基本特征是蒸汽泡从表面产生、生长和离开。界面捕获法是模拟气泡行为的方法之一。在几种界面捕获方法中，流体体积法和水平集法应用最为广泛。VOF和水平集各有优缺点，相互补偿:VOF能保证质量守恒，但不能得到精确的界面曲率，而Level Set刚好相反。因此，发展联合办学成为必然趋势。但这种结合不应该是简单地把两种方法放在一起。我们提出了一种有效的组合方法VOSET[23]:用VOF来确定气泡质量，用Level Set来计算界面曲率;此外，曲率是由几何方法确定的，而不是通过求解微分方程。从而大大提高了界面捕获的效率和精度。在[24]中，VOSET也成功地扩展到3D情况。它现在在解决工程问题中得到了广泛的应用。 界面捕获法模拟沸腾传热的一个缺点是必须先验地设置蒸汽核作为其起点。只有中尺度和微观尺度的数值方法才能模拟气态气泡的生成过程。我们团队通过LBM对此进行了一些调查。[25]提出了结合共轭传热处理的三维多松弛时间(MRT)相变晶格玻尔兹曼方法，在不引入任何人为干扰的情况下模拟粗糙表面的成核。发现表面空腔形状对成核过程有显著影响。[25]中液体密度与气体密度之比为17。[26]中，通过在P-R状态方程中使用较小的系数值，将液气密度比扩展到200。当然，这一比例仍远低于实用价值。例如，大气中的水的这个比值是1600。留了一个大房间供进一步改进。此外，LBM在[27]中也研究了不凝性气体的膜状冷凝HT。 进入21世纪，多尺度模拟成为国际数值模拟界的研究热点。然而，当我们仔细研究需要多尺度模拟的案例时，我们发现不同的研究人员对多尺度问题的含义有不同的理解。从数值模拟的角度来看，我们发现有必要将多尺度问题分为两类:多尺度过程和多尺度系统[28,29]。对于多尺度过程，不同区域有不同的控制方程和数值方法。数值解应在界面处耦合。对于多尺度系统:在不同区域使用相同的控制方程和相同的数值方法。区别在于网格步长。两个典型的例子分别是质子交换膜燃料电池(PEMFC)中的传输过程和数据中心中的传热。据笔者所知，目前文献中所研究的多尺度换热问题均可分为这两类。 12.FVM/LBM/DSMC/MDS之间的耦合方法 对于多尺度过程，不同数值方法下不同区域界面处解的耦合成为一个重要问题。不同区域采用纳米尺度、微/中尺度和宏观尺度的方法，在界面处交换溶液，得到整个全溶液。相邻区域界面处的数值解(信息)交换可以用两个运算符在数学上表示: 
其中Φ和ϕ分别是宏观参数(如速度)和介观或微观参数(如LBM中的分布函数)，C和R分别是压缩和重构算子。运算符C很容易定义，比如用不同平均方法，但算子R很难构造，因为在这里，来自一个解区域的较少信息(例如FVM中某一点的速度)应该转换为邻近区域的更多信息(例如LBM中的密度分布函数)。对于FVM解与LBM解的耦合，我们的团队在[30-32]中提出了成功的重构算子。在多尺度原子连续体方法中，由于原子区域的非均匀性，周期边界条件不再适用。为了在保持MD系统正确的平均压力的同时摆脱MD模拟区域边界处的周期性条件，需要在边界附近的原子上施加一个边界力。[33]中提出了确定这种力的方程。FVM与DSMC之间的耦合研究[34]。 以编织织物为增强相，在合适的基体中制备编织复合材料，具有重量轻、强度高等优点，广泛应用于航空航天、汽车等领域。在编织复合材料等广泛应用的增强纤维材料中，碳纤维的轴向和横向热导率各向异性比较严重，其有效热导率的测定是一个很大的挑战。在适当的边界条件下，研究了三维四向编织复合材料在中观尺度(编织纱线)和宏观尺度(复合材料)上的热传导特性[35]。考虑了不同形状的编织纱线横截面。实验结果验证了该模型的正确性。在[36]中，对平纹编织复合材料的有效导热系数进行了数值预测。利用复合材料的平动对称、反射对称和旋转对称三种不同的对称构造了三种缩小尺寸的单元胞，并推导和验证了相应的热边界条件。 现在让我们对加强传热的研究作一些总结。强化传热是高温研究中一个永恒的课题。1959年春，杨 思敏教授为我们开设了HT课程，他在课程一开始就强调HT增强是HT研究的主要任务之一。板翅管式高温表面被广泛用于增强气侧高温。从几何角度看，板翅管换热面与普通风管相比的主要特点是存在壁筒。尽管已经对不同类型的板翅管高温表面进行了大量的研究，但没有人声称壁筒对高温的影响。该研究在[37]中采用萘升华法进行。研究发现，壁面圆筒起到了湍流促进剂的作用。它可以显著增强低雷诺数区域的高温。而壁面圆筒的存在也使得层流向湍流的过渡时间提前。 尽管对高温增强的研究已经进行了一个多世纪，但直到上世纪末，对流高温增强的基本机制仍不清楚。Guo提出了基于板上层流不可压缩边界层(抛物线流)的高温场协同原理(field synergy principle, FSP)[38]，认为提高速度和流体温度梯度之间的协同作用是其基本机制(由两个矢量之间的夹角表示)。后来FSP扩展到椭圆流[39]、可压缩双曲流[40]和湍流高温流[41]。此外，还表明FSP可以统一现有的增强单相对流高温的机制[42]。 对于涡发生器等复杂的增强技术，也揭示了提高速度和温度梯度之间的协同作用是其基本机理[43]。我们还通过数值算例揭示了FSP和熵理论对一些增强技术的评价结果是完全一致的[44]。FSP被用来指导新的增强技术的发展，如前稀疏-实密集狭缝鳍面[45,46]。 当空气和水是两种高温介质时，通常将空气设置在管外，在那里可以使用一个大房间来开发翅片表面。但由于结构上的要求，有时会在管外设置冷却水，使其在管内横向流动，热空气进入管内。为了增强空气侧高温，可以采用纵翅片。为了固定鳍片，通常在管中同轴地设置一根小管，并在两管的环形空间之间设置鳍片表面。这种换热器在一些压缩机中用作气体冷却器。在管道HT和摩擦系数的实验测量中，我们进行了两种变化，一种是内管堵塞，即没有空气流过，另一种是内管不堵塞。令我们惊讶的是，在一定的外径与内径之比范围内，同样流速下，中心阻塞管的高温速率明显大于未阻塞管，这就是我们所说的内翅片纵波状翅片中心阻塞管[47]。随后，我们对高温和FF进行了数值分析，发现中心阻塞管优越的高温性能是由于外管速度梯度的增加导致协同作用的明显改善。此外，我们根据自己的实验数据和文献中关于水在微鳍管中流动的实验数据，对Gnielinski方程进行了推广，确定了湍流平均高温系数[48]。 历史上，不同作者对微通道摩擦系数和努塞尔数的测量结果存在着严重的差异和偏差。随着制作方法、测量和数据采集技术的进步，造成这种偏差和差异的原因逐渐显露出来。目前人们普遍认为，除了尺寸效应(气体的稀薄性和可压缩性效应、表面粗糙度效应等)外，单相流动的流体流动和传热的常规解和关系式仍然可以预测微通道内的流体流动和传热行为。我们的团队是世界上做出这样声明的研究单位之一[49,50]。在揭示尺寸效应时，我们发现当表面相对粗糙度小于1%时，气液在微通道内的流动符合常规理论。然而，对于直径小于10-20 μm的微通道，存在稀薄效应[51]。在外壁温度恒定的直厚壁管中，二维壁面和流体轴向热传导对同时发展层流和传热的共轭传热问题的影响进行了数值研究[52]。结果表明，壁面轴向导热的基本功能是统一内壁表面热流密度。在[53]中，对水力直径为500微米的波浪形通道中的FF和HT进行了数值研究。文献[54]对多弯曲通道中1.5 mm间隙的三维抛物型湍流强迫对流HT和FF进行了数值模拟。我们团队在气体微通道流体流动和传热方面的主要成果总结于文献之中[55]。 20世纪90年代，为了取代消耗臭氧的工作物质，在世界范围内进行了广泛的研究。我们小组主要研究了152a和134a的替代品，并于1994年发表了第一篇关于152a单翅片管外冷凝实验结果的论文[56]。结果表明，在饱和温度为40◦C时，152a的冷凝HT系数比0f R-12的冷凝HT系数高约20-25%。在此基础上，建立了管内对流强制冷凝沸腾实验装置，并对备选134a进行了实验测量和理论分析。文献[57,58]分别报道了光滑水平管内134a强制对流冷凝和沸腾高温局部换热特性的实验研究。随后，四种制冷剂(包括R134a和R410a)在光滑管和微翅片管内沸腾换热的实验结果发表于[59]。 降膜蒸发器相对于池沸蒸发器具有制冷剂充注量少、体积小、传热温差小、除油容易等固有优点，被认为是冷水机组或热泵系统中池沸蒸发器的潜在替代品。尽管已经进行了大量的研究，但仍然缺乏令人满意的传热相关性，甚至没有一个公认的制冷剂水平单光滑管的传热相关性。基于我们之前的实验研究[60-62]，[63]提出了单个水平光滑管上制冷剂降膜蒸发的高温相关性。建立了五种制冷剂的两个相关关系。全润湿状态的相关性在±30%的偏差范围内拟合96.7%的我们自己的542个数据，而在-30%到+ 15%的范围内拟合73.4%的其他作者的289个数据;部分干枯状态的相关性在±30%的范围内与我们自己的97.5%的数据拟合，在±30%的范围内与其他作者共95个数据的76.8%的数据拟合。 2016年至2017年，在学校的支持下，我校重点实验室建设了微流体实验室。在此基础上进行了微流控研究，摩擦电纳米发电机(TENG)就是其中之一。它是一种革命性的能量收集技术，利用摩擦电将机械能转化为电能。为了适应各种能源和各种应用的适用性和灵活性，已经开发出了种类繁多的具有不同结构设计的teng。我们的团队也对此进行了一些调查[64]。我们提出了一种新的结构，电极和疏水层作为顶板，电极作为底板。我们的实验结果表明，这种结构优于其他器件，这也得到了我们对介电层对输出功率影响的理论分析的支持。此外，还实验研究了液滴数和振动频率对最优负载阻力的影响。 接触热阻(TCR)是航空航天工程中热防护系统设计的关键因素。在过去的几十年里，人们进行了大量的研究来确定不同材料之间的TCR。基于傅立叶热传导定律的一维稳态法是目前应用最广泛的TCR测量方法。为了确定一对材料的TCR，应测量界面温差，并给出测试材料的导热系数。为了呈现TCR的一个主要参数的测量结果，应该提供平均表面粗糙度。因此，界面温差、试样导热系数和试样表面粗糙度的测量构成了TCR的三个决定性要素。这一思路在文献[65]中已经提出，文中也提出了利用ABAQUS对TCR进行数值模拟的方法。对于[65]中所研究的Ti-6Al-4V-Ti-6Al-4V材料对，模拟结果与实验数据进行了比较，最大偏差为9.57%，75%的结果偏差在5%以内。 换热强化技术是提高换热器热性能的有力手段。然而，任何高温强化技术都会导致更大的压降，并且众所周知，压降增加的比例往往大于传热增强的比例。因此，如何对增强技术进行综合评价的问题引起了人们的关注。提出了若干评价标准。随着世界范围内能源短缺危机的出现，强化换热的节能目的变得更加重要。[66]提出了一种基于相同流量、相同压降、相同泵送功率下增强高温表面与参考高温表面对比的综合评价方法，并绘制了性能评价图。在该图中，分别以增强表面的摩擦系数与参考表面的摩擦系数之比和相同雷诺数下相关的传热增强系数之比作为横坐标和纵坐标。由两个坐标包围的第一象限可分为四个区域:在区域1中，高温增强，压降损失较大，每泵送功率相同，高温就会恶化;在相同的泵送功率下，区域2的高温增强，而在相同的压降下，区域3的高温增强，相同的压降下，区域4的高温增强比大于摩擦系数增加比。在该图中，可以很容易和清楚地比较不同的增强技术对同一参考的节能性能。 为了节约能源和减少温室气体排放，中国政府大力采取节能措施。其中一项措施是根据能源效率(EE)来决定新计划的工厂是应该批准还是应该关闭已经存在的工厂。文献[67]提出了一种基于能效指标(EEI)体系的分层指标比较(HIC)方法来评估中国不同行业的能效。HIC方法的基本思想是根据指标将电厂的EE与参考值进行比较。如果比较结果表明该装置比参考装置差，则不允许建立该装置。以化工行业的纯化对苯二甲酸(PTA)为例，说明了HIC法。此外，文献[68]提出了一种用于中国区域能源效率预测与评价的混合预测方法。在这些研究的基础上，[69]对高耗能行业能效评价的方法和政策进行了全面回顾。 24.带有太阳能烟囱和光伏系统的太阳能辅助空气净化系统 太阳能辅助空气净化系统是利用太阳能治理大气污染的一种新尝试。2016年，西安建成了一个示范单元，这是世界上第一个使用太阳能和先进过滤技术进行空气净化的建筑结构[70]。该单元有一个60米高的烟囱和一个2580平方米面积的收集器。我们团队部分参与了该示范装置的设计、测试和评估。利用太阳能产生热上升气流，吸引环境中的污染空气，污染空气被烟囱进风口前设置的过滤器净化，然后通过烟囱排出到环境中。这种设备的重要输出参数是通过烟囱的空气流量。基本上，集热器面积越大，空气流速越大。更大的空气流速是以土地为代价获得的。数值研究了由太阳能烟囱和覆盖在集热器上的光伏(PV)板组成的混合系统和仅由太阳能烟囱组成的系统。结果表明，采用40米宽的光伏板覆盖整个集热器顶部的混合系统，系统容积流量仅降低19%，但系统总输出功率提高了57倍[71]。产生的电力可以用来转动风扇，使更多的空气通过烟囱。因此，太阳能辅助空气净化系统覆盖所有集热器表面的光伏板是一个有效的变种，以节省土地。 对于城市街道污染运移的研究，可以采用数值模拟和风洞试验相结合的方法。为了进行风洞试验，出现了一个问题:如何确定雷诺数?当气流经过街道时，所谓的大气边界层(ABL)就形成了。相似理论要求，要使风洞过程与原型相似，必须满足两个条件:(1)几何相似;(2)流动动力相似性，即风洞试验的Re必须等于原型的Re。如果使用原型的1/ 500的模型，则平均速度为4米/秒。那么风洞中的速度应该是2000米/秒!这在ABL风洞中是无法控制的。我们建立了一个与物理风洞完全相似的数值风洞，在模拟中采用了三种双方程湍流模型[72]。数值结果表明，环境湍流边界层流动可分为低雷诺数范围内的强re依赖区和高雷诺数范围内的弱re依赖区。提出了相对变化比的概念来确定区分两个区域的临界雷诺数。如果原型机雷诺数处于弱Re依赖区，则其风洞试验可以在大于临界Re的任何Re下进行，而临界Re通常远小于原型机雷诺数。因此，确定临界雷诺数将为风洞试验提供极大的方便。 本课题组从2000年开始对质子交换膜燃料电池(PEMFC)的水热管理进行研究。PEMFC最重要的综合性能是其极化曲线(PC)受几十个尺寸参数的影响。为了在数值上预测这个PC，我们开发了几个模型和相应的代码[73-76]。由于受几十个参数的影响，到目前为止，极化曲线都是针对这类参数值的单独一组进行量纲化得到的，并且每条曲线只能适用于这一特定的参数组。文献[77]全面回顾了相似理论在PEMFC各组分研究中的应用，文献[78]对PEMFC三维单相等温模型采用相似分析得出相似准则。得到了7种输入判据，并将无量纲电压和无量纲电流密度定义为两种输出判据。数值验证表明，当7个判据在较大范围内保持各自的值时，无量纲极化曲线保持一致，偏差约为1%，表明了本文分析的有效性和可行性。从对无量纲极化曲线的影响来看，敏感性分析表明，7个标准可分为强、轻度至轻微和可忽略三类。相似度分析方法可以大大节省PEMFC输出特性建模的计算时间。 从上面的简要介绍可以看出，传热是热流体科学中一个非常活跃的子领域。随着世界科学技术的发展，当前的HT研究正在发生着巨大的变化。正如中国谚语所说，活到老学到老，我会跟踪它的发展，尽我所能做研究。 我有幸在杨 思敏教授和E. M. Sparrow教授的指导下学习HT;我很幸运，能和160多名可爱、勤奋的研究生一起工作;我很幸运，能与团队中许多年轻有为的同事一起，认真、和谐地教授和研究HT;此外，我们团队的发展也长期得到了国内外众多同仁的支持。借此机会，我向我的导师、我的学生和我在国内外的同事们表示衷心的感谢。我的团队的研究得到了多项中国政府基金的支持，包括国家重点研发计划和国家自然科学基金。这些基金的支持是非常感谢的。我要特别感谢牛津大学的何李教授、诺丁汉大学的毛雪芮教授、闫玉英教授以及上海交通大学的徐辉教授，感谢你们于2019年9月在诺丁汉大学组织了我80岁生日的庆祝活动，并在备受尊敬的《应用热工程》杂志上发表了这期特刊。由于演讲的限制，我做过的一些课题没有在上面列出。我感谢那些执行这些科目的学生的理解。
著作权归作者所有，欢迎分享，未经许可，不得转载
首次发布时间：2023-12-23
最近编辑：10月前
