Hello,everyone
I'm new to fipy, and i am learning it using the examples. There is weird
problem happened:
Based on the 1-D diffusion example, i changed the boundary conditions.
I need the surface flux is propotional to the square of surface
concentration: J=Kr*C^2.
So i use fick's law to calculate the diffusion flux on surface: J =
-D*dC/dx.
Then i set the boundary conditions to:
cL.faceGrad.constrain(-Kr/D*cL.faceValue**2, mesh.facesRight)
*The problem is:*
*In line 46, depending on whether i print the "cL.faceGrad[0][0]" or
not, the calculated results are quite different. What happends when the
"cL.faceGrad[0][0]" i**s printed?*
At the same time, i also use the change of total inventory in the volume
to calibrate the flux: J=0.5*d(total inventory)/dt. (i use
cL.cellVolumeAverage*Lx to calculate the total inventory). I found that
to get more precise "faceGrad" on the surface,i need very small mesh
grids, is there better ways to calculate the surface flux?
Thank you very much
---------------------------------------
Code:
from fipy import *
from numpy import *
# Mesh
Lx = 0.5e-3 # Thickness
dd = list(2**arange(20))+[2**20]*2
dd+=dd[::-1]
nx=len(dd)
dx=array(dd)*Lx/sum(dd)
mesh = Grid1D(dx=dx) # 1D Meshing
# CellVariables
cI=1.0e5
cL = CellVariable(name="Lattice concentration", mesh=mesh,value=cI)
# Other Variables and constants
time = Variable() # Time (s)
Temp = 700. # Temperature (K)
k = 8.6e-5 # Boltzmann constant (eV/K)
D = 2.0e-7 * numerix.exp(-0.4/k/Temp) # Diffusion Coeffcient (m2/s)
Kr=3e-25/numerix.sqrt(Temp)*numerix.exp(2.06/k/Temp) # Recombination
rate (m4/s)
Cold=float(cL.cellVolumeAverage*Lx)
Flux = 0
# Boundary conditions
cL.faceGrad.constrain(-Kr/D*cL.faceValue**2, mesh.facesRight) # J =
-D*dcL/dx = Kr*cL**2
cL.faceGrad.constrain(Kr/D*cL.faceValue**2, mesh.facesLeft) # J =
-D*dcL/dx = Kr*cL**2
# Diffusion Equation
eqn = TransientTerm(var=cL) == DiffusionTerm(coeff=D, var=cL)
# Time steps
steps = 500
dt = 500./steps
# concentration plots
if __name__ == '__main__':
viewer = Viewer(vars=(cL), datamin=0., datamax=1.5*cI)
viewer.plot()
# Solve the equation
while time() < steps*dt:
time.setValue(time() + dt)
res = 1e25
print time , cL.faceValue[0], Kr/D*cL.faceValue[0]**2, Flux/D
,cL.faceGrad[0][0]
while res/cL[nx/2] > 1e-5:
res = eqn.sweep(dt=dt, var=cL)
if __name__ == '__main__':
viewer.plot() # update the concentration plots
Cnew=float(cL.cellVolumeAverage*Lx) # Total inventory
Flux = -(Cnew-Cold)/dt/2.
Cold=Cnew
if __name__ == '__main__':
raw_input("Calculation finished. Press <return> to proceed...")
from fipy import *
from numpy import *
# Mesh
Lx = 0.5e-3 # Thickness
dd = list(2**arange(20))+[2**20]*2
dd+=dd[::-1]
nx=len(dd)
dx=array(dd)*Lx/sum(dd)
mesh = Grid1D(dx=dx) # 1D Meshing
# CellVariables
cI=1.0e5
cL = CellVariable(name="Lattice concentration", mesh=mesh,value=cI)
# Other Variables and constants
time = Variable() # Time (s)
Temp = 700. # Temperature (K)
k = 8.6e-5 # Boltzmann constant (eV/K)
D = 2.0e-7 * numerix.exp(-0.4/k/Temp) # Diffusion Coeffcient (m2/s)
Kr=3e-25/numerix.sqrt(Temp)*numerix.exp(2.06/k/Temp) # Recombination rate
(m4/s)
Cold=float(cL.cellVolumeAverage*Lx)
Flux = 0
# Boundary conditions
cL.faceGrad.constrain(-Kr/D*cL.faceValue**2, mesh.facesRight) # J = -D*dcL/dx
= Kr*cL**2
cL.faceGrad.constrain(Kr/D*cL.faceValue**2, mesh.facesLeft) # J = -D*dcL/dx
= Kr*cL**2
# Diffusion Equation
eqn = TransientTerm(var=cL) == DiffusionTerm(coeff=D, var=cL)
# Time steps
steps = 500
dt = 500./steps
# concentration plots
if __name__ == '__main__':
viewer = Viewer(vars=(cL), datamin=0., datamax=1.5*cI)
viewer.plot()
# Solve the equation
while time() < steps*dt:
time.setValue(time() + dt)
res = 1e25
print time , cL.faceValue[0], Kr/D*cL.faceValue[0]**2, Flux/D
,cL.faceGrad[0][0]
while res/cL[nx/2] > 1e-5:
res = eqn.sweep(dt=dt, var=cL)
if __name__ == '__main__':
viewer.plot() # update the concentration plots
Cnew=float(cL.cellVolumeAverage*Lx) # Total inventory
Flux = -(Cnew-Cold)/dt/2.
Cold=Cnew
if __name__ == '__main__':
raw_input("Calculation finished. Press <return> to proceed...")
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
