Anderson英文版计算流体力，SIMPLE 算法计算流场程序分享

scott198510

1. 
   
2. % A program of incompressible Couette flow in 2D
   
3. % based on pressure correction (SIMPLE) algorithm.
   
4. % Written by scott , Oct, 2013
   
5. 
   
6. % Based on John D. Anderson, JR., section 9.4, Chapter 9, In
   
7. % Computational Fluid Dynamics, The Basics with Applications,McGraw-Hill, 2002, 4
   
8. clc;clear all; close all;
   
9. nx=21;%x方向网格数
   
10. ny=11;%y方向网格数
    
11. nstop=1000;%迭代步数
    
12. x=zeros(nx,1);%矩阵预分配
    
13. y=zeros(ny,1);%矩阵预分配
    
14. p=zeros(nx,ny);%矩阵预分配
    
15. u=zeros(nx+1,ny);%矩阵预分配
    
16. v=zeros(nx+2,ny+1);%矩阵预分配
    
17. 
    
18. 
    
19. ru=zeros(nx+1,ny);%矩阵预分配
    
20. rv=zeros(nx+2,ny+1);%矩阵预分配
    
21. ux=zeros(nx+1,ny);%矩阵预分配
    
22. vx=zeros(nx+2,ny+1);%矩阵预分配
    
23. 
    
24. 
    
25. px=zeros(nx,ny);%矩阵预分配
    
26. rux=zeros(nx+1,ny);%矩阵预分配
    
27. rvx=zeros(nx+2,ny+1);%矩阵预分配
    
28. d=zeros(nx,ny);%矩阵预分配
    
29. 
    
30. ppm=zeros(nx,ny);%矩阵预分配
    
31. ppb=zeros(nx,ny);%矩阵预分配
    
32. pp=zeros(nx,ny);%矩阵预分配
    
33. cri=zeros(nx,ny);%矩阵预分配
    
34. % initial parameter
    
35. xl=0.5; %x方向长度 单位英尺
    
36. yl=0.01;%y方向长度 单位英尺
    
37. dx=xl/(nx-1); %x方向单位长度
    
38. dy=yl/(ny-1);%y方向单位长度
    
39. dt=0.001;%单位时间
    
40. re=63.6;%雷诺数
    
41. r=0.002377;%密度ρ
    
42. vis=yl*1*r/re;%粘性系数
    
43. fct1=0.5;
    
44. fct2=0.1;
    
45. for i=1:nx
    
46. x(i)=(i-1)*xl/(nx-1);%x方向区域网格化
    
47. end
    
48. for j=1:ny
    
49. y(j)=(j-1)*yl/(ny-1);%x方向区域网格化
    
50. end
    
51. for i=1:nx+1 %u迭代区域 x:1-22, y:1-11
    
52. for j=1:ny
    
53. ux(i,j)=0.0;%u*初始化
    
54. u(i,j)=0.0; %u初始化
    
55. if (j==ny) %u的上边界位置特定值
    
56. u(i,j)=1.0; %u=1
    
57. rux(i,j)=r*ux(i,j);%ρu*
    
58. end
    
59. end
    
60. end
    
61. for i=1:nx+2 %u迭代区域 x:1-23, y:1-12
    
62. for j=1:ny+1
    
63. vx(i,j)=0.0; %v*初始化
    
64. v(i,j)=0.0; %v初始化
    
65. if (i==15&&j==5) %设置特殊点（15,5)速度
    
66. vx(i,j)=0.5;
    
67. rvx(i,j)=r*vx(i,j);%ρv*
    
68. end
    
69. end
    
70. end
    
71. for i=1:nx %压力场p迭代区域 x:1-21, y:1-11
    
72. for j=1:ny
    
73. px(i,j)=0.0; %p*初始化
    
74. p(i,j)=0.0; %p初始化
    
75. pp(i,j)=0.0; % p'初始化
    
76. end
    
77. end
    
78. 
    
79. for nstep=1:nstop; %迭代计算部分
    
80. 
    
81. % Pressure correction (SIMPLE) algorithm
    
82. % step 1. solve at n+1 step ： ρu* and ρv*
    
83. % 1. for ru and u
    
84. 
    
85. for i=1:nx-1 %以（5,3）作为中心点 即（i,j)对应（5,3）
    
86. for j=2:ny-1
    
87. vb=0.5*(v(i+1,j+1)+v(i+2,j+1));%v平
    
88. vbb=0.5*(v(i+1,j)+v(i+2,j));%v2平
    
89. ruu53=r*u(i+2,j)*u(i+2,j);%(ρu^2)53
    
90. ruu33=r*u(i,j)*u(i,j);%(ρu^2)33
    
91. ruv44=r*u(i+1,j+1)*vb; %(ρuvb)44
    
92. ruv42=r*u(i+1,j-1)*vbb;%(ρuvbb)42
    
93. uu1=(u(i+2,j)-2*u(i+1,j)+u(i,j))/dx/dx;%u(i+2,j):u53,u(i+1,j):u43，u(i,j):u33
    
94. uu2=(u(i+1,j+1)-2*u(i+1,j)+u(i+1,j-1))/dy/dy;%u(i+1,j+1):u44,u(i+1,j):u43，u(i+1,j-1):u42
    
95. ax=-((ruu53-ruu33)/2/dx+(ruv44-ruv42)/2/dy)+vis*(uu1+uu2);%A*整项
    
96. rux(i+1,j)=rux(i+1,j)+ax*dt-dt/dx*(px(i+1,j)-px(i,j));%（ρu*)43 in n+1 step 对应式子（9.58）
    
97. ux(i+1,j)=rux(i+1,j)/r; %（u*)43 in n+1 step
    
98. end
    
99. end
    
100. 
     
101. for j=2:ny-1
     
102. ux(1,j)=ux(2,j);%上边界赋值 u1j*=u2j*
     
103. ux(nx+1,j)=ux(nx,j);%下边界赋值 u22j*=u21j*
     
104. end
     
105. 
     
106. % 2. for rv and v
     
107. 
     
108. for i=1:nx %以（3,3）作为中心点 即（i,j)对应（3,3）
     
109. for j=1:ny-1
     
110. ub=0.5*(u(i+1,j)+u(i+1,j+1));%u平
     
111. ubb=0.5*(u(i,j)+u(i,j+1)); %u2平
     
112. rvv45=r*v(i+1,j+2)*v(i+1,j+2); %(ρv^2)45
     
113. rvv43=r*v(i+1,j)*v(i+1,j); %(ρv^2)43
     
114. rvu54=r*v(i+2,j+1)*ub; %(ρvub)54
     
115. rvu34=r*v(i,j+1)*ubb; %(ρvubb)34
     
116. vv1=(v(i+2,j+1)-2*v(i+1,j+1)+v(i,j+1))/dx/dx;%u(i+2,j+1):v54,u(i+1,j+1):v44，u(i,j+1):v34
     
117. vv2=(v(i+1,j+2)-2*v(i+1,j+1)+v(i+1,j))/dy/dy;%u(i+1,j+2):v45,u(i+1,j+1):v44，v(i,j+1):v43
     
118. bx=-((rvv45-rvv43)/2/dy+(rvu54-rvu34)/2/dx)+vis*(vv1+vv2);%B*整项
     
119. rvx(i+1,j+1)=rvx(i+1,j+1)+bx*dt-dt/dy*(px(i,j+1)-px(i,j));%（ρv*)44 in n+1 step 对应式子（9.59）
     
120. vx(i+1,j+1)=rvx(i+1,j+1)/r; %（v*)44 in n+1 step
     
121. end
     
122. end
     
123. 
     
124. for j=2:ny
     
125. vx(nx+2,j)=vx(nx+1,j);%下边界赋值 v23j*=v22j*
     
126. end
     
127. 
     
128. % step 2. solve p' from pressure correction formula
     
129. % by relaxation method (inner iteration)
     
130. a=2*(dt/dx/dx+dt/dy/dy);% 对应（6.104）式子a值
     
131. b=-dt/dx/dx; % 对应（6.104）式子b值
     
132. c=-dt/dy/dy;% 对应（6.104）式子c值
     
133. tiny2=0; %设定值
     
