高等流体力学作业

高等流体力学试题
姓名：岳万凤
学号：20091002061
学院：动力工程学院
专业：动力工程及工程热物理指导老师：何川教授
完成日期：2010.1.10
四、设52/ 1.110/m s μρ-=⨯的气体以10/V m s ∞=的速度以零攻角绕流长度为L ＝1m 的大平板，试用数值解讨论边界层内的流动规律。

解：基本思路：根据绕流平板的大雷诺数边界流动控制方程，运用数量级比较法对控制方程进行简化，然后引入无量纲相似变量将动量方程化为常微分方程。

最后运用龙格库塔方法编写C 程序对该常微分方程进行数值求解，获得沿平板长度方向不同位置处边界层内气体速度在平板法线方向上的变化规律以及边界层厚度沿平板长度方向上的变化规律。

由题意，气体外掠大平板流动的雷诺数：
90909010
1.10.110/Re 5
=⨯⨯=⋅=
-∞ρμL V ＞3
10 1）取平板的左端点O 为坐标原点，x 轴沿着平板，y 轴垂直于平板，建立流向坐标系Oxy ，如图1所示。

边界层控制微分方程组及相应定解条件为：
对二维定常流动，不计重力，其边界层基本方程：
C.F.：
0u v x y
∂∂+=∂∂ M.F.：2222
1()u u p u u
u v
x y x x y μρρ∂∂∂∂∂+=-++∂∂∂∂∂ 22221()v v p v v u v x y y x y
μρρ∂∂∂∂∂+=-++∂∂∂∂∂ 边界条件： 0,0,0
y u v === ,10/y u V m s ∞
→∞== 这里：自变量(0,)x L ∈，(0,)y δ∈；
因变量(0,)u V ∞∈；max (0,)v v ∈。

图1 气体外掠大平板示意图
作变换：*(0,1)x x L =
∈，*(0,1)y y y AL δ==∈，其中A 待定，且A 远小于1 *(0,1)u u V ∞
=
∈，*max (0,1)v v
v v BV ∞=
=∈，其中B 待定,且B 远小于1； *2
(0,1)p
p V ρ∞
=
∈。

代入方程组：
C.F: 0][**
*=∂∂+∂∂*∞y
v A B x u L V 1
1B
A
⋅→A 和B 数量级相当 M.F: ]1//[)(*
*2**222**22
*****
*2
x p y u A L V x u L V L V y u v A B x u u L V ∂∂-∂∂+∂∂=∂∂+∂∂∞∞∞∞ρμρμ 1
1B A ⋅ 11Re ⋅ 211
1Re A
⋅⋅ 1 ]11//[)(**2**222**22
*****
*2
y
p AB y v A L V x v L V L BV y v v A B x v u L BV ∂∂-∂∂+∂∂=∂∂+∂∂∞∞∞∞ρμρμ 1
1B A ⋅ 11Re ⋅ 2111Re A ⋅⋅
1
AB
因为Re 远大于1，所以带1
1Re
⋅的项可忽略，数量级比较后，原控制方程组可
简化为：
C.F:
0u v x y
∂∂+=∂∂ M.F : x
p y u y u v x u u ∂∂-∂∂=∂∂+∂∂ρρμ122 0=∂∂y
p
从上式可以看出，在任一过流断面上，边界层内各个点的压力p 与边界层边界上的压力e p 相等，在边界层外由伯努利方程：
22
1122
e e p V p V ρρ∞∞
+=+，其中2222e e e e u v u V ≈+= 则：2
211111()22
e e p p p V V x x x ρρρρρ∞∞∂∂∂-
=-=-++∂∂∂
= e e e
e dV du
V u dx dx
≈ 对于零攻角绕流平板流动：
(),0e e e du
u x u dx
==
这样控制方程组最终简化为：
C.E ：
0u v x y
∂∂+=∂∂ M.E: 2
2y
u
y u v x u u ∂∂=∂∂+∂∂ρμ B.C:
∞
==∞→===V u u y v u y e :0,0:0
2）引入流函数(,)x y ψ，并作相似性变换：
令,u v y x
ψψ∂∂=
=-∂∂，原方程可表述为以,x y 为自变量，ψ为因变量的方程。

