Runge —Kutta 法求解布拉修斯解摘要薄剪切层方程主要有三种解法，即相似解，非相似条件下对偏微分方程组的数值解和近似解。

布拉修斯解是布拉修斯于1908年求出的，它是零攻角沿平板流动的相似解。

本文用四阶Runge —Kutta 法求解高阶微分方程的方法，并用matlab 编程实现，求得了与实际层流边界层相符合的数值解。

关键词：布拉修斯解，相似解，Runge —Kutta 法，数值解。

1 布拉修斯近似解方程二维定常不可压缩层流边界层的方程为：0=∂∂+∂∂yvx u （1）22yuv dx d y u v x u u u u e e ∂∂+=∂∂+∂∂ （2）边界条件为:0=y )(,0x v u vw==:δ=y )(x u u e =将式（1）和式（2）进行法沃克纳—斯坎变换（简称F —S 变换），将边界层方程无量纲化，即设y x v u e 5.0⎪⎪⎭⎫⎝⎛=η （3）x x = （4）得出F —S 变换后的动量方程()[]()[]⎪⎭⎫ ⎝⎛∂∂''-∂'∂'='-+''++'''+x f f xf f x f m f f m f t k221211 （5）其中k 为流动类型指标，横曲率项t 为212120cos 211⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛⎪⎪⎭⎫ ⎝⎛++-=ηφe u vx L r L t （6） m 是量纲一的压力梯度参数，定义为xd du u x m ee =（7）其边界条件变为:0=η 0='f:∞=η 1='f对于二维平面实壁流动（:0=η0=w f ）可以忽略横曲率项t 的轴对称流动，式（5）成为()[]⎪⎭⎫ ⎝⎛∂∂''-∂'∂'='-+''++'''x f f xf f x f m f f m f 2121 （8） 根据相似解的定义，方程（8）中的函数f 若式相似的，则它应只与η有关而与x 无关，即对x 的偏导数应为零。

于是方程（8）应成为()[]01212='-+''++'''f m f f m f （9） 若f w 为常数，则方程（9）的边界条件为:0=η 常数==w f f ； 0='='wf f :∞=η 1='f2 布拉修斯解布拉修斯于1908年求出了零攻角沿平板流动的解。

这时 0==m u e 常数； 因而方程（9）成为021=''+'''f f f （10） 此即布拉修斯方程。

对于实壁，0=w f ，边界条件成为:0=η 0==w f f ； 0='='wf f :∞=η 1='f3 Runge —Kutta 法求解Runge —Kutta 通过将高阶微分方程化为一阶线性方程组，从而解出高阶方程的数值解。

在方程（10）中令η=0f⎪⎪⎪⎩⎪⎪⎪⎨⎧=====ηηηηd df d f d f d df d df f f f 2221213 （11） 于是方程（10）变为⎪⎪⎪⎩⎪⎪⎪⎨⎧-======313210333321022232101121),,,(),,,(),,,(f f f f f f T d df f f f f f T d df f f f f f T d df ηηη （12） 当区步长为h ，有四阶Runge —Kutta 的形式如下()⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧++++=++++=++++==++++=+),,,()2,2,2,2()2,2,2,2(),,,(2263,3,33,2,23,1,13,0,04,2,3,32,2,22,1,12,0,03,1,3,31,2,21,1,11,0,02,,3,2,1,01,4,3,2,1,,1,hK f hK f hK f hK f T K K h f K h f K h f K h f T K K h f K h f K h f K h f T K f f f f T K K K K K h f f i i i i i i i i i i i i i i i i i i i i i i i i i i i i j i j i （13）使用matlab 软件取步长为0.2，迭代100步视作η→无穷大。

迭代到第40步左右就收敛了，迭代结果如下（本文附录有全程序源代码）表格 1平板层流边界层方程的数值解1.20.237950.393780.316591.40.322990.456270.307871.60.420330.516760.296671.80.529520.574760.2829320.650030.629770.266752.20.78120.681320.248352.40.92230.728990.228092.6 1.07250.772460.206462.8 1.2310.811510.184013 1.39680.846050.161363.2 1.56910.876090.139133.4 1.7470.901770.117883.6 1.92950.923330.0980873.8 2.1160.941120.0801264 2.30580.955520.0642354.2 2.49810.966960.0505214.4 2.69240.975880.0389744.6 2.88830.982690.0294854.8 3.08530.987790.0218735 3.28330.991550.0159085.2 3.48190.994250.0113435.4 3.68090.996160.0079295.6 3.88030.997480.0054345.8 4.07990.998380.003656 4.27960.998980.0024036.2 4.47950.999370.0015516.4 4.67940.999620.0009826.6 4.87930.999770.0006096.8 5.07930.999870.000377 5.27930.999930.0002217.2 5.47930.999960.0001297.4 5.67930.999987.38E-057.6 5.87930.99999 4.15E-057.8 6.07931 2.28E-058 6.27931 1.23E-058.2 6.47931 6.52E-068.4 6.67931 3.38E-068.6 6.87931 1.72E-068.87.079318.59E-0797.27931 4.20E-07由上表可以看出，数值解与布拉修斯解符合程度相当好。

4 结语(1) 布拉修斯用相似变换将N—S方程简化，简化为一维微分方程求解并得到了与实际层流现象相符合的结果。

（2）Runge—Kutta方法用来求解高阶微风方程，有相当高的精度。

参考文献：[1] 陈懋章.粘性流体力学基础.北京.高等教育出版社,2002.[2] 复旦大学数学系《微分方程及其数值解》编写组.《微分方程及其数值解》.上海.复旦大学出版社.2001年[3] 清华大学工程力学系编.流体力学基础：上册，下册.北京：机械工业出版社，1980年附录源程序代码如下% this file is to settle the bulaxiusi function with the method of% Runge-Kutta.% f1 stands for the function f% f2 stands for the function df1/du% f3 stands for the function df2/duf(1:3,1:100)=0; %取A的1,3行，1,100列。

f(3,1)=0.33206;x(1:101)=0;% h stands for the step length;h=0.2;% k1, k2,k3,k4 stands for the coefficient of Rung_Kutta of f1k=[0,0,0,0;0,0,0,0;0,0,0,0];for i=1:100;k(1,1)=f(2,i);k(2,1)=f(3,i);k(3,1)=-f(1,i)*f(3,i)/2k(1,2)=f(2,i)+h*k(2,1)/2;k(2,2)=f(3,i)+h*k(3,1)/2;k(3,2)=-(f(1,i)+h*k(1,1)/2)*(f(3,i)+h*k(3,1)/2)/2;k(1,3)=f(2,i)+h*k(2,2)/2;k(2,3)=f(3,i)+h*k(3,2)/2;k(3,3)=-(f(1,i)+h*k(1,2)/2)*(f(3,i)+h*k(3,2)/2)/2;k(1,4)=f(2,i)+h*k(2,3);k(2,4)=f(3,i)+h*k(3,3);k(3,4)=-(f(1,i)+h*k(1,3))*(f(3,i)+h*k(3,3))/2;f(1,i+1)=f(1,i)+h*(k(1,1)+2*k(1,2)+2*k(1,3)+k(1,4))/6;f(2,i+1)=f(2,i)+h*(k(2,1)+2*k(2,2)+2*k(2,3)+k(2,4))/6;f(3,i+1)=f(3,i)+h*(k(3,1)+2*k(3,2)+2*k(3,3)+k(3,4))/6;x(i+1)=x(i)+h;end。