Wieloskalowe modelowanie tkanki kostnej, Egzaminy z Geometria

Pobierz Wieloskalowe modelowanie tkanki kostnej i więcej Egzaminy w PDF z Geometria tylko na Docsity! POLITECHNIKA ŚLĄSKA WYDZIAŁ MECHANICZNY TECHNOLOGICZNY Instytut Mechaniki i Inżynierii Obliczeniowej mgr inż. Przemysław Makowski Wieloskalowe modelowanie tkanki kostnej Praca doktorska Promotor: dr hab. inż. Wacław Kuś prof. Pol. Śl. – Gliwice 2015 – 2 Spis treści 1 Wstęp .................................................................................................................................. 4 1.1 Przegląd zawartości pracy ........................................................................................... 5 1.2 Teza i cel pracy ............................................................................................................ 6 2 Modelowanie tkanki kostnej ............................................................................................... 8 2.1 Układ kostny człowieka ............................................................................................... 8 2.2 Budowa tkanki kostnej .............................................................................................. 10 2.3 Modele materiału dla tkanki kostnej ......................................................................... 14 2.4 Badania doświadczalne tkanki kostnej ...................................................................... 17 2.5 Obrazowanie medyczne tkanki kostnej ..................................................................... 21 2.6 Tworzenie modeli numerycznych tkanki kostnej na podstawie danych tomograficznych ................................................................................................................... 27 2.7 Metoda elementów skończonych w modelowaniu tkanki kostnej ............................ 30 2.8 Biorusztowania i przebudowa tkanki kostnej ............................................................ 37 3 Modelowanie wieloskalowe w mechanice ........................................................................ 42 3.1 Modelowanie wieloskalowe struktur kostnych ......................................................... 45 3.2 Metoda homogenizacji numerycznej dla zagadnień liniowo sprężystych ................. 48 3.3 Metoda homogenizacji numerycznej dla zagadnień nieliniowych ............................ 53 3.4 Algorytmy ewolucyjne w zastosowaniach wieloskalowych ..................................... 55 4 Algorytm homogenizacji numerycznej dla wieloskalowych modeli tkanki kostnej ........ 57 4.1 Weryfikacja algorytmu homogenizacji numerycznej dla zagadnień liniowych ........ 57 4.2 Weryfikacja algorytmu homogenizacji numerycznej dla zagadnień nieliniowych ... 62 5 Zastosowanie metod wieloskalowych w modelowaniu tkanki kostnej i implantów ........ 68 5.1 Identyfikacja i optymalizacja dla wybranych modeli wieloskalowych tkanek i implantów ............................................................................................................................. 70 5.2 Identyfikacja osobniczych parametrów materiałowych beleczek kostnych na podstawie danych eksperymentalnych i numerycznych ....................................................... 71 5.3 Homogenizacja struktury kości beleczkowej – obliczenie zastępczych parametrów materiałowych tkanki kostnej ............................................................................................... 77 5 prób doświadczalnych oraz modelowania numerycznego pozwala na przeprowadzenie symulacji niemożliwych do realizacji przy użyciu metod eksperymentalnych. Powszechnie stosowane modele numeryczne kości zakładają jednorodność oraz izotropowość parametrów mechanicznych tkanki co często stanowi zbyt daleko idące uproszczenie struktury kości. Tkanka kostna, zwłaszcza kość beleczkowa jest materiałem niejednorodnym i anizotropowym o mikrostrukturze porowatej. Kość jest strukturą hierarchiczną. Parametry tkanki, takie jak np. anizotropia, dla skali całej kości (np. udowej) wynikają z budowy struktury kości w skali mikro. Metody modelowania wieloskalowego, stosowane dla materiałów inżynierskich o budowie niejednorodnej takich jak materiały kompozytowe wzmacniane włóknami, mogą zostać zastosowane w celu modelowania numerycznego hierarchicznej struktury kości oraz implantów tkanki kostnej. Podejście wieloskalowe umożliwia uwzględnienie mikrostruktury kości przy zachowaniu akceptowalnych czasów obliczeń oraz pozwala na obliczanie rozkładów analizowanych pól fizycznych zarówno na poziomach mikro oraz makro. 1.1 Przegląd zawartości pracy Praca składa się z sześciu rozdziałów. W rozdziale pierwszym zamieszczono wstęp pracy oraz sformułowano tezę i cel pracy. W rozdziale drugim opisano strukturę obiektu badań tzn. kości beleczkowej, modelowanej z użyciem modeli liniowych oraz z nieliniowościami fizycznymi. Scharakteryzowano typowe badania doświadczalne prowadzone na próbkach kostnych oraz opisano podstawy obrazowania medycznego tkanki kostnej z użyciem rentgenowskiej tomografii komputerowej. Przedstawiono metodykę tworzenia modeli numerycznych tkanek na podstawie danych tomograficznych. Omówiono metodę elementów skończonych (MES) zastosowaną w pracy do realizacji analiz numerycznych. Opisano biorusztowania kości beleczkowej oraz proces przebudowy kości. Rozdział trzeci zawiera informacje dotyczące metod modelowania wieloskalowego w mechanice oraz biomechanice. Opisano metodę homogenizacji numerycznej dla zagadnień liniowych oraz z nieliniowościami fizycznymi. Opisano zasadę działania oraz zastosowania algorytmów ewolucyjnych dla zagadnień identyfikacji i optymalizacji w układach wieloskalowych. W rozdziale czwartym przedstawiono zadanie weryfikacji zaimplementowanych algorytmów homogenizacji numerycznej z użyciem symulacji próby ściskania dla wieloskalowych modeli 6 uproszczonej struktury tkanki kostnej. Określono błąd dla symulacji numerycznych zrealizowanych z użyciem modeli wieloskalowych oraz porównano czasy obliczeń z modelami dokładnymi. Rozdział piąty przedstawia przykłady zastosowania modeli wieloskalowych tkanek i implantów. Zidentyfikowano parametry materiałowe kości w skali mikro oraz obliczono anizotropowe zastępcze parametry materiałowe kości dla skali makro. Rozwiązano problem optymalizacji topologicznej struktury osobniczego biorusztowania kości beleczkowej z użyciem algorytmu ewolucyjnego i trójskalowego modelu uwzględniającego metodę wytwarzania implantu. Dla optymalnej struktury implantu przeprowadzono wieloskalową analizę układu bliższego końca kości udowej z zaimplantowanym biorusztowaniem kości beleczkowej. Rozdział szósty zawiera wnioski końcowe i podsumowanie zrealizowanych badań. 1.2 Teza i cel pracy Celem pracy jest opracowanie modeli wieloskalowych umożliwiających obliczanie anizotropowych parametrów materiałowych tkanki kostnej z uwzględnieniem hierarchicznej budowy kości oraz projektowanie osobniczych biorusztowań kości z użyciem procedur optymalizacji struktury implantu. Sformułowano następującą tezę pracy: Modelowanie wieloskalowe umożliwia analizę tkanki kostnej uwzględniając jej budowę hierarchiczną oraz własności anizotropowe, jak również może zostać wykorzystane do projektowania biorusztowań kości beleczkowej. 7 Realizacja celu pracy oraz udowodnienie powyższej tezy wymaga rozwiązania następujących zadań cząstkowych:  opracowanie, implementacja komputerowa oraz weryfikacja metody modelowania wieloskalowego tkanek oraz implantów,  budowa modelu numerycznego mikrostruktury kości beleczkowej na podstawie danych mikrotomograficznych,  identyfikacja parametrów materiałowych tkanki kostnej w skali mikro na podstawie danych numerycznych i eksperymentalnych,  opracowanie algorytmu zadawania periodycznych warunków brzegowych dla nieperiodycznych objętościowych siatek MES modeli mikrostruktury tkanki kostnej,  obliczenie anizotropowych parametrów zastępczych tkanki kostnej na podstawie modelu numerycznego mikrostruktury,  weryfikacja poprawności ekstrakcji próbki kostnej z zastosowaniem transformacji układu współrzędnych tensora sprężystości,  opracowanie trójskalowego modelu biorusztowania kości beleczkowej, uwzględniającego metodę wytwarzania implantu,  optymalizacja spersonalizowanego biorusztowania kości beleczkowej,  przeprowadzenie analizy wieloskalowej bliższego końca kości udowej z zaimplantowanym biorusztowaniem. 10 oraz fosforu), funkcję krwiotwórczą (produkcja komórek krwi przez czerwony szpik kostny) oraz funkcję przekazywania informacji poprzez system nerwowy [120]. Rys. 2.2 Struktura kości stawu ramiennego człowieka z wyraźnym podziałem na istotę zbitą i gąbczastą [175] 2.2 Budowa tkanki kostnej Tkanka kostna będąca głównym budulcem kości jest strukturą hierarchiczną o różnorodnej budowie w zależności od skali obserwacji (Rys. 2.3, Rys. 2.4). Rys. 2.3 Hierarchiczna struktura kości beleczkowej [50] 11 Rys. 2.4 Obrazy mikroskopowe struktury kości beleczkowej w skali a) nano b) sub-mikro c) mikro d) meso [50] Na poziomie skali nano, tkanka kostna zbudowana jest białek (głównie włókna kolagenu typu I), składnika mineralnego (głównie kryształy hydroksyapatytu, HA), komórek kostnych, naczyń krwionośnych oraz wody. Włókna kolagenu charakteryzują się ciągłą budową z wieloma nałożeniami co umożliwia przenoszenie obciążeń pomiędzy poszczególnymi włóknami. Kryształy składnika mineralnego w kształcie igieł lub płytek ułożone są wzdłuż oraz pomiędzy włóknami kolagenu. W odniesieniu do kompozytowych materiałów inżynierskich, woda oraz składniki organiczne (osteoid) stanowią matrycę, natomiast składniki mineralne (kryształy) stanowią inkluzje. Wśród komórek kostnych można wyróżnić osteoblasty, osteocyty oraz osteoklasty. Osteoblasty to tzw. komórki kościotwórcze biorące udział w procesie tworzenia oraz gojenia kości, wytwarzając składniki organiczne istoty międzykomórkowej kości. Powstają z komórek macierzystych narastających w szpiku kostnym. Osteoblasty podlegają przemianie w osteocyty – dojrzałe komórki kostne otoczone zmineralizowaną istotą międzykomórkową kości. Osteocyty umiejscowione są w tzw. jamkach kostnych gdzie poprzez kanaliki kostne i wypustki komórek, biorą udział w procesie wymiany substancji odżywczych i metabolitów. Główną funkcją osteoklastów czyli komórek kościogubnych jest niszczenie tkanki kostnej czyli resorpcja. Są to komórki wielojądrowe, znacznie większe od pozostałych komórek kostnych. Osteoblasty i osteoklasty biorą udział w naturalnych procesach tworzenia i niszczenia kości (remodelingu) znajdujących się w stanie równowagi dynamicznej. Zaburzenie tej równowagi może prowadzić np. do nadmiernej resorpcji kości czyli osteoporozy. Kolejny poziomem (skala sub-mikro) w hierarchicznej budowie tkanki kostnej są blaszki kostne (lamellae) zbudowane z równoległych włókien kolagenowych oraz składników mineralnych. W zależności od przestrzennej organizacji i wzajemnego ułożenia blaszek kostnych wyróżniamy kość gąbczastą oraz kość zbitą (skala mikro). Kość zbita zbudowana jest z gęsto ułożonych blaszek kostnych co przekłada się na wysoką wytrzymałość i parametry mechaniczne tkanki. W przypadku kości zbitej, blaszki układają się koncentrycznie wokół kanału Haversa tworząc osteon – podstawową jednostkę struktury 12 kości zbitej w formie współśrodkowych cylindrów znajdujących się jeden w drugim. Ze względu na gęste ułożenie blaszek kostnych, kość zbita jest strukturą stosunkowo jednorodną w porównaniu z kością gąbczastą dla której blaszki kostne uporządkowane są w odmienny sposób tworząc beleczki kostne (trabecula) (Rys. 2.5) [120]. Rys. 2.5 Struktura beleczek kostnych [147] W jamkach kostnych wewnątrz beleczek znajdują się osteocyty połączone z innymi komórkami za pomocą wypustek cytoplazmatycznych biegnących w kanalikach kostnych. Powierzchnia beleczek pokryta jest śródkostną w której znajdują się osteoblasty i mniej liczne osteoklasty [120] [140]. Przestrzenie międzybeleczkowe wypełnione są szpikiem kostnym. Beleczki kostne mogą charakteryzować się geometrią beleczkową, płytkową lub pośrednią w zależności od miejsca występowania w kości (Rys. 2.6). Rys. 2.6 Porowata struktura kości beleczkowej z wyraźną granicą kości zbitej [174] 15 Dla materiału anizotropowego liczba niezależnych stałych sprężystych wynosi 21. Liczba ta ulega zmniejszeniu w zależności od ilości płaszczyzn symetrii sprężystej przechodzących przez dany punkt ciała. Dla materiału anizotropowego nie można wyróżnić żadnej płaszczyzny symetrii sprężystej. W przypadku występowania w materiale trzech wzajemnie prostopadłych płaszczyzn symetrii sprężystej mamy do czynienia z przypadkiem ciała otrotropowego. Krawędzie przecięcia płaszczyzn symetrii wyznaczają tzw. główne kierunki anizotropii (główne kierunki ortotropii, osie materiałowe) x1, x2, x3. Liczba niezależnych elementów macierzy sprężystości redukuje się w tym przypadku do dziewięciu (2.3):                                                            11 11 12 13 11 22 22 23 22 33 33 33 12 44 12 23 55 23 31 66 31 ζ c c c 0 0 0 ε ζ c c 0 0 0 ε ζ c 0 0 0 ε = ζ sym. c 0 0 ε ζ c 0 ε ζ c ε (2.3) Macierzą odwrotną do macierzy sprężystości C jest macierz podatności S (2.4). Wartości elementów macierzy podatności S są powiązane z wartościami stałych inżynierskich materiału ortotropowego: -1= =                                      3121 1 2 3 3212 1 2 3 13 23 1 2 3 12 23 31 -ν-ν1 0 0 0 E E E -ν-ν 1 0 0 0 E E E -ν -ν 1 0 0 0 E E E 1 0 0 0 0 0 G 1 0 0 0 0 0 G 1 0 0 0 0 0 G S C (2.4) gdzie Ei – moduły sprężystości wzdłużnej (moduły Younga) wyznaczone względem głównych kierunków ortotropii, Gij – moduły sprężystości poprzecznej (moduły Kirchhoffa) wyznaczone względem głównych kierunków ortotropii, νij – współczynniki Poissona wyznaczone względem głównych kierunków ortotropii. 16 W przypadku występowania w materiale równoległych do siebie płaszczyzn na których własności ośrodka są identyczne, natomiast ulegają zmianie w kierunku normalnym do tych płaszczyzn, wówczas mamy do czynienia z przypadkiem ciała transwersalnie izotropowego. Płaszczyzny w których własności materiału są identyczne nazywany płaszczyznami izotropii. Można przyjąć, że materiał transwersalnie izotropowy jest szczególnym przypadkiem materiału ortotropowego dla którego własności sprężyste w pewnej płaszczyźnie są sobie równe. Dla materiału transwersalnie izotropowego liczba stałych sprężystych redukuje się do pięciu. W przypadku występowania płaszczyzn izotropii równoległych do płaszczyzny wyznaczonej przez osie x1 oraz x2 zachodzą następujące związki (2.5) w odniesieniu do ciała ortotropowego: E1 = E2 ν12 = ν21 ν31 = ν32 G23 = G31 (2.5) Ciało izotropowe posiada nieskończoną ilość płaszczyzn symetrii sprężystej, tak więc własności mechaniczne materiału izotropowego są identyczne we wszystkich kierunkach. W tym wypadku liczba niezależnych stałych sprężystych redukuję się do dwóch, a macierz sprężystości C związku konstytutywnego (2.2) przyjmuje następującą postać: =                    1 - ν ν ν 0 0 0 1 - ν ν 0 0 0 1 - ν 0 0 0E 0.5 - ν 0 0(1+ ν)(1 - 2ν) sym. 0.5 - ν 0 0.5 - ν C (2.6) gdzie E to moduł Younga, natomiast ν to współczynnik Poissona materiału izotropowego. Wartość modułu Kirchhoffa jest wyznaczana na podstawie zależności E G = 2(1+ ν) [3] [8]. Powszechne stosowany liniowo sprężysty, izotropowy modelu materiału dla tkanki kostnej jest pewnym uproszczeniem. W rzeczywistości kość jest materiałem anizotropowym, nieliniowo lepkosprężystym. Brak jednak wystarczającej ilości porównywalnych wyników badań doświadczalnych, identyfikujących lepkosprężyste parametry kości. Własności lepkosprężyste kości beleczkowej wynikają przede wszystkim z obecności szpiku kostnego w przestrzeniach międzybeleczkowych. Jednakże wpływ własności lepkosprężystych jest pomijalny dla zakresu prędkości odkształcenia odpowiadających kości podczas normalnej aktywności człowieka. Dopiero dla bardzo wysokich prędkości odkształcenia rzędu 10 s -1 17 zaobserwować można wzrost modułu sprężystości tkanki kostnej. Uwzględnienie parametrów lepkosprężyztych wywołanych obecnością szpiku kostnego jest zasadne dla zagadnień badania oddziaływania sił dynamicznych na układ kostny człowieka (np. wypadki komunikacyjne, urazy sportowe). W pozostałych przypadkach szpik kostny i jego oddziaływanie są pomijalne i nie są uwzględniane [73] [21]. Przy pominięciu własności lepkosprężystych, kość może być rozpatrywana jako materiał nieliniowo sprężysty. Klasa materiałów hiposprężystych [160] pozwala na modelowanie ośrodków sprężystych o nieliniowych związkach pomiędzy tensorami naprężenia i odkształcenia. Naprężenie jest w tym wypadku nieliniową funkcją odkształcenia (2.7). Materiał kości beleczkowej może wykazywać nieliniowość w zakresie małych odkształceń [110]. Dla zakresu dużych odkształceń konieczne jest operowanie miarami naprężenia i odkształcenia, takimi jak przykładowo tensor naprężenia Cauchy’ego (ang. true stress tensor) oraz tensor odkształcenia Hencky’ego (ang. true strain tensor, logarithmic strain tensor). Równanie konstytutywne dla klasy materiałów hiposprężystych przyjmuje następującą formę w postaci przyrostowej: ij ijkl kldζ = D dε (2.7) gdzie ijdζ – przyrost naprężenia kldε – przyrost odkształcenia ijklD – tensor stałych materiałowych dla modelu hiposprężystego Tensor stałych materiałowych dla modelu hiposprężystego ijklD jest tensorem czwartego rzędu zależnym od stanu naprężenia. Forma przedstawionego związku konstytutywnego (2.7) wymaga zastosowania implementacji MES w postaci przyrostowej z użyciem uaktualnianego opisu współrzędnych materialnych Lagrange’a (ang. Updated Lagrange additive formulation), np. metody Newtona-Raphsona [121] [113]. 2.4 Badania doświadczalne tkanki kostnej Celem badań doświadczalnych prowadzonych na próbkach kostnych lub całych kościach jest przede wszystkim określanie parametrów materiałowych tkanki. Pomiędzy poszczególnymi ludźmi występują naturalne różnice osobnicze oraz rozbieżności budowy tkanki wynikające z różnicy wieku, płci czy też przebytych stanów chorobowych. Metody 20 Opisane metody są podstawowymi, najpowszechniej stosowanymi badaniami doświadczalnymi przeprowadzanymi w celu określenia parametrów materiałowych kości w skali makro. Metodą eksperymentalną pozwalającą na określenie parametrów kości beleczkowej w skali mikro (pojedyncze beleczki kostne) oraz skalach niższych jest metoda nanoindentacji [135]. Pomiar z zastosowaniem metody nanoindentacji polega na zagłębianiu wgłębnika Berkovitcha lub Vickersa w strukturę badanej próbki kostnej z dokładnym pomiarem wartości siły obciążającej oraz głębokości penetracji. Badanie jest prowadzone do osiągnięcia zadanej wartości maksymalnej siły obciążającej lub maksymalnej wartości penetracji wgłębnika po czym następują fazy utrzymania stałego obciążenia i odciążenia próbki. Na podstawie wykresu siła obciążająca – głębokość penetracji dla fazy odciążania, wyznaczany jest moduł Younga badanej struktury z użyciem zależności teorii kontaktu Hertza oraz teorii sprężystości. Badanie ma charakter punktowy i powinien być przeprowadzony dla co najmniej kilku punktów pomiarowych struktury, ze względu na niejednorodną budowę tkanki kostnej [76]. Parametry materiałowe beleczki kostnej mogą zostać również określone na podstawie próby trójpunktowego zginania pojedynczej beleczki kostnej (Rys. 2.10) oraz symulacji numerycznej MES eksperymentu [63] [90]. Rys. 2.10 Próba trójpunktowego zginania pojedynczej beleczki kostnej [63] Badania doświadczalne służące określaniu parametrów mechanicznych pojedynczej beleczki kostnej wymagają zastosowania specjalistycznego oprzyrządowania oraz zapewnienia bardzo dokładnego pomiaru sił i przemieszczeń z uwagi na niskie zakresy mierzonych wielkości, wynikające ze skali eksperymentu oraz samej próbki. 21 2.5 Obrazowanie medyczne tkanki kostnej Obrazowanie medyczne to grupa metod umożliwiająca utworzenie graficznej reprezentacji wnętrza ciała ludzkiego, czyli struktur znajdujących się pod skórą czy też np. wewnątrz kości. Obrazowanie medyczne stosowane jest przede wszystkim w celach diagnostycznych takich jak wykrywanie torbieli, wad rozwojowych narządów wewnętrznych, obszarów martwicy tkanek czy też zmian nowotworowych. Za początek rozwoju metod obrazowania medycznego przyjąć można rok 1895 w którym niemiecki fizyk Wilhelm Röntgen opublikował wyniki swojej pracy [136] dotyczącej zastosowania elektromagnetycznego promieniowania X (nazwanego później promieniowaniem rentgenowskim), opartej częściowo na wczesnych badaniach Nikola Tesli. Badanie z użyciem aparatu rentgenowskiego wykorzystujące zjawisko zróżnicowania stopnia osłabiania (pochłaniania) promieni X przez różne struktury anatomiczne jest jedną z podstawowych metod obrazowania medycznego. Wynikiem badania rentgenowskiego (RTG) jest dwuwymiarowy obraz prześwietlanych struktur kostnych (Rys. 2.11) znajdujących się między lampą rentgenowską i np. filmem światłoczułym, co skutkuje nakładaniem się na siebie częsci obrazowanych struktur. Rys. 2.11 Zdjęcie rentgenowskie okolicy stawu nadgarstkowego człowieka Metodą badania wykorzystującą w swym działaniu promieniowanie rentgenowskie, jednak umożliwiającą już utworzenie graficznej reprezentacji pojedynczych wybranych przekrojów lub też trójwymiarowej reprezentacji wnętrza ludzkiego ciała jest rentgenowska tomografia komputerowa, nazywana również krócej tomografią komputerową (ang. Computed Tomography ,CT). Pierwszy prototyp tomografu komputerowego został zbudowany w roku 1968 przez Godfreya Hounsfielda w oparciu o podstawy teoretyczne rekonstrukcji obrazów 22 sformułowane przez Allana McLeoda Cormacka. W roku 1971 wykonano pierwsze badanie diagnostyczne z użyciem tomografu komputerowego, natomiast w roku 1973 wprowadzono na rynek pierwszy komercyjny model tomografu EMI CT 1000 [25]. Ogólna zasada działania rentgenowskiego tomografu komputerowego polega na generowaniu wiązki promieniowania X z użyciem lampy rentgenowskiej wykonującej ruch obrotowy wokół skanowanego obszaru (pacjenta). Jednocześnie przez układ detektorów rejestrowane są wartości osłabienia natężenia wiązki promieniowania rentgenowskiego na podstawie której zgodnie z prawem osłabienia promieniowanie rentgenowskiego (2.9) obliczany jest liniowy współczynnik osłabienia μ dla danej struktury. -μx 0I = I e (2.9) gdzie I – natężenie wiązki promieniowania po przejściu przez ośrodek o grubości x [W/sr] I0 – natężenie początkowe wiązki promieniowania padającej na ośrodek [W/sr] μ – liniowy współczynnik osłabienia promieniowania [m -1 ] Wartość liniowego współczynnika osłabienia jest zależna od energii zastosowanego promieniowania oraz gęstości ρ prześwietlanej struktury, tak więc możliwe jest wyznaczenie gęstości struktur anatomicznych na podstawie wyników badania tomograficznego. W trakcie badania tomograficznego, w związku z wykonywaniem ruchu przez lampę rentgenowską rejestrowane są kolejne dwuwymiarowe projekcje prześwietlanej struktury. W obecnie stosowanych, najnowszych tomografach komputerowych czwartej generacji stosowana jest ruchoma lampa rentgenowska generująca wiązkę wachlarzową oraz nieruchomy zamknięty pierścień detektorów składający się nawet z kilku tysięcy elementów. W celu skrócenia czasu badania tomograficznego stosowana jest spiralna tomografia komputerowa polegająca na sprzężeniu ruchu obrotowego lampy rentgenowskiej z liniowym przemieszczeniem stołu tomografu na którym znajduje się pacjent (Rys. 2.12). 25 W przypadku rekonstrukcji trójwymiarowej, prześwietlony obszar zostaje podzielony na sześcienne objętości elementarne, nazywane wokselami (ang. voxel). Podobnie jak ma to miejsce w przypadku rekonstrukcji obrazów dwuwymiarowych, każdemu wokselowi rekonstrukcji trójwymiarowej przypisywana jest wartość gęstości radiologicznej wyrażona w jednostkach skali Hounsfielda [20]. Odmianą tomografii komputerowej jest mikrotomografia komputerowa (μCT). Podobnie jak standardowa metoda rentgenowskiej tomografii komputerowej, metoda mikrotomografii komputerowej jest nieniszczącą metodą obrazowania wewnętrznej struktury prześwietlanych obiektów opierającą się na takich samych zasadach. Główną różnicą między metodami tomografii oraz mikrotomografii komputerowej jest rozdzielczość zrekonstruowanych obrazów prześwietlanych struktur, która w przypadku μCT jest wielkością rzędu mikrometrów (rozdzielczość 1μm lub nawet poniżej 1μm w przypadku najnowszych konstrukcji mikrotomografów). Wyższa rozdzielczość przekłada się na niższą wielkość ogniska (ang. focus size) promieni generowanych przez lampę rentgenowską (<50μm) oraz możliwość skanowania jedynie stosunkowo niewielkich obiektów (wymiary rzędu kilu cm) co ogranicza zastosowania medyczne urządzenia. W przypadku większości mikrotomografów, odwrotnie niż ma to miejsce dla tomografów standardowych, lampa rentgenowska oraz detektory pozostają nieruchome podczas badania, natomiast ruch obrotowy wykonuje skanowana struktura znajdująca się na obrotowym stoliku napędzanym silnikiem krokowym. Pierwszy mikrotomograf został skonstruowany przez Jima Elliotta na początku lat osiemdziesiątych XX wieku [38]. Mikrotomografia znalazła zastosowanie w dziedzinie badań materiałowych czy tez obrazowania wewnętrznej struktury niewielkich zwierząt lub roślin. Na przełomie lat osiemdziesiątych i dziewięćdziesiątych XX wieku, mikrotomografię komputerową zastosowano w celach obrazowania mikrostruktury kości beleczkowej [40]. W przypadku standardowego badania CT, kość beleczkowa jest widoczna na zrekonstruowanych obrazach jako jednorodna struktura o różnej gęstości radiologicznej (zależnej między innymi od porowatości struktury) reprezentowanej w jednostkach skali Hounsfielda. Z kolei w przypadku badania μCT kości beleczkowej, na zrekonstruowanych obrazach widoczna jest już mikrostruktura kości tzn. przestrzenna porowata struktura zbudowana z połączonych beleczek kostnych reprezentowanych przez jaśniejsze obszary (Rys. 2.14). 26 Rys. 2.14 Mikrostruktura kości beleczkowej pozyskana na drodze badania μCT Analiza obrazów tomograficznych mikrostruktury kości beleczkowej pozwala na wyznaczenie parametrów struktury takich jak przykładowo porowatość kości, grubości i długości beleczek kostnych oraz budowę mikrostruktury (struktura beleczkowa, płytkowa lub pośrednia). W połowie lat dziewięćdziesiątych XX wieku wprowadzono na rynek pierwsze komercyjne konstrukcje mikrotomografów [137]. Od tego czasu metoda μCT zyskuje coraz większe zainteresowanie w dziedzinie badań nad kością beleczkową [20]. Oprócz opisanych rentgenowskiej tomografii i mikrotomografii komputerowej, obrazowanie medyczne tkanki kostnej może zostać również zrealizowane z użyciem tomografii i mikrotomografii rezonansu magnetycznego (MRI, μMRI) opartych na zjawisku jądrowego rezonansu magnetycznego. Jednakże ze względu na charakterystykę otrzymywanych obrazów, badanie rezonansem magnetycznym stosowane jest głównie w celach obrazowania tkanek miękkich takich jak mózg, serce czy mięśnie. W celu obrazowania tkanki kostnej zazwyczaj stosowane są dane z rentgenowskiej tomografii i mikrotomografii komputerowej. Dane pozyskane na drodze obrazowania medycznego zapisywane są zgodnie ze standardem cyfrowego zapisu i wymiany danych medycznych DICOM (ang. Digital Imaging and Communications in Medicine). Prace nad standardem rozpoczęto w roku 1983 w celu wprowadzenia standaryzacji i ujednolicenia formatu zapisu danych z obrazowania medycznego i są prowadzone do dnia dzisiejszego. Plik DICOM będący obecnie standardem zapisu danych tomograficznych, oprócz samych obrazów medycznych, zawiera również szereg danych o pacjencie, modelu i parametrach tomografu oraz parametry przeprowadzonego badania. Pliki zawierające informacje medyczne zapisane z użyciem standardu DICOM są odczytywane przez wiele darmowych i komercyjnych oprogramowań służących do przeglądania i analizy danych medycznych, ich segmentacji i przygotowania 27 danych do druku 3D modeli struktur anatomicznych bądź też tworzenia modeli numerycznych poszczególnych struktur anatomicznych, np. kości. 2.6 Tworzenie modeli numerycznych tkanki kostnej na podstawie danych tomograficznych Modele numeryczne tkanek i narządów wewnętrznych mogą zostać budowane na podstawie danych z obrazowania medycznego. W przypadku tkanki kostnej, badaniem pozwalającym najlepiej zobrazować strukturę tkanki jest rentgenowska tomografia komputerowa (CT) oraz rentgenowska mikrotomografia komputerowa (μCT). Modele numeryczne tkanki kostnej dla skali makro (skala całej kości) budowane są na bazie danych CT, natomiast modele struktury kości w skali mikro (porowata mikrostruktura beleczkowa) tworzone są na bazie danych μCT. Bezpośrednia konwersja wokseli na sześcienne elementy skończone odbywa się w sposób automatyczny, jednak otrzymany model struktury może charakteryzować się niską jakością odwzorowania geometrii tkanek i narządów. Ma to szczególne znaczenie w przypadku modeli numerycznych skomplikowanej mikrostruktury kości beleczkowej (Rys. 2.15). Rys. 2.15 Modele numeryczne mikrostruktury kości beleczkowej zbudowane na podstawie konwersji wokseli na elementy sześcienne (góra) oraz modele wygładzone zbudowane z elementów czworościennych (dół) [162] 30 anatomicznej. Do zapisu modelu powierzchniowego stosowany jest format STL (ang. STereoLithography) jednak trójkąty tworzące model mogą być już utożsamiane z powierzchniowymi elementami skończonymi o liniowej funkcji kształtu. Początkowa siatka trójkątów jest modyfikowana (ang. remeshing) w celu redukcji liczby elementów oraz poprawy jakości siatki elementów skończonych. Jakość siatki powierzchniowej wyrażana jest przykładowo znormalizowaną miarą stosunku wysokości elementu trójkątnego do długości jego podstawy (ang. height to base ratio), gdzie wartość 1 powyższego parametru przyporządkowana jest elementowi o najwyższej jakości (trójkąt równoboczny o równych kątach wewnętrznych). Następnie na podstawie zmodyfikowanej siatki powierzchniowej, tworzona jest objętościowa siatka czworościennych elementów skończonych (tetra), a siatka powierzchniowa jest usuwana. Jakość siatki powierzchniowej na podstawie której budowane są czworościenne elementy objętościowe wpływa na początkową jakość elementów objętościowych. Siatka objętościowa może być dalej modyfikowania w celu uzyskania ostatecznej postaci modelu numerycznego. Do tak zbudowanego modelu MES możliwe jest przypisanie niejednorodnego rozkładu parametrów materiałowych (modułów Younga) dla grup elementów skończonych, na podstawie konwersji gęstości radiologicznej danych tomograficznych z użyciem równań konwertujących [134]. Konwersja gęstości radiologicznej nie pozwala jednak na uwzględnienie anizotropii własności mechanicznych tkanek [65]. 2.7 Metoda elementów skończonych w modelowaniu tkanki kostnej Metoda elementów skończonych jest najpowszechniej stosowaną numeryczną metodą obliczania pól przemieszczeń, odkształceń oraz naprężeń dla ośrodków odkształcalnych. W mechanice ciała odkształcalnego poszukuje się pola przemieszczeń spełniającego warunki jednoznaczności. Dla układów o skomplikowanej geometrii i warunkach brzegowych nie możliwe jest wyznaczenie rozwiązania dokładnego na drodze obliczeń analitycznych. Do takich układów niewątpliwie można zaliczyć układ biomechaniczne takie jak struktury kostne czy też mikrostruktury tkanek. Istota metody MES polega na zastąpieniu modelu ciągłego, modelem dyskretnym rozpatrywanego ośrodka, dla którego poszukiwane jest rozwiązanie przybliżone. W przypadku MES, dyskretyzacja to podział obszaru o skomplikowanej geometrii na mniejsze obszary o prostej geometrii nazywane elementami skończonymi. W wierzchołkach elementów skończonych nazywanych węzłami, wyznaczane są wartości poszukiwanych pól fizycznych. W celu spełnienia warunku ciągłości 31 analizowanych pól fizycznych w rozpatrywanym obszarze, wartości węzłowe są interpolowane w obrębie elementów skończonych z użyciem tzw. funkcji kształtu. W wyniku dyskretyzacji początkowy układ równań różniczkowych zastępowany jest układem równań algebraicznych. Dla materiałów modelowanych jako ośrodek liniowo-sprężysty, związek konstytutywny stanowi uogólnione prawo Hooke’a (2.1). W celu obliczenia wartości pól fizycznych w analizowanym układzie sprężystym, poszukiwany jest stan równowagi, któremu odpowiada minimum energii potencjalnej ciała sprężystego. Całkowita energia potencjalna Π (2.13) ciała sprężystego jest sumą energii odkształcenia U (pracy sił wewnętrznych) (2.11) oraz pracy sił zewnętrznych Uzew (2.12), gdzie  T V 1 U = dV 2 σ ε (2.11)   T T T zew i i iV S U = - dV - dS -u f u L u P (2.12) tak więc    T T T T i i iV V S 1 Π = dV - dV - dS - 2 σ ε u f u L u P (2.13) gdzie σ – macierz naprężeń ε – macierz odkształceń u – macierz przemieszczeń f – macierz sił objętościowych L – macierz sił powierzchniowych Pi – macierz sił skupionych przyłożonych w punkcie i V – objętość ciała S – powierzchnia ciała Po przeprowadzeniu dyskretyzacji rozpatrywanego ciała sprężystego, całkowitą energię sprężystą Π można wyrazić jako sumę energii sprężystych elementów skończonych modelu dyskretnego. Dla każdego elementu skończonego definiuje się powiązany z nim lokalny układ współrzędnych oraz funkcje kształtu. W przypadku zastosowania liniowej funkcji interpolacyjnej, mówimy o elementach z liniową funkcją kształtu. Przemieszczenia wewnątrz 32 elementu skończonego (2.14), interpolowane na podstawie wartości węzłowych można wyrazić jako =u Nq (2.14) gdzie u – macierz przemieszczeń danego punktu, interpolowanych wewnątrz elementu skończonego N – macierz funkcji kształtu q – macierz przemieszczeń węzłowych elementu skończonego Dla każdego elementu skończonego wyznaczana jest macierz sztywności elementu zależna od jego geometrii oraz zadanych parametrów materiałowych. Następnie przeprowadzana jest agregacja macierzy sztywności elementów skończonych do globalnej macierzy sztywności całego układu oraz agregacja sił objętościowych, powierzchniowych i skupionych do globalnej macierzy sił . Wówczas równanie całkowitej energii potencjalnej układu sprężystego można zapisać w formie macierzowej (2.15) jako T T1 Π = - 2 Q KQ Q F (2.15) gdzie Q – globalna macierz przemieszczeń węzłowych K – globalna macierz sztywności F – globalna macierz sił węzłowych Warunki brzegowe wprowadzane są z zastosowaniem metody kary bądź metody eliminacji. Następnie obliczane jest ekstremum energii potencjalnej odpowiadające stanowi równowagi. Po zróżniczkowaniu równania (2.15) względem nieznanych przemieszczeń węzłów układu, otrzymujemy układ równań algebraicznych (2.16) w postaci macierzowej =KQ F (2.16) po którego rozwiązaniu otrzymujemy poszukiwane wartości przemieszczeń węzłowych. Odkształcenia obliczane są na podstawie obliczonych przemieszczeń z zastosowaniem zależności (2.17) =ε Bu (2.17) gdzie B jest macierzą geometryczną, zawierającą pochodne funkcji kształtu. Naprężenia obliczane są na podstawie zależności (2.18) = CBu (2.18) gdzie C jest macierzą sprężystości uogólnionego prawa Hooke’a. 35 Ogólna postać równania macierzowego MES dla zagadnień z nieliniowościami fizycznymi przyjmuje formę (2.25) =K(Q)Q F (2.25) Elementy macierzy K są zależne od kolumnowej macierzy przemieszczeń Q co prowadzi do konieczności rozwiązania nieliniowego układu równań. W celu rozwiązania układu równań stosowany jest przyrostowy opis problemu nieliniowego. Przyjmuje się, że obciążenie stopniowo narasta od zera do wartości maksymalnej w funkcji czasu. Czas jest wielkością przyjętą jedynie umownie w celu umożliwienia opisu przyrostu siły w formie inkrementalnej, natomiast rozwiązywane zagadnienie nadal jest statyczne. Rozwiązanie układu równań nieliniowych można rozpatrywać jako wielokrotne rozwiązywanie układu równań liniowych w kolejnych iteracjach dla zmieniającej się macierzy sztywności nazywanej styczną macierzą sztywności KT. Metoda Newtona-Raphsona (Rys. 2.19) jest najczęściej stosowaną metodą inkrementalno- iteracyjną rozwiązywania nieliniowych zagadnień MES. Na podstawie wartości naprężeń w elementach skończonych, w węzłach układu obliczane są wartości sił wewnętrznych. Następnie obliczana jest wartość siły niezrównoważonej będącej różnicą wartości sił zewnętrznych i sił wewnętrznych układu. Poszukiwany jest stan równowagi między siłami wewnętrznymi i zewnętrznymi. W następujących po sobie iteracjach obliczane są kolejne styczne macierze sztywności KT i , gdzie i to numer iteracji. Macierze te mają charakter przybliżony co powoduje lokalne niespełnienie stanu równowagi. Na początku analizy nieliniowej, wobec braku przemieszczeń oraz sił wewnętrznych, macierz styczna jest taka sama jak w zagadnieniu liniowym (KT 0 =K). Dla pierwszego inkrementu sił zewnętrznych, obliczane jest obciążenie niezrównoważone na podstawie którego wyznaczane są dodatkowe przemieszczenia węzłów, czyli zmiana konfiguracji mająca na celu osiągnięcie stanu równowagi. Dla zaktualizowanej konfiguracji obliczana jest kolejna macierz styczna KT 1 i ponownie obliczane jest obciążenie niezrównoważone oraz wyznaczane są dodatkowe przemieszczenia węzłów. Obliczenia przebiegają iteracyjnie w obrębie inkrementuj siły, do momentu osiągnięcia założonej dokładności, po czym następuje kolejny inkrement sił zewnętrznych. 36 Rys. 2.19 Metoda Newtona-Raphsona [32] W przypadku zmodyfikowanej metody Newtona-Raphsona (Rys. 2.20) styczna macierz sztywności obliczana jest tylko raz w obrębie inkrementu. Zmniejsza to koszty obliczeniowe związane z odwracaniem stycznych macierzy sztywności lecz skutkuje zwiększeniem liczby iteracji. Rys. 2.20 Zmodyfikowana metoda Newtona-Raphsona [32] W metodzie Newtona-Raphsona możliwe jest również stosowanie inkrementów przemieszczenia zamiast inkrementów siły. W niniejszej pracy zastosowano komercyjne oprogramowanie metody elementów kończonych – MSC.Marc. Narzędzie to umożliwia sporządzanie własnych podprogramów (ang. user subroutines) w języku Fortran. Zapewnia to duże możliwości modyfikacji sposobu prowadzenia obliczeń lub też analizy określonych wielkości na poziomie każdego inkrementuj czy też iteracji metody Newtona-Raphsona. Podejście to łączy zalety korzystania z profesjonalnego oprogramowania MES (zaawansowane modele materiału, warunki kontaktu 37 itp.) z elastycznością jaką daje możliwość tworzenia własnego kodu (np. definiowanie własnych modeli materiału, niestandardowych warunków brzegowych itp.). 2.8 Biorusztowania i przebudowa tkanki kostnej Inżynieria tkankowa (ang. tissue engineering) to interdyscyplinarna dziedzina nauk będąca połączeniem wiedzy biologicznej oraz inżynierii, mająca na celu projektowanie oraz wytwarzanie funkcjonalnych zamienników uszkodzonych tkanek lub narządów [85]. Badania realizowane w ramach inżynierii tkankowej obejmują między innymi prace dotyczące funkcjonalnych zamienników chrząstki [152], naczyń krwionośnych [143], skóry [107], mięśni [26], pęcherza moczowego [96] oraz tkanki kostnej. Biorusztowania (ang. scaffolds) [55] [61] to przestrzenne struktury stanowiące szkielet implantu (Rys. 2.21) w którego strukturze mogą zostać osadzane komórki pacjenta namnożone poza organizmem (in vitro) oraz czynniki wzrostu mające na celu przyśpieszenie procesu leczenia. Rys. 2.21 Przykładowe struktury biorusztowań [146] Kość charakteryzuje się wysoką zdolnością do regeneracji i odbudowy. Adaptacyjna przebudowa kości (ang. adaptive bone remodelling) jest ciągłym procesem zachodzącym zarówno w tkance gąbczastej i zbitej. W ciągu jednego roku przebudowie podlega około 10% szkieletu człowieka. Na przebudowę kości składają się dwa równoległe procesy tzn. resorpcja kości (ubytek masy tkanki kostnej) wywołana działaniem osteoklastów oraz osteogeneza (odbudowa) związana z działaniem komórek kościotwórczych – osteoblastów. Charakter równowagi między tymi procesami warunkuje dalszy przyrost, ubytek lub 40 implantu. Właściwa topografia powierzchni biorusztowania wspomaga migrację komórek do wnętrza implantu. Występowanie mikroporów (<10μm) ułatwia zagnieżdżenie się w strukturze biorusztowania nowych komórek. Ponadto biorusztowanie powinno charakteryzować się odpowiednimi parametrami materiałowymi, co jest szczególnie istotne w przypadku implantów tkanki kostnej przenoszącej znaczne obciążenia mechaniczne (np. tkanka kostna umiejscowiona w głowie kości udowej). W początkowym okresie leczenia (po implantacji biorusztowania) implant musi być zdolny do przenoszenia całości obciążeń wynikających z aktywności pacjenta. Istotne jest aby pacjent możliwie jak najwcześniej podejmował aktywność fizyczną i choćby w ograniczonym stopniu obciążał mechanicznie obszar w którym zaimplantowano biorusztowanie kości. Aktywność pacjenta stanowi w tym przypadku niezbędny bodziec mechaniczny indukujący proces przebudowy kości i leczenia. W dalszych etapach leczenia i rehabilitacji, nowo utworzona tkanka kostna wrastająca w strukturę biorusztowania, stopniowo przenosi coraz większą część obciążenia mechanicznego. Czas biodegradacji biorusztowania kości beleczkowej jest zależny od miejsca i obszaru implantacji, stopnia aktywności pacjenta, wieku pacjenta oraz własności samego biorusztowania i wynosi do kilku miesięcy. Porowatość struktury biorusztowania kości jest istotnym parametrem warunkującym skuteczność procesu leczenia. W przypadku biorusztowań kości gąbczastej, minimalna porowatość implantu wynosi około 60% [70] [74]. Z kolei minimalny wymiar porów to 100 μm [14], jednak optymalny wymiar porów skutkujący prawidłową przebudową tkanki oraz jej unaczynieniem wynosi co najmniej 300 μm [70] [114]. Górne granice porowatości oraz wymiaru porów biorusztowania uwarunkowane są koniecznością zachowania prawidłowych parametrów materiałowych, zapewniających możliwość przenoszenia obciążeń mechanicznych przez strukturę implantu. Wysoki stopień porowatości oraz parametry materiałowe umożliwiające przenoszenie dużych obciążeń, stanowią przeciwstawne wymagania dotyczące struktury biorusztowania kości beleczkowej. Ponadto parametry mechaniczne biorusztowania powinny być dopasowane do konkretnego pacjenta i miejsca implantacji, imitując parametry materiałowe usuniętej tkanki (biomimetyka) [164] [89]. Zbyt sztywna struktura biorusztowania nie pobudzi procesu przebudowy w tkance kostnej otaczającej implant, natomiast zbyt podatna struktura nie spełni swej funkcji wsporczej [148]. Parametry tkanki kostnej są zmienne dla różnych osób (są zależne od wieku, płci, trybu życia itp.) i nie są jednakowe nawet w obrębie tej samej kości. Biorusztowanie uwzględniające zmienność parametrów materiałowych tkanki kostnej 41 jest implantem spersonalizowanym wpisującym się w koncepcję medycyny spersonalizowanej [128]. Geometria mikrostruktury biorusztowania oraz zastosowany biomateriał wpływają na parametry mechaniczne implantu. Geometria implantu w skali makro, czyli kształt zewnętrzny biorusztowania jest osobnicza i musi być dopasowana do nieregularnego kształtu ubytku kostnego. Konieczność kontroli geometrii mikrostruktury oraz zewnętrznego kształtu biorusztowania warunkują dobór odpowiedniej metody wytwarzania implantu osobniczego. Metody szybkiego prototypowania (ang. Rapid Prototyping – RP) umożliwiają wytwarzanie implantów o precyzyjnie kontrolowanej geometrii z użyciem biodegradowalnych biomateriałów [115]. Biorusztowania wytworzone z zastosowaniem technologii RP stanowią wieloskalowe struktury hierarchiczne, podobnie jak ma to miejsce w przypadku struktury kości beleczkowej. 42 3 Modelowanie wieloskalowe w mechanice W przypadku klasycznych zagadnień mechaniki ośrodków ciągłych zazwyczaj przyjmowany jest model jednorodnego kontinuum materialnego rozpatrywanego materiału, a budowa mikrostruktury ośrodka zostaje pominięta. Często zakładana jest również izotropowość własności rozpatrywanego materiału. Dla materiałów konstrukcyjnych takich jak np. stal o regularnej budowie krystalicznej, powyższe założenia mogą stanowić wystarczające przybliżenie, tak więc izotropowy model ośrodka ciągłego przekłada się na wystarczająco dokładne wyniki analiz stanu przemieszczenia, odkształcenia oraz naprężenia. Jednak w przypadku nowoczesnych materiałów inżynierskich takich jak materiały kompozytowe, pianki metalowe itd., jak również dla struktur biologicznych takich jak tkanka kostna, model jednorodnego kontinuum materialnego może stanowić zbyt daleko idące uproszczenie, przekładające się na nieprawidłowe wyniki analiz strukturalnych. Metody uśredniania parametrów materiałowych i homogenizacji materiałów niejednorodnych początkowo były stosowane głównie w modelowaniu materiałów kompozytowych. Wymagana jest znajomość parametrów materiałowych poszczególnych elementów składowych rozpatrywanej mikrostruktury (np. inkluzji i osnowy). Najprostszą metodą obliczania efektywnych parametrów materiałowych struktur wielofazowych, np. kompozytów składającej się z materiału osnowy oraz umieszczonej w niej inkluzjach jest metoda mieszanin uwzględniająca jedynie procentowy udział powyższych składników mikrostruktury. Z użyciem metod analitycznych czy też półanalitycznych modelowana jest zazwyczaj pojedyncza, elipsoidalna inkluzja otoczona materiałem matrycy przez co często zakładany jest brak interakcji między poszczególnymi inkluzjami. Przy założeniu dla całego analizowanego obszaru mikrostruktury, jednorodnego stanu odkształcenia (model Voigta [172]) lub jednorodnego stanu naprężenia (model Reussa [132]) reprezentującego stan odkształcenia lub naprężenia dla danego punktu w skali makro, otrzymywane są odpowiednio górne i dolne granice oszacowanych parametrów zastępczych. Jedną z szerzej stosowanych analitycznych lub półanalitycznych metod obliczeniowych jest metoda Mori Tanaki [112] należąca do klasy modeli średniego pola (ang. mean field) opartych na rozwiązaniu Eshelby’ego [39] dla elipsoidalnej inkluzji otoczonej nieskończoną matrycą. Metoda ta pozwala na analizę materiałów kompozytowych o udziale inkluzji do około 30% objętości kompozytu. Następną grupą metod służących do analizy materiałów niejednorodnych 45 Wraz ze wzrostem mocy obliczeniowej komputerów coraz większe znaczenie zaczęła odgrywać metoda homogenizacji numerycznej oparta o metodę elementów skończonych [47] lub metodę elementów brzegowych (MEB) [68] [129]. Dla modelu mikrostruktury materiału (RVE) reprezentującego punkt materialny w skali makro, uśrednianiu podlegają naprężenia i odkształcenia co pozwala na wyznaczenie zastępczych parametrów materiałowych ekwiwalentnego materiału jednorodnego. Metody modelowania wieloskalowego są obecnie szeroko stosowane w modelowaniu nanokompozytów [141] lub piezokompozytów [41] czy też nieliniowych zagadnień plastyczności [180] [142]. Zagadnienia modelowania pól sprzężonych są rozpatrywane w ujęciu wieloskalowym [179] [158]. Procesy technologiczne takie jak obróbka plastyczna stali oraz zjawiska wpływu procesów chłodzenia materiałów obrabianych plastycznie na ich mikrostrukturę modelowane są w ujęciu wieloskalowym [126] [127]. Modele wieloskalowe stosowane są przy rozwiązywaniu zadań identyfikacji stochastycznych parametrów materiałowych [82] oraz identyfikacji parametrów układów obciążonych jednocześnie termicznie i mechanicznie [34]. Zastosowania metod wieloskalowych w medycynie dotyczą modelowania tkanek oraz np. modelowania struktury sztucznego serca [78]. Metody wielosklaowe ponadto pozwalają na uwzględnienie w modelu numerycznym struktury materiału na poziomie atomowym [19]. Ze względu na skomplikowaną, niejednorodną geometrię i budowę tkanek w różnych skalach obserwacji, metoda homogenizacji numerycznej oparta o MES stanowi adekwatną technikę wieloskalową dla modelowania struktur biologicznych. 3.1 Modelowanie wieloskalowe struktur kostnych Modelowanie numeryczne tkanki kostnej powiązane jest z metodą elementów skończonych i jej rozwojem. MES pozwala na modelowanie skomplikowanej i nieregularnej geometrii kości w różnych skalach oraz stosowanie różnych modeli materiałów. Zastosowanie MES w obszarze biomechaniki ortopedycznej sięga roku 1972 [15]. W początkowej fazie rozwoju modele numeryczne kości wykorzystywane były w badaniach nad remodelingiem kości oraz projektowaniu endoprotez, jednak charakteryzowały się uproszczoną geometrią oraz liniowością związków konstytutywnych. Wraz z rozwojem MES zaczęto stosować bardziej zaawansowane modele materiału oraz dokładniejsze geometrie 3D kości [60]. Modele numeryczne MES struktur kostnych mogą być wykorzystywane przy wspomaganiu 46 planowania operacji, badania struktur kostnych z zaimplantowanymi endoprotezami czy też w celu badania biomechaniki kręgosłupa człowieka [48]. W powszechnie stosowanych modelach numerycznych kości, własności mechaniczne tkanki charakteryzowane są przez moduł Younga (E) oraz współczynnik Poissona (ν). Wartości parametrów materiałowych można pozyskać na podstawie badań eksperymentalnych, jednak założenie niezmienności modułu Younga w całym obszarze tkanki gąbczastej np. dla struktury bliższego końca kości udowej jest zbyt daleko idącym uproszczeniem. Kość jest tkanką niejednorodną której mikrostruktura (skala beleczek kostnych) dostosowuje się do zewnętrznego obciążenia. Tak więc parametry materiałowe tkanki (moduł Younga) różnią się osobniczo oraz nawet w obrębie tej samej kości. Konwersja gęstości radiologicznej kości na parametry sprężyste (moduł Younga) tkanki z użyciem zależności empirycznych [53] [134] pozwala na uwzględnienie niejednorodnego rozkładu parametrów materiałowych kości beleczkowej i korowej w modelu numerycznym MES [91]. Istnieją próby sformułowania uniwersalnej zależności wiążącej gęstość kości z jej parametrami sprężystymi [31]. Pomiędzy istniejącymi zależnościami występują jednak pewne rozbieżności spowodowane między innymi różnicami w budowie mikrostrukturalnej kości dla różnych obszarów anatomicznych [111], brakiem jednolitych, usystematyzowanych metod badań eksperymentalnych próbek kostnych oraz faktem występowania niedoszacowań sztywności wynikających z badania próbek kostnych już poza organizmem (in vitro) [88] [71]. Niewystarczające rozmiary badanych próbek kostnych również wpływają niekorzystnie na wyniki badań eksperymentalnych, tak więc również na poprawność wyprowadzonych zależności empirycznych. Próbki kostne, jak i również modele numeryczne mikrostruktury kości beleczkowej powinny zawierać co najmniej 5 odległości międzybeleczkowych co przekłada się na minimalny zewnętrzny wymiar rozpatrywanego obszaru kości wynoszący około 4 mm [51]. Powyższe trudności wpływają na szeroki zakres modułów Younga obliczanych dla danej gęstości kości na podstawie różnych zależności empirycznych. Jednakże wykorzystanie ilościowej tomografii komputerowej (qCT) oraz adekwatnych równań konwertujących pozwala na dobre odwzorowanie niejednorodnej dystrybucji parametrów materiałowych w modelu numerycznym struktur kostnych i uzyskanie dokładniejszych wyników analiz numerycznych, które powinny być poddane dalszej weryfikacji. Konwersja gęstości kości na moduły Younga na podstawie danych radiologicznych nie pozwala jednak na uwzględnienie anizotropii kości wynikającej bezpośrednio z budowy mikrostruktury tkanki kostnej. 47 Zbudowanie na podstawie danych mikrotomogarficznych (μCT) modelu numerycznego bliższego końca kości udowej uwzględniającego mikrostrukturę na poziomie beleczek kostnych w całym rozpatrywanym obszarze jest wykonalne [169], jednak przekłada się na bardzo długie czasy obliczeń oraz trudności w modelowaniu geometrii mikrostruktury. Zazwyczaj analizie poddawane są modele numeryczne stosunkowo niewielkich próbek kostnych przy tak wysokim stopniu dokładności odwzorowania mikrostruktury kości (analizy Micro-FE) [167] [66]. Modele wieloskalowe [18] [99] pozwalają na uwzględnienie rzeczywistej mikrostruktury tkanki dla dużych obszarów kości przy zachowaniu akceptowalnych czasów obliczeń. Jedne z pierwszych periodycznych modeli uproszczonej, wyidealizowanej mikrostruktury kości beleczkowej przedstawiono w pracach [45] [46] jako układy o budowie prętowej, tworzące struktury porowate (Rys. 3.3). Rys. 3.3 Model geometryczny uproszczonej mikrostruktury kości beleczkowej [46] Hollister wraz z współautorami w pracy [54] zastosował metodę homogenizacji opartą o MES w celu zbadania wpływu różnych rodzajów uproszczonych topologii mikrostruktury kości na parametry efektywne kości gąbczastej. W pracy [58] przedstawiono koncepcję zastosowania obrazowania cyfrowego (np. CT) i teorii homogenizacji w celu badania struktury kości beleczkowej i jej parametrów. Zbudowano modele rzeczywistej mikrostruktury kości beleczkowej (Rys. 3.4) w oparciu o dane mikrotomograficzne. Rys. 3.4 Model numeryczny MES (RVE) mikrostruktury kości beleczkowej [58] 50 warunek brzegowy skutkuje przeszacowaniem, a naprężeniowy warunek brzegowy niedoszacowaniem zastępczych parametrów materiałowych struktury (Rys. 3.6) [149]. Rys. 3.6 Zbieżność obliczanych zastępczych parametrów materiałowych do parametrów efektywnych w zależności od zastosowanych warunków brzegowych i rozmiaru modelu RVE [79] Zastosowanie periodycznych warunków brzegowych przekłada się na otrzymanie wyników bliższych parametrom efektywnym niż ma to miejsce w przypadku zastosowania innych warunków brzegowych, nawet dla materiałów nie wykazujących periodyczności geometrycznej [156] [118]. Struktury periodyczne dzielimy na lokalnie periodyczne oraz globalnie periodyczne (Rys. 3.7). Z zastosowaniem metody homogenizacji możliwe jest obliczenie zastępczych parametrów materiałowych dla obydwu powyższych rodzajów struktur periodycznych [57]. Rys. 3.7 Przykłady struktury lokalnie periodycznej i globalnie periodycznej Struktury globalnie periodyczne w całym obszarze zbudowane są z powtarzającej się w każdym kierunku jednostki elementarnej (RVE). Charakter periodyczności w całej objętości materiału jest niezmienny, tak więc strukturę można zbudować poprzez kopiowanie 51 danej jednostki elementarnej w dwóch (obszar 2D) lub trzech (obszar 3D) kierunkach. Przykładami struktur globalnie periodycznych są materiały kompozytowe o regularnej budowie mikrostruktury czy też grafen. W przypadku struktur lokalnie periodycznych, charakter periodyczności zmienia się między poszczególnymi podobszarami w skali makro co przekłada się na większą liczbę modeli RVE opisujących mikrostrukturę ośrodka. Każdy z lokalnych podobszarów periodyczności zbudowany jest z odrębnej jednostki elementarnej (RVE). Przykładami struktur lokalnie periodycznych są inżynierskie materiały gradientowe i kość gąbczasta. Wyróżnić można również struktury periodyczne w sensie parametrów materiałowych, nie wykazujące periodyczności geometrycznej. W przypadku braku periodyczności geometrycznej, przeciwległe brzegi modelu RVE nie są identyczne, jednakże model RVE jest otoczony strukturą o takich samych parametrach materiałowych pod względem wartości czy też kierunkowości (anizotropia). Z użyciem periodycznych warunków brzegowych wymuszane są periodyczne przemieszczenia (3.2) oraz antyperiodyczne naprężenia (3.3) na przeciwległych brzegach modelu RVE (Rys. 3.8).  u u (3.2)   t t (3.3) Rys. 3.8 Periodycznie odkształcony model RVE Dla schematu przemieszczeniowego MES, periodyczne warunki brzegowe można zaimplementować w modelu RVE z wykorzystaniem równań liniowych (ang. Multi Point Constraints, MPC) wiążących stopnie swobody węzłów na przeciwległych brzegach modelu RVE. Zastosowanie periodycznych warunków brzegowych z użyciem równań MPC dla schematu przemieszczeniowego MES, zapewnia spełnienie warunku Hilla-Mandela. 52 Dla klasy zagadnień liniowo sprężystych, w celu obliczenia, zastępczych parametrów materiałowych materiału niejednorodnego, wystarczające jest jednokrotne przeprowadzenie sześciu analiz modelu RVE – rozwiązanie sześciu zagadnień brzegowych mechaniki z wykorzystaniem np. MES. W trakcie każdej z sześciu analiz numerycznych MES, do modelu RVE przykładane jest jednostkowe odkształcenie stanowiące jedną z sześciu składowych odkształcenia. Uśrednianiu podlegają obliczone naprężenia (3.4) i odkształcenia (3.5) periodycznie zdeformowanego modelu RVE.  RVE ij ij RVEV RVE 1 < ζ >= ζ dV V (3.4) gdzie <ζij> – uśrednione naprężenia dla skali makro ζij – naprężenia w skali mikro (RVE) VRVE – całkowita objętość modelu RVE  RVE ij ij RVEV RVE 1 < ε >= ε dV V (3.5) gdzie <εij> – uśrednione odkształcenia dla skali makro εij – odkształcenia w skali mikro (RVE) VRVE – całkowita objętość modelu RVE Na podstawie uśrednionych wartości naprężeń i odkształceń obliczane są kolejne kolumny macierzy sprężystości C (3.6) uogólnionego prawa Hookea. ij ij<ζ >= < ε >C (3.6) Macierz sprężystości stanowi w ogólności anizotropowe, zastępcze parametry materiałowe dla skali makro w punkcie materialnym reprezentowanym przez niejednorodny model RVE. Innymi słowy, macierz sprężystości stanowi parametry materiałowe ekwiwalentnego materiału jednorodnego obliczone na podstawie analiz niejednorodnego modelu RVE. Jednorodność materiału ekwiwalentnego w skali makro przekłada się na znaczne skrócenie czasów obliczeń numerycznych MES poprzez redukcję liczby stopni swobody modelu numerycznego. Informacja o budowie mikrostruktury ośrodka jest zachowana poprzez niejednorodne modele RVE. 55 immunologicznych. Zidentyfikowane wielkości zostają wprowadzone jako parametry modelu materiału dla jednorodnej struktury makro, modelowanej z użyciem metody elementów skończonych. W celu analizy stanu odkształcenia i naprężenia w mikrostrukturze (próbce numerycznej) należy stan odkształcenia dla danego punktu lub punktów modelu makro, zadać jako warunki brzegowe w modelu mikro. Jest to procedura analogiczna do procedury lokalizacji dla materiału liniowo-sprężystego. Metoda została zastosowana w celach modelowania wieloskalowego materiałów hipersprężystych [154] czy też materiałów kompozytowych wzmocnionych włóknami [155]. Dla niektórych nieskomplikowanych modeli mikro, możliwe jest określenie nieliniowych modeli zastępczych wprost na podstawie numerycznych testów materiałowych. Zależy to w dużej mierze od przyjętego modelu nieliniowości w skalach mikro i makro. 3.4 Algorytmy ewolucyjne w zastosowaniach wieloskalowych Identyfikacja i optymalizacja dla wybranych modeli wieloskalowych tkanek oraz implantów została zrealizowana z użyciem algorytmu ewolucyjnego [80]. Algorytmy ewolucyjne są szeroko stosowane dla zagadnień optymalizacji mechanicznych układów jednoskalowych m.in. w zagadnieniach sprężystości [77] i termosprężystości [33], optymalizacji struktury materiałów kompozytowych [5] czy też optymalizacji charakterystyk dynamicznych zaawansowanych układów mechanicznych [108]. Algorytmy ewolucyjne [109] oparte są na mechanizmach zaczerpniętych z biologicznej ewolucji gatunków. Generowana losowo populacja początkowa składająca się z osobników, poddawana jest selekcji na podstawie wartości funkcji przystosowania. Wartość funkcji przystosowania jest określana osobno dla każdego osobnika i stanowi stopień jego przystosowania do otaczającego środowiska, które jest reprezentacją przestrzeni zmiennych projektowych. Zazwyczaj każdy osobnik składa się z pojedynczego chromosomu którego składowymi są geny zawierające wartości zmiennych projektowych, natomiast funkcja przystosowania jest określona na podstawie funkcji celu procesu optymalizacji bądź identyfikacji. Proces selekcji wyłania osobniki najlepiej przystosowane, a więc grupę rozwiązań najbliższych poszukiwanemu optimum. Osobniki wyłonione w procesie selekcji poddawane są działaniu operatorów ewolucyjnych: krzyżowaniu oraz mutacji. Dla nowopowstałych osobników ponowne obliczane są wartości funkcji przystosowania. Cały algorytm jest przetwarzany iteracyjnie do momentu spełnienia warunku zatrzymania – braku zmiany wartości 56 przystosowania najlepszego osobnika w określonej liczbie następujących po sobie iteracji, osiągnięcie maksymalnej określonej liczby iteracji lub maksymalnego czasu obliczeń. Po zakończeniu działania algorytmu ewolucyjnego, osobnik o najlepszej wartości funkcji przystosowania przyjmowany jest za rozwiązanie postawionego problemu [1]. Schemat blokowy przedstawiający etapy działania algorytmu został przedstawiony na Rys. 3.11. Rys. 3.11 Schemat blokowy algorytmu ewolucyjnego 57 4 Algorytm homogenizacji numerycznej dla wieloskalowych modeli tkanki kostnej Algorytm homogenizacji numerycznej zastosowano do modelowania wieloskalowego porowatej struktury kości beleczkowej. Zbudowano modele numeryczne MES uproszczonej struktury kości beleczkowej [54] w celu zweryfikowania działania algorytmu i dokładności obliczanych wielkości. Wartości przemieszczeń oraz naprężeń obliczone z użyciem modeli wieloskalowych porównano z wartościami otrzymanymi z jednoskalowych modeli dokładnych i następnie wyznaczono błąd względny zrealizowanych obliczeń wieloskalowych [97]. 4.1 Weryfikacja algorytmu homogenizacji numerycznej dla zagadnień liniowych Zbudowano model numeryczny MES niejednorodnej i porowatej, uproszczonej struktury kości beleczkowej, reprezentujący próbkę kostną o wymiarach 10 x 10 x 10 mm (Rys. 4.1). W modelu zastosowano ośmiowęzłowe elementy skończone HEX8 o liniowej funkcji kształtu. Rys. 4.1 Niejednorodny model makro uproszczonej struktury kości beleczkowej Podczas budowy modelu numerycznego zastosowano szereg uproszczeń geometrii mikrostruktury kości beleczkowej:  beleczka kostna reprezentowana jest przez prostopadłościan o wysokości 1 mm i grubości 0,2 mm – uśrednione wymiary pojedynczej beleczki kostnej [30],  beleczki kostne łączą się pod kątem prostym,  struktura jest w pełni powtarzalna (globalna periodyczność struktury). 60 a) b) Rys. 4.5 Symulacje numeryczne próby ściskania dla modelu a) niejednorodnego b) wieloskalowego Po przeprowadzeniu symulacji prób ściskania, porównano obliczone wartości przemieszczeń maksymalnych. Przemieszczenia górnej płyty ściskającej dla modelu niejednorodnego i wieloskalowego wynoszą odpowiednio 45,72 μm (Rys. 4.6a) oraz 45,70 μm (Rys. 4.6b). Rozpatrując wynik otrzymany dla niejednorodnego modelu jednoskalowego jako dokładną wartość odniesienia, błąd względny obliczonych przemieszczeń dla modelu wieloskalowego wynosi 0,04%. a) b) Rys. 4.6 Wyniki symulacji ściskania - przemieszczenia dla modelu a) niejednorodnego b) niejednorodnego Aby obliczyć wartości naprężeń w mikrostrukturze (beleczkach kostnych) dla modelu wielkoskalowego należy przeprowadzić procedurę lokalizacji – stan odkształcenia dla danego punktu makroskopowego obszaru jednorodnego, przyłożyć do periodycznego modelu RVE. Porównano wartości naprężeń redukowanych zgodnie z hipotezą wytężeniową Hubera dla obszaru naroża modeli makro. Obliczone naprężenia w beleczkach kostnych dla modelu niejednorodnego i wieloskalowego wyniosły odpowiednio 30,21 MPa (Rys. 4.7a) oraz 61 31,17 MPa (Rys. 4.7b). Ponownie rozpatrując wynik otrzymany dla niejednorodnego modelu jednoskalowego jako dokładną wartość odniesienia, błąd naprężeń obliczonych w mikro skali dla modelu wieloskalowego wynosi 3,2%. a) b) Rys. 4.7 Wyniki symulacji ściskania - naprężenia redukowane dla modelu a) niejednorodnego b) wieloskalowego Końcowym etapem weryfikacji było porównanie czasów obliczeń w przypadku symulacji numerycznych próby ściskania dla modeli jednoskalowego oraz wieloskalowego. Zestawienie czasów obliczeń (ang. wall time) oraz liczby stopni swobody analizowanych modeli zamieszczono w Tab. 4.1. Tab. 4.1 Liczba stopni swobody oraz czasy obliczeń dla analizowanych modeli Zaobserwowano znaczne, ponad 20-krotne przyśpieszenie obliczeń w przypadku modelu wieloskalowego, przy jednoczesnym zachowaniu dokładności obliczonych przemieszczeń oraz naprężeń. Zweryfikowano zastosowane procedury obliczeń wieloskalowych oraz ich implementację co umożliwia realizację dalszych obliczeń dla wieloskalowych modeli struktur porowatych takich jak rzeczywista mikrostruktura kości beleczkowej. 62 4.2 Weryfikacja algorytmu homogenizacji numerycznej dla zagadnień nieliniowych W celu zweryfikowania modelu wieloskalowego z nieliniowościami fizycznymi, geometrię modelu uproszczonej mikrostruktury kości beleczkowej poddano modyfikacjom (Rys. 4.8) aby otrzymać materiał o ortortopowych zastępczych parametrach materiałowych w skali wyższej. Rys. 4.8 Model RVE uproszczonej struktury kości beleczkowej o zakładanych ortotropowych parametrach zastępczych W modelu wprowadzono periodyczne warunki brzegowe. Dla materiału beleczek kostnych (skala mikro) zastosowano nieliniowo sprężysty izotropowy model materiału hiposprężystego (2.7). Wartości parametrów materiałowych dla nieliniowego modelu materiału (Rys. 4.9b) określone zostały z wykorzystaniem informacji o spadku modułu Younga dla wartości odkształcenia 0.2% i 0.4% (Rys. 4.9a) zamieszczonej w artykule [110]. a) b) Rys. 4.9 a) Zależność wartości modułu Younga w funkcji odkształcenia b) nieliniowa zależność naprężenia od odkształcenia materiału beleczek 3500 3750 4000 4250 4500 4750 5000 0 0,001 0,002 0,003 0,004 0 4 8 12 16 20 24 0 0,002 0,004 0,006 0,008 0,01 E [MPa] ε ε ε σ [MPa] 65 Przyjęto 18 parametrów określających zależność 9 stałych materiałowych dla ortotropowego modelu hiposprężystego w skali makro od odkształcenia, zgodnie ze wzorem (4.3) p-b ε ij pd = a e (4.3) gdzie ap, bp – parametry krzywej wykładniczej określającej zmienność stałej materiałowej dij w funkcji odkształcenia Wartości parametrów ap oraz bp zamieszczono w Tab. 4.2. Tab. 4.2 Zastępcze parametry materiałowe uproszczonej struktury ortotropowej z nieliniowościami fizycznymi Dla obliczonych parametrów materiałowych przeprowadzono analizy porównawcze z użyciem modelu dokładnego (Rys. 4.11a) który powstał poprzez dziesięciokrotne powielenie mikrostruktury dla każdego z trzech kierunków układu współrzędnych, oraz modelu jednorodnego (Rys. 4.11b) dla którego zastosowano ortotropowe, zastępcze hiposprężyste parametry materiałowe. Z użyciem modeli przeprowadzono symulację numeryczną próby ściskania przy takich samych warunkach brzegowych. 66 a) b) Rys. 4.11 Symulacje numeryczne próby ściskania dla modeli o parametrach ortotropowych z nieliniowościami fizycznymi a) model jednorodny b) model wieloskalowy Po przeprowadzeniu symulacji prób ściskania, porównano obliczone wartości przemieszczeń maksymalnych. Przemieszczenia górnej płyty ściskającej dla modelu z nieliniowościami fizycznymi niejednorodnego jednoskalowego oraz wieloskalowego wynoszą odpowiednio 27,55 μm (Rys. 4.12a) oraz 25,55 μm (Rys. 4.12b). Rozpatrując wynik otrzymany dla niejednorodnego modelu jednoskalowego jako dokładną wartość odniesienia, błąd względny obliczonych przemieszczeń dla modelu wieloskalowego wynosi 7,2%. a) b) Rys. 4.12 Wyniki symulacji ściskania - przemieszczenia dla modelu o parametrach ortotropowych a) niejednorodnego b) niejednorodnego Po przeprowadzeniu lokalizacji dla modelu wieloskalowego, porównano wartości naprężeń redukowanych zgodnie z hipotezą wytężeniową Hubera dla obszaru naroża modelu makro. 67 Obliczone naprężenia w beleczkach kostnych dla niejednorodnego modelu jednoskalowego i wieloskalowego z nieliniowościami fizycznymi wyniosły odpowiednio 16,48 MPa (Rys. 4.13a) oraz 16,16 MPa (Rys. 4.13b). Ponownie rozpatrując wynik otrzymany dla niejednorodnego modelu jednoskalowego jako dokładną wartość odniesienia, błąd naprężeń redukowanych obliczonych w skali mikro dla modelu wieloskalowego wynosi 1,9%. a) b) Rys. 4.13 Wyniki symulacji ściskania - naprężenia redukowane dla modelu z nieliniowościami fizycznymi a) niejednorodnego b) wieloskalowego Porównano czasy obliczeń symulacji numerycznych próby ściskania dla modeli jednoskalowego oraz wieloskalowego z nieliniowościami fizycznymi. Zestawienie czasów obliczeń oraz liczby stopni swobody analizowanych modeli zamieszczono w Tab. 4.3. Tab. 4.3 Liczba stopni swobody oraz czasy obliczeń dla analizowanych modeli z nieliniowościami fizycznymi 70 5.1 Identyfikacja i optymalizacja dla wybranych modeli wieloskalowych tkanek i implantów Zadanie identyfikacji może zostać postawione jako problem optymalizacji. Celem identyfikacji dla modeli wieloskalowych tkanek jest określenie wektora (chromosomu) ch (5.2) zawierającego poszukiwane parametry materiałowe (zmienne decyzyjne) w skali mikro (RVE), które minimalizują funkcję celu FIDNT(ch) (5.1) zależną od wartości pomiarów eksperymentalnych i wyników symulacji numerycznych ( ) IDNT exp numF | w - w ( )|ch ch (5.1) gdzie wexp – wartość zmierzona eksperymentalnie wnum – wartość określona na drodze symulacji numerycznej przeprowadzonego eksperymentu ch=[x1, x2, …, xn] (5.2) xi to poszukiwane parametry materiałowe (geny). W przypadku identyfikacji minimalna wartość funkcji celu jest znana i wynosi zero. Zakres poszukiwanych wartości genów (parametrów materiałowych) dla różnego rodzaju materiałów inżynierskich lub tkanek biologicznych może być znany a priori. Ograniczenia nierównościowe wartości zmiennych projektowych w postaci (5.3) 𝑥𝑖 𝑚𝑖𝑛 ≤ 𝑥𝑖 ≤ 𝑥𝑖 𝑚𝑎𝑥 , 𝑖 = 1,2,… ,𝑛 (5.3) są stosowane w celu poprawy zbieżności algorytmu i przyśpieszenia jego działania, gdzie xi min oraz xi max to odpowiednio minimalna i maksymalna możliwa wartość zmiennej xi w procesie identyfikacji. Celem optymalizacji dla modeli wieloskalowych implantów jest określenie wektora (chromosomu) ch zawierającego parametry opisujące geometrię implantu (zmienne projektowe) w skali mikro (RVE), które minimalizują funkcję celu FOPT(ch) (5.4). Wartość funkcji celu jest zdefiniowana jako suma różnic między zastępczymi parametrami materiałowymi tkanki pacjenta i zastępczymi parametrami materiałowymi implantu zależnymi od wartości chromosomu ch, obliczonymi na podstawie analiz numerycznych modelu RVE. 71 ˆ( )  np OPT i i i=1 F c - cch (5.4) Gdzie ˆ ic to wartość referencyjna parametru materiałowego (tkanki), ic to wartość parametru materiałowego optymalizowanej struktury implantu, natomiast np to liczba porównywanych parametrów materiałowych. Dodatkowe ograniczenia mogą zostać zaimplementowane w postaci (5.5) ograniczeń nierównościowych gk<0 (5.5) gdzie k – numer ograniczenia gk – funkcja ograniczająca 5.2 Identyfikacja osobniczych parametrów materiałowych beleczek kostnych na podstawie danych eksperymentalnych i numerycznych W celu obliczenia efektywnych parametrów materiałowych porowatej struktury kości beleczkowej dla skali makro, konieczne jest wcześniejsze wyznaczenie wartości parametrów materiałowych struktury w skali mikro – modułu Younga oraz współczynnika Poissona beleczek kostnych tworzących mikrostrukturę kości. W skali poniżej mikro (skala sub-mikro) beleczki kostnej zbudowane są z blaszek kostnych (Rys. 2.5). Ze względu na przestrzenne ułożenie blaszek kostnych, można wnioskować o transwersalnej izotropii własności mechanicznych pojedynczej beleczki kostnej [170]. Zastosowanie uproszczenia jakim jest przyjęcie modelu izotropowego dla materiału beleczek przy założeniu jednorodnego rozkładu parametrów materiałowych w skali mikro ma pomijalny wpływ na obliczane parametry zastępcze kości beleczkowej dla skali makro [67] [168]. Powyższe uproszczenie znane w literaturze jako effective isotropic tissue modulus znajduje zastosowanie w całości obliczeń i modeli przestawionych w niniejszej pracy. Modelowanie wieloskalowe kości beleczkowej poczynając od skali niższych (np. struktura włókien kolagenowych, struktura blaszek kostnych) niż skala beleczek kostnych (skala mikro) jest możliwe do zrealizowana [50]. Pozyskanie osobniczych parametrów materiałowych struktury w skalach niższych niż mikro jest jednak problematyczne, a mierzone wielkości mogą być obarczone dużym błędem pomiarowym przekładającym się na nieprawidłowe wartości obliczanych parametrów zastępczych dla skal wyższych. Dwuskalowe modelowanie numeryczne struktury kości 72 beleczkowej (skala beleczek kostnych oraz skala całej kości) stanowi wystarczający stopień dokładności odwzorowania struktury rzeczywistej i pozwala na obliczanie osobniczych, anizotropowych parametrów zastępczych kości oraz dalszą optymalizację struktury spersonalizowanych implantów tkanki kostnej. W celu pozyskania parametrów materiałowych beleczek kostnych przeprowadzane są badania doświadczalne przy założeniu liniowo sprężystego izotropowego modelu materiału beleczek oraz zazwyczaj jednorodnego rozkładu wartości parametrów materiałowych. Najczęściej przeprowadzane badania doświadczalne, prowadzone na beleczkach kostnych w celu estymacji parametrów materiałowych to badania ultradźwiękowe [133], próby rozciągania, ściskania [10] i trójpunktowego zginania [90] oraz próba nanoindentacji [76]. Odnotowany w literaturze zakres modułu Younga dla materiału beleczek kostnych zawiera się w szerokim przedziale 1–15 GPa [92]. Powyższy zakres odzwierciedla trudne do wyeliminowania błędy pomiarów eksperymentalnych, wynikające z niewielkiego rozmiaru próbek oraz ich nieregularnego kształtu. Badania doświadczalne tkanki kostnej w skali mikro znajdują się nadal w fazie rozwoju, nie opracowano jak dotąd spójnych procedur przygotowania próbek beleczek kostnych. Należy również podkreślić osobniczość parametrów kości (również w skali mikro) oraz zmienność tych parametrów w zależności od miejsca występowania badanej beleczki kostnej. Użycie uśrednionych parametrów materiałowych beleczek kostnych, zaczerpniętych np. z danych literaturowych, w modelu numerycznym mikrostruktury kości beleczkowej skutkuje obliczaniem parametrów zastępczych jedynie architektury kości beleczkowej, a nie rzeczywistych osobniczych parametrów zastępczych dla danej próbki kostnej. W celu określenia parametrów materiałowych beleczek kostnych dla danej osoby oraz obszaru kości stosowane jest połączenie wyników badań eksperymentalnych oraz wyników analiz numerycznych modeli zbudowanych na podstawie danych μCT [84]. Rozpatrywana próbka kostna została poddana próbie ściskania na maszynie wytrzymałościowej MTS Insight. Dokładny pomiar przemieszczeń na powierzchni próbki zapewniło zastosowanie metody cyfrowej korelacji obrazu (DIC) w czasie trwania badania eksperymentalnego [76]. Następnie przeprowadzono symulację numeryczną zrealizowanej próby ściskania z użyciem modelu numerycznego MES badanej próbki kostnej (Rys. 5.3). 75 Każdy osobnik algorytmu ewolucyjnego reprezentowany jest przez pojedynczy chromosom zawierający dwa geny (5.7). Pierwszy gen chromosomu reprezentuje moduł Younga beleczek kostnych (ETRAB), natomiast drugi gen współczynnik Poissona materiału beleczek kostnych (νTRAB). Wartości parametrów materiałowych (genów) zapisane są w reprezentacji zmiennoprzecinkowej. ch=[ETRAB, νTRAB] (5.7) Dla wartości identyfikowanych parametrów materiałowych zastosowano ograniczenia nierównościowe (5.8)(5.9). 1000 MPa < ETRAB < 15000 MPa (5.8) 0.2 < νTRAB < 0.4 (5.9) Wartości parametrów materiałowych występujące w powyższych ograniczeniach określono na podstawie danych literaturowych [92]. Obliczenie wartości funkcji przystosowania dla pojedynczego osobnika (chromosomu) wymaga każdorazowego przeprowadzenia symulacji numerycznej próby ściskania (Rys. 5.3) co przekłada się na wysokie koszty i czasy obliczeniowe, zwłaszcza biorąc pod uwagę iteracyjne działanie algorytmu ewolucyjnego. W celu zredukowania czasów obliczeń zastosowano wariant algorytmu ewolucyjnego z dwoma podpopulacjami, realizujący obliczenia w sposób równoległy [17] [81]. W każdej iteracji algorytmu, najlepszy osobnik (o najlepszej wartości funkcji przystosowania) z danej podpopulacji migruje do drugiej podpopulacji w celu zapewnienia szybszej zbieżności obliczeń. Zastosowano następujące parametry algorytmu ewolucyjnego w procesie identyfikacji:  liczba podpopulacji – 2,  liczba osobników w każdej podpopulacji – 10,  prawdopodobieństwo mutacji Gaussa – 0.9,  prawdopodobieństwo mutacji prostej – 0.1,  warunek zatrzymania algorytmu – brak zmiany wartości funkcji celu w 10 kolejnych iteracjach 76 Zidentyfikowane wartości (geny najlepszego osobnika) parametrów materiałowych beleczek kostnych zamieszczono w Tab. 5.1. Tab. 5.1 Zidentyfikowane parametry materiałowe beleczek kostnych Zidentyfikowana wartość modułu Younga materiału beleczek kostnych jest bliska wartości modułu Younga obliczona na drodze eksperymentalnej z zastosowaniem metody nanoindentacji. Średnia wartość modułu Younga tkanki kostnej określoną z użyciem metody nanoindentacji wynosi 8.4 GPa [76]. Różnica 6% pomiędzy wartościami określonymi na drodze eksperymentu oraz na drodze identyfikacji może wynikać między innymi z punktowego charakteru pomiarów metodą nanoindentacji (zbadana została jedna beleczka próbki kostnej). Parametry tkanki określone na drodze identyfikacji z użyciem modelu numerycznego można rozpatrywać jako uśrednione parametry beleczek kostnych całej próbki kostnej. Różnica na poziomie 6%, mając na uwadze również skalę modelu i eksperymentu, jest wartością niską w przypadku materiału biologicznego takiego jak tkanka kostna. Zastosowanie modelowania numerycznego oraz algorytmów ewolucyjnych pozwoliło ponadto na identyfikację wartości współczynnika Poissona materiału kości w skali mikro, co było by znacznie utrudnione do zrealizowania na drodze samego badania eksperymentalnego. Wyniki badania eksperymentalnego w skali makro (próba ściskania) w postaci rozkładu przemieszczeń na powierzchni próbki kostnej (Rys. 5.6a) otrzymanego z zastosowaniem metody cyfrowej korelacji obrazu oraz wynik symulacji numerycznej próby ściskania (Rys. 5.6b) (dla zidentyfikowanych parametrów materiałowych beleczek kostnych) również w postaci rozkładu przemieszczeń na powierzchni próbki są podobne pod względem zarówno wartości przemieszczeń jak i rozkładów przemieszczeń. 77 a) b) Rys. 5.6 Rozkłady przemieszczeń na powierzchni próbki kostnej a) wynik badania eksperymentalnego b) wynik symulacji numerycznej Przeprowadzono identyfikację parametrów materiałowych kości beleczkowej w skali mikro (beleczki kostne) na podstawie pomiarów eksperymentalnych przeprowadzonych w skali makro oraz symulacji numerycznych. Pokonano w ten sposób trudności napotykane w trakcie prowadzenia badań eksperymentalnych w skali mikro oraz zredukowano błędy pomiarowe. Przedstawiona metodologia nie wymaga stosowania specjalistycznego sprzętu pomiarowego oraz ekstrakcji pojedynczych beleczek kostnych z próbki kostnej [90]. Zastosowana metoda cyfrowej korelacji obrazu (DIC) może zostać zastąpiona inną dokładną metodą pomiaru przemieszczeń, np. metodą elektronicznej interferencji obrazów plamkowych (ang. Electronic Speckle Pattern Interferometry, ESPI) [8]. 5.3 Homogenizacja struktury kości beleczkowej – obliczenie zastępczych parametrów materiałowych tkanki kostnej Po przeprowadzeniu identyfikacji parametrów materiałowych w skali mikro (beleczki kostne) dla danego obszaru kości oraz osoby, możliwe jest obliczenie parametrów materiałowych kości dla skali makro (skala całej kości) z zastosowaniem algorytmu homogenizacji numerycznej i periodycznych warunków brzegowych. Model MES mikrostruktury kości (Rys. 5.1) stanowi lokalny, osobniczy model RVE tkanki kostnej. Pomimo, że rozpatrywany model nie jest strukturą periodyczną w sensie geometrycznym, zastosowanie periodycznych warunków brzegowych przekłada się na otrzymanie poprawnych wyników dla wymiarów modelu RVE spełniających założenie kontinuum materialnego dla skali wyższej [86]. Jednakże w przypadku modelu niejednorodnej mikrostruktury kości, wprowadzenie periodycznych warunków brzegowych z zastosowaniem równań MPC nie jest 80 W podejściu self-consistent parametry materiałowe strefy buforowej nie są znane a priori i są określane na podstawie analiz obszaru mikrostruktury czyli rdzenia [13]. Początkowe wartości parametrów materiałowych dla jednorodnej strefy buforowej zostały określone na podstawie wyników próby ściskania próbki kostnej (E=1184 MPa) reprezentowanej przez model numeryczny rdzenia. Dla tak określonych parametrów startowych przeprowadzono obliczenia z zastosowaniem algorytmu homogenizacji numerycznej. W trakcie analiz numerycznych, strefa buforowa jest wykluczana z obliczeń zastępczych parametrów rdzenia [13], tak więc jednostkowe odkształcenie zadawane jest w obszarze rdzenia. Naprężenia również są uśredniane jedynie w obszarze mikrostruktury tkanki. Na podstawie przeprowadzonych analiz numerycznych MES obliczono zastępcze parametry sprężyste (macierzy sprężystości C) rdzenia (5.10) czyli mikrostruktury kości beleczkowej.                    55.171 -9.371 -22.603 35.655 -25.572 0.454 2 1896.975 510.523 566.842 1684.765 513.809 2453.461 = sym. 443.817 466.048 490.968 0.115 -48.886 -35.507 -3.688 -10.968 15.899 C (5.10) Obliczona macierz C jest symetryczna. Elementy macierzy sprężystości oznaczone kolorem szarym nie są równe zeru, jednak są mniejsze co najmniej o około rząd wielkości od pozostałych elementów macierzy. Można więc stwierdzić, że kość beleczkowa jest materiałem anizotropowym bliskim ortotropii. Analizowana sześcienna próbka kostna na podstawie której zbudowano model numeryczny, została wycięta z głowy kości udowej zgodnie z kierunkami anatomicznymi występującymi w tym obszarze (Rys. 2.7). Nie możliwe jest jednak uniknięcie choćby niewielkiego błędu podczas mechanicznej ekstrakcji próbki kostnej z głowy kości udowej. W celu zredukowania błędu ekstrakcji, przeprowadzono transformację początkowego układu współrzędnych związanego z wyciętą próbką kostną, tak aby uzyskać zastępcze parametry materiałowe jak najbliższe ortotropii (5.11) w postaci macierzy transformowanej C ORT [166] 81                    11 12 13 14 15 16 22 23 24 25 26 33 34 35 36 44 45 46 55 56 66 c c c δ δ δ c c δ δ δ c δ δ δ = sym. c δ δ c δ c ORT C (5.11) gdzie δij jest liczbą o jak najmniejszej wartości. Początkowy układ współrzędnych został obrócony (Rys. 5.10) względem osi x,y,z przyjętego kartezjańskiego układu współrzędnych (kąty Eulera) (5.12) z użyciem macierzy transformacji R (5.13). Wyrażenie wiążące macierze sprężystości przed i po transformacji przedstawia się następująco (zastosowana notacja sumacyjna)   ijkl iα jβ kγ lδ αβγδ C = R R R R C (5.12) gdzie  αβγδ C – początkowa macierz sprężystości w notacji tensorowej.  ijkl C – macierz sprężystości w notacji tensorowej po transformacji układu współrzędnych. R – macierz transformująca zawierająca kąty Eulera.          cosαcosχ + sinα sinβ sinχ sinαcosβ -cosα sinχ + sinα sinβ cosχ = -sinαcosχ + cosα sinβ sinχ cosαcosβ sinα sinχ + cosα sinβ cosχ cosβ sinχ -sinβ cosβ cosχ R (5.13) Rys. 5.10 Transformacja układu współrzędnych poprzez rotacje o kąty Eulera α, β, γ [29] 82 W celu obliczenia reprezentacji macierzy sprężystości najbliższej ortotropii, minimalizowano funkcję celu (5.14).   6 6 2 ij i=1 j=1 6 6 2 ij i=1 j=1 δ F( ) = c ch (5.14) gdzie ch = [α, β, γ] jest chromosomem którego geny zawierają kąty transformacji α, β i γ. Minimalizację przeprowadzono z użyciem algorytmu ewolucyjnego [102]. W wyniku transformacji otrzymano macierz (5.15) =                    1894.121 506.982 567.333 1682.723 516.523 2459.042 sym. 445.404 54.833 -25.734 -18.961 39.405 -25.161 2.677 7.634 -10.718 -36.137 -6.089 465. -12.251 15.4482 489. 3 11 5 6 ORT 1 C (5.15) dla zidentyfikowanych kątów transformacji zamieszczonych w Tab. 5.2. Tab. 5.2 Wartości zidentyfikowanych kątów transformacji Niskie wartości zidentyfikowanych kątów Eulera wskazują na prawidłową ekstrakcję próbki kostnej z głowy kości udowej – zgodnie z kierunkami anatomicznymi (łukami gotyckimi). Na podstawie obliczonych niskich wartości kątów transformacji można również wnioskować o zgodności kierunków anatomicznych z głównymi kierunkami ortotropii w głowie kości udowej człowieka. Niskie wartości zidentyfikowanych kątów przekładają się również na niewielkie zmiany wartości elementów macierzy sprężystości próbki kostnej. W przypadku próbki, która nie zostałaby wycięta z głowy kości udowej zgodnie z kierunkami anatomicznymi, zmiany wartości elementów macierzy sprężystości byłyby adekwatnie większe. Można wnioskować, że przedstawiona metodologia umożliwia prowadzenie ekstrakcji sześciennych próbek kości beleczkowej bez konieczności zachowania zgodności z kierunkami anatomicznymi, jednak teza ta wymaga dalszej weryfikacji. 85 Dla modelu geometrycznego przedstawionego na Rys. 5.1 zrealizowano również obliczenia zastępczych parametrów materiałowych przyjmując hiposprężysty model materiału beleczek kostnych z nieliniowościami fizycznymi. Ze względu na brak danych eksperymentalnych dotyczących parametrów materiałowych kości beleczkowej skali mikro uwzględniających nieliniowości fizyczne, zastosowano parametry beleczek kostnych sugerując się danymi zamieszczonymi w artykule [110]. Przyjęto zależność modułu Younga beleczek kostnych od odkształcenia zgodnie z zależnością (5.17): -89.17εE = 7842.2e (5.17) Z użyciem modelu mikrostruktury tkanki kostnej przeprowadzono serię numerycznych testów materiałowych, na podstawie których określono zastępcze parametry materiałowe tkanki kostnej dla modelu materiału hiposprężystego z nieliniowościami fizycznymi. Wykresy zależności stałych materiałowych materiału hiposprężystego od odkształcenia zamieszczono na Rys. 5.11. 86 a) b) c) d) e) f) g) h) i) Rys. 5.11 Wykresy uśrednionych hiposprężystych stałych materiałowych kości beleczkowej w funkcji odkształcenia a) d11 b) d22 c) d33 d) d12 e) d23 f) d31 g) d44 h) d55 i) d66 1000 1200 1400 1600 1800 0,001 0,004 0,007 0,01 800 1000 1200 1400 1600 0,001 0,004 0,007 0,01 1400 1600 1800 2000 2200 2400 0,001 0,004 0,007 0,01 300 350 400 450 500 0,001 0,004 0,007 0,01 320 370 420 470 0,001 0,004 0,007 0,01 300 350 400 450 500 0,001 0,004 0,007 0,01 380 400 420 440 0,01 0,04 0,07 0,1 380 400 420 440 460 0,001 0,004 0,007 0,01 420 440 460 480 0,001 0,004 0,007 0,01 ε ε ε ε ε ε ε ε ε d33 d12 d23 d31 d44 d55 d66 d11 d22 87 Parametry określające zależność 9 stałych materiałowych dla modelu hiposprężystego w skali makro od odkształcenia, przyjętych zgodnie ze wzorem (4.3) zamieszczono w Tab. 5.4. Tab. 5.4 Zastępcze parametry materiałowe struktury kości beleczkowej z nieliniowościami fizycznymi Ze względu na brak danych eksperymentalnych dotyczących rzeczywistych parametrów materiałowych beleczek kostnych uwzględniających nieliniowości fizyczne, w dalszej części pracy (Rozdział 5.5) anizotropowe zastępcze parametry materiałowe dla liniowego modelu materiału (macierz C4 ORT ) zostały użyte w zadaniu optymalizacji topologicznej struktury osobniczego biorusztowania kości [100] [83] jako parametry materiałowe tkanki kostnej chirurgicznie usuniętej z głowy kości udowej. Planowane jest przeprowadzenie badań doświadczalnych pozwalających na wyznaczenie parametrów materiałowych beleczek kostnych z uwzględnieniem nieliniowości fizycznych i następne wyznaczenie zastępczych parametrów materiałowych kości beleczkowej dla skali makro. 5.4 Trójskalowy model biorusztowania kości Biorusztowanie dopasowane do danego pacjenta oraz miejsca implantacji musi charakteryzować się odpowiednią topologią mikrostruktury oraz kształtem dopasowanym do nieregularnej geometrii ubytku kostnego powstałego w wyniku zabiegu chirurgicznego. Uwzględnienie powyższych wymagań przekłada się na wybór adekwatnej metody wytwarzania biorusztowania spersonalizowanego. Metodą umożliwiającą precyzyjną kontrolę geometrii wytwarzanego implantu jest metoda Fused Deposition Modeling (FDM) należąca 90 Rys. 5.14 Struktura biorusztowania w zależności od skali obserwacji Analiza numeryczna MES struktury biorusztowania o rozmiarze ubytku kostnego, z zachowaniem stopnia dokładności na poziomie mikro (osadzone ścieżki) dla całego obszaru implantu nie jest możliwe, szczególnie w przypadku zadań optymalizacji i identyfikacji. Ze względu na metodę wytwarzania (FDM), beleczki biorusztowania zbudowane z ścieżek osadzonego tworzywa termoplastycznego są strukturą globalnie periodyczną. Na przekrojach poprzecznych beleczek implantu, można wyróżnić periodyczne obszary RVE. Na podstawie pomiarów mikroskopowych zgładów (Rys. 5.15a), zbudowano numeryczny model RVE (Rys. 5.15b) mikrostruktury beleczek biorusztowania wytworzonego metodą FDM. a) b) Rys. 5.15 a) Fragment przekroju poprzecznego beleczki biorusztowania b) Model RVE mikrostruktury Z użyciem modeli numerycznych struktury implantu dla różnych skali obserwacji oraz algorytmu homogenizacji numerycznej i periodycznych warunków brzegowych, zbudowano model trójskalowy (Rys. 5.16) struktury biorusztowania wytworzonego metodą FDM. 91 Rys. 5.16 Model trójskalowy biorusztowania wytworzonego metodą FDM Periodyczne warunki brzegowe zostały zaimplementowane w modelu RVE mikrostruktury (Rys. 5.15b) i następnie z użyciem metody homogenizacji numerycznej obliczono transwersalnie izotropowe parametry materiałowe beleczek biorusztowania – analiza mikro- meso. Tak obliczone parametry stosowane są jako parametry materiałowe beleczek biorusztowania (skala meso), które modelowane są jako ośrodek jednorodny o zastępczych transwersalnie izotropowych parametrach materiałowych. W modelu meso, ponownie z użyciem metody homogenizacji numerycznej i periodycznych warunków brzegowych, obliczane są zastępcze parametry materiałowe dla skali wyższej – analiza meso-makro. Po dwukrotnym przeprowadzeniu procedury homogenizacji obliczone zostają anizotropowe parametry materiałowe biorusztowania dla skali makro, które powinny być bliskie anizotropowym parametrom materiałowym chirurgicznie usuniętej tkanki kostnej, której funkcje implant czasowo zastępuje. Szereg parametrów charakteryzujących strukturę biorusztowania w skalach mikro i meso ma wpływ na wartości zastępczych parametrów materiałowych implantu w skali makro. Parametry struktury biorusztowania dla skali mikro to:  parametry materiałowe zastosowanego biomateriału.  geometria osadzanych ścieżek materiału termoplastycznego (średnica dyszy urządzenia FDM) .  głębokość wtopienia kolejnych warstw materiału (ang. negative air gap) zależna od temperatury dyszy oraz liniowej prędkości osadzania,  wzajemna orientacja ścieżek tworzywa w następujących po sobie warstwach materiału – raster. 92 Parametry struktury biorusztowania dla skali meso to:  porowatość,  wymiary porów,  grubości beleczek,  liczba beleczek dla trzech kierunków przyjętego układu współrzędnych. Sterując powyższymi parametrami możliwe jest uzyskanie i następnie wytworzenie zoptymalizowanej struktury biorusztowania kości, dopasowanego dla danego pacjenta i miejsca implantacji. Tak zaprojektowane biorusztowanie spersonalizowane indukuje proces przebudowy kości przy jednoczesnym zapewnieniu funkcji wsporczej dla szkieletu pacjenta. 5.5 Optymalizacja topologiczna struktury spersonalizowanego biorusztowania kości Przeprowadzono optymalizację topologiczną biorusztowania kości beleczkowej z użyciem równoległego algorytmu ewolucyjnego [101]. Modelowano trójskalową strukturę implantu zbudowanego z biozgodnego i biodegradowalnego kompozytu biopolimeru PLGA z bioceramiką HA. Udział procentowy hydroksyapatytu [87] w strukturze materiału stanowi zmienną projektową. Z użyciem modelu mikro struktury biorusztowania (Rys. 5.15b) obliczono stałe sprężyste transwersalnie izotropowego materiału beleczek biorusztowania w zależności od udziału procentowego domieszki hydroksyapatytu (Rys. 5.17). Rys. 5.17 Parametry sprężyste beleczek biorusztowania w funkcji udziału hydroksyapatytu Obliczenia dwuskalowe przeprowadzono z użyciem algorytmu homogenizacji numerycznej. Zastosowano periodyczne warunki brzegowe dla modelu (mikro) ścieżek materiału 0 5000 10000 15000 20000 25000 0 5 10 c11 c12 c13 c22 c23 c33 c44 c55 c66 95 Rys. 5.19 Wartości funkcji celu dla podpopulacji w kolejnych teracjach Wartość funkcji celu wyniosła 886 MPa. Wynikiem optymalizacji jest chromosom (5.22) zawierający geny charakteryzujące geometrię biorusztowania spersonalizowanego o parametrach materiałowych bliskich parametrom materiałowym usuniętej tkanki kostnej. ch OPT = [3, 2, 4, 0.2, 0.14, 0.1] (5.22) Periodyczny model RVE biorusztowania jest generowany na podstawie wartości genów najlepszego chromosomu. Otrzymano model MES (Rys. 5.20) zbudowany z 46148 sześciennych elementów skończonych HEX8 o liniowej funkcji kształtu. Model RVE zoptymalizowanego biorusztowania kości przeskalowano w celu uzyskania minimalnego wymiaru porów równego 0,3 mm umożliwiającego wrastanie tkanki kostnej w strukturę implantu. Skalowanie geometrii nie zmienia wartości obliczonych zastępczych parametrów materiałowych struktury ze względu na liniowość modelu. Grubości beleczek periodycznego modelu RVE biorusztowania wynoszą odpowiednio 0.6 mm, 0.42 mm, oraz 0.3 mm dla kierunków x,y,z. Ilość beleczek zoptymalizowanego biorusztowania dla kierunków x,y,z to odpowiednio 3, 2 oraz 4 beleczki (Rys. 5.20). Ortotropowe parametry materiałowe (5.23) biorusztowania to                    SCAFF 2110.99 241.91 582.16 0 0 0 1479.71 295.17 0 0 0 2629.59 0 0 0 = 118.29 0 0 sym. 121.86 0 504.54 C (5.23) 96 Na Rys. 5.20 przedstawiono model MES zoptymalizowanej struktury biorusztowania kości gąbczastej. Rys. 5.20 Model RVE zoptymalizowanego osobniczego biorusztowania kości a) widok z przodu b) widok z lewej strony c) widok z góry d) widok izometryczny Porowatość zoptymalizowanego biorusztowania kości beleczkowej wynosi 63%. Ze względu na periodyczność struktury, możliwe jest powielenie modelu RVE w trzech kierunkach w celu otrzymania biorusztowania dopasowanego do indywidualnego kształtu ubytku kostnego (Rys. 5.21), przy zachowaniu kierunkowości parametrów materiałowych implantu. d) b) c) a) 97 Rys. 5.21 Przykładowa struktura biorusztowania kości dopasowana do kształtu ubytku kostnego Otrzymany model może zostać następnie eksportowany do formatu pliku STL oraz wytworzony z odpowiedniego biopolimeru przy użyciu urządzenia FDM. Zoptymalizowaną strukturę implantu (Rys. 5.20) zastosowano w wieloskalowej analizie bliższego końca kości udowej z zaimplantowanym biorusztowaniem (układ kość–implant). 5.6 Wieloskalowy model układu kość – implant spersonalizowany Przeprowadzono wieloskalową analizę numeryczną układu bliższego końca kości udowej człowieka (Rys. 5.22) z zaimplantowanym spersonalizowanym biorusztowaniem kości beleczkowej (układ kość–implant). Rys. 5.22 Model numeryczny bliższego końca kości udowej (model makro)