Report

Task 1

1.1 Derive the finite element matrices \(\mathbf{M}\), \(\mathbf{K}\), \(\mathbf{C}\) and \(\mathbf{R}\) from the strong form advection/diffusion/reaction equation. (2.5 points)

Move everything depending on \(c\) to the left hand side (can also be done in a later step), multiply with \(w\) and integrate over the domain

\[ \int_\Omega w\left(\frac{\partial c}{\partial t} - \kappa\nabla\cdot\nabla c + \mathbf{v}\cdot\nabla c + r c\right) \,\mathrm{d}\Omega = 0 \quad \forall \quad w(x,y) \]

Apply integration by parts only on the term with \(\kappa\)

\[ \int_\Omega w\frac{\partial c}{\partial t} + \kappa\nabla w\cdot\nabla c + w\mathbf{v}\cdot\nabla c + w r c \,\mathrm{d}\Omega - \int_\Gamma w\kappa\nabla c\cdot \mathbf{n} \,\mathrm{d}\Gamma = 0 \quad \forall \quad w(x,y) \]

Substitution of Neumann boundary conditions and setting \(w=0\) on the part of the boundary where Dirichlet boundary are applied and moving to right hand side:

\[ \int_\Omega w\frac{\partial c}{\partial t} + \kappa\nabla w\cdot\nabla c + w\mathbf{v}\cdot\nabla c + w r c \,\mathrm{d}\Omega = \int_{\Gamma_N} wq \,\mathrm{d}\Gamma \quad \forall \quad w(x,y) \]

Substitution of \(w = \mathbf{N}\mathbf{w}\), \(\nabla w = \mathbf{Bw}\), \(c = \mathbf{N}\mathbf{c}\), \(\nabla c = \mathbf{B}\mathbf{c}\) gives:

\[ \int_\Omega \mathbf{Nw}\frac{\partial \mathbf{Nc}}{\partial t} + \kappa\mathbf{Bw}\cdot\mathbf{Bc} + \mathbf{Nw}\mathbf{v}\cdot\mathbf{Bc} + \mathbf{Nw} r \mathbf{Nc} \,\mathrm{d}\Omega = \int_{\Gamma_N} \mathbf{Nw}q \,\mathrm{d}\Gamma \quad \forall \quad \mathbf{w} \]

Moving \(\mathbf{w}\) and \(\mathbf{c}\) out of the integral, removing \(\mathbf{w}\), accounting for the fact that \(\mathbf{N}\) is independent of time:

\[ \int_\Omega \mathbf{N}^T\mathbf{N}\,\mathrm{d}\Omega \frac{\partial \mathbf{c}}{t} + \int_\Omega \kappa\mathbf{B}^T\mathbf{B}\,\mathrm{d}\Omega\mathbf{c} + \int_\Omega\mathbf{N}^T\left(\mathbf{v}\cdot\mathbf{B}\right)\,\mathrm{d}\Omega\mathbf{c} + \int_\Omega\mathbf{N}^T r \mathbf{N} \,\mathrm{d}\Omega\mathbf{c} = \int_{\Gamma_N} \mathbf{N}^Th \,\mathrm{d}\Gamma \quad \]

So we have:

\[\begin{split} \mathbf{M} = \int_\Omega \mathbf{N}^T\mathbf{N}\,\mathrm{d}\Omega \\ \mathbf{K} = \int_\Omega \mathbf{B}^T\kappa\mathbf{B}\,\mathrm{d}\Omega \\ \mathbf{C} = \int_\Omega \mathbf{N}^T\left(\mathbf{v}\cdot\mathbf{B}\right)\,\mathrm{d}\Omega \\ \mathbf{R} = \int_\Omega \mathbf{N}^Tr\mathbf{N}\,\mathrm{d}\Omega \\ \mathbf{f} = \int_{\Gamma_N} \mathbf{N}^Th\,\mathrm{d}\Gamma \end{split}\]

