#' One-dimensional diffusion, constant boundary conditions
#'
#' @description Function to solve diffusion equation
#' using numerical implicit finite differences method,
#' applied in the condition where: a) the diffusion occurs
#' inside a slab of material unidimensionally,
#' or only along one axis, b) the coefficient of diffusion is constant,
#' and c) the boundary conditions are constant.
#'
#' @usage diff_1D(Lx, Tt, nt, nx, D, C_ini, C_lim)
#'
#' @param Lx Total length of slab in x direction, usually in mm.
#' @param Tt Total diffusion time, usually in second.
#' @param nt Number of time discretization.
#' @param nx Number of dimension x discretization
#' @param D Coefficient of diffusion. Must be constant value.
#' @param C_ini Initial concentration value inside the slab.
#' @param C_lim Boundary condition, written as a vector
#' with one or two elements. If there is only one element,
#' the other side of the slab will be presented as having
#' Neumann boundary condition with flux = 0. If there are
#' two elements, the first element is the dirichlet
#' concentration on the left side, while the second element
#' is the dirichlet concentration on the right side.
#'
#' @return A matrix with {nx+1} number of row and {nt}
#' number of column, profiling the diffusion on slab with Lx
#' length along the time Tt.
#'
#' @examples
#' Lx <- 5 #Length of slab, in mm
#' Tt <- 20 #Total measured diffusion time, in seconds
#' nt <- 100 #Number of time discretization, wherein the Tt will be divided into 100 equal parts (0, 0.2, 0.4, ..., 20)
#' nx <- 50 #Number of dimension discretization
#' D <- 0.05 #Coefficient of diffusion in mm^2/s.
#' C_ini <- 0.05 #Initial concentration inside the slab.
#' C_lim <- c(0.10,0.25) #Dirichlet boundary concentration diffusing into the slab
#' matC <- diff_1D(Lx, Tt, nt, nx, D, C_ini, C_lim)
diff_1D <- function(Lx,Tt,nt,nx,D,C_ini,C_lim){
#1-Dimensional diffusion solver using implicit method, with border condition on its two extremities/sides
##Basic parameters
if(length(C_lim) == 1){
C_fin <- C_lim
} else {
if(length(C_lim == 2)){
C_left <- C_lim[1]
C_right <- C_lim[2]
} else {
stop('ERROR = Unfortunately, this function cannot take more than two condition limits coming from each side')
}
}
dt <- Tt/nt #Size of one step in t
dx <- Lx/(nx+1) #Size of one step in x
##Calculations
s = D*dt/(dx^2)
##
# Matrix A
A = matrix(0, nrow = nx+1, ncol = nx+1)
for (i in 2:nx) {
A[i, i] <- (1+2*s) #Diagonal elements
A[i, (i+1)] <- -s
A[(i+1), i] <- -s
}
A[1,1] <- (1+2*s); A[nx+1,nx+1] <- A[1,1]
A[2,1] <- -s; A[1,2] <- -s
b <- matrix(0, nrow = nx+1, ncol = 1)
#Conditioning
if(length(C_lim) == 2){
print("Boundary conditions: dirichlet in both sides")
b[1] <- s*C_right #Final, right side
b[nx+1] <- s*C_left #Final, left side
} else {
print("Boundary condition: dirichlet in one side and neumann in another")
A[nx+1,nx] <- -2*s #For neumann boundary condition
b[1] <- s*C_fin #Final, right side
b[nx+1] <- 0 #Zero flux, neumann condition, left side
}
Cmat_i <- matrix(C_ini, nrow = nx+1, ncol = 1)
Cmat_f <- matrix(0, nrow = nx+1, ncol = nt)
t_n <- matrix(0, nrow = nt, ncol = 1)
Cmat_f[,1] <- Cmat_i
t_j <- 0 #Initial time
t_n[1] <- t_j
######
##March solution in time
for (j in 1:nt) {
t_j <- t_j+dt
if(j == 1){
Cmat_j <- Cmat_i
} else {
Cmat_j <- Cmat_j_1
}
c <- Cmat_j + b
Cmat_j_1 = solve(A,c)
Cmat_f[,j] <- Cmat_j_1
#Next time step
t_n[j] <- t_j
}
return(Cmat_f)
}
#' Mean concentration values of one-dimensional diffusion during each time steps.
#'
#' @description A wrapper of \emph{diff_1D} function,
#' returning the mean values of concentration during each
#' time step of the diffusion.
#'
#' @usage diff_1D_Ct(Lx, Tt, nt, nx, D, C_ini, C_lim)
#'
#' @param Lx Total length of slab in x direction, usually in mm.
#' @param Tt Total diffusion time, usually in second.
#' @param nt Number of time discretization.
#' @param nx Number of dimension x discretization
#' @param D Coefficient of diffusion. Must be constant value.
#' @param C_ini Initial concentration value inside the slab.
#' @param C_lim Boundary condition, written as a vector with one or two elements. If there is only one element, the other side of the slab will be presented as having Neumann boundary condition with flux = 0. If there are two elements, the first element is the dirichlet concentration on the left side, while the second element is the dirichlet concentration on the right side.
#'
#' @return A vector with \strong{nt} number of elements with the values of mean concentration in the slab for each T\emph{i}.
#'
#' @examples
#' Lx <- 5 #Length of slab, in mm
#' Tt <- 20 #Total measured diffusion time, in seconds
#' nt <- 100 #Number of time discretization, wherein the Tt will be divided into 100 equal parts (0, 0.2, 0.4, ..., 20)
#' nx <- 50 #Number of dimension discretization
#' D <- 0.05 #Coefficient of diffusion in mm^2/s.
#' C_ini <- 0.05 #Initial concentration inside the slab.
#' C_lim <- c(0.10,0.25) #Dirichlet boundary concentration diffusing into the slab
#' meanCt <- diff_1D_Ct(Lx, Tt, nt, nx, D, C_ini, C_lim)
diff_1D_Ct <- function(Lx,Tt,nt,nx,D,C_ini,C_lim){
#Mean concentration value during each time step of 1-dimensional diffusion.
#Will be useful for plotting
Cmat <- diff_1D(Lx,Tt,nt,nx,D,C_ini,C_lim)
Ct <- vector(mode="numeric", length = nt)
for (i in 1:nt) {
Ct[i] <- mean(Cmat[,i])
}
return(Ct)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.