🔢 mathematics

On the Calculation of the Brinkman Penalization Term in Density-Based Topology Optimization of Fluid-Dependent Problems

This paper investigates how the maximum inverse permeability limit for the Brinkman penalization term in density-based topology optimization of fluid-dependent problems depends on mesh size and flow conditions, aiming to replace trial-and-error selection with a rigorous analysis of the Navier-Stokes equations.

Original authors: Mohamed Abdelhamid, Aleksander Czekanski

Published 2026-08-06

📖 7 min read🧠 Deep dive

CC BY 4.0

Original authors: Mohamed Abdelhamid, Aleksander Czekanski

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are a master architect trying to design the perfect water slide, a super-efficient heart valve, or a filter that cleans oil without clogging. You want to find the absolute best shape for these objects, but you don't know what that shape looks like yet. This is where a field called "topology optimization" comes in. Think of it as a digital sculptor that starts with a block of clay and magically eats away the parts that aren't needed, leaving behind the most efficient shape possible.

However, there's a tricky problem when the "clay" is actually fluid (like water or air) moving through a solid. In the real world, water flows through pipes but stops dead when it hits a wall. In a computer simulation, the software doesn't naturally know where the wall is until the design is finished. To solve this, scientists use a clever trick called "Brinkman penalization." Imagine sprinkling a magical, invisible "thickening powder" into the water. In the areas where you want the water to flow, the powder is absent, and the water moves freely. In the areas you want to be solid walls, the powder is poured in at maximum concentration, turning the water into something as thick as molasses or concrete, effectively freezing it in place. This allows the computer to treat the whole space as one continuous fluid and let the math decide where the walls should be.

But here's the catch: how much of this "thickening powder" do you need? If you use too little, the water leaks through your solid walls like a sieve. If you use too much, the computer gets confused and the numbers go haywire, crashing the simulation. For years, engineers have had to guess the right amount by trial and error, often checking their work against different grid sizes or flow speeds. This paper dives deep into that guessing game to find the exact recipe.

The Paper's Mission: Finding the Perfect "Thickening" Recipe

In this study, authors Mohamed Abdelhamid and Aleksander Czekanski set out to stop the guesswork. They wanted to figure out exactly how the amount of "thickening powder" (which they call the maximum inverse permeability limit, or $\alpha_{max}$ ) should change depending on the size of the computer grid (mesh size) and how fast the water is flowing.

To understand their findings, let's look at how they broke it down. They started by looking at the math behind the fluid's movement, specifically the Navier-Stokes equations (the rules that govern how fluids move). They analyzed these rules in two ways: first, the "strong" form (the raw, continuous math), and second, the "discretized" form (how the computer actually chops the fluid into tiny little boxes or elements to solve the problem).

The Big Discovery: It's All About the Grid Size

The authors found a very clear relationship between the size of the computer's grid and how much "thickening powder" is needed. Imagine your computer screen is made of tiny square tiles. If you make those tiles smaller (a finer mesh), the "thickening powder" needs to be much stronger to stop the water from leaking through the gaps.

Through their mathematical analysis and computer simulations, they discovered that the required strength of the penalty ( $\alpha_{max}$ ) is inversely proportional to the square of the mesh size ( $h$ ). In plain English: if you cut your grid size in half, you need to increase the "thickening" power by a factor of four to keep the walls solid. They also found a secondary relationship where the penalty scales with the inverse of the mesh size itself.

They tested this with a specific design problem: a modified beam in a channel. They ran simulations with different grid sizes (from $h = 1/30$ meters down to $h = 1/190$ meters) and different amounts of "thickening powder" (ranging from 0 up to $10^{20}$ ). The results showed a perfect, straight-line relationship on a log-log scale between the grid size and the penalty needed to keep the water speed in the "solid" areas below a tiny threshold (like $10^{-12}$ meters per second). They even created a specific formula (Equation 33 in the paper) that engineers can use to calculate the exact number needed based on their grid size and how "solid" they want the walls to be.

Flow Speed and Viscosity: The Other Ingredients

The paper also looked at how the speed of the water and its "stickiness" (viscosity) affect the recipe.

