User:Egm6322.s09.three.liu/project

 See my comments below. Egm6322.s09 14:29, 1 May 2009 (UTC)

=The Analytical and Numerical Solutions to the Shallow Water Equations=

Executive Summary
The analytical and numerical solutions to a hyperbolic PDE: linear shallow water equation are derived in this study. The analytical solution, which is solved by the method of separation of variables, consists of a steady part and a transient part. After the analytical solutions were derived, an implicit numerical scheme that applies conjugate gradient method and upwind method to solve the 2-D shallow water equation are developed. The numerical solution by this scheme is obtained and compared with the analytical solution. The comparison indicates that the two solutions are in good agreement.

Introduction
The Shallow Water Equations are widely used to model the propagation of disturbance in incompressible fluids such as tidal hydrodynamics, storm surge and atmospheric circulation. The equations are derived from the depth-integrated Navier-Stokes equation under the assumption that the horizontal length scale is much larger than the vertical length scale. The derivation is based on the equations of conservation of mass and conservation of momentum.

The partial differential equation for linear Shallow Water equation in one dimensional is:

$$\frac {\partial \zeta}{\partial t}+\frac{\partial U}{\partial x}=0  ...(1)$$

$$\frac{\partial U}{\partial t}+gH \frac{\partial \zeta}{\partial t}=0  ...(2)$$

$$\zeta$$: the water surface elevation.

$$U=\int_{-h}^{\zeta}u dz$$: the integral of the velocity with respect to depth.

$$H=h+\zeta$$

If the nonlinear advection terms are considered, then the Equations $$(2)$$goes to:

$$\frac{\partial U}{\partial t}+gH \frac{\partial \zeta}{\partial t}+U\frac{\partial U}{\partial x}=0$$

To solve the shallow water equations, the analytical solutions are usually complicated to obtain, so various numerical schemes have been developed to get the numerical solutions in multi-dimensions. However, the numerical solutions sometimes have to be compared with analytical to evaluate the accuracy of numerical schemes. Therefore, the derivation of the analytical solutions for shallow water equations is always necessary and worthwhile.

Project Description
This project is focused on a tide propagating in a semi-closed rectangular basin, which is shown in Figure 1.



Then length of the basin is 100 km, and the width is 50 km. Three of the boundaries are closed except the entrance to the basin, and a standing wave could be formed in this basin. The tide is a semidiurnal tide with an amplitude 0.5 m and period 12 hours.

 Remember to put a license in each image to prevent their removal. None of the uploaded images had a license; some may be removed soon. Egm6322.s09 17:27, 1 May 2009 (UTC)  The licenses have all been added. Egm6322.s09.three.liu 19:18, 1 May 2009 (UTC)

Project Goals
The goal of this project is to firstly derive the analytical solution for the linear Shallow Water Equation in one dimensional by the method of separation of variables, and then, apply the numerical model to solve the Shallow Water Equation in order to calculate the numerical solution. Finally, compare the analytical and numerical solutions.

Governing Equations
The governing equations for the tidal propagation problem are Shallow Water Equations.

The 1-D linear Shallow Water Equations are:

$$\frac {\partial \zeta}{\partial t}+\frac{\partial U}{\partial x}=0  ...(1)$$

$$\frac{\partial U}{\partial t}+gH \frac{\partial \zeta}{\partial t}=0  ...(2)$$

$$\zeta$$: the water surface elevation.

$$U=\int_{-h}^{\zeta}u dz$$: the integral of the velocity with respect to depth.

$$H=h+\zeta$$

Boundary Condition
The Boundary condition for this problem is:

$$U(x=0,t)=0$$ : the left boundary is closed;

$$\zeta(x=l,t)=a sin\omega t$$ : tide goes into the basin from the entrance;

$$\frac{\partial \zeta}{\partial x}|_{x=0}=0$$ : the elevation gradient in the left boundary is zero.

Initial Condition
The initial conditions are:

$$\zeta(x,t=0)=0$$: there is no water surface elevation in the basin before the tide comes in;

$$U(x,t=0)=0$$: water is stationary before the tide comes in.

Classification of the PDE
Theory Background 1: the Classification of PDE

In the general form of PDE:

$$Au_{xx}+2Bu_{xy}+Cu_{yy}+Du_{x}+Eu_{y}+Fu+G=0$$

