• API
• FAQ
• Tools
• Archive
SHARE
TWEET Untitled a guest Apr 24th, 2019 72 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. clearvars
2. close all
3.
4. OP = [.25 .25 7 7,];
5.
6. xmin = -5;
7. xmax = 5;
8. ymin = -5;
9. ymax = 5;
10. spacing = 0.01;
11.
12. %Define value for vortices and cross flow
13. Gamma1 = -10;
14. Gamma2 = 8;
15.
16. U = 4;
17. a = 2;
18.
19. x = (xmin:spacing:xmax);
20. y = (ymin:spacing:ymax);
21.
22. xsize = numel(x);
23. ysize = numel(y);
24.
25. [x,y] = meshgrid(x,y);
26.
27. %time step
28. ts = 0.000001;
29.
30. %Initial vortex center location
31. yc1 = -1;
32. xc1= -1;
33. yc2 = 1;
34. xc2 = 1;
35.
36. %Initialize zero matrices
37. xloc1 = zeros(1,xsize); %x location of vortex 1
38. yloc1 = zeros(1,ysize); %y location of vortex 1
39. xloc2 = zeros(1,xsize); %x location of vortex 2
40. yloc2 = zeros(1,ysize); %y location of vortex 2
41.
42. u1 = zeros(ysize,xsize);
43. v1 = zeros(ysize,xsize);
44. u2 = zeros(ysize,xsize);
45. v2 = zeros(ysize,xsize);
46. u3 = zeros(ysize,xsize);
47. v3 = zeros(ysize,xsize);
48.
49. dPsidy1 = zeros(ysize,xsize);
50. dPsidx1 = zeros(ysize,xsize);
51. dPsidy2 = zeros(ysize,xsize);
52. dPsidx2 = zeros(ysize,xsize);
53. dPsidy3 = zeros(ysize,xsize);
54. dPsidx3 = zeros(ysize,xsize);
55.
56.
57. %%%%%%%%%%%%Calculate velocity caused by cross flow%%%%%%%%%%%%%%
58. Psi3 = U*((-sind(a).*x) + (cosd(a).*y));
59. for i = 2:(ysize-1)
60.     for j = 2:xsize-1
61.         %The y value changes as you increase the row (so i)
62.         %Essentially this is delta y over delta x
63.         dPsidy3(i,j) = (Psi3(i+1,j) - Psi3(i-1,j) )/ (2*spacing);
64.         dPsidx3(i,j) = (Psi3(i,j+1) - Psi3(i,j-1)) / (2*spacing);
65.     end
66. end
67.
68. u3 = dPsidy3;   %matrix of x velocity
69. v3 = -1* dPsidx3; %matrix of y velocity
70.
71. u3(u1>1000000) = 0;
72. u3(u1<-1000000) = 0;
73. v3(v1>1000000) = 0;
74. v3(v1<-1000000) = 0;
75. %[%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
76.
77.     r1 = sqrt(((x-xc1).^2)+((y-yc1).^2));
78.     r2 = sqrt(((x-xc2).^2)+((y-yc2).^2));
79.
80.         %Psi defined for vortices
81.     Psi1 =(Gamma1*log(r1)./(2*pi));
82.     Psi2 =(Gamma2*log(r2)./(2*pi));
83.
84.         for i = 2:(ysize-1)
85.             for j = 2:xsize-1
86.                 %The y value changes as you increase the row (so i)
87.                 %Essentially this is delta y over delta x
88.                 dPsidy1(i,j) = (Psi1(i+1,j) - Psi1(i-1,j) )/ (2*spacing);
89.                 dPsidx1(i,j) = (Psi1(i,j+1) - Psi1(i,j-1)) / (2*spacing);
90.             end
91.         end
92.
93.         for i = 2:(ysize-1)
94.             for j = 2:xsize-1
95.             %The y value changes as you increase the row (so i)
96.             %Essentially this is delta y over delta x
97.             dPsidy2(i,j) = (Psi2(i+1,j) - Psi2(i-1,j) )/ (2*spacing);
98.             dPsidx2(i,j) = (Psi2(i,j+1) - Psi2(i,j-1)) / (2*spacing);
99.             end
100.         end
101.
102.         u1 = dPsidy1;   %matrix of x velocity
103.         v1 = -1* dPsidx1; %matrix of y velocity
104.         u2 = dPsidy2;   %matrix of x velocity
105.         v2 = -1* dPsidx2; %matrix of y velocity
106.
107.         u1(u1>1000000) = 0;
108.         u1(u1<-1000000) = 0;
109.         v1(v1>1000000) = 0;
110.         v1(v1<-1000000) = 0;
111.
112.         u2(u2>1000000) = 0;
113.         u2(u2<-1000000) = 0;
114.         v2(v2>1000000) = 0;
115.         v2(v2<-1000000) = 0;
116.
117.
118.         xc1 = xc1 + (u2(n1, m1).*ts) + (u3(n1, m1).*ts);
119.         yc1 = yc1 + (v2(n1, m1).*ts) + (v3(n1, m1).*ts);
120.         xc2 = xc2 + (u1(n1, m1).*ts) + (u3(n1, m1).*ts);
121.         yc2 = yc2 + (v1(n1, m1).*ts) + (v3(n1, m1).*ts);
122.
123.
124.         xloc1(1, m1) = xc1;
125.         yloc1(1, m1) = yc1;
126.         xloc2(1, m1) = xc2;
127.         yloc2(1, m1) = yc2;
128.
129.
130.
131. %hold on
132. %contourf(x,y,Psi, 1000, 'LineStyle', 'none');
133. %colormap
134. %colorbar
135. %h = colorbar;
136. %title(h, '\psi')
137. %hold on;
138.
139.
140. %ph = streamslice(x,y,u1,v1);
141.
142. %set(ph,'color','k');
143.
144. %axis([xmin xmax ymin ymax]);
145. %set(gcf,'Units','inches','OuterPosition',OP);
146.
147. %title('Streamlines for vortices');
148. %xlabel('x location');
149. %ylabel('y location');
150. %hold off
151.
152. figure;
153.
154. plot(xloc1, yloc1);
155. hold on;
156. plot(xloc2, yloc2);
157. ylim([-5 5]);
158. xlim([-5 5]);
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top