I'm impressed by the FiPy package, but don't understand how to implement
the convection boundary conditions for my problem.
I can't figure out how to set up the boundary conditions for a steady state
2-D flow problem with a source and radial diffusion. I'd like the inlet on
the left to be at a constant temperature and the outlet on the right to be
what would normally be a constant pressure boundary condition. I want the
top to be at a constant temperature and the bottom to have zero flux.
I have included the code below which I was using for testing and therefore
doesn't include the zero flux BC on the bottom and has several other
inconsistencies, but shows how I'm setting up the problem.
As a side note, I am planning to account for the radial geometry by
multiplying the DiffusionTerm by 1/r and the DiffusionCoeff by r while
using a rectangular grid.
I'd appreciate any input you can provide.
"""
This program solves a 2-d Cylindrical problem with radial conduction and
axial convection with a source.
(k/r)*d/dr(r.dT/dr) - rho.v.cp.dT/dx + q = 0
Boundary conditions:
T(x,R) = Tw
T(0,r) = Tin
dT/dr|r=0 = 0
The system is at steady state.
"""
from fipy import *
R = 0.005 # m, radius of tube
L = 0.1 # m, length of tube
nr = 10 # [], number of radial cells
nx = 10 # [], number of axial cells
k = 0.05e-10 # W/m.K, radial thermal conductivity
rho = 0.5 # kg/m^3, fluid density
v = 1.e1 # m/s, fluid velocity
cp = 1.0 # J/kg.K, fluid heat capacity
q = 0.2e-20 # W/m^3, heat generation
Tw = 5000. # C, wall temperature
Tin = 0. # C, inlet temperature
# FiPy convention conversion, x -> x, y -> r
# T(x, r)
dx = L/float(nx)
dr = R/float(nr)
mesh = Grid2D(dx=dx, dy=dr, nx=nx, ny=nr)
coords = mesh.getCellCenters()
x, y = coords
#print y
# create cell variable and initialize
phi = CellVariable(name = "Temperature", mesh = mesh, value = 0.)
# boundary conditions
# FixedValue at y=R, FixedFlux at y=0
BCs = (FixedValue(faces = mesh.getFacesTop(), value = 1000),\
FixedValue(faces = mesh.getFacesBottom(), value = 2000), \
FixedValue(faces = mesh.getFacesLeft(), value = 3000))
u = ((v,),(0,))
# outflow condition in right side
RHSBC = (u*mesh.getFacesRight()).getDivergence()
# Define diffusion coefficient as vector to only use radial term
D = [[[0.,0.],[0.,k]]]
# Define source term
S = q
# The direction of the velocity vector is only OK if the ConvectionTerm
# is positive although it should really be negative.
eq = (1.)*DiffusionTerm(coeff = D) \
+ ExponentialConvectionTerm(coeff = u) \
+ ImplicitSourceTerm(RHSBC)
# + S \
#eq.solve(phi,boundaryConditions
#print BCs
eq.solve(phi,boundaryConditions=BCs)
print phi
viewer = Viewer(vars = phi)
viewer.plot()
raw_input("finished")
_______________________________________________
fipy mailing list
[email protected]
http://www.ctcms.nist.gov/fipy
[NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy]
