Model matematyczny

Do symulacji numerycznej użyto modelu geometrycznego płytki izolacyjnej składającej się z trzech warstw. Pierwszą warstwę stanowi poszycie kadłuba wykonane z aluminium, drugą polimerowa podkładka usztywniająca SIP, a trzecią materiał izolacyjny LI. Warstwy RCG, oraz RTV w tej fazie obliczeń pominięto. Poszczególne warstwy mają następujące grubości: AL. 2024 - 1,6 mm, NOMEX - 4,394mm, LI900 - 78,486 mm, [34]. Wymiary podstawy płytki wynoszą 76,2 x 76,2 mm. Z płytki został wycięty walec, rys. 2.1, aby w obliczeniach można było rozwiązać zagadnienie dwuwymiarowe, osiowosymetryczne przewodnictwa ciepła Do analizy płytki nieuszkodzonej wykorzystano dwa dodatkowe przypadki różniące się grubością izolacji. W pierwszym z nich grubość wynosiła 98,486 mm, a w drugim 118,486 mm. Ponadto rozpatrzono również przypadek zmiany grubości poszycia i podkładki dla największej średnicy uszkodzenia izolacji (model 5), tabela 2.1. Podział modeli płytki ze względu na rodzaj uszkodzenia przedstawiono w tabeli 2.2.

Tab. 2.1. Dobór grubości warstw dla modelu 5

nr modelu

grubość warstwy

[mm]

zmiana warstwy

5 5

5a 1

poszycie

5b 2

5c 3

5d 3

podkładka

5e 5

5f 6

Rys. 2.1. Sposób wyznaczenia geometrii modelu osiowosymetrycznego płytki izolacyjnej

28 Geometrię uszkodzenia wykonano dla dwóch przypadków. W pierwszym przypadku ubytek materiału izolacji powstaje przez uderzenie w nią małego obiektu poruszającego się z prędkością hiperdźwiękową. W wyniku uderzenia powstaje wyrwa, [34]. Jej kształt i uproszczenie zamieszczono na rysunku 2.2. Przypadek drugi dotyczy możliwości zderzenia z izolacją, która odpadła z innej części pojazdu, rys. 2.3. [12].

Rys. 2.2. Geometria pierwszego rodzaju uszkodzenia: po lewej uszkodzenie rzeczywiste, po prawej uproszczenie

Rys. 2.3. Geometria drugiego rodzaju uszkodzenia: po lewej uszkodzenie rzeczywiste, po prawej uproszczenie

W obszarze płytki rozwiązywane jest równanie przewodnictwa ciepła we współrzędnych cylindrycznych , rys. 2.4:

(

)

(

)

(

)

(

) (2.1)

gdzie:

[kg/m³] – gęstość,

[J/kgK] – ciepło właściwe przy stałym ciśnieniu, [W/m·K] – współczynnik przewodzenia ciepła, [K] – temperatura bezwzględna [K],

z warunkiem brzegowym na powierzchni zewnętrznej osłony postaci:

( ) ( ) ( ) (2.2)

gdzie q(t) jest znaną gęstością strumienia ciepła na powierzchni zewnętrznej osłony, σ oznacza stałą Stefana - Boltzmanna, a  jest emisyjnością osłony. Na pozostałych brzegach płytki przyjęto

D/2

D

D/2

D

D/2

D D/2

D

29 warunek

, oznaczający izolację termiczną tych brzegów. Taki układ warunków brzegowych pozwala rozważać najgorszy wariant przepływu ciepła przez płytkę. Brak wymiany ciepła na izolowanych brzegach płytki powoduje wzrost energii wewnętrznej płytki, a więc jej temperatury. W chwili t = 0 znana jest temperatura w całym obszarze płytki (warunek początkowy) .

Z uwagi na to, że badanie wykonano w oparciu o model dwuwymiarowy, wzór (2.1) przyjmie następującą postać:

( ) ( ) ( ) (2.3) Równanie (2.3) z warunkami brzegowymi oraz warunkiem początkowym, rozwiązywane jest w programie ANSYS metodą elementów skończonych. Po dyskretyzacji obszaru płytki na elementy skończone i przybliżeniu rozwiązania w elemencie za pomocą skończonej kombinacji liniowej funkcji bazowych (w przestrzeni izoparametrycznej – funkcji kształtu), można napisać równanie przewodnictwa cieplnego (2.4) w postaci:

{ } [ ][ ]{ } (2.4)

gdzie:

{ } – wektor gęstości strumienia ciepła [ ] – macierz przewodnictwa cieplnego [ ] – macierz funkcji kształtu

{ } – wektor temperatury w węzłach elementu siatki

Otrzymana temperatura na nagrzanej powierzchni modelu płytki wynika z dopływającego strumienia ciepła do tej powierzchni i promieniowania cieplnego. Dla powierzchni podlegających radiacji (brzeg obszaru płytki, w którym zlokalizowane jest uszkodzenie) zastosowano zależność Siegala oraz Howella, która odnosi się do strat energii i temperatury powierzchni (dla j zmieniającego się od 37 do 51):

∑ ( ) ∑ ( ) (2.5) gdzie:

– liczba radiacyjnych powierzchni,

– delta Kroneckera,

– efektywna emisyjność i – tego elementu powierzchni,

– współczynnik konfiguracji radiacji, – i-ty element powierzchni,

– strumień ciepła tracony przez i-ty element powierzchni, – temperatura bezwzględna i-tego elementu powierzchni

Stan naprężeń dla modelu osiowosymetrycznego w materiałach, z których składa się płytka wynika z prawa Hooke’a (wg rys. 2.4):

{ } [ ] { } (2.6)

gdzie:

