流体流动计算非常复杂,远远超出本教材。最常见的SIMPLE算法编程实现也是比较复杂的,幸运的是MAC算法可以相对容易的实现一些简单算例。Matlab file exchange上一个顶驱方腔流动的例子,使用Matlab计算流体流动,代码如下:
1. clc
2. clearvars
3. %parametersthat can be modified
4. nx=40;%number of points
5. ny=40;
6. xmin=0;xmax=1; %domain dimentions
7. ymin=0;ymax=1;
8. dt=0.0005;
9. nstep=40000;
10. mu=0.00002;
11. un=1;us=0;ve=0;vw=0;% b.c for tangential velocity
12. maxit=850;
13. pit=150;%plot each 'pit' iterations
14. plotvorticity=0;% 1 to plot vorticity 0 to not plot
15. %----------------------------------
16. %--------initializing--------------
17. %----------------------------------
18. time=0.0;
19. dx=(xmax-xmin)/nx;
20. dy=(ymax-ymin)/ny;
21. [X,Y]=meshgrid(xmin:dx:xmax,ymin:dy:ymax);
22. u=zeros(nx+1,ny+2);
23. v=zeros(nx+2,ny+1);
24. p=zeros(nx+2,ny+2);
25. pp=zeros(nx+1,ny+1);
26. ut=zeros(nx+1,ny+2);
27. vt=zeros(nx+2,ny+1);
28. uu=zeros(nx+1,ny+1);
29. vv=zeros(nx+1,ny+1);
30. if plotvorticity==1
31. w=zeros(nx+1,ny+1);
32. end
33.
34. foris=1:nstep
35. %enforcingb.c by interpolation
36. u(1:nx+1,1)=2*us-u(1:nx+1,2);
37. u(1:nx+1,ny+2)=2*un-u(1:nx+1,ny+1);
38. v(1,1:ny+1)=2*vw-v(2,1:ny+1);
39. v(nx+2,1:ny+1)=2*ve-v(nx+1,1:ny+1);
40.
41. i=2:nx;j=2:ny+1; % temporary u-velocity
42. ut(i,j)=u(i,j)+dt*((-0.25)*...
43. (((u(i+1,j)+u(i,j)).^2-(u(i,j)+u(i-1,j)).^2)/dx... %DUUDX
44. +((u(i,j+1)+u(i,j)).*(v(i+1,j)+v(i,j))... %DUVDY 1
45. -(u(i,j)+u(i,j-1)).*(v(i+1,j-1)+v(i,j-1)))/dy)...%DUVDY 2
46. +(mu)*((u(i+1,j)-2*u(i,j)+u(i-1,j))/dx^2... %D2UDX
47. +(u(i,j+1)-2*u(i,j)+u(i,j-1))/dy^2));%D2UDY
48.
49. i=2:nx+1;j=2:ny; % temporary v-velocity
50. vt(i,j)=v(i,j)+dt*((-0.25)*...
51. (((u(i,j+1)+u(i,j)).*(v(i+1,j)+v(i,j))...%DUVDX 1
52. -(u(i-1,j+1)+u(i-1,j)).*(v(i,j)+v(i-1,j)))/dx...%DUVDX 2
53. +((v(i,j+1)+v(i,j)).^2-(v(i,j)+v(i,j-1)).^2)/dy)...%DVVDY
54. +(mu)*((v(i+1,j)-2*v(i,j)+v(i-1,j))/dx^2... %D2VDX
55. +(v(i,j+1)+-2*v(i,j)+v(i,j-1))/dy^2));%D2VDY
56. pt=p;
57. forit=1:maxit % solve for pressure
58. %neumann b.c for pressure
59. pt(1,:)=pt(2,:);
60. pt(nx+2,:)=pt(nx+1,:);
61. pt(:,1)=pt(:,2);
62. pt(:,ny+2)=pt(:,ny+1);
63.
64. i=2:nx+1;j=2:ny+1;
65. pt(i,j)=(0.5/(dx^2+dy^2))*((dy^2)*(pt(i+1,j)+pt(i-1,j))...
66. +(dx^2)*(pt(i,j+1)+pt(i,j-1))...
67. -(dx*dy/dt)*(dy*(ut(i,j)-ut(i-1,j))...
68. +dx*(vt(i,j)-vt(i,j-1))));
69.
70. Er=max(max(pt-p));
71. p=pt;
72. if Er<10^-5
73. break;
74. end
75. end
76. %correct the velocity
77. u(2:nx,2:ny+1)=...
78. ut(2:nx,2:ny+1)-(dt/dx)*(p(3:nx+1,2:ny+1)-p(2:nx,2:ny+1));
79. v(2:nx+1,2:ny)=...
80. vt(2:nx+1,2:ny)-(dt/dy)*(p(2:nx+1,3:ny+1)-p(2:nx+1,2:ny));
81.
82. time=time+dt ; % plotthe results
83. if is==pit*ceil(is/pit)
84.
85. uu(1:nx+1,1:ny+1)=0.5*(u(1:nx+1,2:ny+2)+u(1:nx+1,1:ny+1));
86. vv(1:nx+1,1:ny+1)=0.5*(v(2:nx+2,1:ny+1)+v(1:nx+1,1:ny+1));
87. pp=0.5*(p(1:nx+1,1:ny+1)+p(2:nx+2,2:ny+2));
88. holdoff,h=pcolor(X,Y,flipud(rot90(pp)));
89. set(h,'FaceAlpha',0.3,'EdgeAlpha',0);
90. holdon;
91. quiver(X,Y,flipud(rot90(uu)),flipud(rot90(vv)),2,'color',[00 0]);
92. title(time);
93. if plotvorticity==1
94. w(1:nx+1,1:ny+1)=(u(1:nx+1,2:ny+2)-u(1:nx+1,1:ny+1))/(2*dx)...
95. +(v(1:nx+1,1:ny+1)-v(2:nx+2,1:ny+1))/(2*dy);
96. contour(X,Y,flipud(rot90(w)),10);
97. end
98. axisequal,pause(0.005),drawnow
99. end
100.
101.end
计算结果如下:
此部分内容为选读,仅供参考,有兴趣的读者可以将其移植到HTML5程序。
接下来的课程我们将进入能源与动力工程专业相关的设计类课程,例如换热器的设计、蒸发器、冷凝器、汽包锅炉水位控制等等。
本文分享自 传输过程数值模拟学习笔记 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!