Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

Untitled

By: a guest on May 2nd, 2012  |  syntax: None  |  size: 6.07 KB  |  hits: 17  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. % WAVE SIMULATION
  2.  
  3. Hs=15;  %significant wave height
  4. Tz=15;  % zero upcrossing time
  5. g=9.81;
  6. d=500;  %water depth
  7.  
  8. f=linspace(0.01,.2,100); % frequency
  9. f0=((1/(2*pi))*sqrt(0.161*g/Hs));% peak frequency
  10. q=f/f0;
  11. Sf=(5*Hs^2/(16*f0))./q.^5.*exp(q.^-4.*-1.25); %energy density for unidirectional random wave
  12. w=2*pi*f;
  13. %Sf=(Tz*Hs^2/(8*(pi^2))*((Tz*f).^-5).*exp((-1/pi)*((Tz*f)).^-4));
  14. s=1;% spreading factor
  15. Ns=(gamma(s))/((sqrt(pi))*(gamma(s+0.5)));%normalising factor
  16. gg=1;
  17. thet0=0;
  18. for thet=(-pi:pi/36:pi);% direction
  19. if abs(thet-thet0) <=(pi/2)
  20.     hh=Ns.*((cos(thet)).^2);
  21. else
  22.    hh=0;  
  23. end
  24. disp(hh);
  25. G(gg)=hh;
  26. the(gg)=thet;
  27. gg=gg+1;
  28.  
  29. end
  30. th=(the*180/pi)';
  31. SSf=(Sf'*G);%energy density for multidirectional random wave
  32. figure(1)
  33. surf(th,f,SSf)
  34. xlabel('theta(deg)')
  35. ylabel('f(Hz)')
  36. zlabel('Sf')
  37. figure(1)
  38. grid on
  39. title('Plot of Wave Spectrum')
  40.  
  41. x=0;
  42. y=0;
  43. deltaf=0.001;
  44. deltath=pi/18;
  45. eta=zeros(length(f),length(the)); %wave elevation
  46. u1=zeros(length(f),length(the)); % hoizontal water particle velocity in x direction
  47. v1=zeros(length(f),length(the)); % hoizontal water particle velocity in y direction
  48. w1=zeros(length(f),length(the)); % vertical water particle velocity
  49. ua1=zeros(length(f),length(the));% hoizontal water particle acceleration in x direction
  50. va1=zeros(length(f),length(the));% hoizontal water particle acceleration in y direction
  51. wa1=zeros(length(f),length(the));% vertcal water particle acceleration
  52. Fx=zeros(length(f),length(the));% force in x direction
  53. Fy=zeros(length(f),length(the));% force in y direction
  54. Fz=zeros(length(f),length(the));% force in z direction
  55. Mx=zeros(length(f),length(the));% moment in x direction
  56. My=zeros(length(f),length(the));% moment in y direction
  57. Mz=zeros(length(f),length(the));% moment in z direction
  58. m=1;
  59.  
  60. T = 0;
  61.  
  62. for x=(1:1:100);
  63.     for y = (1:100);
  64.         eta(x,y)=0;
  65.         u1(x,y)=0;
  66.         v1(x,y)=0;
  67.         w1(x,y)=0;
  68.         ua1(x,y)=0;
  69.         va1(x,y)=0;
  70.         wa1(x,y)=0;
  71.     for j=1:(length(f)-1); % loop for frequecy run
  72.            n=1;
  73.            for i=1:(length(the)-1);% loop for direction run
  74.         Sn(j,i)=(SSf(j,i)+SSf((j),(i+1)))/2;% energy densty
  75.         Aa(j,i)=sqrt(2*Sn(j,i)*deltaf*deltath);% amplitude
  76.         Ff(j)=(f(j)+f((j+1)))/2;
  77.         Tpp(j)=1/Ff(j);
  78.         phi(j,i)=rand*2*pi;%phase shift
  79.    
  80.      Lo(j)=1.56.*Tpp(j).^2;%initial wave length
  81.  
  82.     L(j)=(g.*(Tpp(j).^2)/(2*pi)).*tanh(2*pi*d./Lo(j));%wave length
  83.     Lo(j)=L(j);
  84.  
  85.         kj(j)=(2*pi)./L(j);% wave number
  86.        
  87.                
  88.         eta(x,y) = eta(x,y) + Aa(j,i)*(cos((kj(j)*x*cos(the(i)))+(kj(j)*y*sin(the(i)))-(2*pi*Ff(j)*T)+(phi(j,i))));
  89.         u1(x,y)=u1(x,y) +((2*pi)*Aa(j,i)*2*pi*Ff(j)*cosh(kj(j)*(y))*cos(thet0)*(cos((kj(j)*x*cos(the(i)))+(kj(j)*y*sin(the(i)))-(2*pi*Ff(j)*T)+(phi(j,i))))/sinh(kj(j)*(d+eta(j,i))));
  90.         v1(x,y)=v1(x,y) +((2*pi)*Aa(j,i)*2*pi*Ff(j)*cosh(kj(j)*(y))*sin(thet0)*(cos((kj(j)*x*cos(the(i)))+(kj(j)*y*sin(the(i)))-(2*pi*Ff(j)*T)+(phi(j,i))))/sinh(kj(j)*(d+eta(j,i))));
  91.         w1(x,y)=w1(x,y) +((2*pi)*Aa(j,i)*2*pi*Ff(j)*sinh(kj(j)*(y))*(sin((kj(j)*x*cos(the(i)))+(kj(j)*y*sin(the(i)))-(2*pi*Ff(j)*T)+(phi(j,i))))/sinh(kj(j)*(d+eta(j,i))));
  92.         ua1(x,y)=ua1(x,y) +((4*(pi^2))*Aa(j,i)*(2*pi*Ff(j))^2*cosh(kj(j)*(y))*cos(thet0)*(sin((kj(j)*x*cos(the(i)))+(kj(j)*y*sin(the(i)))-(2*pi*Ff(j)*T)+(phi(j,i))))/sinh(kj(j)*(d+eta(j,i))));
  93.         va1(x,y)=va1(x,y) +((4*(pi^2))*Aa(j,i)*(2*pi*Ff(j))^2*cosh(kj(j)*(y))*sin(thet0)*(sin((kj(j)*x*cos(the(i)))+(kj(j)*y*sin(the(i)))-(2*pi*Ff(j)*T)+(phi(j,i))))/sinh(kj(j)*(d+eta(j,i))));
  94.         wa1(x,y)=wa1(x,y) +((4*(pi^2))*Aa(j,i)*(2*pi*Ff(j))^2*sinh(kj(j)*(y))*(sin((kj(j)*x*cos(the(i)))+(kj(j)*y*sin(the(i)))-(2*pi*Ff(j)*T)+(phi(j,i))))/sinh(kj(j)*(d+eta(j,i))));
  95.        Fx(j+1,i+1)=((Vol*rho*Cm*ua1(j,i))+(0.5*Cd*rho*((4*D*dl)+(2*pl*ph))*u1(j,i)*abs(u1(j,i))));
  96.        Fy(j+1,i+1)=((Vol*rho*Cm*va1(j,i))+(0.5*Cd*rho*((4*D*dl)+(2*pl*ph))*v1(j,i)*abs(v1(j,i))));
  97.        Fz(j+1,i+1)=((Vol*rho*Cm*wa1(j,i))+(0.5*Cd*rho*((4*D*dl)+(2*pl*ph))*w1(j,i)*abs(w1(j,i))));
  98.        
  99.        Mx(j+1,i+1)=Fx(j+1,i+1)*(d+z);
  100.        My(j+1,i+1)=Fy(j+1,i+1)*(d+z);
  101.        Mz(j+1,i+1)=0;
  102.         n=n+1;
  103.        
  104.          
  105.        
  106. %    U(m,n)=u1(j+1,i+1);
  107. %    V(m,n)=v1(j+1,i+1);
  108. %    W(m,n)=w1(j+1,i+1);
  109. %    Ua(m,n)=ua1(j+1,i+1);
  110. %    Va(m,n)=va1(j+1,i+1);
  111. %    Wa(m,n)=wa1(j+1,i+1);
  112.    FX(m,n)=Fx(j+1,i+1);
  113.    FY(m,n)=Fy(j+1,i+1);
  114.    FZ(m,n)=Fz(j+1,i+1);
  115.    MX(m,n)=Mx(j+1,i+1);
  116.    MY(m,n)=My(j+1,i+1);
  117.    MZ(m,n)=Mz(j+1,i+1);
  118.    
  119.       end
  120.      end
  121.  
  122.        t(m)=T;
  123.       F=[FX FY FZ MX MY MZ];
  124.        
  125.        
  126.       E(x,y) = eta(x,y);
  127.       U(x,y) = sum(sum(u1));
  128.       V(x,y) = sum(sum(v1));
  129.       W(x,y) = sum(sum(w1));
  130.       Ua(x,y) = sum(sum(ua1));
  131.       Va(x,y) = sum(sum(va1));
  132.       Wa(x,y) = sum(sum(wa1));
  133.       m=m+1;
  134.     end    
  135. end
  136. F=F';
  137. x=[1:100];
  138. y=[1:100];
  139. figure(2)
  140. surf(x,y,E)
  141. xlabel('Theta(deg)')
  142. ylabel('Time(sec)')
  143. zlabel('Surface elevation(m)')
  144. grid on
  145. title('Random elevation of surface wave')
  146.  
  147. x=[1:100];
  148. y=[1:100];
  149. figure(3)
  150. surf(x,y,U)
  151. xlabel('Theta(deg)')
  152. ylabel('Time(sec)')
  153. zlabel('Horizontal velocity X(m)')
  154. grid on
  155. title('Time history of Horizontal velocity, X')
  156.  x=[1:100];
  157. y=[1:100];
  158. figure(4)
  159. surf(x,y,V)
  160. xlabel('Theta(deg)')
  161. ylabel('Time(sec)')
  162. zlabel('Horizontal velocity Y(m)')
  163. grid on
  164. title('Time history of Horizontal velocity, Y')
  165.  x=[1:100];
  166. y=[1:100];
  167. figure(5)
  168. surf(x,y,W)
  169. xlabel('Theta(deg)')
  170. ylabel('Time(sec)')
  171. zlabel('Vertical velocity(m)')
  172. grid on
  173. title('Time history of Vertical velocity')
  174. x=[1:100];
  175. y=[1:100];
  176. figure(6)
  177. surf(x,y,Ua)
  178. xlabel('Theta(deg)')
  179. ylabel('Time(sec)')
  180. zlabel('Horizontal acceleration, X(m)')
  181. grid on
  182. title('Time history of Horizontal acceleration, X')
  183. x=[1:100];
  184. y=[1:100];
  185. figure(7)
  186. surf(x,y,Va)
  187. xlabel('Theta(deg)')
  188. ylabel('Time(sec)')
  189. zlabel('Horizontal acceleration, Y(m)')
  190. grid on
  191. title('Time history of Horizontal acceleration, Y')
  192.  x=[1:100];
  193. y=[1:100];
  194. figure(8)
  195. surf(x,y,Wa)
  196. xlabel('Theta(deg)')
  197. ylabel('Time(sec)')
  198. zlabel('Vertical acceleration(m)')
  199. grid on
  200. title('Time history of Vertical acceleration')
create a new version of this paste RAW Paste Data