it can be classified by the value of $$ac-b^2$$.

$$AC-B^2=det \underline A$$

$$AC-B^2<0 \Rightarrow$$ PDE is hyperbolic.

$$AC-B^2=0\Rightarrow$$ PDE is parabolic.

$$AC-B^2>0\Rightarrow$$ PDE is elliptic.

Combining Equation $$(1)$$ and $$(2)$$ gives the 2nd order wave equation in 1-D:

$$\frac{\partial^2 \zeta}{\partial t^2}=c^2\frac{\partial^2 \zeta}{\partial x^2}$$

$$c=\sqrt{gH}$$ which is the phase velocity.

Therefore, according to the classification of the PDE, for the 2nd order 1-D wave equation,

$$A=c^2;B=0;C=-1;D=0;E=0;F=0;G=0.$$

Then $$AC-B^2=-c^2<0$$,

So it is a hyperbolic PDE.

If the Shallow Water Equations are transformed into curvilinear coordinates which is widely used by ocean models, the characteristic of the hyperbolic Shallow Water Equation still remains the same, since the classification of PDE does not change with the coordinate transformation.

So it can be concluded that the Shallow Water Equations are a set of hyperbolic partial differential equations.

Theory Background 2: the Derivation of the Wave Equation (Hyperbolic) The Derivation of Wave Equation

Consider an infinite small segment of an elastic string shown in the picture, a small amplitude vibration of this string conform to the wave equation.

In the picture,

$$u(x,t)$$ is the vertical displacement of the string from the $$x$$ axis.

$$\theta(x,t$$) is the angle between the string and a horizontal line.

$$T(x,t)$$ is the tension in the string.

$$\rho(x)$$ is the mass density of the string.

$$F(x,t)$$ is the external force on the string.

Applying Newton's law in the vertical direction:

$$\rho(x) \sqrt{\Delta x^2+\Delta u^2} \frac {\partial^2 u}{\partial t^2}(x,t)=T(x+\Delta x,t)sin\theta(x+\Delta x,t)-T(x,t)sin\theta(x,t)+F(x,t)\Delta x$$

Divide $$\Delta x$$ on both RHS and LHS and take the limit as $$\Delta x \rightarrow 0$$ gives:

$$\rho(x)\sqrt{1+(\frac {\partial u}{\partial x})^2}\frac {\partial^2 u}{\partial t^2}(x,t)$$

$$=\frac {\partial}{\partial x}[T(x,t)sin\theta(x,t)]+F(x,t)$$

$$=\frac {\partial T}{\partial x}(x,t)sin\theta(x,t)+T(x,t)cos\theta(x,t)\frac{\partial \theta}{\partial x}{x,t}+F(x,t)$$ ...(1)

From the geometric relation in the picture,

$$sin\theta(x,t)=\frac{\frac{\partial u}{\partial x}(x,t)}{\sqrt{1+(\frac{\partial u}{\partial x}(x,t))^2}}$$

$$cos\theta(x,t)=\frac{1}{\sqrt{1+(\frac{\partial u}{\partial x}(x,t))^2}}$$

$$\theta(x,t)=tan^{-1}\frac{\partial u}{\partial x}(x,t)$$

$$ \frac{\partial \theta}{\partial x}(x,t)=\frac{\frac{\partial^2 u}{\partial x^2}(x,t)}{1+(\frac{\partial u}{\partial x}(x,t))^2}$$

Because $$\left | \theta(x,t) \right \vert<<1$$,

$$\sqrt{1+(\frac{\partial u}{\partial x})^2}\approx 1$$

$$ sin\theta(x,t)\approx\frac{\partial u}{\partial x}(x,t)$$

$$cos\theta(x,t)\approx1$$

$$\frac{\partial \theta}{\partial x}(x,t)\approx\frac{\partial^2 u}{\partial x^2}(x,t)$$

Substituting these into equation (1) gives

$$\rho(x) \frac{\partial^2 u}{\partial t^2}(x,t)=\frac{\partial T}{\partial x}(x,t)\frac{\partial u}{\partial x}(x,t)+T(x,t)\frac{\partial^2 u}{\partial x^2}+F(x,t)$$

Assume that the string only has transverse vibration, which means the horizontal force on it is zero.

