2D heat conduction equation with physical coefficients

variational formulation
1

Hi, beginner on Fenics. I want to solve a basic steady state heat conduction in a 2D plate problem but using PHYSICAL UNITS constants for materials properties (and not constants in non dimensional units as I often see). I want ultimately to parametrize the code to test different materials with different thermal properties.

I started from the code below (that I cleaned and did as clear as possible):
A plate with a flux Psi on left (Neumann) and a temperature theta_0 on right (Dirichlet) for Boundaries.

Originally the code was executed with a non dimensional diffusivity Kappaadim =1 and non dimensional Flux Psiadim (=1). But I started to put physical parameters (conductivity, rho, cp etc…). I have 2 questions:

1-Can I change unity values of Kappaadim and Psiadim to physical values (as I started to do): is the variational formulation still correct ? (the very small physical diffusivity I want to use seems to produce problem to Fenics).

2-How Can I change the Neuman Boudary condition on right by a Robin boundary condition (that would be defined by a Heat Transfer Coefficient H = 1000 W/m2/K and a Text = 20 Celcius for example, if I understood correctly what is a Robin boundary).

Thanks a lot

Nicolas

#Plate with Boundary Conditions (BC): psi (Flux) on left side, theta_0 (temperature) on right side
#C1 Load modules
from dolfin import *
from future import print_function
from fenics import *
import matplotlib.pyplot as plt

#C2 Constants of the problem and BC values
rho = 7800 # kg/m3 (mass per unit volum)
cp = 435 # J/kg/K (heat capacity)
lamda = 50 # W/m/K (thermal conductivity)
#kappaadim = lamda/rho/cp # m2/sec (thermal diffusivity)
#kappaadim = 0.00000025 # m2/sec in physical units
kappaadim = 1 #non dimensional diffusivity
#psiadim= 100 #W/m2, Flux
psiadim = 1 #non dimensional flux
theta_0 = 10. # Celcius, temperature
tol = 1E-14

#C3 Geometry of the problem
X=1
Y=1

#C4 Mesh
nx2=20
ny2=20

mesh2 = RectangleMesh(Point(0.0, 0.0), Point(X, Y), nx2, ny2, diagonal=“left”)

plt.figure(1)
plot(mesh2, title=“Mesh of Plate”)

#C5 Space definition to search the solution
Hh2 = FunctionSpace(mesh2, ‘P’, 1)

#C6 Boundary Conditions

#C6.1 Right side: Dirichlet BC --> Temperature imposed

def RightEdge2(x, on_boundary):
return on_boundary and (abs(x[0] - 1) < tol)
reD2 = DirichletBC(Hh2, Constant(theta_0), RightEdge2)

#C6.2 Left side: Neumann BC --> Flux imposed

class LeftEdge2(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and (abs(x[0] - 0) < tol)
leN2 = LeftEdge2()
boundaries2 = MeshFunction(“size_t”, mesh2, mesh2.topology().dim()-1, 0)
leN2.mark(boundaries2, 1)

ds2 = Measure(“ds”, domain=mesh2, subdomain_data=boundaries2)

C7 Variationnal Formulation

u2 = TrialFunction(Hh2)
v2 = TestFunction(Hh2)
l2 = psiadimv2ds2(1)
a2 = kappaadim*inner(grad(u2), grad(v2))*dx

C8 Find the solution (SOLVER)

u2 = Function(Hh2)
solve(a2 == l2, u2, reD2)

C8 Plot the solution

plt.figure(2)
plu=plot(u2)
plt.title(‘Numerical Solution u2’)
plt.xlabel(‘x’)
plt.ylabel(‘y’)
plt.colorbar(plu)

2

Please format your code using 3x` encapsulation, to make sure indentation and other formatting is correct.

For Robin condition, see for instance the dolfin tutorial: Solving PDEs in Python - <br> The FEniCS Tutorial Volume I

1 Like