四階龍格庫塔求解邊界層問題(C語言)

題目：設(shè)的氣體以的速度以零攻角定常擾流長度為L＝1m的大平板，用數(shù)值解討論邊界層內(nèi)的流動規(guī)律。如下圖，沿板面方向取x坐標(biāo)，板的法線方向取y坐標(biāo)。在流體力學(xué)中，介紹了存在相似性解的二維平面不可壓縮流體的邊界層方程：C.E：M.E:邊界條件為：y=0；u=v=0y=∞；u=v∞(u=u(x)為x點(diǎn)處壁面勢流流速)在外邊界上所以對于平板擾流，那么平板擾流的邊界層方程可簡化為C.E：M.E:其定解的邊界條件為y=0；u=v=0y=∞；u=v∞為了求解邊界層方程，我們可以利用“相似性解〞的方法，對其進(jìn)行轉(zhuǎn)化，由C.E，引進(jìn)流函數(shù)，令由M.E式的偏導(dǎo)階次，可望減少自變量數(shù)目令其中由，得所以，將其都代入M.E，化簡可得這就使原方程變化為：此時(shí)的邊界條件為：η=0；f(0)=0,f’(0)=0η=∞；f’(∞)=1那么，如何利用編輯程序的方法求解這個(gè)變化后的邊界層微分方程呢？解方程的根本思路為了簡化運(yùn)算，此時(shí)邊界層微分方程化簡成：，邊界條件不變。以下為本次計(jì)算的根本思路：步驟1：利用數(shù)學(xué)代換將轉(zhuǎn)化為，使其定解條件為步驟2：利用標(biāo)準(zhǔn)4階龍格-庫塔法，疊代解出的中的的值步驟3:利用，即可計(jì)算出步驟4：通過換算出，即可以將的定解條件轉(zhuǎn)化成步驟5：再次利用標(biāo)準(zhǔn)4階龍格-庫塔法，疊代計(jì)算出定解條件為時(shí)不同η時(shí)的值二.編程前的準(zhǔn)備工作⒈邊界層方程的轉(zhuǎn)換平板邊界層流動問題，那么原方程變?yōu)榱?，引?由那么代入原方程得：〔即〕其中：η=0，得μ=0η→∞，得μ→∞由由此得到新的方程：邊界條件為：運(yùn)用龍格-庫塔法先求解F=F(μ)并且，由推導(dǎo)過程可知：⒉使用龍格-庫塔法時(shí)方程的轉(zhuǎn)換首先，簡單介紹下龍格-庫塔法根本思想：設(shè)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那么利用標(biāo)準(zhǔn)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的誤差函數(shù)為：△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)然后，以下是龍格-庫塔法時(shí)方程的轉(zhuǎn)換：那么由此我們得出了f(t,x,y,z)、g(t,x,y,z)和h(t,x,y,z)的具體函數(shù)，從而容易得到此時(shí)k1.k2.k3.k4,l1.l2，l3.l4和m1.m2.m3.m4的表達(dá)式。即化為平板邊界層流動問題時(shí)，那么：dx/dt=f(x,y,z,t)=ydy/dt=g(x,y,z,t)=zdz/dt=h(x,y,z,t)=-xz三.編程計(jì)算⒈編程的簡單思路利用C++編輯一個(gè)龍格-庫塔法的疊代計(jì)算程序，通過兩次調(diào)用此方法依次計(jì)算出：⑴定解條件為時(shí)的值；⑵求解定解條件為時(shí)方程不同η時(shí)f(η)、f’(η)和f〞(η)的值。⒉龍格-庫塔法的疊代計(jì)算程序的編程流程圖⒊C++程序#include"stdio.h"#include"math.h"#definem0#definet0.05//步長t的設(shè)定//#definel0.00001//精度的設(shè)定//doublef(inti,double*p,intj,doublex)//f(t,x,y,z)和g(t,x,y,z)函數(shù)//{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)函數(shù)//{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)//誤差函數(shù)△y//{doubley;y=*(p+n)*t/6;returny;}doublecharge(double*p)//轉(zhuǎn)化函數(shù)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為程序運(yùn)行次數(shù)////e為疊代次數(shù)//doubleb[]={0,0,1},c[12],d[3];//b[]是x,y,z的存儲數(shù)組，初值0,0,1////c[]為所有ki,li,mi的存儲數(shù)組//TC:tc++;e=0;TURN:e++;for(i=0;i<4;i++)//ki,li,mi的計(jì)算和存儲//{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的計(jì)算//{b[s]=b[s]+differ(s,d);}if(tc==2)//輸出不同η時(shí)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)//轉(zhuǎn)化并再次調(diào)用計(jì)算程序//{A=charge(b); b[0]=0;b[1]=0;b[2]=A; gotoTC;}printf("e=%d\n",e);//輸出疊代次數(shù)//}四.計(jì)算結(jié)果及分析以下是運(yùn)行C++程序后的計(jì)算結(jié)果：由計(jì)算數(shù)據(jù)可以得到〔當(dāng)精度為10-5時(shí)〕：⑴當(dāng)0<η<5.1時(shí)隨著η的增大，