基于matlab的龙格库塔法求解布拉修斯方程

Runge—Kutta法求解布拉修斯解摘要薄剪切层方程主要有三种解法，即相似解，非相似条件下对偏微分方程组的数值解和近似解。布拉修斯解是布拉修斯于1908年求出的，它是零攻角沿平板流动的相似解。本文用四阶Runge—Kutta法求解高阶微分方程的方法，并用matlab编程实现，求得了与实际层流边界层相符合的数值解。关键词：布拉修斯解，相似解，Runge—Kutta法，数值解。1布拉修斯近似解方程二维定常不可压缩层流边界层的方程为：0yvxu（1）22yuvdxdyuvxuuuuee（2）边界条件为:0y)(,0xvuvw:y)(xuue将式（1）和式（2）进行法沃克纳—斯坎变换（简称F—S变换），将边界层方程无量纲化，即设yxvue5.0（3）xx（4）得出F—S变换后的动量方程xffxffxfmffmftk221211（5）其中k为流动类型指标，横曲率项t为212120cos211euvxLrLt（6）m是量纲一的压力梯度参数，定义为xdduuxmee（7）其边界条件变为:00f:1f对于二维平面实壁流动（:00wf）可以忽略横曲率项t的轴对称流动，式（5）成为xffxffxfmffmf2121（8）根据相似解的定义，方程（8）中的函数f若式相似的，则它应只与η有关而与x无关，即对x的偏导数应为零。于是方程（8）应成为01212fmffmf（9）若fw为常数，则方程（9）的边界条件为:0常数wff；0wff:1f2布拉修斯解布拉修斯于1908年求出了零攻角沿平板流动的解。这时0mue常数；因而方程（9）成为021fff（10）此即布拉修斯方程。对于实壁，0wf，边界条件成为:00wff；0wff:1f3Runge—Kutta法求解Runge—Kutta通过将高阶微分方程化为一阶线性方程组，从而解出高阶方程的数值解。在方程（10）中令0fddfdfdfddfddffff2221213（11）于是方程（10）变为313210333321022232101121),,,(),,,(),,,(ffffffTddffffffTddffffffTddf（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,hKfhKfhKfhKfTKKhfKhfKhfKhfTKKhfKhfKhfKhfTKffffTKKKKKhffiiiiiiiiiiiiiiiiiiiiiiiiiiiijiji（13）使用matlab软件取步长为0.2，迭代100步视作η→无穷大。迭代到第40步左右就收敛了，迭代结果如下（本文附录有全程序源代码）表格1平板层流边界层方程的数值解vxuyefeuuff0000.332060.20.0066410.0664080.331990.40.026560.132770.331470.60.0597360.198940.330080.80.106110.264710.3273910.165570.329780.323011.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.61.07250.772460.206462.81.2310.811510.1840131.39680.846050.161363.21.56910.876090.139133.41.7470.901770.117883.61.92950.923330.0980873.82.1160.941120.08012642.30580.955520.0642354.22.49810.966960.0505214.42.69240.975880.0389744.62.88830.982690.0294854.83.08530.987790.02187353.28330.991550.0159085.23.48190.994250.0113435.43.68090.996160.0079295.63.88030.997480.0054345.84.07990.998380.0036564.27960.998980.0024036.24.47950.999370.0015516.44.67940.999620.0009826.64.87930.999770.0006096.85.07930.999870.0003775.27930.999930.0002217.25.47930.999960.0001297.45.67930.999987.38E-057.65.87930.999994.15E-057.86.079312.28E-0586.279311.23E-058.26.479316.52E-068.46.679313.38E-068.66.879311.72E-068.87.079318.59E-0797.279314.20E-07由上表可以看出，数值解与布拉修斯解符合程度相当好。4结语(1)布拉修斯用相似变换将N—S方程简化，简化为一维微分方程求解并得到了与实际层流现象相符合的结果。（2）Runge—Kutta方法用来求解高阶微风方程，有相当高的精度。参考文献：[1]陈懋章.粘性流体力学基础.北京.高等教育出版社,2002.[2]复旦大学数学系《微分方程及其数值解》编写组.《微分方程及其数值解》.上海.复旦大学出版社.2001年[3]清华大学工程力学系编.流体力学基础：上册，下册.北京：机械工业出版社，1980年附录源程序代码如下%thisfileistosettlethebulaxiusifunctionwiththemethodof%Runge-Kutta.%f1standsforthefunctionf%f2standsforthefunctiondf1/du%f3standsforthefunctiondf2/duf(1:3,1:100)=0;%取A的1,3行，1,100列。f(3,1)=0.33206;x(1:101)=0;%hstandsforthesteplength;h=0.2;%k1,k2,k3,k4standsforthecoefficientofRung_Kuttaoff1k=[0,0,0,0;0,0,0,0;0,0,0,0];fori=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