NAMÁHÁNÍ TEPENNÉ STĚNY: LIDSKÁ BŘIŠNÍ AORTA
Kurs: Biomechanika II Obor: Biomechanika a lékařské přístroje Program: Magisterský Fakulta strojní ČVUT v Praze
Lukáš Horný
[email protected]
CÍLE Analytickými metodami získat kvantitativní odhad napjatosti tepenné stěny Tenkostěnná vs. silnostěnná aproximace Elastostatika (nelineární materiál) Konečné deformace Zbytková napětí Jako příklad poslouží lidská břišní aorta
MOTIVACE Numerické metody budou v klinicky podstatných úlohách, optimalizace implantace stentu (z hlediska přetížení) napojení bypassu (opět přetížení) napětí na rozhraní kalcifikovaný plát-stěna interakce tepenné stěny s krví (zejména v místech plátů) disekce stěny a ruptura aneurysmatu atd., vždy hrát prim. Analytické metody neslouží (jak by se v „rakousko-uherské tradici“ výkladu matematiky a fyziky bohužel mohlo někdy zdát) k memorování výrazů, ale k pochopení kvalitativních vlastností řešení, na nichž je testována hodnověrnost výsledků numerických metod.
AD REM
Břišní aorta – kde to je, co to je?
Repro: http://www.doereport.com/ enlargeexhibit.php?ID=15311
Repro: http://my.clevelandclinic.org/heart/heart-blood-vessels/aorta.aspx
AD REM
Břišní aorta – kde to je, co to je? 1 Right lung 2 Right hepatic vein 3 Liver 4 Left hepatic vein 5 Stomach 6 Left colic flexure (splenic flexure of the colon)
7 Spleen 8 Left lung 9 Aorta
Repro: http://www.info-radiologie.ch/en/abdominal_ct.php#
AD REM
Je to elastická tepna
Repro: http://php.med.unsw.edu.au/embryology/ images/a/ae/Artery_histology_16.jpg
Repro: http://php.med.unsw.edu.au/embryology/ index.php?title=File:Blood_vessel_wall_cartoon.jpg
Repro: http://www.lab.anhb.uwa.edu.au/mb140/ corepages/vascular/vascular.htm
MECHANICKÉ INTERAKCE
Pasivní interakce
Přenos tepla a působení vnějšího prostředí
MECHANICKÉ INTERAKCE
Pasivní interakce - externí
Přenos tepla a působení vnějšího prostředí
IN SITU
EXCIZE
EX SITU
MECHANICKÉ INTERAKCE
Pasivní interakce - interní
Po rozříznutí prstýnku se tepna dík uvolnění zbytkových napětí rozevře
MECHANICKÉ INTERAKCE
(nafukování trubice)
Inflační test
Aktivní interakce Maximální kontrakce po podání norepinefrinu céva ztuhne
Bazální tonus
Maximální dilatace Tahová zkouška
po podání papaverinu
Repro: P. Fridez, A. Makino, D. Kakoi, H. Miyazaki, J.-J. Meister, K. Hayashi and N. Stergiopulos 2002 Adaptation of Conduit Artery Vascular Smooth Muscle Tone to Induced Hypertension Annals of Biomedical Engineering 30:7, 905 – 916. http://www.springerlink.com/content/v257812562p17374/
VÝPOČTOVÝ MODEL
Předpoklady pro formulaci úlohy Konstitutivní rovnice Předpoklady o geometrii Předpoklady o zatížení a vazbách Předpoklady o deformaci (konkrétní kinematika)
Repro: http://my.clevelandclinic.org/heart/heart-blood-vessels/aorta.aspx
VÝPOČTOVÝ MODEL
Konstitutivní rovnice Použijeme model Guccione–McCulloch–Waldman, který popisuje cylindricky ortotropní hyperelastický materiál. 2 2 2 c1 c2 EΘΘ + c3 ( EZZ + ERR ) − 1 W = e 2
J. Guccione, A. McCulloch, L. Waldman (1991) Passive material properties of intact ventricular myocardium determined from a cylindrical model. J Biomech Eng 113:42–55 http://dx.doi.org/10.1115/1.2894084
Konkrétní vyčíslení materiálových parametrů pro lidskou břišní aortu převezmeme od kolektivu M.R. Laborsseho, který publikoval výsledky 16 inflačně-extenzních testů... http://www.sciencedirect.com/science/article/pii/S175161611200210X?v=s5 Michel R. Labrosse, Eleanor R. Gerson, John P. Veinot and Carsten J. Beller (2012) Mechanical characterization of human aortas from pressurization testing and a paradigm shift for circumferential residual stress by Journal of the Mechanical Behavior of Biomedical Materials, in press. http://dx.doi.org/10.1016/j.jmbbm.2012.08.004
VÝPOČTOVÝ MODEL
Cylindrická ortotropie ...tři na sebe navzájem kolmé osy materiálové symetrie, které jsou totožné s osami válcového souřadnicového systému
Pasivní odezva materiálu bez aktivace hladkých svalových buněk
Isochorický děj (nestlačitelnost = během deformace materiál nemění objem)
VÝPOČTOVÝ MODEL σ =F
∂W ( F )
c1 W = e 2
∂F
2 c2 EΘΘ + c3
(
Formální konstitutivní rovnice F je deformační gradient W je hustota deformační energie I je jednotkový tenzor 2. řádu p je hydrostatická složka napětí vzniklá reakcí na omezení stlačitelnosti
− pI
2 2 EZZ + ERR
) − 1
zde je použit tenzor E (GreenLagrangeův tenzor deformace), který převedeme na F:
E=
1 2
(F
T
F− I
)
VÝPOČTOVÝ MODEL
Geometrie budeme předpokládat, že břišní aorta je válcová trubice konstantního poloměru R a tloušťky H (po vyjmutí z těla). Tyto údaje převezme opět z literatury: Michel R. Labrosse, Eleanor R. Gerson, John P. Veinot and Carsten J. Beller (2012) Mechanical characterization of human aortas from pressurization testing and a paradigm shift for circumferential residual stress by Journal of the Mechanical Behavior of Biomedical Materials, in press. http://dx.doi.org/10.1016/j.jmbbm.2012.08.004
VÝPOČTOVÝ MODEL
Deformace budeme předpokládat, že během tlakování trubice: • zůstávají všechny řezy rovinné (resp. válcové) • řezy se mohou vzdalovat/přibližovat • řezy se vůči sobě nenatáčejí
Pozor: Jde o pouhou skicu, obrázek nedodržuje konstantní objem deformovaného elementu...
VÝPOČTOVÝ MODEL
Deformace
modelově vylučujeme existenci zkosů...
γ ΘZ ≠ 0 γ RZ ≠ 0
γ RΘ ≠ 0
krut, neboli odklopení od podélné osy (reálně může nastat když (1) je helikální proudění krve, (2) nesymetrie vazeb...)
sklopení ve směru obvodu (přechází-li při nafukování válec ve válec, nemůže nastat)
sklopení do směru podélné osy (reálně musí nastávat jako reakce na tření krve)
VÝPOČTOVÝ MODEL
Deformace budeme si tedy přestavovat, že aorta (válcová trubice) se jen nafukuje a protahuje z tvaru válce do tvaru válce Materiálový bod Q v referenčních a průběžných válcových souřadnicích Q ( r ,θ , z ) Q ( R , Θ, Z )
r
R
Z
z
Q≡Q
⇒ r = r ( R) θ = Θ z = z ( Z)
TENKOSTĚNNÁ SIMULACE Vše (síly, deformace,...) se odehrává jen na úrovni střední plochy σ zz
P
σ θθ
h = λrR H r = λθΘ R z = λzZ Z F=
h
∂x ∂X
F=
l
P
λrR F = 0 0
2r Fred
Fred
σ zz P
Repro: http://nptel.iitm.ac.in/courses/Webcourse-contents/IIT-ROORKEE/strength%20of%20materials/lects%20&%20picts/image/lect15/lecture15.htm
0
λθΘ 0
h 0 H 0 = 0 λzZ 0
0 r R 0
∂h ∂H 0 0
0 0 l L
0 ∂r ∂R 0
0 0 ∂z ∂Z
TENKOSTĚNNÁ SIMULACE Nestlačitelnost (isochorický děj)
dv J= = det F = 1 ⇔ dV
λrR =
λrR det 0 0
1
λθΘλzZ
0
λθΘ 0
0 0 = λrR λθΘ λzZ = 1 λzZ
TENKOSTĚNNÁ SIMULACE Intenzita vnitřních sil (napětí) odpovídá průměrné hodnotě po tloušťce objektu, která je rozprostřena uniformě (membrána) σ zz
P
σ θθ h
l
P 2r Fred
Fred
P σ rr = − 2 rP σ θθ = h Fred rP σ zz = + 2π rh 2h
σ zz P
Repro: http://nptel.iitm.ac.in/courses/Webcourse-contents/IIT-ROORKEE/strength%20of%20materials/lects%20&%20picts/image/lect15/lecture15.htm
TENKOSTĚNNÁ SIMULACE
Konstitutivní rovnice 2 2 c1 c2 EΘΘ + c3 ( EZZ + E2RR ) W = e − 1 2
σ=
∂W ( F ) ∂F
→
E=
FT − p I
1 2
(
FT F− I
)
→
Eij =
(
)
1 2 λ −1 2 ij
∂W σ rr = λrR −p ∂λrR ∂W σ θθ = λθΘ −p ∂λθΘ ∂W σ zz = λzZ −p ∂λzZ
TENKOSTĚNNÁ SIMULACE
Finální soustava elastostatických rovnic ∂W P ∑ Fr : λrR ∂λ − p = − 2 rR ∂W rP ∑ Fθ : λθΘ ∂λ − p = h θΘ Fred ∂W rP ∑ Fz : λzZ ∂λ − p = 2π rh + 2h zZ
TENKOSTĚNNÁ SIMULACE Tuto soustavu vyřešíme ve třech krocích:
(1) Určíme „neurčitou“ reakci na nestlačitelnost p
P ∂W p = − λrR 2 ∂λrR ...a potom ve všech rovnicích substituujeme λrR = 1/(λθΘ· λzZ)
TENKOSTĚNNÁ SIMULACE (2) Vypočteme sílu Fred nutnou k dosažení zvoleného počátečního předepnutí λzZini (ta bude dále konstantní, což odpovídá experimentu, kdy nafukujeme svislou trubku s konstantním axiálním přívažkem) (2a) Pro zvolené λzZini určíme λθΘini za podmínky P = 0 z rovnice:
∂W rP λθΘ −p= ∂λθΘ h
(2b) Pro zvolené λzZini určíme Fred z rovnice:
a vypočtené λθΘ za podmínky P = 0
Fred ∂W rP λzZ −p= + ∂λzZ 2π rh 2 h
TENKOSTĚNNÁ SIMULACE (3) Simulujeme odezvu materiálu na zvolený vnitřní tlak P a předepínací sílu Fred. Čili řešíme soustavu rovnic níže tak, že dosazujeme za P (např. od 0 do 16 kPa) a vypočítáváme λθΘ a λzZ (s tím, že na počátku pro P = 0 se výsledky samozřejmě kryjí s předem vypočtenými λθΘini a λzZini).
∂W rP λθΘ −p= ∂λθΘ h Fred ∂W rP λzZ −p= + ∂λzZ 2π rh 2 h
TENKOSTĚNNÁ SIMULACE Pro další výklad si stáhněte soubor artery-thinwalled-tube.mw z www.biomechanika.cz
TENKOSTĚNNÁ SIMULACE
Ukázka výsledků simulace inflačně-extenzního testu lidských infrarenálních aort z Labrosse et al. (2012)
TENKOSTĚNNÁ SIMULACE
TENKOSTĚNNÁ SIMULACE
Nedostatky modelu σ ij ( x ) = ? x ∈ ( 0 , h ) P σ rr ( r ) = − 2
SILNOSTĚNNÁ SIMULACE
Silnostěnná nádoba integrace po tloušťce stěny...
dσ rr ( r ) σ rr ( r ) − σ θθ ( r ) + =0 dr r ro
Fred + π ri2 P − 2π ∫ σ zz ( r ) rdr = 0 ri
σ rr ( ri ) = − P σ rr ( ro ) = 0
SILNOSTĚNNÁ SIMULACE σ rr ( ri ) = − P
dσ rr σ rr − σ θθ + =0 dr r σ rr ( ro ) = 0 dσ rr σ rr − σ θθ =− dr r dσ rr σ θθ − σ rr = dr r σ θθ − σ rr dσ rr = dr r σ rr ( ro ) ro σ θθ − σ rr dσ rr = σ rr ( ro ) − σ rr ( r ) = 0 − σ rr ( r ) = ∫ dr ∫ r r σ rr ( r ) ro
Integrál jako funkce dolní meze
σ rr ( r ) = − ∫ r
σ θθ − σ rr r
dr
SILNOSTĚNNÁ SIMULACE Dále uvažme jednotkovou krychli, k níž přiložíme napětí σrr, σθθ a σzz, která ji zdeformují na kvádr o hranách λrR, λθΘ a λzZ. Přírůstek deformační energie dW při další diferenciální deformaci o dλkK (tj. např. λrR přejde na λrR +dλrR je: dW = λzZ λrRσ θθ dλθΘ + λzZλθΘσ rr dλrR + λrRλθΘσ zzdλzZ
σzz
Pro nestlačitelný materiál diferencováním: λrRλθΘλzZ = 1 ⇒ dλrRλθΘλzZ + λrRdλθΘλzZ + λrRλθΘdλzZ = 0
λzZ ...dosazením: dW = λrRλzZ (σ θθ − σ rr ) dλθΘ + λrRλθΘ (σ zz − σ rr ) dλzZ σrr když předpokládáme, že nestlačitelnost eliminuje jednu proměnnou...: dW = ∂W dλθΘ + ∂W dλzZ ∂λθΘ
Na závěr porovnáme ∂W σ θθ − σ rr = λθΘ výrazy pro dW: ∂λθΘ
ER e r
∂λzZ
σ zz − σ rr = λzZ
∂W ∂λzZ
σθθ
λθΘ
EZ e z
EΘ e θ
λrR
SILNOSTĚNNÁ SIMULACE To, co jsme udělali, je eliminace parametru p z konstitutivních rovnic...
∂W σ rr = λrR −p ∂λrR ∂W σ θθ = λθΘ −p ∂λθΘ ∂W σ zz = λzZ −p ∂λzZ
∂W σ θθ − σ rr = λθΘ ∂λθΘ ∂W σ zz − σ rr = λzZ ∂λzZ
SILNOSTĚNNÁ SIMULACE
Což dosadíme... σ θθ
ro
∂W − σ rr = λθΘ ∂λθΘ
σ rr ( r ) = − ∫ r
σ θθ − σ rr r
dr
ro
∂W dr σ rr ( r ) = − ∫ λθΘ ∂λθΘ r r Čili vyjádřili jsme radiální napětí na základě kinematiky a konstitutivní rovnice
SILNOSTĚNNÁ SIMULACE ∂W Fred = −π ri P + 2π ∫ σ zz rdr = −π ri P + 2π ∫ σ rr + λzZ ∂λzZ ri ri ro
ro
2
σ zz − σ rr = λzZ
2
∂W ∂λzZ
ro
⇒ σ zz = σ rr + λzZ
ro
rdr
∂W ∂λzZ
ro
( )
d r2
∂W Fred = −π ri P + π ∫ σ rr 2rdr + 2π ∫ λzZ rdr = −π ri2 P + π ∫ σ rr dr + ∂λzZ dr ri ri ri 2
ro
∂W +2π ∫ λzZ rdr ∂λzZ ri
SILNOSTĚNNÁ SIMULACE ro
Fred = −π ri P + π ∫ σ rr 2
ri
( ) dr +2π
d r2 dr
ro
∂W ∫r λzZ ∂λzZ rdr i
o −π ri2 P = π ri2σ rr ( ri ) − π ro2σ rr ( ro ) = −π r 2σ rr ( r )
r
ri
ro
( )
d r2
ro
∂W 2 Fred = −π r σ rr ( r ) + π ∫ σ rr dr + 2π ∫ λzZ rdr = ri dr ∂λzZ ri ri ro
o o o dσ rr 2 σ θθ − σ rr 2 ∂W ∂W = −π ∫ r dr + 2π ∫ λzZ rdr = −π ∫ r dr + 2π ∫ λzZ rdr dr ∂λzZ r ∂λzZ ri ri ri ri
ro
Integrace per partes
r
r
dσ rr σ rr − σ θθ + =0 dr r
r
SILNOSTĚNNÁ SIMULACE ro
Fred = −π ∫
σ θθ − σ rr
ri
r
ro
∂W r dr + 2π ∫ λzZ rdr ∂λzZ ri 2
σ θθ − σ rr = λθΘ ro
Fred = −π ∫ λθΘ ri
∂W ∂λθΘ
ro
∂W ∂W rdr + 2π ∫ λzZ rdr ∂λθΘ ∂λzZ ri ∂W ∂W Fred = π ∫ 2λzZ − λθΘ rdr ∂λzZ ∂λθΘ ri ro
SILNOSTĚNNÁ SIMULACE
Finální model σ rr ( ri ) = − P σ rr ( ro ) = 0
ro
∂W dr σ rr ( r ) = − ∫ λθΘ ∂λθΘ r r
W = W ( λθΘ , λzZ )
∂W σ θθ = λθΘ + σ rr ∂λθΘ ∂W σ zz = λzZ + σ rr ∂λzZ
∂W ∂W Fred = π ∫ 2λzZ − λθΘ rdr ∂λzZ ∂λθΘ ri ro
λzZ ≠ λzZ ( z )
SILNOSTĚNNÁ SIMULACE
Jak zahrnout zbytková napětí Rozevírání prstýnku můžeme modelovat jako otevírání uzavřeného kruhového prutu... Úhel rozevření
α
SILNOSTĚNNÁ SIMULACE Mo
2α
Mo
Rozstřiženou konfiguraci budeme modelovat jako mezikruhovou výseč
ρi ρ ρo Ri R
Ro
α
SILNOSTĚNNÁ SIMULACE Rozstřiženou konfiguraci budeme považovat za beznapěťovou, a tak za referenční Uvažujeme-li tepnu, jako trubici, pak budeme zavírání modelovat jako přechod válcové výseče 2π–2α do válce za podmínky zachování objemu
SILNOSTĚNNÁ SIMULACE
Kinematika – 3 konfigurace ∂R ( ρ ) ∂ρ F1 = 0 0
0
πR (π − α ) ρ 0
0 0 δ
∂r ( R ) ∂R F2 = 0 0
ro
0 r 0 R 0 λ 0
r ri
Otevřená = referenční Uzavřená = zbytková napětí Uzavřená = zbytková napětí Natlakovaná-napnutá
(ρ
ψ
ξ ) → ( R Θ Z)
(R
R = R( ρ )
ξ ∈ 0 , Ξ
Θ=
Z ∈ 0 , L
πψ π −α Z = δξ
δ=
L Ξ
z ∈ 0 , l
Θ Z ) → (r θ
z)
r = r ( R)
λ=
l L
θ =Θ z = λZ
SILNOSTĚNNÁ SIMULACE
Kinematika ∂r ( R ) ∂R F = F2 F1 = 0 0
∂R ( ρ ) 0 0 ∂ρ r 0 0 R 0 λ 0
0
πR (π − α ) ρ 0
∂r 0 ∂ρ 0 = 0 δ 0
0
πr (π − α ) ρ 0
0 λr ρ 0 = 0 0 λδ
0
λθψ 0
0 0 λzξ
SILNOSTĚNNÁ SIMULACE
Kinematika Pro integraci je třeba vyjádřit všechny funkce jakožto závislé na zdeformovaném (finálním) poloměru r, což provedeme z podmínky zachování objemu... λ
(
)
(
π l ro2 − r 2 = π L Ro2 − R2
)
2π − 2α 2 ρ o − ρ 2 Ξ = π L Ro2 − R2 2
(
)
(
)
2π − 2α 2 ρ o − ρ 2 Ξ = π l ro2 − r 2 2
(
)
(
→ R = Ro2 −
→ ρ = ρ o2 −
)
(
l 2 2 ro − r L
π
) δ
(
L 2 Ro − R2 π −α Ξ
→ ρ = ρ o2 −
π π −α
(
)
λzξ ro2 − r 2
λzξ =
)
l = δλ Ξ
SILNOSTĚNNÁ SIMULACE Pro další výklad si stáhněte soubor artery-thickwalled-tube.mw z www.biomechanika.cz
SILNOSTĚNNÁ SIMULACE
Zbytková deformace (parametrizovaná úhlem rozevření α) má jen malý vliv na tvar křivky P – ro a Fred. Fred (α = 0° , λz = 1.15 , P = 16 kPa ) = 1.95 N Fred (α = 80° , λz = 1.15 , P = 16 kPa ) = 1.99 N
SILNOSTĚNNÁ SIMULACE Zbytková deformace homogenizuje průběh napětí
