简单三维面元法的实现
很多物理现象或者工程方法，二维的时候，我们经常能够想到很多巧妙方法来处理。而一旦我们把问题转向三维之后，往往就发现原来的方法失效了，甚至是错误的。二维中行之有效，全面封闭所有可能性的奇思妙想，面对三维之后，处处都是漏洞。这些漏洞单靠简单的升维，或者东拼西补，竟完全堵不完。
记得《时间简史》里面提到的一个想法，这个想法是尝试解释为啥生物是三维的存在，而不是二维的存在。霍金的想法是这样的，说一个生物通常要具备一个贯穿通道，用来吸收能量，排泄废物。如果这个生物是三维，比如一个球形，加了一个贯穿通道后，还是一个整体。如果是二维，加了贯穿后，这就变成了两个部分，再不是一个整体。
不得不说，大科学家就是大科学家，多有趣的想法。这也间接反映了二维和三维的差别可能远超乎我们的想象，这之间的界限甚至可以决定生死。
二维到三维，难度可能是指数级的增加。比如二维本构，刚度矩阵9个参数，到三维就是36个参数。我们之前解释的Delaunay网格生成方法，二维的时候，考虑一个圆的位置关系，勉强可以定义清楚。三维之后，上下左右，前前后后的球，这个时候的位置关系变得及其复杂，也因此虽然Delaunay生成非结构网格的思路非常清晰，但是要彻底实现，难度极大。
本次聊的三维面元法，甚至是后续的三维结冰算法，要做的就是这样一个升维的工作。
1 三维面元法理论
仍旧参考徐华舫版《空气动力学基础》。思路如下：
（1）面网格离散。考虑四边形网格，以每个单元中心点作为控制点。计算并存储好每个单元的局部坐标系。
（2） 面元强度计算。对每个单元（i）进行循环，考虑第j个单元对它的扰动：
a. 考虑在第j个单元局部坐标系下，第i个单元的法向向量；
b. 考虑在第j个单元局部坐标系下，第j个单元的坐标；
c. 循环计算第j个单元四个边对第j个控制点的扰动；
d. 形成扰动矩阵K(I,J)。
e. 求解面元。
（3） 流场计算。针对给定的点，循环计算每个单元对其的扰动：
a. 考虑给定点在第j个单元下的局部坐标；
b. 循环每个单元四个边对给定点的扰动；
c. 将局部坐标系下的扰动转换到总体坐标系下。
2 一些注意事项
单看上面的思路好像和二维区别不大。实际上，首先每个单元的局部坐标系的处理方法就很不一样。
另外循环计算每个面元时，循环的方向是要特别注意的，搞不好就反了。
最后是程序细节的问题，比如对分母可能出现的为0情况的处理，对反tan函数，可能出现的无限的处理。对判断一个值是否为0是，是否考虑容差的处理等等。任何一个细节，都会导致出来的结果有问题。
3 一些结果的展示
局部坐标系的显示
流场的显示
4 程序调试的一点思考
写了这么多程序，我的感觉，一个程序能不能成，关键是调试方法的好坏。
找到合适的算例，找到合适的判据，根据判据来定位问题，然后各个击破。
这次调试这个三维面元法，我先用了平直机翼，然后后掠翼，再用球。我的判据主要是面元强度的分布，这个要感谢前面的二维面元法，正是现有了二维面元法，能够快速得到二维面元分布值，然后才能这么快的帮我定位出各种程序bug。
做仿真找我们！
欢迎大家关注公众号:320科技工作室
