数值计算方法在风洞翼型实验中的应用

前言

​ 夏学期选课阶段，莫得课上的我选了一堆奇怪的个性课，其中就有控制学院的数值计算方法。选中而错过退课的时候很虚，一来我没有修过常微分方程，常微分是这门课的预修课程，二来计院的相关课程数值分析在大二上开设，我还没有修过。硬着头皮上的时候发现，内容是很有意义的，作业和报告是很多的，控院的选修课是很水的…

​ Numerical Method 是研究分析用计算机求解数学计算问题的数值计算方法及其理论的学科。数值计算方法一词是从苏联教材翻译而来，在美国，被称为 numerical analysis，译作数值分析。这是一个奇妙的学科分支，它和离散数学格格不入，专门研究连续的世界，广泛应用于实际工程中，因此，数值计算方法被称为工科生都应该学一学的课程。 本文是我由数值计算方法这门课的大作业整理而来，大作业代码于6月13日完成，报告于6月15日完成，虽然非常水，但还是解决了一些小问题，回顾了一些有趣的东西。

问题背景

​ 数值计算方法经常用于飞机外形的设计，尤其是翼型。飞机的制造耗资巨大，因此，设计之初的模型实验必不可少。如何保证模型的数据可以应用到真实飞机上呢？空气动力学给出了解答，只有无量纲化之后的升力系数可用。问题因此变得简单明了，我们只需要测出相应攻角、迎角、动压下翼型各离散点的压强，对此进行曲线积分，即可求出翼型的升力系数。

问题分析

​ 飞机要实现飞行，首先依靠机翼的升力。那么升力是怎样产生的呢？这就是人们熟知的伯努利原理：水与空气等流体，流速大的地方，压强小；流体流速小的地方，压强大。把机翼纵向剖开，会形成一个翼截面或翼剖面，在航空上称翼型。当空气流过机翼时，气流会沿上下表面分开，并在后缘处汇合。上表面弯曲，气流流过时走的路程较长，下表面较平坦，气流的行程较短。上下气流最后要在一处汇合，因而上表面的气流必须速度较快，才能与下表面气流同时到达后缘。根据伯努利原理，上表面高速气流对机翼的压力较小，下表面低速气流对机翼压力较大，这就产生了一个压力差，也就是向上的升力。

​ 机翼的翼型选择是机翼设计的关键，优秀的翼型能为飞机在不同飞行状态下提供良好的升力。翼型对机翼升阻比特性有重要影响，在机翼设计的不同历史时期出现了不同类型的翼型。

​ 而笔者获得的数据是DU 93-210翼型的0.8米模型数据，根据测压节点坐标绘制翼型图如下：

理论背景

压强系数

​ 我们所用的模型是0.8米模型，尺寸显然与真实的机翼相差甚远，而如何衡量翼型的性能呢？在流体力学中，一般将压强用无量纲化的参数压强系数来表示各个测点的压强系数值：

​ 式中，分别是测点压强、来流压强、驻点压强（总压），公式由伯努利方程得来。

升力系数

​ 气流给予翼型的总合力在y轴上的分量成为升力，记作。根据资料，在紊流绕流中，粘性切应力对总合力的贡献仅占很小份额，因此，通常仅考虑压强作用，升力系数定义为：

​ 式中，是升力作用面的面积，对于二元翼型，升力的作用面等于弦长C乘以单位宽度。

计算方法

​ 设上表面的微面积ds，设该面积上的压强为p，则压力为pds，投影到y轴得，负号表示压力方向为y轴负向。对于下表面，压力应取负值，因此，升力是上下表面压力的代数和：

​ 升力系数为

​ 还可以用速度环量法计算升力系数，我们在这里不做考虑。可见，根据已有数据条件，升力系数的计算可以转化为数值积分问题。

数据说明

​ 文件会放在GitHub相应项目下。

​ 另有洞壁上的压强分布数据，根据牛顿第三定律，洞壁上下表面的压力差应该与翼型的升力非常接近，我们可以以此来验证实验。

​ 通过这些数据我们可以做两件事，一是对翼型坐标进行插值、拟合，求出轮廓曲线；二是对压强进行数值积分，求出翼型的升力系数，也可以求出升力，跟洞壁压力差比较。

曲线拟合

