Free surface flows 123
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)
*hr-u(i,j)*hl+v(i,j+1)*ho-v(i,j)*hu);
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;
surfx=fs*wx*sqrt(wx^2+wy^2)/hm;
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*(surfx+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
un(2,j)=-z(2,j)*sqrt(g/h8);
end
h9=fm(nx-2,j,df);
if h9>0.1
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;