令无因次变量：
()
η=
=
=
y y
y
g x ； 即
=()g x ；
则()
()ψ
η∞=
=⋅f f V g x ；
故而可得：
()()
()()
1
f V
g x u V g x f V f y y g x ψηη∞∞∞∂⋅⋅∂∂''
====∂∂∂()22V u f V f V V f f x x x x x η
ηηη∞∞∞∞
'⎛⎫∂∂∂∂'''''===-=- ⎪∂∂∂∂⎝⎭
； ()()()
1
V f V u f V V f f y y y g x g x ηη∞∞∞∞'∂'∂∂∂''''====∂∂∂∂； ()()()()()
2221
V V V V u f f f f y y g x g x g x g x y g x ηη∞∞∞∞⎛⎫''∂∂∂∂''''''''==== ⎪ ⎪∂∂∂∂⎝⎭
()()
()()
()()()()22V g x f f v V g x f g x x x
x g x V g x f g x f V f fg x x x ψ
ηηηη∞∞∞∞
∂⎡⎤
∂∂∂'=-=-
=-+⎢⎥∂∂∂∂⎣⎦
⎡⎤⎡⎤
⎛⎫''''=-+-=-⎢⎥⎢⎥ ⎪⎝⎭⎢⎥⎣
⎦⎣⎦
代入方程，有：
()()()()
()2
22g x V V V V f f V f fg x f f x x g x g x μηηρ∞∞
∞∞∞⎡⎤⎛⎫
''''''''''-+-=⎢⎥ ⎪⎢⎥⎝⎭⎣⎦； 又 (
)()22g x g x x x ''===
； 则有：
()
()2222
2222g x V V V V f f f f f f f x x x x g x μηημρρ∞∞∞∞''''''''''''-+-=； 整理为：
f f f f ff ff f ηη'''''''''''''-+-=-=；
即：0f ff '''''+=。

边界条件：
η∞''=⇒====00,0;故(0)0y u V f f ；
ηη∞'=⇒==-+=='
00,[()'()]0;故(0)0x y v V f g x g x f f ；
()()η∞∞''→∞⇒→∞=∞=∞=,;故1y u V f V f 。

因此，原定解问题可表示如下：
()
()()0
00,00,1f ff f f f '''''⎧+=⎪⎨
''==∞=⎪⎩ 将上述定解问题中的高阶常微分方程表示为一阶常微分方程组：
η⎧'==⎪⎪⎪
'''==⎨⎪
''''''''===-=-⎪⎪⎩
112
122df
f f d f f f f f f ff ff
即有：1
1222
f f f f f ff '⎧=⎪⎪
'=⎨⎪
'=-⎪⎩因此可令：⎧=⎪⎪
⎪=⎨⎪⎪=-⎪⎩dx
y dt dy z dt dz xz dt
相应边界条件为：()()()00,00,1x y y ==∞=。

3) 将上述常微分方程的边值问题变换为初值问题
令：()()()1133
0;;f A A f A F ξηηξ''=== 则有：
()()12
33
22334
3f A F A F f A F A F AF f AF A F ξξηξξξηηξξη⎧∂'''==⎪∂⎪
⎪⎛⎫∂∂∂⎪'''''''=== ⎪⎨
∂∂∂⎝⎭⎪
⎪∂∂⎪''''''''==⎪∂∂⎩
代入原方程有：
()
414333
0A F A F AF A
F FF ''''''''''+⋅=+=
由于A 为非零常数，则有：0F FF '''''+=。

初始条件：
1
3
232332(0)(0)0(0)0
'(0)'(0)0'(0)0
''(0)''(0)''(0)1'()'()1['()]f A F F f A F F f AF A F f A F A F -⎧==⇒=⎪⎪
⎪==⇒=⎨
⎪==⇒=⎪⎪∞=∞=⇒=∞⎩
即转化的初值问题为：
()()()0
00,00,01
F FF F F F '''''⎧+=⎪⎨
'''===⎪⎩ 这样原问题转化为：
()
()()0
00,00,0f ff f f f A '''''⎧+=⎪⎨
'''===⎪⎩ 为求得A 值，须计算()F '∞．
4)、将该常微分方程表示为一阶常微分方程组：
令：ξ⎧'==⎪⎪⎪
'''==⎨⎪
''''''''===-=-⎪⎪⎩112122dF
F F d F F F F F F FF FF
即有：11222
F F F F F FF '⎧=⎪⎪
'=⎨⎪'=-⎪⎩
相应初值为：()()()'''===00,00,01F F F
利用四阶龙格库塔法求解此微分方程初值问题，令F x =，1F y =，2F z =，
则有：ξξξ
⎧=⎪⎪⎪=⎨⎪⎪=-⎪⎩dx
y d dy z d dz
xz d
初始条件：()()()00,00,01x y z ===。

