tran.2D: General Two-Dimensional Advective-Diffusive Transport

In ReacTran: Reactive Transport Modelling in 1d, 2d and 3d

Description

Estimates the transport term (i.e. the rate of change of a concentration due to diffusion and advection) in a two-dimensional model domain.

Usage

tran.2D (C, C.x.up = C[1,], C.x.down = C[nrow(C),],
         C.y.up = C[,1], C.y.down = C[ ,ncol(C)],
         flux.x.up = NULL, flux.x.down = NULL, 
         flux.y.up = NULL, flux.y.down = NULL,
         a.bl.x.up = NULL, a.bl.x.down = NULL, 
         a.bl.y.up = NULL, a.bl.y.down = NULL, 
         D.grid = NULL, D.x = NULL, D.y = D.x,
         v.grid = NULL, v.x = 0, v.y = 0,
         AFDW.grid = NULL, AFDW.x = 1, AFDW.y = AFDW.x,
         VF.grid = NULL, VF.x = 1, VF.y = VF.x,
         A.grid = NULL, A.x = 1, A.y = 1,
         grid = NULL, dx = NULL, dy = NULL,
         full.check = FALSE, full.output = FALSE)

Arguments

C	concentration, expressed per unit volume, defined at the centre of each grid cell; Nx*Ny matrix [M/L3].
C.x.up	concentration at upstream boundary in x-direction; vector of length Ny [M/L3].
C.x.down	concentration at downstream boundary in x-direction; vector of length Ny [M/L3].
C.y.up	concentration at upstream boundary in y-direction; vector of length Nx [M/L3].
C.y.down	concentration at downstream boundary in y-direction; vector of length Nx [M/L3].
flux.x.up	flux across the upstream boundary in x-direction, positive = INTO model domain; vector of length Ny [M/L2/T].
flux.x.down	flux across the downstream boundary in x-direction, positive = OUT of model domain; vector of length Ny [M/L2/T].
flux.y.up	flux across the upstream boundary in y-direction, positive = INTO model domain; vector of length Nx [M/L2/T].
flux.y.down	flux across the downstream boundary in y-direction, positive = OUT of model domain; vector of length Nx [M/L2/T].
a.bl.x.up	transfer coefficient across the upstream boundary layer. in x-direction; Flux=a.bl.x.up*(C.x.up-C[1,]). One value [L/T].
a.bl.x.down	transfer coefficient across the downstream boundary layer in x-direction; Flux=a.bl.x.down*(C[Nx,]-C.x.down). One value [L/T].
a.bl.y.up	transfer coefficient across the upstream boundary layer. in y-direction; Flux=a.bl.y.up*(C.y.up-C[,1]). One value [L/T].
a.bl.y.down	transfer coefficient across the downstream boundary layer in y-direction; Flux=a.bl.y.down*(C[,Ny]-C.y.down). One value [L/T].
D.grid	diffusion coefficient defined on all grid cell interfaces. A prop.2D list created by setup.prop.2D [L2/T]. See last example for creating spatially-varying diffusion coefficients.
D.x	diffusion coefficient in x-direction, defined on grid cell interfaces. One value, a vector of length (Nx+1), a prop.1D list created by setup.prop.1D, or a (Nx+1)* Ny matrix [L2/T].
D.y	diffusion coefficient in y-direction, defined on grid cell interfaces. One value, a vector of length (Ny+1), a prop.1D list created by setup.prop.1D, or a Nx*(Ny+1) matrix [L2/T].
v.grid	advective velocity defined on all grid cell interfaces. Can be positive (downstream flow) or negative (upstream flow). A prop.2D list created by setup.prop.2D [L/T].
v.x	advective velocity in the x-direction, defined on grid cell interfaces. Can be positive (downstream flow) or negative (upstream flow). One value, a vector of length (Nx+1), a prop.1D list created by setup.prop.1D, or a (Nx+1)*Ny matrix [L/T].
v.y	advective velocity in the y-direction, defined on grid cell interfaces. Can be positive (downstream flow) or negative (upstream flow). One value, a vector of length (Ny+1), a prop.1D list created by setup.prop.1D, or a Nx*(Ny+1) matrix [L/T].
AFDW.grid	weight used in the finite difference scheme for advection in the x- and y- direction, defined on grid cell interfaces; backward = 1, centred = 0.5, forward = 0; default is backward. A prop.2D list created by setup.prop.2D [-].
AFDW.x	weight used in the finite difference scheme for advection in the x-direction, defined on grid cell interfaces; backward = 1, centred = 0.5, forward = 0; default is backward. One value, a vector of length (Nx+1), a prop.1D list created by setup.prop.1D, or a (Nx+1)*Ny matrix [-].
AFDW.y	weight used in the finite difference scheme for advection in the y-direction, defined on grid cell interfaces; backward = 1, centred = 0.5, forward = 0; default is backward. One value, a vector of length (Ny+1), a prop.1D list created by setup.prop.1D, or a Nx*(Ny+1) matrix [-].
VF.grid	Volume fraction. A prop.2D list created by setup.prop.2D [-].
VF.x	Volume fraction at the grid cell interfaces in the x-direction. One value, a vector of length (Nx+1), a prop.1D list created by setup.prop.1D, or a (Nx+1)*Ny matrix [-].
VF.y	Volume fraction at the grid cell interfaces in the y-direction. One value, a vector of length (Ny+1), a prop.1D list created by setup.prop.1D, or a Nx*(Ny+1) matrix [-].
A.grid	Interface area. A prop.2D list created by setup.prop.2D [L2].
A.x	Interface area defined at the grid cell interfaces in the x-direction. One value, a vector of length (Nx+1), a prop.1D list created by setup.prop.1D, or a (Nx+1)*Ny matrix [L2].
A.y	Interface area defined at the grid cell interfaces in the y-direction. One value, a vector of length (Ny+1), a prop.1D list created by setup.prop.1D, or a Nx*(Ny+1) matrix [L2].
dx	distance between adjacent cell interfaces in the x-direction (thickness of grid cells). One value or vector of length Nx [L].
dy	distance between adjacent cell interfaces in the y-direction (thickness of grid cells). One value or vector of length Ny [L].
grid	discretization grid, a list containing at least elements dx, dx.aux, dy, dy.aux (see setup.grid.2D) [L].
full.check	logical flag enabling a full check of the consistency of the arguments (default = FALSE; TRUE slows down execution by 50 percent).
full.output	logical flag enabling a full return of the output (default = FALSE; TRUE slows down execution by 20 percent).

Details

The boundary conditions are either

• (1) zero-gradient

• (2) fixed concentration

• (3) convective boundary layer

• (4) fixed flux

This is also the order of priority. The zero gradient is the default, the fixed flux overrules all other.

Value

a list containing:

dC	the rate of change of the concentration C due to transport, defined in the centre of each grid cell, a Nx*Ny matrix. [M/L3/T].
C.x.up	concentration at the upstream interface in x-direction. A vector of length Ny [M/L3]. Only when full.output = TRUE.
C.x.down	concentration at the downstream interface in x-direction. A vector of length Ny [M/L3]. Only when full.output = TRUE.
C.y.up	concentration at the the upstream interface in y-direction. A vector of length Nx [M/L3]. Only when full.output = TRUE.
C.y.down	concentration at the downstream interface in y-direction. A vector of length Nx [M/L3]. Only when full.output = TRUE.
x.flux	flux across the interfaces in x-direction of the grid cells. A (Nx+1)*Ny matrix [M/L2/T]. Only when full.output = TRUE.
y.flux	flux across the interfaces in y-direction of the grid cells. A Nx*(Ny+1) matrix [M/L2/T]. Only when full.output = TRUE.
flux.x.up	flux across the upstream boundary in x-direction, positive = INTO model domain. A vector of length Ny [M/L2/T].
flux.x.down	flux across the downstream boundary in x-direction, positive = OUT of model domain. A vector of length Ny [M/L2/T].
flux.y.up	flux across the upstream boundary in y-direction, positive = INTO model domain. A vector of length Nx [M/L2/T].
flux.y.down	flux across the downstream boundary in y-direction, positive = OUT of model domain. A vector of length Nx [M/L2/T].

Note

It is much more efficient to use the grid input rather than vectors or single numbers.

Thus: to optimise the code, use setup.grid.2D to create the grid, and use setup.prop.2D to create D.grid, v.grid, AFDW.grid, VF.grid, and A.grid, even if the values are 1 or remain constant.

There is no provision (yet) to deal with cross-diffusion. Set D.x and D.y different only if cross-diffusion effects are unimportant.

Author(s)

Filip Meysman <[e-mail]>, Karline Soetaert <[e-mail]>

References

Soetaert and Herman, 2009. a practical guide to ecological modelling - using R as a simulation platform. Springer

See Also

tran.polar for a discretisation of 2-D transport equations in polar coordinates

tran.1D, tran.3D

Examples

## =============================================================================
## Testing the functions
## =============================================================================
# Parameters
F        <- 100             # input flux [micromol cm-2 yr-1]
por      <- 0.8             # constant porosity
D        <- 400             # mixing coefficient [cm2 yr-1]
v        <- 1               # advective velocity [cm yr-1]