​ 当插值节点数目较多时，普通的拉格朗日、牛顿插值法会出现龙格现象，即在区间边缘出现的振荡现象，而采用分段线性插值时，函数光滑性又较差。我们对两种方法做出尝试，特别地，对于当前问题，分段线性插值可以简单地看作将各测压点连线，运行结果如下图。可以看出，由于翼型自身的光滑性较好，即使是线性插值肉眼也看不出光滑性的损失，而由于插值节点过多，在第十个节点左右即出现了非常明显的龙格现象，并且无法收敛。

​ 因此，为了兼顾消除龙格现象和提高函数的光滑性，我们采取分段插值，并且要求分点处若干阶导数值也相等，实际应用中常采用分段三次Hermite插值，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

function varargout=pchinterp(xd,yd,ydot,xi)
% PCHINTERP  分段三次Hermite插值
n=length(xd);  
if length(unique(xd))<n  
    error('互异')
end
if nargin==2 || isempty(ydot)
    ydot=gradient(yd,xd); 
end
h=diff(xd); 
pp=zeros(n-1,6);  
y=@(x)zeros(size(x));  
for k=1:n-1
    pp(k,1:4)=2*yd(k)/h(k)^3*poly([xd(k+1),xd(k+1),(3*xd(k)-xd(k+1))/2])+...
        -2*yd(k+1)/h(k)^3*poly([xd(k),xd(k),(3*xd(k+1)-xd(k))/2])+...
        ydot(k)/h(k)^2*poly([xd(k),xd(k+1),xd(k+1)])+...
        ydot(k+1)/h(k)^2*poly([xd(k),xd(k),xd(k+1)]);
    pp(k,5:6)=xd(k:k+1);  
    y=@(x)y(x)+polyval(pp(k,1:4),x).*(x>=xd(k) & x<xd(k+1));
end
if nargin==4
    y=feval(y,xi); 
    y(xi==xd(end))=yd(end);
end
pp=array2table(pp,'variablenames',{'Ai','Bi','Ci','Di','xi','xi_1'});
[varargout{1:2}]=deal(y,...  
    pp);

​ 在大程序中添加如下代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

xdata = wing_x(1:48);
ydata = wing_y(1:48)
xi = xdata
yi = pchinterp(xdata,ydata,[],xi);
plot(xdata,ydata,'ko',xi,yi)
axis equal
axis([-100,900,-200,200]) 
hold on
xdata = wing_x(49:96);
ydata = wing_y(49:96);
xi = xdata
yi = pchinterp(xdata,ydata,[],xi);
plot(xdata,ydata,'-o',xi,yi)
axis equal
axis([-20,820,-150,150]) 
title({'分段三次Hermite插值'})

​ 运行结果如下：

​ 可以看出，由于翼型本身光滑性极好，因此分段插值得到的曲线连续性、光滑度都极好，完全可以用于曲线积分时求导。

洞壁压力差

​ 对于洞壁测得的离散点压强数据，我们进行复化梯形积分，得出上下表面压力差。将积分区间[a,b]划分为n等份，分点为,其中,在每个小区间上利用梯形公式，则有

​ 代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

function z = trapequad(x,y,dim)
%TRAPE
    perm = []; nshifts = 0;
      dim = min(ndims(y)+1, dim);
      perm = [dim:max(ndims(y),dim) 1:dim-1];
      y = permute(y,perm);
      m = size(y,1);
    elseif nargin == 2 && isscalar(y) 
      dim = y;
      y = x;
      x = 1;

      end
      dim = min(ndims(y)+1, dim);
      perm = [dim:max(ndims(y),dim) 1:dim-1];
      y = permute(y,perm);
      m = size(y,1);
    else 
      if nargin < 2
          y = x;
          x = 1;
      end
      [y,nshifts] = shiftdim(y);
      m = size(y,1);
    end
    
    x = x(:);
    

    if isempty(perm) && isequal(y,[])
      z = zeros(1,class(y));
      return;
    end
 
    if isscalar(x)
        z = x * sum((y(1:end-1,:) + y(2:end,:)), 1)/2;
    else
        z = diff(x,1,1).' * (y(1:end-1,:) + y(2:end,:))/2;
    end
 
    siz = size(y);
    siz(1) = 1;
    z = reshape(z,[ones(1,nshifts),siz]);
    if ~isempty(perm) && ~isscalar(z)
        z = ipermute(z,perm);
    end
