116 Computational Modelling in Hydraulic and Coastal Engineering
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*(cor+pres+frict+diff+adv);
end
end
end
% Open sea boundary conditions for velocity v:
for i=2:nx-2
h6=fm(i,2,df);
h7=fm(i,ny-2,df);
if h6<0.1
if h7<0.1
continue;
else
vn(i,ny-1)=z(i,ny-2)*sqrt(g/h7);
end
end
vn(i,2)=-z(i,2)*sqrt(g/h6);
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);
% Zero flux boundary condition;
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
KinE(k)=ke;
%
index=k
end
plot(1:k,KinE,'Color','k','Linewidth',1.5)
xlabel('Number of time steps')
ylabel('Total kinetic energy (u^2+v^2 on all grid points)')