5）分别计算x =0.01m 、0.1m 、0.5m 和1.0m 处，边界层内速度随平板法向距离变化规律的变化曲线。

源程序见附录。

水平方向速度u 随平板法向距离y 的变化曲线如图2所示；法向速度v 随平板法向距离y 的变化曲线如图3所示。

从图2可以看出，在相同的y 处，若该位置离平板前端越远，其水平方向速度越小。

而从图3可以看出，在离平板起始端越远的位置，对应的边界层内方向速度越小，在离平板起始端附近，边界层内法向速度较大。

6）计算边界层厚度沿流动方向的变化曲线如图4所示。

若要将边界厚度的值拟合成无因次的函数关系，需要这样处理：
令：X x
δ
=
，Y =
=可见X 和Y 均是无量纲的量。

通过程序
计算得到的数据计算X 和Y 的值，再拟合就可得到边界层厚度的无因次关系式， 计算结果如附表1所示。

拟合得到的无因次关系式：
Y X 0503.9=
表1 边界层内速度u、v沿平板法向距离的计算值
图2 平板不同位置处边界层内水平方向速度沿法向距离的变化曲线
图3 平板不同位置处边界层内法向速度沿法向距离的变化曲线
即：x
x
Re 0503.9=
δ
，无因次曲线关系曲线如图5所示。

图4 边界层厚度沿流动方向的变化曲线
δ（m ）
图5 边界层厚度无量纲关系式曲线
附源程序：
#include<stdio.h>
#include<math.h>
main()
{double k1,k2,k3,k4,l1,l2,l3,l4,m1,m2,m3,m4;
double a,b=0,u,v,m,n,y0,y1,x=0,y=0,z=1,h=0.1,h1=0; //定义变量和步长，赋初值// do{ y0=y;
k1=y; l1=z; m1=-x*z;
k2=y+l1*h/2; l2=z+m1*h/2; m2=-(x+k1*h/2)*(z+m1*h/2);
k3=y+l2*h/2; l3=z+m2*h/2; m3=-(x+k2*h/2)*(z+m2*h/2);
k4=y+l3*h; l4=z+m3*h; m4=-(x+k3*h)*(z+m3*h);
x=x+h*(k1+2*k2+2*k3+k4)/6;
y=y+h*(l1+2*l2+2*l3+l4)/6;
z=z+h*(m1+2*m2+2*m3+m4)/6; //四阶龙格库塔法求'()
F //
} while(fabs((y0-y)/y)>0.0000001);
printf("%f\n", y);
a=pow(y,-1.5);
printf("%f\n", a); //求A='''(0)
f//
x=0; y=0; z=a;
m=0.000011*10/2/0.01;
n=2*0.000011*0.01/10; //对x=0.1m、0.5m和1m处，直接将0.01替换即可// do { y0=y;
u=10*y;
v= (h1*y-x)*pow (m, 0.5);
y1=h1*pow(n,0.5);
printf ("%f %f %f %f\n", h1, y1, u, v);
l1=z; k1=y; m1=-x*z;
k2=y+l1*h/2; l2=z+m1*h/2; m2=-(x+k1*h/2)*(z+m1*h/2);
k3=y+l2*h/2; l3=z+m2*h/2; m3=-(x+k2*h/2)*(z+m2*h/2);
k4=y+l3*h; l4=z+m3*h; m4=-(x+k3*h)*(z+m3*h);
x=x+h*(k1+2*k2+2*k3+k4)/6;
z=z+h*(m1+2*m2+2*m3+m4)/6;
y=y+h*(l1+2*l2+2*l3+l4)/6; //四阶龙格库塔法求x=0.01m处的u, v//
h1=h1+h;
} while (fabs ((y-y0)/y)>0.0000001);
x=0;
do { n=2*0.000011*x/10; //求边界层厚度// h=h1*pow (n, 0.5);
printf ("%f %f\n", x, h);
x=x+0.01;
} while(x<=1.0);
}
附表1
五、谈谈自己所关心的流动问题，并给出自己对该问题的看法。

希门茨流动问题求解
——对非线性N-S 方程精确解的认识
希门茨流动是在平面滞止区域的一种流动，通常称平面驻点流动。

绕流问题总存在流动的滞流区域。

