### Fenics Heat conduction-convection problem

848
views
0
10 months ago by
Hello everyone,

I am new to Fenics and I am trying to solve the simple heat conduction - convection problem of a  cube with a prescribed temperature at left end and convective heat transfer in the right end. Using the fenics tutorial chap 3 I have tried to frame the heat equation pde in variational form.

To verify the pde I have used a simple square plate with a prescribed temperature of 873K at left end and specified a convective zone with htc = 0.01 at the right end of the plate.

The free stream temperature is set at 293K.

All the other constants e.g. thermal conductivity, density, are set at 1.00

I am running the  model under steady state conditions by removing the ‘dt’ from the pde.

On solving the problem for T ( temperature at the right edge) I should get the following :    K*A * dT/dx = htc * A *(T-To)  : T = 873+293 / 2 = 583 K
However fenics provides an answer of 129.9 K

If anyone has any idea where I am going wrong, please let me know. I have pasted the code for reference.
Thank you very much.
### Heat Conduction and Convection problem - simple 2d plate

from __future__ import print_function
from fenics import *

rho = 1.0 # density
k =  1.0  # thermal conductivity coefficient
C_p = 1.0  # specific heat capacity
source = 0
T_left = 873.0
T_amb = 293.0
h_c = 0.01   # convection heat transfer coefficient
flux = 0.0

# Create mesh and define function space
nx = ny = 4
mesh = RectangleMesh(Point(0, 0),Point(100, 100), nx, ny,"right/left")
V = FunctionSpace(mesh, 'P', 2)

# Define boundary conditions

class LeftEdge(SubDomain):
def inside(self, x, on_boundary):
tol = 1e-6
return on_boundary and abs(x[0]) < tol

class RightEdge(SubDomain):
def inside(self, x, on_boundary):
tol = 1e-6
return on_boundary and near(x[0], 100, tol)

left_edge = LeftEdge()
right_edge = RightEdge()
sub_domains = FacetFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
right_edge.mark(sub_domains, 1)
left_edge.mark(sub_domains, 2)

ds = Measure('ds', domain = mesh, subdomain_data = sub_domains)
dx = Measure('dx', domain = mesh, subdomain_data = sub_domains)

# Applying a prescribed temperature at left edge

bc = DirichletBC(V, T_left, left_edge)

# Define initial value
u_n = project(T_amb,V)

# Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
K_o = k / (rho*C_p)

#  Define the pde

#F_tran = (u*v*dx) + (K_o*inner(grad(u), grad(v))*dx) + ((u_n + dt*source)*v*dx) + (dt*h_c*u*v*ds(2)) +(dt*h_c*T_amb*v*ds(2))

F = (u*v*dx) + (K_o*inner(grad(u), grad(v))*dx) - ((u_n + source)*v*dx) + (h_c*u*v*ds(2)) + (h_c*T_amb*v*ds(2))      #   removing the dt term for steady state soln.

a, L = lhs(F), rhs(F)

# Compute solution

u = Function(V)
t = 0
solve(a == L, u, bc)

print(u.vector().array())
print(u.vector().array().min(),u.vector().array().max())

interactive()
​
Community: FEniCS Project
2
I'm not sure what you're trying to solve, but my guess is a stationary heat equation without sources, i.e., the Laplace equation
$\Delta u = 0$
with boundary conditions
$\left. u \right|_{\Gamma_{left}} = 873, \qquad \left. [\nabla u] n \right|_{\Gamma_{right}} = h_c (u-u_{amb}).$

If this is the case, the weak formulation reads
$\langle \Delta u, v \rangle = - \langle \nabla u, \nabla v \rangle + \langle [\nabla u] n, v \rangle_{\Gamma_{right}} = - \langle \nabla u, \nabla v \rangle + \langle h_c (u-u_{amb}), v \rangle_{\Gamma_{right}}.$

In your formulation, you didn't remove dt rather than set it to dt = 1.0. To find the steady state solution, drop the time derivative in the original (strong) formulation.
written 10 months ago by Adam Janecka