算例收集3：NACA0012机翼绕流
1.问题描述
通过计算流体力学（Computational-Fluid-Dynamics,CFD）模拟机翼绕流是一个经典的问题，这也是计算流体力学这一学科起源的原因之一。
通过CFD模拟，可以快速评估不同翼型设计的升阻比，识别最优几何形状或攻角，从而在减少阻力的同时最大化升力，显著提升飞行器的整体性能；同时，CFD可在虚拟环境中模拟复杂流动，快速迭代设计方案，减少物理实验次数；而且，CFD还可以输出流场细节（如压力分布、边界层分离、涡流结构等），这些信息帮助工程师理解流动现象的物理机制，例如失速原因或激波位置，从而针对性优化翼型设计。
NACA0012是一种经典的对称翼型，属于NACA四位数字系列。其名称中的“00”表示无弯度（中弧线为直线），“12”代表最大厚度为弦长的12%。
1、对称结构
上下表面完全对称，零攻角时不产生升力，适用于双向流动或需要对称气动性能的场景，如直升机旋翼、螺旋桨和飞机尾翼。
2、气动特性
升力随攻角线性增加，失速前具有稳定的气动性能；
阻力特性在中低攻角下较为平缓，适合需要低阻设计的部件；
常用于验证CFD模拟和风洞实验，作为基准案例测试湍流模型与计算精度。
1.1控制方程
目前，实际流动问题的解多是通过求解N-S方程来获得的，因为N-S方程能够反映流体的质量守恒、动量守恒和能量守恒的规律。而由于本文针对的是机翼升阻特性，可以不考虑粘性影响，即二维无粘可压缩流，因此本文研究的控制方程是Euler方程。
Euler方程可以改写为：
接回上文，W是守恒变量：
1.2计算网格
本文所选择的进行算例分析的是NACA0012翼型，采用以三角形为基本单元的非结构网格作为计算网格。网格总节点数为5306，总边数为15688,总单元数为10382。
fid=fopen(('naca0012.grd'),'r'); % READ GRIDS
nnodes=fscanf(fid,'%d',1); %!nnodes : node numbers
fseek(fid,4,1);
nedges=fscanf(fid,'%d',1); %!nedges : edge numbers
fseek(fid,4,1);
ncells=fscanf(fid,'%d',1); %!ncells : cell numbers
fscanf(fid,'\n');
xy(2,nnodes)=0;
iedge(4,nedges)=0;
icell(3,ncells)=0;
vol(ncells)=0;
cell_v(4,ncells)=0;
node_v(nnodes)=0;
node_P(nnodes)=0;
for i=1:nnodes
xy(1,i)=fscanf(fid,'%f',1); %!xy : cartesian coordinates xy(1,i) for x xy(2,i) for y
xy(2,i)=fscanf(fid,'%f',1);
end
for i=1:nedges %!iedge : matrix for edge
iedge(1,i)=fscanf(fid,'%f',1); %!iedge(1,i) : edge's start point 'a'
iedge(2,i)=fscanf(fid,'%f',1); %!iedge(2,i) : edge's end point 'b'
iedge(3,i)=fscanf(fid,'%f',1); %!iedge(3,i) : edge's left cell (all positive)
iedge(4,i)=fscanf(fid,'%f',1); %!iedge(4,i) : positive for right cell
% !iedge(4,i) -1 for wall boundary
% !iedge(4,i) -2 for farfiled boundary
end
for i=1:ncells %!icell : triangle cell's three nodes index
icell(1,i)=fscanf(fid,'%f',1);
icell(2,i)=fscanf(fid,'%f',1);
icell(3,i)=fscanf(fid,'%f',1);
end
for i=1:ncells %!vol : cell's volume(area in 2d)
vol(i)=fscanf(fid,'%f',1);
end
fclose(fid);
2.数值方法
2.1中心格式 + 人工粘性(jamson格式)
将计算域划分为有限个互不重叠的单元，并将积分守恒方程应用于每个单元，由于各个单元的面积不随时间变化，所以可将控制方程改写为：
U = 0.5 * (cell_v(1,L) + cell_v(1,r));
V = 0.5 * (cell_v(2,L) + cell_v(2,r));
pa = 0.5 * (cell_v(3,L) + cell_v(3,r));
rho = 0.5 * (cell_v(4,L) + cell_v(4,r));
H = 0.5 * (cell_v(6,L) + cell_v(6,r));
alpha(i)=abs(U*dy-V*dx)+sqrt(gama*pa/rho)*ds;
Z=U*dy-V*dx;
flux(1) = Z*rho;
flux(2)= Z*rho*U+dy*pa;
flux(3) = Z*rho*V-dx*pa;
flux(4)= Z*rho*H;
for j=1:4
cell_s(j,L)=cell_s(j,L)+flux(j);
cell_s(j,r)=cell_s(j,r)-flux(j);
end
如上所述的中心格式是不包含耗散项的，所以任何误差（离散误差，循环误差等等）都不会衰减，最后的定常解可能会出现振荡。为了减少这些振荡，在方程的右边加了人工耗散项，对于k单元，方程变为：
由此确定的自适应系数变为：
function [cell_D]=compute_artificialdissipation(nedges,ncells,iedge,cell_w,cell_v,k_2,k_4,alpha)
d2_w(4,ncells)=0;
cell_D(4,ncells)=0;
d2(4)=0;
d4(4)=0;
for i=1:nedges
L=iedge(3,i);
r=iedge(4,i);
if(r~=-1&&r~=-2)
for j=1:4
average=0.5*(cell_w(j,L)+cell_w(j,r));
d2_w(j,L)=d2_w(j,L)+average-cell_w(j,L);
d2_w(j,r)=d2_w(j,r)+average-cell_w(j,r);
end
end
end
for i=1:nedges
L=iedge(3,i);
r=iedge(4,i);
if(r~=-1&&r~=-2)
V=abs(cell_v(3,r)-cell_v(3,L))/(abs(cell_v(3,r)+cell_v(3,L)));
e2=k_2*V;
e4=max(0,k_4-e2);
for j=1:4
d2(j)=alpha(i)*e2*(cell_w(j,r)-cell_w(j,L));
d4(j)=-alpha(i)*e4*(d2_w(j,r)-d2_w(j,L));
cell_D(j,L)=cell_D(j,L)+(d2(j)+d4(j));
cell_D(j,r)=cell_D(j,r)-(d2(j)+d4(j));
end
end
end
2.2边界条件和流场初始条件
2.2.1---边界条件
物面边界条件
对于无粘流体，在物面边界上，应该满足流动方向与物面相切，即物面的法向速度分量为零，可知
这就是物面无穿透边界条件。因此，X,Y方向动量的通量不是零，而与压强有关。
远场边界条件
用有限体积法进行空间离散时，只能取一个有限远的边界作为远场边界，而在绕物体的实际流动中并不存在这个边界，物体所产生的扰动可以认为被传到无穷远而没有反射，即无反射边界条件。这里，采用一维Riemann不变量来处理远场边界条件。
可以将远场边界条件分为如下四种情况：
（1）亚音速入流条件
它有三条入流特征线，需规定三个条件：
其中下标t表示切向，n代表法向, \(()_{\infty}\)代表自由流，\(()_{e}\)代表从流场内到边界的外插值，上式的右端皆为已知，可解出边界上的qn、c、qt、s的值，再由s和c求出p和 \(\rho\)。
(2)亚音速出流条件
它有一条入流特征线，需规定一个条件:
（3）超音速出流条件
无入流特征线，不需在边界上规定边界条件:
其中的W表示边上的守恒变量，\(W_{l}\)表示与此边相邻元素的守恒变量值。
（4）超音速入流条件
它有四条入流特征线，需规定全部四个条件:
其中的W表示边上的守恒变量，\(W_{\infty}\)表示来流值。
qn(1)=Ma * clocal * cos(Theta)*y_s+Ma * clocal * sin(Theta)*x_s; %qn_远场
qn(2)=cell_v(1,L)*y_s+cell_v(2,L)*x_s; %qn_内场
qt(1)=-Ma * clocal * cos(Theta)*x_s+Ma * clocal * sin(Theta)*y_s; %qt_远场
qt(2)=-cell_v(1,L)*x_s+cell_v(2,L)*y_s; %qt_内场
c(1)=clocal; %远场声速
c(2)=sqrt(gama*cell_v(3,L)/cell_v(4,L)); %内场声速
qn(3)=0.5*(qn(1)+qn(2)+2*(c(2)-c(1))/(gama-1)); %qn_边界（公式计算）
c(3)=(gama-1)/4*(qn(2)-qn(1)+2*(c(1)+c(2))/(gama-1)); %声速_边界（公式计算）
if(Ma<=1) %亚声速远场条件
if(qn(3)<=0) %入流条件
qt(3)=qt(1);
rho=(c(3)/c(1))^(2/(gama-1))*flow(4);
pa=rho*(c(3)^2)/gama;
else %出流条件
qt(3)=qt(2);
rho=(c(3)/c(2))^(2/(gama-1))*cell_v(4,L);
pa=rho*(c(3)^2)/gama;
end
else %超声速远场条件
if(qn(1)<=0) %入流条件
qn(3)=qn(1);
qt(3)=qt(1);
rho=flow(4);
pa=flow(3);
else %出流条件
qn(3)=qn(2);
qt(3)=qt(2);
rho=cell_v(4,L);
pa=cell_v(3,L);
end
end
U=qn(3)*y_s-qt(3)*x_s;
V=qn(3)*x_s+qt(3)*y_s;
H=pa/(rho*(gama-1))+(U^2+V^2)/2+pa/rho;
alpha(i)=abs(U*dy-V*dx)+sqrt(gama*pa/rho)*ds;
Z=U*dy-V*dx;
flux(1) = Z*rho;
flux(2)= Z*rho*U + dy*pa;
flux(3) = Z*rho*V-dx*pa;
flux(4)= Z*rho*H;
for j=1:4
cell_s(j,L)=cell_s(j,L)+flux(j);
end
2.3时间推进
定常解是由对常微分方程进行时间积分获得的，此方程可以写成：
方程的右侧的值表示每个单元k中心的残值：
采用显式的四步格式来完成，由于对于定常解时间的精确并不重要，所以选择此格式只是因为它的稳定和衰减的特性。目前采用的是以下的格式：
式中n是当时的时间步数，n+1是新的时间步数：
为了减少计算时间，只在第一步时计算耗散项D，然后在接下来的各步中D是一个常值。对于一个标量模型方程，这种做法改变了格式的稳定区，但是精度和收敛特性被保留了下来。以上格式在应用于非定常Euler方程时，能保持稳定的CFL数最大可以取到2*sqrt(2)
这个显式格式最大的缺点是最大的时间步长受到限制，因为稳定域受到限制。而且，对于多维的方程组，最大时间步长只可以用近似的方法获得。对于固定形状的网格，采用以下的表达式：
下面通过时间推进直接展示了整个计算程序的框架：
1、先通过单元守恒变量计算界面通量
2、在4步Runge-Kutta中的第一步计算人工粘性
3、计算残值，更新单元守恒变量，最后直至收敛
4、根据需要通过守恒变量计算出独立的变量
a=[1/4,1/3,1/2,1];
for T=1:15000
for j=1:4
for i=1:ncells
cell_wR(j,i)=cell_w(j,i);
end
end
for m=1:4 %4步Runge-Kutta
[alpha,cell_s]=compute_flux(nedges,ncells,iedge,xy,flow,cell_v,Ma,gama,clocal,Theta);
if(m==1)
cell_D=compute_artificialdissipation(nedges,ncells,iedge,cell_w,cell_v,k_2,k_4,alpha);
cell_dt=compute_timestep(nedges,ncells,iedge,alpha,vol,CFL);
end
for i=1:ncells
for l=1:4
Rk=-(cell_s(l,i)-cell_D(l,i))/vol(i);
cell_w(l,i)=cell_wR(l,i)+a(m)*cell_dt(i)*Rk;
end
end
for i=1:ncells
cell_v(4,i)=cell_w(1,i); %rho
cell_v(1,i)=cell_w(2,i)/cell_w(1,i); %U
cell_v(2,i)=cell_w(3,i)/cell_w(1,i); %V
cell_v(5,i)=cell_w(4,i)/cell_w(1,i); %E
cell_v(3,i)=(gama-1)*cell_v(4,i)*(cell_v(5,i)-0.5*(cell_v(1,i)^2+cell_v(2,i)^2));%Pa
cell_v(6,i)=cell_v(5,i)+cell_v(3,i)/cell_v(4,i); %H
end
end
end
3.计算结果
计算状态如表所示：
浙公网安备 33010602011771号