252 Computational Modelling in Hydraulic and Coastal Engineering
for j=2:ny-2
for i=3:nx-2
un(i,j)=u(i,j)-Dt/Dd*cp^2*(h(i,j)-h(i-1,j))+ed*5
*Dt/Dd^2*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)4*u(i,j));
end
end
for j=1:ny
un(2,j)=-(h(2,j)-1)*cp/h(2,j);
un(nx-1,j)=(h(nx-2,j)-1)*cp/h(nx-2,j);
end
for i=1:nx
un(i,1)=un(i,2);
un(i,ny-1)=un(i,ny-2);
end
% Computation of velocities along the y-axis;
for j=3:ny-2
for i=2:nx-2
1))+ed*5*Dt/Dd^2*(v(i,j+1)+v(i,j-1)+v(i-1,j)+v(i+1,j)4*v(i,j))+Dt*(c(i,j)+c(i,j-1))/2/Dd^2*tau;
end
end
for j=1:ny
vn(1,j)=vn(2,j);
vn(nx-1,j)=vn(nx-2,j);
end
% Updating of the velocity values;
for i=1:nx
for j=1:ny
u(i,j)=un(i,j);
v(i,j)=vn(i,j);
end
end
% Particle relocation by advection, buoyancy and
diffusion;
for i=1:ipt
X=fix(x(i)/Dd)+1;
Y=fix(y(i)/Dd)+1;
if X1 && Y1
uu=u(X,Y)+(u(X+1,Y)-u(X,Y))*(x(i)-(X-1)*Dd)/Dd;
vv=v(X,Y)+(v(X,Y+1)-v(X,Y))*(y(i)-(Y-1)*Dd)/Dd;
end
tempx=x(i);
tempy=y(i);
ax=rand;
x(i)=x(i)+Dt*(uu+(2*ax-1)*D);
ay=rand;