SIMPLE算法处理流场分享

scott198510

1. clear all;clc;close all;
   
2. %initialization
   
3. % Guess a temporary pressure and velocity fields
   
4. p_star = zeros(41);
   
5. u_star = zeros(41);u_star(1,:) = 5;u_star(2,:) = 5;
   
6. v_star = zeros(41);
   
7. 
   
8. 
   
9. p_prime = zeros(41);
   
10. p_prime2= zeros(41);
    
11. u_prime = zeros(41);
    
12. v_prime = zeros(41);
    
13. 
    
14. 
    
15. p = p_star+p_prime;
    
16. u = u_star+u_prime;
    
17. v = v_star+v_prime;
    
18. 
    
19. alpha = 0.4;
    
20. dt = 0.001;
    
21. dx = 0.01;
    
22. dy = 0.01;
    
23. ro = 1000;
    
24. MU = 0.0025;
    
25. 
    
26. for itr = 1:2000
    
27. for i = 2:40
    
28. for j = 2:40
    
29. A_star = MU*((u(i+1,j)-2*u(i,j)+u(i-1,j))/dx^2 + (u(i,j+1)-2*u(i,j)+u(i,j-1))/dy^2);
    
30. B_star = MU*((v(i+1,j)-2*v(i,j)+v(i-1,j))/dx^2 + (v(i,j+1)-2*v(i,j)+v(i,j-1))/dy^2);
    
31. u_star(i,j) = u_star(i,j) + A_star*dt - dt/dx*(p_star(i-1,j) - p_star(i,j))/ro;
    
32. v_star(i,j) = v_star(i,j) + B_star*dt - dt/dx*(p_star(i,j-1) - p_star(i,j))/ro;
    
33. end
    
34. end
    
35. u_star(41,:) = 0;u_star(40,:) = 0;u_star(:,1) = 0;u_star(:,2) = 0;u_star(:,41) = 0;u_star(:,40) = 0;u_star(2,:) = 5;u_star(1,:) = 5;
    
36. v_star(41,:) = 0;v_star(40,:) = 0;v_star(:,1) = 0;v_star(:,2) = 0;v_star(:,41) = 0;v_star(:,40) = 0;v_star(2,:) = 0;v_star(1,:) = 0;
    
37. 
    
38. u = u_star;
    
39. v = v_star;
    
40. 
    
41. for itrate = 1:1000
    
42. for i = 2:40
    
43. for j = 2:40
    
44. p_prime2(i,j) = 0.25*(p_prime(i-1,j)+p_prime(i+1,j)+p_prime(i,j-1)+p_prime(i,j+1)) + dx*(u_star(i,j)-u_star(i+1,j)+v_star(i,j)-v_star(i,j+1));
    
45. end
    
46. end
    
47. mx = max(max(abs(p_prime-p_prime2)));
    
48. p_prime = p_prime2;
    
49. if mx<0.001, break, end
    
50. end
    
51. 
    
52. p = p_star + alpha*p_prime;
    
53. p_star = p;
    
54. 
    
55. end
    
56. 
    
57. figure,surf(p) % pressure distribution in cavity
    
58. figure,contour(p,20) % pressure distribution in cavity
    
59. figure,surf(u) % u component of velocity
    
60. figure,contour(u,20) %u component of velocity
    
61. figure,surf(v) %v component of velocity
    
62. figure,contour(v,20)%v component of velocity
    

复制代码

鉴于有不少学习SIMPLE算法处理流场（耦合方程中的速度场）的迭代问题，附上一个小程序，便于新手学习。

采用SIMPLE算.法实施速度分量和压力场耦合方程的离散求解时，计算步骤如下：
（1） 假定一个速度分布，记为u0,v0,w0 ，以此计算动量离散方程中的系数及常数项；
（2） 假设一个压力场p0 ；
（3） 依次求解动量方程，得 u1,v1,w1；
（4） 对压力加以修正，得p1 ；
（5） 根据p1 改进速度值；
（6） 利用改进后的速度场求解那些通过源项物性等与速度场耦合的 Φ变量，如果Φ 变量并不影响流场，则应在速度场收敛后再求解；
（7） 利用改进后的速度场重新计算动量离散方程的系数，并利用改进后的压力场作为下一层次迭代计算的初值。重复上述步骤，直到获得收敛的解。

