四阶龙格库塔求解边界层问题.doc

第四题：设x的气体以的速度以零攻角定常绕流长度为L=1m的大平板，试用数值解讨论边界层内的流动规律。解：如图1所示，由于，所以可以抽象为高雷诺数气流绕平板的边界层定常流动问题。沿板面取气体流向为x轴正向，板的法线方向取y坐标，记边界层厚度，因厚度很薄，既有。该问题可按二维模型处理，其求解过程如下： 图1 绕流大平板的边界层模型一、流体动力学基本方程根据问题的特性，忽略质量力，设边界层内的速度矢量为，压强，则流体动力学基本方程具体化为：C.E （4.1）M.E. （4.2） （4.3）下面，用数量级比较的方法分析方程中各项的数量级关系：在此问题中，自变量为：x （0，L），y （0，）；因变量为：u （0，ue）；v （0，vmax）； 无量纲作变换： ； ；将各无量纲变量代入连续性方程4.1，则：1 1 可知A与B应属于同一个数量级;将各无量纲变量代入动量方程M.E(x方向)4.2，则： 1 1 1 1 1将各无量纲变量代入动量方程M.E(y方向)4.3，则： 1 1 1 1 1在连续性方程中，去掉任何一项都缺乏物理依据和工程事实，故应考虑两项具有相同的数量级，即有A=B；将些关系带入分析x方向的动量方程，发现左端两项与压力梯度项具有相同的数量级，两个粘性扩散项中，第一项明显具有很低的数量级，可以去掉，由于边界层中粘性的作用不容忽视，第二项应予保留，即认为该项也具有压力梯度项相同的数量级，则有；分析y方向的动量方程，发现对流项与粘性扩散项具有相同的数量级，而压力梯度项却具有高得多的数量级。去掉各方程中具有明显较低数量级的项，得到简化的近似方程组（边界层方程）：C.E M.E 由上式可知，在任一过流断面上，边界层内各点的压力与其边界上的势流压力相等，即，而边界层外势流区满足：既有：因此，边界层方程即可表示为：C.E M.E. 对于零攻角绕流平板流动：Ve不随x变化，则平板绕流边界层方程可以简化为：C.E M.E. 其定解的边界条件为：上述方程及边界条件即为零攻角定常绕流平板的边界层控制方程及定解条件。二、将控制方程化为常微分方程下面采用无量纲相似性解法：（1）引入流函数（x，y），将未知函数的数量由两个化为一个：令；使连续性方程自动满足，动量方程变形为：（2）作相似性变换：引入无量纲相似性变量令，。则有： 上式可变为常微分方程：相应的定解条件为：（3）将初边值问题变换为初值问题：令（这里，A为非零常数）引入新变量，；则上式可变为： ；定解条件相应的变换为： ；其中，三、数值计算以及编写程序（1）数值计算运用标准四阶龙格库塔法先求解,得到后，再次运用标准四阶龙格库塔法求解，并得到及后，可以直接得到关于相似性自变量的流函数，再通过转换，即可求得所需的物理量的函数关系。x方向速度：；y方向速度：；边界层厚度：，；板面切应力：；运用标准四阶龙格库塔法求解和的具体步骤如下：首先令 其中，求解时C=1；求解时C=；若令，原方程进一步化为： 取步长为，由上可知：，带入四阶龙格-库塔公式得： 式中，（），这样就可以求出和的值。（2）程序源代码本作业采取Matlab7.5作为计算软件，共编制五个文件：RGKT.m，MyFunction1.m，f.m，u.m,v.m为主程序，主要负责迭代计算和输出显示；F、G、H经转换成用龙格-库塔方法计算的标准形式后，储存于MyFunction1.m中供前者调用。由于方程形式上相同，故采用同一套迭代程序求解。主程序RGKT.m的内容为:function result = RGKT(varargin)%=% 调用格式:RGKT(MyFunction,f(1)_0,f(2)_0,f(3)_0,.,t_0,t_1,t_span);% MyFunction为函数名% f(1)_0,f(2)_0,f(3)_0,.为变量f(1),f(2),f(3),.的初始值% t_0为求解区间初始值% t_1为求解区间终值% t_span为步长%=F=varargin1;%=把调用函数名赋给Ff0=varargin2;%=把初值赋给f0t_0=varargin3;%=把求解区间起点赋给t_0t_1=varargin4;%=把求解区间终点赋给t_1t_span=varargin5;%=把求解步长赋给t_spant=t_0;%=把求解区间起点值赋给循环变量tyita=(t:t_span:t_1);%=保存自变量序列，以供结果显示n=length(t:t_span:t_1);%=求出迭代次数m=length(feval(F,f0,t);%=m表示方程个数(也就是因变量个数)x=zeros(n,m);%=生成并初始化解矩阵x(1,:)=f0;%=把初值赋给解矩阵第一行x_temp=x(1,:);%=把第一行(初值)赋给循环中变量for i=2:nx_0=x_temp;%=得x_0的值k1=feval(F,x_temp,t);%=根据初值和时间求K1k1=reshape(k1,1,m);%=把K1转换成m列的行向量x_temp=x_0+k1*(t_span/2);%=重新计算自变量；t=t+t_span/2;%=重新计算时间k2=feval(F,x_temp,t);%=依照重新计算的自变量和时间计算K2k2=reshape(k2,1,m);%=把K2转换成m列的行向量x_temp=x_0+k2*(t_span/2);%=重新计算自变量k3=feval(F,x_temp,t);%=依照重新计算的自变量K2k3=reshape(k3,1,m);%=把K2转换成m列的行向量x_temp=x_0+k3*(t_span/2);%=重新计算自变量；t=t+t_span/2;%=重新计算时间；k4=feval(F,x_temp,t);%=依照重新计算的自变量和时间计算K4k4=reshape(k4,1,m);%=把K4转换成1行m列的列向量x(i,:)=x_0+t_span*(k1+2*k2+2*k3+k4)/6;%=计算解矩阵中第i行的值x_temp=x(i,:);%=把本次计算结果赋给循球中变量作为下次计算的初值endresult=yita,x;%=把自变量和结果重整成一个矩阵作为结果返回end%=结果显示t=t_0:t_span:t_1;length(t);h=plot(t,x(:,1),r,t,x(:,2),g,t,x(:,3),k);grid on;h1=legend(h,第一个变量,第二个变量,第三个变量),2);set(h1,FontSize,10);set(h1,FontName,仿宋); 原方程可以通过代换转换成用龙格-库塔方法求解的标准形式：在程序中，令,，,（f和x为向量，,),）,采用矩阵形式表示该方程组，即有：。在程序中，单独建立一文件MyFunction1.m储存表示该方程组的函数。MyFunction1.m内容为：function f=MyFunction1(x,yita)f= zeros(3,1);%初始化列向量f(1) =x(2);f(2) = x(3);f(3) = -x(1)*x(3);四、计算过程及结果分析（1），和的变化在命令窗口中输入：f(MyFunction1,0,0,1,0,9,0.05)第一次解算得出的图形如1所示（红色为，绿色为，黑色为）:图1绿色线表示的迭代结果，由图中绿线变化趋势以及由迭代结果矩阵中数据可认为=1.6639，由此可算出=0.4659。图2在命令窗口中输入：f(MyFunction1,0,0,0.4659,0,9,0.05)第二次求解结果如下图（红色为，绿色为，黑色为）如图2。（2）X方向速度，其中。在命令窗口中输入：u(MyFunction1,0,0,0.4659,0,9,0.05)如上图中绿色线示，在大于3之后逐渐趋于1；图3为的分布示意图（横坐标为速度大小，纵坐标为无纲量距离）：图3 的分布示意图从图中可以看，边界层内上升的速率很快。在，也就是如果的话，只要约等于就能达到主流速度。（3）Y方向速度用相应离散点的值可计算出处的分布,令x=0.25m的v分布（如图4，横轴表示速度值，纵坐标表示无纲量距离）。