Thanks for the reply. Those are very useful suggestions. I have some new
findings which may help shedding some light on this problem.
I have added TransientTerm and tried to see how the solution is developed
over time. I also remove the datamin and datamax input on the viewer, so
that I can see the full range of the solution variable data.
When I see the temporal development of the solution (attached in the
email), It looks stable. However, I noticed that the right outlet boundary
has extreme high concentration and did not outflow the flow domain. I think
It has due to the improper setup of my right outlet condition. Since the
materials keep building up on the right outlet boundary, it appears to
affect the solution in the actual flow domain(also attached in email) and
drives the solution unstable. This could be the main reason that solving
steady state equation directly never converges.
I then begin to look into how to properly set up "open outlet" boundary,
the closest I could find is following:
https://www.mail-archive.com/fipy@nist.gov/msg02797.html
and in the section of "Applying fixed flux boundary conditions" of official
website
http://www.ctcms.nist.gov/fipy/documentation/USAGE.html
I tried to adapt the official method listed on the official website, but no
success is achieved. Thus, I wonder is it possible to point out some
examples that used this method to define fixed flux outlet boundary
condition ? You may find my current progress on this part in the code I
attached.
Another thing I noticed is that, there is a reported issue on robin
boundary condition:
https://github.com/usnistgov/fipy/issues/426
Would it be possible to explain more about this ? I don't think the up &
bottom boundary condition in my case causes the non convergent problem,
because from the Transient equation, the segregation of the material starts
to develop at the top&bottom boundary, which is consistent with my
expectation.
I have attached my current progress of the code, and the temporal
development of the solution.
Best,
Zhekai
On Thu, Sep 15, 2016 at 9:33 AM, Daniel Wheeler <daniel.wheel...@gmail.com>
wrote:
> On Wed, Sep 14, 2016 at 7:52 PM, Guyer, Jonathan E. Dr. (Fed)
> <jonathan.gu...@nist.gov> wrote:
> >
> > I don't know offhand. With a Péclet number of 100, your problem is
> almost completely hyperbolic, which FiPy (and cell-centered Finite Volume)
> isn't very good at. Daniel knows more about this and may have some
> suggestions.
> >
> > You might consider adding a TransientTerm for artificial relaxation and
> trying the VanLeerConvectionTerm.
>
> Definitely use a TransientTerm and see the time step you can get away
> with. Don't sweep, just use time steps. Using the VanLeerConvetionTerm
> isn't necessary as the accuracy is not important in a steady state
> problem.
>
> --
> Daniel Wheeler
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
from fipy import *
# mesh set up
nx = 100.
nz = 100.
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]]]
#This is the new part to account for outflow fixed flux boundary (this does not
set up correctly)
outflowflux = FaceVariable(mesh = mesh, value = phi.faceValue*velocityVector)
velocityVector.constrain(0.,mesh.facesRight)
# equation set up
eq = (TransientTerm() + UpwindConvectionTerm(coeff = velocityVector)\
+ UpwindConvectionTerm(coeff =local_shear_rate*(1-phi).faceValue)\
- (mesh.facesRight * outflowflux).divergence \
== 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.constrain(0, where = mesh.facesRight) #not working outlet boundary
condition
#phi.faceGrad.constrain(velocityVector*phi.faceValue,where = mesh.facesRight) #
problematic outlet boundary conditions
# Viewer set up
viewer = Matplotlib2DViewer(vars=phi, title="final solution")
# Begin the calculation of the pde
timeStepDuration = 10*0.9*dx**2/(2*0.1)
steps = 5000
for step in range(steps):
eq.solve(var = phi,dt = timeStepDuration)
viewer.plot()
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
