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

[Karline Soetaert]

Chris, the problem is that it multiplies the flux with 'A', the surface, and in your case, A at the origin = 0, so the flux vanishes (there is no surface to cross).
Karline

-----Original Message-----
From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of cfriedalek .
Sent: zaterdag 7 februari 2015 9:43
To: r-sig-dynamic-models at r-project.org
Subject: [R-sig-dyn-mod] Radial Diffusion with prescribed flux

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]]

_______________________________________________
R-sig-dynamic-models mailing list
R-sig-dynamic-models at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models