Rozwiązywanie równań różniczkowych cząstkowych wymaga znajomości zaawansowanego aparatu matematycznego. Dla prostych przypadków jesteśmy w stanie uzyskać rozwiązanie analityczne, dla przypadków złożonych najczęściej jest to niemożliwe. Jedną z metod numerycznego rozwiązywania równań cząstkowych jest metoda objętości skończonych, która polega na podziale analizowanego obszaru na małe elementy kontrolne i przeprowadzenia po nich całkowania. Kształt elementu kontrolnego może być dowolny, co jest jedną z przyczyn dużej popularności tej metody, zwłaszcza do rozwiązywania równań mechaniki płynów.
Adam Piechna
W niniejszym artykule zastosujemy metody objętości skończonych do zamodelowania dwóch zagadnień: przewodzenia ciepła w pręcie i rozkładu ciśnień w szczelinie pomiędzy pierścieniem tłokowym a ścianą cylindra.
Wykorzystanie metody objętości skończonych przedstawimy na przykładzie równania dyfuzji. Ogólna postać równania jest następująca:
gdzie :
Φ – zmienna, której rozkładu szukamy
Γ – współczynnik dyfuzji; w ogólnym przypadku może być zależny od x
S – człony źródłowe
Ideą metody objętości skończonych jest podział analizowanego obszaru na skończenie małe objętości tzw. objętości kontrolne i przeprowadzenie po nich całkowania. Schematyczny podział obiektu jednowymiarowego przedstawiono na rysunku 1.
Rys. 1 Dyskretyzacja (podział) obiektu 1D na objętości kontrolne. Analizowana objętość kontrolna ma środek w punkcie C. Indeksami L i P oznaczono kolejno: lewą komórkę i prawą komórkę. Indeksami l i p oznaczono powierzchnie pomiędzy komórką centralną, a komórkami po lewej i prawej stronie.
Po dokonaniu podziału interesującego nas obszaru możemy wykonać całkowanie po wydzielonych objętościach:
Stosując twierdzenie Ostrogradskiego-Gaussa możemy zamienić całki objętościowe na całki powierzchniowe. W najprostszym przypadku 1D, stosując oznaczenia jak na rysunku 1 otrzymamy:
Przyjmując Γ = const., A = const., dx = const. (W ogólnym przypadku Γ , A, dx mogą być funkcją x) oraz zamieniając różniczkę na różnicę otrzymamy:
Podstawiając otrzymujemy:
Porządkując otrzymujemy:
Jest to równanie algebraiczne liniowe opisujące zależności pomiędzy trzema elementami. Dzieląc analizowany obszar na przykład na pięć elementów musimy sformułować układ równań składający się z pięciu równań – po jednym dla każdego elementu kontrolnego. Równania dla elementu pierwszego i ostatniego będą różnić się od pozostałych. Elementy te są elementami brzegowymi i konieczne jest zdefiniowanie tzw. warunku brzegowego. Dwa typowe warunki to tzw. warunek Dirichlet’a (wartość zmiennej na brzegu) i Neumann’a (wartość pochodnej na brzegu). Rozpisując równania dla przykładu pięciu elementów otrzymamy układ równań :
W postaci macierzowej :
Równania w takiej postaci są bardzo proste do rozwiązania np. korzystając ze środowisk obliczeniowych GNU Octave (darmowe) lub Matlab.
Analiza rozkładu temperatury w pręcie
Stosując przedstawioną metodykę rozwiążemy zagadnienie: Obliczyć rozkład temperatury w aluminiowym pręcie. Pręt zanurzony jest w przepływającej wodzie o temperaturze 300 K. Oba końce pręta stykają się z płytą o temperaturze 500 K.
Rys. 2 Analizowany przypadek: rozkład temperatury w pręcie zanurzonym w płynącej wodzie
Pręt oddaje ciepło do powietrza i wody zgodnie z zależnością:
Równanie opisujące niniejsze zagadnienie, po uwzględnieniu członu źródłowego ma postać:
Współczynnik przewodzenia dla aluminium wynosi: k = 250 [W/K/m].
Współczynnik przejmowania ciepła α pręt – płynąca woda = 100 [W/m2/K]
Współczynnik przejmowania ciepła α pręt – powietrze = 5 [W/m2/K]
Po analogicznych do opisanych wcześniej przekształceniach otrzymujemy układ równań:
Porządkując otrzymujemy:
Poniżej kod rozwiązujący tak zdefiniowane zagadnienie (w programie GNU Octave):
Poniżej wyniki przykładowej symulacji:
Rys. 3 Rozwiązanie postawionego zadania. Zakropkowano obszar kontaktu z powietrzem. Falki oznaczają obszar kontaktu z płynącą wodą. Przy tak dużych różnicach we współczynniku przejmowania temperatura w obszarze kontaktu z powietrzem zmienia się praktycznie liniowo, jak przy warunku izolacji.
Analiza rozkładu ciśnień w szczelinie pomiędzy pierścieniem tłokowym a ścianą cylindra
Takie same podejście możemy zastosować do analizy rozkładu ciśnień w wypełnionej olejem szczelinie pomiędzy pierścieniem tłokowym a ścianą cylindra.
Rys. 4 Schemat analizowanego zagadnienia. Szukany rozkład ciśnienia w szczelinie pomiędzy pierścieniem tłokowym a ścianą cylindra.
Rozkład ciśnienia w szczelinie wypełnionej cieczą opisany jest równaniem Reynold'sa:
gdzie :
p – ciśnienie
h – lokalna wysokość szczeliny
U – prędkość pierścienia
μ – lepkość dynamiczna oleju
Przeprowadźmy analogiczny tok postępowania :
Tym razem współczynnik przy pochodnej jest również funkcją x:
Poniżej rozwiązanie problemu w środowisku GNU Octave. Konieczne było stworzenie dodatkowego skryptu-funkcji opisującej kształt pierścienia.
Poniżej otrzymane wyniki:
Rys. 5 Wyniki uzyskane w symulacji. Na górze – analizowane kształty pierścienia. Na dole – rozkład ciśnienia w szczelinie. Przerywaną linią oznaczony jest pierścień z defektem; obserwujemy znaczne zwiększenie panujących w szczelinie ciśnień.
Podsumowanie
W niniejszym artykule położyliśmy większy nacisk na pokazanie praktycznego wykorzystania metody niż na opis symulowanych zjawisk fizycznych.
Zachęcam do poeksperymentowania z zaproponowanymi modelami we własnym zakresie i do poświęcenia chwili na refleksję nad samą fizyką zjawisk. Przedstawiona metoda może zostać wykorzystana do rozwiązywania bardziej skomplikowanych dwu i trzywymiarowych zagadnień. Komplikuje się jednak wtedy proces dyskretyzacji analizowanego obszaru. Również liczba równań, którą trzeba rozwiązać wymaga zastosowania innego niż przedstawione podejścia.
Adam Piechna
artykuł pochodzi z wydania Czerwiec 6 (57) 2012