基于simple算法的流场模拟计算Word格式文档下载.docx

1、nenenene,pu,（）（）,uudxdyvudxdydxdyudxdy，,，,xyxxx,swswswswnene,u ，udxdySdxdyu,yy,swsw上式积分得到:1FuFuyFuFuxppyDuuDuu,，,，,）-）（）（）（，eewwnnssweeEPwPW2，,，,DuuDuuSSuxy）-）+（，nNPsPSCPP对流项采用一阶迎风格式，扩散项采用中心差分，,yywe（） , （） , FuyDFuyD,wwweeexx,WPPE ,xy,ns（） , （） , FvxDFuyD,nnnsssyx,PNSP1,auauauauauppyb,，,，（）PPWWEESSN

2、NWE,2,aDFaD,，, , WwwEe,11, aDFaDF,，, , ,SssNnn22,aaaaaFFFF,，,，,（）（）,PWESNewns,1,bPP,y（）WE,2,开始 变量的初场是否需要视情况而定） 假设一个速度初场（其它,假设一个压力场，即给定压力猜测值 根据当前的已知量，计算动量离散方程等方程中的系数和常数项 ,0赋值: pp, 步骤1:依次求解动量离散方程 * auauAppb,，,，（）,eenbnbePE *avavAppb,，,，（）, nnnbnbnPN赋值: ,0uv, 步骤2:根据速度求解压力修正方程 pp,0uu,apapapapapb,， PPWWE

3、ESSNN 0* ,buuyvvx,，,（）（）（）（）,wesn 步骤3:对压力和速度进行修正 *ppp,， * （）uuuudpp,，,，,eeeeePE * （）vvvvdpp,，,，,nnnnnPNaab,，,步骤4:求解其他变量的离散方程（视需要进行） ,ppnbnb 否 收敛否, 是 结束 TDMA算法 ,C11,，,DC,2122232 ,，,DC,3233343,，,DC4344454,，,DCnnnnnnn,，11 ,C,nn，11在上式中，假定和是边界上的值，为己知。上式中任一方程都可写成: ,1n，1,，,DCjjjjjjj,，11 除第一及最后一个方程外，其余方程可写为

4、:,C222,，,231DDD222C,333 ,，,342DDD333C444,，453,DDD444,Cnnn,， ,nnn，,11DDDnnn这些方程可通过消元和回代两个过程求解。,C,22+C,313,DD,22,3, ,，,34,22,DD,3333,DD,22,现引入记号:,C222AC,+ ,221DDD222,CC+3323则， ,，,34DADA,332332,CC+3323,AC, 33,DADA,332332,，AC即: ，jjjj1,CC+jjjj,1 ,AC,jj,DADA,jjjjjj,11在边界点，j=1与j=n+1， AC,0,111 AC,0,nnn，111为

5、了求解方程组，首先要对方程组按的形式编排，并明确其中的系数和。从D,Cjjjjj=2起，计算出和，直到j=n。由于在边界位置（n+1）的数值是已知的，因此，C,Ajj可连续计算出。 ,j3、计算结果对比 本文基于simple算法用编写MATLAB程序对整个流场进行计算，另外借助fluent计算流体力学软件数值模拟，通过对比分析计算结果，得出整个流场的流速矢量图和速度云图。0.30数值计算0.25 matlab计算0.200.15u0.100.050.000.00000.00250.00500.00750.01000.0125y 上图为两种计算方法的计算结果，发现两者速度分布趋势一致，编程计算结

6、果数值偏大，原因在于计算采用一阶迎风格式，节点数值趋近于其上游节点值。下图为速度云图，入口段效应的影响导致流动未达到充分发展。MATLAB程序:% dx=5e-4;dy=5e-4;den=998;dyna=1001.6e-6;a=0.01/dy+3;%x节点% b=0.15/dx+1;%y节点% % %边界条件的设置、初始值的设置% % T=0.2*ones（a,b）;U1_old=zeros（a,b）;P1_old=zeros（a,b）;%初值% U1_old（3:a-2,:）=T（3:）;%加速收敛，边界条件% uw_sum（1,b）=0;ue_sum（1,b）=0;de1_sum（1,b

