自研软件介绍| HFEM 含随机夹杂的传热问题的建模与仿真

关键词：飞机结冰，水滴轨迹，结冰模型，朗缪尔，龙格库塔法，人工降雨，流体力学Ps：昨天版本有些问题，修改重新发一下科学与技术的发展，是无数人智慧接力的结果。所谓愚公移山，子又生孙，孙又生子；子又有子，子又有孙，一代代前仆后继，才有今日之昌明科技。本篇尝试从人物切入，介绍一下飞机结冰计算中经典的水滴轨迹求解问题。这个问题目前在国内外已不是什么难题，相关专业的研究生毕业论文大都有相关介绍，站在前人肩膀上，只要通晓原理，很快就可以编写出程序求解。但是就是这样一个问题，涉及了数学方法、力学模型、试验技术、计算机技术，多位天才接力才得以解决。我在多年前开发结冰软件的时候，曾仔细做过相关研究，因此本打算写一篇算法介绍的文章。但在准备阶段，查阅资料的时候竟发现很多有趣的东西。为此才准备转变写法，以期分享这些发现。1. 龙格自牛顿建立经典力学与微积分手段以后，大量工程问题有了统一的分析框架，即受力分析、构建平衡方程、求解。然而很多问题构建方程容易，求解可就难了。有些有解析解的，数学家们发挥聪明才智给出了答案。但是没有解析解的，该如何处理呢？如果这些问题实际生活中碰不到也就罢了，可偏偏有的问题又很重要。尤其工业革命以后，大家需要研究天体运行以解释宇宙，研究大炮的弹道以提升军事。于是离散化的数值求解方法应运而生，最初欧拉提出了欧拉法。以弹道问题为例，通过牛顿第二定律，建立考虑重力、阻力、升力、推力、科里奥利力等等的控制方程后，可以知晓弹体的速度控制方程。这样弹道轨迹就变成了一个已知初值（初始位置）的预测未来的问题：如何计算出它下一秒在哪个位置，下下秒又在哪里。欧拉给出来的思路通俗易懂：把整个时间历程分成多步求解，假设这一秒（步）内的速度是恒定不变的，那么根据速度矢量，乘以单步的时间，再和初始位置矢量相加，就能预测下一秒的位置。这种思想（将时间分割为小步长，用当前点的切线斜率近似下一时刻的值）简单是简单了，可是精度很差，如盲人开车，容易越开越偏。自欧拉后，又过了将近150年，德国的数学家、物理学家卡尔·龙格决定解决这个问题。当时的德国数学界由克莱因等人推动的“纯数学至上”风气主导，对应用性的数学研究不够重视。龙格当时已经从纯数学研究转向光谱分析，欧拉法的一阶精度导致误差累积过快，无法满足光谱数据对精度的苛刻要求。卡尔·龙格 Carl Runge，1856 – 1927龙格的思路是：既然欧拉法是基于“当前”的斜率往前推一步，越走越偏，属于一锤子买卖。那么如果我假装往前迈个半步，然后用那个半步的斜率再回来修我当前斜率，如此来来回 回，四渡赤水呢？结果还真行，龙格的方法综合多个局部信息加权平均，得到比单点外推更精确的结果。龙格这种多阶段试探+综合修正的思想竟然暗合兵法，神似四渡赤水。不直接突围，而是将战略转移拆成多次渡河机动，每次渡河都是一次试探和一次修正，获取敌军部署、地形和时机信息，综合四次机动后的全局态势，确定战略方向，最终选择最优路线，避免了硬冲导致的被围歼风险。2. 库塔龙格的方法发表不久后，引起了马丁·威廉·库塔的注意。库塔当时正在攻读博士论文，其方向就是微分方程的数值解法。在龙格的基础上，库塔将这种“多点试探”的方法进行了高度的数学泛化。他不仅将方法扩展到了方程组，更通过严密的数学证明，推导出了今天工程界奉为圭臬的经典四阶龙格-库塔法（RK4）的系数矩阵。龙格库塔法在工程中应用极广泛，作为求解常微分方程初值问题的黄金标准，其应用场景几乎覆盖了所有需要随时间演化模拟的科学与工程领域。从弹道轨迹到病毒传播，从电路瞬态到游戏物理，它无处不在。值得一提的是，库塔不仅仅是数学家，也是空气动力学名宿。航空史上有著名的“库塔-儒可夫斯基定理”（解释了机翼是如何产生升力的），库塔就是其中的关键贡献者。库塔所处的时代正是飞机刚刚诞生（1903）的时候，数学、力学与航空工程互相融合促进。库塔一定想不到，他完成的龙格库塔法 会在半个多世纪以后用在飞机结冰上。马丁·威廉·库塔 Martin Wilhelm Kutta，1867 – 19443. 朗缪尔3.1 粒子调控就在库塔完善了龙格库塔法的同时，远在美国的20岁年轻人欧文·朗缪尔正在潜心电灯泡的研究。欧文·朗缪尔 Irving Langmuir 1881 – 1957做科研的人大都有自己的三板斧，不管问题是什么，都有相对固定的几个套路去处理。朗缪尔是研究物理化学的，两者的交叉天然的就会出现在微观粒子上。朗缪尔研究灯泡也是从这个角度切入，他发现，通过将灯泡填充以惰性气体(如氩气)可以显著延长钨丝的使用寿命。随后，这种在一个系统里面引入额外“粒子”观察反应，成了朗缪尔研究的主要手段。自此，他开始涉足表面化学领域，并发现当分子氢被引入到钨丝灯泡中时，会分解为原子氢并在灯泡表面形成一层厚度仅为一个原子的物质，这相当于早期的“原子制造”。单纯引入额外粒子已经不能满足他，还必须要让它们“听话”。继续沿着这个思路，朗缪尔又研究了等离子体，并成为等离子体的命名人。终于在1932年，由于表面化学的工作，朗缪尔被授予诺贝尔奖。有趣的是，朗缪尔荣获诺贝尔奖的工作是关于油膜的，通过控制界面分子排布来调控宏观界面性质，把粒子层面的应用算是玩明白了。3.2 数学模型与试验二次世界大战促使航空飞行器技术快速发展，轰炸机、运输机尺寸越来越大。在著名的驼峰航线，美军动用大量C-46/C-47运输机飞跃喜马拉雅山向中国重庆抗战后方运输物资。C-46运输机该航线气象条件恶劣，结冰事故频发。根据库塔-儒可夫斯基定理，飞机升力和阻力大小受到外形影响极大，机翼形状一定设计好，最忌“整形”。然而飞机机翼结冰后，外形发生改变，升力下降，阻力增加，造成飞机失速迎角提前，再加复杂地形和气流影响，因此事故不断。现在我们知道，为了应对机翼结冰，需要空气动力学设计和除防冰系统协同。当时美军首先想到的是如何让机翼表面不结冰，即防冰。这一下撞到了朗缪尔的研究领域，“通过控制界面分子排布来调控宏观界面性质”，如果在飞机表面镀膜呢？我想最初军方和朗缪尔的思路可能是这样。就这样朗缪尔开始了飞机结冰的研究。但是对于当时的朗缪尔来说，还有个更关键的问题要搞明白。高空之中，尤其巡航高度，温度通常在零下数十度。在地面这个温度环境中，水都结成了冰，冰撞在飞机上要么弹开要么破碎，怎会积聚呢？除非，它们在零下温度情况下，仍未结冰。虽然已年逾六十，朗缪尔研究手段和物理直觉仍旧狠辣，经过对高空云层的研究，他很快确定，在没有外来粒子（灰尘等）作为凝结核的情况下，云层中漂浮的水滴果然以液态形式存在，这种低于0℃的水滴被称作过冷水滴。飞机机翼表面刚好也是提供了凝结核，才产生了结冰。随后朗缪尔分两路推进研究，一路是理论。根据前面提到的经典力学的分析框架：受力分析、构建平衡方程、求解。朗缪尔把水滴等效成球体，同时考虑流场作用，由于机翼前缘等局部位置可以看成是圆柱的一部分，为此经过简化和处理，建立了圆柱体流场中水滴运动的控制方程。随后又通过一系列推导，由于圆柱体流场求解相对简易，在不借助龙格库塔法的基础上，终于实现了水滴轨迹的艰难求解。另一路则是试验，水滴轨迹和其尺寸关系极大。可是云层中的水滴尺寸究竟是多少？它们又是否一样大，如果不一样大，又有何分布规律？科研就是如此，一环扣一环。朗缪尔根据他推导的圆柱体绕流流场水滴轨迹规律发现，不同直径圆柱体获取的水滴收集量不同，如果设计一个含有多个直径的圆柱体进行测试，就可以根据结冰量的差异反推出云层里的水滴数据。这就是旋转多圆柱法。朗缪尔通过在华盛顿山结冰观测站的大量实验，归纳出了著名的 Langmuir A 到 G 等一系列水滴尺寸分布谱（其中 Langmuir D 谱最为经典）。直到今天，民航局的适航认证标准（如 FAA 的 Appendix C）以及各大结冰仿真软件（LEWICE, FENSAP-ICE）在设置云层过冷水滴的水动力学直径分布（MVD）时，依然在使用 Langmuir 分布 作为标准输入之一。旋转多圆柱法3.3 “科学狂人”经此一役，朗缪尔对大气科学产生了极大兴趣。他开始思考，高空中的过冷水滴没有依附的微粒，就无法打破表面张力结成冰晶，也就无法变重掉落。如果人为向云层提供凝结核，或者强行让它降温产生冰晶，就能触发连锁反应，逼迫整片云下雨。这就是人工降雨的思维起点了。为了验证这个设想，他带领团队开启了试验。引入外来粒子，这是朗缪尔从博士就用惯的套路。朗缪尔助手文森特·谢弗在实验室里用一个改装的家用冷柜来模拟高空低温云层。他在冷柜里吹气，让呼出的水分在低温下形成一片微型的“过冷云雾”。谢弗尝试往里面扔各种粉末（粉尘、滑石粉、土壤）充当凝结核，但过冷云雾毫无反应。直到 1946 年 7 月的一天，纽约天气酷热，实验室里的冷柜温度不够低。为了让冷柜快速降温，谢弗顺手把一大块干冰（固态二氧化碳，温度为 -78.5°C）扔进了冷柜。神迹发生了。就在干冰进入冷柜的瞬间，原本雾蒙蒙的过冷水滴瞬间消失，取而代之的是无数的微小冰晶在光线下闪烁。干冰强大的局部降温能力，强行把局部空气降到了 -40°C 以下。在这种极寒下，过冷水滴不需要任何凝结核，自身就会自发结晶（自发成核）。这时候团队意识到，飞机结冰的原因可能不止凝结核那么简单。过冷水滴本就不稳定，水滴与机翼的碰撞会改变张力状态，同样可以造成结冰的发生。后来朗缪尔的另一位同事伯纳德·冯内古特从化学文献中找到了更优方案。他发现碘化银（AgI）的晶体结构与冰相似，可以作为"冰核"让过冷水滴在其表面冻结，且不需要干冰那样的极低温。不管怎么说，这个冰箱里的试验验证了原理可行性。随后朗缪尔与美军开始策划真正人工降雨试验。1946年11月13日，谢弗乘坐一架小型飞机，从纽约州斯克内克塔迪起飞，在附近寻找到了一片过冷的层云。他从飞机上向约3英里长的云带播撒了约3磅（约1.4公斤）碎干冰，干冰撒入后不到5分钟，云滴就转化为冰晶，几分钟后，雨雪从云层底部倾泻而下，试验成功。现代无人机增雨1947年，在美军的支持下，朗缪尔主持了著名的“卷云计划”。他甚至曾疯狂地派遣 B-17 轰炸机飞入一个正在肆虐的大西洋飓风中，倾倒了 180 磅干冰，试图用人工降雨人为改变飓风的能量分布。虽然那次试验导致飓风诡异地大拐弯并席卷了佐治亚州，尽管有人认为是巧合，但是也引发了巨大的法律和政 治风波，同时也证明了朗缪尔这套理论对大气的恐怖影响力。正所谓“天有不测风云”，自古以来天气从来不以人的意志为转移。但是作为战争中非常关键的要素，任何作战计划都必须考虑天气。近年一直传闻美军已具备相当水平的气象武器，从朗缪尔的气象试验中可以看出，这些传闻绝非空穴来风。真是蝴蝶效应，本身只是一项飞机结冰研究，最后竟然演变成了气象武器的雏形。4. NACA与NASA与此同时，NACA（NASA的前身）在俄亥俄州组建了Lewis结冰研究团队，Rinaldo J. Brun、William Lewis等多位科学家，将朗缪尔圆柱流场结果扩展到翼型流场上。他们发表了多篇具有里程碑意义的 NACA 技术报告（如著名的 NACA Report 1215）。他们利用当时极其原始的模拟计算机（如机械微分分析仪）或者手算，画出了水滴在 NACA 0012 等经典翼型周围的运动轨迹。又过去近30年，此时CFD、数值模拟、计算机技术蓬勃发展，已经具备了大规模求解水滴轨迹的基本能力。NASA组建了 LEWICE 软件团队，试图真正把拉格朗日轨迹法、面元法与现代网格技术完美结合并写成工程软件。所谓的拉格朗日轨迹法就是前面的朗缪尔的粒子追踪技术。根据相对原理，假定机翼不动，远方的水滴粒子撞击过来。根据粒子受力平衡，创建微分方程组。再通过龙格库塔法时间推进求解。水滴在流场中，会受到流场速度的扰动。为简化，采用势流作为控制方程，面元法作为求解技术。面元法的基本思路是在翼面每个单元中心（控制点）布置源强，流场中的每个点的速度由这些单元和对应的源强叠加确定。面元法面元法有个很大好处，就是一旦求解出来翼型表面的源分布，空间中任意一个点的流速，都可以通过一个矩阵乘法完成计算。如果采用更复杂的CFD技术，只能获取流场网格点的速度，这样一来，需要时刻判断水滴所处单元位置，再通过单元点的速度进行水滴位置插值获取速度，如此计算量大增。LEWICE软件于80年代问世，90年代逐步公开，该软件通过一系列精妙设计，可以在数秒内完成冰形计算。LEWICE软件5. 后续我在《我是如何开发飞机结冰软件的》中详述了自己研发飞机结冰软件的历程。人做一个事情前，结果是未知的。越是难、越是耗时的事情，越需要发大愿行之。尽管前人已经做了诸多工作，很多资料也都详实完备，轮到自己而为，仍然走了很多弯路。在多年以后重新整理这个研究的历史，故事之曲折精彩令我瞠目结舌。随着国内民机的发展，属于我们的“结冰”传奇，或许才刚刚开始。来源：静界有限元