表达式如图片中所示，前两式为动量守恒方程，第三式为质量守恒方程。

liuyalong008

好久不见scott198510，其实matlab做流场是有局限性的，很多情况下进行差分求解动量方程时用的是迎风格式，下一个step是受上一个step的影响的，matlab所擅长的向量化运算的优势是无从体现啊，这也就是现在好多流场模拟都是fortran等其他语言来完成的原因。不过simple算法确实是一个经典的案例

scott198510

liuyalong008 发表于 2013-9-18 14:04
好久不见scott198510，其实matlab做流场是有局限性的，很多情况下进行差分求解动量方程时用的是迎风格式， ...

表达式如上。

scott198510

1. clear all;clc;close all;
   
2. %initialization
   
3. nx=100;
   
4. ny=30;
   
5. %速度场
   
6. u=zeros(ny+2,nx+1,3);
   
7. v=zeros(ny+1,nx+2,3);
   
8. %压力场
   
9. p=zeros(ny,nx,3);
   
10. %===========动量守恒方程=============
    
11. aeu=zeros(ny+2,nx+1);
    
12. awu=zeros(ny+2,nx+1);
    
13. asu=zeros(ny+2,nx+1);
    
14. anu=zeros(ny+2,nx+1);
    
15. apu=zeros(ny+2,nx+1);
    
16. aev=zeros(ny+1,nx+2);
    
17. awv=zeros(ny+1,nx+2);
    
18. asv=zeros(ny+1,nx+2);
    
19. anv=zeros(ny+1,nx+2);
    
20. apv=zeros(ny+1,nx+2);
    
21. %============压力修正================
    
22. pce=zeros(ny,nx);
    
23. pcw=zeros(ny,nx);
    
24. pcs=zeros(ny,nx);
    
25. pcn=zeros(ny,nx);
    
26. aep=zeros(ny,nx);
    
27. awp=zeros(ny,nx);
    
28. asp=zeros(ny,nx);
    
29. anp=zeros(ny,nx);
    
30. app=zeros(ny,nx);
    
31. aepc=zeros(ny,nx);
    
32. awpc=zeros(ny,nx);
    
33. anpc=zeros(ny,nx);
    
34. aspc=zeros(ny,nx);
    
35. appc=zeros(ny,nx);
    
36. pc=zeros(ny+2,nx+2,3);
    
37. %=========质量不守恒，需要通过压力修正===========
    
38. mp=zeros(ny,nx);
    
39. errp=zeros(ny,nx);
    
40. raw=1000;
    
41. dx=0.001/ny;
    
42. dy=0.001/ny;
    
43. d=0.001;
    
44. u(2:ny+1,:,:)=0.01;
    
45. for q=1:5000
    
46. %===========动量守恒方程=============
    
47. for ll=1:10
    
48. for i=2:ny+1
    
49. for j=2:nx
    
50. if i==2
    
51. aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
    
52. awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
    
53. asu(i,j)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*v(i-1,j+1,1)+0.25*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
    
54. anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
    
55. if j==nx
    
56. aeu(i,j)=0;
    
57. end
    
58. apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
    
59. u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy)/apu(i,j);
    
60. else if i==ny+1
    
61. aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
    
62. awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
    
63. asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
    
64. anu(i,j)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*v(i,j+1,1)+0.25*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
    
65. if j==nx
    
66. aeu(i,j)=0;
    
67. end
    
68. apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
    
69. u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy)/apu(i,j);
    
70. else
    
71. aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
    
72. awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
    
73. asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
    
74. anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
    
75. if j==nx
    
76. aeu(i,j)=0;
    
77. end
    
78. apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
    
79. u(i,j,1)=(aeu(i,j)*u(i,j+1,1)+awu(i,j)*u(i,j-1,1)+anu(i,j)*u(i+1,j,1)+asu(i,j)*u(i-1,j,1)+(p(i-1,j-1,1)-p(i-1,j,1))*dy)/apu(i,j);
    
