Hi Zhekai,
Try adding a `TransientTerm()` to the equation and just using the
`faceCenters` as the velocity so it looks like
eq = (TransientTerm() == DiffusionTerm(coeff=D) -
ExponentialConvectionTerm(coeff = mesh.faceCenters))
and making the time step large. I used 1e4 to test. Using the
faceCenters ensures that the inlet/outlet condition has the same value
as the mesh is resized. FiPy doesn't interpolate those values near the
boundary from cell to face. I think that using a TransientTerm() makes
the solution unique. It's probably always a good idea to use a
TransientTerm() anyway and it never hurts to do so.
Cheers,
Daniel
On Tue, Mar 8, 2016 at 1:04 PM, Zhekai Deng
<zhekaideng2...@u.northwestern.edu> wrote:
> Dear all,
>
> I have a probably simple problem: I tried to solve the 2D
> Convection-Diffusion problem with a velocity field that is spatially
> dependent. I initialize my solution variable, phi, as 0.5. I found out that
> If I don't implement the inlet boundary condition on the left (Dirichlet on
> the left, on the line 24) and just impose velocity field, the solution
> changes a lot when I increase mesh size. However, If I implement inlet
> boundary condition, the solution seems stays the same as the mesh size
> increases. I am not completely sure why mesh sizes will play a role in the
> first case. Could someone clarify this problem to me ? Thank you very much.
> I have copied and pasted the code in the following, and attached the source
> code in the email.
>
> -----------------------------------
>
> from fipy import *
>
> # when I change the nx from 400 to 200 and to 100 and ny from 400 to 200 and
> to 100.
>
> # The solutions changes a lot. I don't know exactly why the mesh would
> change
>
> # the solution this much.
>
> nx = 200.
>
> ny = 200.
>
> dx = 1/nx
>
> dy = 1/ny
>
> L = 1.
>
> mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
>
> #initial phi is equal to 0.5 in all the space
>
> phi = CellVariable(name = "solution variable",mesh = mesh, value = 0.5)
>
>
> # the follwing set up the velocity in vector form.
>
> # for those direction that I don't need it, I set its value as zero.
>
> velocityX = CellVariable(value = mesh.y,mesh=mesh, name=r'$u_x$') # linear
> velocity profile
>
> velocityY = CellVariable(value = 0,mesh=mesh, name=r'$u_y$') # zero velocity
> in y direction
>
> velocityVector = CellVariable(mesh=mesh, name=r'$\vec{u}$', rank=1)
>
> velocityVector[0] = velocityX
>
> velocityVector[1] = velocityY
>
> D = 0.1
>
>
> #This is the inlet boundary conditions that I tried to implement.
>
> #phi.constrain(1, where = mesh.facesLeft)
>
> viewer = Matplotlib2DViewer(vars=phi, title="final solution")
>
>
>
> eq = (0 == DiffusionTerm(coeff=D)
>
> - ExponentialConvectionTerm(coeff = velocityVector.faceValue))
>
> eq.solve(var = phi)
>
> viewer.plot()
>
> --------------------------------------------
>
>
>
> Thank you very much!
>
> Best,
>
> Zhekai
>
>
>
> _______________________________________________
> fipy mailing list
> fipy@nist.gov
> http://www.ctcms.nist.gov/fipy
> [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
>
--
Daniel Wheeler
_______________________________________________
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