7、）=0;m=1;A（1:a-4,2:b-1）=0.2;%初值为1% e22=1;while e220.01%整个流场迭代% if m4 break end m=m+1;for j=2:1:b-1;if j2 n=1;%迭代次数% e=1;%初次迭代% B=0;0;0.0823455;0.149211;0.196652;0.224194;0.23674; 0.241087;0.242138;0.242214;0.242125;0.242083; 0.242125;0.241087; 0.224195;0.149212;0.0823459;0;U1_old（3:a-2,1）=B（3:a-2,1）;

8、a-2,2）=B（3:while e0.001%迭代，算两列，其他为已知% % %对（i,j）点求uxin% for i=3:a-2 Fe1（i,j）=den*（U1_old（i,j）+U1_old（i,j+1）*dy/2;%对流强度% Fw1（i,j）=den*（U1_old（i,j-1）+U1_old（i,j）*dy/2;ap_u1（i,j）=max（Fe1（i,j）,0）+max（-Fw1（i,j）,0）+4*dyna;%离散方程对应系数% ae_u1（i,j）=dyna+max（-Fe1（i,j）,0）;aw_u1（i,j）=dyna+max（Fw1（i,j）,0）;an_u1（i,j

9、）=dyna;as_u1（i,j）=dyna;A1_u（i,i-1）=-an_u1（i,j）;A1_u（i,i）=ap_u1（i,j）;A1_u（i,i+1）=-as_u1（i,j）;b1_u（i,j）=aw_u1（i,j）*U1_old（i,j-1）+ae_u1（i,j）*U1_old（i,j+1）. +（P1_old（i,j）-P1_old（i,j+1）*dy/2;%（i,j+1）点的源项% A1_u（1,1）=1;%附加点% A1_u（2,2）=1;%边界点% A1_u（a-1,a-1）=1;A1_u（a,a）=1;b1_u（1,j）=0;b1_u（2,j）=0;b1_u（a-1,j）=

10、0;b1_u（a,j）=0;U1_new（:,j）=inv（A1_u）*b1_u（:,j）;%得到（i,j+1）点的速度计算值% %对（i，j）点求p% uw_sum（1,j）=0;ue_sum（1,j）=0;de1_sum（1,j）=dy/ap_u1（2,j）;uw_sum（1,j）=uw_sum（1,j）+U1_old（i,j-1）;ue_sum（1,j）=ue_sum（1,j）+U1_new（i,j）;de1_sum（1,j）=de1_sum（1,j）+dy/ap_u1（i,j）;P1_fix（1,j）=（uw_sum（1,j）-ue_sum（1,j）/de1_sum（1,j）; %对（

11、i，j）点进行速度修正% U1_fix（i,j）=（dy/ap_u1（i,j）*P1_fix（1,j）;%速度修正% a-2,j）=U1_new（3:a-2,j）+U1_fix（3:a-2,j）;P1_old（2:a-1,j）=P1_old（2:a-1,j）+0.4*P1_fix（1,j）;ee=uw_sum（1,j）-ue_sum（1,j）;e=max（max（abs（ee）;%判断的是否合理,% n=n+1;N（1,j）=n;if n1000 A（a-4）*m+1:（a-4）*（m+1）,2:b-1）=U1_old（3:a-2,2:b-1）;%赋值，即m从1开始% E（1:j）=A（a-4）*m+1:j）-A（a-4）*（m-1）+1:（a-4）*m,2:j）;e22=max（max（abs（E）;end%整个流场迭代结束% y=0:10;%结果的输出% u=U1_old（2:2:a-1）;u plot（y,u）