So $$T(x+\Delta x,t)cos\theta(x+\Delta x,t)-T(x,t)cos\theta(x,t)=0$$

As $$\Delta x \rightarrow 0$$

$$\frac{\partial}{\partial x}[T(x,t)cos\theta(x,t)]=0$$

For small vertical vibration,$$cos\theta \rightarrow 1$$, and$$\frac{\partial T}{\partial x}(x,t)\rightarrow 0$$

Therefore, Equation (1) can be further simplified to:

$$\rho(x)\frac{\partial^2 u}{\partial t^2}(x,t)=T(t)\frac{\partial^2 u}{\partial x^2}(x,t)+F(x,t)$$

Assume the density is constant and no external forces, the second order 1-D wave equation finally goes to

$$\frac{\partial^2 u}{\partial t^2}(x,t)=c^2 \frac{\partial^2 u}{\partial x^2}(x,t), c=\sqrt{\frac{T}{\rho}}$$

In water waves, $$c$$ is the phase speed.

Derivation of the Analytical Solution to Shallow Water Equation
Based on the method of Separation of Variables, we may seek a solution with the following form :

$$\zeta(x,t)=\frac{a cos(kx)}{cos kl}sin \omega t+\sum_{k=0}^{\infty} A_n sin [k_n (l-x)] sin \omega_n t...(3)$$

$$k=\frac{\omega}{c}$$ which is the wave number.

Applying the boundary condition: $$\frac{\partial \zeta}{\partial x}|_{x=0}=0$$,

$$ \sum_{k=0}^{\infty} A_n cos k_n l sin \omega_n t=0$$

$$\Rightarrow cos k_n l=0$$

$$\Rightarrow k_n=\frac{(2n+1)\pi}{2l}$$

In Equation $$(3)$$ $$\omega_n=k_n c$$

Applying the initial condition: $$U(x,t=0)=0$$

$$\sum_{k=0}^{\infty} A_n sin [k_n (l-x)]=\frac{-a\omega cos(kx)}{coskl}$$

$$\Rightarrow A_n=\frac{-2a\omega}{l\omega_n}\int_{0}^{l}\frac{cos (kx)}{cos kl}sin [k_n (l-x)]dx=\frac{-2a\omega}{lc({k_n}^2-k^2)}$$

Therefore, the analytical solutions to the 1-D linear Shallow Water Equations are:





$$\zeta(x,t)=\frac{acos(kx)}{cos(kl)}sin\omega t+\sum_{k=0}^{\infty}\frac{-2a\omega}{lc({k_n}^2-k^2)}sin[k_n (l-x)]sin\omega_nt ...(4)$$

$$U(x,t)=\frac{accos(kx)}{cos(kl)}cos\omega t+\sum_{k=0}^{\infty}\frac{-a\omega}{l({k_n}^2-k^2)}cos[k_n (l-x)]cos\omega_nt ...(5)$$

The solutions can be divided into two parts:

the steady part:

$$\zeta_s=\frac{acos(kx)}{cos(kl)}sin\omega t$$

$$U_s=\frac{accos(kx)}{cos(kl)}cos\omega t$$

the transient part:

$$\zeta_t=\sum_{k=0}^{\infty}\frac{-2a\omega}{lc({k_n}^2-k^2)}sin[k_n (l-x)]sin\omega_nt$$

$$U_t=\sum_{k=0}^{\infty}\frac{-a\omega}{l({k_n}^2-k^2)}cos[k_n (l-x)]cos\omega_nt$$

The plot for the analytical solutions $$(4)$$ and $$(5)$$ at the location x=50 km is shown in Figure 2 and Figure 3.

Description of the Numerical Model
This model developed by the author uses conjugate gradient method to solve the Shallow Water Equations in 2-D. Under the condition that the bottom is flat, and the open boundary condition $$\zeta(x=l_x,y,t)=asin(wt)$$ is applied, the solutions can be compared with the 1-D analytical solutions.

 Nice animation. Egm6322.s09 17:27, 1 May 2009 (UTC)

The Shallow Water Equations in 2-D is:

$$\frac{\partial \zeta}{\partial t}+\frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}=0...(6)$$

$$\frac{\partial U}{\partial t}+gH\frac{\partial \zeta}{\partial x}=0...(7)$$



$$\frac{\partial V}{\partial t}+gH\frac{\partial \zeta}{\partial y}=0...(8)$$

