《中国科学院大学 计算流体力学 作业 11.1》由会员分享，可在线阅读，更多相关《中国科学院大学 计算流体力学 作业 11.1（8页珍藏版）》请在金锄头文库上搜索。
1、U计算流体力学作业11.1习题11.1求解方腔问题问题描述：如图示边K为L的方腔，匕表血流体以常速度U运 动,求解里面的流场（假设流动定常）。考虑 Re = =100,400,1000 二种情况要求：数值方法不限（人匚技缩性方法、投影法、制收流函数方法及SIMPLE方 法均可）；空间离散采用差分法,建议采用较高阶精度的方法o绘制出定常解的流线图。请详细写明方程及公式的推导过程及计算流程，切勿只上交计算结果。首先，进行公式推导。涓第加伊力号作划I岑祐依刈忖）羽小功残十为二。东十嗫+V为二-节十古内 。 我十联叶埼二啮十古加奇-聂依5二奇一会，有攀十啜二右血 Q加的根，u二喙v二一头,有便无靴碱-
2、痴淞缸刚心格闻、柄抨判谭晶（G版:不富依：球网5恼网棒衣 _，3W _呜（做）-6奶匠岫广诩加 啊（W吐一峭3啊cW*D.X/hj 弧64K（I ?吆）；- Yij（W珂dWirj十,电彳十刈卬j） +均（助坎十那么尸刈叫）弟枝尿：物中心为今以勾川十助人一刈可以勾川十助人一刈可二 W%，j 十助fj-2WqA犬加跖-就无无吊漏翱：御务克多白外十7厂2%上H十知T -W可4X 十折二外的可悔遗：?毕炎河,JJ42.=用 一）七4七令 Ax = Ay = A/? o次边界涡量采用中心差分格式。迭代采用Jacobi迭代形式。例如：ij=R-+ +k,H 一（八力）22,）当速度u、V与上一步的最大误
3、差小于10-5时结束循环。边界条件速度边界条件：上壁面u=l,v=0;其他壁面u=v=0；流函数的边界条件：h=o涡量的边界条件：可按照定义由速度给出I du 切尸。=而Ax解题思路:设置初始边界u,v得出边界一一进行空间迭代一一进行时间推 进一一得到内部一一时间推进一一反求得到流函数反求得到 最终速度u, Vo注：在设置节点的时候，参考同学思路，令1、3和104、105为固壁 边界处，3-103为所求方腔及其内部。流线图结果:(1) Re=100(2) Re=4000.80.70.60.50.40.30.20.10.20.40.60.81(3) Re=1000由图可知，上壁面有一个向右的速度
4、，带动下层流体运动，形成一个涡。旋向为顺时针。随着Re数的增大，方腔的涡中心不断向下 移动。除此之外，在边角处逐渐形成角涡，旋向与中心涡相反，为逆 时针，且随Re增加而逐渐长大。附上源代码：clc;clear;dh=0.01;%空间步长dt=0.001;%时间步长R=1000; %修改 Re=100,400, 1000之初始化变量u=zeros (200,200) ; %x 方向速度 uv=zeros (200,200) ; %y 方向速度 vu_nxt=z;ros (200,200) ; %u 迭代后结果v_nxt=zros (200,200) ; %v 迭代后结果w=ones (200,2
5、00) ; %涡量streamf=zeros (200,200) ; %流函数%定义误差函数用于存放迭代前后U、V的误差值的最大值error_um=l;error vm=l;%迭代求解k=l;超迭代次数%上边界速度为1u(103:105,3:103)=1;%当最大误差小于10人(-5)时结束循环while error_um10A(-5)|error_vm10A(-5)茗求边界的涡量for i=3:102w(i, 103) = (u(i + 1,103)-u(i,103)/dh-(v (i,103)-v (if102) /dh;w(i,3) = (u(i,3) - u(i-l,3) ) /dh-
6、 (v (i, 4 ) -v (i, 3) ) /dh; endfor j=3:102w(3, j) = (u(4,j)-u(3,j)/dh-(v(3,j+1)-v (3,j)/dh;w(103,j)=(u(103,j)-u(102,j)/dh-(v(103zj)-v(103,j-l) /dh;end%求次边界涡量(采用中心差分)for j=4:102w(4,j)=l/4*(w(5,j)+w(3,j)+w(4,j+l)+w(4,j-l)-R*dh/8* (u (4, j )夫(w (5, j ) -w (3, j ) )+v(4,j)*(w(4,j+l)-w(4,j-l);w(102,j)=1
7、/4* (w(103,j)+w(101,j)+w(102,j+1)+w(102, j-l) )-Re*dh/8*(u(102,j)*(w(103,j)-w(101,j)+v(102,j)*(w(1 02,j+l)-w(102,j-1);endfor i=4:102w (i,4)=1/4*(w(i + 1,2)+w(i-l,4)+w(i,5)+w(i,3)-R*dh/8* (u(i,4) * (w(i+l,4)-w(i-l,4) )+v(i,4) * (w (i, 5)-w (i, 3);w(i, 102) =1/4*(w(i + 1,102)+w(i-lz102)+w(i,103)+w(iA1
8、01) )-R-dh/8* (u (i, 102) * (w (i + 1,102) -w(i-l,102) )+v(i,102)* (w(i,103)-w(i,101);end%对内部涡量进行迭代for i=4:1:102for j=4:1:102%速度正负U_pos(i,j) = (u(i,j)+abs(u(i, j ) ) ) /2;u_neg(i,j) = (u(i,j)-abs(u(iz j)/2;v_pos(i,j)=(v(i,j)+abs(v(i,j)/2;v_neg(i,j) = (v(i,j)-abs(v(i, j ) ) ) /2;用对流项采用三阶迎风格式差分diffwx_
9、pos(i,j)=(w(i,j-2)-6*w(i,j-l)+3*w(i,j)+2*w(i,j +1)/(6*dh);diffwx_neg(i,j)=-(w(i,j+2)-6*w(i,j+l)+3*w(i,j)+2*w(i, j-1)八6夫dh);diffwy_pos(i,j)=(w(i-2,j)-6*w(i-l,j)+3*w(i,j)+2*w(i+l ,j)/(6*dh);diffwy_ng(i,j)=-(w(i+2,j)-6*w(i+l,j)+3*w(i,j)+2*w(i- 1, j) )/?6*dh);%粘性项采用中心差分diff2wx(i,j)=(w(i,j+l)+w(i,j-l)-2*
10、w(i,j)/(dh*dh);diff2wy(i,j)=(w(i+l,j)+w(i-l,j)-2*w(i,j)/(dh*dh);%迭代涡量ww(i,j)=dt*(diff2wx(i,j)+diff2wy(iA j)/Re-(u_pos(i,j) *diffwx_neg(i,j)+u_neg(i,j)*diffwx_pos(i,j)-(v_pos(i ,j)*diffwy_neg(i,j)+v_neg(i,j)*diffwy_pos(i,j)+w(i, j)；%反解流函数、迭代后速度 streamf(i,j) = (streamf(i + lz j)+streamf(i-1,j)+streamf
11、(i ,j+1)+streamf(iA j-1)-w(i,j)*dhA2)/4;u_next(i,j)=(streamf(i+1,j)-stramf(i-1,j)/2/dh;v_next(i,j)=-(streamf(i,j+1)-streamf(iz j-1)/2/dh; ndndu_next(103:105z 3:103)=1;for i=3:1:103for j=3:1:103国定义误差error_u(i-2z j-2)=u_next(i,j)-u (i,j);error_v(i-2,j-2)=v_next(i,j)-v(i,j);endnd告迭代后速度u_next替换原速度u进行下一步迭代u=u_next;v=v_next;e r r o r_um=ma x(max(rror_u);e r r o r_vm=ma x(max(error_v);k=k+l;end%绘制流线x=0:0.01:1;y=0:0.01:1;figurestreamslice(x,y,u(3:103,3:103),v(3:103A 3:103);axis equal;axis ( 0 1 0 1);box on;