80. end
    
81. end
    
82. end
    
83. end
    
84. u(:,nx,1)=u(:,nx,1)*(sum(u(:,1,1))/sum(u(:,nx,1)));
    
85. u(:,nx+1,1)=u(:,nx,1);
    
86. for i=2:ny
    
87. for j=2:nx+1
    
88. if j==2
    
89. aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
    
90. awv(i,j)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j-1,1)+0.25*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
    
91. asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
    
92. anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
    
93. if j==nx+1
    
94. aev(i,j)=0;
    
95. end
    
96. apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
    
97. v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy)/apv(i,j);
    
98. else if j==nx+1
    
99. aev(i,j)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j,1)+0.25*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
    
100. awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
     
101. asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
     
102. anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
     
103. if j==nx+1
     
104. aev(i,j)=0;
     
105. end
     
106. apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
     
107. v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+(p(i-1,j-1,1)-p(i,j-1,1))*dy)/apv(i,j);
     
108. else
     
109. aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
     
110. awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
     
111. asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
     
112. anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
     
113. if j==nx+1
     
114. aev(i,j)=0;
     
115. end
     
116. apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
     
117. v(i,j,1)=(aev(i,j)*v(i,j+1,1)+awv(i,j)*v(i,j-1,1)+anv(i,j)*v(i+1,j,1)+asv(i,j)*v(i-1,j,1)+...
     
118. (p(i-1,j-1,1)-p(i,j-1,1))*dy)/apv(i,j);
     
119. end
     
120. end
     
121. end
     
122. end
     
123. end
     
124. %================压力场修正============
     
125. for i=2:ny+1
     
126. for j=2:nx
     
127. if i==2
     
128. aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
     
129. awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
     
130. asu(i,j)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*v(i-1,j+1,1)+0.25*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
     
131. anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
     
132. if j==nx
     
133. aeu(i,j)=0;
     
134. end
     
135. apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
     
136. else if i==ny+1
     
137. aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
     
138. awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
     
139. asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
     
140. anu(i,j)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*v(i,j+1,1)+0.25*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
     
141. if j==nx
     
142. aeu(i,j)=0;
     
143. end
     
144. apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
     
145. else
     
146. aeu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j+1,1)+0.5*u(i,j,1)));
     
147. awu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i,j,1)));
     
148. asu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j+1,1)+0.5*v(i-1,j,1)));
     
149. anu(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i,j+1,1)+0.5*v(i,j,1)));
     
150. if j==nx
     
151. aeu(i,j)=0;
     
152. end
     
153. apu(i,j)=aeu(i,j)+awu(i,j)+asu(i,j)+anu(i,j);
     
154. end
     
155. 
     
156. end
     
157. end
     
158. end
     
159. for i=2:ny
     
160. for j=2:nx+1
     
161. if j==2
     
162. aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
     
163. awv(i,j)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j-1,1)+0.25*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
     
164. asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
     
165. anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
     
166. apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
     
167. else if j==nx+1
     
168. aev(i,j)=2*d*max(0,(1-0.1*abs(raw*dy*(0.25*u(i,j,1)+0.25*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
     
169. awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
     
170. asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
     
171. anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
     
172. apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
     
173. else
     
174. aev(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1))/d))^5)+max(0,-raw*dy*(0.5*u(i,j,1)+0.5*u(i+1,j,1)));
     
175. awv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1))/d))^5)+max(0,raw*dy*(0.5*u(i,j-1,1)+0.5*u(i+1,j-1,1)));
     
176. asv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,raw*dy*(0.5*v(i-1,j,1)+0.5*v(i,j,1)));
     
177. anv(i,j)=d*max(0,(1-0.1*abs(raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1))/d))^5)+max(0,-raw*dy*(0.5*v(i+1,j,1)+0.5*v(i,j,1)));
     
178. apv(i,j)=aev(i,j)+awv(i,j)+asv(i,j)+anv(i,j);
     
179. end
     
180. end
     
181. end
     
182. end
     
183. 
     
184. for i=1:ny
     
185. for j=2:nx-1
     
