Topologieoptimierung mit der SIMP-Methode

Diese Seite fasst zusammen, was das Werkzeug berechnet. Bezeichnungen und Gleichungsnummern stimmen mit dem Skript überein (Kap. 5 und 8); die Nummern rechts neben den Gleichungen sind die des Skripts. Skalare sind kursiv, diskrete Vektoren und Matrizen aufrecht fett ($\vu$, $\vy$, $\mK$).

1Aufgabe

Gesucht ist die Verteilung von Material im Entwurfsraum $\Omega_\mathrm{D}$, die unter den gegebenen Lasten und Lagern möglichst steif ist, also die kleinste Nachgiebigkeit hat, und dabei nur das Volumen $V_0$ verwendet:

Die Nachgiebigkeit $J$ ist die Arbeit der äußeren Lasten. Kleine Nachgiebigkeit bedeutet kleine Verschiebungen in Lastrichtung, also eine steife Struktur. Die spannungsbasierte Variante folgt in Abschnitt 8.

2Diskretisierung und gezeichnete Felder

Der Entwurfsraum wird in $n_\mathrm{e} = n_x \times n_y$ gleiche, quadratische, bilineare Elemente der Kantenlänge $h = 1$ und der Dicke $1$ zerlegt (ebener Spannungszustand). Jedes Element $e$ hat eine Dichte $u_e$, $0 < u_\mathrm{min} \le u_e \le 1$; das sind die Entwurfsvariablen. Der Vektor $\vy^{(e)}$ enthält die acht Knotenverschiebungen von Element $e$. Die gezeichneten Felder bedeuten:

• Dichte $u$: der Startentwurf $\vu_0$. Die Iterationen schreiben den aktuellen Entwurf $\vu_i$ in dieses Feld zurück.
• Lager: markierte Knoten (blaue Punkte) sind in $x$ und $y$ fest. Es sind mindestens zwei Knoten nötig, damit keine Starrkörperbewegung möglich ist.
• Last $x$, Last $y$: Kraft je Fläche $\bb = (\bar{b}_x, \bar{b}_y)$ im Element. Jeder der vier Knoten erhält die konsistente Knotenlast $\bb \, A_e / 4$ mit $A_e = h^2 = 1$; positives $\bar{b}_y$ zeigt nach oben. Belastete Elemente sind passiv fest ($u_e = 1$, Abschn. 5.1 des Skripts): sie werden nicht optimiert, sonst könnte die Optimierung das Material unter der Last entfernen.
• Ausschluss (schraffiert): diese Elemente gehören nicht zum Entwurfsraum. Anders als passive leere Elemente ($u_e = u_\mathrm{min}$) werden sie ganz aus dem Modell entfernt: keine Steifigkeit, kein Volumen, kein Filter, und Knoten, die nur an ausgeschlossene Elemente grenzen, werden aus dem Gleichungssystem genommen. $V_\mathrm{D}$ bezieht sich auf die übrigen Elemente. Ein belastetes Element wird nie ausgeschlossen.

3SIMP-Interpolation

Der Elastizitätsmodul jedes Elements wird mit seiner Dichte skaliert,

und die Elementsteifigkeitsmatrix ist $\mK^{(e)} = E(u_e) \, \tilde{\mK}^{(e)}$, wobei $\tilde{\mK}^{(e)}$ die Elementsteifigkeitsmatrix für $E = 1$ ist. Aus den $\mK^{(e)}$ wird $\mK(\vu)$ assembliert.

4Sensitivitäten

Das adjungierte Problem der Nachgiebigkeit hat die Lösung $\vyh = -2 \, \vy$ (Abschn. 4.5 des Skripts); es muss also kein zweites Gleichungssystem gelöst werden. Mit $\partial \mK^{(e)} / \partial u_e = (E_0 - E_\mathrm{min}) \, p \, u_e^{p-1} \, \tilde{\mK}^{(e)}$ folgen die Sensitivitäten

Da ${\vy^{(e)}}^\T \tilde{\mK}^{(e)} \vy^{(e)} \ge 0$ proportional zur Formänderungsenergie des Elements ist, sind alle Sensitivitäten negativ oder null: mehr Material macht die Struktur immer steifer.

5Filter

Ohne Glättung entstehen Schachbrettmuster, und das Ergebnis hängt vom Netz ab. Der Filter mittelt über die Nachbarschaft $\mathcal{N}_i$ aller Elemente $e$, deren Mittelpunkt näher als $R$ am Element $i$ liegt:

• Sensitivitätsfilter (99-Zeilen-Code): die Sensitivitäten werden mit dem zusätzlichen Gewicht $u_e$ geglättet, $g_i \leftarrow \sum_e w(r_{ie}) \, u_e \, g_e \big/ \big( u_i \sum_e w(r_{ie}) \big)$; im Nenner steht im Werkzeug $\max(u_\mathrm{min}, u_i)$. Das ist eine Heuristik: die geglätteten $g_e$ sind keine exakten Ableitungen mehr.
• Dichtefilter (88-Zeilen-Code): Steifigkeit und Volumen werden mit den gefilterten („physikalischen“) Dichten berechnet, und die Ableitungen folgen mit der Kettenregel:

