Optymalizacja topologii z użyciem metody SIMP

Optymalizacja topologii jest najbardziej powszechnym rodzajem optymalizacji struktury. Jest stosowana w początkowej fazie projektu do przewidywania optymalnego rozkładu materiału w danej początkowej przestrzeni projektu struktury i uwzględnia specyfikację funkcjonalną oraz ograniczenia produkcyjne.

Najbardziej popularną metodą matematycznej optymalizacji topologii jest metoda litego materiału izotropowego z penalizacją (SIMP). Metodę SIMP jako pierwsi zaproponowali Bendsoe i Kikuchi (1988) oraz Rozvany i Zhou (1992). Metoda SIMP przewiduje optymalny rozkład materiału w danej przestrzeni projektu dla określonych przypadków obciążenia, warunków brzegowych, ograniczeń produkcyjnych i wymagań wydajnościowych.

Według Bendsoe'a (1989): „optymalizacja kształtu w najbardziej ogólnym ujęciu powinna polegać na ustaleniu dla każdego punktu w przestrzeni, czy jest w nim materiał, czy go nie ma”. Tradycyjnym podejściem do optymalizacji topologii jest dyskretyzacja domeny w sieć elementów skończonych, zwanych izotropowymi mikrostrukturami litymi. Każdy element jest wypełniony materiałem w obszarach, które wymagają materiału, albo opróżniony z materiału w obszarach, w których można usunąć materiał (reprezentujących puste przestrzenie). Rozkład gęstości materiału w obrębie domeny projektu ρ jest dyskretny, a każdy element ma przypisaną wartość binarną:
  • ρ(e) = 1 w przypadku, gdy materiał jest wymagany (czarny)
  • ρ(e) = 0 w przypadku, gdy materiał jest usuwany (biały)

Na przykładowej ilustracji przedstawiono zoptymalizowany układ materiału w obciążonej belce. Elementy lite o gęstości ρ(e) = 1 są czarne, natomiast elementy puste o ρ(e) = 0 są usuwane.



Wprowadzenie funkcji ciągłego względnego rozkładu gęstości pozwala uniknąć binarnego, nieciągłego charakteru problemu. Dla każdego elementu przypisana gęstość względna może się zmieniać pomiędzy wartością minimalną ρmin i wartością 1, co umożliwia przypisanie elementom gęstości pośrednich (elementy porowate):

ρmin to minimalna dopuszczalna wartość gęstości względnej dla elementów pustych, które są większe od zera. Ta wartość gęstości zapewnia stabilność numeryczną w analizie metodą elementów skończonych.

Ponieważ gęstość względna materiału może zmieniać się w sposób ciągły, moduł Younga materiału w każdym elemencie może również zmieniać się w sposób ciągły. Dla każdego elementu e zależność pomiędzy współczynnikiem gęstości względnej materiału ρe a modułem Younga sprężystości przypisanego modelu materiału izotropowego E0 określa prawo potęgowe:

Współczynnik kary p zmniejsza wpływ elementów o gęstościach pośrednich (elementy szare) na sztywność całkowitą. Współczynnik kary ukierunkowuje rozwiązanie optymalizacyjne na elementy lite czarne (ρe = 1) lub puste białe (ρe = ρmin). Z doświadczeń numerycznych wynika, że odpowiednia jest wartość współczynnika kary p = 3.

Zmniejszenie modułu sprężystości materiału elementu powoduje zmniejszenie sztywności elementu. Zgodnie z metodą SIMP sztywność globalna jest modulowana według wzoru:

gdzie to macierz sztywności elementu, ρmin to minimalna gęstość względna, ρe to gęstość względna elementu, p to współczynnik kary, a N to liczba elementów w domenie projektu.

Na przykład dla elementu o przypisanej gęstości względnej ρe = 0,5, współczynniku kary = 3 i ρmin = 0,001 globalna macierz sztywności jest skalowana według współczynnika (0,001 + (1-0,001)* 0,5 ^3) = 0,12587.

Funkcja celu: Maksymalizacja sztywności

Celem optymalizacji jest często zmaksymalizowanie ogólnej sztywności konstrukcji lub zminimalizowanie jej podatności po usunięciu określonej ilości masy.

Podatność jest miarą ogólnej elastyczności lub miękkości konstrukcji i jest odwrotnością sztywności. Globalna podatność jest równa sumie energii sprężystości lub naprężenia. Minimalizacja globalnej podatności C jest równoznaczna z maksymalizacją globalnej sztywności. Algorytm optymalizacji ma na celu znalezienie — za pomocą iteracyjnego procesu — takich gęstości elementów (będących zmiennymi projektowymi optymalizacji), które minimalizują globalną podatność konstrukcji.



[ue] jest węzłowym wektorem przemieszczenia elementu e, [Ke] jest sztywnością elementu e, a {ρ} wektora zawiera względne gęstości elementów ρe.

W trakcie każdej iteracji optymalizacji muszą być zachowane powiązania masy obiektu docelowego, równowaga globalna siła-sztywność oraz wymagane ograniczenia funkcjonalne:

ve jest objętością elementu, a Mtarget jest docelową masą optymalizacji.


[K{ρ}] jest globalną macierzą sztywności modulowaną wektorem gęstości względnych, {u} jest wektorem przesunięcia, a {F} jest wektorem siły zewnętrznej.


Powyższy wzór zawiera ograniczenia reakcji projektu, takie jak ograniczenia w zakresie naprężeń, przesunięć, częstotliwości drgań własnych itd.

Analiza czułości

Podczas każdej iteracji algorytm optymalizacji przeprowadza analizę czułości, aby ocenić wpływ zmiany gęstości materiału na funkcję celu i zmaksymalizować sztywność.

W ujęciu matematycznym analiza czułości jest wyrażana jako pochodna funkcji celu w odniesieniu do gęstości materiału:



Podczas analizy czułości elementy ważone niskimi współczynnikami gęstości materiału ostatecznie tracą swoje znaczenie strukturalne i są eliminowane podczas kolejnych iteracji.

Kiedy dla każdego elementu czułość obliczana jest niezależnie, bez uwzględniania połączenia między elementami, skutkiem może być nieciągłość materiału i brak połączenia objętości z główną geometrią. To tzw. efekt szachownicy. W celu ograniczenia efektu szachownicy schemat filtrowania stosuje promień wpływu elementu i uśrednia czułość każdego elementu wewnątrz jego obszaru wpływu.

Iteracje optymalizacji są kontynuowane do momentu, gdy zmiany funkcji celu staną się zbieżne i iteracje osiągną kryteria zbieżności.