Tutorial 27 – Shallow Water Equations

Posted on July 17, 2017 - Boundary Conditions, Postprocessing, Shallow Water, Transient Problems, Tutorials

This tutorial is intended for the full version of the toolbox.

Introduction

Problem definition The purpose of this tutorial is to show how to solve simplified, reduced to two dimensions Navier-Stokes Equations called shallow water or Saint-Venant equations. In order to use this simplification domain of phenomenon that we want to simulate has to be significantly smaller in vertical direction.

    $$\frac{\partial h}{\partial t}+H\nabla\cdot\vec{u}=0$$

    $$\frac{\partial\vec{u}}{\partial t}+f\left[\begin{array}{cc}0 & -1\\1 & 0\end{array}\right]\vec{u}=-g\nabla h-b\vec{u}$$

Where:

$h$ – Height of free surface (solution field)

$H$ – Mean height of free surface (constant)

$f$ – Coriolis coefficient

$g$ – Acceleration due to gravity

$b$ – Viscous drag coefficient

Shallow water equations are widely used in computational fluid dynamics, for example, to simulate propagation of tsunami waves.

Geometry

Since we are solving 3-D problem that is reduced to 2-D our geometry will be represented by surface shown below.

Code

Clear MATLAB Command Window and Workspace

clc;
clear;

% Import and display mesh
[p, e, t] = importMeshGmsh('shallowWater.msh');
displayMesh2D(p,t);


nnodes = size(p,2);

% Define problem parameters
swe.b = 0;
swe.f = 0;
swe.g = 1;
swe.H = 0.4;

% Define time-stepping parameters
dt = 0.005;
nstep = 1050;
% Init and plot solution
vx = zeros(nnodes,1);
vy = zeros(nnodes,1);
x = p(1,:); y = p(2,:);
h = 2*exp(-((x-2.5).^2 + (y-0.5).^2)/0.25)';
% movieObj = initMovie('shallowWater2D_vis_fine.avi',60);

% Begin time-stepping
for step = 1:nstep
% Assemble matrix
[M,F] = assembleShallowWaterMatrix2D(p,t,vx,vy,h,swe,dt);
% Impose boundary conditions at rigid walls
[M,F] = imposeShallowWaterBoundaryCondition2D(p,e,M,F,9,'wall');
% Solve equations and update solution
[h, vx, vy] = solveShallowWaterEquations(M,F,h,vx,vy);
% Plot current solution
figure(1)
trisurf(t(1:3,:)',p(1,:)',p(2,:)',h)
shading('interp')
colormap('jet')
colorbar
axis equal
zlim([-0.4 2])
% view([0, 90])
drawnow;
% writeVideo(movieObj,getframe(gcf));
disp(step)
end
% close(movieObj);

Tags