该流动实际过程可以简化为：一平壁面无限大，距壁面无穷远处流动速度垂直壁面，在流向壁面的过程中，有一条始终垂直壁面的流线，该流线上的速度逐渐减小至达到壁面后速度降为零，该流线与壁面的交点即为驻点，滞止流线两侧的流体分别向左右分流，如图1所示。

图1 希门茨流动
1、数学描述
由不可压缩流动的连续方程
0=∂∂+∂∂y
v
x u (1) 定义流函数),(y x ψ
y u ∂∂=
ψ , x
v ∂∂-=ψ
(2) 设无粘性的流函数axy y x id =),(ψ 在如图所示的平面有势流动中有
⎩
⎨
⎧-==ay V ax
U (3) 于是由伯努利方程，可得压强分布为
22201
p p a (x y )2
ρ=-+ (4)
P 0 为驻点压强。

（4）式中的流速与压强均满足势流方程，也满足不可压缩粘性流动方程。

而对于势流有ϕ∇=u ，所以0)(22=∇∇=∇ϕu 。

这就是说N-S 方程中的粘性项对于势流而言恒等于零。

所以，尽管它们也满足不可压缩粘性流动的运动
方程，但却不满足粘性流动的无滑移边界条件。

为此，对适用于粘性流体的流函数做出如下假设)(),(y axf y x =ψ，则：
u xf '(y)v f (y)=⎧⎨
=-⎩ (3)
22
01p p a [x F(y)]2ρ-=
+ (4)
并代入恒定的N —S 方程,并简化为
22f 'ff ''f '''a ν-=+ (5)
2
1ff 'F'f ''2a ν=
- (6)
原边界条件：
⎪⎩
⎪
⎨⎧=====∞====0p
p 0;y 0,x ;y 0v ,0;0y ax U u u (7) 边界条件化简为：
⎪⎩
⎪
⎨⎧====∞====0F 0;y 0,x 1
'f ;y 0
f ,0'f ;0y (8) 对方程（8）作变量置换，令：
)(A )y (f ,y ηϕαη== (9)
代入（8）式得：
22223A (''')a A '''αϕϕϕναϕ-=+ (10)
如果使3222A a A ανα==，则使方程简化。

由此，（9）式可改写为
)
(a )y (f ,a
y
ηϕνν
η== (11)
于是方程（5）可简化为
2''''''10ϕϕϕϕ+-+= (12)
其边界条件变为：
1';0',0;0=∞====ϕηϕϕη (13)
(12)式为一个非线性方程，要得到它的其解析解仍十分困难,我们采用数值方法求解。

2、微分方程衍化
微分方程的数值解，首先由希门茨给出，之后由霍华斯加以改进。

在此，笔者采用解常微分方程的常用数值解法—龙格-库塔(Runge-Kutta)法求解的方法，对边界层进行试探性求解。

在无限远处的边界条件对编写程序很不方便的,因为在程序中每一个规定或计算的值必须是有限的。

为此，为了避免出现计算1',=∞=ϕη的编程困难，进一步作变量置换，令：
β
φηϕβη)
g ()(,g =
= (14) 引入的β为待定常数，在此处起调节作用。

于是有：
n n n
n+1n d 1d d dg ϕφ
ηβ= (15)
这时方程和边界条件变为：
242
''''''0g 0;(0)0,'(0)0g ;'()φφφφβφφφβ⎫
+-+=⎪
===⎬
⎪
=∞∞=⎭ (16)
由于)0('')0(''3ϕβφ=，因待定常数β可以起到调节作用，所以，尽管)0(''ϕ未知，我们仍可假定一个''(0)φ，以起到假定)0(''ϕ的目的。

不妨假定1)0(''=φ，看)('∞φ是否趋向1，如果不趋近1则改变)0(''φ，最后试算出 1.232533)0(''=φ。

这样，方程（16）就可为化为：
⎪
⎭
⎪
⎬
⎫
======+-+232533.1)0(;00)0(,0)0(;00'''42''''''φφφβφφφφg g (17)
此时，对(17)编程求解就可以避免边界条件中对无穷远处定解条件的处理。

为了方便计算和编程，取调节常数为β=1，方程（17）可化为下面的形式
⎪
⎭
⎪
⎬
⎫
======+-+232533.1)0(;00)0(',0)0(;001''2''''''φφφφφφφg g (18)
3、微分方程数值解
3.1用MATLAB 求常微分方程的数值解
针对方程的特征，采用以Runge-Kutta 为基础的常微分方程初值问题数值解的求解函数—ode 函数进行求解。

方程（18）可写为如下形式：
⎪
⎭
⎪
⎬⎫
======+-+232533.1;00,0;001'''2''''''y x y y x y yy y (19)
下面用MATLAB 的ode 函数求方程(19)的数值解。

⑴、将方程化为ode 函数的标准方程
MATLAB 提供的函数命令只适用于求解一阶的常微分方程的初值问题，对于高阶微分方程，必须先用替换法化为形如y'=f(x,y)的一阶微分方程组后，才能用ode 函数进行求解。

将微分方程的导数降阶，可令：123,',''y y y y y y ===，则原方程变为：
12232
3132'''1y y y y y y y y ⎧=⎪
=⎨⎪=-+-⎩ (20) 其初值条件为：
⎪⎩⎪
⎨⎧===232533.1)0(0
)0(0)0(3
21y y y (21)
⑵、建立函数文件odex.m function ydot=odex(x,y)
ydot=[y(2);y(3);y(2)^2-1-y(1)*y(3)]
⑶、解微分方程:在MATLAB 命令窗口中运行函数[Beta,phy]=ode45(@fun, [0:0.2:4.6],[0;0;1.232533])。

⑷、画图观察其变化趋势 在MATLAB 命令窗口中运行函数
plot(Beta,phy(:,1),'+r',Beta,phy(:,2),'*g',Beta,phy(:,3),'-b') 得出数值，绘成图象如下。

3.2 利用excel 求解常微分方程的数值解
在excel 界面里，Excel 对单元格数值的默认精度是15位（包括小数点在内），会进行“四舍五入”。

所以计算数值会出现误差。

故在该计算过程中，我们选取较小步长为x ∆=0.01进行计算与求解，以尽量减少误差。

利用欧拉公式，
x y y y n n n ∆+=--'11，x y y y n n n ∆+=--''1'1'，x y y y n n n ∆+=--'''1''1''。

计算过程：由已知232533.1''0=y ，012''''''=+-+y yy y 可知，'
''0y =-1 由欧拉公式，x y y y ∆+='
01 ⇒ 001.0001=⨯+=y x y y y ∆+=''0'0'1 ⇒ 1232533.001.0232533.10'
1=⨯+=y x y y y ∆+='''0''0''1 ⇒ 222533.101.01232533.1'
'1=⨯-=y 根据原方程，可知
99985.01'
'112
'2'
''1-=--=y y y y
然后，以x ∆=0.01为步长，在excel 界面里往下拉动即可得到其他解的值。

将计算的数据绘成的图像如下。

4、结果对比与分析
将上述两种方法计算的数值与霍华斯求解的数值进行对比。

分析对比图如下。

采用以上方法进行的计算结果与前人霍华斯求解的结果能够很好的符合，由于目前笔者对计算流体力学的学习不够深入，对计算流体力学的思想认识很肤浅，在求解过程中遇到了很多困难，比如在边值问题转化为初值问题时存在编程上的问题，所以在今后的学习过程中还要更深入的学习数值计算方面的内容。

5、对非线性N-S方程精确解求解的认识
经过对希门茨流动问题的求解，笔者对流动问题的非线性N-S方程的精确解有了更深的认识，非线性常微分方程得出数值解，较偏微分方程的数值解精度要高，在这个意义下，该数值解一般认为是精确解。

N-S方程至今尚未得到其通解，只有在一些特定的流动条件下，如流场中对流项等于零或可线性化，或偏微分方程可转化为常微分方程，才有可能得到其精确解。

一般说来，存在精确解的情况和实际流动情况相比，都做了一定程序的简化。

讨论精确解的意义主要在于：1、如果实际流动与精确解的流动与精确解的流动情况相近，可用摄动法求解流动问题，精确解构成这种方法的基础；2、用来校验数值计算的结果，如编写计算机程序进行数值计算时，以精确解作为算例，来验证程序的正确程度；3、也可以用来校核测试仪器的精准度。

参考文献
[1] 董曾南.章梓雄.粘性流体力学.北京:清华大学出版社.2001.
[2] 陈矛章.粘性动力学基础.北京:高等教育出版社.1993.
[3] 易大义.计算方法.李有法.浙江:浙江出版社.2002
[4] 苏金明.阮沈勇. MA TLAB实用教程.电子工业出版社.2008。