Ćwiczenie 5a: Adwekcja w dwóch wymiarach

Ćwiczenie 5a: Adwekcja w dwóch wymiarach ^∗

Na poprzednich zajęciach (laboratorium nr 4) wyznaczyliśmy pole prędkości cieczy płynącej przez rurę. Znając prędkości, zasymulujemy zachowanie plamy oleju rozlanej na niej. Przyjmiemy pudło obliczeniowe takie, jak w zestawie o przepływie lepkim nieściśliwym: (i, j)∈ [0, 300] × [0, 90] z dx = dy = 0.01, gdzie (x, y) = (i· dx, j · dy) - zgodnie z ilustracją 1 z treści laboratorium nr 4.

Gęstość początkowa plamy oleju dana jest wzorem:

ρ⁰(x, y) = ρ(x, y, t = 0) = exp[

−25(

(x− 0.4)²+ (y− 0.45)²)]

. (1)

ρ(x, y, t) oznacza gęstość, natomiast prędkość w kierunku osi OX to u(x, y), a w kierunku osi OY : v(x, y).

Zbadamy unoszenie plamy ρ(x, y, t) w czasie. Zjawisko opisuje równanie ad- wekcji, które dla cieczy nieściśliwej (znikająca dywergencja prędkości) przyjmuje postać

∂ρ

∂t + u(x, y)∂ρ

∂x + v(x, y)∂ρ

∂y = 0. (2)

Równanie (3) zdyskretyzowane zgodnie ze schematem leapfrog ma postać ρⁿ⁺¹_ij = ρⁿ_ij⁻¹− ∆t

( uij

ρⁿ_i+1,j− ρⁿi−1,j

dx + vij

ρⁿ_i,j+1− ρⁿi,j−1

dy

)

. (3)

Metoda leapfrog nie jest iteracyjna - pętlę po całym pudle obliczeniowym prze- prowadzamy tylko raz dla każdego t (indeks górny n oznacza chwilę czasową).

Schemat wymaga zadania rozkładu gęstości nie tylko w pierwszej, ale i w drugiej chwili czasowej (potrzebne są więc trzy macierze: jedna dla obliczanej w danej chwili gęstości ρⁿ⁺¹ i dwie dla dwóch poprzednich: ρⁿ, ρⁿ⁻¹). Przyjmiemy na- stępujący przepis na drugą chwilę czasową:

ρ¹_ij= ρ⁰(x_i− uij∆t, y_j− vij∆t) (4) (gęstość w każdym punkcie x_i, y_iunoszona przez czas ∆t przez lokalną prędkość - jest to dobre przybliżenie, gdy prędkość słabo zmienia się tam, gdzie plama jest wstawiona jako warunek początkowy).

∗Laboratorium z inżynierskich metod numerycznych, Wydział Fizyki i Informatyki Sto- sowanej AGH 2014/2015. Bartłomiej Szafran (bszafran@agh.edu.pl), Krzysztof Kolasiński (kolasinski@fis.agh.edu.pl), Elżbieta Wach (Elzbieta.Wach@fis.agh.edu.pl), Dariusz Żebrowski (Dariusz.Zebrowski@fis.agh.edu.pl)

1

Teraz i w następnych zadaniach: krok czasowy przyjmiemy: ∆t = _4·max^dx , gdzie max =√

u²+ v² to maksymalna prędkość cieczy w całym pudle (należy ją najpierw znaleźć). Zastosujemy periodyczne warunki brzegowe: wzór (3) trzeba zaprogramować osobno dla krańców pudła - tak, aby lewym sąsiadem punktu (i = 0, j) był punkt (i = 300, j) i vice versa.

W symulacjach jako test poprawności w każdym kroku czasowym śledzimy całkę z gęstości:

I(t) =

∫

ρ(x, y, t)dxdy. (5)

Ma się nie zmieniać w sposób znaczący. Śledzimy również środek pakietu:

⟨x⟩ =

∫xρ(x, y, t)dxdy

I(t) . (6)

Symulację przerywamy po t≃ 100.

Zadanie 1: Adwekcja w rurze bez zastawki - schemat leapfrog.

Zbadamy adwekcję w przepływie lepkim przez rurę. Zadajemy paraboliczny rozkład prędkości u(x, y) = _2µ^Q(y− ymin)(y− ymax), v(x, y) = 0 z ymin= 0.0, ymax= 0.9, gradient ciśnienia Q = −1, i współczynnik lepkości µ = 1. Wpro- wadzamy plamę oleju w dwóch pierwszych chwilach czasowych: ρ⁰_ij, ρ¹_ij i sto- sujemy schemat (3) z pominięciem górnego i dolnego brzegu rury. Narysować I(t) i⟨x⟩(t):20 pkt. Zrobić 5 zdjęć plamie w czasie ruchu (tj. narysować mapę gęstości ρ(x, y, t) w pięciu wybranych chwilach czasowych t):20 pkt.

Zadanie 2: Adwekcja w rurze z przegrodą - schemat leapfrog.

Do rury wstawimy przegrodę o wymiarach identycznych, jak w laborato- rium nr 4: i1= 85, i2= 101, i3= 116, j1= 50, j2= 70. Przyjmujemy Q =

−1. Prędkości należy wyliczyć ze swojego programu z zeszłego tygodnia (przypa- dek pierwszy zadania 2) lub użyć pliku predkosc.dat, do którego link znajduje się na stronie. W pliku dane są ułożone w następujący sposób:

x y u(x, y) v(x, y)

Dla liczących samodzielnie: proszę zadbać o to, żeby obydwie składowe prędkości na przeszkodzie były równe dokładnie zeru.

Zbadać, jak plama opływa przeszkodę. Należy narysować całki I(t), ⟨x⟩(t) (20 pkt) oraz gęstości w pięciu wybranych chwilach czasowych (20 pkt).

Zadanie 3: Adwekcja w rurze z przegrodą - schemat Laxa-Friedrichsa.

Zbadać problem przepływów w rurze z przeszkodą z zadania drugiego sche- matem Laxa-Friedrichsa, który wprowadza dyfuzję (numeryczną)

ρⁿ⁺¹_ij = ρⁿ_i+1,j+ ρⁿ_i−1,j+ ρⁿ_i,j−1+ ρⁿ_i,j+1 4

− ∆t (

u_ijρⁿ_i+1,j− ρⁿi−1,j

2dx + v_ijρⁿ_i,j+1− ρⁿi,j−1

2dy )

. (7)

Wstawić dwukrotnie mniejszy krok czasowy, niż w metodzie leapfrog. Wzór (3) gwarantował, że plama nie wpłynie na przeszkodę (dokładna adwekcja pro- wadzi plamę wzdłuż linii strumienia, nigdy więc plama nie trafi w przeszkodę).

2

Dla schematu Laxa trzeba o to zadbać - nie przeprowadzamy obliczeń na brze- gach ani wewnątrz przeszkód.20 pktza 5 zdjęć w wybranych chwilach t.

3