Viscosity ( $\mu$ ): They found that if the fluid is stickier (higher viscosity), you need a stronger penalty to stop it. Their simulations showed a linear relationship: as viscosity goes up, the required $\alpha_{max}$ goes up.
Flow Speed ( $v_c$ ): Similarly, if the water is rushing faster, you need a stronger penalty to hold it back. The relationship here is also linear.
Fluid Density ( $\rho_f$ ): Interestingly, they found that the density of the fluid (how heavy it is) barely matters at all, as long as the flow isn't moving at extreme speeds. Whether the fluid is light or heavy, the "thickening powder" requirement stays roughly the same.

The Tricky Part: The Size of the Channel

One part of the puzzle didn't fit the simple math they expected. They looked at the characteristic length of the channel ( $L_c$ ), which is basically the width of the entrance. They initially thought the penalty would scale with the square of this length (like the grid size), but the computer simulations said "nope." The relationship was more complex. Instead of a simple square rule, they found that the penalty scales with the length raised to a power of about 0.6. They suspect this is because the "thickening" effect depends on the tiny, microscopic structure of the porous material, which relates to the big channel size in a complicated way.

What They Ruled Out

The authors explicitly ruled out a few common assumptions.

Trial and Error is Overkill: They showed that you don't need to blindly guess and check. You can calculate the number directly if you know your grid size and flow conditions.
Density Doesn't Matter: They proved that for most practical engineering problems, you don't need to worry about the fluid's density when setting this penalty.
No "Accuracy Loss" at High Numbers: Unlike some previous studies that suggested using a huge penalty number would crash the simulation or lose accuracy, the authors found that in their setup, you can crank the number up to $10^{20}$ without the solution falling apart, as long as the grid is fine enough.

How Sure Are They?

The authors are very confident in their findings, but with a specific caveat. They derived these relationships mathematically and then proved them through computer simulations using a specific design problem (the modified beam in a channel). They are not claiming these numbers are universal laws for every single fluid problem in the universe. Instead, they suggest that these proportional relationships are generally true for this type of problem, and that engineers only need a few data points to tweak the specific numbers for their own unique designs.

In short, this paper hands engineers a new tool: a way to stop guessing and start calculating exactly how much "thickening powder" they need to freeze the water in their digital designs, ensuring the final shape is perfect whether they are using a coarse grid for a quick test or a fine grid for the final product.

Technical Summary: On the Calculation of the Brinkman Penalization Term in Density-Based Topology Optimization of Fluid-Dependent Problems

Problem Statement
In density-based topology optimization (TO) of fluid-dependent problems, the Brinkman penalization method is the standard approach for interpolating between fluid and solid domains. This method introduces a volumetric force term, proportional to fluid velocity, to enforce near-zero velocity in solid regions. The magnitude of this force is governed by a maximum inverse permeability limit, denoted as $\alpha_{max}$ .

Currently, the selection of $\alpha_{max}$ relies heavily on trial-and-error analysis. Practitioners typically choose a value high enough to suppress flow in solid regions but low enough to avoid numerical instability. While this approach suffices for single design problems, it presents significant challenges when calibrating interpolation parameters across different mesh sizes or flow conditions. Specifically, parameters calibrated on a coarse mesh often fail to produce consistent behavior on the finer meshes required for final results, particularly in complex multiphysics problems like fluid-structure interaction (TOFSI). The paper addresses the lack of a systematic, analytical method to determine $\alpha_{max}$ based on mesh size and flow conditions.

Methodology
The authors investigate the dependency of $\alpha_{max}$ on mesh size ( $h$ ) and flow conditions (Reynolds number components: characteristic velocity $v_c$ , dynamic viscosity $\mu$ , characteristic length $L_c$ , and fluid density $\rho_f$ ) through a two-pronged approach:

Analytical Derivation: The study begins by analyzing the Navier-Stokes equations in both their strong (PDE) and discretized finite element (FE) forms. By examining the diagonal terms of the resulting system matrices, the authors derive proportionality relations between $\alpha_{max}$ and the governing parameters.
- Mesh Size: The analysis reveals that the original Navier-Stokes terms (viscous and convective) scale with derivatives of shape functions, which are inversely proportional to $h$ and $h^2$ . In contrast, the Brinkman penalization term matrix is independent of these derivatives. Consequently, $\alpha_{max}$ must scale inversely with $h$ and $h^2$ to maintain balance.
- Flow Conditions: Through non-dimensionalization, the authors deduce that $\alpha_{max}$ is dependent on viscosity ( $\mu$ ) and characteristic length ( $L_c$ ), and dependent on characteristic velocity ( $v_c$ ) due to the convective term. They hypothesize independence from fluid density ( $\rho_f$ ) within reasonable Reynolds number ranges.
Numerical Verification and Curve Fitting: To validate the analytical derivations and determine specific proportionality coefficients, the authors conducted a series of numerical experiments using a modified "beam in a channel" benchmark problem.
- Setup: The problem was solved using both a proprietary MATLAB code and COMSOL Multiphysics. A parametric sweep was performed over a wide range of $\alpha_{max}$ values ($0$ to $10^{20}$ ) and varying parameters ( $h$ , $\mu$ , $L_c$ , $v_c$ , $\rho_f$ ).
- Metric: The validity of the penalization was assessed by measuring the maximum velocity in the solid regions ( $max(|v_{solid}|)$ ) and the relative errors in fluid state variables.
- Fitting: Using log-log plots, the authors performed curve fitting on selected data points to derive explicit algebraic expressions for $\alpha_{max}$ as a function of the input parameters for a specific target solid velocity.

Key Contributions and Results
The paper provides the following specific findings and contributions:

Analytical Dependencies: The study confirms that $\alpha_{max}$ is inversely proportional to the square of the mesh size ( $h^2$ ) and the mesh size ( $h$ ), and is dependent on $\mu$ , $v_c$ , and $L_c$ . It is largely independent of $\rho_f$ .
Empirical Formulas: The authors derived specific equations to calculate $\alpha_{max}$ $α_{ma x}$ for a target solid velocity ( $10^q$ $1 0^{q}$ ):
- Mesh Size: $\alpha_{max} \propto \frac{1}{h^2} + \frac{1}{h}$ . The fitted relation (Eq. 33) shows a strong linear log-log correlation.
- Viscosity ( $\mu$ ): A linear relationship was found: $\alpha_{max} \propto \mu$ (Eq. 34).
- Characteristic Velocity ( $v_c$ ): A linear relationship was found: $\alpha_{max} \propto v_c$ (Eq. 36).
- Characteristic Length ( $L_c$ ): Contrary to the initial hypothesis of a simple $1/L_c^2$ scaling derived from non-dimensionalization, curve fitting revealed a more complex power-law relationship: $\alpha_{max} \propto L_c^{-0.6073}$ (Eq. 35). The authors suggest this discrepancy arises because the relevant length scale for the penalization is a microscale characteristic of the porous microstructure, which relates to the macroscale $L_c$ in a non-trivial way.
- Density ( $\rho_f$ ): Numerical results confirmed that $\alpha_{max}$ is effectively independent of $\rho_f$ for variations within an order of magnitude, as viscous forces dominate inertia in the penalized solid regions.
Accuracy: The derived formulas demonstrated high accuracy, with maximum errors generally below 3.4% when compared to numerical data points, provided the target solid velocity is kept within the linear regime of the log-log relationship (typically $max(|v_{solid}|) \leq 10^{-2}$ ).

Significance and Claims
The authors claim that this work provides a systematic, physics-based method for calculating the Brinkman penalization limit, moving beyond the traditional trial-and-error approach. The primary significance lies in the ability to calibrate $\alpha_{max}$ such that different mesh sizes and flow conditions produce consistent physical behavior. This is particularly valuable for:

Multi-scale Optimization: Allowing the use of coarse meshes for initial parameter calibration (interpolation and projection) before refining to computationally expensive fine meshes without re-tuning parameters.
Robustness in TOFSI: Addressing the strong nonlinearity and nonconvexity in fluid-structure interaction problems by ensuring the penalization term behaves predictably across varying discretizations.

The paper modestly notes that while the proportionality relations are generally true for the physics involved, the specific coefficients are problem-dependent. However, the authors demonstrate that only a small number of data points are required to calculate these coefficients for new problems, provided the analytical dependencies are known. The work does not propose new applications but rather offers a foundational tool to improve the reliability and efficiency of existing density-based topology optimization workflows for fluid problems.

More like this