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

[cfriedalek .]

Thanks again Karline. I'm on vacation now so I'll try this out when I get
back. Hopefully I can post a successful result to this list in a couple of
weeks.
Regards Chris

On Thu, Feb 12, 2015 at 10:33 PM, Karline Soetaert <Karline.Soetaert at nioz.nl
> 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]]
> >
> > _______________________________________________
> > 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
>
> _______________________________________________
> 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]]