# Diffusive Transport in cylindrical (r, theta, z) and spherical (r, theta, phi) coordinates. 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. ```r tran.cylindrical (C, C.r.up = NULL, C.r.down = NULL, C.theta.up = NULL, C.theta.down = NULL, C.z.up = NULL, C.z.down = NULL, flux.r.up = NULL, flux.r.down = NULL, flux.theta.up = NULL, flux.theta.down = NULL, flux.z.up = NULL, flux.z.down = NULL, cyclicBnd = NULL, D.r = NULL, D.theta = D.r, D.z = D.r, r = NULL, theta = NULL, z = NULL) tran.spherical (C, C.r.up = NULL, C.r.down = NULL, C.theta.up = NULL, C.theta.down = NULL, C.phi.up = NULL, C.phi.down = NULL, flux.r.up = NULL, flux.r.down = NULL, flux.theta.up = NULL, flux.theta.down = NULL, flux.phi.up = NULL, flux.phi.down = NULL, cyclicBnd = NULL, D.r = NULL, D.theta = D.r, D.phi = D.r, r = NULL, theta = NULL, phi = NULL) ``` ## Arguments - `C`: 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] ## Returns 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]. ## Examples ```r ## ============================================================================= ## 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) ``` ## 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 ** ## See Also `tran.polar` for a discretisation of 2-D transport equations in polar coordinates `tran.1D`, `tran.2D`, `tran.3D`

tran.cyl_spher function