Bonjour,
J'aimerais utiliser tran.1D afin de calculer la concentrations d'espèces chimiques autour d'une sphère. J'ai écrit le code ci-dessous, mais lorsque je veux ajouter des conditions limites et ajouter un coefficient de transfert de la couche limite (a.bl.down) , il y a un message d'erreur qui apparaît : Error in if (v$int[N + 1] >= 0) { : missing value where TRUE/FALSE needed
Comment pourrais-je corriger le problème?
Merci beaucoup,
Michel
Voilà le code :
--------------------------------
library(ReacTran)
#------------------------------#
# MODEL #
#------------------------------#
#--- Partial derivative equations
boundary_layer <- function(t, state, parms) {
{ with( as.list( c(t, state, parms) ), {
# Reshape state variable as a 2D matrix
S <- matrix(nrow = X.grid$N, ncol = 5, data = state)
# Initialize dC/dt matrix
dCdt <- 0*S
# Rate of chage for CO2
tran <- tran.1D(C = S[,1], D = D, a.bl.down = P1, C.down = C[1], A = A.grid, dx = X.grid)
dCdt[,1] <- tran$dC
# Rate of change of SO4
tran <- tran.1D(C = S[,2], D = D, flux.up = F1, C.down = C[1], A = A.grid, dx = X.grid)
dCdt[,2] <- tran$dC
return(list(dCdt))
} )
}
}
#------------------------------#
# Model grid definition #
#------------------------------#
# Number of grid layer
N <- 100
# Radius of the cell or boundary layer thickness (m)
R <- 4E-06
# Model grid setup
X.grid <- setup.grid.1D(x.up = R, L = R, N = N) # x.up = radius of the cell; radius = boundary layer thickness
# Interface area
A.grid <- setup.prop.1D(grid = X.grid, func = function(r) 4*pi*r^2)
#------------------------------#
# Model parameters #
#------------------------------#
# Concentration of CO2,HCO3,CO3,H,OH in bulk solution (mol m^-3)
C <- c(1.05070044705462E-02, 4.69330649964981E-02)
# Growth rate (u in d-1)
u <- 0.79 * ((R*1E+06)^-0.21)
P1 <- 0.01
# Uptake flux of CO2 by the cell (mol m-^2 s^-1)
F1 <- -(23230/3) * R * u / 86400
# Diffusion coefficient of CO2,... (m^2 s^-1)
D <- 1E-09
# Define parameters + grid definition vector
parms <- list(C=C, F1=F1, D=D, P1=P1, X.grid=X.grid, A.grid=A.grid)
#------------------------------#
# Model solution #
#------------------------------#
# Numerical solution at steady state
Cini <- matrix(C, nrow=N, ncol=2, byrow=T)
# I added pos = TRUE
boundary <- steady.1D(y = Cini, func = boundary_layer, pos = TRUE, parms = parms, nspec = 2, names = c('CO2','HCO3'))
#------------------------------#
# Plotting output #
#------------------------------#
# Using S3 plot method of package rootSolve'
plot(boundary, grid = X.grid$x.mid, xlab = 'distance from centre, m', ylab = 'mol/m3', main = c('CO2', 'HCO3'))