114 Computational Modelling in Hydraulic and Coastal Engineering
cf=0.0001;
fb=0.001;
Dd=5;
Dt=0.5;
nx=70;
ny=40;
nt=14400;
load depths.txt
df=depths;
% Subroutine h=fm(i,j,df): For definition of the marina
breakwater;
% Initial conditions;
for i=1:nx
for j=1:ny
u(i,j)=0;
un(i,j)=0;
v(i,j)=0;
vn(i,j)=0;
z(i,j)=0;
end
end
% Main program;
for k=1:nt
% Estimation of the eddy viscosity;
for i=2:nx-1
for j=2:ny-1
t1=(u(i+1,j)-u(i,j))/Dd-(v(i,j+1)-v(i,j))/Dd;
t2=(u(i,j+1)+u(i+1,j+1)-u(i,j-1)-u(i+1,j-1))/2/
Dd+(v(i+1,j+1)+v(i+1,j)-v(i-1,j+1)-v(i-1,j))/4/Dd;
ev(i,j)=sp*Dd^2*sqrt(t1^2+t2^2);
if ev(i,j) < 0.05;
ev(i,j)=0.05;
end
end
end
% Solution of the continuity equation;
for i=2:nx-2
for j=2:ny-2
h1=fm(i,j,df);
h2=fm(i-1,j,df);
h3=fm(i+1,j,df);
h4=fm(i,j+1,df);
h5=fm(i,j-1,df);
if h1>0.1
hl=(h1+h2)/2;
hr=(h1+h3)/2;
ho=(h4+h1)/2;
hu=(h1+h5)/2;
z(i,j)=z(i,j)-Dt/Dd*(u(i+1,j)*hru(i,j)*hl+v(i,j+1)*ho-v (i,j)*hu);