四阶龙格库塔求解边界层问题.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) %=====================================================