A simple space and time staggered finite difference mesh (50 by 25) is used to discretize the domain as depicted in Figure 1. The governing equations are discretize by the Upwind Method and solved implicitly to guarantee stability.

Applying the upwind method, the finite difference equations for the governing equations are:

$$\frac{\zeta_{i,j}^{n+1}-\zeta_{i,j}^n}{\Delta t}+\frac{U_{i,j}^{n+1}-U_{i-1,j}^{n+1}}{\Delta x}+\frac{V_{i,j}^{n+1}-V_{i,j-1}^{n+1}}{\Delta y}=0 ...(9)$$

$$\frac{U_{i,j}^{n+1}-U_{i,j}^n}{\Delta t}+gH\frac{\zeta_{i+1,j}^{n+1}-\zeta_{i,j}^{n+1}}{\Delta x}=0...(10)$$

$$\frac{V_{i,j}^{n+1}-V_{i,j}^n}{\Delta t}+gH\frac{\zeta_{i+1,j}^{n+1}-\zeta_{i,j}^{n+1}}{\Delta x}=0...(11)$$

Substituting Equation $$(10)$$ and $$(11)$$ into $$(9)$$ gives:

$$[1+2gH(\frac{\Delta t}{\Delta x})^2+2gH(\frac{\Delta t}{\Delta y})^2]\zeta_{i,j}^{n+1}-gH(\frac{\Delta t}{\Delta x})^2\zeta_{i-1,j}^{n+1}-gH(\frac{\Delta t}{\Delta x})^2\zeta_{i+1,j}^{n+1}-gH(\frac{\Delta t}{\Delta y})^2\zeta_{i,j-1}^{n+1}-gH(\frac{\Delta t}{\Delta y})^2\zeta_{i,j+1}^{n+1}$$

$$=\zeta_{i,j}^n-\frac{\Delta t}{\Delta x}(U_{i,j}^n-U_{i-1,j}^n)-\frac{\Delta t}{\Delta y}(V_{i,j}^n-V_{i,j-1}^n)...(12)$$

The open boundary condition is:

$$\zeta(i=nx,j,t)=asin(wt)$$

The Equation $$(12)$$ generates a sparse matrix which is mainly five diagonal. The Conjugate Gradient Method is an efficient way to solve the sparse matrix problem. Applying the conjugate gradient method, the values for $$\zeta_{i,j}^{n+1}$$ are all solved. Substituting the values of $$\zeta_{i,j}$$ into Equation $$(10)$$ and $$(11)$$ gives the values for $$U$$ and $$V$$. For a flat bottom, the velocity $$V$$ should be zero.

The propagation of the tide in a semi-closed basin simulated by this model is shown in Figure 4(1) (in this animation, the amplitude of the tide is 1 and the period is 1000 seconds), and the semidiurnal tide with a amplitude 0.5 m is shown in Figure 4(2). It is seen that a standing wave is formed in this basin.

Comparison between the analytical and numerical solutions
The comparison between the analytical and numerical solutions to the water surface elevation $$\zeta$$ and the longitudinal velocity $$U$$ at the location x=50 km is shown in Figure 5 and Figure 6.





From the comparison for the water surface elevation in Figure 5, it is seen that the analytical and numerical solutions are fit well. However, the numerical model is based on implicit method which may cause some numerical dissipation. So after 80 hours, there are some difference between the analytical and numerical solutions.

From the comparison for the longitudinal velocity in Figure 6, it is seen that the two solutions are generally fit well. The numerical model is developed for 2-D problem, and for a flat bottom, the velocity in y direction should be zero. However, in fact, the velocity in y direction is not zero, although very small (10-3 scale), due to some numerical errors. Hence, besides the implicit method that may cause dissipation, the non-zero velocity in y direction may also cause some errors to the longitudinal velocities.

Conclusion
The analytical and numerical solutions to a hyperbolic partial differential equation: linear Shallow Water Equation for tidal propagation in a semi-closed basin are derived and compared. The analytical solution consists of a steady part and a transient part. The numerical solution is obtained from an implicit numerical model which applies upwind method and conjugate gradient method. The comparison shows that the analytical and numerical solutions are generally in good agreement, and the analytical solution derived from the method of separation of variables can be used to evaluate the accuracy of the numerical solution.