Topology Optimization and homogenization problem
This post documents my journey learning topology optimization, homogenization, and inverse homogenization problems. As a computer science student with zero background in graphics or material physics, I'm focusing on understanding the implementation rather than diving deep into mathematical proofs. I'm writing this primarily as my own study notes, so it won't cover every detail - definitely check out the original papers for the complete picture. Here are the three key articles I'm working through:
- Sigmund, O. (2001). A 99 line topology optimization code written in Matlab. Structural and Multidisciplinary Optimization, 21(2), 120-127. https://doi.org/10.1007/s001580050176
- Xia, L., Breitkopf, P. Design of materials using topology optimization and energy-based homogenization approach in Matlab. Struct Multidisc Optim 52, 1229–1241 (2015). https://doi.org/10.1007/s00158-015-1294-0
- Zhang, D., Zhai, X., Liu, L. et al. An optimized, easy-to-use, open-source GPU solver for large-scale inverse homogenization problems. Struct Multidisc Optim 66, 207 (2023). https://doi.org/10.1007/s00158-023-03657-y
The first article covers the fundamentals of topology optimization from an engineering perspective and includes a complete 99-line MATLAB implementation. It's the perfect entry point if you're new to the field - the code is surprisingly compact yet fully functional.
The second article extends topology optimization into material design using homogenization theory. Originally, my advisor suggested starting here, but I quickly hit a wall with unfamiliar concepts. After working through Paper #1 first, everything clicked much better. Learn from my mistake and read these in order.
The third article tackles the same material design problem as Paper #2 but moves all computations to GPU for dramatic speedup. Once you understand the concepts from Paper #2, this reads smoothly since it's mainly about clever implementation optimizations rather than new theory.
Topology Optimization Introduction
First, what is topology optimization? Imagine you're designing a chair with limited material. How do you create a structure that's both durable and uses the least amount of material possible? That's essentially what topology optimization solves.
In the example above, I gave the computer two key inputs: where forces will be applied (people sitting) and boundary conditions (chair legs fixed to the ground - otherwise the chair would just slide away due to Newton's laws). I also specified a volume fraction (how much total material I'm willing to use). The algorithm then calculates the optimal structure that satisfies all these constraints. The resulting chair design is pretty striking - I'd definitely consider buying one if it hit the market.
There are multitude of approaches for the solving of topology optimization problems. In the original paper Bendsøe and Kikuchi (1988) a so-called microstructure or homogenization based approach was used. An alternative approach to topology optimization is the so-called “power-law approach” or SIMP approach (Solid Isotropic Material with Penalization) (Bendsøe 1989; Zhou and Rozvany 1991; Mlejnek 1992).
$$ \begin{align} \min_{\mathbf{x}} : \quad & c(\mathbf{x}) = \mathbf{U}^T \mathbf{K} \mathbf{U} = \sum_{e=1}^{N} (x_e)^p \mathbf{u}_e^T \mathbf{k}_0 \mathbf{u}_e \\ \text{subject to} : \quad & \frac{V(\mathbf{x})}{V_0} = f \\ & \mathbf{K}\mathbf{U} = \mathbf{F} \\ & 0 < x_{\min} \leq \mathbf{x} \leq 1 \end{align} $$
The foundation of topology optimization can be summarized to be the function above. Where $U$ and $F$ are the global displacement and force vectors. $K$ is the global stiffness matrix, $u_e, k_0$ are the element displacement vector and stiffness matrix. $p$ is the penalization power (typically $p=3$). $x_min$ is a fixed small number and used to ensures every element has a tiny bit of "fake material" so the math works, but it's so small it doesn't affect the final design. $U^T K U$ is the mathematical way to compute total strain energy or compliance. For the optimization problem above, our goal is that with the condition volume fraction is fixed: $\frac{V(x)}{V_0} =f$, structure is in static equilibrium: $KU=F$ This optimization asks: "Given a fixed amount of material and applied loads, where should I place material to create the stiffest possible structure?",or in mathematical representation: find our designed variable $x$ which contains each element density information ( $x_{min}\lt x \le 1$, where $1$ represent solid material, $x_{min}$ means essentially void material, and grey material if in the middle)。
The optimization problem can be solved using serval different approaches such as Optimality Criteria (OC) methods, Sequential Linear Programming (SLP) methods or the Method of Moving Asymptotes (MMA by Svanberg 1987) and others. In this article, for simplicity, they use OC-method. Following Bendsøe (1995) a heuristic updating scheme for the design variables can be formulated as:
$$ x_e^{\text{new}} = \begin{cases} \max(x_{\min}, x_e - m) & \text{if } x_e B_e^{\eta} \leq \max(x_{\min}, x_e - m) \\ x_e B_e^{\eta} & \text{if } \max(x_{\min}, x_e - m) < x_e B_e^{\eta} < \min(1, x_e + m) \\ \min(1, x_e + m) & \text{if } \min(1, x_e + m) \leq x_e B_e^{\eta} \end{cases} $$
where $m$ (move) is a positive move limit, $\eta$ ($=1/2$) is a numerical damping coefficient and $B_e$ can be found from the optimality conditions (the mathematical condition that must be satisfied at the optimal solution) as:
$$ B_e = \frac{-\frac{\partial c}{\partial x_e}}{\lambda \frac{\partial V}{\partial x_e}} $$
Where $\lambda$ is the Lagrange multiplier - it balances the trade-off between compliance and volume, and it can be found by a bi-sectioning algorithm.
The sensitivity of the objective function is found as:
$$ \frac{\partial c}{\partial x_e} = -p(x_e)^{p-1} \mathbf{u}_e^T \mathbf{k}_0 \mathbf{u}_e $$
In order to ensure existence of solutions to the topology optimization problem, some sort of restriction on the resulting design must be introduced. Here we use a filtering technique (Sigmund 1994, 1997). It must be emphasized that this filter has not yet been proven to ensure existence of solutions, but numerous applications by the author have proven the the filter produces mesh-independent designs in practice.
The mesh-independency filter works by modifying the element sensitivities as follows:
$$ \frac{\overline{\partial c}}{\partial x_e} = \frac{1}{x_e \sum_{f=1}^{N} \tilde{H}_f} \sum_{f=1}^{N} \tilde{H}_f x_f \frac{\partial c}{\partial x_f} $$
$$ \tilde{H}_f = r_{\min} - \text{dist}(e,f) $$
$$ \{f \in N | \text{dist}(e,f) \leq r_{\min}\}, \quad e = 1,\ldots,N $$
where the operator $\text{dist}(e, f)$ is defined as the distance between center of element $e$ and center of element $f$. Please looks really careful to this formula, its goal is to apply a filter to smooths the sensitivities by averaging them over neighboring elements. let's analyze how such filtering works step by step.
- For each element e, look at all neighboring elements $f$ within radius $r_{\min} $
- Weight each neighbor's sensitivity by $\tilde{H}_f $ (closer neighbors have more influence)
- Also weight by density $x_f $ (solid elements have more influence than voids)
- Normalize by the total weights
Now, after understanding all the prior formula to do topology optimization, we can then write code. However, I will not dive into the code now because we read this article only for having prior knowledge for later articles.
Design Material using Topology Optimization and Homogenization approach
Now, let's come to the second article. This guy extends topology optimization from designing structures (like the chair example) to designing the internal microstructure of materials themselves. It also provide an educational MATLAB code that can automatically design materials with specific mechanical properties by optimizing how solid and void regions are arranged at the microscale. The key innovation behind this paper is that it use energy-based homogenization (simpler than asymptotic homogenization) to predict how a complex microstructure behaves as a bulk material. Basically, You can input desired material properties (like "I want maximum stiffness" or "I want negative Poisson's ratio") and the algorithm designs the optimal microstructure pattern!
Below are some definitions for terms in material science in cases you don't know.
Homogenization: A mathematical technique to find the "effective" or "average" properties of complex materials. Think of a sponge - instead of analyzing every tiny pore, homogenization gives you the overall behavior as if it were a single uniform material.
Stress $\sigma$: Internal force per unit area within a material when external forces are applied. Like pressure, but inside the material. Measured in force/area (Pa, MPa).
Strain $\varepsilon $: How much a material deforms relative to its original size. If a 10cm bar stretches to 11cm, the strain is 10%. Dimensionless (just a ratio).
Poisson's Ratio: When you stretch something in one direction, it usually gets thinner in the perpendicular directions. Poisson's ratio measures this effect. Most materials have positive ratios (get thinner when stretched). "Negative Poisson's ratio" materials get thicker when stretched - very unusual and useful.
Shear Modulus: Resistance to shearing forces (like sliding one layer over another). Think of how hard it is to twist or slide parts of a material relative to each other.
Young's Modulus: Specific measure of stiffness in tension/compression. The slope of the stress-strain curve in the elastic region. Higher values = stiffer material.
Bulk Modulus: Resistance to uniform compression in all directions. High bulk modulus = hard to compress (like water). Related to how materials behave under pressure.
Elasticity Matrix $E_{ijkl}$: measures how much stress component $\sigma_{ij}$ results from strain component $\varepsilon_{kl}$, that is also the reason why stiff matrix is a order 4 tensor. The dimension for $E_{ijkl}$ is $R^{6\times 6}$ for 3D case because stress and strain matrix only have 6 degree of freedom.
Homogenization
Within the scope of linear elasticity, the equivalent constitutive behavior of periodically patterned microstructures can be evaluated using the homogenization method.
$$ u^{\epsilon}(x) = u_0(x,y) + \epsilon u_1(x,y) + \epsilon^2 u_2(x,y) + \ldots, \quad y = x/\epsilon $$
In this formula, $\epsilon$ measures the aspect ratio between macro and micro scale, and $y=x/ \epsilon$. This formula illustrate that the displacement field has components at different scales - some vary slowly (macroscale) and some vary rapidly (microscale patterns).
For the Homogenized Stiffness Tensor, it can be calculate by averaging the integral over the the base cell $Y$ as
$$ E^H_{ijkl} = \frac{1}{|Y|} \int_Y E_{ijpq}(\varepsilon^{0(kl)}_{pq} - \varepsilon^{*(kl)}_{pq})dY $$
where:
- $E^H_{ijkl}$: The homogenized (effective) stiffness tensor we want to find
- $|Y|$: Volume of the unit cell
- $E_{ijpq}$: Local material stiffness (varies within the unit cell)
- $ε^{0(kl)}_{pq}$: Unit test strain fields (prescribed strains we apply)
- $ε^{*(kl)}_{pq}$: Fluctuation strain fields (how the material responds)
This formula averaging the local material response over the entire unit cell to get effective properties. In physical meaning, this formula describe how does the material internally deform when you apply a unit strain to the unit cell.
More specifically, for each component $E^H_{ijkl}$, this formula:
Run test case ($kl$)
- Apply unit strain in direction $kl$
- Solve for resulting displacement field
At each point in the unit cell
- Calculate: $ε^{0(kl)}*{pq} - ε^{*(kl)}*{pq}$
- Multiply by local stiffness: $E_{ijpq}$
- This gives local stress contribution
Integrate over unit cell and normalize
- Sum all local contributions
- Divide by unit cell volume
The relation between stress and strain also satisfied: stress is the derivative of energy density with respect to strain. Hence, with the intension to favor effective existing algorithms used in topology optimization, the formula above can be rewritten in an equivalent form in terms of element mutual energies
$$ E^H_{ijkl} = \frac{1}{|Y|} \int_Y E_{pqrs}\varepsilon^{A(ij)}_{pq} \varepsilon^{A(kl)}_{rs} dY $$
The finite element implementation for the energy form above is:
$$ E^H_{ijkl} = \frac{1}{|Y|} \sum_{e=1}^N (u^{A(ij)}_e)^T k_e u^{A(kl)}_e $$
Where $k_e$ is the element stiffness matrix. $(u^{A(ij)}_e)^T k_e u^{A(kl)}_e$ calculate the energy contribution from each element. This formula basically done the following steps:
- Apply unit test strains to the unit cell boundaries
- Solve the FEM system to get displacements $u^{A(ij)}_e$
- Compute element energies ($u^T k u$) for each element
- Sum over all elements to get the homogenized property
In 2D case, we can rewrite the above formula to the following expanded form:
$$ \begin{bmatrix} E^H_{11} & E^H_{12} & E^H_{13} \\ E^H_{21} & E^H_{22} & E^H_{23} \\ E^H_{31} & E^H_{32} & E^H_{33} \end{bmatrix} = \begin{bmatrix} Q_{11} & Q_{12} & Q_{13} \\ Q_{21} & Q_{22} & Q_{23} \\ Q_{31} & Q_{32} & Q_{33} \end{bmatrix} $$
where $Q_{ij}=\frac{1}{|Y|}\sum^{N}_{e=1}q_e^{ij}$ are the sums of element mutual energies $q_e^{ij}=(u^{A(i)}_e)^T k_e u^{A(j)}_e$
Periodic Boundary Condition (PBC)
Similar to the chair example, we need to have a boundary condition to ensure the unit cell deforms with the desired macroscopic strain. in microscope homogenization problem, our boundary condition is set to be:
$$ u_i^{k+}-u_i^{k-}=\varepsilon_{ij}^0(y_j^{k+}-y_j^{k-}=\varepsilon_{ij}^0\Delta y_j^k) $$
which means that The displacement difference between opposite boundaries equals the macroscopic strain times the cell dimension.
Inverse Homogenization Problem
I just directly jump to the article since the article 2 is painful to read. Article 3 use one paper long to introduce how to apply inverse homogenization problem clearly. Even though it does not touch the deeper principle like why those formula make sense etc., but come on, we are doing engineering, knowing how to use is enough for us. For the following section, I will introduce how to apply the inverse homogenization (IHP) step by step.
What IHP try to solve: Given some expected structure characteristic, like Poisson ratio, and constrains include volume, volume fraction; generate a structure satisfied such characteristic and constrain
above is a structure generated by inverse homogenization, with set Poisson ratio=-0.6644 and volume fraction=0.2. The inverse homogenization problem is an optimization problem as following:
[TODO]
Multigrid V-cycle Solver
To solve a system of equation:
$$ Ax=b $$
of course we can use Gauss-Seidel approach, which iteratively update updates each variable in $x$ by one using the most recent values. However, Gauss-Seidel has a major problem: it's very slow for large systems, especially those arising from discretized PDEs. For a 1000×1000 grid, Gauss-Seidel might need 10,000+ iterations to converge! In addition, it eliminates local errors quickly but global errors very slowly. For example:
Exact: ———————————————————— (smooth curve)
Current: ∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿∿ (wiggly approximation)For local error (high frequency one), the wiggles neighboring points disagree Global (Low-Frequency) Errors: The overall level is wrong - all points shifted up/down. Multigrid V-cycle Solver solves the global error problem by using a hierarchy of coarser grids.
https://www.youtube.com/watch?v=4UoJpazu_fU&ab_channel=FluidMechanics101
