122 Computational Modelling in Hydraulic and Coastal Engineering
% nt = Number of time steps;
clc; clear all; close all;
% Input data;
g=9.81;
ho=0;
wx=7;
wy=-7;
dif=1;
cf=0.0001;
fs=0.000003;
fb=0.01;
sp=0.3;
Dd=2000;
Dt=30;
nx=22;
ny=23;
nt=21600;
load ThermD.txt
df=ThermD;
% Import h=fm(i,j,df): Subroutine that defines depths:
% Initial conditions;
for i=1:nx
for j=1:ny
u(i,j)=0;
un(i,j)=0;
v(i,j)=0;
vn(i,j)=0;
z(i,j)=0;
end
end
% Main program;
for k=1:nt
% Estimation of the eddy viscosity;
for i=2:nx-1
for j=2:ny-1
t1=(u(i+1,j)-u(i,j))/Dd-(v(i,j+1)-v(i,j))/Dd;
t2=(u(i,j+1)+u(i+1,j+1)-u(i,j-1)-u(i+1,j-1))/2/Dd+
(v(i+1,j+1)+v(i+1,j)-v(i-1,j+1)-v(i-1,j))/4/Dd;
ev(i,j)=sp*Dd^2*sqrt(t1^2+t2^2);
if ev(i,j) < 1;
ev(i,j)=1;
end
end
end
% Solution of the continuity equation;
for i=2:nx-2
for j=2:ny-2
h1=fm(i,j,df);
h2=fm(i-1,j,df);
h3=fm(i+1,j,df);