Table of Contents

  • 1  Finite volumes 6 - flux boundary conditions

  • 2  Learning goals

  • 3  Prescribed flux boundary conditions

    • 3.1  Summary of finite-volume equations (1-d, steady-state, no sources within the volume)

      • 3.1.1  Fundamental 1-d steady-state FV equation

      • 3.1.2  In terms of specific fluxes and areas of finite-volume faces (for regular retangular gridblocks):

      • 3.1.3  Location of flux boundary condition

      • 3.1.4  Implications for grid and boundary condition locations

      • 3.1.5  Example

      • 3.1.6  Prescribed flux boundary condition in a discrete approximation

      • 3.1.7  In terms of concentrations in gridblocks:

    • 3.2  Your turn

    • 3.3  Your turn - general case

    • 3.4  Your turn

  • 4  2-D steady-state

    • 4.1  Your turn

      • 4.1.1  Question

  • 5  Gridblock example - flux boundary conditions

    • 5.1  Your turn

Finite volumes 6 - flux boundary conditions

Learning goals

  1. Be able to apply flux boundary conditions to a finite volume approximation.

  2. To be able to compute local and global mass balance from a finite-volume discrete approximation.

  3. Justify the use of harmonic averaging for heterogeneous diffusion coefficients to maintain mass balance.

  4. Be able to express a finite-volume formulation in integral notation and identify each term in the integral expression.

  5. Be able to describe the steps required to develop a general finite volume approximation for arbitrarily shaped elements.

Prescribed flux boundary conditions

For many problems, a prescribed flux is a natural boundary condition: * for solute diffusion - it prescribes the rate at which solute is fluxing into the problem domain (eg sulfate). * for water flow, a flux boundary condition describes the rate at which water enters into the problem domain. * perhaps the most common is a zero flux boundary condition - where the boundary is prescribed to not allow flux of a quantity across the boundary.

Well first review the fundamental FV equations and show that flux boundary conditions can be easily incorporated into our stencils.

Summary of finite-volume equations (1-d, steady-state, no sources within the volume)

Fundamental 1-d steady-state FV equation

We developed the finite-volume approximations starting from the fundamental conservation equation which is already written in terms of fluxes:

\begin{align} \left(J_{EC}+J_{WC}\right) &= 0 \label{eq761}\\ \end{align}

In terms of specific fluxes and areas of finite-volume faces (for regular retangular gridblocks):

\begin{align} &\left(j_{EC}+j_{WC}\right) (\Delta y) (\Delta z) = 0 \label{eq762}\\ \end{align}

If our gridblock C lies on the boundary, a prescribed flux boundary condition means that one fluxes entering the gridblock is prescribed (known).

Location of flux boundary condition

For our 1-d problem, with the fluxes aligned with the \(x\) direction, we have something that looks like

The fluxes \(j_{EC}\) and \(j_{WC}\) enter at the boundary of gridblock C - between EC and WC:

It is therefore best to conceive of the fluxes as being located at the boundary between gridblocks.

Implications for grid and boundary condition locations

When we are modeling a physical system: * we place the edge of the gridblock on the boundary if the boundary is a prescribed flux boundary; * we place the node on the boundary if it is a prescribed concentration boundary.

Example

Consider our simple 5-node example. In the physical system to be modeled, a flux of \(j=0.37\times 10^{-8}~mg/(s\cdot m^2)\) of sulfate is entering from the left of the domain and the concentration is prescribed (fixed) at \(c=9 ~mg/L\) on the right boundary:

Prescribed flux boundary condition in a discrete approximation

Prescribing a flux boundary condition is equivalent to prescribing the value of \(j_{WC}\) or \(j_{EC}\) in a discrete approximation. So in our example above, for gridblock 1, the equation would look like: \begin{align} &\left(j_{EC}+j_{WC}\right) (\Delta y) (\Delta z) = 0 \label{eq763}\\ &\left(j_{21}+0.37\times 10^{-8}\right) (\Delta y) (\Delta z) = 0 \label{eq764}\\ \end{align}

Be careful with the signs! Is the sign correct for the flux if it is entering gridblock 1?

In terms of concentrations in gridblocks:

Below is the equation for a generic gridblock C in the middle of the mesh. Let’s see below how to modify it for a flux boundary condition. \begin{align} &\left(D\theta {c_E - c_C \over \Delta x} + D\theta {c_W - c_C \over \Delta x} \right) (\Delta y) (\Delta z) = 0 \label{eq765}\\ \end{align}

Your turn

For the 5-node example above, write the equation for gridblock 1, where \(\Delta x=1~m\), \(\Delta y = 1.5~m\), \(\Delta z = 2.0~m\), \(D\theta=10^{-9}~m^2/s\) and the flux on the left edge of gridblock 1 is \(j = 0.37\times 10^{-8}~ mg/(s\cdot m^2)\): \begin{align} &\left(D\theta {c_E - c_C \over \Delta x} + D\theta {c_W - c_C \over \Delta x} \right) (\Delta y) (\Delta z) = 0 \label{eq766}\\\end{align}

Your turn - general case

