1D Diffusion Equation: Boundary Condition

[beantsja]

Solving a 1D diffusion equation:

with the boundary condition:

I implemented the equation of the first picture, but I want to make sure about setting up the boundary equation.

eqI = TransientTerm() == DiffusionTerm(coeff=D) after setting the coefficient as D,

X.faceGrad.constrain(bcRight * mesh.faceNormals, mesh.facesRight) for this line, does bcRight should be J or J/D? Do I have to bring the diffusion coefficient to the RHS?

[guyer]

X.faceGrad is the gradient of X, or $\partial C_A/\partial x$, so constrain it to $J\cdot\hat{n}/D$. Alternatively, constrain the flux at the boundary, in which case, exteriorFlux would be $\vec{J}$.

[beantsja]

`from fipy import CellVariable, Grid1D, TransientTerm, DiffusionTerm, Viewer
import math
import matplotlib.pyplot as plt

da = 5e-5 # 1D length.
nx = 100 # Spatial Grid Points.
dx = da / nx # Each Spatial Step Length in Domain.
initial_X = ((1)*(5e-5))/(nx)
D = 1e-7
timeStepDuration = 0.00001
step = 10 #timeStepDuration * second
X_values = []
exteriorFlux = 1e-7
bcLeft = 0
mesh = Grid1D(nx=nx, dx=dx) # Creating the Grid in nx points with dx step length.
X = CellVariable(name="X", mesh=mesh, value=float(initial_X))

D.constrain(0., mesh.exteriorFaces)
X.faceGrad.constrain(bcLeft * mesh.faceNormals, mesh.facesLeft)

eqI = (TransientTerm() == DiffusionTerm(coeff=D) + (mesh.exteriorFaces * exteriorFlux).divergence)

for i in range(step):
eqI.solve(var=X, dt=timeStepDuration)
X_values.append(float(X[-1]))`

Here I am trying to implementing the situation in my original problem. Instead of using the FaceGrad, I am trying to use the alternative way from your first reply (exterior.Flux). "D.constrain(0., mesh.exteriorFaces)" this line is giving me an error with AttributeError: 'float' object has no attribute 'constrain'.

What would be the issue?

[guyer]

D must be a FaceVariable in order to constrain its values at mesh.exteriorFaces.

D = FaceVariable(mesh=mesh, name="D", value=1e-7)
D.constrain(0., mesh.exteriorFaces)

[beantsja]

from fipy import CellVariable, Grid1D, TransientTerm, DiffusionTerm, Viewer
import math

da = 5e-5   # 1D length. cm
nx = 25 # Spatial Grid Points. unitless
dx = da / nx  # Each Spatial Step Length in Domain. cm
initial_X = (2e-6)#((1)*(5e-5))/(nx) #mole
D = 1e-7 
timeStepDuration = 1
step = 100 #timeStepDuration * second
flow = 1e-6
A = 5
n= 2
F =96485
X_values = []


mesh = Grid1D(nx=nx, dx=dx)  # Creating the Grid in nx points with dx step length.
X = CellVariable(name="X", mesh=mesh, value=float(initial_X))  
eqI = TransientTerm() == DiffusionTerm(coeff=D)
bcRight = -((flow*A)) / (n * F*D)
flux = bcRight*D
bcLeft = 0

Xi=sum(X)

X.faceGrad.constrain(bcRight* mesh.faceNormals, mesh.facesRight)
X.faceGrad.constrain(bcLeft * mesh.faceNormals, mesh.facesLeft)

for i in range(step):
    eqI.solve(var=X, dt=timeStepDuration)
    X_values.append(float(X[-1]))
    if sum(X)<0:
        # print(i)
        break

Xf=sum(X)
print("flux:", flux)
print("Total mole:",Xi)
print((Xi-Xf)) #Total mol Change
print(flux*(i+1)) #Total Flux Amount
print((Xi-Xf)*(dx)) #Unit Conversion

Here, I am trying to check the total concentration change is equal to the total flux out. I expect print((Xi-Xf)) #Total mol Change to be equal to print(flux*(i+1)) #Total Flux Amount. However, it is only equal when I do this unit conversion. print((Xi-Xf)*(dx)) #Unit Conversion. I am confused with how Fipy takes the input of the initial variables at X = CellVariable(name="X", mesh=mesh, value=float(initial_X)). There must be something I am misunderstanding how Fipy is taking control of the variables, and I believe it happens when I locate the initial variable values in each index when mesh is created. But, I am not quite sure what is going on, and why the unit conversion above is necessary.

[guyer]

The amount of "stuff" in a cell is the value at the cell center times the volume of the cell. You need to compare the before and after amount of substance (extensive) rather than the concentration (intensive). The preferred way to do that in FiPy (which accounts for nonuniform meshes) is

Xi = (X.cellVolumeAverage * mesh.cellVolumes.sum()).value
  :
  :
Xf = (X.cellVolumeAverage * mesh.cellVolumes.sum()).value

X.cellVolumeAverage represents $\sum_\text{cells} (X V_\text{cell}) / \sum_\text{cells} V_\text{cell}$.

note: if you do not use .value, then FiPy's lazy evaluation will result in Xi having the same value as Xf. You need .value to hold the value at that time.

[beantsja]

Thank you. The volume of the cell is (the length of the index)^3. Right?

[guyer]

In 1D, the volume of the cell is ${\tt dx} \times 1 \times 1$

[beantsja]

When I print X, the whole grids, it shows the X values in each index. However, what matters is the X times volume in Fipy, right? Is there a way to print X, but present X*volume values instead of just X? It is because I want to see the results of the actual amount of stuffs in the cell, rather than the value X.

Another thing I want to confirm is that X*volume is the value that Fipy is using to solve the diffusion equation. Is that correct?

[guyer]

X * mesh.cellVolumes

Another thing I want to confirm is that X*volume is the value that Fipy is using to solve the diffusion equation. Is that correct?

Yes. Please see the documentation on the discretizations that FiPy uses.