Hi all,
I am trying to use Fipy to solve convection only problem for the
concentration moved only by solid body rotation in a "circular" shape
geometry.
By looking at the examples online, I found out that
http://www.ctcms.nist.gov/fipy/examples/convection/
generated/examples.convection.source.html#module-examples.convection.source
and some of the level set example appears to allow me do it.
I implemented the approach from convection example. However, the solution
still looks has diffusion ( or maybe artificial smoothness ) as the
concentration move with velocity field. I have attached my example code,
that concentration enter from the top right side of the geometry, and
undergo solid body rotation eventually to the left side, and flow out of
the domain. So my question is that is there any way to further reduce the
diffusion? Also, does anyone know where this "diffusion" is coming from ?
The approach that I have tried but did not work are following:
1. Solve the equation with very small diffusion coefficient (1e-8)
2. Reduce the timestep or refine the mesh size does not seem to help very
much
I have attached my example code in this email. Thank you very much.
Best,
Zhekai
from fipy import *
import numpy as np
import matplotlib as matplotlib_plot
mesh_2 = Gmsh2D('''
lc = 1e-2;
Point(1) = {-1, 0, 0, lc};
Point(2) = {0, 0, 0, lc};
Point(3) = {1, 0, 0, lc};
Point(4) = {0,-0.1,0,lc};
Ellipse(5) = {1,2,2,4};
Ellipse(6) = {4,2,2,3};
Line(7) = {3,2};
Line(8) = {2,1};
Circle(9) = {1,2,3};
Line Loop(9) = {6,7, 8,5};
Line Loop(10) = {-6,-5,9};
Plane Surface (12) = {10};
Physical Surface("bulk_flow") = {12};
Physical Line("Bottom_Left") = {5};
Physical Line("Bottom_Right") = {6};
''' % locals())
phi_2 = CellVariable(name = "solution variable",
mesh = mesh_2,
value = 0.) # doctest: +GMSH
xFace_2,yFace_2=mesh_2.faceCenters
velocityX_2 = yFace_2
velocityY_2 = -xFace_2
velocityVector_2 = FaceVariable(mesh=mesh_2, rank=1)
velocityVector_2[0] = velocityX_2
velocityVector_2[1] = velocityY_2
exit_velocityX_2 = yFace_2
exit_velocityY_2 = -xFace_2
exit_velocityVector_2 = FaceVariable(mesh=mesh_2, rank=1)
exit_velocityVector_2[0] = exit_velocityX_2
exit_velocityVector_2[1] = exit_velocityY_2
out_flow_mask_2 = ((xFace_2 <= 0.) & (yFace_2 < 0.))
exteriorCoeff_2 = FaceVariable(mesh_2, value=exit_velocityVector_2, rank=1)
exteriorCoeff_2.setValue([[0.],[0.]], where= ~(mesh_2.exteriorFaces &
mesh_2.physicalFaces["Bottom_Left"]))
velocityVector_2.setValue([[0.],[0.]], where= mesh_2.exteriorFaces &
mesh_2.physicalFaces["Bottom_Left"])
eq = TransientTerm() \
+ PowerLawConvectionTerm(coeff = velocityVector_2)\
+ ImplicitSourceTerm(exteriorCoeff_2.divergence) == 0
phi_2.constrain(1., where = mesh_2.exteriorFaces &
mesh_2.physicalFaces["Bottom_Right"]) # doctest: +GMSH
viewer_2 = Matplotlib2DViewer(vars=phi_2, title="final solution", cmap =
matplotlib_plot.cm.hot,datamin = 0., datamax = 1.)
timeStepDuration = 1/200.
steps = 800
for step in range(steps):
eq.sweep(var=phi_2,dt=timeStepDuration) # doctest: +GMSH
viewer_2.plot()___
fipy mailing list
fipy@nist.gov
http://www.ctcms.nist.gov/fipy
[ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
