Table of Contents
1 Finite volumes 7 - transient problems and time integration
2 Learning goals
3 The two Eulers
3.1 Time stepping
3.1.1 Forward Euler - explicit
3.1.2 Backward Euler - implicit
4 Transient finite-volume equation for diffusion in porous media
4.1 Your turn
5 Implicit and explicit methods
6 Explicit 1-d stencil for diffusion in porous media
7 Implicit 1-d stencil for diffusion in porous media
7.1 Your turn
7.2 Your turn
8 Summary
8.1 Explicit method
8.2 Implicit method
9 Accuracy, stability and time-step constraints
9.1 Accuracy
9.2 Stability
9.3 Explicit methods are conditionally stable
9.4 Your turn
9.5 Implicit methods are unconditionally stable
10 Reflection
Finite volumes 7 - transient problems and time integration¶
Learning goals¶
Be able to formulate a transient finite-volume method.
Be able to formulate a stencil for the Euler method for a finite - volume approximation.
To recognize the advantages and disadvantages of explicit, implicit and mixed time-stepping schemes.
To be able to code a transient problem in python.
To be able to plot the results in python.
The two Eulers¶
We’ve covered transient ordinary differential equations of the form \begin{align} {dy\over dt} &= f(y)\label{7fv71}\\ \end{align} The idea of the Euler method is to approximate the time derivative as finite difference: \begin{align} {dy\over dt} \approx {y(t+\Delta t) - y(t)\over \Delta t} \label{7fv72}\\ \end{align} Then we can approximate the ordinary differential equation with this algebraic equation: \begin{align} {y(t+\Delta t) - y(t)\over \Delta t} &= f(y)\label{7fv73}\\ \end{align} This finite-difference approximation is accurate if \(\Delta t\) is small. The approach behind most transient problems is to time step with time steps that are sufficiently small to maintain accuracy and stability. We start with the initial condition (knowledge of the dependent variable (\(y\)) at the start of the period of interest (\(t_0\))) and then find \(y\) at time \(t+\Delta t\), …: \begin{align*} \text{Step 1}\quad y(t_0)\quad &\underrightarrow{\Delta t} \quad y(t_0+\Delta t)\\ \text{Step 2}\quad y(t_1)\quad &\underrightarrow{\Delta t} \quad y(t_1+\Delta t)\\ \text{Step 3}\quad y(t_2)\quad &\underrightarrow{\Delta t} \quad y(t_2+\Delta t)\\ \text{Step 4}\quad y(t_3)\quad &\underrightarrow{\Delta t} \quad y(t_3+\Delta t)\\ \end{align*}
etc The value of the dependent variable at the start of the simulation period, usually \(t=0\) is \(y(0)\) and is called the initial condition. The initial condition is part of the problem specification and, along with boundary conditions, is information required before the problem can be solved.
Time stepping¶
In a finite-volume approach with gridblocks in space, the approach is: 1. Start with an initial condition at \(t = t_0\) 2. Determine solution at every node at \(t = t_0 +\Delta t\), where \(\Delta t\) is the time step. 3. Repeat, optionally with a different time step size, and march forward in time, one time step at a time.
The Euler approach then can be written: \begin{align} {y(t+\Delta t)} &= y(t) + \Delta t f(y)\label{7fv74}\\ \end{align} The two Euler methods differ in the time level at which we evaluate \(f(y)\); either the old time level or the future time level:
Forward Euler - explicit¶
\begin{align} {y(t+\Delta t)} &= y(t) + \Delta t \underbrace{f(y(t))}_\text{evaluated at present time}\label{7fv75}\\ \end{align} #### Backward Euler - implicit \begin{align} {y(t+\Delta t)} &= y(t) + \Delta t \underbrace{f(y(t+\Delta t))}_\text{evaluated at future time}\label{7fv76}\\ \end{align}
In this notebook, we’ll develop the transient 1-D finite volume stencil and then show how we can apply explicit and implicit time stepping schemes. It will then be clear why forward Euler is called explicit and backward Euler implicit (I prefer the terms implicit and explicit).
Transient finite-volume equation for diffusion in porous media¶
We’ll develop the equations for 1-D. The 2- and 3-D equations follow the same pattern. Recall from notebook 5.1finite_volume_2.ipynb the following transient mass conservation statement \begin{align}
m_C(t+\Delta t) - m_C(t) & = \left(J_{EC}+J_{WC}\right)\Delta t \\
\end{align} where \(m_C \left[M\right]\) is the mass in control volume C and \(J_{EC}, J_{WC}~\left[M/T\right]\) are the mass fluxes into C. If we take \(\Delta t\) to the other side we get: \begin{align}
{m_C(t+\Delta t) - m_C(t)\over \Delta t} & = J_{EC}+J_{WC} \label{7fv77} \\
\end{align} We developed a stencil for the right side of this equation \(0=J_{EC}+J_{WC}\) for diffusion in porous media in notebooks 6_finite_volume_4.ipynb and 6_finite_volume_5.ipynb. The stencil is (1-D): \begin{align}
&J_{EC}+J_{WC}=\left(D\theta {c_E - c_C \over \Delta x} + D\theta {c_W - c_C \over \Delta x} \right) (\Delta y) (\Delta z) \label{7fv78}\\
\end{align}
Or with terms that multiply \(c_E\), \(c_C\) and \(c_W\) collected together: \begin{align} &J_{EC}+J_{WC} = \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_E - \left(2 D\theta {\Delta y \Delta z \over \Delta x}\right) c_C + \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_W \label{7fv79}\\ \end{align} So, substituting \ref{7fv79} into \ref{7fv77} we can write the transient conservation law for the case of 1-D diffusion in porous media as: \begin{align} &{m_C(t+\Delta t) - m_C(t)\over \Delta t}= \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_E - \left(2 D\theta {\Delta y \Delta z \over \Delta x}\right) c_C + \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_W \label{7fv710}\\ \end{align}
Your turn¶
Show that the dimensions of both sides of transient 1-D diffusion in porous media stencil \ref{7fv710} are the same.
The left side has dimensions \(\left[\text{Left side dimension here}\right]\)
The right side has dimensions \(\left[\text{Right side dimension here}\right]\)
As a last step, we convert the mass in the control volume C, \(m_C\), into a concentration in the control volume \(c_C\).
In porous media, when the pores are fully saturated with water, the volume of water in the gridblock is porosity \(\times\) total volume of the gridblock: \(\theta (\Delta x)(\Delta y)(\Delta z)\), so the mass of solute in the water in the gridblock is the volume of water \(\times\) the concentration in the water: \(c \theta (\Delta x)(\Delta y)(\Delta z)\). So for \(m_C\) we substitute \(c\theta (\Delta x)(\Delta y)(\Delta z)\). \begin{align} &{c_C(t+\Delta t) - c_C(t)\over \Delta t}\theta (\Delta x)(\Delta y)(\Delta z)= \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_E - \left(2 D\theta {\Delta y \Delta z \over \Delta x}\right) c_C + \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_W \label{7fv711}\\ \end{align} We’re not quite done. We have to determine at what time level we compute the right side of this expression.
Implicit and explicit methods¶
From the Euler perspective, the question we ask, is at what time level do we treat the spatial terms. First some notation.
\(c_C^t\) = the concentration in gridblock C at the present time level \(t\). These values are known. At the first time step, \(c_C^{t_0}\) is the initial condition.
\(c_C^{(t+\Delta t)}\) = the concentration in gridblock C at the next time step at time \(t+\Delta t\). These values are not known. Conceptually, we can think of our transient discrete approximation like the simple ODEs in equation \ref{7fv74} that we used to introduce the Euler method: \begin{align} &\quad\quad c_C(t+\Delta t) = c_C(t) + \Delta t f(c_W, c_C, c_E)\\ \nonumber \text{where}&\\ &f(c_W, c_C, c_E)= {1\over\theta (\Delta x)(\Delta y)(\Delta z)}\left\{ \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_E - \left(2 D\theta {\Delta y \Delta z \over \Delta x}\right) c_C + \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_W \right\} \label{7fv712}\\ \end{align} In an explicit method (Forward Euler), we evaluate \(f\) in \ref{7fv712} at the present (known) time level \(t\): \begin{align} &c_C^{(t+\Delta t)} = c_C^t + \Delta t \overbrace{f(c_W^t, c_C^t, c_E^t)}^\text{evaluated at known time}\label{7fv713}\\ \end{align} In an implicit method (Backward Euler), we evaluate \(f\) in \ref{7fv712} at the future (unknown) time level \(t+\Delta t\): \begin{align} &c_C^{(t+\Delta t)} = c_C^t + \Delta t \overbrace{f(c_W^{(t+\Delta t)}, c_C^{(t+\Delta t)}, c_E^{(t+\Delta t)})}^\text{evaluated at future time $t+\Delta t$}\label{7fv714}\\ \end{align} ## Explicit 1-d stencil for diffusion in porous media Here, the spatial terms on the right are evaluated using the concentrations at the known time level, \(c_W^t\), \(c_C^t\), \(c_E^t\): \begin{align} &{c_C^{(t+\Delta t)} - c_C^t\over \Delta t}\theta (\Delta x)(\Delta y)(\Delta z)= \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_E^t - \left(2 D\theta {\Delta y \Delta z \over \Delta x}\right) c_C^t + \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_W^t \label{7fv715}\\ \end{align}
Implicit 1-d stencil for diffusion in porous media¶
Here, the spatial terms on the right are evaluated using the concentrations at the future \(t+\Delta t\) time level, \(c_W^{(t+\Delta t)}\), \(c_C^{(t+\Delta t)}\), \(c_E^{(t+\Delta t)}\): \begin{align} &{c_C^{(t+\Delta t)} - c_C^t\over \Delta t}\theta (\Delta x)(\Delta y)(\Delta z) = \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_E^{(t+\Delta t)} - \left(2 D\theta {\Delta y \Delta z \over \Delta x}\right) c_C^{(t+\Delta t)} + \left(D\theta {\Delta y \Delta z \over \Delta x}\right) c_W^{(t+\Delta t)}\label{7fv716}\\ \end{align}
Your turn¶
Re-arrange the explicit equation \ref{7fv715} to collect together the terms that multiply each variable, with the known variables on the right side (those at time level \(t\)) and the unknowns on the left (those at time level \(t+\Delta t\).
Put the result here
Your turn¶
Re-arrange the implicit equation \ref{7fv716} to collect together the terms that multiply each variable, with the known variables on the right side (those at time level \(t\)) and the unknowns on the left (those at time level \(t+\Delta t\):
Put the result here
Summary¶
Explicit method¶
initial conditions provide solution for the first \(t_0\) time step
each individual unknown at the next time level \(t+\Delta t\) can be computed from an explicit formula containing values from the known time level \(t\).
no matrix equation is necessary - the value of each unknown can be written as an explicit expression of values from a known time level.
each unknown at time level \(t+\Delta t\) is independent of other values at the \(t+\Delta t\) time level.
Conceptual model of the explicit method: For a 1-D diffusion problem, the C node at time level \(t+\Delta t\) depends on the values of C, W, E nodes at time level \(t\):
Implicit method¶
initial conditions provide solution for the first \(t_0\) time step
spatial derivatives evaluated at time level \(t+\Delta t\).
cannot rearrange equation with one unknown at \(t+\Delta t\) in terms of values from the known time level \(t\).
BUT can write one equation for each gridblock - the solution is IMPLICIT in the system of equations.
must solve a system of equations at each time step.
each unknown at time level \(t+\Delta t\) is dependent on the other values at the \(t+\Delta t\) time level.
Conceptual model for the 1-D implicit diffusion stencil:
## Accuracy, stability and time-step constraints
Accuracy¶
For both explicit and implicit methods, the accuracy will improve with smaller time steps. That is good. What is the cost of using smaller time steps?
Stability¶
A stable time discretization is one where round-off errors and small errors in initial conditions do not grow large with each time step. In an unstable scheme, roundoff errors quickly pollute the solution beyond recognition.
Explicit methods are conditionally stable¶
Explicit methods are conditionally stable, with conditions that depend upon the gridblock size, the parameters of the system and the time step. The conditions on stability are (for one-dimensional diffusion): \begin{align} {\Delta t\over (\Delta x)^2}< {1\over 2 \max D} \label{7fv717}\\ \end{align} where \(D\) is the diffusion coefficient.
Your turn¶
What is the largest time step (in seconds) that can be used for a grid where \(\Delta x=10~cm\) and \(D = 10^{-9}~m^2/s\).
[1]:
### BEGIN SOLUTION
dt = (0.1) * (0.1) / (2.0 * 1e-9)
### END SOLUTION
Implicit methods are unconditionally stable¶
Implicit methods are unconditionally stable. While large time steps may lead to inaccurate solutions, the solutions will not become unstable. For this reason, implicit methods are preferred for their stability.
Reflection¶
What are the trade-offs between an explicit scheme and an implicit scheme?
What considerations should you take into account when choosing a time-integration scheme and a time-step size?
[ ]: