SHARE
TWEET

Untitled

a guest Apr 24th, 2019 60 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. OK, I Understand
 
Top