ÚVOD DO MODELOVÁNÍ V MECHANICE Přednáška č. 9
BIOMECHANICKÉ MODELY ČLOVĚKA Ing. Luděk Hynčík, Ph.D..
OBSAH • Modely na bázi tuhých těles (multi-body) • Jednoduchý 1D model • 2D model v časové oblasti • Komplexní 3D model • Aktivní modely • Deformovatelné modely • Vývoj geometrického modelu • Výpočetní síť • Materiálové parametry • Škálování modelu člověka • Aplikace v dopravě
MOTIVACE Moderní matematické metody umožňují rozsáhlou analýzu technických problémů numerickou cestou Kombinace výpočtu a experimentu šetří náklady a vytváří optimální cestu při vývoji nových výrobků – virtuální prototyping HMI = Human Machine Interaction – nutnost adekvátního biomechanického modelu člověka Různé stupně zjednodušení pro různé aplikace – ne vždy je třeba mít časově náročný model
Vývoj a validace s důrazem na aplikace Metody mechaniky vázaných systémů Metoda konečných prvků Jednoduchý model na bázi tuhých těles Deformovatelný model vybraných částí Aktivní a pasivní model Úzká spolupráce s ESI Group Využití komerčního výpočtového prostředí PAM
VERTIKÁLNÍ 1D MODEL Absolutní a relativní z a = z st z 0 + z r výchylka
Boileau P.E., Rakheja S.: Whole-body vertical biodynamic response characteristics of the seated vehicle driver: Measurement and Model development, International Journal of Industrial Ergonomics 22 (1998) 449-472
Energetický přístup 1 n E k = ∑ mi zɺ i2 2 i =1
1 n E p = ∑ mi zi 2 i =1
Lagrangeovy rovnice
ɺɺ(t ) + Bqɺ (t ) + Kq(t ) = f (t ) Mq f (t ) = k 4 z 0 (t ) + b4 zɺ 0 (t ) Frekvenční oblast;
(K − ω
2
)
M + iωB Q(iω ) = F (iω )
Boileau P.E., Rakheja S.: Whole-body vertical biodynamic response characteristics of the seated vehicle driver: Measurement and Model development, International rychlosti Journal of Industrial Ergonomics 22 (1998) 449-472
Převod do frekvenční oblasti
(
)
Q(iω ) = K − ω 2 M + iωB P(iω ) Impedance = poměr síly a Z (iω ) =
(b4 + iωb4 )(z 0 (ω ) − z 4 (iω )) iωz 0 (ω )
Přenos T (iω ) =
z1 (iω ) z 0 (ω )
Identifikace parametrů m1 = 5.31kg ,
k1 = 310kN / m,
b1 = 400 Ns / M
m2 = 28.49kg ,
k 2 = 183kN / m,
b2 = 4750 Ns / M
m3 = 8.62kg ,
k 3 = 162.8kN / m,
b3 = 4585 Ns / M
m4 = 12.78kg ,
k 4 = 90kN / m,
b4 = 2064 Ns / M
MODEL NA BÁZI TUHÝCH TĚLES Úroveň zjednodušení na korektní popis globální artikulace Z hlediska mechaniky se jedná o vázaný mechanický systém a využívá poznatků o dynamice tuhých těles a kloubových spojení Tuhé těleso = hmotnost, střed hmotnosti a matice setrvačnosti Kloubový prvek = sférický, posuvný, rotační, …
JEDNODUCHÝ 2D MODEL Vývojové prostředí MATLAB Nízká časová náročnost výpočtu
Otevřený řetězec elips Univerzálnost algoritmu
Základní kinematika a dynamika pohybu ve 2D x1 = −a1 sin ϕ1 x2 = x1 − a1 sin ϕ1 − a 2 sin ϕ 2 x3 = x2 − a 2 sin ϕ 2 − a3 sin ϕ 3 ⋮ y1 = a1 cos ϕ1 y 2 = y1 + a1 cos ϕ1 + a 2 cos ϕ 2 y3 = y 2 + a 2 cos ϕ 2 + a3 cos ϕ 3 ⋮
VALIDACE MODELU Pád z výšky a jeho příčiny Ztráta kontaktu Základní kinematika a dynamika lidského těla Podmínky dopadu Experimenty na policejní akademii v Praze
Model těla jako jednoho tělesa Vázaný Mechanický systém
Experimentální data (dobrovolníci)
Srovnání s experimentem
KOMPLEXNÍ 3D MODEL NA BÁZI TUHÝCH TĚLES Komerční software PAM Realistický tvar (povrch vytvořen na základě triangulace, databáze ViewPoint DataLabs) Všechna tuhá tělesa spojena pomocí vazeb (kloubů) do globálního modelu a přidány vzájemné kontakty
SYSTÉM HYPERMESH Grafické prostředí pro „Pre-“ i „Post-Processing“ Geometrické úpravy dat a definice vlastností
VIRTUÁLNÍ VÝVOJOVÉ PROSTŘEDÍ PAM Grafické prostředí zahrnující základní „Pre-Processing“, „Řešič“ i „Post-Processing“ Drobné geometrické úpravy, definice materiálových vlastností, počátečních i okrajových podmínek a vnějšího zatížení Komplexní virtuální prostředí pro virtuální prototyping
MODEL NA BÁZI TUHÝCH TĚLES Reakce modelu na vnější zatížení Síla (náraz), zpomalení, … Definice kontaktních ploch
Kontakty Model kontaktu Časová náročnost
MODEL NA BÁZI TUHÝCH TĚLES – SHRNUTÍ Pasivní chování Kostra a kůže Pouhý povrchový popis „Nulový model materiálu“ 44 tuhých těles Anatomické klouby 13 ohybově-torzních 14 posuvných 18 sférických Fyziologické klouby Klouzání lopatky po hrudním koši – kontakt Kluzné kontakty „Nulové“ proužky kůže pro artikulaci
MODEL NA BÁZI TUHÝCH TĚLES – KOLENO Validace samotného kolenního kloubu
VALIDACE MODELU S HLEDISKA ČELNÍHO NÁRAZU Definované testy zmapované experimenty Zjednodušení externích prvků Sáňová zkouška, počáteční rychlost, definované zpomalení 3-bodový pás a 2-komorový airbag Akcelerometry v hlavě, hrudníku, pánvi, … Head Acceleration Compared to Experiment 60 simulation experiment 50
Srovnání s experimentem
a [g]
40
30
20
10
0 0
20
40
60
80
100 t [ms]
120
140
160
180
200
VALIDACE MODELU S HLEDISKA CHODCE Střet s dodávkou, srovnání s experimentem
VNITŘNÍ ORGÁNY Dynamické vlastnosti vypočtené na základě hustoty Gauss-Ostrogradského věta ∂f dΩ = ∫∫ fni dA ∫∫∫ ∂xi Ω ∂Ω
Kontakty
JEDNODUCHÝ MODEL SVALU
šlacha
Deformovatelný prutový prvek
šlacha
sval
šlacha
sval KP
PP
PP
TP
TP
length L
kontraktilní prvek KP => aktivní síla FKP(L,v,t) = Na(t) Fv(v) FL(L) visko-elastický prvek PP/TP => pasivní síla PP nelineární pružina FPP = k(dL) dL TP lineární tlumič FTP = b v
šlacha
KP PP
TP KP kontraktilní kontraktilní prvek PP pasivní pasivní prvek TP tlumí tlumící prvek
PP
TP
JEDNODUCHÝ MODEL SVALU Náhrada prutovým prvkem – metoda konečných prvků zkrá zkrácený
optimum
nataž natažený
šlacha 1,5 1.5
Fsval/Fmax úroveň roveň aktivace 100% 100%
1.0
pasivn pasivníí část
aktivní aktivní část
sval
50%
FKP+FPP FKP
25%
FPP
0 0
šlacha
Fmax = σAstřední
1
L/L L/Lopt
2
σ = 0.001 GPa
JEDNODUCHÝ MODEL SVALU – AKTIVNÍ MODEL Všechny významné svalové skupiny Každý svalový snopec nahrazen prutovým prvkem
ERGONOMIE – komfort Optimální poloha z hlediska minimální napjatosti 1 2
2 f = ∑ γ i (α i − c ) → min i
Cílová funkce
METODA KONEČNÝCH PRVKŮ (MKP, an. FEM) Moderní, vysoce efektivní numerická metoda pro řešení technických a vědeckých úloh Kontinuum rozdělíme na konečný počet jednotlivých částí (prvků, elementů), čímž vznikne síť (mesh). Posuv libovolného vnitřního bodu prvku vyjádříme pomocí náhradních funkcí (aproximační, tvarové funkce, s kompaktním nosičem) a zobecněných posuvů uzlů prvků
PROSTOROVÁ DISKRETIZACE Ve všech konfiguracích kontinua, resp. u konečných prvků těchto konfigurací, se aplikují stejné interpolační (izoparametrické) funkce pro vyjádření polohového vektoru i posuvu bodu 1D, 2D a 3D prvky úsečka pruty, nosníky
trojúhelníky, čtyřúhelníky
čtyřstěny, dorty, cihly
membrány, skořepiny
Transformační vztahy pro převod funkcí a jejich derivací na izoparametrický prvek Rovnice popisující rovnováhu kontinua v maticové formě
3D REKONSTRUKCE Vychází z řezů (CT, MRI) Hranice jednotlivých orgánů Poloautomatický proces
MKP MODEL Spojení měkkých a tvrdých tkání Vhodné elementy pro MKP Definice materiálových vlastností Definice vzájemných vazeb Fyziologické podmínky
DEFORMOVATELNÝ MODEL
elastická
Tvrdé tkáně – kosti • Tuhá těles • Skořepinová kompakta a 3D spongióza Měkké tkáně • Nelineární viskoelastický materiál • Hyperelastický materiál • Elastický materiál Ligamenty – vázaný kontakt Vzájemné kontakty Předpětí svalů (aktivní model) Tlak v cévách
neelastická pěna
EXPERIMENTÁLÍ MĚŘENÍ Etika Viskoelastické vlastnosti
Ep [kPa]
Es [kPa]
Ŋ [kPa s]
Liver (n=7)
124.1 ± 11.9
3459.9 ± 1566.9
779.4 ± 199.6
Liver
175.9
Spleen (n=2)
5476.1
3029796.9
16969.7
Spleen
119
Critical tension [kPa]
BIOBAG Převod stavové rovnice plynu na kapalinu p = ρRT0
⇒
ρ ∆ρ p = ρ 0 RT0 + ρ 0 RT0 − 1 = p0 + ρ 0 RT0 ρ0 ρ0 ∆ρ ∆ρ p − p0 = ρ 0 RT0 ∆ p = ρ RT ⇒ 0 0
ρ0 Hookeův zákon pro pevné látky
ρ0
E ∆V ∆V = −K p=− 3(1 − 2ν ) V0 V0
Kapaliny ∆V p = p0 − K V0
⇒
∆V p − p0 = − K V0
⇒
∆V ∆p = − K V0
BIOBAG Zákon zachování hmoty
m = ρ 0V = ρV = ( ρ 0 + ∆ρ )(V0 + ∆V )
m = ρ 0V0 + ρ 0 ∆V + ∆ρV0 + ∆ρ∆V ≈ m + ρ 0 ∆V + ∆ρV ∆V ∆ρ =− ⇒ ρ 0 ∆V + ∆ρV0 = 0 ⇒ V0 ρ ∆ρ ∆ρ Plyn ∆p = ρ 0 RT0 Kapalina ∆p = K ρ0 ρ0 K R= ρ 0T0 cp κ= cV
VALIDACE MODELU Standardní testy Nízké zátěže – dobrovolníci Vysoké zátěže – mrtvá těla Etika Experimentální koridory Kroellův test – impaktor dané hmotnosti a rychlosti (energie) naráží do dané oblasti lidského těla, měří se deformace, síla, … Čelní test
VALIDACE MODELU INRETS test – impaktorem je 23.4 kg těžká deska o rozměrech 200 x 100 mm předepsanou rychlostí
S ac rum A cc eleration - High S peed
Contact Force - High S peed
0
15 lower c orridor upper corridor variant 1 variant 2 variant 3 variant 4
-20
lower c orridor upper corridor variant 1 variant 2 variant 3 variant 4 10
F [k N]
a [g]
-40
-60
-80
5
-100
-120
0 0
2
4
6
8
10 t [m s]
12
14
16
18
20
0
2
4
6
8
10 t [m s]
12
14
16
18
20
Contact Force - Low S peed
VALIDACE MODELU
8 lower c orridor upper corridor s im ulation
7
Test podle Viana
Impaktorem je tuhý válec hmotnosti 23.4 kg a poloměru průřezu 150 mm
5 F [k N]
Boční test
6
4
3
2
1
0 0
2
4
6
8
10 t [m s ]
12
14
16
18
20
VALIDACE MODELU Test podle Cavanaugha Čelní test – volant Impaktorem je tuhá tyč o hmotnosti 48 kg a poloměru průřezu 25 mm
0 ms
20 ms
40 ms
60 ms
TĚHOTNÁ ŽENA Zranitelný účastník dopravy 30. týden těhotenství
Analýza účinku stávajících bezpečnostních prvků Vývoj a optimalizace
VALIDACE MODELU Čelní test podle Cavanaugha
Boční test podle INRETS
ŠKÁLOVÁNÍ MODELŮ LIDSKÉHO TĚLA Tuhé nebo deformovatelné Celé modely nebo modely částí Modely pro speciální účely Pasivní a aktivní modely Modely pro legislativu
MOTIVACE 5%-ní, 50%-ní a 95%-ní velikosti Každý jsme jiný Biologická rozmanitost lidského těla Pohlaví, rasa, geografické podmínky, posilování atd. Škálování a morfování referenčního modelu • Škálování = globální rozměry a materiálové vlastnosti apod. • Morfování = detaily (např. rysy nebo patologické změny apod.
REŠERŠE Populační skupiny Vliv chování člověka a přírodní vlivy Izolace a adaptace k určitému území Ovlivní i vývoj populace do budoucna Míry lidské rozmanitosti Statistika, interpretace ”průměrného člověka” Výška, hmotnost, forma
ANTROPOMETRIE Výška • Normální mezi 1.4 m a 2 m • Výjimky existují • Vyšší lidé dále od rovníku
Hmotnost • Normální mezi 30 kg a 100 kg
Etnické rozdíly • Ženy kmene Khoi a San nebo !Kung v jižní Africe • Těžší lidé v chladnějším klimatu
Forma • Poměr stojící a sedící výšky (cormic index) • Poměr délky trupu k délce těla • Kratší tělo + delší nohy < 50
DEFINICE ROZMĚRŮ
BODY MASS INDEX
body mass BMI = body height ⋅ body height
SHRNUTÍ Vybrat množinu antropometrických parametrů Definovat cílovou skupinu lidí Najít vhodná antropometrická data jako funkci nějakého parametru Vyvinout algoritmus škálování antropometrických dat jako funkci parametru
CÍLOVÁ SKUPINA Československá spartakiáda 1985 5000 mužů a žen Jedinečná množina lidí od 6 do 55 let Interpolace Parametry – věk a pohlaví Děti – nutné morfování
REFGERENČNÍ MODELY Na základě spolupráce s Rodina modelů ROBBY ROBBY2 50% male HARB model ROBINA 5% female HARB model Výhody a nevýhody Segmentace těla Deformovatelné modely
ROBBY2 Hmotnost = 76.59 kg Stojící výška = 178.72 cm Sedící výška 98.44 cm BMI = 24
ROBINA Hmotnost = 48.22 kg Stojící výška = 141.71 cm Sedící výška 76.75 cm BMI = 24
BOBBY Hmotnost = 24 kg Stojící výška = 112.87 cm Sedící výška 64.65 cm BMI = 19
SEGMENTACE Hlava a krk Hrudník Abdomen Pánev Paže Předloktí Dlaň Stehno Lýtko Chodidlo
ALGORITMUS 1. Škálování všech segmentů v globálním souřadnicovém systému 2. Přepočet hmotnosti a matice setrvačnosti 3. Spojení modelu 4. Pozice modelu
EVOLUCE Algoritmus škálování muže a ženy Československé populace z roku 1985 v rozmezí od 5 do 55 let Software pro škálování základního modelu
VALIDACE Univerzita v Heidelbergu 13-ti letý chlapec Sáňová zkouška 3-bodový pás a airbag Počáteční rychlost = 41 km/h Předepsané zpomalení Akcelerometry a silové senzory
ANALÝZA VÝSLEDKŮ
POPIS PORANĚNÍ Kritérium poranění aorty (síla, tlak, prodloužení) Studium poranění aorty během kritického přistání nebo „pádu“ letadla (vrtulníku) Využití normovaného ATD místo modelu člověka Spolupráce s AMSAFE (USA) Reálné případy
založené na databázích ATD sedící v kokpitu letadla
PODPORA VÝVOJE SEDAČKY KLUZÁKU Kritické přistání kluzáku – pravidla FAR „Nosnost sedačky“ během nárazu
OCHRANA CHODCŮ Optimalizace kapoty automobilu, pasivní a aktivní prvky Korektní biomechanický model kolene
Čas = 0 ms
50 ms
100 ms
250 ms
OCHRANA CHODCE VŮČI TRAMVAJI CDV identifikovalo 3 scénáře nehody Chodec přecházející (přebíhající) před tramvají Aktivní chodec Analýza rizika poranění při primárním nárazu
POPIS NEHODY AUTOMOBILU Přímý popis poranění jako demonstrátor
OCHRANA ŘIDIČE A CESTUJÍCÍCH V TRAMVAJI Osazení modelu kolejového vozidla řidičem a pasažéry Analýza vlivu nárazníku (geometrie, materiál) Vyhodnocení zrychlení na sedadle Vyhodnocení kritérií poranění Návrhy na zlepšení bezpečnosti
DĚKUJI VÁM ZA POZORNOST
[email protected]
