---
# System prepended metadata

title: '***Computational Hydraulics***'
tags: [期中作業]

---

---
tags: 期中作業
---
# ***Computational Hydraulics***
## **Group 01:** N88071053 洪豪男 N88091011 黃彥鈞
## Shallow Water Equations (2D) on horizontal plane

$\quad\quad\dfrac{\partial h}{\partial t} + 
\dfrac{\partial}{\partial x} \left(h\overline{u}\right)+\dfrac{\partial}{\partial y} \left(h\overline{v}\right) = 0$
$\quad\quad\dfrac{\partial}{\partial t}\left(h\overline{u}\right) + 
\dfrac{\partial}{\partial x} \left(h\overline{u}^2 + \dfrac{1}{2}gh^2\right)+ 
\dfrac{\partial}{\partial y} \left(h\overline{u}\,\overline{v}\right)  = -hg_z\dfrac{\partial z_B}{\partial x}- \dfrac{\tau_{bx}}{\rho}$
$\quad\quad\dfrac{\partial}{\partial t}\left(h\overline{v}\right) + 
\dfrac{\partial}{\partial x} \left(h\overline{u}\,\overline{v}\right)+ 
\dfrac{\partial}{\partial y} \left(h\overline{v}^2 + \dfrac{1}{2}gh^2\right)  = -hg_z\dfrac{\partial z_B}{\partial y}- \dfrac{\tau_{by}}{\rho}$

---
1. on horizontal plane -> $-hg_z\dfrac{\partial z_B}{\partial x}$ and $-hg_z\dfrac{\partial z_B}{\partial y}$ = 0
2. assume : drag can ignore -> $\dfrac{\tau_{bx}}{\rho}$ and $\dfrac{\tau_{by}}{\rho}$= 0

----