end

​ 也可以直接调用matlab自带的trapz函数，我们在大程序中加入下列代码：

1
2
3

sum_up_wall_pressure =trapz(up_wall,up_wall_pressure);
sum_low_wall_pressure =trapz(low_wall,low_wall_pressure);
sum_wall_pressure_difference = sum_low_wall_pressure - sum_up_wall_pressure;

​ 得到上下洞壁压力差约为-112.8941，单位为PSI·mm。

升力

​ 计算升力，即是对翼型上下表面压力求积，做法类似洞壁压力差的计算，我们先采用梯形积分，在大程序中加入下列代码：

1
2

sum_wing_pressure_difference = sum(wing_difference);
sum_wing_pressure = trapz(wing_x,wing_pressure);

​ 得到上下洞壁压力差约为120.8783，单位为PSI·mm。

​ 当然，梯形公式虽然简单易用，但是精度较低，可以推导出复化科斯特公式如下：

对于当前问题这种全为离散点的情况，意味着每个区间用于积分的函数值由其附近几个点的函数值加权得到，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

function q=cotesquad(varargin)
% COTESQUAD  复化柯特斯求积
args=varargin;  % 将输入参数赋给args
if isnumeric(args{1})  
    if mod(length(args{1}),4)~=1
        error('数据长度不满足柯特斯公式')
    end
    if nargin==1  
        y=args{1};  % 函数值
        x=1:length(y);  % 自变量
    else 
        [x,y]=deal(args{1:2}); 
    end
    h=x(5:4:end)-x(1:4:end-4); 
else
    if nargin==3
        [fun,a,b]=deal(args{1:3});  
        n=1000; 
    else
        [fun,a,b,n]=deal(args{1:4});  % 被积函数，积分区间和区间等分数
    end
    h=(b-a)/n;  % 等距节点的步长
    x=a+(0:4*n)*h/4;  % 4n+1个等距节点
    y=feval(fun,x);  % 计算节点处的函数值
end
q=sum(h/90.*(7*y(1:4:end-4)+32*y(2:4:end-3)+12*y(3:4:end-2)+...
    32*y(4:4:end-1)+7*y(5:4:end)));

​ 之前的问题中，我们的重点在于曲线插值和数值积分，并没有考虑量纲的问题。

升力系数

​ 对于刚才所求的压力，进行单位换算后，再除以动压和面积，即可得到该翼型的升力系数。根据空气动力学的理论知识，升力系数与雷诺数、迎角等相关，我得到了从-15°到25°四十一个迎角和两个不同雷诺数下的数据，数据说明如下：

原始数据：

第1行

1 -15.0 0.0 191.477 56.517 191.980 56.591 973.20 14.40 32.50 →

序号 攻角 未用 动压 风速 动压 风速 大气压 温度 湿度

（落差法）（落差法）（风速管）（风速管）（100帕）（℃）

-0.02657 -0.29635 0.00002 -0.02653 -0.29521

P1 P2 未用 风速管总压 风速管静压

（落差法稳定段压力）（落差法试验段入口压力）

（以下各行说明）

共1-96行，模型表面压力 ，第一列是序号，二、三列不用，第四列是各点压力

共1-187行，尾迹排管的总压 ，第一列是序号，二、三列不用，第四列是各点压力

共1-135行，风洞上壁面的静压 ，第一列是序号，二、三列不用，第四列是各点压力

共1-135行，风洞下壁面的静压 ，第一列是序号，二、三列不用，第四列是各点压力

共1-9行，尾迹排管的静压 ，第一列是序号，二、三列不用，第四列是各点压力

2行， 不用

（如上重复开始下一个攻角,直到全部41个迎角数据完）

注：

1、大气压的单位是百帕（100Pa），即所给值乘以100单位为pa，且大气压是绝对压力（绝压）；

2、除大气压和动压外，压力单位为PSI，1PSI=6894.76Pa；

2、动压单位为kgf/m2 , 1kgf/m2=9.8Pa

4、除大气压和动压外，其余所有的压力均为表压，即该处压力减去大气压的差值。对于直流风洞（如NF-3），根据伯奴利原理，随风速增加，风洞内的静压降低，用表压表示时即为负值；由于存在能量损失，故总压也为一个绝对值较小的负值。所有压力加上大气压才是相应的绝对压力；

