124 Computational Modelling in Hydraulic and Coastal Engineering
for i=2:nx-2
for j=3:ny-1
edd=(ev(i,j)+ev(i,j-1))/2;
h1=fm(i,j,df);
h5=fm(i,j-1,df);
if h1>0 && h5>0
hm=(h1+h5)/2;
um=(u(i,j)+u(i+1,j)+u(i,j-1)+u(i+1,j-1))/4;
surfy=fs*wy*sqrt(wx^2+wy^2)/hm;
adv=-v(i,j)*(v(i,j+1)-v(i,j-1))/2/
Dd-um*(v(i+1,j)-v(i-1,j))/2/Dd;
cor=-um*cf;
frict=-fb*sqrt(v(i,j)^2+um^2)*v(i,j)/hm;
pres=-g*(z(i,j)-z(i,j-1))/Dd;
diff=evd*(-4*v(i,j)+v(i,j-1)+v(i,j+1)+v(i1,j)+v(i+1,j))/Dd^2;
vn(i,j)=v(i,j)+Dt*(surfy+cor+pres+frict+diff
+adv);
end
end
end
% Boundary conditions for velocity v:
for i=2:nx-2
h6=fm(i,2,df);
if h6>0.1
vn(i,2)=-z(i,2)*sqrt(g/h6);
end
h7=fm(i,ny-2,df);
if h7>0.1
vn(i,ny-1)=z(i,ny-2)*sqrt(g/h7);
end
end
for j=2:ny-2
vn(1,j)=vn(2,j);
vn(nx-1,j)=vn(nx-2,j);
end
% Updating the velocities u and v;
for i=1:nx
for j=1:ny
u(i,j)=un(i,j);
v(i,j)=vn(i,j);
end
end
% Calculation of the kinetic energy;
ke=0;
for i=1:nx
for j=1:ny
ke=ke+u(i,j)^2+v(i,j)^2;
end
end