【游戏流体力学基础及Unity代码（十一）】理想流体机翼绕流和升力原理
势函数和流函数质量越大的物体，下落的距离越长，具有的势能就越大，想阻止这个物体继续下落就需要越大的能量。速度越快的流体，流过的路径的越长，具有的势能也越大。流体的势能可以使用势函数(Potential Function)来计算。势函数只是个标量场，仅在无旋流中存在，因此无旋流也叫势流。符号为phi
https://www.zhihu.com/equation?tex=%E7%9B%B4%E8%A7%92%E5%9D%90%E6%A0%87%E7%B3%BB%3Au+%3D+%5Cfrac%7B%5Cpartial+%5Cphi%7D%7B%5Cpartial+x%7D+%5Cquad+v+%3D+%5Cfrac%7B%5Cpartial+%5Cphi%7D%7B%5Cpartial+y%7D+%5Cquad+w+%3D+%5Cfrac%7B%5Cpartial+%5Cphi%7D%7B%5Cpartial+z%7D%5C%5C%E5%9C%86%E6%9F%B1%E5%9D%90%E6%A0%87%E7%B3%BB%3AV_r+%3D+%5Cfrac%7B%5Cpartial+%5Cphi%7D%7B%5Cpartial+r%7D+%5Cquad+V_%7B%5Ctheta%7D+%3D+%5Cfrac%7B1%7D%7Br%7D%5Cfrac%7B%5Cpartial+%5Cphi%7D%7B%5Cpartial+%5Ctheta%7D+%5Cquad+V_z+%3D+%5Cfrac%7B%5Cpartial+%5Cphi%7D%7B%5Cpartial+z%7D 同样，我们也可以用一个函数来表示二维速度，这就是流函数(Stream Function)
https://www.zhihu.com/equation?tex=%E7%9B%B4%E8%A7%92%E5%9D%90%E6%A0%87%E7%B3%BB%EF%BC%9Au+%3D+%5Cfrac%7B%5Cpartial+%5Cpsi%7D%7B%5Cpartial+y%7D+%5Cqquad+v+%3D+-%5Cfrac%7B%5Cpartial+%5Cpsi%7D%7B%5Cpartial+x%7D%5C%5C%E6%9E%81%E5%9D%90%E6%A0%87%E7%B3%BB%EF%BC%9AV_r+%3D+%5Cfrac%7B1%7D%7Br%7D%5Cfrac%7B%5Cpartial+%5Cpsi%7D%7B%5Cpartial+%5Ctheta%7D+%5Cqquad+V_%7B%5Ctheta%7D-%5Cfrac%7B%5Cpartial+%5Cpsi%7D%7B%5Cpartial+r%7D 有了流函数，就很容易绘制流线(Stream Line)，流线可以很好地表示速度的方向。 三维中的流函数有点特殊，本系列不会讨论。
在所有势流中，最简单的就是均匀流了(Unifrom Flow)。均匀流的流线就是一条直线，其势函数和流函数也很简单。
此处应该有一张配图，但要找一张好看的势函数/流函数图很费时间，好不容易找到一张，想微调又不可能。而且我想配上的图实在太多了，而且不少都是动图。所以本篇文章的配图主要靠咱自己画~
画图当然是用着色器了啦。代码在此clatterrr/FluidSimulationTutorialsUnity。首先计算xy方向的速度，然后是势函数和流函数。这个均匀流如果速度为正数，则是由左向右流。。首先计算xy方向的速度，然后是势函数和流函数。这个均匀流如果速度为正数，则是由左向右流。
Vx += UniformFlowSpeed * cos(UniformFlowAngle);
Vy += UniformFlowSpeed * sin(UniformFlowAngle);
phi += Vx * i.uv.x + Vy * i.uv.y;
psi += Vx * i.uv.y - Vy * i.uv.x;
python和matlab展示流线很简单，着色器里要绘制也不难
float4 col = float4(1.0f, 1.0f, 1.0f, 1.0f);
if (floor(phi * 100) % 10 == 0) col -= float4(0.0f, 0.2f, 1.0f, 0.5f);//橙黄色的线是等势线
if (floor(psi * 100) % 10 == 0)col -= float4(1.0f, 1.0f 0.0f, 0.5f);//蓝色的线是流线
return col;
最后均匀流的势函数和流函数效果如下，我为了让绘出的流线和势线好看一些，将屏幕大小设置成了1000x1000。试着调节UniformFlowSpeed和UniformFlowAngle会有什么效果吧。
均匀流很简单但也很无聊。来加个Source/Sink吧。对于源和汇来说，它们的角速度为零，而径向速度定义为Gamma/2pi r，其中Lambda是源的强度，大于则说明它是源，小于零则说明它是汇。之后通过简单的推导就能算出势函数和流函数。
写成代码如下：
dx = i.uv.x - SourcePosX1;
dy = i.uv.y - SourcePosY1;
r = sqrt(dx * dx + dy * dy);
Vx += (SouceStrength1 * dx) / (2 * 3.14159 * r * r);
Vy += (SouceStrength1 * dy) / (2 * 3.14159 * r * r);
theta = atan2(dy,-dx);
phi += (SouceStrength1 / (2 * 3.14159) * log(r));
psi += (SouceStrength1 / (2 * 3.14159) * theta);
如果均匀的速度为零，将源的强度设置为正一，那么效果如下，顺着蓝色流线，无数的流体从源离开。
将均匀流的速度由0逐渐增加到1。如果谁要写时空扭曲的小说，这张动图说不定就挺合适。
上图势函数会变成这样的原因是，由于均匀流从左向右流，在源的左边速度很小，而源的右边速度很大，因此左侧势能变化小，势能线稀疏，右侧势能变化大，势能线密集。至于流线嘛，喷泉也是这么玩的。
更直观的解释，可以用彩虹颜色来代表标量值。代码也在同一个着色器中，主要部分如下：
float x = (value + 2)/4;//value的值在-2到2之间，x的值在0到1之间
if (x < 0.25)col = fixed3(0.0, 4.0 * x, 1.0);
else if (x < 0.5)col = fixed3(0.0, 1.0, 1.0 + 4.0 * (0.25 - x));
else if (x < 0.75) col = fixed3(4.0 * (x - 0.5), 1.0, 0.0);
elsecol = fixed3(1.0, 1.0 + 4.0 * (0.75 - x), 0.0);
下图就是X轴的速度，蓝色越深表示越向-2靠近，红色越深表示越向+2靠近，绿色代表零。很显然，红色圆与蓝色圆的交点就是源的位置。
如果均匀流从左向右，在左侧放置一个源，即强度大于零，在右侧放置一个汇，即强度小于零，最终效果如下：
如果把左侧的源看成磁铁的南极，右侧的汇看成磁铁的南极，那么上面的中间那块流线完全看成是磁感线。方向由左向右。对于漩涡来说，速度则是下面这样，原理是根据上篇的Bio-Sawart Law。
代码如下：
dx = i.uv.x - VortexPosX;
dy = i.uv.y - VortexPosY;
r = sqrt(dx * dx + dy * dy);
theta = atan2(dy,dx);
Vx += (VortexStrength * dy) / (2 * 3.14159 * r * r);
Vy += (- VortexStrength * dx) / (2 * 3.14159 * r * r);
phi -= VortexStrength / (2 * 3.14159) * (theta);
psi += VortexStrength / (2 * 3.14159) * log(r);
下图是只有漩涡时候的图像
同样，保持漩涡强度为一，也就是顺时针流动。让后均匀流的速度由零逐渐增加到一。
同样，漩涡上方速度高，所以势线密集。下方则相反。Doublet则由一个源和汇组成，源在上流，汇在小流，就是我们刚才模拟的那个例子。如下图
那么
这是因为
最后速度计算公式如下：
代码如下：
dx = i.uv.x - DoubletPosX;
dy = i.uv.y - DoubletPosY;
r = sqrt(dx * dx + dy * dy);
theta = atan2(dy,dx);
Vx += (-DoubletStrength * (dx * dx - dy * dy)) / (3.14159 * r * r * r * r);
Vy += (-2.0f * DoubletStrength * dx * dy) / (3.14159 * r * r * r * r);
phi += (DoubletStrength * dx) / (2 * 3.14159 * r * r);
psi += (-DoubletStrength * dy) / (2 * 3.14159 * r * r);
然后是Doublet强度设置为一，均匀流从零逐渐增大到一。
上图虽然是doublet，但可以用于模拟理想无粘性流体的圆柱绕流，圆柱并未画出来。好了，我们已经学会圆柱绕流了，现在可以开始设计机翼来制作飞机了。【真的】
升力理论
把一个静止的圆柱放在理想流体中，就像上面那样DoubletStrength设置为0.3，均匀流速度为一，那么没什么奇怪的。但如果继续增加一个漩涡，也继续将VortexStrength设置为1，那么流线最终变成下面这样
这时候流线弯曲了。看看一下X轴速度图，上面圆柱上面气流的速度要大于气流流体的速度。
这是因为圆柱在顺时针旋转，圆柱上方的流体，由于摩擦力的作用，本来就是从左向右流的流体，会被圆柱从左边甩到右边，速度更快了。而下方一开始虽然是从左向右流，但会被圆柱从右边甩到左边，速度变慢。再根据伯努利原理，速度快的流体，所在地方压力小，速度慢的流体所在地方压力大。因此上面这个圆柱的上方压力小，下方压力大，自然就会被顶起来，有了升力。虽然我们在本章讨论的是无粘性流体，但升力的产生与摩擦力紧密相连。
圆柱体旋转，同时也产生了环量。升力与密度，空气速度以及环量有关
这就是Kutta–Joukowski定理，物体上的升力与气流密度，气流速度与物体周围的环量都成正比例。
下图则是机翼的侧面图。机翼在随着飞机移动过程中，主要会受到两个力的作用，与来流方向一致的是阻力(Drag)，与来流方向垂直的是升力(Lift)。
升力产生的直接原因是机翼上部压力比下部低，为此上部气流速度要比下部快，这样也就产生了正的环量。为了让固定翼飞机的机翼能弯曲气流，就需要合适翼型加上一定的攻角。为了增大升力和减小阻力，美国航天咨询委员会(he National Advisory Committee for Aeronautics,NACA)从1938开始设计一系列翼型来探索规律。很多翼型的数据都可以在UIUC的翼型数据库查到https://m-selig.ae.illinois.edu/ads/coord_database.html，不过你也可以直接用油管博主JoshTheEngineer在github上准备好的https://github.com/jte0419/Panel_Methods/blob/master/Airfoil_DAT_Selig.zip。本篇接下来的机翼绕流用的是NACA1412。
Source Panel Method
在之前绘制势函数和流函数的时候，我们都是预先就设定好了Source/Sink/Vortex的强度，还可以使用一对Source/Sink来模拟圆柱绕流。但对于各种各样的机翼而言，则需要另外的方法。还是首先看看圆柱绕流。要在计算机上准确画个圆是不可能的，但是可以用有限的线段去近似。这就是Panel Method，也就是用多边形来代替真正的圆柱。
多边形的边叫平板(Panel)，端点(Boundary Point)的坐标为XB/YC，中心点或称控制点(Control Point)的坐标为XC/YC。
首先假设多边形的一平板上有无数源/汇，它们排布类似于下图这样，称为Source Sheet。这种排布与涡线上漩涡排布方式不同。
一条边上所有源/汇的强度都相同，并且将它们强度加起来后来强度为lambda。而对于不同边来说，强度则各不相同。而且所有边的源强度加起来要等于零。然后我们只要求出每条边上的源强度就行了。如果知道了某条边的总强度，那么某处的势函数应该用如下公式，也就是从这条边的左端点积分到右端点，s代表边上无限小的一段
但强度应该怎么求呢？对于的平板中心点来说，可以将那里的速度向量分为两个部分，分别是沿着物体表面法向量的速度Vn和沿着物体切线的速度Vt。由于流体不能穿过这个物体，所以Vn一定等于零。有了这个等式，我们就可以求出每个平板上源强度，然后再将这个强度去求切向速度Vt。要求解的方程如下，也就是均匀流的速度加上其它源/汇的速度影响
上式要写成代码需要一点技巧，想直接看结论可以往后翻到写代码的地方。令
https://www.zhihu.com/equation?tex=I_%7Bi%2Cj%7D+%3D+%5Cint%5Cfrac%7B%5Cpartial+%7D%7B%5Cpartial+n_i%7D%5Cln%28r_%7Bi%2Cj%7D%29ds_%7Bi%7D+%3D++%5Cint%5Cfrac%7B1%7D%7Br_%7Bi%2Cj%7D%7D%5Cfrac%7B%5Cpartial+%28r_%7Bi%2Cj%7D%29%7D%7B%5Cpartial+n_i%7Dds_%7Bi%7D%5C%5C 又
https://www.zhihu.com/equation?tex=r_%7Bi%2Cj%7D+%3D+%5Csqrt%7B%28x_i+-+x_j%29%5E2+%2B+%28y_i+-+y_j%29%5E2%7D%5C%5C 那么
https://www.zhihu.com/equation?tex=%5Cfrac%7B%5Cpartial+%28r_%7Bi%2Cj%7D%29%7D%7B%5Cpartial+n_i%7D+%3D+%5Cfrac%7B1%7D%7B2%7D%28%28x_i+-+x_j%29%5E2+%2B+%28y_i+-+y_j%29%5E2%29%5E%7B-1%2F2%7D%282%28x_i-x_j%29%5Cfrac%7Bx_i%7D%7Bn_i%7D+%2B+2%28x_i-x_j%29%5Cfrac%7By_i%7D%7Bn_i%7D%29
得
https://www.zhihu.com/equation?tex=I_%7Bi%2Cj%7D+%3D+%5Cint_j%5Cfrac%7B%28x_i-x_j%29%5Cfrac%7Bx_i%7D%7Bn_i%7D+%2B+%28x_i-x_j%29%5Cfrac%7By_i%7D%7Bn_i%7D%7D%7B%28x_i+-+x_j%29%5E2+%2B+%28y_i+-+y_j%29%5E2%7Dds_i+%3D+%5Cint_j%5Cfrac%7B%28x_i-x_j%29cos%28%5Cdelta_i%29+%2B+%28y_i-y_j%29sin%28%5Cdelta_i%29%7D%7B%28x_i+-+x_j%29%5E2+%2B+%28y_i+-+y_j%29%5E2%7Dds_i 如下图，橙色直线为平板，它的角度为phi，蓝色直线为平板向外的法向量，它的角度为delta。这里将delta转变为phi更好计算。黑色直线为x轴。左侧为攻角为零。如果有攻角，那么就成了下图右侧情况，
这一段写成代码是这样的
phi = Mathf.Atan2(dy, dx);//平板角度
if (phi < 0) phi += 2.0f * Mathf.PI;
delta = phi + Mathf.PI / 2.0f;//平板法向量
beta = delta - AngleOfAttack * Mathf.Deg2Rad;//受到攻角影响
所以
https://www.zhihu.com/equation?tex=I_%7Bi%2Cj%7D++%3D+%5Cint_j%5Cfrac%7B%28x_i-x_j%29%28-sin%28%5Cvarphi_i%29%29+%2B+%28y_i-y_j%29cos%28%5Cvarphi_i%29%7D%7B%28x_i+-+x_j%29%5E2+%2B+%28y_i+-+y_j%29%5E2%7Dds_i%5C%5C 对于上式积分符号内分式的分子，可以作如下变换，小xy是平板上某点，大XY是平板的端点，角度则是平板的角度，s则是这块平板的长度的一半。而且到最后会用中心控制点来代替xi/yi的。
https://www.zhihu.com/equation?tex=x_j+%3D+X_j+%2B+s_jcos%28%5Cvarphi_j%29+%5Cqquad+y_j+%3D+Y_j+%2B+s_jsin%28%5Cvarphi_j%29%5C%5C 那么对于上式积分符号内中分式的分子来说，
https://www.zhihu.com/equation?tex=%7B%28x_i-x_j%29%28-sin%28%5Cvarphi_i%29%29+%2B+%28y_i-y_j%29cos%28%5Cvarphi_i%29%7D%5C%5C+%3D%7B%28x_i-X_j-s_jcos%28%5Cvarphi_j%29%29%28-sin%28%5Cvarphi_i%29%29+%2B+%28y_i-Y_j+-+s_jsin%28%5Cvarphi_j%29%29cos%28%5Cvarphi_i%29%7D%5C%5C+%3D+%28x_i-X_j%29%28-sin%28%5Cvarphi_i%29%29+%2B+%28y_i-Y_j%29cos%28%5Cvarphi_i%29+%2B+s_j%28sin%28%5Cvarphi_i%29cos%28%5Cvarphi_j%29+-+cos%28%5Cvarphi_i%29sin%28%5Cvarphi_j%29%29%5C%5C+%3D+%28x_i-X_j%29%28-sin%28%5Cvarphi_i%29%29+%2B+%28y_i-Y_j%29cos%28%5Cvarphi_i%29+%2B+s_jsin%28%5Cvarphi_i+-+%5Cvarphi_j%29 继续简化这个式子
https://www.zhihu.com/equation?tex=%28x_i-X_j%29%28-sin%28%5Cvarphi_i%29%29+%2B+%28y_i-Y_j%29cos%28%5Cvarphi_i%29+%2B+s_jsin%28%5Cvarphi_i+-+%5Cvarphi_j%29+%3D+Cs_j+%2B+D%5C%5C+%E5%85%B6%E4%B8%AD%EF%BC%9AC+%3D+sin%28%5Cvarphi_i+-+%5Cvarphi_j%29%5C%5C+D+%3D+%28x_i-X_j%29%28-sin%28%5Cvarphi_i%29%29+%2B+%28y_i-Y_j%29cos%28%5Cvarphi_i%29%5C%5C对于分式中的分母分母
https://www.zhihu.com/equation?tex=%7B%28x_i+-+x_j%29%5E2+%2B+%28y_i+-+y_j%29%5E2%7D+%3D+x_i%5E2+-+2x_ix_j+%2B+x_j%5E2+%2B+y_i%5E2+-+2y_iy_j+%2B+y_j%5E2%5C%5C 同样用平板的端点来代替
https://www.zhihu.com/equation?tex=x_i%5E2+-+2x_i%28X_j+%2B+s_jcos%28%5Cvarphi_j%29%29+%2B+%28X_j+%2B+s_jcos%28%5Cvarphi_j%29%29%5E2+%2B+y_i%5E2+-+2y_i%28X_j+%2B+s_jsin%28%5Cvarphi_j%29%29+%2B+%28X_j+%2B+s_jsin%28%5Cvarphi_j%29%29%5E2%5C%5C+%3D+x_i%5E2+-+2x_iX_j+-+2x_is_jcos%28%5Cvarphi_j%29+%2B+X_j%5E2+%2B+2X_js_jcos%28%5Cvarphi_j%29%2Bs_j%5E2cos%5E2%28%5Cvarphi_j%29++%5C%5C%2B+y_i%5E2+-+2y_iY_j+-+2y_is_jsin%28%5Cvarphi_j%29%2BY_j%5E2%2B2Y_js_jsin%28%5Cvarphi_j%29%2Bs_j%5E2sin%5E2%28%5Cvarphi_j%29%5C%5C+%3D+s_j%5E2+%2B+2As_j+%2B+B%5C%5C+%E5%85%B6%E4%B8%AD%EF%BC%9AA+%3D+-%28x_i-X_j%29cos%28%5Cvarphi_j%29-%28y_i-Y_j%29sin%28%5Cvarphi_j%29%5C%5C+B+%3D+%28x_i-X_j%29%5E2+%2B+%28y_i-Y_j%29%5E2 最后得到如下化简的积分式子
https://www.zhihu.com/equation?tex=I_%7Bi%2Cj%7D+%3D+%5Cint_j%5Cfrac%7BCs_j+%2B+D%7D%7Bs_j%5E2+%2B+2As_j+%2B+B%7Dds_j+%3D+%5Cint_j%5Cfrac%7BCx+%2B+D%7D%7Bx%5E2+%2B+2Ax+%2B+B%7Ddx%5C%5C 计算这个积分需要用到一些技巧。令
https://www.zhihu.com/equation?tex=u%3D+x+%2B+A+%5Cquad+dx+%3D+du%5C%5C 那么上式
https://www.zhihu.com/equation?tex=%5Cint%5Cfrac%7BCx+%2B+D%7D%7Bx%5E2+%2B+2Ax+%2B+B%7Ddx+%3D+%5Cint%5Cfrac%7BCx+%2B+D%7D%7B%28x%2BA%29%5E2+%2B+E%5E2%7Ddx+%3D+%5Cint%5Cfrac%7BCu+%2B+D+-+AC%7D%7Bu%5E2+%2B+E%5E2%7Ddu+%5C%5C%3D+C%5Cint%5Cfrac%7Bu%7D%7Bu%5E2+%2B+E%5E2%7Ddu+%2B+%28D+-+AC%29%5Cint%5Cfrac%7B1%7D%7Bu%5E2+%2B+E%5E2%7Ddu 上面最后两个积分式的第一个积分结果如下
https://www.zhihu.com/equation?tex=%5Cint%5Cfrac%7Bu%7D%7Bu%5E2+%2B+E%5E2%7Ddu+%3D+%5Cfrac%7B1%7D%7BE%7Dtan%5E%7B-1%7D%28%5Cfrac%7Bu%7D%7BE%7D%29%5C%5C+%3D+%5Cfrac%7B1%7D%7BE%7Dtan%5E%7B-1%7D%28%5Cfrac%7Bx+%2B+A%7D%7BE%7D%29%7C_0%5E%7Bs_j%7D+%5C%5C+%3D%5Cfrac%7B1%7D%7BE%7D%28tan%5E%7B-1%7D%28%5Cfrac%7Bs_j+%2B+A%7D%7BE%7D%29+-+tan%5E%7B-1%7D%28%5Cfrac%7BA%7D%7BE%7D%29%29 第二个是这样计算的
https://www.zhihu.com/equation?tex=%5Cint%5Cfrac%7B1%7D%7Bu%5E2+%2B+E%5E2%7Ddu+%3D+%5Cfrac%7B1%7D%7B2%7D%5Cln%28u%5E2+%2B+E%5E2%29+%3D+%5Cfrac%7B1%7D%7B2%7D%5Cln%28u%5E2+%2B+E%5E2%29%5C%5C+%3D+%5Cfrac%7B1%7D%7B2%7D%5Cln%28%28x+%2B+2AX%2BB%29+-+%5Cln+%28B%29%29+%3D+%5Cfrac%7B1%7D%7B2%7D%5Cln%28%5Cfrac%7Bs_j%5E2+%2B+2As_j%2BB%7D%7BB%7D%29 最后得到
https://www.zhihu.com/equation?tex=I_%7Bi%2Cj%7D+%3D+%5Cfrac%7BC%7D%7B2%7D%5Cln%28%5Cfrac%7Bs_j%5E2+%2B+2As_j%2BB%7D%7BB%7D%29+%2B+%5Cfrac%7BD+-+AC%7D%7BE%7D%28tan%5E%7B-1%7D%28%5Cfrac%7Bs_j+%2B+A%7D%7BE%7D%29+-+tan%5E%7B-1%7D%28%5Cfrac%7BA%7D%7BE%7D%29%29%5C%5C+A+%3D+-%28x_i-X_j%29cos%28%5Cvarphi_j%29-%28y_i-Y_j%29sin%28%5Cvarphi_j%29%5C%5C+B+%3D+%28x_i-X_j%29%5E2+%2B+%28y_i-Y_j%29%5E2%5C%5C+C+%3D+sin%28%5Cvarphi_i+-+%5Cvarphi_j%29%5C%5C+D+%3D+-%28x_i-X_j%29sin%28%5Cvarphi_i%29+%2B+%28y_i-Y_j%29cos%28%5Cvarphi_i%29%5C%5C+E+%3D+%5Csqrt%7BB+-+A%5E2%7D 上式推导过程来源于https://www.youtube.com/watch?v=76vPudNET6U&t=324s。这个频道的博主JoshTheEngineer主要就是讲Panel Method的，并且有非常全的配套代码库。其中https://github.com/jte0419/Panel_Methods/blob/master/COMPUTE_IJ_SPM.py展示了如何把上面那一坨公式写成python代码。其中XC/YC是平板中心点的坐标，XB/YB是平板起始点的坐标，phi是平板的角度。
A= -(XC-XB)*np.cos(phi)-(YC-YB)*np.sin(phi) # A term
B= (XC-XB)**2 + (YC-YB)**2 # B term
Cn = np.sin(phi-phi) # C term (normal)
Dn = -(XC-XB)*np.sin(phi)+(YC-YB)*np.cos(phi) # D term (normal)
Ct = -np.cos(phi-phi) # C term (tangential)
Dt = (XC-XB)*np.cos(phi)+(YC-YB)*np.sin(phi) # D term (tangential)
E= np.sqrt(B-A**2) # E term
# Compute I (needed for normal velocity), Ref
term1= 0.5*Cn*np.log((S**2 + 2*A*S + B)/B) # First term in I equation
term2= ((Dn-A*Cn)/E)*(math.atan2((S+A),E)-math.atan2(A,E)) # Second term in I equation
I = term1 + term2 # Compute I integral
# Compute J (needed for tangential velocity), Ref
term1= 0.5*Ct*np.log((S**2 + 2*A*S + B)/B) # First term in I equation
term2= ((Dt-A*Ct)/E)*(math.atan2((S+A),E)-math.atan2(A,E)) # Second term in I equation
J = term1 + term2 经过上面一堆计算，法向速度最终表示如下
切向速度推导虽然有一点不同，但基本方法一样，这里就不再展示了。如果需要详细推导，可以看看我刚才推荐的油管博主JoshTheEngineer。但到目前为止仅仅是搞定了积分符号，强度则要通过矩阵的形式算出来。如果有三个平板，并且均匀流为x方向，对于三个平板来说，那么三个法向速度公式如下
写成矩阵就是下面这个形式的。注意A的对角线上的元素，也就是平板上的源对自身中心点造成的影响为0.5，乘以2pi后结果为pi。
然后就是愉快地算矩阵啦，先求A的逆矩阵，然后乘上矩阵b即可求得强度。我求逆矩阵用的是约旦消元法，不太熟悉的可参考https://www.luogu.com.cn/problem/solution/P4783或是csdn上各种博客，大概看起来就是这样的。
for(int i = 0;i < PanelNumber;i++)//求A的逆矩阵
{
float temp = Amatrix;//获取A对角线上的元素
for(int j = 0; j < PanelNumber;j++)
{
Amatrix = Amatrix/temp;
AInverse = AInverse/temp;
}
for(int k = 0; k < PanelNumber;k++)
{
if (i == k) continue;
float mul = Amatrix;//其它行对应的元素，现在要将它减成零
for(int j = 0;j < PanelNumber;j++)
{
Amatrix -= mul * Amatrix;
AInverse -= mul * AInverse;
}
}
}
算出了每块的平板的总强度，然后我们再把平板上无限的源和汇看成一个位置处在平板中心点的大源/汇，而强度就是平板的总强度，就可以方便地用着色器绘制流线图了。对应代码在clatterrr/FluidSimulationTutorialsUnit。最后用Source Panel Method模拟的圆柱绕流效果如下：
中间的那段黑线代表的模拟的圆柱体。下面则是攻角为17度的NACA1412翼型
然后计算切向速度以及压力，也很简单
for (int i = 0; i < PanelNumber; i++)
{
float addval = 0.0f;
for (int j = 0; j < PanelNumber; j++)
addval += (SourceStrength / (2.0f * Mathf.PI)) * J;
Vt = Vinf * Mathf.Sin(beta) + addval;
Cp = 1 - (Vt / Vinf) * (Vt / Vinf);
}
VortexPanelMethod
然而SourcePanelMethod可以模拟一般的绕流，但最大的缺点就是没有环量，因此根本不可能有升力。解决方法就是源/汇换成漩涡，这就是Vortex Panel Method。它的公式与代码与SourcePanelMethod非常接近，因此接下来只说不同之处，第一处不同是计算漩涡强度积分的时候。
float Cn = -Mathf.Cos(phi - phi);
float Dn = (XC - XB) * Mathf.Cos(phi) + (YC - YB) * Mathf.Sin(phi);
float Ct = Mathf.Sin(phi - phi);
float Dt = (XC - XB) * Mathf.Sin(phi) - (YC - YB) * Mathf.Cos(phi);
第二处不同是求解漩涡强度矩阵的时候，矩阵形式本来看起来是这个样子的，对角线上的元素是零。
但是根据实验观察，如果机翼的尾部上边和下边都很平滑，但上边和下边的交界点很尖的画，那么在尾部处，上下边的气流速度是一样的，这就是Kutta条件。这样要求机翼尾部的漩涡强度相加为零。按照我的代码中的方法，机翼尾部下边编号为1，顺时针旋转到尾部上部的编号为n，也就是Gamma_1 + Gamma_n = 0。改写矩阵A最后一行第一个和第一个元素为一，其它都为零。矩阵的最后一行为零。最终矩阵形式如下。
然而上面矩阵A对角线上的元素都是零，直接使用之前的求逆矩阵方法会出错。解决方法就是使用初等行变换，将每一行都下移一行，将最后一行换成第一行即可。
第三处不同是计算速度Vt的时候
for (int i = 0; i < PanelNumber; i++)
{
float addval = 0.0f;
for (int j = 0; j < PanelNumber; j++)
addval -= (VortexStrength / (2.0f * Mathf.PI)) * J;
Vt = Vinf * Mathf.Sin(beta) + addval + VortexStrength / 2.0f;
}
代码在clatterrr/FluidSimulationTutorialsUnity。Vortex Panel Method结果如下，用八个平板模拟圆柱绕流。嗯，下图确实有八个平板，只是为了让流线好看一些，我将它的Scale设置为0.1，导致最终画出来像个正方形...
尽管如此，原来的流线还是太密了，于是我更改了着色器的如下设置，让流线不那么密
if (floor(psi * 20) % 10 == 0)col -= float4(1.0f, 1.0f, 0.0f, 0.5f);//原来是floor(psi*100)%10
最后是NACA1412翼型从10度攻角逐渐增加到30度的时候。
而且原本算出的漩涡强度太大，都是几十上百，绘出的流线乱成一团。所以我在着色器里将漩涡强度缩小为原来的百分之一。而且为了让流线连续，又将漩涡强度的精度改为0.1f。
VortexStrength /= 100.0f;
VortexStrength = floor(VortexStrength * 10.0f) / 10.0f;
VortexPanelMethod和SourcePanelMethod各有优缺点，因此更好的做法是将它们结合起来。用源模拟厚度，用漩涡模拟环量。
最后，本篇中的所有方法都假设流体是理想无粘性流体，用来模拟一些粘度低而速度快的流体可以得到合理的结果。但某些时候，在流体的物体的边界处，流体受到摩擦力的影响，速度变慢，粘性的作用开始变大，形成了需要特别处理的边界层(Boundary Layer)。下篇的内容就是用边界层理论，更准确地分析各种机翼周围的压力分布，计算升力和阻力，敬请期待！
上一篇：光影帽子：【游戏流体力学基础及Unity代码（十）】漩涡和模拟二维烟雾
推荐阅读
《Fundamentals of Aerodynamics》by John D. Anderson已经推荐这本书三遍了。这次主要是上面有关升力原理的第3.13至3.20节，有关PanelMethod的3.17节，以及有关机翼绕流的第四章内容。
《Foundations of Aerodynamics_ Bases of Aerodynamic Design》by chuen yen chow里面也有关于PanelMethod的内容
油管博主JoshTheEngineer，讲授PanelMethod最好的教程。油管主页https://www.youtube.com/channel/UCQq2_Dk2MvGGqi8-eBVwRww，配套代码https://github.com/jte0419/Panel_Methods
https://github.com/shohei/simulation，这是日本教授Shohei Aoki的一个开源库，主要是用势函数和流函数在webgl上模拟流体，从简单的均匀流到复杂的卡门涡街有差不多40份不同代码，简洁易懂。
MIT的网络开放课堂上也有很多不错的流体讲义http://ocw.abu.edu.ng/search/ocwsearch.htm?q=fluids搜索引擎很难搜索到。比如下面两篇
http://web.mit.edu/16.unified/www/SPRING/fluids/Lecture_Notes/f18.pdf
https://ocw.mit.edu/courses/aeronautics-and-astronautics/16-01-unified-engineering-i-ii-iii-iv-fall-2005-spring-2006/fluid-mechanics/f17_fall.pdf
页:
[1]