mdfnewtondiffu: 1D diffusion, boundary constant, FD newtonian method,...

In ahmad-alkadri/Rdiffsolver: A simple R package for solving diffusion equation with numerical methods.

Description

Function to solve diffusion equation using finite difference method, backward euler or implicit scheme, applied in the condition where: a) the diffusion occurs inside a slab of material unidimensionally, or only along one axis, b) the coefficient of diffusion is dependent on the concentration (D = f(c)), and c) the boundary conditions are constant.

Usage

mdfnewtondiffu(Nt,Nx,T,C_f,C_i,l,Dfunstr)

Arguments

Nt	number of points in time (t)
Nx	number of points in space (x)
T	Total calculated diffusion time
C_f	Dirichlet boundary condition, final concentration coming from the outside of the slab, with one or two elements. If there is only one element, the two sides of the slab (x = -l and x = l) will have the same C_f. If there are two elements, the first element (C_f[1]) will be on x = -l while the second one will be on x = l.
C_i	Initial concentration value inside the slab
l	half-length of the slab, usually in cm
Dfunstr	the coefficient of diffusion written as a function of concentration (c) as string; please check the example for details

Value

A matrix with Nx number of row and Nt number of column, profiling the diffusion.

Examples

C_i = 0.00 # Initial concentration inside the slab
C_f = c(0.10,0.07)	# Final concentration coming from outside
Nt = 30 # points in time
Nx = 20 # points in space
l = 0.25 # half-thickness of the slab, in cm
T = 432000 # Total measured time in seconds (~5 days)

# Coefficient of diffusion as a linear function of c
penta = 2.21e-7 #cm^2/s
bro = 1.91e-7 #cm^2/s
Dfunstr = paste("function(c){",penta,"*c+",bro,"}",sep="")
print(Dfunstr)

# Calculation
U <- mdfnewtondiffu(Nt,Nx,T,C_f,C_i,l,Dfunstr)

# Visualization
x_s = seq(-l,l,length.out=Nx)
y_s = seq(0,T,length.out=Nt)
persp(x=x_s, y=y_s, z=U,
      theta = -35, phi = 25,
      ticktype = "detailed",
      xlab = "x (mm)",
      ylab = "t (seconds)",
      zlab = expression(c[x,t]))

ahmad-alkadri/Rdiffsolver documentation built on Feb. 4, 2020, 9:45 p.m.