### Fenics Heat conduction-convection problem

423

views

0

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.

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

Please login to add an answer/comment or follow this question.

\[ \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.That solved the problem