Univerzita Karlova v Praze Matematicko-fyzikální fakulta
DIPLOMOVÁ PRÁCE
Pavel Čížek Matematická analýza a numerická simulace remodelačních procesů v kostech Matematický institut Univerzity Karlovy
Vedoucí diplomové práce: Prof. Ing. František Maršík, DrSc. Studijní program: Matematické modelování ve fyzice a technice 2008
Děkuji panu Profesorovi Františkovi Maršíkovi za odborné vedení práce, za trpělivost a za podnětné připomínky. Děkuji panu Magistrovi Martinovi Mádlíkovi za pomoc při výpočtech a práci s programem Fstrin.
Prohlašuji, že jsem svou diplomovou práci napsal(a) samostatně a výhradně s použitím citovaných pramenů. Souhlasím se zapůjčováním práce a jejím zveřejňováním. V Praze dne 1.4.2008
Pavel Čížek
2
Obsah Seznam značení
7
1 Úvod 1.1 Stavba kosti . . . . . . . . . . . . . . . . . . . . . . . . . . .
9 9
2 Deformace kosti 2.1 Bilanční zákony . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Materiálové vlastnosti kosti . . . . . . . . . . . . . . . . . .
13 13 14
3 Slabá formulace stacionárního problému deformace kosti 3.1 Slabá formulace elastostatiky . . . . . . . . . . . . . . . . . 3.2 Problém s monotonnií . . . . . . . . . . . . . . . . . . . . . 3.3 Linearizace . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Specifikace okrajové podmínky . . . . . . . . . . . . . . . . . 3.5 Aplikace Lax - Milgramovy věty . . . . . . . . . . . . . . . . 3.5.1 Omezenost formy A˜ . . . . . . . . . . . . . . . . . . . 3.5.2 V -elipticita formy A˜ . . . . . . . . . . . . . . . . . .
18 18 19 21 22 23 24 26
4 Metody řešení 4.1 Ritzova metoda . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Triangulace . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Volba báze . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29 29 30 32
5 Implementace metody konečných prvků v jazyku C++ 5.1 Základní členění programu . . . . . . . . . . . . . . . . . . 5.2 Inicializace triangulace . . . . . . . . . . . . . . . . . . . . 5.3 Inicializace hraničních podmínek . . . . . . . . . . . . . . . 5.4 Inicializace bázových funkcí . . . . . . . . . . . . . . . . . 5.5 Výpočet přibližného řešení . . . . . . . . . . . . . . . . . .
36 36 37 39 40 41
3
. . . . .
5.6
Výhody a nevýhody objektového přístupu . . . . . . . . . .
6 Výsledky numerického řešení 6.1 Zatížení kosti . . . . . . . . . . . . 6.2 Zatížení kosti a chrupavky . . . . . 6.3 Deformace kosti uchycené v trubce 6.4 Testování modelu . . . . . . . . . . 6.5 Vztah deformace a remodelace kosti
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
42 44 44 46 48 50 55
7 Závěr
56
Literatura
57
4
Název práce: Matematická analýza a numerická simulace remodelačních procesů v kostech Autor: Pavel Čížek Katedra (ústav): Matematický institut Univerzity Karlovy Vedoucí diplomové práce: Prof. Ing. František Maršík, DrSc. e-mail vedoucího:
[email protected] Abstrakt: V předložené práci studujeme model deformace kosti v důsledku působení vnějších sil. Po biologickém úvodu, kde popisujeme remodelaci kosti a její spojitost s deformací, se zabýváme bilančními zákony a materiálovými vlastnostmi kosti. Problém formulujeme slabě. Dále rovnice linearizujeme a dokazujeme existenci a jednoznačnost řešení. Následně popisujeme metodu konečných prvků a naprogramování této metody pro výpočet přibližného řešení problému s využitím logiky objektově orientovaného programování. Uvádíme také popis využívaných datových struktur. V závěru prezentujeme numerické výsledky dosažené popisovaným programem zachycující deformaci stehenní kosti a deformaci obratlů a meziobratlové ploténky. Klíčová slova: remodelace kosti, deformace, metoda konečných prvků, objektově orientované programování Title: Mathematical analysis and numerical simulation of remodelling processes in bones Author: Pavel Čížek Department: Mathematical Institute of Charles University Supervisor: Prof. Ing. František Maršík, DrSc. Supervisor’s e-mail address:
[email protected] Abstract: In the presented work, we study a model of bone deformation caused by action of external forces. After the biological introduction, in which we describe a bone remodelling process and its connection with bone deformation, we deal with conservation laws and material characteristics of a bone. We define weak formulation. Further, we linearize the equations and prove existence and uniqueness of the solution. Consecutively, we describe the finite element method and the programming of this method for the computing of the approximate solution of the problem using the logic of object-oriented programming. We also adduce a description of the used 5
data structures. In the conclusion, we present numerical results reached by the described program that represent deformation of femur and deformation of vertebras and intervertebral foramina. Keywords: bone remodelling, deformation, finite element method, objectoriented programming
6
Seznam značení Tučné symboly značí vektory nebo body z prostoru R3 . a aij a ˜ A A˜ cijkl (n) ci Cij C 0 (Ω) ekl e˜kl E1 , E2 , E3 f f g G12 , G23 , G31 h hi H 1 (Ω) Ki n n Pk (Ki ) s suppψ Sym (Ω) tli
zobrazení odvozené z Hookova zákona ij-tá složka zobrazení a linearizované zobrazení a operátor odvozený ze zobrazení a bilineární forma odvozená ze zobrazení a˜ ijkl-tá složka tenzoru elastických konstanta i-tá konstanta lineární kombinace bázových funkcí ϕi prostoru Vn materiálová konstanta prostor spojitých funkcí na Ω kl-tá složka Eurelova tenzoru konečných deformací kl-tá složka tenzoru malých deformací materiálové konstanty lineární funkcionál odvozený z funkce g objemová síla funkce definující okrajovou podmínku materiálové konstanty maximum průměrů hi průměr i-tého elementu triangulace Th Hilbertův prostor funkcí na Ω i-tý element triangulace Th počet vrcholů triangulace Th jednotková vnější normála k hranici ∂Ω prostor polynomů k-tého stupně na Ki počet částí hranice ∂Ω nosič funkce ψ prostor symetrických matic, jejichž prvky jsou funkce na Ω li-tá složka Cauchyho tenzoru napětí 7
T Th u un v V Vn , Vh w W 1,2 (Ω) x
operátor stop triangulace polygonální oblasti Ωh posunutí přibližné řešení unášivá rychlost 3 podprostor prostoru [H 1 (Ω)] konečně dimenzionální podprostor prostoru V testovací funkce Sobolevův prostor funkcí na Ω bod prostoru R3 v prostorových souřadnicích
Γ1 Γ2 Γ3 δij λij νij ρ ρi ϕi Ω Ωh
část hranice ∂Ω s nulovou Dirichletovou okrajovou podmínkou část hranice ∂Ω s Neumannovou okrajovou podmínkou část hranice ∂Ω s nulovou Neumannovou okrajovou podmínkou Kroneckerovo delta j-tá barycentrická souřadnice i-tého elementu triangulace Th materiálové konstanty hustota maximum průměrů vepsaných koulí do i-tého elementu triangulace Th i-tá bázová funkce prostoru Vn , resp. Vh pracovní oblast polygonální aproximace oblati Ω
v˙
materiálová derivace funkce v časová derivace derivace podle l-té složky prostorové proměnné hranice oblasti Ω gradient funkce v
∂ ∂t ∂ ∂xl
∂Ω ∇v
8
Kapitola 1 Úvod Deformace kosti patří k základním fyzickým vlastnostem zdravé kosti. Od velikosti deformace se odvíjí její přestavba. Úvodem popíšeme biologickou stavbu kosti a dále budeme studovat její deformaci.
1.1
Stavba kosti
Kosti jsou základním stavebním orgánem opěrné soustavy lidského těla. Tvoří kostru, jejíž funkcí je podpírání (nesení) těla, ochrana ostatních orgánů a také krvetvorba. Kosti se rovněž podílejí na metabolismu minerálů, v kostře se ukládá několik důležitých prvků, např. vápník, fosfor, sodík a draslík. Kosti se skládají z vláken, která jsou orientována ve směru největšího a nejčastějšího zatěžování. Tak je zajišťována co největší nosnost a pevnost kosti. Rozlišujeme dva druhy kostní tkáně - kompaktní a houbovitou. Kompaktní kostní tkáň je na povrchu kostí, je pevná a hustá a tvoří 80% váhy kostry. Houbovitá kostní tkáň je uvnitř kostí, především kloubních hlavic, je řídká a má složitou strukturu s mnoha dutinami. Struktura tenkých vláken (trámců), která houbovitou kost tvoří, závisí na zatěžování kosti a vyvíjí se celý život vlivem tlaků a napětí tak, aby pevnost kosti byla co největší. Funkce systému trámců je zatížení částečně absorbovat a částečně rovnoměrně přenést na kompaktní kost, která je pevnější. Jsou-li trámce vhodně uspořádány, mají schopnost absorbovat poměrně vysoké množství energie, což je důležité např. při pádu.
9
Kompaktní i houbovitá kostní tkáň je tvořena třemi typy kostních buněk - osteoblasty, osteocyty a osteoklasty a mezibuněčnou hmotou zvanou kostní matrix. Osteoblasty jsou hlavní stavební buňkou kosti. Vytváří vrstvy složené z jedné řady buněk. Tyto vrstvy se tvoří na površích, kde se ukládá nová kost. Osteoblasty také produkují kostní matrix. Osteocyty vznikají z osteoblastů uvíznutím v dutinách mineralizované kostní matrix (lakunách). Dotýkají se navzájem dlouhými tenkými výběžky sloužícími k předávání živin a kyslíku a k odvodu odpadních látek. Přetvářejí a udržují okolní matrix. Osteoklasty resorbují starou kost, rozkládají kostní matrix a uvolňují z ní minerály. Kostní tkáň je schématicky znázorněna na obrázku 1.1, který je převzat z [6]. Fotografie mikroskopické struktury kostní tkáně je na obrázku 1.2. Písmeno A značí kostní matrix s osteocyty, písmeno B vrstvu osteoblastů a písmeno C resorbující osteoklasty.
Obrázek 1.1: Nákres kostní tkáně 10
Obrázek 1.2: Fotografie mikroskopické struktury kostní tkáně Kost má schopnost přetvářet architektoniku trámců v reakci na dlouhodobější změny zatěžování. Přetváření probíhá odstraňováním staré kosti a jejím nahrazováním kostí novou, která má strukturu lépe odolávající probíhajícímu zatěžování. Dochází k řídnutí kosti v oblastech s nízkou mechanickou zátěží a ke zpevňování v místech opakované vyšší zátěže. Tento proces se nazývá remodelace a je jednou ze základních vlastností zdravé kosti. Kromě mechanické zátěže má na remodelaci vliv mnoho dalších faktorů, např. frekvence zátěže, pohlaví, věk nebo přísun vitamínů, minerálů či hormonů. Remodelace kosti probíhá celý život, je to reakce organismu na změny některých z výše uvedených faktorů. 11
Remodelace kosti je složitý proces, na kterém se podílejí osteoklasty a osteoblasty. Nejprve osteoklasty resorbováním staré kosti vyhloubí tunel. Stěny tohoto tunelu pokrývají tenkou vrstvou osteoblasty, které začínají produkovat kostní matrix. Dále se osteoblasty podílejí na postupné mineralizaci kostní matrix, kterou se pomalu obklopují. Nakonec osteoblasty obklopené mineralizovanou kostní matrix dozrávají a přetvářejí se v osteocyty. Na následujících stránkách budeme studovat deformaci kosti v důsledku působení vnějších sil, neboli zatížení kosti chůzí či doskokem. Budeme sledovat posunutí materiálových elementů popisující absorbci zatížení v kosti. Při dlouhodobějším působení stejné zátěžové síly je způsob absorbce zatížení, tedy posunutí vnitřních elementů, jedním z hlavních faktorů ovlivňujících remodelační proces kosti. Studium vlivu zatížení kosti na její remodelaci je důležité např. v lékařství. U pacientů s protézou kyčelního kloubu je hlavním důvodem uvolňování implantátu převládající resorbce kostní tkáně v jeho okolí způsobená změnou absorbce zatížení kosti váhou těla. V této oblasti by mohly matematické výsledky pomoci lékařům predikovat nejvhodnější mechanické zatěžování kosti s implantátem.
12
Kapitola 2 Deformace kosti V biologickém úvodu jsme popsali mj. důvody studia deformace kosti v důsledku působení vnějších sil. Nyní uvedeme matematický model, který tuto situaci popisuje.
2.1
Bilanční zákony
Nejprve formulujeme zákony bilance v prostorových souřadnicích. Zanedbáním některých členů je postupně zjednodušíme a dojdeme tak k výsledné parciální diferenciální rovnici, kterou budeme dále studovat. Úlohu budeme řešit na množině Ω ⊂ R3 , o které předpokládáme, že je to omezená oblast s lipschitzovskou hranicí. Bilance hmotnosti:
Bilance hybnosti:
∂ρ ∂ + l ρv l = 0, ∂t ∂x ρv˙i +
∂tli = ρf i , ∂xl
x ∈ Ω.
x ∈ Ω.
Objemové síly f (gravitace) zanedbáme, neboť jejich působení deformaci kosti ovlivňuje minimálně. Hustotu ρ > 0 budeme uvažovat konstatní. Takto
13
se nám bilanční zákony zjednoduší na následující rovnice: ∂ l v = 0, x ∈ Ω, ∂xl ∂tli ρv˙i + l = 0, x ∈ Ω. ∂x Jak ukážeme níže, problém se transformuje na výpočet posunutí u. Unášivá rychlost v je v porovnání s posunutím u zanedbatelná, což demonstruje následující, pro jednoduchost jednorozměrný, příklad: Uvažujme periodické zatěžování kosti ve směru osy x, neboli chůzi, posunutí je tvaru u = u0 ei(kx−ωt) , kde u0 , k a ω jsou konstanty. Užitím vztahu v = ∂u = −iωu můžeme vyjádřit materiálovou derivaci rychlosti v násle∂t dovně: ∂v ∂v ∂ 2u ∂ 2u ∂u ∂u v˙ = +v = 2 − iωu = −iω − ω2u = −ω 2 u − ikω 2 u2 . ∂t ∂x ∂t ∂x∂t ∂t ∂x Nyní porovnejme absolutní hodnoty prvního a druhého členu: 2 ω u 1 |−ω 2 u| = 2 2 = . 2 2 |−ikω u | kω u k | u|
π Označíme-li L délku kosti, lze konstantu k vyjádřit vztahem k = 2L . Uvážíme-li dále, že délka L je mnohem větší než posunutí u, tedy že platí L ≫ | u|, dostáváme odhad: 1 2L = ≫ 1. k | u| π | u| Tento odhad ilustruje, že lze derivaci rychlosti v podle prostorové proměnné zcela zanedbat a materiálovou derivaci nahradit druhou časovou derivací posunutí u. My budeme ovšem uvažovat pouze stacionární případ, takže lze materiálovou derivaci rychlosti zcela zanedbat. Bilanční zákony jsme zjednodušili na rovnici:
∂tli = 0, ∂xl
2.2
x ∈ Ω.
(2.1)
Materiálové vlastnosti kosti
Protože uvažujeme isotermický proces lze Cauchyho tenzor napětí tli vyjádřit pomocí Hookova zákana následovně: tli = cijkl ekl , 14
kde ekl je Eulerův tenzor konečných deformací a cijkl je tenzor elastických konstant, který lze ze čtvrtého řádu převézt na symetrickou matici elastických konstant o rozměrech 6 × 6 využitím symetrií cijkl = cjikl = cijlk = cjilk ,
cijkl = cklij ,
které plynou ze symterie tenzoru napětí a ze záměnnosti pořadí parciálních derivací. Zaveďme označení: C11 = c1111 , C22 = c2222 , C44 = c2323 ,
C12 = C21 = c1122 , C23 = C32 = c2233 , C55 = c1313 ,
C13 = C31 = c1133 , C33 = c3333 , C66 = c1212 ,
jenž nám umožní psát Hookův zákon ve tvaru: 11 t C11 C12 C13 0 0 0 t22 C21 C22 C23 0 0 0 33 t C31 C32 C33 0 0 0 23 = t 0 0 0 C44 0 0 13 t 0 0 0 0 C55 0 t12 0 0 0 0 0 C66
e11 e22 e33 e23 e13 e12
.
(2.2)
Kost, jakožto příčně izotropní materiál, je charakterizována následujícími materiálovými konstantami: E1 = E2 = 15 GPa, E3 = 18 GPa, ν12 = ν21 = 0.2,
ν13 = ν23 = 0.25,
G23 = G31 = 6 GPa,
ν31 = ν32 = 0.3,
(2.3)
G12 = 4.5 GPa.
E1 , E2 a E3 jsou moduly pružnosti v tahu (Youngovy moduly) v materiálových směrech (x1 , x2 , x3 ), νij Poissonova čísla pro příčnou deformaci v j-tém směru při jednoosém napětí v i-tém směru a G23 , G31 a G12 jsou smykové moduly v rovinnách (x2 , x3 ), (x3 , x1 ) a (x1 , x2 ). S využitím označení ∆=
1 − ν12 ν21 − ν23 ν32 − ν31 ν13 − 2ν21 ν32 ν13 , E1 E2 E3
a symetrií ν21 ν12 = , E1 E2
ν13 ν31 = , E1 E3 15
ν23 ν32 = , E2 E3
jsou nenulové prvky matice elastických konstant v následujícím vztahu k uvedeným materiálovým konstantám Cij : C11 =
1 − ν23 ν32 , E2 E3 ∆
C22 =
C44 = G23 ,
1 − ν13 ν31 , E1 E3 ∆
C55 = G31 ,
C33 =
1 − ν12 ν21 , E1 E2 ∆
C66 = G12 ,
C23 = C32 =
ν23 + ν21 ν13 ν32 + ν12 ν31 = , E1 E3 ∆ E1 E2 ∆
C31 = C13 =
ν31 + ν21 ν32 ν13 + ν23 ν12 = , E2 E3 ∆ E1 E2 ∆
(2.4)
ν12 + ν32 ν13 ν21 + ν31 ν23 = . E2 E3 ∆ E1 E3 ∆ Eulerův deformační tenzor ekl lze vyjádřit pomocí vektoru posunutí u takto (podrobněji viz [4]): 1 ∂uk ∂ul ∂um ∂um . (2.5) ekl = + − 2 ∂xl ∂xk ∂xk ∂xl C21 = C12 =
Dosazením tohoto vztahu do rovnice (2.2) dostáváme: ∂uj ∂um ∂um 1 ii , Cij 2 − t = 2 ∂xj ∂xj ∂xj 1 ∂u2 ∂u3 ∂um ∂um 23 32 t =t = , C44 + − 2 ∂x3 ∂x2 ∂x2 ∂x3 1 ∂u1 ∂u3 ∂um ∂um 13 31 t =t = , C55 + − 2 ∂x3 ∂x1 ∂x1 ∂x3 1 ∂u1 ∂u2 ∂um ∂um 12 21 t =t = . C66 + − 2 ∂x2 ∂x1 ∂x1 ∂x2 Takto zapsaný Cauchyho tenzor napětí lze považovat za zobrazení def
3
aij (u (x)) = tij (x) , 3 a : H 1 (Ω) → Sym (Ω) ,
z prostoru [H 1 (Ω)] do prostoru Sym (Ω), kde symbolem Sym (Ω) myslíme prostor všech symetrických matic, jejichž prvky jsou funkce bodů oblasti Ω. 16
Bilančí rovnici (2.1) lze pak vyjádřit ve tvaru: ∂ali (u (x)) = 0, ∂xl
x ∈ Ω.
(2.6)
Přidejme okrajovou podmínku, která je dána funkcí g (x, t) udávající zatížení na povrchu kosti v bodě x a v čase t. Čas t budeme uvažovat pouze jako parametr funkce g a úlohu budeme řešit stacionárně. Označíme-li dále n jednotkovou vnější normálu k hranici ∂Ω má okrajová podmínka tvar (Neumannova okrajová podmínka): tli (x) nl = g i (x, t) ,
x ∈ ∂Ω.
Použijeme-li tvar rovnice (2.6) okrajová podmínka vypadá následovně: ali (u (x)) nl = g i (x, t) ,
17
x ∈ ∂Ω.
(2.7)
Kapitola 3 Slabá formulace stacionárního problému deformace kosti V této kapitole nejprve odvodíme slabou formulaci stacionární deformace kosti a po té rozebereme existenci a jednoznačnost řešení toho problému.
3.1
Slabá formulace elastostatiky
Nyní odvodíme slabou formulaci rovnice (2.6) s okrajovou podmínkou (2.7). 3 Hledáme posunutí u, o kterém předpokládáme, že je z prostoru [H 1 (Ω)] , 3 tedy Sobolevova prostoru [W 1,2 (Ω)] se skalárním součinem indukovaným normou, proto rovnici (2.6) skalárně vynásobíme libovolnou testovací funkcí 3 w ∈ [H 1 (Ω)] ∂ali (u (x)) i w = 0, ∂xl
3 ∀ w ∈ H 1 (Ω) .
Tuto rovnici zintegrujeme přes oblast Ω Z ∂ali (u (x)) i w dx = 0, ∂xl Ω
3 ∀ w ∈ H 1 (Ω) .
Na prostorový integrál použijeme Greenovu větu Z Z ∂wi li − a (u (x)) l dx + ali (u (x)) nl wi dx = 0, ∂x Ω
∂Ω
18
3 ∀ w ∈ H 1 (Ω) .
Ještě použijeme okrajovou podmínku. Neboť ve funkci g (x, t) je čas t pouze parametr, budeme pro přehlednost psát pouze g (x). Dostáváme slabou formulaci problému: Z Z 3 ∂wi li a (u (x)) l dx = g i (x) wi dx, ∀ w ∈ H 1 (Ω) . (3.1) ∂x Ω
3.2
∂Ω
Problém s monotonnií
Ve zbytku této kapitoly se budeme věnovat důkazu existence a jednoznačnosti řešení úlohy (3.1). Označme: Z ∂wi def hA (u) , wi = ali (u (x)) l dx. ∂x Ω
∗ 3 3 Takto je definován nelineární operátor A : [H 1 (Ω)] → [H 1 (Ω)] . Pro důkaz existence řešení rovnice (3.1) je třeba, aby tento operátor byl monotonní (podrobněji viz [2]). Definice (Monotonnie). Nechť V je Banachův prostor a nechť A : V → V ∗ je zobrazení. Pak toto zobrazení je monotonní právě když splňuje nerovnost hA (u) − A (v) , u − vi ≥ 0,
∀u, v ∈ V.
My však ukážeme, že operátor A monotonní není, že existují taková 3 u, v ∈ [H 1 (Ω)] , že platí opačná nerovnost než v uvedené definici, tedy nerovnost hA (u) − A (v) , u − vi < 0. Volme: u = (x, y, z) , v = c (x, y, z) , kde c ∈ R je konstanta. Ukážeme, že pro vhodně zvolené c platí nerovnost: Z ∂ui ∂v i li li dx < 0. (3.2) − a (u) − a (v) ∂xl ∂xl Ω
Integrujeme skalární součin dvou matic. Druhý člen součinu má tvar: 19
1 0 0 1−c 0 0 1−c 0 = − (c − 1) 0 1 0 . ∇u − ∇v = 0 0 0 1 0 0 1−c
Nyní rozepišme člen ali (u) C11 + C12 + C13 0 0 1 0 C21 + C22 + C23 0 ali (u) = 2 0 0 C31 + C32 + C33
a člen ali (v) ali (v) =
C11 + C12 + C13 0 0 1 . 0 C21 + C22 + C23 0 = 2c − c2 2 0 0 C31 + C32 + C33
První člen skalárního součinu je tedy tvaru: ali (u) − ali (v) =
C + C + C 0 0 11 12 13 1 0 C21 + C22 + C23 0 = 1 − 2c + c2 2 0 0 C31 + C32 + C33 C11 + C12 + C13 0 0 1 . 0 C21 + C22 + C23 0 (c − 1)2 = 2 0 0 C31 + C32 + C33
S využitím těchto výpočtů můžeme psát (”:” značí maticový skalární součin): Z Ω
a (u) − a (v) li
li
∂ui ∂v i − ∂xl ∂xl
dx =
C11 + C12 + C13 0 0 1 : 0 C21 + C22 + C23 0 = (c − 1)2 2 0 0 C31 + C32 + C33 Ω 3 1 0 0 X 1 3 dx = − (c − 1) |Ω| Cij . : −1 (c − 1) 0 1 0 2 i,j=1 0 0 1 Z
20
Prostým dosazením konstant (2.3) do vztahů (2.4) dostáváme Cij > 0, i, j = P3 1, 2, 3 , proto je také i,j=1 Cij > 0, zřejmě rovněž |Ω| > 0. Volíme-li tedy c > 1 platí nerovnost: 3 X 1 3 − (c − 1) |Ω| Cij < 0. 2 i,j=1
Takto jsme dokázali nerovnost (3.2), z čehož plyne, že operátor A není monotonní. Jedna z možností, jak dokázat existenci řešení úlohy (3.1) je zobrazení a linearizovat. Tento postup se pokusíme aplikovat.
3.3
Linearizace
Protože nelineární členy zobrazení a se znatelně projeví až v 15% deformacích, ke kterým v kosti při obvyklých situacích nedochází, tyto členy zanedbáme. Definujeme zobrazení a ˜, k čemuž využijeme tenzor malých deformací ∂ul 1 ∂uk e˜kl = 2 ∂xl + ∂xk . Hookův zákon (2.2) pak dostává tvar:
t11 t22 t33 t23 t13 t12
=
C11 C12 C13 0 0 0 C21 C22 C23 0 0 0 C31 C32 C33 0 0 0 0 0 0 C44 0 0 0 0 0 0 C55 0 0 0 0 0 0 C66
e˜11 e˜22 e˜33 e˜23 e˜13 e˜12
.
(3.3)
Dosadíme-li opět vyjádření tenzoru malých deformací do rovince (3.3), dostáváme vztahy: ∂uj , ∂xj ∂u2 ∂u3 , C44 + ∂x3 ∂x2 ∂u1 ∂u3 , C55 + ∂x3 ∂x1 ∂u1 ∂u2 , C66 + ∂x2 ∂x1
tii = Cij 1 2 1 = 2 1 = 2
t23 = t32 = t13 = t31 t12 = t21
21
3
definující zobrazení a ˜ z prostoru [H 1 (Ω)] do prostoru Sym (Ω) def
a ˜ij (u (x)) = tij (x) , 3 a ˜ : H 1 (Ω) → Sym (Ω) . Bilanční rovnice (2.6) pak dostává tvar: ∂˜ ali (u (x)) = 0, x ∈ Ω. ∂xl Podobně slabá formulace (3.1) dostává tvar: Z Z 3 ∂wi li a ˜ (u (x)) l dx = g i (x) wi dx, ∀ w ∈ H 1 (Ω) . ∂x Ω
3.4
(3.4)
(3.5)
∂Ω
Specifikace okrajové podmínky
Pokud více specifikujeme okrajovou podmínku lze u linearizovaného problému dokázat kromě existence také jednoznačnost řešení. Ke specifikaci okrajové podmínky využijeme Větu o operátoru stop. Věta (O operátoru stop). Buď Ω ⊂ Rn omezená oblast s lipschitzovskou hranicí. Pak existuje lineární omezený operátor T : W 1,p (Ω) → Lp (∂Ω) tak, že • T u = u⌊∂Ω , ∀u ∈ C Ω ∩ W 1,p (Ω), • existuje konstanta c > 0 tak, že kT ukp,∂Ω ≤ c kuk1,p,Ω ,
∀u ∈ W 1,p (Ω).
Oblast Ω představuje část kosti. Při výpočtech budeme za tuto oblast volit válec stojící (ve směru osy z) na jedné z podstav. Okrajová podmínka (2.7) popisuje zatížení působící na povrch kosti při chůzi. Je zřejmé, že toto zatížení působí pouze na svrchní podstavu. Předpokládejme zjednodušeně, že kost stojí na pevné podložce. Rozdělíme tedy hranici ∂Ω na části Γ1 , Γ2 a Γ3 . Přičemž platí Γ1 ∩ Γ 2 Γ1 ∩ Γ 3 Γ2 ∩ Γ 3 Γ1 ∪ Γ2 ∪ Γ3 22
= = = =
∅, ∅, ∅, ∂Ω.
Množina Γ1 představuje spodní podstavu, na které válec stojí. Zde předepíšeme nulovou Dirichletovu okrajovou podmínku. Část Γ2 je svrchní podstava, na které předepisujeme Neumannovu okrajovou podmínku danou funkcí g z původní okrajové podmínky (2.7). Na plášti válce, tedy na množině Γ3 , chování posunutí neznáme, tedy předepíšeme nulovou Neumannovu okrajovou podmínku. Okrajová podmínka pak dostává tvar u (x) = 0, x ∈ Γ1 , a ˜ (u (x)) nl = g i (x, t) , x ∈ Γ2 , a ˜li (u (x)) nl = 0, x ∈ Γ3 . li
(3.6)
3
Definujeme prostor V ⊂ [H 1 (Ω)] : n o 1 3 def V = w ∈ H (Ω) , T w |Γ1 = 0 ,
kde T je operátor stop z výše uvedené věty, podrobněji viz [2], přičemž předpokládáme platnost |Γ1 | > 0. Požadujeme-li, aby řešení u bylo z prostoru V , je zaručeno splnění okrajové podmínky na hranici Γ1 . Proto i při odvozování slabé formulace násobíme rovnici (3.4) testovací funkcí w z prostoru V, která má nulovou stopu na hranici Γ1 . Slabá formulace (3.5) pak dostává tvar Z Z ∂wi li a ˜ (u (x)) l dx = g i (x) wi dx, ∀ w ∈ V. (3.7) ∂x Ω
3.5
Γ2
Aplikace Lax - Milgramovy věty
K důkazu existence a jednoznačnosti řešení rovnice (3.7) s okrajovou podmínkou (3.6) použijeme Lax - Milgramovu větu. Věta (Lax-Milgramova). Nechť X je reálný Hilbertův prostor s normou k.k a prostor X ∗ je jeho duál. Nechť a : X × X → R je bilineární forma splňující (omezenost) (X-elipticita)
∃M > 0 : |a (u, v)| ≤ M kuk kvk , 2
∃α > 0 : |a (u, u)| ≥ α kuk , 23
∀u, v ∈ X,
∀u ∈ X,
a nechť dále f ∈ X ∗ , aplikaci prvku f na prvek u značíme (f, u). Pak má rovnice pro neznámou u ∈ X a (u, v) = (f, v)
∀v ∈ X
právě jedno řešení, které navíc splňuje kuk ≤
1 α
kf kX ∗ .
Důkaz věty je možno nalézt např. v [3]. Pracujeme na prostoru V, což je zřejme reálný vektorový prostor s nor3 mou a skalárním součinem shodným jako v [H 1 (Ω)] . Aby to byl Hilbertův prostor je však ještě třeba dokázat, že je úplný, což přímo plyne ze spojitosti operátoru stop T, která je dokázána ve Větě o operátoru stop (podrobněji viz [2]). Podobně jako v části 3.2 definujeme zobrazení A˜ (u, w): Z ∂wi def ˜ (3.8) a ˜li (u (x)) l dx. A (u, w) = ∂x Ω
Protože jsme ale v sekci 3.3 provedli linearizaci zobrazení a, je předpisem (3.8) definována bilineární forma A˜ : V × V → R. Aby jsme mohli aplikovat Lax-Milgramovu větu je třeba dokázat, že je tato bilineární forma omezená a 3 V -eliptická. Připomeňme ještě definici normy k.k prostoru [H 1 (Ω)] a tedy i prostoru V :
def kuk =
3 Z X i=1 Ω
1 2 2 Z ∂ui 2 |ui | dx + . ∂xj dx i,j=1 3 X
Výrazem
∂ui ∂xj
3.5.1
Omezenost formy A˜
Ω
myslíme slabou derivaci funkce u.
˜ Nechť u, w ∈ V jsou liboNejprve studujme omezenost bilineární formy A. volné funkce. Z Z Z ˜ ˜ (u) : ∇w dx ≤ |˜ a (u) : ∇w| dx ≤ k˜ a (u)k2 k∇wk2 dx, A (u, w) = a Ω
Ω
Ω
24
a ˜ (u) značí matici s prvky a ˜ij (u), takže integrujeme skalární součin dvou matic, proto jsme v poslední nerovnosti mohli použít Cauchyho nerovnost. Zavedeme označení: 1 1 1 def CM = max C44 , C55 , C66 , Cij , i, j = 1, 2, 3 , 2 2 2 a rozepíšeme člen k˜ a (u)k2 : k˜ a (u)k2 = 2 3 2 2 3 X X 1 2 ∂u2 ∂u3 1 2 ∂u1 ∂u3 ∂uj = Cij j + C44 3 + 2 + C55 3 + 1 + ∂x 2 ∂x ∂x 2 ∂x ∂x i=1 j=1 2 1 2 ∂u1 ∂u2 + C66 2 + 1 2 ∂x ∂x !2 3 X ∂u1 ∂u3 2 ∂u2 ∂u3 2 ∂u j 2 ≤ CM 3 ∂xj + ∂x3 + ∂x2 + ∂x3 + ∂x1 + j=1 # ∂u1 ∂u2 2 . + 2 + 1 ∂x ∂x
∂u ∂u k i , i 6= k, j 6= l, Po umocnění závorek dostaneme také členy tvaru 2 ∂x j ∂xl na které použijeme Youngovu nerovnost: ∂ui ∂uk ∂ui 2 ∂uk 2 2 j l ≤ j + l . ∂x ∂x ∂x ∂x Pak lze psát:
2
k˜ a (u)k ≤
2 6CM
3 X ∂ui 2 2 2 ∂xj = 6CM k∇uk2 .
i,j=1
25
Nyní pokračujme v odhadování celého integrálu: Z k˜ a (u)k2 k∇wk2 dx ≤ Ω
≤
√
6CM
Z
k∇uk2 k∇wk2 dx
Ω
1 1 2 2 Z Z √ 2 2 6CM k∇ukL2 (Ω) dx k∇wkL2 (Ω) dx ≤ ≤
√
Ω
Ω
6CM kukH 1 (Ω) kwkH 1 (Ω) .
Při přechodu od druhého ke třetímu řádku jsme použili Hölderovu nerovnost. Výsledná nerovnost má tvar √ ˜ A (u, w) ≤ 6CM kukH 1 (Ω) kwkH 1 (Ω) ,
˜ dokazující omezenost formy A.
3.5.2
V -elipticita formy A˜
K důkazu V -elipticity je třeba znát přesné hodnoty parametrů Cij . Ty dostaneme dosazením elastických konstant (2.3) do vztahů (2.4). Zaokrouhleno na jeden desetinný řád vypadají parametry Cij následovně: C11 =17.8, C22 =17.8, C33 =22.2,
C12 =5.3, C13 =6.9, C23 =6.9,
C44 =6, C55 =6, C66 =4.5.
Nyní již můžeme začít s důkazem C11 def B = C12 C13
elipticity. Zaveďmě označení: C12 C13 C22 C23 . C23 C33
(3.9)
O symetrické matici B dokážeme, že je pozitivně definitní. K tomu využijeme tvrzení, že matice s reálnými koeficienty je pozitivně definitní, právě když jsou kladné všechny její hlavní subdeterminanty, podrobněji viz [1]. Hlavní 26
subdeterminanty označme po řadě K, L, M . Z vyčíslených konstant (3.9) přímo vypočítáme, že hlavní subdeterminanty jsou kladné K = C11 > 0, 2 L = C11 C22 − C12 > 0, 2 2 2 M = C11 C22 C33 + 2C12 C13 C23 − C11 C23 − C12 C33 − C13 C22 > 0. Matice B je tedy symetrická pozitivně definitní, proto existuje konstanta α > 0 tak, že platí nerovnost: u · Bu ≥ α |u|2 ,
∀u ∈ R3 .
Rozepišme aplikaci bilineární formy A˜ na prvek u: ∂u1 2 Z ∂xl ∂u ∂u C ∂u ∂u ∂u 1 1 66 2 3 2 ∂u 2 ˜ A (u, u) = · B ∂x2 + , , + ∂xl ∂x2 ∂x3 2 ∂x2 ∂x1 ∂u3 Ω ∂x3 2 2 # C44 ∂u2 ∂u3 C55 ∂u1 ∂u3 dx. + + + + 2 ∂x3 ∂x1 2 ∂x3 ∂x2 Zavedeme konstantu Cm : def
Cm = min{α, C44 , C55 , C66 }. S využitím pozitivní definitnosti matice B lze psát: A˜ (u, u) ≥ 2 2 Z 3 3 X X ∂ui 1 ∂ui ∂uj ≥ + Cm + i α dx i j ∂x 4 ∂x ∂x i,j=1 i=1 Ω
≥ Cm
i6=j
2 X 2 3 ∂u ∂u 1 ∂u i j i i + dx + + j i ∂xi ∂xi 2 ∂x ∂x i,j=1
Z 3 X 1 ∂u Ω
2 i=1
i6=j
2 Z X 3 1 ∂ui ∂uj = Cm + i dx 2 ∂xj ∂x i,j=1 Ω Z Z 2 ≥ CK Cm k∇ukL2 (Ω) dx ≥ CF CK Cm kuk2H 1 (Ω) dx. Ω
Ω
27
V posledním řádku jsme po řadě využili Kornovu (konstanta CK ) (viz [5]) a Friedrichsovu nerovnost (konstanta CF ), které můžeme použít, protože pracujeme nad prostorem funkcí, které mají nulovou stopu na části hranice ∂Ω. Tím je dokončen důkaz V -elipticity bilineární formy A˜ (u, v). Aby jsme mohli aplikovat Lax-Milgramovu větu musíme ještě dokázat, že pravá strana rovnice (3.7) lze zapsat jako prvek prostoru V ∗ aplikovaný na prvek prostoru V . Definujeme lineární funkcionál f : Z def (f, w) = g i (x) wi dx, ∀ w ∈ V. Γ2
f ∈ V ∗ je třeba, aby pro všechny w ∈ V existoval integrál RK platnosti 3 i i g (x) w dx, což zaručíme přidáním předpokladu g ∈ [L2 (Γ2 )] , podrobΓ2 něji viz [3]. Nyní jsou splněny všechny předpoklady Lax-Milgramovy věty. Dokázali jsme tedy existenci a jednoznačnost řešení u rovnice (3.7) s okrajo3 vou podmínkou (3.6) a předpokladem g ∈ [L2 (Γ2 )] , o kterém navíc platí kuk ≤ CF C1K Cm kf kV ∗ . 3 Platnost g ∈ [L2 (Γ2 )] budeme dále předpokládat, aniž bychom to explicitně zmiňovali.
28
Kapitola 4 Metody řešení Pro výpočet řešení rovnice (3.7) s okrajovou podmínkou (3.6) použijeme metodu konečných prvků, která je odvozena z Ritzovy metody. Tuto metodu krátce teoreticky rozebereme.
4.1
Ritzova metoda
Úlohu řešíme na prostoru V, který je separabilní (podrobněji viz [7]), tedy existuje jeho spočetná báze. Zvolme některou bázi a označme ji ϕ1 (x1 , x2 , x3 ) , ϕ2 (x1 , x2 , x3 ) , . . . Zvolme dále pevné n a řešme úlohu (3.7) na konečně dimenzionálním prostoru Vn s bází ϕ1 (x1 , x2 , x3 ) , ϕ2 (x1 , x2 , x3 ) , . . . , ϕn (x1 , x2 , x3 ) . Dostaneme tak přibližné řešení un , které splňuje rovnici: A˜ (un , w) = (f, w) ,
(4.1)
∀ w ∈ Vn .
Protože přibližné řešení un leží v prostoru Vn , lze vyjádřit jako lineární kombinace bázových funkcí: (n)
un (x1 , x2 , x3 ) = c1 ϕ1 (x1 , x2 , x3 ) + · · · + cn(n) ϕn (x1 , x2 , x3 ) .
(4.2)
Dosadíme-li toto vyjádření přibližného řešení un do rovnice (4.1) a využijemeli linearitu formy A˜ v první složce, dostaneme rovnici: (n) ˜ c1 A˜ (ϕ1 , w) + · · · + c(n) n A (ϕn , w) = (f, w) ,
29
∀ w ∈ Vn .
Nyní dosadíme za testovací funkci w po řadě bázové funkce ϕ1 , . . . , ϕn . Zís(n) (n) káme tak soustavu lineárních algebraických rovnic pro neznámé c1 , . . . , cn . (n) (n) (n) A˜ (ϕ1 , ϕ1 ) c1 + A˜ (ϕ2 , ϕ1 ) c2 + · · · + A˜ (ϕn , ϕ1 ) cn = (f, ϕ1 ) (n) (n) (n) A˜ (ϕ1 , ϕ2 ) c1 + A˜ (ϕ2 , ϕ2 ) c2 + · · · + A˜ (ϕn , ϕ2 ) cn = (f, ϕ2 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . (4.3) (n) (n) (n) A˜ (ϕ1 , ϕn ) c1 + A˜ (ϕ2 , ϕn ) c2 + · · · + A˜ (ϕn , ϕn ) cn = (f, ϕn )
Řešení rovnice (4.1) se tak redukuje na vyřešení soustavy lineárních al(n) (n) gebraických rovnic (4.3). Řešením této soustavy jsou konstanty c1 , . . . , cn , které dosadíme do vztahu (4.2) a tak získáme přibližné řešení un . Algebraickým metodám na řešení soustavy (4.3) se věnovat nebudeme. Metoda konečných prvků se zabývá volbou vhodného prostoru Vn a jeho báze. Snažíme se zajistit, aby matice výše uvedené soustavy byla řídká nebo pásová. Tomuto problému se budeme věnovat ve zbytku této kapitoly. Při výběru konečně dimenzionálního prostoru Vn jsme omezeni pouze dvěmi požadavky: • Vn ⊂ V ∀n = 1, 2, . . . , což v praxi znamená, že bázové funkce ϕ1 , . . . , ϕn musí vyhovovat okrajovým podmínkám • systém {Vn } , n = 1, 2, . . . musí být úplný ve V , tj. musí platit: [
Vn = V.
n=1,2,...
Jestliže je bilineární forma A˜ omezená a V -eliptická, jak tomu v našem případě je, pak výše uvedené podmínky zaručují konvergenci metody konečných prvků, tedy konvergenci přibližného řešení un k přesnému řešení u, jak je popsáno např. v [3].
4.2
Triangulace
Nejprve oblast Ω ⊂ R3 nahradíme polygonální oblastí Ωh , kterou lze rozdělit na množinu čtyřstěnů. Jednotlivé čtyřstěny budeme nazývat elementy a označíme je Ki . Množinu všech elementů označíme Th , tedy Th = {Ki }ni=1 . 30
Množinu Th nazveme triangulace oblasti Ωh , jestliže platí následující vlastnosti: • pro polygonální oblast Ωh platí: Ωh =
n [
Ki ,
i=1
• nastane právě jedna z možností: (i) Ki ∩ Kj = ∅
pro
i 6= j
(ii) jestliže pro i 6= j platí:
Ki ∩ Kj 6= ∅, pak Ki ∩ Kj je t-dimenzionální stěna, t = 0, 1, 2, tedy Ki ∩ Kj je vrchol, celá hrana nebo celá stěna O triangulaci Th dále předpokládáme, že je konzistentní s dělením ∂Ω podle okrajových podmínek, tedy pro ∂Ω = Γ1 ∪ Γ2 ∪ · · · ∪ Γs , požadujeme, aby ∀i = 1, . . . , s byla Γi sjednocením celých dvoudimenzionálních stěn elementů K ∈ Th . def Pro Ki , element triangulaci Th , označme symbolem hi = diam Ki a symbolem ρi maximum průměrů vepsaných koulí do Ki . Dále zaveďme symbol h, který značí maximální průměr elementů triangulace Th , tedy def
h = max hi . i
(4.4)
Máme-li systém triangulací {Th } , h → 0+ , chceme zamezit tomu, aby nastala degenerace elementů, tedy aby pro některá i nastala situace: hi −−−−→ ∞. ρi hi →0+ Proto budeme předpokládat, že existuje konstanta α > 0 nezávislá na h splňující: hi ≤ α, ρi
∀Ki ∈ Th , ∀h → 0+ . 31
(4.5)
Systém triangulací {Th } , h → 0+ splňujících podmínku (4.5) nazveme regulárním dělením oblasti Ωh . Výše popsané vlastnosti triangulace jsou nutné k tomu, aby byla metoda konečných prvků konvergentní, podrobněji viz [3], proto budeme-li dále mluvit o triangulaci oblasti Ωh , budeme mít na mysli triangulaci splňující tyto vlastnosti. Pojmem vrchol triangulace Th myslíme nějaký vrchol nějakého elementu triangulace. Tento vrchol může být (a ve většině případů také je) společným vrcholem více elementů.
4.3
Volba báze
Předpokládejme triangulaci Th a popišme konstrukci báze prostoru Vn . Nejprve tento prostor přeznačme na Vh , jak je zvykem, neboť platí n = n (h), tedy dimenze prostoru závisí na triangulaci Th . Prostor Vh definujeme jako prostor vektorových funkcí, jejichž složky jsou po částech polynomiální spojité funkce: n o 3 def Vh = ϕ ∈ C 0 (Ωh ) : ϕ⌊K ∈ [Pk (K)]3 , ∀K ∈ Th . Každá funkce restringovaná na jednotlivý element triangulace je polynomem stupně k. Podívejme se nyní na ij-tý prvek levé strany soustavy (4.3). A˜ ϕi , ϕj =
Z
a ˜lk (ϕi (x))
∂ϕkj dx. ∂xl
Ω
Integrál přes oblat Ω nejprve aproximujeme integrálem přes polygonální oblast Ωh a pak rozložíme na součet integrálů přes jednotlivé elementy triangulace Th : Z Ω
∂ϕkj dx ≈ a ˜ (ϕi (x)) ∂xl lk
Z
Ωh
X Z ∂ϕkj ∂ϕkj lk a ˜ (ϕ (x)) a ˜ (ϕi (x)) dx = dx. i ∂xl ∂xl K∈T lk
h
K
Snažíme se volit bázové funkce ϕi tak, aby bylo co nejvíce členů sumy integrálů nulových. Dále usilujeme o to, aby bylo co nejvíce prvků matice soustavy (4.3) nulových. Uvedené požadavky splníme volbou bázových funkcí s malými nosiči vzhledem k velikosti Ωh . 32
Bázové funkce budeme volit následujícím způsobem: pro každý vrchol triangulace Th vytvoříme tři bázové funkce. Každá bude v jedné složce nabývat v tomto vrcholu hodnoty jedna a v ostatních dvou složkách bude nulová. V nenulové složce bude na elementech, jejichž vrcholem je i vrchol, na kterém právě pracujeme, polynomem k-tého stupně a to takovým, že mimo tyto elementy bude možné funkci spojitě dodefinovat nulou. Tak dosáhneme toho, že funkce bude nulová na všech elementech, které neobsahují daný vrchol. Vzhledem k tomu, že předpokládáme, že h je velmi malé, tedy platí |K| ≪ |Ωh | , ∀K ∈ Th , dostáváme i |supp ϕi | ≪ |Ωh | pro všechny bázové funkce. Pro lepší představu se podívejme na konkrétní bázové funkce v nejjednodušším případě. Máme jednodimenzionální skalární problém na intervalu (0, 1) a volíme po částech lineární bázové funkce. Oblast Ω = Ωh = (0, 1) rozdělíme rovnoměrně na podintervaly délky h. Dělící body označíme 0 = x0 , x1 , . . . , xn = 1, přičemž platí n = h1 . Tyto body jsou vrcholy triangulace Th , její i-tý element je uzavřený interval [xi−1 , xi ]. Bázové funkce ϕ1 , ϕ2 , . . . , ϕn jsou definovány následovně: 1 − xx1 : x ∈ [x0 , x1 ] def , ϕ1 = 0 : x ∈ (x1 , xn ] def
ϕ2 = def
ϕn =
: x ∈ [x0 , x1 ] : x ∈ (x1 , x2 ] , 0 : x ∈ (x2 , xn ]
x x1 x−x2 x1 −x2
x−xn−1 xn −xn−1
: x ∈ (xn−1 , xn ] . 0 : x ∈ [x0 , xn−1 ]
Prostor Vh je definován jako lineární obal funkcí ϕ1 , ϕ2 , . . . , ϕn . Dále je třeba definovat bázi konzistentní s okrajovými podmínkami. Je-li na některé části hranice ∂Ω definována nulová Dirichletova okrajová podmínka, musíme i na stěnách elementů, které jsou součástí této části hranice, definovat bázové funkce nulové. (O nenulových Dirichletových okrajových podmínkách se záměrně nezmiňujeme, protože se v našem problému nevyskytují.)
33
Definujeme-li v našem pomocném ilustračním příkladě nulovou Dirichletovu okrajovou podmínku v bodě 0. Pak bázi prostoru Vh tvoří pouze funkce ϕ2 , . . . , ϕ n . Podívejme se na konstrukci báze pro náš problém, tedy rovnici (3.7). Silná formulace je třídimenzionální vektorová rovnice, bázové funkce budeme volit po částech lineární v každé složce a jejich předpis odvodíme z barycentrických souřadnic jednotlivých elementů. Předpokládejme triangulaci Th oblasti Ω ⊂ R3 . Nechť Ki je její i-tý element s vrcholy xi1 , xi2 , xi3 , xi4 , pak jsou jeho barycentrické souřadnice dány obecným předpisem: λij (x1 , x2 , x3 ) = aij + bij x1 + cij x2 + dij x3 ,
j = 1, 2, 3, 4.
Z vlastnosti λij (xik ) = δjk (δjk značí Kroneckerovo delta) vyplývá, že koeficienty aij , bij , cij a dij jsou řešením soustavy algebraických rovnic: 1 1 1 1
xi11 xi21 xi31 xi41
xi12 xi22 xi32 xi42
i aj xi13 δj1 i i x23 bj δj2 . = xi33 cij δj3 xi43 δj4 dij
Nyní vezměme libovolný, např. l-tý, vrchol triangulace Th , nechť je tento vrchol zároveň k-tým vrcholem elementu Ki , takže xl = xik . Nechť je tento vrchol průnikem právě ml elementů triangulace Th . Pro jednodušší zápis přepokládejme, že je vrchol xl v každém z těchto elementů právě k-tý. Pak definujeme funkci ψl : i λk (x) : x ∈ Ki , i = 1, . . . , ml def ψl = . S Ki 0 : x ∈ Ωh \ i=1,...,ml
Takto definovaná funkce má ve vrcholu xl hodnotu 1, na elementech, jejichž vrchol je i vrchol xl , lineárně klesá a mimo tyto elementy je spojitě dodefinována nulou. P l VzhledemPk tomu, že zřejmě platí m i=1 |Ki | ≪ |Ωh |, neboť ml je malé, l a supp ψ l = m K , dostáváme také i i=1 |supp ψ l | ≪ |Ωh | . 34
Takto budeme postupovat pro všechny vrcholy triangulace Th , které neleží na hranici Γ1 , na které je definována nulová Dirichletova okrajová podmínka. Označme n počet vrcholů triangulace a p počet vrcholů triangulace ležících na Γ1 . Prostor Vh pak definujeme následovně: * + def
Vh =
(ψ , 0, 0), (0, ψ , 0), (0, 0, ψ ), . . . , (0, 0, ψn−p ) . | 1{z } | {z1 } | {z 1} | {z } ϕ1
ϕ2
ϕ3
(4.6)
ϕ3(n−p)
Při řešení našeho problému budeme pracovat s tímto konečně dimenzionálním prostorem Vh nebo s prostorem, kde jsou bázové funkce polynomy stupně vyššího než 1, který se konstruuje obdobně jako je tomu v lineárním případě.
35
Kapitola 5 Implementace metody konečných prvků v jazyku C++ V této kapitole popíšeme způsob naprogramování metody konečný prvků s využitím logiky objektově oreintovaného programování.
5.1
Základní členění programu
Program je rozčleněn do několika bloků, jejichž základní příkazy jsou volány z metody main, která je v souboru main.cpp. Nejprve stručně popíšeme tuto základní strukturu programu a pak se podrobněji zmíníme o jednotlivých blocích. Program musí dostat jako argument příkazového řádku vstupní soubor, kde je uložena triangulace Th oblasti Ωh ve formátu neutral. V prvním bloku je triangulace ze vstupního souboru načtena, uložena a jsou inicializovány její parametry - barycentrické souřadnice, objemy a těžiště elementů. Ve druhém bloku jsou inicializovány okrajové podmínky na jednotlivých částech hranice ∂Ωh . Ve třetím bloku jsou inicializovány bázové funkce. V dalším bloku je vypočteno přibližné řešení rovnice (3.7). Metodou konečných prvků je naplněna matice soustavy algebraických rovnic (4.3), která je vyřešena s využitím lineárního řešiče umfpack. V posledním bloku jsou vytvořeny výstupní soubory ve formátu vtk a dat, kterých lze použít k vizualizaci výsledků např. v pogramu MayaVi 1.5 36
nebo Tecplot. Celý program je navržen v režimu výjimek, který zaručuje větší ochranu před chybami programátora. V tomto režimu jsou navrženy i všechny datové struktury popsané níže.
5.2
Inicializace triangulace
Všechny datové struktury pro práci s triangulací Th jsou v souboru mesh.cpp. Bod oblasti Ωh je reprezentován instancí třídy point c. Pro tuto třídu jsou definovány operátory ==, =, + a - jako standartní porovnávání, přiřazování, sčítání a odčítání bodů prostoru R3 a operátor * pro násobení vektoru skalárem, resp. vektorový součin dvou bodů (vektorů) - v závislosti na tom, jestli jsou argumenty při volání operátoru * typu double a point c nebo point c a point c. Dále jsou implementovány funkce scalar product, norm a make vector a operátor ≪ pro výpis souřadnic bodů. Všechny tyto funkce výrazně usnadňují práci se třídou point c, ke které dochází při běhu programu velmi často. Vrchol triangulace Th je reprezentován instancí třídy vertex c, která obsahuje jako privátní datovou položku třídu point c, která se inicializuje konstruktorem, dále je možno jí pouze číst konstantní veřejnou metodou třídy vertex c. Všechny vrcholy triangulace jsou uloženy v kontejneru vector. Každý vrchol je v tomto kontejneru uložen pod stejným indexem, pod kterým je uveden ve vstupním souboru. Prostorový element triangulace Th reprezentuje instance třídy element c, která má jako privátní datové položky index každého ze čtyř vrcholů (jak jsou uloženy v kontejneru obsahujícím vrcholy triangulace - tedy i jak jsou zapsány ve vstupním souboru), typ (v naší rovnici (3.7) jsou všechny elementy stejného typu), objem a těžiště elementu a barycentrické souřadnice každého ze čtyř vrcholů. Indexy vrcholů a typ elementu se nastaví konstruktorem, ostatní privátní položky se nastaví metodou třídy mesh c, jak bude popsáno níže. Dále je možno privátní položky pouze číst konstantními veřejnými metodami třídy element c. Všechny elementy triangulace jsou uloženy v dalším kontejneru vector, který zachovává stejná pravidla jako kontejner, ve kterém jsou uloženy vrcholy triangulace. Hraniční element triangulace Th oblasti Ωh je reprezentován instancí třídy boundary element c, která má opět jako privátní datové položky index 37
každého ze tří vrcholů (jak jsou uloženy v kontejneru obsahujícím vrcholy triangulace) a obsah a těžiště elementu. S daty se pracuje shodně jako v případě prostorového elementu. Hraniční elementy patřící do stejné části hranice jsou uloženy stejným způsobem jako prostorové elementy vždy v novém kontejneru vector, pro který je zaveden datový typ boundary elements t. Tento kontejner reprezentuje jednu část hranice se stejným indexem. Všechny části hranice jsou uloženy v dalším kontejneru vector, který reprezentuje celou hranici ∂Ωh , a obsahuje jednotlivé kontejnery obsahující hraniční elementy příslušející konkrétní části hranice. Jednotlivé části hranice jsou přístupné pod stejnými indexy, kterými jsou popsány ve vstupním souboru. Pro třídy vertex c, element c a boundary element c je definován operátor ≪, který u vrcholu triangulace vypíše prostorové souřadnice a u prostorového a hraničního elementu vypíše indexy vrcholů. Pro práci s triangulací Th je vytvořena třída mesh c, o které se předpokládá, že bude právě jednou instanciována (tedy je typu singleton). Jako privátní položky obsahuje kontejnery obsahující vrcholy, prostorové elementy a časti hranice. K jednotlivým položkám těchto tří kontejnerů je umožněn přístup konstantními veřejnými metodami. Datovou strukturu mesh c znázorňuje obrázek 5.1. Třída mesh c obsahuje metodu read, která má jeden z argumentů proměnnou typu input stream, tedy typicky vstupní soubor, ze kterého načte a uloží triangulaci, přičemž jednotlivé vrcholy, elementy a hraniční elementy vytváří jako instance tříd popsaných výše. Dále tato metoda nastaví u všech elementů a hraničních elementů zbývající privátní datové položky (barycentrické souřadnice,. . .), což je umožněno tím, že je ve třídách element c a boundary element c nastavena jako friend. Datovou strukturu mesh c jsme se pokusili navrhnout tak, aby splňovala požadavky programátora obvyklého numerického softwaru. Tedy aby měl programátor inkludováním souboru mesh.h a zavoláním příkazů mesh c mesh; mesh.read(input file); k dispozici triangulaci s informacemi potřebnými k naprogramování metody konečných prvků. Navíc jsou všechny datové položky této třídy privátní, čímž je zajištěna ochrana před chybným přepsáním dat.
38
mesh c
(vertex c, vertex c, vertex c, . . . ) (element c, element c, element c, . . . ) (boundary elements t, boundary elements t, . . . )
boundary element c, boundary element c, ...................
boundary element c, boundary element c, ................... Obrázek 5.1: Struktura třídy mesh c
5.3
Inicializace hraničních podmínek
Datové struktury pro uložení hraničních podmínek jsou implementovány v souboru boundary.cpp. Indexy částí hranice ∂Ωh , na kterých je předepsána nulová Dirichletova okrajová podmínka, jsou uloženy jako privátní položky třídy zero dirichlets c. Přístup k těmto položkám je opět zajištěn konstantními veřejnými metodami. Konstruktor této třídy je funkce s volitelným počtem argumentů. První argument je povinný a obsahuje počet částí hranice, kde je předepsána nulová Dirichletova okrajová podmínka (tedy může být i nula). Další argumenty jsou indexy těchto částí hranice, jak jsou uvedeny ve vstupním souboru. Je-li např. nastavena nulová Dirichletova okrajová podmínka na částech hranice s indexy 1 a 4, je tato hraniční podmíka reprezentována vytvořením instance třídy zero dirichlets c příkazem 39
zero dirichlets c zero dirichlets(2,1,4); Indexy částí hranice ∂Ωh , na kterých je předepsána nenulová Dirichletova okrajová podmínka jsou uloženy totožným způsobem jako privátní položky třídy non zero dirichlets c a reprezentovány vytvořením její instance. Tato třída navíc obsahuje metodu u0, která popisuje Dirichletovu okrajovou podmínku. Obdobným způsobem je v programu reprezentována nenulová Neumannova okrajová podmínka vytvořením instance třídy non zero neumanns c a její metodou g, která představuje funkci popisující Neumannovu okrajovou podmínku, jako je tomu např. ve (2.7). U částí hranice ∂Ωh , na kterých není předepsána ani jedna z výše popsaných okrajových podmínek, implicitně předepisujeme nulovou Neumannovu okrajovou podmínku. Třídy zero dirichlets c, non zero dirichlets c a non zero neumanns c musí být vždy právě jednou instanciovány (tedy jsou opět typu singleton), protože některé z funkcí programu mají jako argumenty konstantní reference na tyto třídy. Proto např. v případě, že se v řešené úloze nevyskytuje nulová Dirichletova okrajová podmínka, musíme vytvořit instanci třídy zero dirichlets c příkazem zero dirichlets c zero dirichlets(0); Tato instance reprezentuje nulovou část hranice oblasti ∂Ωh . Třídy zero dirichlets c, non zero dirichlets c a non zero neumanns c jsou implementovány jako potomci třídy part boundary c, která nese datové položky s informacemi o částech hranice. Tyto položky jsou typu protected. Budeme-li chtít později naprogramovat třídu reprezentující jinou okrajovou podmínku, navrhneme jí jako potomka třídy part boundary c, čímž budeme mít usnadněnu její implementaci.
5.4
Inicializace bázových funkcí
Datové struktury a funkce potřebné k incializaci bázových funkcí jsou v souboru bases.cpp. Nejprve je implementována abstraktní třída abase c, která obsahuje jako čistě virtuální metody vracející hodnotu bázové funkce v daném bodě a hodnotu derivací i-té složky bázové funkce podle j-té proměnné. Dále obsahuje data typu protected, kde je v kontejneru set uložen nosič bázové funkce jako
40
seznam indexů prostorových elementů triangulace, index vrcholu triangulace, ke kterému byla bázová funkce zkonstruována (viz algoritmus v sekci 4.3) a parametr p nesoucí informaci o tom, která složka bázové funkce je nenulová (viz níže). Potomkem této abstraktní třídy je třída baselin c, která reprezentuje lineární bázovou funkci. Její nosič, index vrcholu a parametr p jsou inicializovány konstruktorem. Dále je v souboru bases.cpp implementována globální funkce set bases, která pro konkrétní instance tříd mesh c a zero dirichlets c vytvoří bázové funkce algoritmem popsaným v sekci 4.3. Parametr p udávající, která složka bázové funkce je nenulová, je např. pro funkci ϕ1 v definici (4.6) roven 1, pro funkci ϕ3(n−p) roven 3, atd. Bázové funkce jsou uloženy v kontejneru vector, který obsahuje ukazatele na třídu abase c. Budeme-li potřebovat naprogramovat jiné bázové funkce než lineární, navrhneme je jako potomka třídy abase c. Protože veškerá manipulace programu s bázovými funkcemi probíhá přes ukazatele na třídu abase c (která je reprezentována instancí potomkové třídy, v našem případě třídy baselin c), bude po naprogramování dalšího typu bázových funkcí třeba v kódu programu změnit pouze tvorbu bázových funkcí ve funkci set bases. Jeden z parametrů této funkce (parametr degree) udává stupeň polynomu bázových funkcí. V této verzi programu je přípustná pouze hodnota 1. Tvorba bázových funkcí v závislosti na stupni polynomu je řešena konstrukcí switch(degree). Navrhneme-li tedy třídu reprezentující bázové funkce vyššího stupně, je třeba pouze naprogramovat další větev příkazu switch, kde se budou vytvářet instance této třídy.
5.5
Výpočet přibližného řešení
Pro výpočet přibližného řešení problému je implementována globální fuknce fem lin baselin, která je v souboru algorithm.cpp. Tato funkce uloží přibližné řešení do kontejneru vector. Každý prvek tohoto kontejneru je hodnota přibližného řešení ve vrcholu se shodným indexem jako je index prvku. Jak je patrné z názvu, funkce fem lin baselin počítá přibližné řešení pouze pro lineární bázové funkce. Funkce fem lin baselin má jako argumenty mj. dva ukazatele na funkci. Jako první z těchto argumetnů je třeba dosadit ukazatel na funkci reprezentující levou stranu rovnice, kterou programem řešíme, v našem případě ukazatel na funkci reprezentující formu a ˜. Druhý argument představuje uka41
zatel na funkci popisující objemový integrál na levé straně řešené rovnice. V našem případě by tento integrál byl dán skalárním součinem bázové funkce a objemové síly f , tu jsme však zanedbali, proto dosazujeme ukazatel na funkci konstantně vracející nulu. Povrchový integrál na levé straně řešené rovnice je reprezentován metodou g třídy non zero neumanns c. Konstantní reference na tuto třídu je dalším z argumetnů funkce fem lin baselin. Funkce reprezentující formu a ˜ je zapsána v souboru equation.cpp. Výpočet objemových integrálů pravé strany řešené rovnice není ve funkci fem lin baselin implementován, protože tyto integrály při řešení našeho problému počítat nepotřebujeme. Funkce fem lin baselin počítá pouze přibližné řešení lineární rovnice, ve které se nevyskytují objemové síly, s využitím lineárních bázových funkcí. Hlavička této funkce je navržena co nejobecněji, tedy aby vyhovovala i jiným tvarům řešené rovnice a jiným bázovým funkcím. Proto obsahuje jako argument i ukazatel na funkci popisující objemové integrály pravé strany rovnice, které se v našem problému nevyskytují. Pro implementaci funkcí reprezentujících jiné tvary rovnice lze ještě využít pomocné funkce use umfpack a create solution, které jsou také součástí souboru algorithm.cpp. Funkce use umfpack vyřeší soustavu lineárních algebraických rovnic (4.3) s využitím lineárního řešiče umfpack. Funkce create solution z řešení soustavy (4.3), vypočte hodnoty přibližného řešení uh ve vrcholech triangulace Th a uloží je do kontejneru vector.
5.6
Výhody a nevýhody objektového přístupu
Výhodou objektového přístupu k implementaci metody konečných prvků je logické uspořádání dat do tříd, které kromě samotných dat obsahují i metody pro práci s nimi. Takto je např. zaručeno, že veškeré informace o triangulaci oblasti Ωh jsou uloženy v položkách třídy mesh c. Tyto datové struktury mohou být využívány i samostatně v jiných numerických softwarech inkludováním příslušných hlavičkových souborů. Další výhodou objektového přístupu je částečná ochrana před chybami programátora, která je zajištěna privátností datových položek. Tyto položky se obvykle inicializují konstruktorem třídy, která je obsahuje, a dále je možné je pouze číst konstantními metodami této třídy. Tak je zamezeno nechtěnému přepsání dat, která se v programu inicializují pouze jednou a dále mají být pouze čtena.
42
Nevýhodou je, že program využívající logiku objektového programování je v porovnání s obvyklým přístupem výrazně pomalejší. Jednou z možností jak program zrychlit a zároveň alespoň částečně zachovat jeho objektovost je zveřejnění datových položek jednotlivých tříd. Změníme-li typ datových položek všech tříd, se kterými pracujeme, z private nebo protected na public a přistupujeme-li k nim přímo, nikoli přes metody tříd, zrychlíme program zhruba o 20%. Touto změnou zachováme logické uspořádání dat do tříd, včetně metod pro práci s nimi, ale přijdeme o ochranu dat. Programátor pracující s takto změněnými třídami má možnost kdykoli přepsat jejich datové položky, což může vést k chybám.
43
Kapitola 6 Výsledky numerického řešení V 2. kapitole jsme popsali matematický model posunutí v kosti při jejím zatížení. Nyní uvedeme několik numerických výsledků řešících tento model, které byly dosaženy softwary využívajícími metodu konečných prvků.
6.1
Zatížení kosti
V nejjednodušším případě používáme linearizovaný model (3.7). Kost stojící na pevné podložce zatěžujeme silou 50 000 N ve směru osy z, tedy g = (0, 0, −50 000). Velikost síly volíme s ohledem na dynamické zatěžování kostí, které bývá zvlášť významné v okmžiku doskoku či pádu. Za oblast Ω volíme válec o výšce 40 cm a průměru 2 cm, který na jedné podstavě stojí na pevné podložce (množina Γ1 ) a na druhé je zatěžován silou (množina Γ2 ). Na povrchu válce předepíšeme okrajové podmínky (3.6). Nejprve jsme použili program Netgen, který nahradil oblast Ω polygonální oblastí Ωh a vygeneroval triangulaci Th . Informace o triangulaci jsme použili jako vstupní soubor pro náš program, který numericky vyřešil rovnici (3.5). Výsledky jsme graficky zpracovali v programu Tecplot a jsou zachyceny na obrázeku 6.1. Pro větší názornost výsledků jsme při zobrazování přeškálovali osu z substitucí z ′ = 0.5z. Šipky na povrchu válce znázorňují posunutí u. Použitá triangulace obsahuje 6276 elementů a parametr h definovaný vztahem (4.4) má hodnotu 0.0199.
44
Obrázek 6.1: Posunutí v kosti délky 40 cm, s velikostí Youngova modulu E3 = 18 GPa a zatížené silou 50 000 N. 45
6.2
Zatížení kosti a chrupavky
Jestliže hodnoty Youngova modulu E1 , E2 a E3 snížíme o jeden řád, získá popisovaný materiál vlastnosti chrupavky. Pokusíme se takto simulovat zatížení páteře. Za oblast Ω opět volme válec. Fyziologické rozměry obratle jsou průměr 3 cm a výška 2 cm, meziobratlová ploténka má průměr 3 cm a výšku 2 mm. Oblast Ω je tedy válec o průměru 3 cm a výšce 4.2 cm. Na hranici ∂Ω předepíšeme stejné okrajové podmínky jako v předchozím případě s tím rozdílem, že velikost působící síly stanovíme na 5 000 N. V části válce ve výšce od 2 cm do 2.2 cm nastavíme hodnoty Youngova modulu nižší o jeden řád. Reálná geometrie bohužel nedávám uspokojivé výsledky. Výška chrupavky je v poměru k výšce obratlů příliš malá, takže část výpočetní oblasti, kde jsou menší Youngovy moduly, neobsahuje dostatek elementů triangulace Th , aby na vypočtených výsledcích byla patrná změna vlastností materiálu. Pro dosažení demonstrativních výsledků zvětšíme výšku meziobratlové ploténky pětkrát a celou pracovní oblast dále zvětšíme třicetkrát, adekvátně zvětšíme i působící sílu. Řešíme-li úlohu stejným způsobem jako v předchozím případě, narazíme na problém, neboť na celé oblasti Ω už materiálové konstanty (2.4) nejsou konstantami ale nespojitými funkcemi, takže metoda konečných prvků není konvergentní. Tento problém lze obejít vhodnou volbou bodů nespojitosti těchto materiálových funkcí. Body nespojitosti tvoří dvoudimenzionální stěny příčně přetínající oblast Ω. Jestliže tyto stěny zvolíme tak, aby neprotínaly žádné elementy triangulace Th , ale ležely pouze ve sjednocení stěn elementů, je probléme odstraněn, neboť materiálové funkce jsou pak spojité na každém elementu triangulace. Oblast Ω definujeme jako tři válce položené na sobě (obratel, meziobratlová ploténka, obratel). Triangulace vygenerovaná programem Netgen pro tuto geometrii má požadované vlastnosti s tou vyjímkou, že mezi části hranice ∂Ω jsou zahrnuty i roviny mezi jednotlivými válci. Protože jsme ovšem volili válce těsně přiléhající k sobě, můžeme přebytečné části hranice z triangulace odstranit, aniž bysme tím způsobili, že triangulace přestane vyhovovat podmínkám popsaným v sekci 4.2. Výpočet provedeme obdobným způsobem jako v minulém případě. Použitá triangulace obsahuje 3357 elementů a parametr h má hodnotu 0.353. Výsledek je znázorněn na obrázku 6.2. Pro větší patrnost rozdílu mezi deformovanou a nedeformovanou oblastí 46
Obrázek 6.2: Deformace obratlů a meziobratlové ploténky zatížených silou 5 000 N. Youngovy moduly pružnosti obratlů a chrupavky jsou postupně E3,vert = 18 GPa a E3,car = 1.8 GPa.
47
Obrázek 6.3: Porovnání deformované a nedeformované oblasti, posunutí dvacetkrát zvětšeno. vypočtené posunutí zvětšíme dvacetkrát. Výsledek je zachycen na obrázku 6.3, nedeformovanou oblast znázorňuje příčný řez rovinou osy y. Aby byla lépe patrná absorpce posunutí v meziobratlové ploténce, nastavíme v této části Youngovy moduly nižší o tři řády nikoli pouze o řád, jak jsme to učinili v minulém případě, a celý problém vyřešíme shodným způsobem. Výsledek je zachycen na obrázku 6.4.
6.3
Deformace kosti uchycené v trubce
Pro demonstraci materiálových vlastností kosti proveďme výpočet řešení modelu kosti pevně uchycené v nedeformovatelné trubce. Za oblast Ω opět 48
Obrázek 6.4: Deformace obratlů a extrémě měkké meziobratlové ploténky zatížených silou 5 000 N. Youngovy moduly pružnosti obratlů a chrupavky jsou postupně E3,vert = 18 GPa a E3,car = 0.018 GPa.
49
volme válec, na svrchní podstavu tlačíme silou 1 000 000 N, předepíšeme zde tedy Neumannovu okrajovou podmínku s funkcí g = (0, 0, −1 000 000). Plášť válce je pevně uchycen v trubce, takže na něm definujeme nulovou Dirichletovu okrajovou podmínku. Spodní podstava je volná, proto zde předepíšeme nulovou Neumannovu okrajovou podmínku. Fyziologické rozměry kosti použité v části 6.1 nedávají ilustrativní výsledky. Průměr kosti je v poměru k délce kosti příliš malý, takže je na průřezu kosti málo elementů triangulace Th na to, aby byl charakter deformace dobře patrný. Proto do trubky uchytíme pouze část kosti. Pracovní oblast Ω bude mít průměr 2 cm a výšku 8 cm. Výpočet provedeme obdobným způsobem jako v minulých případech. Použitá triangulace obsahuje 1744 elementů a paramter h má hodnotu 0.0125. Výsledek je zachycen na obrázku 6.5.
6.4
Testování modelu
Za oblast Ω opět volme válec. Ve třetí pětině válce rovněž nastavme materiálové konstanty o jeden řád nižší než ve zbývajích částech. Válec stojící na pevné podložce znovu zatěžujme silou. Jestliže nastavíme sílu a rozměry válce tak, že u chrupavky nastává větší než 15% deformace, projeví se zanedbání nelineráních členů zobrazení a (u), proto proveďme porovnání řešení linearizovaného a nelinearizovaného problému. Pro řešení linearizované rovnice (3.5) s okrajovými podmínkami (3.6) i pro řešení nelinearizovaného modelu (3.1) opět s okrajovými podmínkami (3.6), tedy rovnice Z Z ∂wi li a (u (x)) l dx = g i (x) wi dx, ∀ w ∈ V, ∂x Ω
Γ2
použijeme program Fstrin. Výsledky vizualizujeme v programu Tecplot. Použitá triangulace obsahuje 10142 elementů a parametr h má hodnotu 0.00326. Výsledky linearizovaného problému zachycuje obrázek 6.6, nelinearizovaného problému obrázek 6.7. Rozdíl deformace válce při použití linearizovaného a nelinearizovaného modelu zachycuje obrázek 6.8. Při použití nelinearizovaného modelu je válec méně deformovaný. Červené plošky zobrazují místa, kde jsou patrné rozdíly mezi jednotlivými modely.
50
Obrázek 6.5: Deformace kosti délky 8 cm, s velikostí Youngova modulu E3 = 18 GPa, uchycené v trubce a zatížené silou 1 000 000 N.
51
Obrázek 6.6: Posunutí při zatížení obratlů a meziobratlové ploténky, linearizovaný model. Youngovy moduly pružnosti obratlů a chrupavky jsou postupně E3,vert = 18 GPa a E3,car = 1.8 GPa. 52
Obrázek 6.7: Posunutí při zatížení obratlů a meziobratlové ploténky, nelinearizovaný model. Youngovy moduly pružnosti obratlů a chrupavky jsou postupně E3,vert = 18 GPa a E3,car = 1.8 GPa. 53
Obrázek 6.8: Rozdíl deformace při použití linearizovaného a nelinearizovaného modelu. 54
6.5
Vztah deformace a remodelace kosti
Jak bylo popsáno v úvodní kapitole, remodelační proces zahajují osteoklasty resobováním staré kosti. Aktivita osteoklastů je stimulována retězcem látek RANKL-RANK-OPG, které jsou důležité jako kontrolní mechanismus kostní remodelace. Distribuce těchto látek po objemu kosti probíhá v závislosti na rozložení dlouhodobějšího zatěžování kosti (způsobeného např. chůzí), tedy na rozložení posunutí po objemu kosti. V místech, kde byla resorbována stará kost, začínají osteoblasty budovat kost novou, předtím ale potřebují být aktivovány. Aktivační látky jsou v určitém množství uvolňovány při resorbci staré kosti. Toto množství opět závisí na dlouhodobějším zatížení daného místa. Je-li určité místo v kosti méně zatěžováno probíhá větší resorbce kostní tkáně a dochází k jejímu řídnutí, je-li naopak některé místo více zatěžováno probíhá masivnější tvorba nové tkáně a v daném místě dochází k jejímu zhoustnutí. Toto chování působí uvolňování implantátu kyčelního kloubu ve stehenní kosti. V důsledku zatížení implantátu při chůzi působí zatížení na stehenní kost, která se přetváří. V některých místech dochází k jejímu řídnutí, v jiných naopak k jejímu zhoustnutí. Rozložení posunutí ve stehenní kosti s implantátem kyčelního kloubu bude předmětem našeho dalšího studia.
55
Kapitola 7 Závěr V této práci jsme odvodili model deformace kosti působením vnějších sil a krátce popsali metodu konečných prvků, kterou používáme k numerickému řešení modelu. Dále jsme metodu konečných prvků naprogramovali v jazyku C++ s využitím logiky objektově orientovaného programování a popsali využívané datové struktury. V poslední části práce jsme uvedli několik numerických výsledků dosažených naším programem, které znázorňují deformaci stehenní kosti a deformaci dvou obratlů a meziobratlové ploténky. V kapitole 5 je popsán náš program, který je navržen tak, aby bylo snadné jej rozšířit na řešení složitějších problémů. Cílem naší další práce bude např. naprogramovat nelineární bázové prvky a adaptovat program na řešení nelineárních rovnic s nenulovými objemovými silami nebo implementovat programy umožňující přenos složité geometrie kosti , např. z CT (Computer Tomography) nebo z MRI (Magnetic Resonance Imaging). Dalším objektem našeho studia bude také remodelace kosti a její užší spojitost s dynamickou zátěží (deformací) kosti. Proces remodelace kosti je popsán dodatečnou soustavou obyčejných diferenciálních rovnic, které je nutno řešit v každém elementu kosti - konečném prvku. Neposledním cílem pak bude modelovat deformaci kosti na přechodu kost-implantát, především v případě implantátu kyčelního kloubu.
56
Literatura [1] Bican L.: Lineární algebra a geometrie, ACADEMIA, Praha, 2002. [2] Fučík S., Kufner A.: Nelineární diferenciální rovnice, Nakladatelství technické literatury, Praha, 1978. [3] Haslinger J.: Metoda konečných prvků pro řešení eliptických rovnic a nerovnic, Státní pedagogické nakladatelství, Praha, 1980. [4] Maršík F.: Termodynamika kontinua, Academia, Praha, 1999. [5] Nečas J., Hlaváček I.: Úvod do matematické teorie pružných a pružně plastických těles, Nakladatelství technické literatury, Praha, 1983. [6] Paulsen D. F.: Histologie a buněčná biologie, Nakladatelství H&H Vyšehradská, Praha, 2004. [7] Rektorys K.: Variační metody v inženýrských problémech a v problémech matematické fyziky, Nakladatelství technické literatury, Praha, 1974.
57
