User:EGM6322.S09.TIAN/report

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

=Numerical Modeling of Shallow Water Equation=

Summary
A one-dimensional model of the shallow water equation is developed. The model is based on the one-dimensional continuity equation and momentum equation. A staggered-grid and theta method scheme is utilized in this finite-difference model. The model is used to test two PDE cases, and the results obtained from these numerical experiments are then compared with their analytical solutions. Results indicate that the model was able to give a relatively precise simulation of the solution to the one dimensional shallow water equation. The numerical simulation is affected by initial conditions and boundary conditions while formulating the model, but the results show promise for usefulness of the numerical simulation with future exploration.

Introduction
Currently, partial differential equations (PDEs) are widely used in research of coastal engineering. The fundamental equations of hydrodynamic processes, such as the continuity equation and momentum equation, are PDEs. Analytical solutions to some PDEs, which have simple initial conditions and boundary conditions, are given by mathematicians. However, some of the simple PDEs with complicated initial and boundary conditions are difficult to solve. And even worse, some PDEs are so complex that mathematicians are not able to develop analytical solutions for them. Fortunately, with the development of computer technology, numerical simulation methods are applied to give numerical solutions to PDEs.

The goal of the research is to develop an accurate and efficient model of the one-dimensional shallow water equation, which can be readily expanded into a multi-dimensional model. In the following, the formulation of this one-dimensional hydrodynamics model is described first. The model is then used to examine two shallow-water equation cases with different initial conditions. Comparing them with their analytical solutions, the numerical results obtained from these tests would be used to evaluate whether the one-dimensional model is sufficiently precise to simulate shallow water hydrodynamics process.

On the basis of the material from a paper, “Numerical Modeling of Tidal Hydrodynamics and Salinity Transport in the Indian River Lagoon”, which is given by Y. Peter Sheng, S. Peene, and Y. M. Liu (June 6, 1989), the code of the model is written in FORTRAN. Moreover, the PDE part of this report is based on two reference books as follows, “Partial Differential Equations for Scientists and Engineers”, “Partial Differential Equations in Mechanics 1”. The form of the shallow-water equation was derived in “Introduction to Physical Oceanography”. Last but not least, all the figures used in this report were plotted by MATLAB. The animations were generated by Macromedia.

Basic Partial Differential Equations
The shallow-water equation is given as follows:

$$\frac {\partial \eta}{\partial t}+ \frac {\partial}{\partial x}[u(h+ \eta)]=0 (1)$$

$$\frac {\partial u}{\partial t}+ u \frac {\partial u}{\partial x}=-g \frac {\partial \eta}{\partial x}(2)$$

Neglecting the non-linear terms in (2), assuming a flat bottomed basin $$(h=constant)$$, and assuming $$ \eta \ll h $$ , the one-dimensional equations of motion for the free surface elevation $$\eta$$, vertically and laterally averaged velocity $$u$$ are simplified to:

$$\frac {\partial \eta}{\partial t}+h \frac {\partial u}{\partial x}=0 (3)$$

$$\frac {\partial u}{\partial t}+g \frac {\partial \eta}{\partial x}=0 (4)$$

By differentiating (3) with respect to $$t$$ and (4) with respect to $$x$$ and subtracting the two results we obtain:

$$c^2 \frac {\partial^2 \eta}{\partial x^2}-\frac {\partial^2 \eta}{\partial t^2}=0(5)$$

Or by reversing the order of differentiation we obtain:

$$c^2 \frac {\partial^2 u}{\partial x^2}-\frac {\partial^2 u}{\partial t^2}=0(6)$$

where $$c= \sqrt {gh}$$.

According to the wave theory, $$c \equiv \sqrt {gh}$$is the wave celerity in shallow water.

Applying the method of classifying PDEs, which we have learned in the class, to this case,

$$B^2-AC=c^2>0$$

we learn that this is a hyperbolic PDE.

Moreover, we know that (5) or (6) is the one-dimensional wave equation. From the class notes, we know that the solution of one-dimensional wave equation with initial condition was given as d'Alembert solution, which we will give in the following section.

d'Alembert solution

$$PDE:$$

