Contaminant and sediment transport by advection and diffusion 253
y(i)=y(i)+Dt*(vv+(2*ay-1)*D+wb);
end
% Estimation of particle concentration in each cell;
for j=1:ny-1
for i=1:nx-1
c(i,j)=0;
end
end
for i=1:ipt
X=fix(x(i)/Dd)+1;
Y=fix(y(i)/Dd)+1;
if X1 && Y1
c(X,Y)=c(X,Y)+1;
end
end
% Estimation of the kinetic energy;
KE=0;
for i=2:nx-2
for j=2:ny-2
KE=KE+u(i,j)^2+v(i,j)^2;
end
end
KinE(k)=KE;
index=k
end
plot(1:k,KinE,'Linewidth',1.5);
xlabel('Time x 0.05 [s]')
ylabel('Kinetic energy u^2+v^2')
figure
[i,j]=meshgrid(1:1:nx,1:1:ny);
quiver(i,j,u',v');
xlabel('Horizontal distance [m]')
ylabel('Vertical distance [m]')
figure
contour(i,j,c')
xlabel('Horizontal distance [m]')
ylabel('Vertical distance [m]')
PROBLEM 8.9
Solve the same problem by making the suggested modifications while keeping the rest of the data constant:
1. Change the air pseudo-compressibility constant from 4 m/s to 1 and
8 m/s.
2. Change the eddy viscosity from 0.1 m 2 /s to 0.005, 0.05 and
0.2 m 2 /s.