🧠 Building a Topology Optimization Solver from Scratch: A Journey in FEM, SIMP, and Structural Efficiency

By Anmol Shrestha , Piyush K.C., and Nepal Keshav

---

✨ Introduction

Topology optimization is revolutionizing how we design mechanical structures. By mimicking the efficiency of natural forms—like bone structures or tree branches—this method helps us place material only where it is structurally necessary.

But most engineers interact with topology optimization through black-box software tools like ANSYS or Abaqus. What if we could break free from the black box and build our own solver from the ground up?

That's exactly what we did. This blog details our journey of developing a 2D topology optimization solver using the SIMP method, implemented entirely in C++. It walks through the mathematics, programming decisions, FEM formulation, and performance comparison with ANSYS.

---

🎯 The Goal

We aimed to develop a topology optimization solver that:

• Uses the SIMP (Solid Isotropic Material with Penalization) method

• Solves 2D plane stress problems using Finite Element Method (FEM)

• Optimizes structural stiffness (compliance minimization)

• Compares favorably with industry-standard tools like ANSYS

---

🧱 Understanding the Foundations

Before optimization, we had to build the underlying FEM infrastructure.

📐 Plane Stress Assumption

For thin 2D structures, we use the plane stress condition, which simplifies the stress-strain relationship:

σ = D * ε

Where:

• σ is the stress vector [σx, σy, τxy]

• ε is the strain vector [εx, εy, γxy]

• D is the material stiffness matrix for isotropic material under plane stress:

ini
D = (E / (1 - ν²)) * 
    [  1    ν    0
       ν    1    0
       0    0  (1 - ν)/2 ]

E is the Young’s modulus, and ν is Poisson’s ratio.

---

🧮 Element Stiffness Using CST Elements

We used Constant Strain Triangle (CST) elements due to their simplicity and efficiency. The stiffness matrix for each element is given by:

Ke = t * A * Bᵀ * D * B

Where:

• t is the thickness

• A is the area of the triangle

• B is the strain-displacement matrix

• D is the material stiffness matrix

These element matrices are assembled into a global stiffness matrix (K) using standard FEM procedures.

---

🧠 SIMP Topology Optimization: The Core Algorithm

Topology optimization modifies the structure by varying material density (denoted as Xe) for each element. The objective is to minimize compliance (maximize stiffness) while using a limited amount of material.

🎯 Objective Function

The compliance (C) is given by:

C = Uᵀ * K * U

Where:

• U is the displacement vector

• K is the global stiffness matrix

Our goal is:

Minimize C, subject to: Σ Xe ≤ V₀ × volume fraction

Each element stiffness is penalized to avoid intermediate densities using:

Ke = Xeᵖ * K₀

Where:

• p is the penalization factor (typically 3)

• K₀ is the stiffness of a fully solid element

---

🔁 Optimization Loop

The solver iterates through these steps:

1. Initialize mesh and densities

2. Assemble global stiffness matrix

3. Apply boundary conditions and solve K * U = F

4. Compute compliance and sensitivities

5. Update densities using optimality criteria and a Lagrange multiplier

6. Repeat until convergence

The Lagrange multiplier ensures the total volume constraint is satisfied through a bisection method.

---

🧪 Validation and Results

To validate the solver, we used a benchmark test: a cantilever beam loaded at the tip. We compared the optimized structure with an I-beam of the same mass simulated in ANSYS.

Metric	Optimized Structure	I-Beam (ANSYS)
Deflection (mm)	0.41191	0.50641
Volume / Mass	Same	Same
Compliance (Uᵀ * K * U)	Lower	Higher

This means the optimized structure was ~18.6% stiffer than the standard I-beam, using the same amount of material.

---

📤 From Data to Design: CAD Reconstruction

After exporting the final material distribution as a CSV file, we:

• Applied image thresholding to identify solid vs void regions

• Cleaned the geometry using morphological operations

• Rebuilt the structure in SolidWorks using extrusion

This allows the result of the topology optimization to be converted into a fabricable 3D model.

---

⚠️ Challenges and Fixes

We faced several technical challenges:

• Checkerboarding effects due to the lack of filtering

• Singular stiffness matrices, fixed using regularization

• Convergence instability, resolved via volume constraint bisection

Despite these issues, the solver remained stable and produced meaningful output.

---

✅ Conclusion

This project demonstrates the feasibility of creating a custom topology optimization solver that matches the performance of commercial tools for basic structural cases. More importantly, it provides full control and visibility over every step in the simulation and optimization process.

Our solver:

• Successfully minimized compliance within volume constraints

• Showed superior stiffness compared to traditional I-beams

• Was validated through FEM simulations in ANSYS

• Produced exportable data for 3D CAD reconstruction

---

🚀 What’s Next?

We plan to:

• Integrate density filtering to avoid checkerboarding

• Extend the solver to 3D topology optimization

• Add stress-based failure constraints

• Build a simple GUI for educational users

• Explore additive manufacturing constraints

---

🙏 Acknowledgements

Special thanks to:

• Assoc. Prof. Dr. Krishna P. Shrestha (Supervisor)

• Dr. Surendra Sujakhu (Co-supervisor)

• Department of Mechanical Engineering, Kathmandu University