四阶龙格库塔求解边界层问题(C语言)

题目：设的气体以的速度以零攻角定常扰流长度为L＝1m的大平板，用数值解讨论边界层内的流动规律。如下图，沿板面方向取x坐标，板的法线方向取y坐标。在流体力学中，介绍了存在相似性解的二维平面不可压缩流体的边界层方程：C.E：M.E:边界条件为：y=0；u=v=0y=∞；u=v∞(u=u(x)为x点处壁面势流流速)在外边界上所以对于平板扰流，那么平板扰流的边界层方程可简化为C.E：M.E:其定解的边界条件为y=0；u=v=0y=∞；u=v∞为了求解边界层方程，我们可以利用“相似性解〞的方法，对其进行转化，由C.E，引进流函数，令由M.E式的偏导阶次，可望减少自变量数目令其中由，得所以，将其都代入M.E，化简可得这就使原方程变化为：此时的边界条件为：η=0；f(0)=0,f’(0)=0η=∞；f’(∞)=1那么，如何利用编辑程序的方法求解这个变化后的边界层微分方程呢？解方程的根本思路为了简化运算，此时边界层微分方程化简成：，边界条件不变。以下为本次计算的根本思路：步骤1：利用数学代换将转化为，使其定解条件为步骤2：利用标准4阶龙格-库塔法，叠代解出的中的的值步骤3:利用，即可计算出步骤4：通过换算出，即可以将的定解条件转化成步骤5：再次利用标准4阶龙格-库塔法，叠代计算出定解条件为时不同η时的值二.编程前的准备工作⒈边界层方程的转换平板边界层流动问题，那么原方程变为令，引入,由那么代入原方程得：〔即〕其中：η=0，得µ=0η→∞，得µ→∞由由此得到新的方程：边界条件为：运用龙格-库塔法先求解F=F(µ)并且，由推导过程可知：⒉使用龙格-库塔法时方程的转换首先，简单介绍下龙格-库塔法根本思想：设dx/dt=f(t,x,y,z)dy/dt=g(t,x,y,z)dz/dt=h(t,x,y,z)初值t=t0,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+1-yn=1/6(l1+2l2+2l3+l4)其中：k1=△t*f(tn,xn,yn,zn)k2=△t*f(tn+△t/2,xn+k1/2,yn+l1/2,zn+m1/2)k3=△t*f(tn+△t/2,xn+k2/2,yn+l2/2,zn+m2/2)k4=△t*f(tn+△t,xn+k3,yn+l3,zn+m3)l1=△t*h(tn,xn,yn,zn)l2=△t*h(tn+△t/2,xn+k1/2,yn+l1/2,zn+m1/2)l3=△t*h(tn+△t/2,xn+k2/2,yn+l2/2,zn+m2/2)l4=△t*h(tn+△t,xn+k3,yn+l3,zn+m3)m1=△t*g(tn,xn,yn,zn)m2=△t*g(tn+△t/2,xn+k1/2,yn+l1/2,zn+m1/2)m3=△t*g(tn+△t/2,xn+k2/2,yn+l2/2,zn+m2/2)m4=△t*g(tn+△t,xn+k3,yn+l3,zn+m3)然后，以下是龙格-库塔法时方程的转换：那么由此我们得出了f(t,x,y,z)、g(t,x,y,z)和h(t,x,y,z)的具体函数，从而容易得到此时k1.k2.k3.k4,l1.l2，l3.l4和m1.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三.编程计算⒈编程的简单思路利用C++编辑一个龙格-库塔法的叠代计算程序，通过两次调用此方法依次计算出：⑴定解条件为时的值；⑵求解定解条件为时方程不同η时f(η)、f’(η)和f〞(η)的值。⒉龙格-库塔法的叠代计算程序的编程流程图⒊C++程序#include"stdio.h"#include"math.h"#definem0#definet0.05//步长t的设定//#definel0.00001//精度的设定//doublef(inti,double*p,intj,doublex)//f(t,x,y,z)和g(t,x,y,z)函数//{doubley;if(i==0)y=*(p+j);elseif(i==1||i==2)y=x*t/2+*(p+j);elsey=x*t+*(p+j); returny;}doubleg(inti,doublek,double*p,doublex,doubley,doublez)//h(t,x,y,z)函数//{doubleg,a,b,c; if(i==0) {a=*(p+0); b=*(p+1); c=*(p+2);} elseif(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-a*c; returng;}doublediffer(intn,double*p)//误差函数△y//{doubley;y=*(p+n)*t/6;returny;}doublecharge(double*p)//转化函数f’’(0)=A=[F’(∞)]-3/2//{doubley,a; a=-1.5;y=pow(*(p+1),a); returny;}voidmain(){doublek1=0,k2=0,k3=0;//k1,k2,k3分别为ki,li,mi// doubleq,v,w,A; inti,j,s,u,h,e,tc=0;//tc为程序运行次数////e为叠代次数//doubleb[]={0,0,1},c[12],d[3];//b[]是x,y,z的存储数组，初值0,0,1////c[]为所有ki,li,mi的存储数组//TC:tc++;e=0;TURN:e++;for(i=0;i<4;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);c[3*i]=k1;c[3*i+1]=k2;c[3*i+2]=k3;}for(j=0;j<3;j++){d[j]=c[j]+2*c[j+3]+2*c[j+6]+c[j+9];}for(s=0;s<3;s++)//x,y,z的计算//{b[s]=b[s]+differ(s,d);}if(tc==2)//输出不同η时x,y,z值//{printf("n=%lf",t*e);printf("f=%lf",*b);printf("f`=%lf",*(b+1));printf("f``=%lf\n",*(b+2));}h=1;q=differ(h,d);//精度控制//if(fabs(q)>l)gotoTURN;if(tc==1)//转化并再次调用计算程序//{A=charge(b); b[0]=0;b[1]=0;b[2]=A; gotoTC;}printf("e=%d\n",e);//输出叠代次数//}四.计算结果及分析以下是运行C++程序后的计算结果：由计算数据可以得到〔当精度为10-5时〕：⑴当0<η<5.1时随着η的增大，