《四阶龙格库塔求解边界层问题C语言》由会员分享，可在线阅读，更多相关《四阶龙格库塔求解边界层问题C语言（9页珍藏版）》请在金锄头文库上搜索。
1、题目：设/1.1 105m2/s的气体以v 10m/s的速度以零攻角定常扰流长度为L = 1m的大平板，用数值解讨论边界层内的流动规律。如图所示，沿板面方向取x坐标，板的法线方向取y坐标。在流体力学中，介绍了存在相似性解的二维平面不可压缩流体的边界层方程:C.E： 0 x y 2 u u 1 dp uM.E： u2x y dx y上0 y边界条件为：y = 0; u = v = 0为x点处壁面势流流速cy =8； u =v-( u = u (x)在外边界上pe 1 v2pe 1 v222所以dpe 0 0dx对于平板扰流，则平板扰流的边界层方程可简化为C.E： 0yuM.E: u x2 u u
2、2 y y其定解的边界条件为y = 0; u = v = 0y =8； u =v-为了求解边界层方程，我们可以利用“相似性解”的方法，对其进行转化, 由C.E,引进流函数 （x,y）,令由M.E式的偏导阶次，可望减少自变量数目y _y_2x / g(X) v-f()V g(x)其中一一 x 2xy g(x)由 vg(x)f,得 u v g(x) f V y yv v g(x)f x x所以，x u y(v f)fx2x. v .(v f )f yg(x)f)v g2(x)Ifv g(x) ( f2xf)将其都代入M.E,化简可得f ff 0这就使原方程变化为:f ff 0此时的边界条件为：4=
3、0； f(0) = 0 , f (0= 04=oo； f ( c=1那么，如何利用编辑程序的方法求解这个变化后的边界层微分方程呢?一.解方程的基本思路为了简化运算，此时边界层微分方程化简成：f ff 0,边界条件不变。以下为本次计算的基本思路：步骤1:利用数学代换将f ff 0转化为F FF 0,使其定解条件为F(0) 0, F (0) 0, F (0) 1步骤2:利用标准4阶龙格-库塔法，叠代解出的F FF 0中的F()的值步骤3:利用f( ) A2/3F( ) 1,即可计算出A F ( ) 3/2ff 0的定解条件步骤4:通过f (0) AF(0)换算出f(0),即可以将ff (0)0,
4、f (0)0, f ( )1 转化成f(0)0, f (0) 0, f (0) A步骤5:再次利用标准4阶龙格-库塔法，叠代计算出定解条件为f (0)0, f (0)0, f (0)A 时 fff 0不同时 f( ),f( )ffif ()的值二.编程前的准备工作1 .边界层方程的转换平板边界层流动问题，则原方程变为 f ff 0令 f(0) A,引入A1/3 ,F( ) f( )A1/3,由 fFA1/3则 f A1/3 F A1/3 F A1/3 A2/3 F f A2/3F AFAFa4/3F代入原方程得：A4/3 F (A1/3 F )(AF )其中：”=0得产04-oo得亚8由 f
5、(0)A2/3 F (0) 0 F (0) 0_1/3 _f (0) A F(0) 0F(0) 0f (0) AF (0) A F (0)1_2/3 _ _ 3/2f ( ) A F ( ) 1 A F()由此得到新的方程：F FF 00 (即 F” FF 0)边界条件为：F(0) 0, F (0) 0, F (0)运用龙格-库塔法先求解F=F（由并且，由推导过程可知:_ _ 3/2f (0) A F()2 .使用龙格-库塔法时方程的转换首先，简单介绍下龙格-库塔法基本思想:设 dx/dt=f(t,x,y,z)dy/dt= g(t,x,y,z)dz/dt= h(t,x,y,z)已知初值 t=t
6、0,x(t0)=x0,y(t0)=y0,z(t0)=z0那么利用标准4阶龙格-库塔法：xn+1=xn+1/6(k1+2k2+2k3+k4)yn+1 =yn+1/6(l1+2l2+2l3+l4)zn+1 =zn+1/6(m1+2m2+2m3+m4)则其中 y 的误差函数为： y= yn+i- yn=1/6(l 1+212+213+14)其中：kl = At*f(t n ,Xn,yn,Zn)k2= t*f(tn+At/2,Xn+ki/2,yn+1i/2,Zn+mi/2)k3= t*f(tn+At/2,Xn+k2/2,yn+12/2,Zn+m2/2)k4= At*f(t n+At ,Xn+k3,yn
7、+13,Zn+m3)1l=At*h(tn ,Xn,yn,Zn)12=At*h(tn+At/2,Xn+ki/2,yn+1i/2,Zn+mi/2)13=At*h(tn+At/2,Xn + k2/2,yn+12/2,Zn+m2/2)14=At*h(tn+At ,Xn+k3,yn+13,Zn+m3)mi = t*g(tn ,Xn,yn,Zn)m2= At*g(t n+ t/2,xn+ki/2,yn+1 i/2,zn+mi/2)m3= At*g(t n+ t/2,Xn+k2/2,yn+12/2,Zn+m2/2)m4= At*g(t n+ At ,Xn+k3,yn+13,Zn+m3)然后，以下是龙格-库塔
8、法时方程F FF 0的转换：则 F y xF y zF y z由此我们得出了 f(t,X,y,Z)、g(t,X,y,Z)和h(t,X,y,Z)的具体函数，从而容易得到 此时 ki. k2 .k3 .k4,1i.12, 13 .14 和 mi. m2 .m3 .m4的表达式。即化为平板边界层流动问题时，则：dX/dt=f(X,y,Z,t)=ydy/dt=g(X,y,Z,t) =zdZ/dt=h(X,y,Z,t) =-xz三.编程计算1,编程的简单思路利用C+编辑一个龙格-库塔法的叠代计算程序，通过两次调用此方法依次 计算出：定解条件为 F(0) 0, F (0) 0, F (0) HF，FF 0
9、的F()值；求解定解条件为f(0) 0, f (0) 0, f (0) A F ( ) 3/2时方程 f ff 0不同时共小)f(和f (的值。2 .龙格-库塔法的叠代计算程序的编程流程图/步长t的设定/精度的设定/ f(t,x,y,z)和 g(t,x,y,z)函数3 .C+程序#include stdio.h#include math.h#define m 0#define t 0.05#define l 0.00001double f(int i,double *p,int j,double x) double y;if(i =0) y=*(p+j);else if(i=1|i=2) y=
10、 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-
11、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=F (冏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,c12,
12、d3;/ b是 x,y,z 的存储数组，初值 0,0,1/c口为所有ki,li,mi的存储数组 TC: tc+;e=0;TURN: e+;for(i=0;i4;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);c3*i=k1;c3*i+1=k2; c3*i+2=k3;for(j=0;j3;j+)/x,y,z的计算输出不同刀时x,y,z值/dj=cj+2*cj+3+2*cj+6+cj+9;for(s=0;sl) goto TURN;if(tc=1) A=charge(b);/转化并再次调用计算程序/b0=0;b1=0;b2=A; goto TC; /输出叠代次数/printf(e=%dn,e);四.计算结果及分析以下是运行C+程序后的计算结果:,091081 .ftBBdllS015081276i3344 4444 一二二i=.1 = 一二-= nnnMnnnftnn14333=8.083322Mn?5一模 097G9&时59储.64477 ,083861g?56437Z31IE099- -33333333333333 JUUL二二二
