I am using FEniCS to perform thermal simulation. And I set one of boundary surfaces is adiabatic. However, I use the temperature solution to compute the heat flux of that boundary surface and it is not zero, which is conflicted with the adiabatic boundary condition. Any suggestions or ideas about this?
The standard implementation of an adiabatic/insulated boundary using Galerkin’s method is to enforce it weakly, through the variational form (i.e., dropping the boundary term after integrating the flux divergence by parts). For a finite mesh, this does not satisfy the boundary condition in an exact pointwise sense (in contrast with strong boundary conditions, built into the trial function space and implemented in FEniCS via DirichletBC). However, the flux should converge to zero as the mesh is refined.
Thanks for your valuable reply and the boundary condition is implemented as what you said. You are right that the heat flux will be approach to zero with finer and finer mesh. I will try the strong boundary condition as what you said.
If you want to apply flux boundary conditions strongly, you would typically use a mixed formulation, as demonstrated with FEniCS for the Poisson equation (i.e., steady heat conduction) here. Of course, in that formulation, you’ll find that your flux sigma is not exactly equal to grad(u) pointwise, and u does not exactly satisfy its Dirichlet boundary condition; these two conditions are enforced weakly in the mixed formulation.