41, s. 283-290, Gliwice 2011
MODEL NUMERYCZNY ODDZIAŁYWAŃ TERMICZNYCH MIĘDZY TKANKĄ SKÓRNĄ I ZEWNĘTRZNYM ŹRÓDŁEM CIEPŁA
BOHDAN MOCHNACKI
MARIUSZ CIESIELSKI
Wydział Inżynierii Mechanicznej i Informatyki, Politechnika Częstochowska e-mail: bohdan.mochnacki@im.pcz.pl, mariusz.ciesielski@icis.pcz.pl
Streszczenie. W pracy rozpatruje się niejednorodny obszar tkanki skórnej (naskórek, skóra właściwa, tkanka podskórna) poddany oddziaływaniu zewnętrznego strumienia ciepła. Przyjęta postać funkcji opisującej przestrzenny rozkład zewnętrznego źródła ciepła determinuje osiowo-symetryczny charakter zadania i zorientowanie rozpatrywanego fragmentu tkanki w układzie współrzędnych walcowych. Procesy cieplne zachodzące w obszarze tkanki opisuje układ cząstkowych równań parabolicznych (tzw. równań Pennesa) uzupełniony odpowiednimi warunkami brzegowo-początkowymi. W obliczenia numerycznych zastosowano schemat jawny metody różnic skończonych dla zadań nieliniowych. W końcowej części pracy przedstawiono wyniki symulacji numerycznych.
1. WSTĘP
Tkanka skórna jest obszarem niejednorodnym (warstwowym), kolejne warstwy odpowiadają podobszarom naskórka, skóry właściwej i tkanki podskórnej – rys. 1.
Rys. 1. Tkanka skórna
Grubości warstw są cechą osobniczą zależną od płci, wieku, a nawet wykonywanego zawodu, w obliczeniach wprowadzono wartości średnie zaczerpnięte z [1].
Procesy cieplne zachodzące w rozważanym obszarze opisuje układ równań różniczkowych cząstkowych w postaci [1, 2, 3]
( ) ( )
T T xt t[ ( )T T( )
xt]
Q ( )
T Q ( )
T
c
x e e e = λe e + pe + me
∂ Ω ∂
∈ , div grad ,
: , (1)
gdzie e = 1, 2, 3 dotyczy kolejnych warstw tkanki, c jest objętościowym ciepłem właściwym, λ jest współczynnikiem przewodzenia ciepła, Qp, Qm są objętościowymi wewnętrznymi źródłami ciepła związanymi z perfuzją i metabolizmem, T, x, t oznaczają temperaturę, współrzędne geometryczne i czas. Składnik perfuzyjny dany jest zależnością
( )
T k( )
T[
T T( )
xt]
Qpe = e b− e , , (2)
gdzie ke(T) = Gbe(T)⋅cb, Gb jest perfuzją wyrażoną w [m3 krwi/s/m3 tkanki], cb jest objętościo- wym ciepłem właściwym krwi, Tb jest temperaturą krwi zasilającej tkankę. Zależność (2) dotyczy podobszaru tkanki zasilanego dużą liczbą włoskowatych naczyń krwionośnych, w przypadku potrzeby uwzględnienia par lub pojedynczych dużych naczyń równanie Pennesa (1) pozostaje w mocy, ale uzupełnia się je równaniami dotyczącymi naczyń i warunkami brzegowymi determinującymi wymianę ciepła między naczyniem i tkanką (por. [3]). Składnik metaboliczny Qme traktuje się z reguły jako stały, przy czym w zależności od warunków fizycznych (wysiłek, odpoczynek) wartość jego może zmieniać się nawet o rząd.
Warunki brzegowe na powierzchniach kontaktu między naskórkiem a skórą właściwą oraz między skórą właściwą a tkanką podskórną są postaci
( ) ( )
( ) ( )
⎪⎩
⎪⎨
⎧
= ∂
λ ∂
−
∂ = λ ∂ Γ −
∈
+ + +
t x T t x T
n t x T n
t x T x
e e
e e e
e e
, ,
, ,
:
1
1
1 , (3)
przy czym wskaźniki e i e+1 identyfikują dwa pozostające w kontakcie podobszary, n oznacza kierunek normalny do brzegu obszaru (w przypadku omawianego zagadnienia n odpowiada kierunkowi osi pionowej z).
Na powierzchni naskórka pozostającej w bezpośrednim kontakcie ze źródłem zewnętrznym zadany jest brzegowy strumień ciepła (warunek Neumanna)
( )
nxt q( )
xtx T , b ,
: 1 1
0 =
∂ λ ∂
− Γ
∈ , (4)
a w szczególności
{ } ( ) [ ]
2
0 0 0
0
, : b , 1 r , 0, , p
r z q x t q r r t t
r
⎡ ⎛ ⎞ ⎤
⎢ ⎥
∈Γ = −⎜ ⎟ ∈ ≤
⎢ ⎝ ⎠ ⎥
⎣ ⎦
, (5)
gdzie r jest współrzędną promieniową w walcowym układzie współrzędnych x = {r, z}, r0 jest promieniem działania źródła, q0 jest maksymalną wartością zewnętrznego strumienia ciepła, tp jest czasem ekspozycji. Dla r > r0 , a dla t > tp na całej górnej powierzchni ciepło oddawane jest do otoczenia zgodnie z prawem Newtona, czyli
( ) [T( )
xt Ta]
n t x
x T =α −
∂ λ ∂
− Γ
∈ , ,
: 1
1 1
0 , (6)
gdzie α, Ta oznaczają współczynnik wymiany ciepła i temperaturę otoczenia. Na powierzchni bocznej umownie przyjętego obszaru walcowego oraz na dolnej powierzchni ograniczającej ten obszar założono warunki adiabatyczne
{ } ( )
n r zn t x z T
r e e = = ∪
∂ λ ∂
− Γ
∈ ∞ , 0,
:
, . (7)
Założono w tym miejscu, że wymiary rozpatrywanego obszaru tkanki są na tyle duże, że na jego peryferiach można przyjąć zerowe strumienie ciepła.
Znany jest również rozkład temperatury w chwili t = 0, czyli
( )
xt T( )
xT
t=0: e , = e0 . (8)
Przedstawiony wyżej układ równań i warunków brzegowo-początkowych stanowi opis matematyczny procesów zachodzących w obszarze tkanki skórnej pokazanym na rys. 2.
Można w tym miejscu zaznaczyć, że działanie bardziej intensywnych źródeł zewnętrznych wymaga pewnej modyfikacji przedstawianego wyżej modelu [4]. Na zakończenie należy podkreślić, że równanie Pennesa dla podobszaru naskórka nie zawiera składników odpowiadających wewnętrznym źródłom ciepła [5].
Rys. 2. Osiowo-symetryczny podobszar tkanki
2. MODEL NUMERYCZNY
Model numeryczny procesów cieplnych w nagrzewanej tkance bazuje na metodzie różnic skończonych w wersji omówionej szczegółowo w [6, 7]. W pierwszej kolejności wprowadza się siatkę czasu
0 1 2 1
... f f f ... F
t < < <t t − <t − < < < < ∞t t (9) o stałym kroku Δt .
Różnicową siatkę geometryczną pokazano na rys. 3. Można zauważyć, że węzły ,,brzegowe’’
umieszczone są w odległości 0.5h lub 0.5k od brzegu (h, k oznaczają kroki regularnej siatki geometrycznej w kierunku r i z). Taka lokalizacja węzłów ,,brzegowych’’ poprawia dokładność rozwiązania numerycznego [6].
Rys. 3. Siatka geometryczna
Operator ∇(λ∇T) dla zadania osiowo symetrycznego w węźle wewnętrznym (i, j) w chwili f-1 (rozpatrujemy jawny schemat różnicowy) jest postaci
( )
1 1
1 ,
, ,
1 f f
f i j
i j i j
T T
T r
r r r z z
− −
− ⎡ ∂⎛ ∂ ⎞⎤ ⎡∂⎛ ∂ ⎞⎤
⎡∇ λ∇ ⎤ =⎢ ⎜ λ ⎟⎥ +⎢ ⎜λ ⎟⎥
⎣ ⎦ ⎣ ∂ ⎝ ∂ ⎠⎦ ⎣∂ ⎝ ∂ ⎠⎦ . (10)
Sposób tworzenia aproksymacji różnicowej operatora (10) dla siatki geometrycznej w otoczeniu punktu (i, j) pokazanego na rys. 4 jest następujący.
Rys. 4. Otoczenie punktu (i, j) siatki geometrycznej
W odległościach 0.5h i 0.5k od punktu (i, j) wprowadza się dodatkowe punkty zaznaczone na rys. 4. Wartości pochodnych w węzłach pomocniczych przybliżamy za pomocą ilorazów różnicowych średnich, czyli
1 1 1 1 1
, 1 , , 1 ,
1
, 0.5 , 0.5 , 1
, 1
, 0.5 2
f f f f f
i j i j i j i j
f
i j i j i j f
i j i j
T T T T
T h
r r r
r h R
− − − − −
+ +
+ −+ −
+ +
− −
⎛ λ∂ ⎞ = λ =⎛ + ⎞
⎜ ∂ ⎟ ⎜⎝ ⎟⎠
⎝ ⎠ , (11)
1 1 1 1 1
, , 1 , , 1
1
, 0.5 , 0.5 , 1
, 0.5 2 1,
f f f f f
i j i j i j i j
f
i j i j i j f
i j
i j
T T T T
T h
r r r
r h R
− − − − −
− −
− −− −
− −
− −
⎛ λ∂ ⎞ = λ =⎛ − ⎞
⎜ ∂ ⎟ ⎜⎝ ⎟⎠
⎝ ⎠ . (12)
W tym miejscu współczynniki przewodzenia λi jf,−+10.5 oraz λi jf,−−10.5 przybliżono średnimi harmo- nicznymi współczynników przewodzenia w tych węzłach, to znaczy
1 1 1 1
, , 1 , , 1
1 1
, 0.5 1 1 , 0.5 1 1
, , 1 , , 1
2 2
,
f f f f
i j i j i j i j
f f
i j f f i j f f
i j i j i j i j
− − − −
+ −
− −
+ − − − − −
+ −
λ λ λ λ
λ = λ =
λ + λ λ + λ . (13)
Wówczas występujące we wzorach (11), (12) wielkości
1 1
, 1 1 1 , 1 1 1
, , 1 , , 1
0.5 0.5 0.5 0.5
f , f
i j f f i j f f
i j i j i j i j
h h h h
R −+ − − R −− − −
+ −
= + = +
λ λ λ λ (14)
odpowiadają oporom cieplnym między węzłem (i, j) i węzłami (i, j+1) oraz (i, j−1). Na podkreślenie zasługuje fakt, że z formalnego punktu widzenia takie same definicje oporów należy przyjąć dla węzłów leżących w pobliżu powierzchni granicznych oddzielających podobszar naskórka od podobszaru skóry właściwej i podobszar skóry właściwej od podobszaru tkanki podskórnej (pamiętając oczywiście, że współczynniki przewodzenia dotyczą dwóch różnych materiałów). Dowód poprawności takiej aproksymacji pochodnych w węzłach ,,brzegowych’’ można znaleźć w monografii [6].
Wykorzystując powtórnie iloraz różnicowy średni, otrzymuje się
1 1 1 1 1
, 1 , , 1 ,
, 1 , 1
, , 1 , 1
,
1 1
2 2
f f f f f
i j i j i j i j
i j f i j f
i j i j i j
i j
T T T T
T h h
r r r
r r r r h R R
− − − − −
+ −
− −
+ −
⎡ − − ⎤
⎡ ∂⎛ λ∂ ⎞⎤ = ⎢⎛ + ⎞ +⎛ − ⎞ ⎥
⎢ ∂ ⎜⎝ ∂ ⎠⎟⎥ ⎢⎜⎝ ⎟⎠ ⎜⎝ ⎟⎠ ⎥
⎣ ⎦ ⎣ ⎦. (15)
Podobnie postępuje się przy konstrukcji wyrażenia różnicowego przybliżającego drugi składnik operatora (10)
1 1 1 1 1
1, , 1, ,
1
0.5, 1
0.5, 1,
f f f f f
i j i j i j i j
f
i j f
i j
i j
T T T T
T
z k R
− − − − −
+ +
+− −
+ +
− −
⎛λ∂ ⎞ = λ =
⎜ ∂ ⎟
⎝ ⎠ , (16)
1 1 1 1 1
, 1, , 1,
1
0.5, 1
0.5, 1,
f f f f f
i j i j i j i j
f
i j f
i j
i j
T T T T
T
z k R
− − − − −
− −
−− −
− −
− −
⎛λ∂ ⎞ = λ =
⎜ ∂ ⎟
⎝ ⎠ . (17)
Następnie
1 1 1 1 1
1, , 1, ,
1 1
1, 1,
,
1
f f f f f
i j i j i j i j
f f
i j i j
i j
T T T T
T
z z k R R
− − − − −
+ −
− −
+ −
⎡ − − ⎤
⎡∂⎛λ∂ ⎞⎤ = ⎢ + ⎥
⎢∂ ⎜⎝ ∂ ⎟⎠⎥ ⎢ ⎥
⎣ ⎦ ⎣ ⎦, (18)
gdzie
1 1
1, 1 1 1, 1 1
, 1, , 1,
0.5 0.5 0.5 0.5
f , f
i j f f i j f f
i j i j i j i j
k k k k
R+− − − R −− − −
+ −
= + = +
λ λ λ λ (19)
są oporami cieplnymi między węzłami (i, j ) i (i+1, j ) oraz (i−1, j).
Ostatecznie
( )
( ) ( )
( ) ( )
1 , 1 1 1 , 1 1 1
, 1 , , 1 ,
1 1
,
, 1 , 1
1, 1 1 1, 1 1
1, , 1, ,
1 1
1, 1,
f i j f f i j f f
i j i j i j i j
f f
i j
i j i j
i j f f i j f f
i j i j i j i j
f f
i j i j
T T T T T
R R
T T T T
R R
− − − − + − −
− +
− −
− +
− − − + − −
− +
− −
− +
Φ Φ
⎡∇ λ∇ ⎤ = − + − +
⎣ ⎦
Φ Φ
− + −
, (20)
gdzie
, ,
, 1 , 1 1, 1,
, ,
0.5 0.5 1
, ,
i j i j
i j i j i j i j
i j i j
r h r h
r h r h k
− + − +
− +
Φ = Φ = Φ = Φ = . (21)
Zgodnie z zasadami dotyczącymi konstrukcji schematu jawnego dla przejścia t f−1 → t f otrzy- muje się równanie różnicowe w postaci
( ) ( )
( ) ( ) ( )
( )1
, , , 1 , 1
1 1 1 1 1
, 1 , 1 , 1 , 1 ,
, 1 , 1
1 1
1, 1 1 1, 1 1
1, , 1, ,
1 1 , ,
1, 1,
f f
i j i j i j i j
f f f f f
i j f i j i j f i j i j
i j i j
f f
i j f f i j f f
i j i j i j i j p m
f f i j i j
i j i j
T T
c T T T T
t R R
T T T T Q Q
R R
− − +
− − − − −
− +
− −
− +
− −
− − − + − −
− +
− −
− +
− Φ Φ
= − + − +
Δ
Φ Φ
− + − + +
, (22)
z którego wyznacza się wartość temperatury w węźle wewnętrznym (i, j) w chwili f.
Założony krok czasu należy tak dobrać, aby zapewnić warunki stabilności rozwiązania numerycznego. Problemy aproksymacji warunków brzegowych (w omawianym zagadnieniu warunków Neumanna i Robina) dla przyjętej siatki geometrycznej omówiono szczegółowo w [6].
3. PRZYKŁADY OBLICZEŃ NUMERYCZNYCH
Omówiony wyżej wariant MRS zastosowano do opracowania programu autorskiego realizującego obliczenia numeryczne związane z nagrzewaniem tkanki. Rozważano podobszar tkanki skórnej w kształcie walca o wymiarach R = 20 mm, Z = 12.1 mm, przy czym grubości poszczególnych warstw tkanki wynosiły: L1 = 0.1 mm, L2 = 2 mm, L3 = 10 mm (indeksy dolne oznaczają kolejno: 1 – naskórek, 2 – skórę właściwą, 3 – tkankę podskórną). Przyjęto następujące uśrednione parametry termofizyczne poszczególnych warstw: λ1 = 0.235 W/(mK), λ2 = 0.445 W/(mK) , λ3 = 0.185 W/(mK) , c1 = 4.3068·106 J/(m3K) , c2 = 3.96·106 J/(m3K), c3 = 2.674·106 J/(m3K) oraz cb = 3.9962·106 J/(m3K), Tb = 37 ºC, Gb1 = 0, Gb2 = Gb3 = 0.00125 (m3 krwi/s/m3 tkanki), Qm1 = 0, Qm2 = Qm3 = 245 W/m3 (warunki odpoczynku). Warunki początkowe wynosiły T1 0 = T2 0 = T3 0 = 37 ºC, temperatura otoczenia Ta = 20 ºC, α = 10 W/(m2K).
Celem symulacji było określenie zależności między wartościami q0 oraz tp, przy których może dojść do oparzenia I lub II stopnia (tj. gdy temperatura powierzchni pomiędzy warstwą naskórka a skórą właściwą osiągnie wartość powyżej 44 ºC). W symulacjach przyjęto stały promień działania źródła ciepła r0 = R/2 = 10 mm.
Na rys. 5 przedstawiono krzywe nagrzewania (dla t ≤ tp) i stygnięcia (dla t > tp) w wybranych punktach analizowanego obszaru o współrzędnych: „0” – (0, 0), „A” – (0, L1),
„B” – (r0/2, L1) oraz „C” – (r0, L1) – punkty A, B i C umieszczono pod warstwą naskórka.
Pokazano wyniki dwóch symulacji, dla których oszacowano: q0 = 5 kW/m2 i tp = 3.726 s oraz q0 = 1.5 kW/m2 i tp = 34.670 s.
Rys. 5. Krzywe nagrzewania i stygnięcia w wybranych punktach analizowanego obszaru (dwa przypadki)
Na rys. 6 zebrano wyniki symulacji numerycznych, które określają zależność q0 od tp. Analizując tę zależność, można np. oszacować dla określonej wartości q0, po jakim czasie ekspozycji tp może nastąpić oparzenie I lub II stopnia. Dla otrzymanej zależności można dobrać funkcję aproksymującą typu potęgowego.
Rys. 6. Zależność q0 w funkcji czasu ekspozycji tp (wyniki numeryczne)
4. WNIOSKI
W pracy dokonano analizy przebiegu procesów cieplnych, które zachodzą w tkance skórnej poddanej oddziaływaniom zewnętrznego źródła ciepła o kształcie paraboloidalnym, które może powodować oparzenia. Rozpatrywano zadanie osiowo-symetryczne. Do symulacji numerycznej zastosowano metodę różnic skończonych, która w stosunkowo prosty sposób może uwzględniać niejednorodność poszczególnych warstw tkanki skórnej. Zaprezentowane w pracy wyniki pozwalają oszacować wpływ natężenia zewnętrznego strumienia ciepła
i czasu ekspozycji na możliwość wystąpienia oparzenia I lub II stopnia. Dalsze etapy prac w tym zakresie dotyczyć będą uzupełnienia algorytmu o procedurę wyznaczania tzw. całki oparzeń [2], której wartość decyduje o stopniu oparzenia.
LITERATURA
1. Torvi D.A., Dale J.D.: A finite element model of skin subjected to a flash fire. “Journal of Biomechanical Engineering” 1994, 116, p. 250-255.
2. Majchrzak E., Jasiński M:: Sensitivity analysis of burns integrals. “Computer Assisted Mechanics and Engineering Sciences” 2004, 11, 2/3, p. 125-136.
3. Majchrzak E., Mochnacki B.: Analysis of bio-heat transfer in the system of blood vessel- biological tissue. KLUWER Academic Publishers, 2001, p. 201-211.
4. Abraham J.P., Sparrow E.M.: A thermal-ablation bioheat model including liquid-to vapour phase change, pressure- and necrosis-dependent perfusion, and moisture- dependent properties. “International Journal of Heat and Mass Transfer” 2007, 50, p.
2537-2544.
5. Majchrzak E., Mochnacki B., Jasiński M.: Numerical modelling of bioheat transfer in multi-layer skin tissue domain subjected to a flash fire. W: Computational Fluid and Solid Mechanics. Elsevier, Vol. II, 2003, p. 1766-1770.
6. Mochnacki B, Suchy J.S.: Numerical methods in computations of foundry processes.
Cracow: PFTA, 1995.
7. Majchrzak E., Mochnacki B.: Metody numeryczne : podstawy teoretyczne, aspekty praktyczne, algorytmy. Gliwice: Wyd. Pol. Śl., 2004.
NUMERICAL MODEL OF THERMAL INTERACTIONS BETWEEN SKIN TISSUE AND EXTERNAL HEAT SOURCE
Summary. Non-homogeneous skin tissue domain subjected to an external heat flux is considered. The form of external heating determines the axially- symmetrical type of task considered and the sub-domain of tissue is oriented in cylindrical co-ordinate system. The thermal processes proceeding in the tissue domain are described by the system of partial differential equations (the Pennes equations) supplemented by the adequate boundary and initial conditions. At the stage of numerical computations the finite difference method for transient and non-linear problems is used. In the final part of the paper the examples of computations and also the conclusions are presented.
Pracę wykonano w ramach realizacji projektu NR 13-0124-10.
