資源描述：
《計算流體驅(qū)動方腔程序.doc》由會員上傳分享，免費在線閱讀，更多相關(guān)內(nèi)容在行業(yè)資料-天天文庫。
1、驅(qū)動方腔（粘性不可壓流），方腔示意圖如左下：物理背景：如左圖所示,裝滿黏性不可壓流體的方腔上底面以速度U運動,其它壁面固定,其內(nèi)部流體將作類似湍流的運動.流體運動復(fù)雜度將依賴于以下要素:(1)方腔形狀；(2)U的大??；(3)流體的黏性,即雷諾數(shù)Re的大小。本題只討論理想流體在二維正方形空腔中的流動，如左圖，作無量綱化處理后,方腔長度和上底面運動速度均取為1,通過差分?jǐn)?shù)值模擬計算出運動的流函數(shù)，渦量函數(shù)，畫出Re=100和200時的流線圖和渦量圖。并標(biāo)出窩心位置。采用流函數(shù)-渦方法：對于二維不可壓
2、縮的理想流體,流速應(yīng)滿足如下方程組:：（1）連續(xù)方程：（1）(2)N-S方程：（2）（3）引入流函數(shù)：和：，，（4）（4）代入（1）的流函數(shù)的方程：；（5）把（5）代入，的渦量方程：（6）以上方法稱之為流函數(shù)－渦方法。數(shù)值方法：將區(qū)域分為N×N個方格單元，只有（N+1）×（N+1）格網(wǎng)格點要計算，令h=1/N為空間步長，為時間步長。1：流函數(shù)與速度的邊界條件：2渦量的邊界條件：固壁上的渦量：計算：1，方腔內(nèi)的流函數(shù)：采用中心五點差分格式對（5）進(jìn)行離散并簡化得：算出后，對（4）采用中心差分離散可
3、以求得u,v.2,方腔內(nèi)得渦量：采用迎風(fēng)格式對渦量方程（6）進(jìn)行離散：程序?qū)崿F(xiàn)過程如下：clearallclcjishi=cputime;timestep=2000;[x,y]=meshgrid(0:1/100:1);Nx=length(x);Ny=length(y);dh=1/100;dt=0.6*dh;renold=200;fluidfun=zeros(Nx,Ny);vortex=zeros(Nx,Ny);velocityx=zeros(Nx,Ny);velocityy=zeros(Nx,N
4、y);%%%初始條件velocityx(:,Ny)=1;vortex(:,Ny)=-2/dh;%fluidfun(Ny,:)=0;%velocityy(Ny,:)=0;fort=1:timestepfori=2:Nx-1forj=2:Ny-1fluidfun(i,j)=(fluidfun(i+1,j)+fluidfun(i-1,j)+fluidfun(i,j+1)+fluidfun(i,j-1)+dh^2*vortex(i,j))/4;endendfori=2:Nx-1forj=2:Ny-1ve
5、locityx(i,j)=(fluidfun(i,j+1)-fluidfun(i,j-1))/(2*dh);velocityy(i,j)=(fluidfun(i-1,j)-fluidfun(i+1,j))/(2*dh);endendfori=2:Nx-1forj=2:Ny-1a=renold*(abs(velocityx(i,j))+abs(velocityy(i,j)))/dh+4/dh^2+renold/(dt);b=renold*abs(velocityx(i,j))/dh+1/dh^2;
6、c=renold*abs(velocityy(i,j))/dh+1/dh^2;if(velocityx(i,j)>=0)&&(velocityy(i,j)>=0)vortex(i,j)=((vortex(i+1,j)+vortex(i,j+1))/dh^2+b*vortex(i-1,j)+c*vortex(i,j-1)+renold*vortex(i,j)/(dt))/a;elseif(velocityx(i,j)>=0)&&(velocityy(i,j)<0)vortex(i,j)=((vor
7、tex(i+1,j)+vortex(i,j-1))/dh^2+b*vortex(i-1,j)+c*vortex(i,j+1)+renold*vortex(i,j)/(dt))/a;elseif(velocityx(i,j)<0)&&(velocityy(i,j)>=0)vortex(i,j)=((vortex(i-1,j)+vortex(i,j+1))/dh^2+b*vortex(i+1,j)+c*vortex(i,j-1)+renold*vortex(i,j)/(dt))/a;elseif(v
8、elocityx(i,j)<0)&&(velocityy(i,j)<0)vortex(i,j)=((vortex(i-1,j)+vortex(i,j-1))/dh^2+b*vortex(i+1,j)+c*vortex(i,j+1)+renold*vortex(i,j)/(dt))/a;endendendendc1=max(abs(fluidfun));b1=find(c1==max(c1));%%%在第幾列d1=find(abs(fluidfun(:,b1))==max(c1));%%%在第幾行