$$c^2 \frac {\partial^2 u}{\partial x^2}-\frac {\partial^2 u}{\partial t^2}=0$$, $$- \infty <x<\infty $$; $$0 <t<\infty $$

Initial Conditions:

$$\begin{cases} u(x,0)=f(x) \\ u_t(x,0)=g(0) \end{cases} $$, $$- \infty <x<\infty $$.

The d'Alembert solution is:

$$u(x,t)= \frac {1}{2} [f(x-ct)+f(x+ct)]+ \frac {1}{2c}\int_{x-ct}^{x+ct} g(\xi)\, d \xi$$

Numerical scheme


In the numerical model, we use a stagger grid. The staggered grid means that the water surface elevation $$\eta$$ (p points in Figure 1)and velocity points $$u$$ are staggered points on an axis.

Theta time-stepping method is applied to this model. The following scheme is explicit time-stepping when theta is equal to 0, implicit when theta is equal to 1, and semi-implicit when theta is equal to 0.55.

The discrete form of (3) and (4) are given as follows:

$$\frac {\eta^{n+1}_i-\eta^{n}_i}{\Delta t}+(1-\theta) h\frac {u^{n}_i-u^{n}_{i-1}}{\Delta t}+\theta h\frac {u^{n+1}_i-u^{n+1}_{i-1}}{\Delta t}=0 (7)$$

$$\frac {u^{n+1}_i-u^{n}_i}{\Delta t}+(1-\theta)g\frac {\eta^{n}_{i+1}-u^{n}_{i}}{\Delta t}+\theta g\frac {\eta^{n+1}_{i+1}-\eta^{n+1}_{i}}{\Delta t}=0 (8)$$

Based on the equations above, a one-dimensional shallow water model is developed, and then examined by two cases.

First Test
As the first step, the one-dimensional model is used to examine the following case :

$$ \frac {\partial^2 \eta}{\partial x^2}-\frac {\partial^2 \eta}{\partial t^2}=0$$

Initial conditions:

$$f(x)=\begin{cases} \eta(x,0)=1;31\le x \le 40 \\ \eta(x,0)=0;otherwise \end{cases} $$,

$$g(x)=\eta_t(x,0)=0; -\infty  44 \end{cases} $$



when $$t=5$$

$$f(x)=\begin{cases} \eta(x,4)=0;x < 26 \\ \eta(x,4)=\frac {1}{2};26\le x < 35 \\ \eta(x,4)=\frac {1}{2};35< x \le 45 \\ \eta(x,4)=0;x > 45 \end{cases} $$



when $$t=15$$

$$f(x)=\begin{cases} \eta(x,4)=0;x < 16 \\ \eta(x,4)=\frac {1}{2};16\le x < 25 \\ \eta(x,4)= 0; x=36 \\ \eta(x,4)=\frac {1}{2};46< x \le 55 \\ \eta(x,4)=0;x > 55 \end{cases} $$

The above solutions， which are the water surface elevation( $$\eta$$ )at four certain moments, are plotted on the right hand side as Figure2.

In numerical experiment, I set $$\Delta x =1, \Delta t =1, c=1, \theta =0. $$

The numerical results, which are the water surface elevation( $$\eta$$ )at four certain moments, are plotted on the right hand side as Figure3.

And an animation to describe the process of this case are given as Figure 4.

 If you had used matlab, you could also produce an animation with matlab. The background of this animation is a bit grainy. Egm6322.s09 18:49, 1 May 2009 (UTC)

 I'm sorry for the animation, actually I produce the movie by MATLAB, and as I transform it from movie format to gif format by Macromedia, the grainy thing happened. I notice that you appreciate Liu's animation. However, he developed his animation by Tecplot,which is able to produce better animation than MATLAB. --EGM6322.S09.TIAN 20:34, 1 May 2009 (UTC)Egm6322.s09 18:49, 1 May 2009 (UTC)

The physical meaning of this case is that, in the beginning (t=0), the water surface between point 31 and 40 are 1 meter higher than average water surface. According to the analytical solution, as t increases, the initial surface elevation would be split into two half-heighted waves which move in an opposite direction. Comparing the numerical results with the analytical results, we find out that the numerical results showed the general idea of this case. However, we notice that the amplitudes of the numerical waves are much lower than they are predicted. In addition, the shapes of these two new waves are not consistent with their original shapes. We will discuss this phenomenon in the next section.

Second Test
$$ \frac {\partial^2 \eta}{\partial x^2}-\frac {\partial^2 \eta}{\partial t^2}=0$$

Initial Conditions:

$$\begin{cases} \eta(x,0)=sin[(x-55) \pi /10];56\le x \le 65 \\ \eta(x,0)=0;otherwise \end{cases} $$,

$$\eta_t(x,0)=0; -\infty  75 \end{cases} $$

when $$t=15$$

$$f(x)=\begin{cases} \eta(x,10)=0;x < 41 \\ \eta(x,10)=\frac {1}{2}sin \frac {(x-40)\pi}{10};41\le x \le 50 \\ \eta(x,10)=0;50 80 \end{cases} $$

when $$t=20$$

$$f(x)=\begin{cases} \eta(x,20)=0;x < 36 \\ \eta(x,20)=\frac {1}{2}sin \frac {(x-35)\pi}{10};36\le x \le 45 \\ \eta(x,20)=0;45 85 \end{cases} $$

The above solutions, which are the water surface elevation $$(\eta)$$ at four certain moments, are plotted on the right hand side as Figure5.

In numerical experiment, I set $$\Delta x =1, \Delta t =1, c=1, \theta =0. $$ first. The results are plotted in Figure 6. We find that the shapes of the results are similar to the original one. Nonetheless, we face the same problem as we have in the first case, that is the amplitudes of the two waves are less then half of the original wave height as they propagates.

Then I set $$\Delta x =1, \Delta t =1, c=1, \theta =0.55 $$. The results are showed in Figure 7.

As far as I know, $$\theta$$ is related to numerical dissipation. Hence the great damping we observed in Figure 6 is due to numerical dissipation, when the numerical scheme is fully-implicit. In addition, we notice that some of the values of u are negative in Figure 7, which is not in accordance with the analytical solutions. The reason is given as follows. When $$\theta=1$$, the numerical scheme is fully-implicit. And when $$\theta=0.55$$, the numerical scheme is semi-implicit. The model would not have as much numerical dissipation in semi-implicit scheme as it has in fully-implicit scheme, hence the results of Figure 7 is more conservative than that of Figure 6. However, it lose stability in semi-implicit scheme. In other words, the fully-implicit scheme is more stable than semi-implicit scheme. By the way, the fully-explicit scheme( $$\theta=0$$ ) is not stable, which would give us unreasonable results. The negative values of u in Figure 7 are caused by instability of semi-implicit scheme. Apparently, the results in Figure 7 are more consistent with the analytical results in Figure 5. Hence we usually apply semi-implicit scheme, that is $$\theta=0.55$$, to a numerical experiment. And an animation to describe the process of this case when $$\theta=0.55$$ are given as Figure 8.

Conclusion
A semi-implicit finite difference technique for the shallow-water equations have been presented and analyzed. As seen in the previous part of the report, there appear to be differences between the analytical solutions and the values obtained by the model. However, the results given in Figure 7 simulate the second case accurately. Although the initial conditions are simple and there is no boundary condition, the case is physically meaningful. To some degree, we can imagine this as a dam breaking case. Hence we can also apply boundary conditions, such as open or radiation boundary conditions, to this case. If we do this, the analytical solutions seem to be hard to obtain, while the one-dimensional shallow-water model can still give us reasonable results. The model is meant to be a simplified representation of the actual conditions. To make it more powerful, we are able to take into account wind stresses, bottom friction, and Coriolis forcing. These factors would affect the results greatly, and make the model a more precise simulation of the actual conditions. Moreover, we can expand the one-dimensional model to a two-dimensional one, which would give us more exiting results. The steps described above were unable to be taken during this project because of the time restrain, but if a longer time had been available, more complex solutions could have been accomplished.

List of References
--EGM6322.S09.TIAN 19:30, 30 April 2009 (UTC)