Feat: Core Physics Engine for Fluid-Thermal Topology Optimization (Fixes #86)

[rainerrodrigues]

Description

This PR introduces the core finite element physics architecture required for Fluid-Thermal Topology Optimization, addressing the foundational requirements of Issue #86.

It implements a monolithic, one-way coupled solver tracking Velocity (:u), Pressure (:p), and Temperature (:T) using MixedDofHandler in Ferrite.jl.

Key Mathematical & Architectural Implementations:

1. Mixed Element Formulation (P1-P1-P1): Implements equal-order linear interpolations for velocity, pressure, and temperature.
2. LBB Stabilization: Injects a standard penalty stabilization term (1e-6) into the pressure block to bypass the Ladyzhenskaya-Babuška-Brezzi (LBB) checkerboard instability inherent to P1-P1 Stokes flow.
3. Brinkman Penalization: Implements fluid resistance mapping $\alpha_e = \alpha_{max}(1 - x_e^3)$ directly into the velocity stiffness matrix block (based on Borrvall & Petersson, 2003, Topology optimization of fluids in Stokes flow).
4. Thermal Advection: Extracts the local element velocity vector and constructs the asymmetric advection matrix $\mathbf{K}_{adv}(\mathbf{u}) = \rho c_p \mathbf{N}^T (\mathbf{u} \cdot \nabla \mathbf{N})$ to physically transport heat.
5. Sequential Solver (solve_fluid_thermal): Executes a two-step Picard iteration (solving the pure fluid domain first, then passing the velocity field into the coupled advection-diffusion matrix).

Testing

• test_fluid_thermal.jl added to the test suite.
• Validates DoF allocations (484 DoFs for a 10x10 P1-P1-P1 grid).
• Validates physical boundary constraints (Hot wall at 100°C).
• Literature Validation: The test explicitly verifies the presence of standard Galerkin undershoots (T < 0.0) characteristic of high-Péclet number advection-dominated flows without SUPG stabilization, as mathematically expected (Brooks & Hughes, 1982).
• Validates Brinkman penalty matrix norms.

Next Steps (Out of Scope for this PR)

Wiring this sequential solver into Nonconvex.jl via an objective function (e.g., Power Dissipation) and providing the adjoint sensitivities via ChainRulesCore. I will open a separate tracking issue for this integration.