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 = psiadim*v2*ds2(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)