***Computational Hydraulics***

--- 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)