134. for i=2:nx-1
     
135. for j=2:ny-1
     
136. pp(i,j)=0;
     
137. d(i,j)=(rux(i+1,j)-rux(i,j))/dx+(rvx(i+1,j+1)-rvx(i+1,j))/dy;% 对应（6.104）式子d值
     
138. end
     
139. end
     
140. for nii=1:1000
     
141. tiny1=0;
     
142. for i=2:nx-1
     
143. for j=2:ny-1 %再以（3,3）作为中心点 即（i,j)对应（3,3）
     
144. ppb(i,j)=-1/a*(b*pp(i+1,j)+b*pp(i-1,j)+c*pp(i,j+1)+c*pp(i,j-1)+d(i,j));% 式（9.60) p33'
     
145. ppm(i,j)=pp(i,j)+fct1*(ppb(i,j)-pp(i,j));%迭代后的pp
     
146. cri(i,j)=abs(ppm(i,j)-pp(i,j))/abs(ppm(i,j));%计算增量
     
147. pp(i,j)=ppm(i,j);%迭代后赋值
     
148. if (cri(i,j)>=tiny1)
     
149. tiny1=cri(i,j);
     
150. end
     
151. if (abs(d(i,j))>=tiny2)
     
152. tiny2=d(i,j);
     
153. end
     
154. end
     
155. end
     
156. if (tiny1<=0.0001); break,end %小于设定值时，即认为已经收敛
     
157. end
     
158. % step 3. calculate p at n+1 step
     
159. for i=2:nx-1
     
160. for j=2:ny-1
     
161. p(i,j)=px(i,j)+fct2*pp(i,j);%对应 式子（6.106)
     
162. end
     
163. end
     
164. % step 4. calculate u and v at n+1 step
     
165. for i=1:nx-1
     
166. for j=1:ny
     
167. u(i+1,j)=ux(i+1,j)-dt/dx/r*(pp(i+1,j)-pp(i,j)); % u更新
     
168. end
     
169. end
     
170. for i=1:nx
     
171. for j=1:ny-1
     
172. v(i+1,j+1)=vx(i+1,j+1)-dt/dy/r*(pp(i,j+1)-pp(i,j));% v更新
     
173. end
     
174. end
     
175. % update boundary conditions
     
176. for i=1:nx
     
177. pp(i,ny)=0;
     
178. pp(i,1)=0;
     
179. end
     
180. for i=1:nx+1
     
181. u(i,ny)=1;
     
182. u(i,1)=0;
     
183. end
     
184. for i=1:nx+2
     
185. v(i,ny+1)=0;
     
186. v(i,1)=0;
     
187. end
     
188. for j=1:ny
     
189. pp(1,j)=0;
     
190. pp(nx,j)=0;
     
191. end
     
192. for j=1:ny
     
193. u(1,j)=u(2,j);
     
194. u(nx+1,j)=u(nx,j);
     
195. end
     
196. for j=1:ny+1
     
197. v(1,j)=0;
     
198. v(nx+2,j)=v(nx+1,j);
     
199. end
     
200. 
     
201. % step 5. check the primary iteration
     
202. for i=2:nx-1
     
203. for j=2:ny-1
     
204. px(i,j)=p(i,j);
     
205. end
     
206. end
     
207. end
     
208. figure,surf(p) % pressure distribution in cavity
     
209. figure,contour(p,15) % pressure distribution in cavity
     
210. figure,surf(u) % u component of velocity
     
211. figure,contour(u,15) %u component of velocity
     
212. figure,surf(v) %v component of velocity
     
213. figure,contour(v,15)%v component of velocity</P>
     
214. set(gcf,'color','w');
     
215. 
     
216. 
     

复制代码

SIMPLE 算法在不可压缩流场中应用广泛，看到有不少做这方面的开发的，附上 John D. Anderson  英文版计算流体力学，
section 9.4, Chapter 9，采用SIMPLE 算法的算例的程序，程序采用matlab语言编写，属原创。

小七工作室

好资料 谢谢分享

pank

不错，好好学习一下，一直先做这个开发

scott198510

pank 发表于 2013-10-19 11:24
不错，好好学习一下，一直先做这个开发

1.       
   
2.      clc;clear all;close all;
   
3.       u=zeros(32,1);uf=zeros(33,1);
   
4.       a=zeros(33,1);p=zeros(32,1);
   
5.       pp=zeros(32,1);aa=zeros(32,1);
   
6.       su=zeros(32,1); aup=zeros(33,1);
   
7.       aue=zeros(33,1);auw=zeros(33,1);
   
8.       ap=zeros(32,1);ae=zeros(32,1);
   
9.       aw=zeros(32,1);s=zeros(32,1);
   
10.       
    
11.      dx=0.01;
    
12.       for i=1:33
    
13.                 aaa=0.1/16^2;
    
14.                 bb=0.1;
    
15.                 a(i)=aaa*((i-17).^2)+bb;
    
16.       end
    
17.                 for i=1:32
    
18.                    aa(i)=0.5*(a(i)+a(i+1));
    
19.         end
    
20.   % 密度参数设置
    
21.                 den=1000;
    
22.   % 边界条件设置
    
23.                 u(1)=10;
    
24.                 u(32)=10;
    
25.                 viscos=0.001;
    
26.    %面对速度的边界条件
    
27.                 uf(1)=10;
    
28.                 uf(33)=10;
    
29.    % 初始条件设置
    
30.           for i=2:32
    
31.              u(i)=10;
    
32.              uf(i)=10;
    
33.              p(i)=0;
    
34.       end
    
35.    % 循环 (迭代开始)
    
36. for N=1:1000
    
37.   % 二维网格解动量守恒方程
    
38.           for i=2:32
    
39.           aue(i)=max(abs(den*u(i)*aa(i)/2),viscos*aa(i)/dx)-den*u(i)*aa(i)/2;
    
40.           auw(i)=max(abs(den*u(i-1)*aa(i-1)/2),viscos*aa(i-1)/dx) +den*u(i-1)*aa(i-1)/2;
    
41.           su(i)=a(i)*(p(i-1)-p(i));
    
42.           aup(i)=auw(i)+aue(i);
    
43.       end
    
44. %==== write(6,*)aue(10),auw(10),aup(10),su(10),' aue,auw,aup,su'
    
45. %  small iterations
    
46.           for i=2:32
    
47.             uf(i)=(aue(i)*uf(i+1)+auw(i)*uf(i-1)+su(i))/aup(i);
    
48.           end
    
49. % extrapolation
    
50.           aue(1)=aue(2);
    
51.           auw(1)=auw(2);
    
52.           aue(33)=aue(32);
    
53.           auw(33)=auw(32);
    
54. %  主网格解压力修正方程
    
55.           for i=2:31
    
56.           aw(i)=den*a(i).^2/aup(i);
    
57.           ae(i)=den*a(i+1).^2/aup(i+1);
    
58.           ap(i)=aw(i)+ae(i);
    
59.           s(i)=den*a(i)*uf(i)-den*a(i+1)*uf(i+1);
    
60.       end
    
61. % set pp=0
    
62.           for i=1:32
    
63.              pp(i)=0;
    
64.       end
    
65. % small iterations
    
66.           for its=1:30
    
67.           for i=2:31
    
68.           pp(i)=(ae(i)*pp(i+1)+aw(i)*pp(i-1)+s(i))/ap(i);
    
69.           end
    
70.           end  
    
71. % correct velocity
    
72.           for i=2:32
    
73.           uf(i)=uf(i)+(pp(i-1)-pp(i))*a(i)/aup(i);
    
74.           end
    
75. % interpolation
    
76.           for i=1:32
    
77.           u(i)=0.5*(uf(i+1)+uf(i));
    
78.           end
    
79. % correct pressure
    
80.           for i=2:31
    
81.             p(i)=p(i)+0.1*pp(i);
    
82.       end         
    
83. % extrapolation
    
84.           p(1)=2.*p(2)-p(3);
    
85.           p(32)=2.*p(31)-p(30);
    
86. end

复制代码

再发一个。

pank

非常感谢scott198510的无私共享！

ybli

感谢楼主，先收藏了，细细拜读