132 Computational Modelling in Hydraulic and Coastal Engineering
hb(k)=hbo;
hun(k)=huo;
hbn(k)=hbo;
end
for k1=1:nx
zb(k1)=0;
uu(k1)=0;
ub(k1)=0;
uun(k1)=0;
ubn(k1)=0;
end
% Main program;
for i=1:nt
% Continuity equation for the upper fluid layer;
for j=2:nx-2
qul(j)=uu(j)*(hu(j)+hu(j-1))/2;
qur(j)=uu(j+1)*(hu(j)+hu(j+1))/2;
hun(j)=hu(j)-(Dt/Dx)*(qur(j)-qul(j));
end
% Continuity equation for the lower fluid layer;
for j=2:nx-2
qbl(j)=ub(j)*(hb(j)+hb(j-1))/2;
qbr(j)=ub(j+1)*(hb(j)+hb(j+1))/2;
hbn(j)=hb(j)-(Dt/Dx)*(qbr(j)-qbl(j));
end
% Updating of layer depth values;
for j=2:nx-2
hu(j)=hun(j);
hb(j)=hbn(j);
end
% Force equilibrium equation for the upper layer;
for j=3:nx-2
hum(j)=(hu(j)+hu(j-1))/2;
ts(j)=fs*w*abs(w)/hum(j);
ti(j)=fi*(uu(j)-ub(j))*abs(uu(j)-ub(j))/hum(j);
ed(j)=ev*(uu(j-1)-2*uu(j)+uu(j+1))/Dx^2;
uun(j)=uu(j)-Dt*g*(hu(j)+hb(j)+zb(j)-hu(j-1)-hb(j1)-zb(j-1))/Dx+Dt*(ed(j)+ts(j)-ti(j));
end
% Force equilibrium equation for the lower layer;
for j=3:nx-2
hbm(j)=(hb(j)+hb(j-1))/2;
ti(j)=fi*(uu(j)-ub(j))*abs(uu(j)-ub(j))/hbm(j);
tb(j)=fb*ub(j)*abs(ub(j))/hbm(j);
ed(j)=ev*(ub(j-1)-2*ub(j)+ub(j+1))/Dx^2;
ubn(j)=ub(j)-Dt*g*((hu(j)+hb(j)+zb(j)-hu(j-1)-hb(j1)-zb(j-1))/Dx-rdd*(hu(j)-hu(j-1))/
Dx)+Dt*(ed(j)+ti(j)-tb(j));
end
% Updating of fluid velocity values;