186. aepc(i,j)=raw*dy^2/apu(i+1,j+1);
     
187. awpc(i,j)=raw*dy^2/apu(i+1,j);
     
188. end
     
189. end
     
190. for i=2:ny-1
     
191. for j=1:nx
     
192. anpc(i,j)=raw*dx^2/apv(i+1,j+1);
     
193. aspc(i,j)=raw*dx^2/apv(i,j+1);
     
194. end
     
195. end
     
196. for i=1:ny
     
197. awpc(i,1)=0;
     
198. aepc(i,1)=raw*dy^2/apu(i+1,2);
     
199. awpc(i,nx)=raw*dy^2/apu(i+1,nx);
     
200. aepc(i,nx)=0;
     
201. end
     
202. for j=1:nx
     
203. anpc(1,j)=raw*dx^2/apv(2,j+1);
     
204. aspc(1,j)=0;
     
205. anpc(ny,j)=0;
     
206. aspc(ny,j)=raw*dx^2/apv(ny,j+1);
     
207. end
     
208. for i=1:ny
     
209. for j=1:nx
     
210. appc(i,j)=aepc(i,j)+awpc(i,j)+anpc(i,j)+aspc(i,j);
     
211. end
     
212. end
     
213. 
     
214. for i=1:ny
     
215. for j=1:nx
     
216. mp(i,j)=raw*dx*((u(i+1,j+1,1)-u(i+1,j,1))+(v(i+1,j+1,1)-v(i,j+1,1)));
     
217. end
     
218. end
     
219. 
     
220. for ll=1:1000
     
221. for i=2:ny+1
     
222. for j=2:nx+1
     
223. pc(i,j,1)=(pc(i,j+1,1)*aepc(i-1,j-1)+pc(i,j-1,1)*awpc(i-1,j-1)+pc(i+1,j,1)*anpc(i-1,j-1)+...
     
224. pc(i-1,j,1)*aspc(i-1,j-1)-mp(i-1,j-1))/appc(i-1,j-1); %求出p'
     
225. end
     
226. end
     
227. end
     
228. ref=pc(2,2,1);
     
229. for i=2:ny+1
     
230. for j=2:nx+1
     
231. pc(i,j,1)=pc(i,j,1)-ref;
     
232. end
     
233. end
     
234. for i=1:ny
     
235. for j=1:nx
     
236. p(i,j,1)=p(i,j,1)+0.01*pc(i+1,j+1,1);
     
237. end
     
238. end
     
239. 
     
240. for i=2:ny+1
     
241. for j=2:nx
     
242. u(i,j,1)=u(i,j,1)+dx/apu(i,j)*(pc(i,j,1)-pc(i,j+1,1));
     
243. end
     
244. end
     
245. for i=2:ny
     
246. for j=2:nx
     
247. v(i,j,1)=v(i,j,1)+dy/apv(i,j)*(pc(i,j,1)-pc(i+1,j,1));
     
248. end
     
249. end
     
250. err=norm(mp(:));
     
251. if err<1e-5;break, end
     
252. for i=2:ny-1
     
253. for j=2:nx-1
     
254. mp(i,j)=raw*dx*((u(i,j+1,1)-u(i,j,1))+(v(i+1,j,1)-v(i,j,1)));
     
255. end
     
256. end
     
257. end
     
258. 
     
259. contourf(u(1:31,:,1),30);axis equal
     
260. set(gcf,'color','w')

复制代码

上面的方程比较简单，再附上一个，其实通过矩阵整体迭代求解可以适当加快求解过程，
matlab的向量迭代优势可以进一步体现出来，流场不管是否迎风模式，都可以采用。
看到有考虑多重求解雷诺方程的，雷诺方程本质上和Navier-Stokes是一样的，只是一个变种。
做这方面的都可以参与讨论。

scott198510

1. clear all;clc;close all;
   
2. %initialization
   
3. % Guess a temporary pressure and velocity fields
   
4. N=100;
   
5. a=0.0525;
   
6. R=15;
   
7. p_star = zeros(N);
   
8. u_star =a*ones(N);
   
9. v_star = zeros(N);
   
10. p_prime = zeros(N);
    
