四阶龙格库塔求解边界层问题(C语言)

题目： 设52
/ 1.110/m s μρ-=⨯的气体以10/v m s ∞=的速度以零攻角定常扰流长
度为L ＝1m 的大平板，用数值解讨论边界层内的流动规律。

如图所示，沿板面方向取x 坐标，板的法线方向取y 坐标。

在流体力学中，介绍了存在相似性解的二维平面不可压缩流体的边界层方程：
C.E ：
0=∂∂+∂∂y
x u υ
M.E : 221y u
dx dp y u x u u ∂∂+-=∂∂+∂∂νρυ
0p y
∂=∂ 边界条件为：y = 0；u = v = 0
y =∞；u =v ∞ ( u = u (x) 为x 点处壁面势流流速 )
在外边界上22
1122e e p v p v c ρρ∞+=+=
所以 00e dp
dx
+=
对于平板扰流，则平板扰流的边界层方程可简化为 C.E ：
0=∂∂+∂∂y
x u υ M.E : 22u u u
u x y y
υ
ν∂∂∂+=∂∂∂ 其定解的边界条件为
y = 0；u = v = 0
y =∞；u =v ∞
为了求解边界层方程，我们可以利用“相似性解”的方法，对其进行转化，由C.E ，引进流函数),(y x ψ，令
x
y u ∂∂-=∂∂=
ψυψ, 由M.E 式的偏导阶次，可望减少自变量数目
令()y
g x η==
()()f f v g x η∞=
= 其中
2x x ηη∂=-∂
1()
y g x η∂=∂ 由()v g x f ψ∞=，得()'v g x f
u v f y y
ψ∞∞∂∂=
==∂∂ ()()
(')2v g x f v g x v f f x x x
ψη∞∞∂∂=-
==-∂∂ 所以，
(')"2v u v f f x x x η∞∞∂∂
==-∂∂
(')"()
v u v f f y y g x ∞∞∂∂
==∂∂
222(")"'()()
v v u f f y y g x g x ∞
∞∂∂==∂∂ 将其都代入M.E ，化简可得"'"0f ff += 这就使原方程变化为：
"'"0f ff +=
此时的边界条件为：η = 0；f(0) = 0 , f’(0) = 0
η = ∞；f’(∞) = 1
那么，如何利用编辑程序的方法求解这个变化后的边界层微分方程呢？
一. 解方程的基本思路
为了简化运算，此时边界层微分方程化简成："'"0f ff +=，边界条件不变。

以下为本次计算的基本思路：
步骤1：利用数学代换将"'"0f ff +=转化为"'"0F FF +=，使其定解条件为
(0)0,'(0)0,"(0)1F F F ===
步骤2：利用标准4阶龙格-库塔法，叠代解出的"'"0F FF +=中的'()F ∞的值 步骤3: 利用2/3'()'()1f A F ∞=∞=，即可计算出3/2['()]A F -=∞
步骤4：通过"(0)"(0)f AF =换算出"(0)f ，即可以将"'"0f ff +=的定解条件
(0)0,'(0)0,'()1f f f ==∞=转化成(0)0,'(0)0,"(0)f f f A ===
步骤5：再次利用标准4阶龙格-库塔法，叠代计算出定解条件为
(0)0,'(0)0,"(0)f f f A ===时"'"0f ff +=不同η时(),'()"()
f f f ηηη和的值
二. 编程前的准备工作
⒈ 边界层方程的转换
平板边界层流动问题，则原方程变为"'"0f ff +=
令"(0)f A =，引入1/31/3,()()A F f A μημη-==, 由1/3f FA = 则1/31/31/32/3'''''f A F A F A A F ημ===
2/3""'"f A F AF ημ== 4/3'""'''"f AF A F ημ==
代入原方程得：4/31/3'"()(")0A F A F AF +=（即'""0F FF +=） 其中：η=0，得µ=0
η→∞，得µ→ ∞
由2/3'(0)'(0)0'(0)0f A F F ==→=
1/3(0)(0)0(0)0f A F F ==→=
"(0)"(0)"(0)1f AF A F ==→= 2/33/2'()'()1['()]f A F A F -∞=∞=→=∞
由此得到新的方程：'""0F FF += 边界条件为：(0)0,'(0)0,"(0)1F F F === 运用龙格-库塔法先求解F=F(µ)
并且，由推导过程可知：3/2"(0)['()]f A F -==∞ ⒉ 使用龙格-库塔法时方程的转换 首先，简单介绍下龙格-库塔法基本思想： 设 dx/dt=f(t,x,y,z) dy/dt= g(t,x,y,z) dz/dt= h(t,x,y,z)
已知初值 t=t 0,x(t 0)=x 0,y(t 0)=y 0,z(t 0)=z 0 那么利用标准4阶龙格-库塔法： x n+1=x n +1/6(k 1+2k 2+2k 3+k 4) y n+1 =y n +1/6(l 1+2l 2+2l 3+l 4) z n+1 =z n +1/6(m 1+2m 2+2m 3+m 4)
则其中y 的误差函数为：△y= y n+1- y n =1/6(l 1+2l 2+2l 3+l 4) 其中：k 1=△t*f(t n ,x n ,y n ,z n )
k 2=△t*f(t n +△t/2,x n +k 1/2,y n +l 1/2,z n +m 1/2) k 3=△t*f(t n +△t/2,x n +k 2/2,y n +l 2/2,z n +m 2/2) k 4=△t*f(t n +△t ,x n +k 3,y n +l 3,z n +m 3) l 1=△t*h(t n ,x n ,y n ,z n )
l 2=△t*h(t n +△t/2,x n +k 1/2,y n +l 1/2,z n +m 1/2) l 3=△t*h(t n +△t/2,x n +k 2/2,y n +l 2/2,z n +m 2/2) l 4=△t*h(t n +△t ,x n +k 3,y n +l 3,z n +m 3) m 1=△t*g(t n ,x n ,y n ,z n )
m 2=△t*g(t n +△t/2,x n +k 1/2,y n +l 1/2,z n +m 1/2) m 3=△t*g(t n +△t/2,x n +k 2/2,y n +l 2/2,z n +m 2/2) m 4=△t*g(t n +△t ,x n +k 3,y n +l 3,z n +m 3)
然后，以下是龙格-库塔法时方程'""0F FF += 的转换： 则''F y x ==
"'F y z ==
"'"'F y z ==
由此我们得出了f(t,x,y,z) 、g(t,x,y,z) 和h(t,x,y,z)的具体函数，从而容易得到此时k 1. k 2 .k 3 .k 4,l 1 .l 2，l 3 .l 4 和m 1. m 2 .m 3 .m 4的表达式。

即化为平板边界层流动问题时， 则： dx/dt=f(x,y,z,t)=y dy/dt=g(x,y,z,t) =z
dz/dt=h(x,y,z,t) =-xz 三. 编程计算 ⒈ 编程的简单思路
利用C++编辑一个龙格-库塔法的叠代计算程序，通过两次调用此方法依次计算出：⑴定解条件为(0)0,'(0)0,"(0)1F F F ===时'""0F FF +=的'()F ∞值；
⑵求解定解条件为3/2(0)0,'(0)0,"(0)['()]f f f A F -====∞时方程
"'"0f ff +=不同η时f(η)、f’(η)和f” (η)的值。

⒉ 龙格-库塔法的叠代计算程序的编程流程图
⒊C++程序
#include "stdio.h"
#include "math.h"
#define m 0
#define t 0.05 //步长t的设定//
#define l 0.00001 //精度的设定//
double f(int i,double *p,int j,double x) // f(t,x,y,z) 和g(t,x,y,z)函数// { double y;
if(i==0) y=*(p+j);
else if(i==1||i==2)
y= x*t/2+*(p+j);
else y= x*t +*(p+j);
return y;
}
double g(int i,double k,double *p,double x,double y,double z) // h(t,x,y,z)函数// { double g,a,b,c;
if(i==0)
{ a=*(p+0);
b=*(p+1);
c=*(p+2); }
else if(i==1||i==2)
{ a= x*t/2+*(p+0);
b= y*t/2+*(p+1);
c= z*t/2+*(p+2);
}
else
{ a= x*t +*(p+0);
b= y*t +*(p+1);
c= z*t +*(p+2); }
g=k*b*b-k-a*c;
return g;
}
double differ(int n,double *p) //误差函数△y//
{ double y;
y=*(p+n)*t/6;
return y;
}
double charge(double *p) //转化函数f’’(0)=A=[F’(∞)]-3/2// { double y,a;
a=-1.5;y=pow(*(p+1),a);
return y;
}
void main()
{ double k1=0,k2=0,k3=0; //k1,k2,k3分别为ki,li,mi// double q,v,w,A;
int i,j,s,u,h,e,tc=0; //tc为程序运行次数//
//e为叠代次数//
double b[]={0,0,1},c[12],d[3]; // b[]是x,y,z的存储数组，初值0,0,1//
//c[]为所有ki,li,mi的存储数组// TC: tc++;
e=0;
TURN: e++;
for(i=0;i<4;i++) //ki,li,mi的计算和存储//
{ w=k3;
v=2*m/(m+1);k3=g(i,v,b,k1,k2,k3);
h=1;k1=f(i,b,h,k2);
h=2;k2=f(i,b,h,w);
c[3*i]=k1;
c[3*i+1]=k2;
c[3*i+2]=k3;
}
for(j=0;j<3;j++)
{ d[j]=c[j]+2*c[j+3]+2*c[j+6]+c[j+9]; }
for(s=0;s<3;s++) //x,y,z的计算//
{ b[s]=b[s]+differ(s,d); }
if(tc==2) //输出不同η时x,y,z值// { printf("n=%lf ",t*e);
printf("f=%lf ",*b);
printf("f `=%lf ",*(b+1));
printf("f ``=%lf\n",*(b+2));
}
h=1;q=differ(h,d); //精度控制//
if(fabs(q)>l) goto TURN;
if(tc==1) //转化并再次调用计算程序// { A=charge(b);
b[0]=0;b[1]=0;b[2]=A;
goto TC; }
printf("e=%d\n",e); //输出叠代次数//
}
四. 计算结果及分析
以下是运行C++程序后的计算结果：
由计算数据可以得到（当精度为10-5时）：
⑴当0<η<5.1时随着η的增大，f `(η)是逐渐增大的。

⑵当η>5.1以后，f `(η)已经收敛于1。

这是可以用流体力学的知识来解释的。

在流体绕壁面流动时，边界层内沿x 方向的速度u=U f `(η)，其中U为势流速度，即f `(η)=u/U。

我们知道，当η→∞时，u=U，即u/U=1。

所以，当η→∞时，f `(η)应该收敛于1。

另外，由实际程序计算可以得到，当η>5.1时绝对误差已经小于10-5，因此可以认为此时f `(η)已经收敛。