6 不可压缩牛顿流体流动

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