5、x-actual文件中是风洞上下壁测压孔坐标；

6、翼型测压孔坐标.dat文件中是翼型表面测压点的x、y坐标

​ 对不同条件下得到的数据进行处理，加入代码段如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77

mm = 1000;
 s=0.8;
 
Wall = load('洞壁坐标.txt');  % 执行后，会得到一个N行3列的矩阵A。
wall_i = Wall(:,1);  % 这样就可以取出A的第1列。
up_wall = Wall(:, 2)/mm;
low_wall = Wall(:, 3)/mm;
 
Wing = load('翼型坐标.txt');
wing_x = Wing(:, 1)/mm;
wing_y = Wing(:, 2)/mm;
 
xdata = wing_x(1:48);
ydata = wing_y(1:48)
xi = xdata
yi = pchinterp(xdata,ydata,[],xi);
plot(xdata,ydata,'ko',xi,yi)
axis equal
axis([-0.020,0.820,-0.150,0.150]) 
hold on
xdata = wing_x(49:96);
ydata = wing_y(49:96);
xi = xdata
yi = pchinterp(xdata,ydata,[],xi);
plot(xdata,ydata,'-o',xi,yi)
axis equal
axis([-0.020,0.820,-0.150,0.150]) 
title({'分段三次Hermite插值'})
 
figure(2);
 
for leinuoshu=1:2
 
    if leinuoshu == 1
        data = textread('CJ0130.jin');  
    else
        data = textread('CJ0140.jin');
    end
    num_line = 1+9+96+2*135+187+2;
 
    for num=1:41
 
        angle(num) = data( (num-1).*num_line+1 ,2);
 
 
        PSI=6894.76;
        p = data( (num-1).*num_line+1 ,6)*9.8;
 
        for i=1:135
            up_wall_pressure(i) = data((num-1).*num_line+1+96+187 +i,4 )*PSI; 
            low_wall_pressure(i) = data((num-1).*num_line+1+96+187+135 +i,4 )*PSI;
        end
 
        for i = 1:96
            wing_pressure(i)= data((num-1).*num_line+1 +i,4 )*PSI;
        end
 
        sum_up_wall_pressure(num) =trapequad(up_wall,up_wall_pressure);
        sum_low_wall_pressure(num) =trapequad(low_wall,low_wall_pressure);
        sum_wall_pressure_difference(num) = sum_low_wall_pressure(num) - sum_up_wall_pressure(num);
 
        sum_wing_pressure(num) =-trapequad(wing_x,wing_pressure);
 
        sum_wall_pressure_difference(num) = sum_wall_pressure_difference(num)/p/s;
        sum_wing_pressure(num) = sum_wing_pressure(num)/p/s;
    end
 
    n = 1:1:41;
    plot(angle,sum_wall_pressure_difference,'-o',angle,sum_wing_pressure,'-*');
    title('不同迎角下DU 93-210翼型的升力系数');
    if leinuoshu == 1
        legend('洞壁压力差','升力','Location','best');
    else 
        legend('洞壁压力差','升力','Location','best');    
    end
    hold on
end

​ 我们得到不同迎角下的升力系数如图所示，其中上图为雷诺数为2.0×106时的结果，下图为雷诺数为3.0×106时的结果：

​ 可以看出，随着迎角的增大，升力系数增大。雷诺数对翼型升力系数的影响较小，我掌握的数据较少，总的来看，雷诺数小时，升力系数的变化曲线较为平缓。升力系数的绝对值在0-1.5之间，符合流体力学常识。我们看到，DU 93-210 翼型在迎角为10°时升力系数达到了1.4以上，可以说是一款非常优秀的翼型，达到了风力机的翼型所能达到的较好水平。

​ 当然，翼型的升力系数并不能完全衡量飞机的升力，坊间传言歼-20的升力系数高达2以上，那是全机系数。在増升装置的作用下，高速飞行的飞机升力系数最大甚至超过4。

翼型实验的重要性

​ 众所周知，为了保持民用飞机的安全性和经济性，飞行马赫数不宜超过临界马赫数。想要提高飞行速度就要设法提高机翼临界马赫数。减小机翼厚度或采用后掠机翼可以提高临界马赫数，但是这样会增加机翼重量。采用超临界机翼则可提高临界马赫数，同时不必付出增加机翼重量的代价。