We are not very strict about having a consistent notation for the products in \(\mathbf{C}\). A good alternative is the one below, but the one above, with or without parentheses, also receives a full score.

\[ \mathbf{C} = \int_\Omega \mathbf{N}^T \mathbf{v}^T \mathbf{B} \]

Task 2

2.1 In which part of the code are Dirichlet boundary conditions applied? Include the function, code line(s) and any parameters that need to be defined (0.5 point)

In the simulate function itself, when only the [free_nodes,free_nodes] part of the matrix is considered for solving the system of equations, while values for the constrained_nodes are set to 0.

    Kmodff = Kmod[np.ix_(free_nodes, free_nodes)]
        :
        concentrations[i+1, constrained_nodes] = 0

The part where constrained_nodes are selected is also relevant but not essential. It is not where the boundary conditions are applied.

2.2 In which part of the code are Neumann boundary conditions applied? Include the function, code line(s) and any parameters that need to be defined (0.5 point)

In the get_global_f function

2.3 What are the effective boundary conditions for the parts of the boundaries on which the simulation does not do anything? Motivate your answer with an observation about the simulation results. (1 point)

Homogeneous Neumann boundary conditions \(q=0\). Observed by contours remaining normal to the boundary.

2.4 Give a mathematical expression for the boundary conditions that are applied on the different boundaries of the domain (0.5 point)

\(c=0\) at the left boundary

\(q=0\) at the right boundary, bottom boundary and interior boundary

\(q=1\cdot10^9\left(1+\cos\left(\frac{2\pi t}{25000}\right)\right)\) while \(t<25000\) and \(q=0\) afterwards at the top boundary for \(x\in[50,80]\)

\(q=0\) at the remainder of the top boundary

Task 3

3.1 What happens if you assume the contaminant is inert by setting \(r=0\)? Include any quantitative values that change (1 point)

See solution notebook for graphs. They should be included in the student reports.

The peak is much higher (about ~100 and ~200 times for the two measurement nodes) and arrives later. With nonzero \(r\) the peak is happening earlier because concentrations get reduced progressively over time.

3.2 What happens if the diffusivity is reduced with a factor 10 (use the original value for \(r\))? (1 point)

For node 776 (closer to the surface) the peak is similar of similar height to that from the reference case. The peak arrives later. Reduced diffusion keeps the contimanant cloud more concentrated, but because it takes longer for the contaminant to reach the measurement point, there is also more time for reactive decay. For the node 1176 (deeper underground) the peak is much lower. As the flow is mostly parallel to the water table and the spill is at the surface, diffusion is needed for the contaminant to reach depth.

Task 4

4.1 If you want to code a simulate function to compute the hydraulic head distribution on the domain, which of the finite element functions in this notebook could you reuse? Also mention any modifications you would make, if any. (0.5 point)

The simulate function can use get_global_K and get_global_f in unmodified form.

4.2 How would you deal with the boundary conditions for this problem? (1.0 point)

Neumann boundary conditions should be passed to the get_global_f function. Nonzero Dirichlet boundary conditions need to be applied by modifying the right-hand side vector when eliminating the degrees of freedom on which \(h\) is prescribed.

4.3 In this notebook, we have used velocity vectors inside each element. After solving for \(h\) on the nodes, how would you compute the velocity at the center of a given element in the mesh? The correct answer here involves the finite element shape functions. (0.5 point)

In every element, we can evaluate the \(\mathbf{B}\) matrix with the derivative of the finite element shape functions. We can then compute \(\nabla h=\mathbf{B} \mathbf{h}_\mathrm{elem}\) where \(\mathbf{h}_\mathrm{elem}\) contains the three nodal values of \(h\) for the element. Finally, we can use

\[v=-(\alpha/\mu)\nabla h \]

to compute the velocity. For a given element with nodes elnodes and nodal coordinates n1, n2 and n3 all this could be done in one line as:

v_local = -alpha / porosity * getBMatrix(n1, n2, n3) @ head[elnodes]

See solution notebook for a complete implementation.

By Frans van der Meer, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.