Flow Beneath An Embankment Using Laplace & PDE

Sunday. May 12, 2019

This project implements the Laplace Equation with Partial Different Equations (PDE) to determine the flow beneath an embankment.

Method

1. Discretization

The solution domain must be broken down into a set of representative points (discretized). This will allow the domain below the embankment to be visualized as a mesh.

2. Set of Governing Equations

Laplace equation is used to represent the physical process through the embankment, expressed as:

$\nabla^2 H=\frac{\partial^2 H}{\partial x^2}+\frac{\partial^2 H}{\partial y^2}$

3. Numerical Representation of Governing Equations

Using partial derivatives that are accurate to second order in the mesh spacing. The 2nd order finite difference equation is:

$-−2(h^2+k^2)H_{i,j}+k^2 H_{i-1,j}+k^2 H_{i+1,j}+h^2 H_{i,j-1}+h^2 H_{i,j+1}=0$

where $H_{i,j}$ is the head at any internal grid point i, j.

4. Boundary Conditions

The domain is bounded and at each boundary there are potential effects of the boundary on the interior of the domain. The mesh must be constructed so that all boundary points coincide with the boundary of the mesh. An equation must be specified at every boundary point.

At the inflow and outflow boundaries, the equation is trivial and the value of $H_{i,j}$ at these point is simply the total head. For example, at each inflow point:

$H_{i,j}=h_1$

At the other boundaries, a finite difference expression must be constructed by combining the governing equation with another equation stating that the gradient in $H$ normal to the boundary is zero. We need four expressions for the left, right, upper (that is, beneath the embankment) and lower boundaries which are respectively:

$−-2(h^2+k^2)H_{i,j}+2k^2 H_{i+1,j}+h^2 H_{i,j-1}+h^2 H_{i,j+1}=0$
$−-2(h^2+k^2)H_{i,j}+2k^2 H_{i-1,j}+h^2 H_{i,j-1}+h^2 H_{i,j+1}=0$
$-−2(h^2+k^2)H_{i,j}+k^2 H_{i-1,j}+k^2 H_{i+1,j}+2h^2 H_{i,j-1}=0$
$−-2(h^2+k^2)H_{i,j}+k^2 H_{i-1,j}+k^2 H_{i+1,j}+2h^2 H_{i,j+1}=0$

There are two special boundary condition points: the lower left grid point and the lower right grid point. The applicable boundary conditions at these are respectively:

$-−2(h^2+k^2)H_{i,j}+2k^2 H_{i+1,j}+2h^2 H_{i,j+1}=0$
$−-2(h^2+k^2)H_{i,j}+2k^2 H_{i-1,j}+2h^2 H_{i,j+1}=0$

Solution Method for Numerical Equations

Constructs the full set of linear equations that incorporates an equation for each boundary and internal point. This will form a matrix that maps each point i,j to a row and column number in the matrix.

The solution will be:

$A \hat H =\hat b$

where $A$ is the coefficient matrix, $\hat H$ is the unknown column vector of head values at each grid point and $\hat b$ is the column vector containing the constant value from the equation corresponding to each grid point.

Thus in order to solve for $\hat H$ , the equation is rearranged to $\hat H=A^{-1}\hat b$ </p>

Once $\hat H$ is solved, the velocity components can be computed using 2nd order finite differences again as:

$u=-K \frac{\partial^2 H}{\partial x^2}=\frac{H_{i-1,j}+H_{i+1,j}}{2h}$
$v=-K \frac{\partial^2 H}{\partial y^2}=\frac{H_{i,j-1}+H_{i,j+1}}{2k}$

The total volume flux of water beneath the embankment can be determined by integrating the velocity normal to a defined surface. The total volume flux per unit length along the crest of the embankment $q$ can be computed as:

$q=\oint \vec{u} \cdot \hat n \partial A$

$u$ is the velocity vector at the defined surface, $\hat n$ is the local unit normal vector of the defined surface and $\partial $$A$$ is an elemental area of the defined surface.