Zhekai -
1. I think what you have is fine. Explicitly, I'd write:
phi.faceGrad.constrain(..., where=mesh.facesTop | mesh.facesBottom)
See
http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-spatially-varying-boundary-conditions
2. For inlet/outlet conditions, please see the comment at
http://www.ctcms.nist.gov/fipy/documentation/USAGE.html#applying-outlet-or-inlet-boundary-conditions
If you impose constraints, then FiPy is inlet/outlet, not no-flux.
I've made some changes to your script that result in the initial sweeps looking
much more like your matlab solution, but it is still not convergent. I wonder
at this point if a velocity is backwards someplace.
I don't understand why $\tilde{\dot{\gamma}} = 1 \hat{\j}$.
Shouldn't it be $\tilde{\dot{\gamma}} = 1 \hat{\i}$, given that $\tilde{u} =
\tilde{z} \hat{\i}$?
Is, perhaps, the definition
$\tilde{\dot{\gamma}} \equiv \frac{\partial \tilde{u}_x}{\partial \tilde{z}}
\hat{\j}$
where
$\tilde{u} \equiv \tilde{u}_x \hat{\i} + \tilde{u}_y \hat{\j}$
(I *really* strongly recommend not writing vectors like scalars)
- Jon
> On Sep 13, 2016, at 5:45 PM, Zhekai Deng <[email protected]>
> wrote:
>
> 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
> <Advection_Diffusion_Segregation.py><matlab_demo.m><equation_explain.pdf>_______________________________________________
> fipy mailing list
> [email protected]
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
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
D_Matrix = [[[0, 0],
[0, 1/Pe]]]
eq = (ExponentialConvectionTerm(coeff = velocityVector) + ExponentialConvectionTerm(coeff = shear_rate*(1-phi).faceValue) == DiffusionTerm (coeff=D_Matrix))
# 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
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!
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
