R/FuncHolCyl.R

In ahmad-alkadri/AnalyticDiffusion: Analytical Model of Fick's Second Law Diffusion

#' @title Analytical Model of Fick's Second Law Diffusion on a Hollow Cylindrical Plane
#'
#' @description This package contains several functions for determining the flux of diffusion, estimating the coefficient of diffusion, and basically performing analytical modelization of nonsteady state of diffusion based on Fick's second law of diffusion FOR hollow cylindrical plane.
#'
#' @return NULL
#'
#' @examples
#'
#' @export
func_besselU = function(b,a,v){
  U_v = besselJ(nu=0,x=a*v)*besselY(nu=0,x=b*v)-besselJ(nu=0,x=b*v)*besselY(nu=0,x=a*v)
  return(U_v)
}

root_besselU = function(i,j,N,k){
  n = rootSolve::uniroot.all(func_besselU,
                  a=i,b=j,
                  lower=0.001,
                  upper=N)
  return(n[1:k])
}

a_alpha = function(a,b,N,k){
  a_n = root_besselU(i=a,j=b,N=N,k=k)*a
  return(a_n)
}

Flux_cyl = function(a,b,D,t,N,k){
  alpha = a_alpha(a=a,b=b,N=N,k=k)/a
  Flux = 1-sum((besselJ(x=alpha*a,nu=0)-besselJ(x=alpha*b,nu=0))/(besselJ(x=alpha*a,nu=0)+besselJ(x=alpha*b,nu=0))/alpha^2/exp(D*alpha^2*t))*(4/(b^2-a^2))
  return(Flux)
}

Flux_cyl_val = function(a,b,D,T,N,k){
  temp_val<-c()
  for(x in T){
    temp_val<-c(temp_val,Flux_cyl(a=a,
                                  b=b,
                                  D=D,
                                  t=x,
                                  N=N,
                                  k=k))
  }
  return(temp_val)
}

Flux_exp_val = function(DATA){
  F_val=(DATA$Mass-DATA$Mass[1])/(DATA$Mass[nrow(DATA)]-DATA$Mass[1])
  return(F_val)
}

DiffCoeffCyl = function(DATA,ratio_b_a,N,k){
  sum_temp = summary(minpack.lm::nlsLM(formula = F~Flux_cyl_val(a=1,
                                                    b=ratio_b_a*1,
                                                    D=K,
                                                    T=dT,
                                                    N=N,
                                                    k=k),
                           data = DATA,
                           control = nls.lm.control(maxiter=1024),
                           start = list(K=0.1)))
  DiffCoef_Cyl = as.numeric(sum_temp$coefficients[1])
  return(DiffCoef_Cyl)
}