# Grid definition
x.N <- 4   # number of cells in x-direction
y.N <- 6   # number of cells in y-direction
x.L <- 8   # domain size x-direction [cm]
y.L <- 24  # domain size y-direction [cm]
dx  <- x.L/x.N             # cell size x-direction [cm]
dy  <- y.L/y.N             # cell size y-direction [cm]
 
# Intial conditions 
C <- matrix(nrow = x.N, ncol = y.N, data = 0, byrow = FALSE)

# Boundary conditions: fixed concentration  
C.x.up   <- rep(1, times = y.N)
C.x.down <- rep(0, times = y.N)
C.y.up   <- rep(1, times = x.N)
C.y.down <- rep(0, times = x.N)

# Only diffusion 
tran.2D(C = C, D.x = D, D.y = D, v.x = 0, v.y = 0,
  VF.x = por, VF.y = por, dx = dx, dy = dy,
  C.x.up = C.x.up, C.x.down = C.x.down,
  C.y.up = C.y.up, C.y.down = C.y.down, full.output = TRUE)

# Strong advection, backward (default), central and forward 
#finite difference schemes 
tran.2D(C = C, D.x = D, v.x = 100*v, 
  VF.x = por, dx = dx, dy = dy,
  C.x.up = C.x.up, C.x.down = C.x.down, 
  C.y.up = C.y.up, C.y.down = C.y.down)
  
tran.2D(AFDW.x = 0.5, C = C, D.x = D, v.x = 100*v, 
  VF.x = por, dx = dx, dy = dy,
  C.x.up = C.x.up, C.x.down = C.x.down, 
  C.y.up = C.y.up, C.y.down = C.y.down)

tran.2D(AFDW.x = 0, C = C, D.x = D, v.x = 100*v, 
  VF.x = por, dx = dx, dy = dy,
  C.x.up = C.x.up, C.x.down = C.x.down, 
  C.y.up = C.y.up, C.y.down = C.y.down)

# Boundary conditions: fixed fluxes 

flux.x.up <- rep(200, times = y.N)
flux.x.down <- rep(-200, times = y.N)
flux.y.up <- rep(200, times = x.N)
flux.y.down <- rep(-200, times = x.N)
tran.2D(C = C, D.x = D, v.x = 0, 
  VF.x = por, dx = dx, dy = dy,
  flux.x.up = flux.x.up, flux.x.down = flux.x.down,
  flux.y.up = flux.y.up, flux.y.down = flux.y.down)

# Boundary conditions: convective boundary layer on all sides

a.bl <- 800   # transfer coefficient
C.x.up <- rep(1, times = (y.N)) # fixed conc at boundary layer
C.y.up <- rep(1, times = (x.N)) # fixed conc at boundary layer
tran.2D(full.output = TRUE, C = C, D.x = D, v.x = 0, 
  VF.x = por, dx = dx, dy = dy, 
  C.x.up   = C.x.up, a.bl.x.up = a.bl, 
  C.x.down = C.x.up, a.bl.x.down = a.bl, 
  C.y.up   = C.y.up, a.bl.y.up = a.bl,
  C.y.down = C.y.up, a.bl.y.down = a.bl)

# Runtime test with and without argument checking

n.iterate <-500

test1 <- function() {
  for (i in 1:n.iterate )
    ST <- tran.2D(full.check = TRUE, C = C, D.x = D, 
      v.x = 0, VF.x = por, dx = dx, dy = dy,
      C.x.up = C.x.up, a.bl.x.up = a.bl, C.x.down = C.x.down)
} 
system.time(test1())

test2 <- function() {
  for (i in 1:n.iterate )
    ST <- tran.2D(full.output = TRUE, C = C, D.x = D, 
      v.x = 0, VF.x = por, dx = dx, dy = dy,
      C.x.up = C.x.up, a.bl.x.up = a.bl, C.x.down = C.x.down)
} 
system.time(test2())

test3 <- function() {
  for (i in 1:n.iterate )
    ST <- tran.2D(full.output = TRUE, full.check = TRUE, C = C,
      D.x = D, v.x = 0, VF.x = por, dx = dx, dy = dy,
      C.x.up = C.x.up, a.bl.x.up = a.bl, C.x.down = C.x.down)
} 
system.time(test3())

## =============================================================================
## A 2-D model with diffusion in x- and y direction and first-order
## consumption - unefficient implementation
## =============================================================================

N     <- 51          # number of grid cells
XX    <- 10           # total size
dy    <- dx <- XX/N  # grid size
Dy    <- Dx <- 0.1   # diffusion coeff, X- and Y-direction
r     <- 0.005       # consumption rate
ini   <- 1           # initial value at x=0

N2  <- ceiling(N/2)
X   <- seq (dx, by = dx, len = (N2-1))
X   <- c(-rev(X), 0, X)

# The model equations

