R: Diffusive Transport in cylindrical (r, theta, z) and...
tran.cylindrical
R Documentation
Diffusive Transport in cylindrical (r, theta, z) and spherical (r, theta, phi)
coordinates.
Description
Estimates the transport term (i.e. the rate of change of a concentration
due to diffusion) in a cylindrical (r, theta, z) or spherical (r, theta, phi)
coordinate system.
concentration, expressed per unit volume, defined at the centre
of each grid cell; Nr*Nteta*Nz (cylindrica) or Nr*Ntheta*Nphi (spherical
coordinates) array [M/L3].
C.r.up
concentration at upstream boundary in r(x)-direction;
one value [M/L3].
C.r.down
concentration at downstream boundary in r(x)-direction;
one value [M/L3].
C.theta.up
concentration at upstream boundary in theta-direction;
one value [M/L3].
C.theta.down
concentration at downstream boundary in theta-direction;
one value [M/L3].
C.z.up
concentration at upstream boundary in z-direction (cylindrical
coordinates); one value [M/L3].
C.z.down
concentration at downstream boundary in z-direction(cylindrical
coordinates); one value [M/L3].
C.phi.up
concentration at upstream boundary in phi-direction (spherical
coordinates); one value [M/L3].
C.phi.down
concentration at downstream boundary in phi-direction(spherical
coordinates); one value [M/L3].
flux.r.up
flux across the upstream boundary in r-direction,
positive = INTO model domain; one value [M/L2/T].
flux.r.down
flux across the downstream boundary in r-direction,
positive = OUT of model domain; one value [M/L2/T].
flux.theta.up
flux across the upstream boundary in theta-direction,
positive = INTO model domain; one value [M/L2/T].
flux.theta.down
flux across the downstream boundary in theta-direction,
positive = OUT of model domain; one value [M/L2/T].
flux.z.up
flux across the upstream boundary in z-direction(cylindrical
coordinates); positive = INTO model domain; one value [M/L2/T].
flux.z.down
flux across the downstream boundary in z-direction,
(cylindrical coordinates); positive = OUT of model domain; one value [M/L2/T].
flux.phi.up
flux across the upstream boundary in phi-direction(spherical
coordinates); positive = INTO model domain; one value [M/L2/T].
flux.phi.down
flux across the downstream boundary in phi-direction,
(spherical coordinates); positive = OUT of model domain; one value [M/L2/T].
cyclicBnd
If not NULL, the direction in which a cyclic
boundary is defined, i.e. cyclicBnd = 1 for the r direction,
cyclicBnd = 2 for the theta direction and
cyclicBnd = c(1,2) for both the r and theta direction.
D.r
diffusion coefficient in r-direction, defined on grid cell
interfaces. One value or a vector of length (Nr+1), [L2/T].
D.theta
diffusion coefficient in theta-direction, defined on grid cell
interfaces. One value or or a vector of length (Ntheta+1), [L2/T].
D.z
diffusion coefficient in z-direction, defined on grid cell
interfaces for cylindrical coordinates. One value or a vector of length
(Nz+1) [L2/T].
D.phi
diffusion coefficient in phi-direction, defined on grid cell
interfaces for cylindrical coordinates. One value or a vector of length
(Nphi+1) [L2/T].
r
position of adjacent cell interfaces in the r-direction.
A vector of length Nr+1 [L].
theta
position of adjacent cell interfaces in the theta-direction.
A vector of length Ntheta+1 [L]. Theta should be within [0,2 pi]
z
position of adjacent cell interfaces in the z-direction (cylindrical
coordinates). A vector of length Nz+1 [L].
phi
position of adjacent cell interfaces in the phi-direction (spherical
coordinates). A vector of length Nphi+1 [L]. Phi should be within [0,2 pi]
Details
tran.cylindrical performs (diffusive) transport in cylindrical coordinates
tran.spherical performs (diffusive) transport in spherical coordinates
The boundary conditions are either
(1) zero gradient
(2) fixed concentration
(3) fixed flux
(4) cyclic boundary
This is also the order of priority. The cyclic boundary overrules the other.
If fixed concentration, fixed flux, and cyclicBnd are NULL then
the boundary is zero-gradient
A cyclic boundary condition has concentration and flux at upstream and
downstream boundary the same. It is useful mainly for the theta and
phi direction.
** Do not expect too much of this equation: do not try to use it with
many boxes **
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 Nr*Nteta*Nz (cylindrical) or
Nr*Ntheta*Nphi (spherical coordinates) array. [M/L3/T].
flux.r.up
flux across the upstream boundary in r-direction,
positive = INTO model domain. A matrix of dimension Nteta*Nz (cylindrical)
or Ntheta*Nphi (spherical) [M/L2/T].
flux.r.down
flux across the downstream boundary in r-direction,
positive = OUT of model domain. A matrix of dimension Nteta*Nz (cylindrical)
or Ntheta*Nphi (spherical) [M/L2/T].
flux.theta.up
flux across the upstream boundary in theta-direction,
positive = INTO model domain. A matrix of dimension Nr*Nz (cylindrical) or
or Nr*Nphi (spherical) [M/L2/T].
flux.theta.down
flux across the downstream boundary in theta-direction,
positive = OUT of model domain. A matrix of dimension Nr*Nz (cylindrical) or
Nr*Nphi (spherical) [M/L2/T].
flux.z.up
flux across the upstream boundary in z-direction,
for cylindrical coordinates;
positive = OUT of model domain. A matrix of dimension Nr*Nteta [M/L2/T].
flux.z.down
flux across the downstream boundary in z-direction
for cylindrical coordinates;
positive = OUT of model domain. A matrix of dimension Nr*Nteta [M/L2/T].
flux.phi.up
flux across the upstream boundary in phi-direction,
for spherical coordinates;
positive = OUT of model domain. A matrix of dimension Nr*Nteta [M/L2/T].
flux.phi.down
flux across the downstream boundary in phi-direction,
for spherical coordinates;
positive = OUT of model domain. A matrix of dimension Nr*Nteta [M/L2/T].
See Also
tran.polar
for a discretisation of 2-D transport equations in polar coordinates
tran.1D, tran.2D, tran.3D
Examples
## =============================================================================
## Testing the functions
## =============================================================================
# Grid definition
r.N <- 4 # number of cells in r-direction
theta.N <- 6 # number of cells in theta-direction
z.N <- 3 # number of cells in z-direction
D <- 100 # diffusion coefficient
r <- seq(0, 8, len = r.N+1) # cell size r-direction [cm]
theta <- seq(0,2*pi, len = theta.N+1) # theta-direction - theta: from 0, 2pi
phi <- seq(0,2*pi, len = z.N+1) # phi-direction (0,2pi)
z <- seq(0,5, len = z.N+1) # cell size z-direction [cm]
# Intial conditions
C <- array(dim = c(r.N, theta.N, z.N), data = 0)
# Concentration boundary conditions
tran.cylindrical (C = C, D.r = D, D.theta = D,
C.r.up = 1, C.r.down = 1,
C.theta.up = 1, C.theta.down = 1,
C.z.up = 1, C.z.down = 1,
r = r, theta = theta, z = z )
tran.spherical (C = C, D.r = D, D.theta = D,
C.r.up = 1, C.r.down = 1, C.theta.up = 1, C.theta.down = 1,
C.phi.up = 1, C.phi.down = 1,
r = r, theta = theta, phi = phi)
# Flux boundary conditions
tran.cylindrical(C = C, D.r = D, r = r, theta = theta, z = z,
flux.r.up = 10, flux.r.down = 10,
flux.theta.up = 10, flux.theta.down = 10,
flux.z.up = 10, flux.z.down = 10)
tran.spherical(C = C, D.r = D, r = r, theta = theta, phi = phi,
flux.r.up = 10, flux.r.down = 10,
flux.theta.up = 10, flux.theta.down = 10,
flux.phi.up = 10, flux.phi.down = 10)
# cyclic boundary conditions
tran.cylindrical(C = C, D.r = D, r = r, theta = theta, z = z,
cyclicBnd = 1:3)
tran.spherical(C = C, D.r = D, r = r, theta = theta, phi = phi,
cyclicBnd = 1:3)
# zero-gradient boundary conditions
tran.cylindrical(C = C, D.r = D, r = r, theta = theta, z = z)
tran.spherical(C = C, D.r = D, r = r, theta = theta, phi = phi)
## =============================================================================
## A model with diffusion and first-order consumption
## =============================================================================
N <- 10 # number of grid cells
rr <- 0.005 # consumption rate
D <- 400
r <- seq (2, 4, len = N+1)
theta <- seq (0, 2*pi, len = N+1)
z <- seq (0, 3, len = N+1)
phi <- seq (0, 2*pi, len = N+1)
# The model equations
Diffcylin <- function (t, y, parms) {
CONC <- array(dim = c(N, N, N), data = y)
tran <- tran.cylindrical(CONC,
D.r = D, D.theta = D, D.z = D,
r = r, theta = theta, z = z,
C.r.up = 0, C.r.down = 1,
cyclicBnd = 2)
dCONC <- tran$dC - rr * CONC
return (list(dCONC))
}
Diffspher <- function (t, y, parms) {
CONC <- array(dim = c(N, N, N), data = y)
tran <- tran.spherical (CONC,
D.r = D, D.theta = D, D.phi = D,
r = r, theta = theta, phi = phi,
C.r.up = 0, C.r.down = 1,
cyclicBnd = 2:3)
dCONC <- tran$dC - rr * CONC
return (list(dCONC))
}
# initial condition: 0 everywhere, except in central point
y <- array(dim = c(N, N, N), data = 0)
N2 <- ceiling(N/2)
y[N2, N2, N2] <- 100 # initial concentration in the central point...
# solve to steady-state; cyclicBnd = 2,
outcyl <- steady.3D (y = y, func = Diffcylin, parms = NULL,
dim = c(N, N, N), lrw = 1e6, cyclicBnd = 2)
STDcyl <- array(dim = c(N, N, N), data = outcyl$y)
image(STDcyl[,,1])
# For spherical coordinates, cyclic Bnd = 2, 3
outspher <- steady.3D (y = y, func = Diffspher, parms = NULL, pos=TRUE,
dim = c(N, N, N), lrw = 1e6, cyclicBnd = 2:3)
#STDspher <- array(dim = c(N, N, N), data = outspher$y)
#image(STDspher[,,1])
## Not run:
image(outspher)
## End(Not run)