Dear All.
I am trying to solve a heat transfer equation with a conductive and convective term with the form of
k*grad^2*T-Cp*u*grad(T)=0
over a cube space with 2 boundary conditions, a Dirichlet boundary on one surface and a numan boundary conditions on an other where -dT/dx=Q(T-Text)
Here k is conductive heat transfer coefficient, Cp is heat capacity and u is velocity profile defined by an function for laminar flow so I do not solve Navia stokes here.
The problem above is linear, but my final goal is to solve a non-linear problem (as Cp and k are functions of T) so the example code will use a nonlinera solver in it, and the function is written as such as well.
I did integration by parts on the first term of the function above and got the following variational problem as written in Fenics UFL format
Here T is the trial function, and v is the test function
F=inner(k*grad(T),grad(v))*dx-inner(Cp*dot(grad(T),u),v)*dx+Q*T*v*ds(1)-Q*T_ext*v*ds(1)
Here T_ext, Cp, and k are constants, and the dot product is needed as velocity function is a vector as well as grad(T) should be a vector, and needs to be multiplied out, and ds(1) is the surface perpendicular to the flow direction.
I am able to solve this now with the code bellow, and would appreciate if some one could double check my implementation for any errors. Thank you.
from dolfin import *
# Create mesh and define function space
class lf_function(Expression):
def setup(self,FH=0.25,FL=0.25,Fxoff=0,Fyoff=0):
self.l=FH #y length
self.L=FL #x length
self.A=self.l/self.L
self.n=100
self.xOffset=Fxoff
self.yOffset=Fyoff
def m_0(self,n):
a=lambda n:(128.0*(-1.0)**n)/(self.A*pi**5.0*(2.0*n+1.0)**5.0)
b=lambda n:np.tanh((2.0*n+1.0)*pi*self.A/2.0)
c=lambda n:np.sin((2.0*n+1.0)*pi/2.0)
sumabc=lambda n:a(n)*b(n)*c(n)
return -2.0/3.0+np.sum(sumabc(n))
def U(self,x):
n=np.array(range(self.n))
X=(x[0]-self.xOffset)/self.L
Y=(x[1]-self.yOffset)/self.l
sumTop =lambda n:32.0*((-1.0)**n)*np.cosh((2.0*n+1.0)*pi*self.A*Y/2.0)
sumBottom=lambda n:((2.0*n+1.0)**3.0)*(pi**3.0)*np.cosh((2.0*n+1)*pi*self.A/2.0)
sumMultiplier=lambda n:np.cos((2.0*n+1)*pi*X/2.0)
sumTotal=lambda n:(sumTop(n)/sumBottom(n))*sumMultiplier(n)
m_0=self.m_0(n)
return (1.0/m_0)*((X**2.0-1.0)+np.sum(sumTotal(n)))
def eval(self, value, x):
if x[0]<=self.L+self.xOffset and x[0]>=self.xOffset-self.L \
and x[1]<=self.l+self.yOffset and \
x[1]>=self.yOffset-self.l:
value[2]=self.U(x)
value[0]=0
value[1]=0
else:
value[2]=0.0
value[1]=0.0
value[0]=0.0
def value_shape(self):
return (3,)
cellMesh = UnitCubeMesh(5, 5,5)
V = FunctionSpace(cellMesh, "Lagrange", 2)
# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
return x[2] < DOLFIN_EPS
class LowerRobinBoundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
#if abs(x[1]) <= tol:
#print x[0],x[1],x[2]
return on_boundary and abs(x[1]) <= tol
# Define DiricletBC condition
u0 = Constant(300)
bc = DirichletBC(FunctionSpace(cellMesh, "Lagrange", 2), u0, boundary)
bcs=[bc]
# Defining Numan Boundary
boundary_parts = \
MeshFunction("size_t", cellMesh, cellMesh.topology().dim()-1)
boundary_parts.set_all(0)
Gamma_R = LowerRobinBoundary()
Gamma_R.mark(boundary_parts, 1)
uflow=lf_function(degree=2)
uflow.setup(FH=0.5,FL=0.5,Fyoff=0.5,Fxoff=0.5)
# Define variational problem
ds=Measure('ds')[boundary_parts]
T = TrialFunction(V)
v = TestFunction(V)
dT = TrialFunction(V)
k=Constant(0.6)
Cp=Constant(200)
u=uflow
Q=Constant(100)
T_ext=Constant(270)
#uvec=interpolate(u,VectorFunctionSpace(cellMesh, "Lagrange", 2))
F=inner(k*nabla_grad(T),nabla_grad(v))*dx\
-inner(Cp*dot(grad(T),u),v)*dx(2)\
+Q*T*v*ds(1,subdomain_data = boundary_parts)\
-Q*T_ext*v*ds(1,subdomain_data = boundary_parts)
T_=Function(V)
F=action(F,T_)
J=derivative(F,T_,dT)
problem = NonlinearVariationalProblem(F,T_,J=J,bcs=bcs) #,bcs=BC3d
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-8
prm['newton_solver']['maximum_iterations'] = 10000
prm['newton_solver']['relaxation_parameter'] = 1.0
set_log_level(PROGRESS)
solver.solve()
# Save solution in VTK format
file = File("convectionTest.pvd")
file << T_