Diff2D <- function (t, y, parms)  {

 CONC  <- matrix(nrow = N, ncol = N, y)
 dCONC <- tran.2D(CONC, D.x = Dx, D.y = Dy, dx = dx, dy = dy)$dC + r * CONC

 return (list(dCONC))

}

# initial condition: 0 everywhere, except in central point
y <- matrix(nrow = N, ncol = N, data = 0)
y[N2, N2] <- ini  # initial concentration in the central point...

# solve for 10 time units
times <- 0:10
out <- ode.2D (y = y, func = Diff2D, t = times, parms = NULL,
                dim = c(N,N), lrw = 160000)

pm <- par (mfrow = c(2, 2))

# Compare solution with analytical solution...
for (i in seq(2, 11, by = 3))  {
  tt   <- times[i]
  mat  <-  matrix(nrow = N, ncol = N, 
                  data = subset(out, time == tt))
  plot(X, mat[N2,], type = "l", main = paste("time=", times[i]),
       ylab = "Conc", col = "red")
  ana <- ini*dx^2/(4*pi*Dx*tt)*exp(r*tt-X^2/(4*Dx*tt))
  points(X, ana, pch = "+")
}

legend ("bottom", col = c("red","black"), lty = c(1, NA), 
  pch = c(NA, "+"), c("tran.2D", "exact"))
par("mfrow" = pm )



## =============================================================================
## A 2-D model with diffusion in x- and y direction and first-order
## consumption - more efficient implementation, specifying ALL 2-D grids
## =============================================================================

N     <- 51          # number of grid cells
Dy    <- Dx <- 0.1   # diffusion coeff, X- and Y-direction
r     <- 0.005       # consumption rate
ini   <- 1           # initial value at x=0

x.grid    <- setup.grid.1D(x.up = -5, x.down = 5, N = N)
y.grid    <- setup.grid.1D(x.up = -5, x.down = 5, N = N)
grid2D    <- setup.grid.2D(x.grid, y.grid)

D.grid    <- setup.prop.2D(value = Dx, y.value = Dy, grid = grid2D)
v.grid    <- setup.prop.2D(value = 0, grid = grid2D)
A.grid    <- setup.prop.2D(value = 1, grid = grid2D)
AFDW.grid <- setup.prop.2D(value = 1, grid = grid2D)
VF.grid   <- setup.prop.2D(value = 1, grid = grid2D)

# The model equations - using the grids

Diff2Db <- function (t, y, parms)  {

   CONC  <- matrix(nrow = N, ncol = N, data = y)

   dCONC <- tran.2D(CONC, grid = grid2D, D.grid = D.grid, 
      A.grid = A.grid, VF.grid = VF.grid, AFDW.grid = AFDW.grid, 
      v.grid = v.grid)$dC + r * CONC
  
  return (list(dCONC))
}

# initial condition: 0 everywhere, except in central point
y <- matrix(nrow = N, ncol = N, data = 0)
y[N2,N2] <- ini  # initial concentration in the central point...

# solve for 8 time units
times <- 0:8
outb <- ode.2D (y = y, func = Diff2Db, t = times, parms = NULL,
                dim = c(N, N), lrw = 160000)

image(outb, ask = FALSE, mfrow = c(3, 3), main = paste("time", times))

## =============================================================================
## Same 2-D model, but now with spatially-variable diffusion coefficients
## =============================================================================

N     <- 51          # number of grid cells
r     <- 0.005       # consumption rate
ini   <- 1           # initial value at x=0
N2    <- ceiling(N/2)

D.grid <- list()

# Diffusion on x-interfaces
D.grid$x.int <- matrix(nrow = N+1, ncol = N, data = runif(N*(N+1)))

# Diffusion on y-interfaces
D.grid$y.int <- matrix(nrow = N, ncol = N+1, data = runif(N*(N+1)))

dx <- 10/N
dy <- 10/N

# The model equations

Diff2Dc <- function (t, y, parms)  {

   CONC  <- matrix(nrow = N, ncol = N, data = y)

   dCONC <- tran.2D(CONC, dx = dx, dy = dy, D.grid = D.grid)$dC + r * CONC

  return (list(dCONC))
}

# initial condition: 0 everywhere, except in central point
y <- matrix(nrow = N, ncol = N, data = 0)
y[N2, N2] <- ini  # initial concentration in the central point...

# solve for 8 time units
times <- 0:8
outc <- ode.2D (y = y, func = Diff2Dc, t = times, parms = NULL,
                dim = c(N, N), lrw = 160000)

outtimes <- c(1, 3, 5, 7)
image(outc, ask = FALSE, mfrow = c(2, 2), main = paste("time", outtimes),
      legend = TRUE, add.contour = TRUE, subset = time %in% outtimes)

ReacTran documentation built on May 2, 2019, 9:38 a.m.