Spring-ICE 结冰算法述评-（2）水滴轨迹计算
Spring-ICE系列前述两篇文章：
今天要聊的内容，既常规又不常规。
说其常规，是因为水滴轨迹的计算是所有结冰计算的基础算法，不可或缺，是必备技能。说其非常规是因为，尽管目前文献对这个计算过程描述的已经足够详细，但是还有很多的技术细节，非实际操作而不可揣摩，正是这些具体的技术细节决定着计算的准确性、效率，甚至是整个结冰算法的设计逻辑。
我这段时间感悟，任何工程的东西，越涉及到基础的和实践的，自己必须要亲自去试一试，做一做，且不能是简单的感受下，而要做出效果，做出感悟。要一竿子插到底，摸到边。
如果我做复材，我就要亲自去做一做材料成型，研究试件加工，然后试验测试，最后是结果处理。反复的想，整个过程和我仿真或者理论，哪些是对应的，哪些是不同的，这些不同会有什么影响，有些试验现象的机理是什么，能不能从理论层面找到依据。
如果我们现在的身份不是一个学生，不是面向自己的课题。而是一个工程师面对一个项目，整个项目不是一个人能够完成的，那么作为组织者，是不是给别人分个工，等着大家的结果就完事了？我想不是，即便没有精力面面俱到，那么至少应该精通某一方面，然后瞄准项目的目标，搞清楚大家的工作和目标的联系，然后对重点问题深入研究，要搞懂要把控。总而言之，必须实践，必须讲求细节。
1 水滴轨迹计算到底要干啥？
从结冰计算整个层面来说，所谓结冰就是水滴撞在机翼表面然后冻结的现象。水滴轨迹的计算，就是要搞清楚水滴会撞击到机翼表面的哪些区域，每个区域水滴撞击的量（即水滴收集量）是多少，为下一步的结冰物理计算提供输入。
从算法的层面，水滴轨迹计算需要做以下工作：
（1）控制方程建立。说白了，水滴在流场中的运动，满足什么方程。
（2）求解方程，计算轨迹。水滴在流场中运动，肯定要有流场的计算，然后是水滴运动方程的计算。
（3）水滴运动状态的判断。这个意思是说，要判断水滴在运动中是否已经撞击到了翼面，是否已经冲出了我们预设的计算边界，如果是，这个时候需要结束计算，然后存储这个结果，为水滴收集量的计算作输入。
2 水滴轨迹计算控制方程
首先是水滴运动控制方程，这个文献里都有，通常的拉格朗日法基本是下面这种形式：
这个方程，是考虑了水滴的浮力、重力、阻力、升力，从而根据牛顿第二定律建立的。说白了，就是高中的“受力分析”，然后列方程。
说着简单，要求解出来却要费点劲。好在文献里都说了，这个方程进一步化简为一阶常微分方程组，根据初值，采用龙格库塔法进行时间推进求解。这里涉及的是本科数值计算的内容。
转眼间，我们就用到了牛顿第二定律和数值计算方法。别急，还有重头戏。上面的方程中，需要用到流场速度。流场计算用什么算法？
流场算法的选择直接影响到需不需要做网格，结冰后网格需不需要写更新算法，这是结冰算法大框架的问题，十分重要。我就在这个问题上徘徊很久。最早我是准备写NS方程求解，为此我还摸索一套自己想的网格自动生成和更新方法。最后卡在NS方程上。
回过头看，当时NS方程应该是写出来了，但是网格域的问题，还有别的基本认知问题导致我认为自己走不下去了。被迫放弃了NS方程，在师兄的指导下，选择了面元法。
从现在看，面元法的选择无疑是对的，首先它足够简单，计算效率高。早期的空气动力学还管它加鳞片布源法，就是根据空气动力里“源”和“汇”的概念，在翼面上布置“源”，每个源对其他位置的流场都有一个影响，通过叠加这些影响再综合其在翼面上法向速度为0这样的边界条件，运用线性代数方法，求解这组方程，得到每个源的强度。然后再将每个源对流场中特定点影响求出来，叠加就是这个点的流场速度。
最简单的面元法，不考虑涡，两组线性代数方程，简单易上手。
除了简单，通过我上面的描述可知，面元法不需要划分流场网格，那么冰形长出来后，只需要把边界往外推进即可，不用重构网格，无疑解决了冰形更新带来的很多问题。
有人立马就要问了，那你这么算出来的结果准不准？
要回答这个问题，首先我们要搞清楚水滴的计算是两相流的问题，水滴自身的惯性使得它不会完全按照流线运动，也就是说，水滴的运动对流场是敏感的，但是不是我们想象中那么敏感，这就给我们的流场计算留出了裕度。另外一个方面，冰形的冰角位置和水滴收集量大的地方不是绝对的对应，因为水滴在翼面上会流动，这个流动对冰角位置的影响是很大的，而这个流动由结冰热力学方程确定的。也就说，即便前面流场的计算存在一定的误差，导致水滴的运动出现一定的偏差，但是最后的结冰热力学计算，可以把冰角的位置修正过来。
这个对我很有启发，作为一个涉及多个学科的平台开发，在能不能算的准的问题上，一定要首先考虑好，最后的结果有没有兜底的修正项，找到这个兜底的，前面遇到一些影响效率的东西或者很复杂的东西，可以适当舍弃。
3 水滴状态的判断
水滴状态无外乎三种：
（1）处于运动区；
（2）冲出设定区域；
（3）撞到翼面。
前两种状态容易判断，关键是如何判断水滴撞到的翼面。这个说起来也很容易，用到的是初中的代数知识。只要连接水滴当前位置和上一个时间的位置，和翼面各段逐个求交，有交点的就是撞到了。如果撞到了，还要给出来撞到的位置。
4 算法设计
程序的设计结构如下：
（1）初始条件定义。水滴的初始位置，不能太远，不能太近。远了，计算效率低；近了，算的不准。具体多少，看你经验了，多试几次看看。
（2）当前水滴位置的流场速度计算；
（3）求解水滴运动方程，得到水滴在当前时间间隔内运动的距离，然后更新水滴位置；
（4）判断水滴状态，如果冲出防线，或者撞到翼面，停止计算；
（5）如果处于运动状态，重复步骤（2）-（4）。
需要注意的是，时间推进求解水滴轨迹，理论上时间间隔越小，肯定越准，但是呢这就意味着算的慢。太粗了也不行，后面水滴接近翼面后，时间间隔大，算的不准。我们折中下，定义一个范围，当水滴接近前缘一定范围内，加密计算，其他的区域就适当粗一点。下面的图很明显的反映了这一设置特点：
5 小结
上面的方法是Spring-ICE水滴轨迹计算算法思想。有时候看文献的时候常常想，前辈们真聪明，这些控制方程明明白白，还想到求解方法。其实仔细了解下，早期的计算模型往往也很简单，我们常人去做大部分也能想到。随着时间的推移，无数学者的努力，这些模型被不断完善，不断变的复杂高深。
我们还是要从这些繁复的理论里面，探寻些朴素且基础的设计思想，必要时候还要融入自己的理解，起码做到这个算法的灵活运用。
最后 欢迎大家关注我们的微信公众号。
微信公众号：320科技工作室。