Modify the equation below for the general case where the flux from W is prescribed to be \(j_{WC} = f\), \(f>0\). \begin{align} &\left(D\theta {c_E - c_C \over \Delta x} + D\theta {c_W - c_C \over \Delta x} \right) (\Delta y) (\Delta z) = 0 \label{eq767}\\\end{align}

Your turn

Modify the generic equation below for the case where the flux from E is prescribed to be \(j_{EC} = 0\), (no-flux boundary condition): \begin{align} &\left(D\theta {c_E - c_C \over \Delta x} + D\theta {c_W - c_C \over \Delta x} \right) (\Delta y) (\Delta z) = 0 \label{eq768}\\ \end{align}

2-D steady-state

In 2-D, the general stencils for an interior gridblock C in terms of total flux, specific flux and concentrations are: \begin{align} &\left(J_{EC}+J_{WC} + J_{NC}+J_{SC}\right) = 0 \label{eq769}\\ \end{align}

\begin{align} &\left(j_{EC}+j_{WC}\right) (\Delta y) (\Delta z) +\left(j_{NC}+j_{SC}\right) (\Delta x) (\Delta z) = 0 \label{eq7610}\\ \end{align}

\begin{align} &\left(D\theta {c_E - c_C \over \Delta x} + D\theta {c_W - c_C \over \Delta x} \right) (\Delta y) (\Delta z) +\left(D\theta {c_N - c_C \over \Delta y} + D\theta {c_S - c_C \over \Delta y} \right) (\Delta x) (\Delta z)= 0 \label{eq7611}\\ \end{align}

Your turn

Modify each of the three equations below for the situation where there is no flux to or from gridblocks N and S:

\begin{align} \left(J_{EC}+J_{WC} + J_{NC}+J_{SC}\right) = 0 \label{eq7612}\\ \end{align}

\begin{align} \left(j_{EC}+j_{WC}\right) (\Delta y) (\Delta z) +\left(j_{NC}+j_{SC}\right) (\Delta x) (\Delta z) = 0 \label{eq7613}\\ \end{align}

\begin{align} &\left(D\theta {c_E - c_C \over \Delta x} + D\theta {c_W - c_C \over \Delta x} \right) (\Delta y) (\Delta z) +\left(D\theta {c_N - c_C \over \Delta y} + D\theta {c_S - c_C \over \Delta y} \right) (\Delta x) (\Delta z)= 0\label{eq7614}\\ \end{align}

Question

The result should be a familiar stencil. What is it?

Gridblock example - flux boundary conditions

Let’s write the equations for the 5-gridblock example from before, now with a flux boundary condition as shown (note \(x\) and \(y\) coordinate directions are as indicated,with \(z\) into the page).

The total length of the domain is \(2~m\) in the \(x\) direction, each gridblock is the same size and porosity and diffusion coefficient are spatially homogeneous: * \(\Delta x=2/4.5 = 0.44~ m\) * \(\Delta y=0.3~m\) * \(\Delta z=0.5~m\) * \(D=2\times 10^{-10}~m^2/s\) * \(\theta = 0.25\)

Note that the edge of gridblock 1 is on the physical boundary whereas the node of gridblock 5 is on the physical boundary. Accordingly, there are 4.5 gridblocks across the 2-m domain, so each gridblock is 0.44 m in the x direction.

Your turn

  1. Write the equations for the unknowns in the problem.

  2. Solve the system of equations for the unknown concentrations. (If you are doing this by hand, note that you can rescale any equation by dividing or multiplying all terms in the equation by the same number. Often it is convenient to divide all terms by the coefficient of the unknown on the “diagonal” - for equation number :math:`n`, the diagonal is the coefficient that multiplies :math:`c_n`.)

  3. In the cell below, report the concentration you computed in gridblock 1.

  4. In the cell after that, use concentration values to compute the magnitude (i.e. a positive number) of the specific flux (in \(mg/(s\cdot m^2)\)) between gridblocks 4 and 5.

[ ]:
# YOUR CODE HERE
raise NotImplementedError()

2a. Fill out the following A and b matrices for this problem, overwriting the non-zero values with your coefficients:

A=

 [[0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0.]]

b = [0. 0. 0. 0. 0.]

YOUR ANSWER HERE

[ ]:
# 2b. solve the equations you just wrote down for the concentration c in steadcy state
#
# Our solution in this cell follows the same pattern as the week6_assign.  That is,
# we defined a function called build_1d_matrix_source that produce a 5x5 matrix A
# and a 5x1 vector b, such that the 5 cell concentrations c can be found with the
# following python code
#
# A,b = build_1d_matrix_source(...)
# c = np.linalg.solve(A, b)
#
# You can also show your hand calculations here.  If you do it by hand
# leave the concentrations in a list named c, so we can check the values
# like this
#
# c=[c1,c2,c3,c4,c5]
#
# we will check your c[0]=c1 value against ours in a test below.
#

# YOUR CODE HERE
raise NotImplementedError()
[ ]:
# In this cell show that for you value of the c vector, you do in fact get:
# A@c=b  (either by hand or with python code in this cell)

# YOUR CODE HERE
raise NotImplementedError()
[ ]:
# 3. We check your c[0] for cell 1 here

[ ]:
# 4. We check your value for the flux between grid blocks 4 and 5 (c[3] and c[4])