11. p_prime2= zeros(N);
    
12. u_prime = zeros(N);
    
13. v_prime = zeros(N);
    
14. p = p_star+p_prime;
    
15. u = u_star+u_prime;
    
16. v = v_star+v_prime;
    
17. for i=1:N
    
18.   for j=1:N
    
19.       if (i-N/2)^2+(j-N/2)^2<=R^2
    
20.         u_star(i,j)=0;
    
21.         v_star(i,j)=0;
    
22.       end
    
23.   end
    
24. end
    
25. alpha=0.4;
    
26. dt=0.004;
    
27. dx=0.4;
    
28. dy=0.4;
    
29. ro=1000;
    
30. MU=10;
    
31. for itr = 1:4000
    
32.     for i = 2:N-1
    
33.         for j = 2:N-1
    
34.            if (i-N/2)^2+(j-N/2)^2>R^2
    
35.             A_star1 = MU*((u(i+1,j)-2*u(i,j)+u(i-1,j))/dx^2+(u(i,j+1)-2*u(i,j)+u(i,j-1))/dy^2);
    
36.             A_star2 = u(i,j)*(u(i+1,j)-u(i-1,j))/2/dx+v(i,j)*(u(i,j+1)-u(i,j-1))/2/dy;
    
37.             A_star=A_star1-A_star2;
    
38.             B_star1 =  MU*((v(i+1,j)-2*v(i,j)+v(i-1,j))/dx^2 + (v(i,j+1)-2*v(i,j)+v(i,j-1))/dy^2);
    
39.             B_star2 =  u(i,j)*(v(i+1,j)-v(i-1,j))/2/dx+v(i,j)*(v(i,j+1)-v(i-1,j))/2/dy;
    
40.             B_star=B_star1-B_star2;
    
41.             u_star(i,j) = u_star(i,j) + A_star*dt-dt/dx*(p_star(i+1,j) - p_star(i,j))/ro;
    
42.             v_star(i,j) = v_star(i,j) + B_star*dt-dt/dx*(p_star(i,j+1) - p_star(i,j))/ro;
    
43.            end
    
44.         end
    
45.     end
    
46.     u = u_star;
    
47.     v = v_star;
    
48. 
    
49.     for itrate = 1:2000
    
50.         for i = 2:N-1
    
51.             for j = 2:N-1
    
52.                 p_prime2(i,j) = 0.25*(p_prime(i-1,j)+p_prime(i+1,j)+p_prime(i,j-1)+p_prime(i,j+1))+...
    
53.                     dx*(u_star(i-1,j)-u_star(i+1,j)+v_star(i,j-1)-v_star(i,j+1))/2;
    
54.             end
    
55.         end
    
56.         mx = max(max(abs(p_prime-p_prime2)));
    
57.         p_prime = p_prime2;
    
58.         if mx<0.001, break, end
    
59.     end
    
60. 
    
61.     p = p_star + alpha*p_prime;
    
62.     p_star = p;
    
63. 
    
64. end
    
65. 
    
66. figure,surf(p) %  pressure distribution in cavity
    
67. figure,contour(p,100) % pressure distribution in cavity
    
68. figure,surf(u) % u component of velocity
    
69. figure,contour(u,100) %u component of velocity
    
70. figure,surf(v) %v component of velocity
    
71. figure,contour(v,100)%v component of velocity
    
72. 
    
73. S=ones(N,N);
    
74. for i=1:N
    
75.   for j=1:N
    
76.       if (i-N/2)^2+(j-N/2)^2<=R^2
    
77.         S(i,j)=0;
    
78.       end
    
79.   end
    
80. end
    
81. figure,imshow(S);
    
82. hold on
    
83. xmax=N;
    
84. ymax=N;
    
85. nx=2;
    
86. ny=2;
    
87. [x,y]=meshgrid(1:nx:xmax,1:ny:ymax);
    
88. quiver(x,y,u(1:nx:xmax,1:ny:ymax)',v(1:nx:xmax,1:ny:ymax)',3);
    
89. axis equal;
    

复制代码

控制方程如图中所示，忽略重力作用gx,gy，模拟流体中含有固体障碍物时，流场分布