Hello All, I am trying to use Fipy to solve advection, segregation and diffusion equation in 2D. I have attached the pdf to explain in details about my equation, and corresponding fipy code that I have problem with it.
The solution I got has "kind of similar shape" to the same equation that I solved using MATLAB(the code is also attached, I have verified that MATLAB is solving correctly). However, the code I have in Fipy does not converge, and the residual jump between large and small number. Also,I have a few things that I am not entire sure. 1. In 2D, how could we specify the which direction of the "faceGrad.constrain" that we want to constrain ? 2. In my problem specifically, I would like to implement convection inlet flux and convection outlet flux. This seems is against the Fipy default no flux condition. Any suggestion to on how to implement it ? Best, Zhekai
from fipy import *
# mesh set up
nx = 300.
nz = 300.
dx = 1/nx
dz = 1/nz
mesh = Grid2D(dx=dx, dy=dz, nx=nx, ny=nz)
phi = CellVariable(name = "solution variable",mesh = mesh, value = 0.5)
# the follwing set up the velocity and shear rate in vector form.
# for those direction that I don't need it, I set its value as zero.
velocityX = mesh.faceCenters[1] # linear velocity profile
velocityY = FaceVariable(mesh=mesh, value = 0) # zero velocity in y direction
velocityVector = FaceVariable(mesh=mesh, rank=1)
velocityVector[0] = velocityX
velocityVector[1] = velocityY
shear_rate_X = FaceVariable(mesh=mesh, value = 0) # linear velocity profile
shear_rate_Y = FaceVariable(mesh=mesh, value = 1) # zero velocity in y direction
shear_rate = FaceVariable(mesh=mesh, rank=1)
shear_rate[0] = shear_rate_X
shear_rate[1] = shear_rate_Y
Pe = 100. # some constants
#Setting up the diffusion part
Diffusion_X = FaceVariable(mesh=mesh, value = 0)
Diffusion_Y = FaceVariable(mesh=mesh, value = 1/Pe)
D_Vector = FaceVariable(mesh=mesh, rank=1)
D_Vector[0] = Diffusion_X
D_Vector[1] = Diffusion_Y
# mask represents where the boundary condition is applied
mask = mesh.facesTop | mesh.facesBottom
phi.constrain(0.5, where = mesh.facesLeft) # fixed the inlet as 0.5
concentration
phi.faceGrad.constrain([shear_rate*Pe*(1-phi.faceValue)*phi.faceValue],mask)
#boundary conditions at the top and bottom (Not sure whether it is correct or
not)
# Viewer set up
viewer = Matplotlib2DViewer(vars=phi,datamin=0.,datamax = 1., title="final
solution")
# Begin the calculation of the pde
res = 1e+10
while res >1e-6: #wait until it converges
eq = (ExponentialConvectionTerm(coeff = velocityVector) +
ExponentialConvectionTerm(coeff = shear_rate*(1-phi).faceValue) ==
DiffusionTerm (coeff=D_Vector))
res = eq.sweep(var = phi)
viewer.plot()
print res
#I implmented both advection and segregation term using
ExponentialConvectionTerm.
#However, I would want there is a inlet flux on the left, and outlet flux on
the right due to only advection.
#Any suggestions on how to implement is very appreciated!
matlab_demo.m
Description: application/vnd.wolfram.mathematica.package
equation_explain.pdf
Description: Adobe PDF document
_______________________________________________ fipy mailing list [email protected] http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
