Difficulty computing smooth local Nusselt number on inner circular boundary in FEniCSx

General
1

I’m trying to compute the local Nusselt number (dimensionless temperature gradient) along an inner circular boundary of a 2D natural convection problem using FEniCSx (DOLFINx + UFL).

Specifically, I need

image

where T is temperature, n is the unit normal pointing outward from the fluid region, and D is the cylinder diameter.

My setup:

  • Mesh from Gmsh with physical groups: outer walls and a circular inner boundary.

  • facet_markers are available and correct.

  • Temperature field T is a second-order Function(Q).

  • I can compute average Nusselt number using

n = ufl.FacetNormal(domain)
ds_cyl = ufl.Measure(“ds”, domain=domain, subdomain_data=facet_markers, subdomain_id=2)
Nu_avg = fem.assemble_scalar(fem.form(-dot(grad(T), n) * ds_cyl))

However, I’d like to get local values of -∂T/∂n per boundary facet (or vertex) for post-processing (e.g., plotting Nu vs θ around the cylinder).

I’ve tried looping over facets and assembling the flux using a facet-specific measure like:

ds_tag = ufl.Measure(“ds”, domain=domain, subdomain_data=my_tags, subdomain_id=facet_id)
flux = fem.assemble_scalar(fem.form(inner(grad(T), n) * ds_tag))

but I’m unsure if this gives a consistent local gradient or if I should use a projection approach (e.g., project -dot(grad(T), n) onto a function space defined on the boundary).

Questions:

  1. What’s the correct or recommended way in FEniCSx to compute local ∂T/∂n on a curved internal boundary?

  2. How can one obtain facet-wise or vertex-wise values that are numerically smooth (e.g., suitable for plotting Nu vs θ)?

(I can provide a minimal example if needed — but essentially, it’s a square cavity with an inner circular boundary, steady-state temperature field from a diffusion problem.)

2

Note that the n in question, the inwards pointing normal, from a mesh perspective is not continuous.
You are using the facet normal, which is a piecewise constant vector per facet.
Similarly, the gradient grad(T) is not continuous at element boundaries.

Therefore, I would suggest evaluating the field at the midpoint of each facet, where the field is continuous and well defined.

You can for instance do that with the functions in:

specifically; move_to_facet_quadrature. You can choose the number of points by controlling the degree of the quadrature space on the interface.