​ 国产大型客机C919采用的正是超临界翼型，它是国内第一次完全自主设计的超临界机翼，采用“一体化，弱激波”的设计理念，代表我国民机技术的时代性和先进性。采用这种翼型设计，使得C919在跨音速范围内的气动性能大幅增强，提高机翼的临界马赫数，使高亚音速时阻力急剧增大的现象推迟发生，并提高姿态可控性。由于采用超临界机翼，C919的飞机结构重量也相应减少，机翼内的容积增加，可以放置更多的燃油或其他设备，获得更好的经济性。

​ 经过无数工程师的勤奋探索，已经发展出了许多性能优异且使用成熟的翼型，但遗憾的是，至今为止仍然没有找到一款可以适合所有飞机的超级翼型。究其原因，机翼设计实际上也反映了飞机设计的系统性和复杂性，任何一款翼型，如果拥有较高的升力系数，其阻力也会随之增大，对发动机的要求就会提高，而发动机性能的提升又会带来重量的增加等一系列问题；如果超音速性能优异，则其低速性能特别是起降过程中的气动效率就会有损失。

​ 为此，飞机设计过程中必须确定该型飞机的最主要飞行状态，设计师根据其飞行状态再具体为其选择最合适的翼型。换言之，飞机设计并不是单独追求某个性能的最优，而是在反复权衡比较之后的整体性能最优，是一项关联复杂、逻辑严密的系统工程。

​ 因此，翼型实验的意义在于，通过实验手段测量可用的翼型无量纲化升力系数、阻力系数，而飞机总设计师根据这些系数结合飞机整体需求选择最合适的翼型。由此可见，通过数值计算方法提高数值积分精度，进一步提高翼型系数的测量精度，在工程中是非常有意义的。

参考

1. 大飞机报《机翼那些事儿》http://www.comac.cc/mjkp/xzs/201612/21/t20161221_4665263.shtml?from=timeline
2. 数值分析_维基百科 https://zh.wikipedia.org/zh/数值分析

备注

​ 2019年夏学期上课初期较为痛苦，当时还在打中控杯，第一次作业很仓促地用c++应付了(老师没有限制语言和平台)，后来的本想用matplotlib和NumPy两个Python库来完成，然而助教一脸你们应该用Matlab对叭我要统一跑程序的表情，转而学了Matlab这样一个不在我计划内的工具。数值计算方法内容广泛，远非我们课堂所能覆盖，这门课程又莫得期末考，我就抱着划水的心态完成了若干个数十页报告，以至于到现在对那些算法也没有什么深入的理解，想要真正掌握，日后还得还债。

​ 数值计算方法这门课启发了我的思考，一来是引发的我对混沌的好奇，在这里不介绍这个概念，众所周知的蝴蝶效应就是对非线性系统的形象描述。我们认为，在线性系统中，通过数值计算方法提高数值的精确度，就可以提升运算结果的精确度，而生活中存在着大量非线性系统，最早意识到它的重要性的是一位气象学家，他反复向同一个天气预测系统中输入同一组数据，得到了迥乎不同的结果，其诱因就是输入存在着微小的误差。

​ 二来是让我感受到了迭代的魅力。所谓迭代法，在计算数学中，迭代是通过从一个初始估计出发寻找一系列近似解来解决问题（一般是解方程或者方程组）的数学过程，为实现这一过程所使用的方法统称。迭代法用于求解非线性方程组，令我觉得稀松平常。但是迭代法用于求解大型线性代数方程组，却令我觉得非常迷人。迭代法有两个重要的点，一是如何量化误差；二是如何从当前点向精确解迭代，第二点又包含两个细节，一是搜索方向，二是确定步长。在求解非线性方程时，我们的关注点都放在迭代过程上，因为误差易于量化，取差的绝对值即可。而对于线性代数方程组，我们引入了范数来量化误差。在人工智能的神经网络中也有类似的问题，去通过学习确定各层神经网络的参数，也涉及量化误差和迭代这两步。常见的损失函数有均方误差和交叉熵误差，也常常通过梯度来迭代，求得损失函数的极值。在大数据时代，根据公式求得精确解的方法愈发受到限制，在更多的问题当中，我们都是通过迭代的思想逐步逼近真实解，求得一个满足精度要求的近似解。

• ← Previous Post
• Next Post →