Fenics Heat conduction-convection problem

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)
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)


Community: FEniCS Project
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  
Thank you very much Adam.
That solved the problem
written 10 months ago by Subhro  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »