[R-sig-dyn-mod] Radial Diffusion with prescribed flux

[cfriedalek .]

Hi

Thanks for ReacTran and deSolve. To help my understanding I'm trying to
simulate a line source heating problem in 1-D. Heat is generated along the
central axis of an infinitely long cylinder of finite radius and diffuses
radially to the exterior wall of the cylinder of radius R. I'm modelling
this as a 1-D line along the radial (diffusion) direction with a
prescribed, positive heat flux at the origin and an insulated exterior wall
(zero heat flux). My code is shown below.

The problem I'm having is that when I set a positive heat flux at r=0 and
zero heat flux at r=R the solver returns zero for the temperature (C) field
at all times.
However if I set a negative heat flux at r=R (heat flow into the cylinder)
and zero at r=0 I get a sensible result (the wall temperature rises in time
and eventually so too does the center temperature). I don't understand why
this doesn't work for the case with heating in the center.

Have I set this up properly?

I'm really enjoying working with R and R-studio so even though I have
solved this problem with other programming languages, I would really like
to migrate to working with R. Getting my head around this is one step of
the way. Ultimately I want to extend this to modelling a real line-source
device with zones of different materials and electrical heating in the
core, so any help is greatly appreciated.

regards

Chris

# 1-D radial diffusion
library(ReacTran)
library(deSolve)
N <- 10
xgrid <- setup.grid.1D(N = N, L = 1)
r1 <- setup.prop.1D(grid = xgrid, func = function(r) r)
Yini <- rep(0, N)
times <- seq(from = 0, to = 10, length.out = 101)
D.coeff <- 0.05
radial_diffusion <- function (t, Y, parms, A){
  tran.1D(C = Y, flux.up = flux.up, flux.down = flux.down, D = D.coeff, A =
A, dx = xgrid)
}
# heating from outside - THIS WORKS
flux.up <- 0
flux.down <- -1
# # heating from inside - THIS FAILS
# flux.up <- 1
# flux.down <- 0
# solve
out1 <- ode.1D(y = Yini, times = times, func = radial_diffusion, parms =
NULL, nspec = 1, A = r1)
# Plot Results
par (mfrow=c(1, 2))
# time evolution at center "1" and wall "N"
plot(out1[1:nrow(out1), "time"], out1[1:nrow(out1), as.character(N)], type
= "l", lwd = 2,
     xlab = "time", ylab = "C", main = "Time Evolution")
lines(out1[1:nrow(out1), "time"], out1[1:nrow(out1), as.character(1)])

# radial distribution every 10 time steps
plot(xgrid$x.mid, out1[1, 2:(N+1)], type = "l", lwd = 2, xlab = "r/R", ylab
= "C",
     main = "Radial Distribution", ylim=c(0, max(out1[1:nrow(out1),
as.character(N)])))
time_indicies <- seq(1, length(times), by=10)
for (i in time_indicies) lines(xgrid$x.mid, out1[i, 2:(N+1)])
par (mfrow=c(1, 1))

	[[alternative HTML version deleted]]