Thank you for the reply. I try to answer any confusion in below.
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
Thanks. It has been implemented in this way now.
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.
So, in my previous one, I only have constraints on left (fixed 0.5 inlet
concentration), top and bottom ( the boundary conditions that I applied),
but not boundary condition on the right. I have implemented an additional
one on the right to account for the outlet:
phi.faceGrad.constrain(velocityVector*phi.faceValue,where = mesh.facesRight)
but it seems does not help very much.
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.
Thank you for the change. Well, even after I added the right side outlet
boundary condition, the sweep solution still looks similar to what you
have. In the first initial sweeps, it looks like the matlab solution, but
it is not convergent. However, the residual seems does not explore as
quickly as before by adding the outlet condition.
I am not sure what does "velocity is backwards someplace" mean. From the
way I defined it, the $u \hat{i} = \tilde{z}$ and $u \hat{j} = 0$. Both are
independent of solution variable. Would it be possible to clarify what does
"velocity is backwards" mean?
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}$?
Thanks for bring it up. After the discussion with my colleague, we think it
is better to understand this as a vector term that only act on the normal
direction ($\hat{j}$) of the flow domain. The reason for this is that the
shear rate term is gradient of velocity field which is a tensor. However,
we are not entire sure the functional dependence of segregation velocity on
shear rate tensor. That's why we usually call it as local shear rate (du/dz
\hat{j} ), and I should have been clear on this issue in the beginning.
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}$
Yes, this is the definition we are currently following.
In summary, it does seems that the initial sweeps give the right "trend"
toward the solution. However, I am still not sure about why it is not
convergent. Any suggestion on how to proceed to the next step ? I have
attached the newest version of the code.
Best,
Zhekai
On Tue, Sep 13, 2016 at 6:35 PM, Guyer, Jonathan E. Dr. (Fed) <
[email protected]> wrote:
> 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 <zhekaideng2014@u.
> northwestern.edu> 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 ]
>
> _______________________________________________
> 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
local_shear_rate_X = FaceVariable(mesh=mesh, value = 0) # this term only acts
on normal direction
local_shear_rate_Y = FaceVariable(mesh=mesh, value = 1) # linear velocity
profile du/dz
local_shear_rate = FaceVariable(mesh=mesh, rank=1)
local_shear_rate[0] = local_shear_rate_X
local_shear_rate[1] = local_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 =local_shear_rate*(1-phi).faceValue) ==
DiffusionTerm (coeff=D_Matrix))
# set up the boundary conditions
phi.constrain(0.5, where = mesh.facesLeft) # fixed the inlet as 0.5
concentration
phi.faceGrad.constrain(local_shear_rate*Pe*(1-phi.faceValue)*phi.faceValue,
where = mesh.facesTop | mesh.facesBottom) #boundary conditions at the top and
bottom
phi.faceGrad.constrain(velocityVector*phi.faceValue,where = mesh.facesRight) #
outlet boundary conditions
# 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
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
