Harish Ragul — CFD & Computational Engineering

The starting point of my M.E. thesis was a deceptively simple question: can you solve a PDE for heat conduction in a material with spatially varying conductivity through a web browser? The answer turned out to be yes — but getting there required understanding the finite element method from the ground up.

The Strong Form

We want to find the temperature field $T(\mathbf{x})$ in a domain $\Omega$ that satisfies:

\nabla \cdot (k(\mathbf{x})\,\nabla T) = 0 \quad \text{in } \Omega

with boundary conditions:

T = T_D \quad \text{on } \Gamma_D \qquad (\text{Dirichlet})

k\,\frac{\partial T}{\partial n} = q_n \quad \text{on } \Gamma_N \qquad (\text{Neumann})

k\,\frac{\partial T}{\partial n} = h(T_\infty - T) \quad \text{on } \Gamma_R \qquad (\text{Robin})

The conductivity $k(\mathbf{x})$ can vary spatially — this is what makes the problem interesting and practically relevant (composite materials, non-uniform heating).

Deriving the Weak Form

The finite element method does not solve the PDE directly in this strong form. Instead, we multiply both sides by a test function $v \in V_0$ (which vanishes on $\Gamma_D$ ) and integrate over the domain:

\int_\Omega \nabla \cdot (k\,\nabla T)\,v\,d\Omega = 0

Applying the divergence theorem (Green's first identity):

\int_\Omega k\,\nabla T \cdot \nabla v\,d\Omega = \int_{\partial\Omega} k\,\frac{\partial T}{\partial n}\,v\,d\Gamma

Splitting the boundary integral and substituting the Neumann and Robin conditions:

\int_\Omega k\,\nabla T \cdot \nabla v\,d\Omega = \int_{\Gamma_N} q_n\,v\,d\Gamma + \int_{\Gamma_R} h(T_\infty - T)\,v\,d\Gamma

This is the weak form: find $T \in V$ such that for all $v \in V_0$ :

a(T, v) = L(v)

where:

a(T,v) = \int_\Omega k\,\nabla T \cdot \nabla v\,d\Omega + \int_{\Gamma_R} h\,T\,v\,d\Gamma

L(v) = \int_{\Gamma_N} q_n\,v\,d\Gamma + \int_{\Gamma_R} h\,T_\infty\,v\,d\Gamma

Implementing in FEniCS

The elegance of FEniCS is that the weak form maps almost directly to Python code. Here is a simplified version:

from fenics import *

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, 'P', 1)

# Boundary conditions
bc = DirichletBC(V, Constant(0.0), left_boundary)

# Spatially varying conductivity
k = Expression('1 + x[0]*x[0]', degree=2)

T = TrialFunction(V)
v = TestFunction(V)

a = k * dot(grad(T), grad(v)) * dx
L = Constant(0.0) * v * dx

T_sol = Function(V)
solve(a == L, T_sol, bc)

The dx integration measure automatically handles Gaussian quadrature over the mesh elements. The ds measure handles boundary integrals for the Neumann and Robin terms.

Key Insights from the Implementation

Mesh refinement matters. At first I used a coarse 8×8 mesh — the solution looked smooth but was wrong near the boundaries. Refining to 32×32 and checking the $L^2$ error against an analytic solution revealed the convergence rate matched the theoretical $\mathcal{O}(h^2)$ for linear elements.

Weak BCs vs strong BCs. Neumann and Robin conditions enter naturally through the boundary integral in the weak form — you do not apply them explicitly. Only Dirichlet conditions are imposed strongly (the DirichletBC object modifies the stiffness matrix rows).

Function spaces define accuracy. Using 'P', 1 (linear Lagrange elements) is sufficient for smooth conductivity fields. For sharp gradients (e.g., a crack or interface), 'P', 2 (quadratic) or adaptive refinement would be needed.

The Web Application Layer

The solver is wrapped in a Django view that accepts material parameters and boundary condition values via a JSON POST request, runs the FEniCS solve, and returns the temperature field data as a Matplotlib contour plot encoded in base64. The front-end displays it inline — no page reload required.

This architecture demonstrates a principle I think is underappreciated: numerical solvers can be first-class web services. The same approach could serve interactive CFD tools for students or remote engineering consultations.

What Comes Next