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

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

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

wrote:

What I would do, is to implement this heat source as a regular source
term, outside of the transport routine. Of course you then need to probably
divide by the thickness of the first cell.

So,

dY <- tran.1D(C = Y, flux.down = flux.down, D = D.coeff, A = A, dx =
xgrid)$dC  # No flux.up!

dY[1] <-dY[1] + flux.up/xgrid$dx.int[1]   # something like that...

 Karline

-----Original Message-----
From: R-sig-dynamic-models [mailto:
r-sig-dynamic-models-bounces at r-project.org] On Behalf Of cfriedalek .
Sent: dinsdag 10 februari 2015 11:25
To: Special Interest Group for Dynamic Simulation Models in R
Subject: Re: [R-sig-dyn-mod] Radial Diffusion with prescribed flux

Thanks Karline. So I think I need to consider the heated central zone as
an infinitely long cylinder of finite radius rather than a line of zero
radius. This would mimic the real situation of a heated slender wire in a
domain.

I understand how I can specify a two zone grid where, say, the inner 1% of
the total radius is the "wire" as follows:

xgrid <- setup.grid.1D(x.up = 0, x.down = c(0.01, 1), L=c(0.01, 0.99), N =
c(1, 99))

Then setup the properties as before.

For constant diffusion coefficient of value 1:

r1 <- setup.prop.1D(grid = xgrid, func = function(r) r)

For different diffusion coefficients for wire and domain, say, 10 and 1
respectively:

D.coeffs <- c(10, 1)
wire_radius <- 0.01
r1 <- setup.prop.1D(grid = xgrid, func = function(r) r*D.coeffs[[(r >
wire_radius)+1]])


However now I am unclear how to proceed.

How do I specify a heat source that applies only to the inner zone of the
grid? The heat source is a electrical resistance wire so the power
generation is Volt^2/Resistance (Watts). I guess that must be done in the
transport function but I'm not sure how. Also, do I need to specify the D
parameter since it appears in the properties?

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) # Add heat generation term for the wire zone }

I hope you can help.

Regards

Chris

PS: how to a achieve nice code layout in these posts. Compared to other
posts it looks like I'm missing some part of CR/LF. I'm using R on Windows.

On Tue, Feb 10, 2015 at 5:50 AM, Karline Soetaert <
Karline.Soetaert at nioz.nl>
wrote:

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

        [[alternative HTML version deleted]]