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

[Karline Soetaert]

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]]
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>

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