{ } ⌊ ⌋ - wektor naprężeń

30 [ ]

[

]

– macierz sprężystości

{ } { } { } – odkształcenie sprężyste;

{ } ⌊ ⌋ – wektor odkształcenia całkowitego

{ } ⌊ ⌋ – wektor odkształcenia wynikającego z obciążeń termicznych

Rys. 2.4. Oznaczenia w globalnym układzie współrzędnych

Używając funkcji kształtu dla dowolnego odkształcenia wykorzystano następującą zależność:

[ ]{ } (2.7)

gdzie:

[ ] – macierz przemieszczeń bazująca na funkcji kształtu { } – wektor przemieszczenia węzła.

Naprężenia Hubera (von Mises’a) wyliczono zgodnie z następującym wzorem:

[ ( ) ( ) ( ) ( )] (2.8) Istotnym ograniczeniem metody elementów skończonych jest to, że powoduje niefizyczne „związanie” ze sobą wszystkich węzłów siatki w danym kroku czasowym. Oznacza to, że mała ilość elementów siatki powoduje istotne błędy w obliczeniach. Efekt ten zmniejsza się przez zwiększenie liczby węzłów pokrywających płytkę, [38]. Symulacja numeryczna dotyczyła określenia nie tylko temperatury w modelu płytki izolacyjnej, ale również naprężeń.

Wymagało to zastosowania takiej siatki elementów skończonych, która umożliwiłaby przeprowadzenie takich obliczeń w sposób prawidłowy. Przy tworzeniu tej siatki wykorzystano typ elementów siatki ”plane223”, określanych w skrócie jako P223. W środowisku programu ANSYS obszary te są powierzchniami (elementami) o ośmiu lub sześciu węzłach. Każdy z węzłów posiada do 4 stopni swobody. Mają zastosowanie w symulacjach modeli osiowosymetrycznych, a taki typ rozpatrywano w pracy. W analizie termomechanicznej na podstawie przyjętego strumienia ciepła i właściwości fizycznych materiałów, program generuje temperatury w każdym elemencie. Na podstawie otrzymanych temperatur obliczane są przemieszczenia oraz naprężenia. Elementy P223 stosuje się m. in. w analizie termicznej w stanie nieustalonym.

𝑧

𝑦 𝑥

𝑟 𝜃

31 Elementy P223 dla modelu osiowosymetrycznego wykorzystują do przybliżenia rozwiązania równania różniczkowego kombinację liniową funkcji kształtu, [29]:

1) dla macierzy sztywności naprężeń i wektora rozszerzalności termicznej:

a) w elemencie siatki czworokątnej:

( ( )( )( ) ( )( )( ) a) w elemencie siatki trójkątnej:

( ) ( ) ( )

( ) ( ) ( ) (2.11)

( ) ( ) ( )

( ) ( ) ( ) (2.12)

2) dla macierzy przewodnictwa cieplnego i wektora pojemności cieplnej:

a) w elemencie siatki czworokątnej:

( ( )( )( ) ( )( )( ) ( )( )( ) ( )( )( ))

( ( )( ) ( )( )

( )( ) ( )( )) (2.13) a) w elemencie siatki trójkątnej:

( ) ( ) ( )

( ) ( ) ( ) (2.14)

Oznaczenia stosowane we wzorach (2.9) – (2.14) są zgodne z rysunkiem 2.5.

Rys. 2.5. Oznaczenia węzłów w elemencie siatki w globalnym układzie współrzędnych Dane dotyczące modeli płytek wraz ze średnicami uszkodzeń przedstawiono w tabeli 2.2.

𝑦 𝑣

32 Tab. 2.2. Numeracja modeli z liczbą elementów i węzłów

nr modelu średnica uszkodzenia [cm] liczba elementów liczba węzłów

1 0 310 1009

2 1 308 999

3 2 328 1071

4 3 321 1064

5 4 361 1186

6 4 373 1226

7 3 337 1106

8 2 318 1043

9 1 297 980

We wszystkich rozważanych przypadkach dla osiowosymetrycznego modelu 2D oś „z”

stanowi oś symetrii. W obliczeniach przyjęto następujące warunki (wg rys. 2.6.):

 gęstość strumienia ciepła q na ściankach S1 oraz S2 jest taka sama (rys. 2.7), wyniki obliczeń wykonane zostały dla dwóch przebiegów gęstości strumienia ciepła, rys. 2.7;

 gęstość strumienia ciepła q = 0 dla ścianek S3, S4, S5;

 emisyjność cieplna ε = 0,85 dla ścianek S1 oraz S2

Rys. 2.6. Numeracja ścianek modelu 𝑆

𝑆

𝑆

𝑆

𝑆₅ 𝑧

𝑥 𝑟

𝑆

𝑆

𝑆

𝑆

𝑆5

𝑧

𝑥 𝑟

33 Rys. 2.7. Gęstość strumienia ciepła w funkcji czasu na brzegu zewnętrznym izolacji (ścianki S1 i S2) [28]

W obliczeniach wytrzymałościowych założono, że płytka posiada górne oraz dolne przytwierdzenie. Poszycie może być rozciągane lub ściskane. W związku z tym zastosowano określony typ podpór rys. 2.8.

Rys. 2.8. Umiejscowienie podpór względem przykładowego modelu: podpora 1 (z lewej), podpora 2 (z prawej)

Statek kosmiczny wykonał lądowanie po upływie 2100 s od momentu wlotu w atmosferę ziemską, [28]. W przypadku pojazdów kosmicznych wyposażonych w taką izolację najwyższe temperatury w poszyciu występują po pewnym upływie czasu od lądowania. W związku z tym do symulacji numerycznej przyjęto czas końcowy równy 8000 sekund.