---

Data Generation

This FEniCSx script simulates fluid flow and heat transfer around a 2D circuit board with heat-generating electronic components.It first solves the transient, incompressible Navier-Stokes equations using a Discontinuous Galerkin (DG) method to find the velocity and pressure fields. Then, it uses this computed velocity field to solve the steady-state convection-diffusion heat transfer equation (with SUPG stabilization) to determine the temperature distribution.

Installation

This project uses conda for environment management. Environment Setup Create a file named environment.yml with the following contents. This specifies all necessary dependencies.

name: circuit_board_ht
channels:
  - conda-forge
dependencies:
  - python=3.10
  - fenics-dolfinx
  - mpi4py
  - petsc4py
  - numpy
  - gmsh
  - adios2

Create and activate the conda environment by running the following commands in your terminal: conda env create -f environment.yml

How to Run the Script

First, make sure you have activated the conda environment first conda activate circuit_board_ht

Then run the simulation python src/data_generation/main.py

The script will:

1. Generate the mesh file using gmsh.
2. Solve the initial Stokes problem.
3. Time-step through the transient Navier-Stokes problem.
4. Solve the steady-state heat transfer problem.
5. Save the results (velocity u.bp, pressure p.bp, and temperature T.bp) in the output/ folder in VTX format, viewable with tools like ParaView.

Formulation and Equations

This simulation solves two main physical problems sequentially.

Part 1: Incompressible Navier-Stokes Equations (DG)

The fluid flow is modeled by the transient, incompressible Navier-Stokes equations, which enforce conservation of momentum and mass.

Strong Form: The equations for the velocity vector $\mathbf{u}$ and pressure $p$ in the fluid domain $\Omega_f$ are:

1. Momentum Equation: $$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{f}$$
2. Continuity Equation (Incompressibility): $$\nabla \cdot \mathbf{u} = 0$$ Where:
   1. $\rho$ is the fluid density.
   2. $\nu$ is the kinematic viscosity.
   3. $\mathbf{f}$ is an external body force (assumed zero).

Weak Form (Discontinuous Galerkin): The script uses a DG method, which is solved on a mixed function space $V \times Q$. The method relies on fluxes and penalty terms across element facets ($\mathcal{F}$) to enforce continuity and stability.

1. Initial Stokes Problem (t=0): The initial Stokes problem bilinear form $a((\mathbf{u}, p), (\mathbf{v}, q))$ combines the standard continuous Galerkin terms with specific DG flux terms. The standard part is: $$a_{std} = \int_{\Omega} \nu \nabla \mathbf{u} : \nabla \mathbf{v} , dx - \int_{\Omega} p (\nabla \cdot \mathbf{v}) , dx - \int_{\Omega} (\nabla \cdot \mathbf{u}) q , dx$$ Boundary Condition Terms: Similar terms are added on the boundary facets $\mathcal{F}_{ext}$ (ufl.ds) to weakly impose the Dirichlet boundary conditions. The crucial DG Terms enforce continuity and penalize jumps with Symmetry Terms (enforce consistency) and Penalty Term (enforce stability for the discontinuous solution - This term penalizes large jumps in velocity across interior facets. $\alpha$ is a penalty parameter and $h$ is the cell diameter):

$$ - \int_{\mathcal{F}_{int}} \nu \left( \mathrm{avg}(\nabla \mathbf{u}) : \text{jump}(\mathbf{v}, \mathbf{n}) + \mathrm{jump}(\mathbf{u}, \mathbf{n}) : \text{avg}(\nabla \mathbf{v}) \right) , dS $$

$$ \ + \int_{\mathcal{F}_{int}} \frac{\alpha \nu}{\text{avg}(h)} \text{jump}(\mathbf{u}, \mathbf{n}) : \text{jump}(\mathbf{v}, \mathbf{n}) , dS$$

2. Transient Navier-Stokes Problem:

   The formulation is modified with three additional terms to account for time and advection:

   (1) Time Derivative (Backward-Euler): This term is split between the LHS and RHS for the implicit time-stepping scheme.

$$\int_{\Omega} \frac{\mathbf{u} - \mathbf{u}_n}{\Delta t} \cdot \mathbf{v} , dx$$

($\mathbf{u}_n$ is the velocity from the previous time step, $\Delta t$ is the time step size). (2.) Advection Term (Upwind Stabilized): The non-linear advection $(\mathbf{u}_n \cdot \nabla) \mathbf{u}$ is linearized using the previous time step velocity $\mathbf{u}_n$. After integration by parts, the contribution is:

$$ \ - \int_{\Omega} \mathbf{u} \cdot (\nabla \cdot (\mathbf{v} \otimes \mathbf{u}_n)) , dx + \text{Upwind Flux Terms}$$

(The Upwind Flux Terms use the operator $\lambda$, which selects the upwind value of $\mathbf{u}$ at internal and external facets to ensure stability for convection).

Steady-State Convection-Diffusion (Heat)

After solving for the flow, the steady-state temperature $T$ is found. Strong Form

The governing equation is the steady-state convection-diffusion equation:

$$\underbrace{\rho c_p (\mathbf{u} \cdot \nabla T)}_{\text{Convection (Fluid)}} - \underbrace{\nabla \cdot (k \nabla T)}_{\text{Diffusion (All)}} = \underbrace{Q}_{\text{Source (Components)}} $$

Where:

• $\rho c_p$ is volumetric heat capacity (only relevant in the fluid domain).
• $k$ is the thermal conductivity, which is piecewise constant (fluid, PCB, components).
• $Q$ is the volumetric heat source (non-zero only in components).

Weak Form (SUPG Stabilized)

The total weak form $F_{total}$ consists of the standard Galerkin form ($F_{std}$) and the Streamline Upwind Petrov-Galerkin stabilization term ($F_{supg}$):

$$F_{total} = F_{std} + F_{supg} = 0$$

Standard Weak Form ($F_{std}$):

$$F_{std} = \int_{\Omega_{fluid}} \rho c_p (\mathbf{u} \cdot \nabla T) T' , dx + \int_{\Omega} k \nabla T \cdot \nabla T' , dx - \int_{\Omega_{comp}} Q T' , dx = 0$$

(Where $T'$ is the test function). This form accounts for convection in the fluid domain and diffusion (conduction) across all domains, balancing them against the heat source.

SUPG Stabilization Term ($F_{supg}$): This term is added only in the fluid domain to prevent non-physical oscillations in convection-dominated flows.

$$F_{supg} = \int_{\Omega_{fluid}} \tau \left( \rho c_p (\mathbf{u} \cdot \nabla T') \right) R(T) , dx$$

Where $R(T)$ is the residual of the strong-form equation, and $\tau$ is the stabilization parameter, calculated based on the Péclet number ($Pe$):

$$\tau = \begin{cases} \frac{h}{2 |\mathbf{u}|} & \text{if } Pe \ge 1 \text{ (convection-dominated)} \ \frac{h^2}{12 \alpha} & \text{if } Pe < 1 \text{ (diffusion-dominated)} \end{cases}$$

($\alpha = k / \rho c_p$ is the thermal diffusivity). This $\tau$ parameter ensures that artificial diffusion is added only where necessary and in the streamline direction.