```
%% Midterm HW#  -  2D shallow water equation : 
% 
% NOC scheme
% 2D Shallow Water Equation : 
%
% Assumptions : 
% Horizontal plane --> -hgz(zB_x) and -hgz(zB_y) = 0
% Ignore Drag force --> tau/tho = 0
%
% h_t + hu_x + hv_y = 0
% hu_t + (hu^2 + 0.5*g*h^2)_x + huv_y = 0 
% hv_t + huv_x + (hv^2 + 0.5*g*h^2)_y = 0
% 
% U(:,:,1) = h
% U(:,:,2) = hu
% U(:,:,3) = hv
% U(:,:,4) = u = U(:,:,2)/U(:,:,1)
% U(:,:,5) = v = U(:,:,3)/U(:,:,1)


close all ;
clc ;
clear all ;
mainFolder = pwd ;
%
% Sam Y. J. Huang, 2020.11.05
% Sam Y. J. Huang, 2020.12.11 modified
% Sam Y. J. Huang, 2021.01.02 modified

%% setting

x = -20:1:20; 
y = -20:1:20; 
dx = x(2)-x(1);
dy = y(2)-y(1);
[xx,yy] = meshgrid(x,y);
nx = length(x) ;
ny = length(y) ;
U0 = zeros(nx,ny,5);
g = 9.8; 
CFL = 0.1;
xpos0 = 0.0;
ypos0 = 0.0;
vex0 = 0.0;
vey0 = 0.0;
Rwave = 5.0 ;
H0 = 1.0 ; 
Hamp = 1.5;
Vel_bgn = 0.01;

% initial condition : sinus wave

for i = 1:length(x)
    for j = 1:length(y)
        xtmp = x(i) - xpos0;
        ytmp = y(j) - ypos0;
        dist = sqrt(xtmp^2+ytmp^2);
    if  dist <= Rwave
        tmp = dist/Rwave;
        U0(i,j,1) =  Hamp*(1.0-tmp^2)+H0; % h
        U0(i,j,4) =  vex0; % u
        U0(i,j,5) =  vey0; % v
    else
        U0(i,j,1) =  0.0 + H0; % h
        U0(i,j,4) =  0.1; % u
        U0(i,j,5) =  0.1; % v
    end
    U0(i,j,2) = U0(i,j,1)*U0(i,j,4); % hu
    U0(i,j,3) = U0(i,j,1)*U0(i,j,5); % hv
    end
end
hcut = 1.0E-8; % E-12

figure(1);
mesh(x,y,U0(:,:,1)), colormap gray, axis([-20 20 -20 20 0 3])
title('initial condition')
xlabel x, ylabel y; zlabel h;

%% Calculation

U = U0;
u_two = zeros(size(U));

H_x  = zeros(nx,ny);
H_y  = zeros(nx,ny);
HU_x = zeros(nx,ny);
HU_y = zeros(nx,ny);
HV_x = zeros(nx,ny);
HV_y = zeros(nx,ny);
U_x  = zeros(nx,ny);
U_y  = zeros(nx,ny);
V_x  = zeros(nx,ny);
V_y  = zeros(nx,ny);

stepNo = 0;
TimeT = 0.0;

figure(2);
for nstep = 1:200
    disp(nstep)
    stepNo = stepNo + 1 ;
    lmbdx1 = abs(U(:,:,4));
    lmbdy1 = abs(U(:,:,5));
    lmb2 = sqrt(U(:,:,1)*g);
    lmbdx = lmbdx1+lmb2;
    lmbdy = lmbdy1+lmb2;
    lmbdx = abs(U(:,:,4)) + sqrt(U(:,:,1)*g);
    lmbdy = abs(U(:,:,5)) + sqrt(U(:,:,1)*g);
    ax_max = max(max(lmbdx));
    ay_max = max(max(lmbdy));
    a = max(ax_max,ay_max);
    dt = CFL*dx/a;
    TimeT = TimeT + dt;
    
   
    for jj = 1:nx-1
         for ii = 1:ny-1
           
            H_x(jj,ii)  = 0.5*(U(jj+1,ii+1,1)+U(jj,ii+1,1)) - ...
                          dt/dx*(U(jj+1,ii+1,2)-U(jj  ,ii+1,2));
            
            H_y(jj,ii)  = 0.5*(U(jj+1,ii+1,1)+U(jj+1,ii,1)) - ...
                          dt/dx*(U(jj+1,ii+1,3)-U(jj+1,ii  ,3));
            
            HU_x(jj,ii) = 0.5*(U(jj+1,ii+1,2)+U(jj,ii+1,2)) - ...
                          dt/dx*(U(jj+1,ii+1,2).*U(jj+1,ii+1,4)+...
                                 0.5*g*U(jj+1,ii+1,1).^2-U(jj,ii+1,2).*...
                                 U(jj,ii+1,4)-0.5*g*U(jj,ii+1,1).^2);
            
            HU_y(jj,ii) = 0.5*(U(jj+1,ii+1,2)+U(jj+1,ii,2)) - ...
                          dt/dy*(U(jj+1,ii+1,2).*U(jj+1,ii+1,5)-...
                                 U(jj+1,ii,2).*U(jj+1,ii,5));
            
            HV_x(jj,ii) = 0.5*(U(jj+1,ii+1,3)+U(jj,ii+1,3)) - ...
                          dt/dx*(U(jj+1,ii+1,2).*U(jj+1,ii+1,5)-...
                                 U(jj,ii+1,2).*U(jj,ii+1,5));
            
            HV_y(jj,ii) = 0.5*(U(jj+1,ii+1,3)+U(jj+1,ii,3)) - ...
                          dt/dy*(U(jj+1,ii+1,3).*U(jj+1,ii+1,5)+...
                                 0.5*g*U(jj+1,ii+1,1).^2-U(jj+1,ii,3).*...
                                 U(jj+1,ii,5)-0.5*g*U(jj+1,ii,1).^2);
            
            U_x = HU_x(jj,ii)./H_x(jj,ii);
            U_y = HU_y(jj,ii)./H_y(jj,ii);
            V_x = HV_x(jj,ii)./H_x(jj,ii);
            V_y = HV_y(jj,ii)./H_y(jj,ii);

         end
    end

    lmbdx = abs(sqrt(U_x^2+U_y^2)) + sqrt(U(:,:,1)*g);
    lmbdy = abs(sqrt(V_x^2+V_y^2)) + sqrt(U(:,:,1)*g);
    ax_max = max(max(lmbdx));
    ay_max = max(max(lmbdy));
    a = max(ax_max,ay_max);
    dt = CFL*dx/a;
    TimeT = TimeT + dt;
    
    % Second step u_two(jj) value at jj
    for jj = 2:nx-1
         for ii = 2:ny-1
u_two(jj,ii,1) = U(jj,ii,1)-dt/dx*(HU_x(jj,ii-1)-HU_x(jj-1,ii-1))...
                           -dt/dy*(HV_y(jj-1,ii)-HV_y(jj-1,ii-1));
           
u_two(jj,ii,2) = U(jj,ii,2) - dt/dx*((HU_x(jj,ii-1).^2./H_x(jj,ii-1)+...
                 0.5*g*H_x(jj,ii-1).^2)-(HU_x(jj-1,ii-1).^2./...
                 H_x(jj-1,ii-1)+0.5*g*H_x(jj-1,ii-1).^2))- dt/dy*...
                 (HU_y(jj-1,ii).*HV_y(jj-1,ii)./H_y(jj-1,ii)-...
                 HU_y(jj-1,ii-1).*HV_y(jj-1,ii-1)./H_y(jj-1,ii-1));
                        
u_two(jj,ii,3) = U(jj,ii,3) - dt/dx*(HV_x(jj,ii-1).*HU_x(jj,ii-1)./...
                 H_x(jj,ii-1)-HV_x(jj-1,ii-1).*HU_x(jj-1,ii-1)./...
                 H_x(jj-1,ii-1))-dt/dy*((HV_y(jj-1,ii).^2./H_y(jj-1,ii)+...
                 0.5*g*H_y(jj-1,ii).^2)-(HV_y(jj-1,ii-1).^2./...
                 H_y(jj-1,ii-1) + 0.5*g*H_y(jj-1,ii-1).^2));
         
            u_two(jj,ii,4) = u_two(jj,ii,2)./(u_two(jj,ii,1)+hcut);
            u_two(jj,ii,5) = u_two(jj,ii,3)./(u_two(jj,ii,1)+hcut);
            
         end
    end
    
% B.C.

    % h
    u_two(:,1,:)   = u_two(:,end-1,:);
    u_two(:,end,:) = u_two(:,end-1,:);
    u_two(1,:,:)   = u_two(end-1,:,:);
    u_two(end,:,:) = u_two(end-1,:,:);    
    
    % hu
    u_two(:,1,2)   =  u_two(:,end-1,2);
    u_two(:,end,2) =  u_two(:,end-1,2);    
    u_two(1,:,2)   = u_two(end-1,:,2);
    u_two(end,:,2) = u_two(end-1,:,2);    
    
    % hv
    u_two(:,1,3)   = u_two(:,end-1,3);
    u_two(:,end,3) = u_two(:,end-1,3);    
    u_two(1,:,3)   =  u_two(end-1,:,3);
    u_two(end,:,3) =  u_two(end-1,:,3);  
    
    U = u_two;
    
    mesh(x,y,U(:,:,1)) ;
    colormap gray
    axis([-20 20 -20 20 0 3]);
    xlabel x ;
    ylabel y ;
    zlabel('h') ;
    text(-10,10,3,['Step Number = ', num2str(stepNo)])
    text(-10,10,2.6,['Time = ', num2str(TimeT)])
    pause(0.01);
    
   
    
end





 ```


![](https://i.imgur.com/fCUWJvK.jpg)
![](https://i.imgur.com/GmNG07l.jpg)
![](https://i.imgur.com/KvJqwdm.jpg)








































