Skip to content

Commit fb6ec8c

Browse files
jerabekjakjerabekjak
authored
Merge pull request #12 from storm-fsv-cvut/numerical_solution
Numerical solution
2 parents 5bd928e + df78b08 commit fb6ec8c

File tree

2 files changed

+234
-165
lines changed

2 files changed

+234
-165
lines changed

‎reference_manual/text_en/hillhyd.tex

-164
Original file line numberDiff line numberDiff line change
@@ -275,170 +275,6 @@
275275
changes during the calculation.
276276

277277

278-
\subsection{Implicit /Explicit approach}
279-
280-
\subsubsection{Explicit}
281-
The time derivative in Equation \ref{eq:bilance} is calculated using an
282-
explicit method. The computation is therefore sensitive to the size of the time
283-
step. The size of the time step is controlled by the Courant criterion, which
284-
needs to be kept below the theoretical maximum value of 1, while the maximum
285-
value in practise is lower than 1
286-
\cite{zhang1989modeling, esteves2000overland}.
287-
288-
289-
The sheet flow water level of the next time t + 1 step in Equation
290-
\ref{eq:bilance} which incorporates the sum \ref{eq:routing} is calculated with the
291-
explicit time discretisation scheme for cell i as:
292-
\begin{equation}
293-
h_{i,t} =h_{i,t-1} + \mathrm{d}t (es_{i,t-1} + \sum_j^m q^{out}_{j,t-1}-
294-
inf_{i,t-1} - q^{out}_{i,t-1}),
295-
\label{eq:bilanceexpl}
296-
\end{equation}
297-
298-
299-
After inserting the equations~(\ref{eq:infiltration})
300-
and~(\ref{eq:powerlaw2}) into equation~(\ref{eq:bilanceexpl}) the
301-
final explicit form balance equation can be written as:
302-
\begin{dmath}
303-
h_{i,t} = h_{i,t-1} +
304-
es_{i,t-1} + exf_{i,t-1} + \sum_{j}^{inflows}
305-
1/n_{sheet,j}I_j^{y_j} h_{j,t-1}^{b_j} -
306-
1/n_{sheet,i}I_i^{y_i}h_{i,t-1}^{b_i} - (
307-
\frac{1}{2}S_it_1^{-1/2}+K_{s,i}).
308-
\end{dmath}
309-
310-
311-
\subsubsection{Implicit}
312-
Alternatively the equation (\ref{eq:bilance}) solved using implicit
313-
scheme. The right hand side can be expressed in for time $t$ as:
314-
315-
$$
316-
\frac{h_{i,t} - h_{i,t-1} }{\mathrm{d} t} = \left(es_{i,t} +
317-
\sum_j^{inflows} q^{in}_{j,t} - inf_{i,t} - q^{out}_{i,t}\right).
318-
$$
319-
320-
After inserting the power-law equation (\ref{eq:powerlaw}) and several rearrangements the formula above can be written
321-
as follows:
322-
323-
324-
325-
$$
326-
h_{i,t} + \mathrm{d} t a_ih^{b_{i}-1}_{i,t} h_{i,t} - \mathrm{d} t \sum_j^{inflows}
327-
a_jh^{b_{j}-1}_{j,t} h_{j,t} = h_{i,t-1} + \mathrm{d} t es_{i,t} - \mathrm{d}
328-
t inf_{i,t}.
329-
$$
330-
Here the known part of the balance equation (rainfall and infiltration) are at the right hand side. To
331-
extract the unknown $h$ form the power law a following rearrangement was used:
332-
333-
334-
335-
\begin{equation}
336-
(1+\mathrm{d} t a_ih^{b_{i}-1}_{i,t})h_{i,t} - \mathrm{d} t \sum_j^{inflows}
337-
a_jh^{b_{j}-1}_{j,t} h_{j,t} = h_{i,t-1} + \mathrm{d} t es_{i,t} -
338-
\mathrm{d} t inf_{i,t}
339-
\label{eq:impl}
340-
\end{equation}
341-
342-
From the equation (\ref{eq:impl}) the system of linear equations can be constructed.
343-
344-
As example, lets show the equation (\ref{eq:impl}) for
345-
$
346-
i=5\ \mathrm{and}\ inflows\in\{7,8,9\}
347-
$:
348-
349-
350-
\begin{dmath}
351-
(1+\Delta T a_5h^{b_{5}-1}_{5,t})h_{5,t} - \Delta T a_7h^{b_{7}-1}_{7,t} h_{7,t} + \Delta T a_8h^{b_{8}-1}_{8,t} h_{8,t} + \Delta T a_9 h^{b_{9}-1}_{9,t} h_{9,t} = h_{5,t-1} + \Delta T es_{5,t} - \Delta T inf_{5,t}.
352-
\end{dmath}
353-
354-
In matrix form the above equation can be written as:
355-
356-
357-
358-
\begin{dmath}
359-
\begin{bmatrix}
360-
\hdots& & & & & & \\
361-
\hdots & 1+\Delta T a_5h^{b_{5}-1}_{5,t} & \hdots & \Delta T a_7h^{b_{7}-1}_{7,t} & \Delta T a_8h^{b_{8}-1}_{8,t} & \Delta T a_9 h^{b_{9}-1}_{9,t} & \hdots \\
362-
\hdots& & & & & & \\
363-
\hdots& && & & & \\
364-
\hdots& && & & & \\
365-
\hdots& && & & & \\
366-
367-
\end{bmatrix}
368-
%
369-
\begin{bmatrix}
370-
\vdots \\
371-
h_{5,t} \\
372-
\vdots \\
373-
h_{7,t} \\
374-
h_{8,t} \\
375-
h_{9,t} \\
376-
\vdots \\
377-
%
378-
\end{bmatrix}
379-
=
380-
\begin{bmatrix}
381-
\vdots \\
382-
h_{5,t-1} + \Delta T es_{5,t} - \Delta T inf_{5,t} \\
383-
\vdots \\
384-
\vdots \\
385-
\vdots \\
386-
\vdots \\
387-
\vdots \\
388-
%
389-
\end{bmatrix}
390-
\label{eq:sys}
391-
\end{dmath}
392-
393-
The system (\ref{eq:sys}) is solved using \texttt{scipy} package method \texttt{root}
394-
using \texttt{krigov} methof that shown the best convergence.
395-
396-
397-
\paragraph{Implicit solution with rill flow} To calculated rill flow the water level in rill and needs to be calculated as
398-
\begin{equation}
399-
h_{rill} = \text{max}(h-h_{crit},0),
400-
\label{eq:hrill2}
401-
\end{equation}
402-
where
403-
\begin{equation}
404-
h_{sheet} = \text{min}(h,h_{crit}).
405-
\label{eq:hsheet2}
406-
\end{equation}
407-
408-
Mannings formula is used to calculated the steady state flow in
409-
the rill assuming rill is a channel (description in section \ref{sec:rill}.
410-
$$
411-
q_{rill} = 1/n R(h_{rill})^{2/3} i^{1/2}
412-
$$
413-
414-
415-
The balance equation becomes
416-
\begin{dmath}
417-
\frac{h_{i,t} - h_{i,t-1} }{\Delta T} =
418-
es_{i,t} + \sum_j^{inflows} a_jh^{b_{j}}_{sh,j,t} + \sum_k^{inflows} 1/n_k R_k(h_{rl,k,t})^{2/3} i_k^{1/2} - inf_{i,t} - a_ih^{b_{i}}_{sh,i,t} - 1/n R(h_{rl,i,t})^{2/3} i^{1/2},
419-
\end{dmath}
420-
where $h_{sheet}$ and $h_{rill}$ are defined as (\ref{eq:hsheet2}) and (\ref{eq:hrill2}).
421-
422-
This formula is rearranged and, for practical purposes, the linear equation is the system is calculated with condition \\
423-
for $h<=h_{crit}$ as
424-
\begin{equation}
425-
\left(\frac{1}{\Delta T}+a_ih^{b_{i}-1}_{i,t}\right)h_{i,t} - \sum_j^{inflows} a_jh^{b_{j}-1}_{j,t} h_{j,t} = \frac{h_{i,t-1}}{\Delta} + es_{i,t} - inf_{i,t}
426-
\end{equation}
427-
%
428-
%
429-
and for $h>h_{crit}$ as
430-
\begin{multline}
431-
\left(\frac{1}{\Delta T}
432-
+ 1/n_i R_i(h_{i,t}-h_{crit,i})^{2/3} i_i^{1/2} \frac{1}{h_{i,t}}\right)h_{i,t} \\
433-
- \sum_k^{inflows} \left( 1/n_k R_k(h_{k,t}-h_{crit,k})^{2/3} i_k^{1/2} \frac{1}{h_{k,t}}\right)h_{k,t}
434-
= \\
435-
\frac{h_{i,t-1} }{\Delta T}
436-
+ es_{i,t}
437-
+ \sum_j^{inflows} a_j h^{b_{j}}_{crit,j}
438-
- inf_{i,t}
439-
- a_i h^{b_{i}}_{crit,i}
440-
\end{multline}
441-
The later system of equations is constructed in similar fashion as the case of (\ref{eq:sys}) and solved with the same solver.
442278

443279

444280
\FloatBarrier

0 commit comments

Comments
 (0)