Der Radius $R$ (in Elementlängen) bestimmt die kleinste Stabdicke. Ausgeschlossene Elemente nehmen am Filter nicht teil.

6Update und Bisektion

Das Volumen wird als Volumenanteil gemessen,

und $V_0$ ist der vorgegebene Anteil am Volumen $V_\mathrm{D}$ des Entwurfsraums (im Werkzeug $V_0/V_\mathrm{D}$). Die Dichten werden multiplikativ aktualisiert,

mit dem Dämpfungsexponenten $\beta$, der Schrittweite $|(\vu_{i+1})_e - (\vu_i)_e| \le m$ und der Schranke $u_\mathrm{min} \le (\vu_{i+1})_e \le 1$. Der Lagrange-Multiplikator $\lambda > 0$ der Volumennebenbedingung wird durch Bisektion bestimmt: $V(\lambda)$ fällt monoton, das Intervall $[\lambda_0^-, \lambda_0^+] = [10^{-9}, 10^{9}]$ wird halbiert, bis $(\lambda^+ - \lambda^-) / (\lambda^+ + \lambda^-) < 10^{-4}$ ist.

7Algorithmus

„1 Iteration“ führt die Schritte 2–9 einmal aus, „Optimieren“ wiederholt sie bis zum Abbruch. Wie im 99-Zeilen-Code bricht das Werkzeug ab, wenn sich keine Dichte mehr um mehr als $\mathrm{tol}$ ändert, $\max_e |(\vu_{i+1})_e - (\vu_i)_e| < \mathrm{tol}$, oder wenn die maximale Iterationszahl erreicht ist. Weicht das Volumen der gezeichneten Startdichte von $V_0$ ab, wird es im ersten Update (im Rahmen von $m$) angepasst.

8Spannungsbasierte Optimierung

Statt der Steifigkeit wird die Festigkeit betrachtet: „Minimiere das Volumen und halte die Spannungen unter der zulässigen Spannung.“

Statt der Sensitivität wird der Aktualisierungsfaktor

verwendet, wobei $\sigma_e$ die Vergleichsspannung im Element $e$ ist. Weil die Spannung eines grauen Elements keine physikalische Spannung ist, wird eine „relaxierte“ Spannung $\sigma_e = u_e^q \, \sigma^0_{\mathrm{v},e}$ verwendet. Dabei ist $\sigma^0_{\mathrm{v},e}$ die Von-Mises-Spannung im Elementmittelpunkt, berechnet mit den Dehnungen des Elements und dem Elastizitätsmodul $E_0$ des Vollmaterials, und $q \le p$ (mit Dichtefilter wird $\tilde{u}_e$ verwendet). Die Dichten werden mit $(\vu_{i+1})_e = (\vu_i)_e \, g_e^\beta$ (Schrittweite $m$) aktualisiert: Material wird hinzugefügt, wo die Spannung über der zulässigen liegt ($g_e > 1$), und entfernt, wo sie darunter liegt ($g_e < 1$).

Wie im Skript werden die Faktoren wie die Sensitivitäten geglättet, also mit dem zusätzlichen Gewicht $u_e$: $g_i \leftarrow \sum_e w(r_{ie}) \, u_e \, g_e \big/ \big( \max(u_\mathrm{min}, u_i) \sum_e w(r_{ie}) \big)$. Sonst verkleinern die spannungsfreien leeren Elemente neben einem Stab dessen geglätteten Faktor, und die Stäbe bleiben dünner als nötig.

Die Standardwerte $q = 1$ und $\sigma_\mathrm{allow} = 1$ sind die des Skripts. Für das Kragträger-Beispiel konvergiert das Verfahren damit mit denselben übrigen Parametern wie bei der Nachgiebigkeitsminimierung in rund 50 Iterationen; das Volumen stellt sich auf etwa $V/V_\mathrm{D} = 0{,}31$ ein, und das Material ist mit etwa $0{,}9$ bis $1{,}5 \, \sigma_\mathrm{allow}$ ausgelastet. Ein kleineres $\sigma_\mathrm{allow}$ ergibt mehr Material.

Grenzen. An einspringenden Ecken, an Punktlagern und an Lasteinleitungen sind die Spannungen singulär (Abschn. 8.1 des Skripts). Dort bleibt $\sigma_e > \sigma_\mathrm{allow}$ auch bei $u_e = 1$, und die angezeigte größte Auslastung $\max_e \sigma_e / \sigma_\mathrm{allow}$ wächst bei Netzverfeinerung. Belastete Elemente sind passiv und werden dabei nicht berücksichtigt.

9Parameter