# Untitled

By: a guest on May 2nd, 2012  |  syntax: None  |  size: 6.07 KB  |  hits: 17  |  expires: Never
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));
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')