Free surface flows 115
end
end
end
% Solution x-axis equilibrium equation;
for i=3:nx-2
for j=2:ny-2
evd=(ev(i,j)+ev(i-1,j))/2;
h1=fm(i,j,df);
h2=fm(i-1,j,df);
if h1>0 && h2>0
hm=(h1+h2)/2;
vm=(v(i,j)+v(i,j+1)+v(i-1,j)+v(i-1,j+1))/4;
adv=-u(i,j)*(u(i+1,j)-u(i-1,j))/2/
Dd-vm*(u(i,j+1)-u(i,j-1))/2/Dd;
cor=vm*cf;
frict=-fb*sqrt(u(i,j)^2+vm^2)*u(i,j)/hm;
pres=-g*(z(i,j)-z(i-1,j))/Dd;
diff=evd*(-4*u(i,j)+u(i+1,j)+u(i1,j)+u(i,j+1)+u(i,j-1))/Dd^2;
un(i,j)=u(i,j)+Dt*(cor+pres+frict+diff+adv);
end
end
end
% Boundary conditions for velocity u;
for j=2:ny-2
h8=fm(2,j,df);
if h8>0.1
% Constant flow boundary condition;
un(2,j)=cv;
end
h9=fm(nx-2,j,df);
if h9>0.1
% Open sea boundary condition;
un(nx-1,j)=z(nx-2,j)*sqrt(g/h9);
end
end
for i=2:nx-2
un(i,1)=un(i,2);
un(i,ny-1)=un(i,ny-2);
end
% Solution of the y-axis equilibrium condition;
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;