Obliczenie rozkładu pola prędkości

3.5 Obliczenie rozkładu pola prędkości

Algorytm stosowany do wyznaczenia rozkładu pola prędkości w metodzie DPIV opiera się na statystycznej analizie zarejestrowanych obrazów. Porównując ze sobą poszczególne zdjęcia, na których zapisany został obraz cząstek znacznikowych, poszukuje się podobieństwa między nimi, za co w analizie statystycznej odpowiada funkcja korelacji. W obliczeniach, w zależności od sposobu zapisu danych, wykorzystywane są dwa typy korelacji: autokorelacja i korelacja krzyżowa. Pierwszą z nich aplikuje się w przypadku, gdy pojedyncze zdjęcie naświetlane jest kilkakrotnie, zawierając jednocześnie obraz cząstek znacznikowych w kilku chwilach czasowych. W odwrotnej sytuacji, gdy na poszczególnej

klatce obrazu zapisany jest rozkład cząstek tylko i wyłącznie z danej chwili czasowej, stosowana jest funkcja korelacji krzyżowej. W prezentowanej rozprawie doktorskiej, poszczególne pola prędkości wyznaczane były w oparciu o algorytm wykorzystujący funkcję korelacji krzyżowej, dlatego też poniżej omówiona została metoda obliczeniowa oparta na jej zastosowaniu. Algorytmy wykorzystujące metody autokorelacji zostały szczegółowo opisane przez Westerweel’a [66] i Raffel’a i innych [21].

W ogólnym założeniu, celem metody DPIV jest obliczenie rozkładu wektorowego pola prędkości w oparciu o zarejestrowane obrazy cząstek znacznikowych. Aby było to możliwe, każde ze zdjęć w początkowym etapie analizy dzielone jest na tzw. okna interrogacji, czyli niewielkie (kwadratowe) obszary, w których zapisany jest pewien ściśle określony wzór (rozkład) cząstek znacznikowych (rys. 3.10a). Wielkość takich obszarów dostosowywana jest do warunków aktualnego eksperymentu i zależy głównie od gęstości posiewu (ilości cząstek znacznikowych przypadających na jednostkę objętości), szybkości przepływu badanego medium oraz rozdzielczości zapisanego obrazu. Typowe wymiary stosowanych okien interrogacji są równe 4x4, 6x6, 12x12, 16x16, 32x32, 64x64, 128x128, 256x256, 512x512 i 1024x1024 pikseli.

(a) (b)

Rys. 3.10 (a) Podział pojedynczej klatki obrazu na okna interrogacji, (b) Schemat wyznaczenia prędkości średniej na podstawie n par zdjęć zawierających chwilowe wartości prędkości cząstek

W każdym z takich obszarów wyznaczany jest jeden wektor prędkości. Ich całkowita ilość zależy od wielkości obranego okna interrogacji. Przykładowo, na zdjęciu o rozdzielczości

2048x2048 piksela, wybierając okno interrogacji o rozmiarze 32x32 piksela można wyznaczyć maksymalnie 4096 wektorów prędkości.

Prędkość płynu w danym punkcie, w metodzie DPIV wyznaczana jest ze wzoru (3.25) na prędkość ciała w ruchu jednostajnie prostoliniowym.

dr dU

dt

= (3.25)

Zgodnie z nim do wyznaczenia poszczególnego wektora prędkości konieczna jest znajomość przemieszczenia cząstki dr oraz czasu dt, w którym ono nastąpiło. W najprostszym przypadku tylko dwa zdjęcia wystarczają do wyznaczenia poszukiwanego, chwilowego rozkładu wektorowego pola prędkości. W takiej sytuacji, na pierwszym z obrazów zapisany jest rozkład cząstek znacznikowych w czasie t0.

Rys. 3.11 Schematyczna reprezentacja geometrii pomiarowej w metodzie DPIV [21]

Drugie zdjęcie przedstawia rozkład tych samych cząstek znacznikowych w czasie t1. Czas jaki upłynął pomiędzy rejestracją obydwu zdjęć to właśnie dt. Jest to wielkość znana i ustalana indywidualnie dla danego układu pomiarowego przez osobę wykonującą pomiar. Jego wartość nie przekracza kilku mikrosekund, stąd prędkość cząstki z dobrym przybliżeniem można określić wzorem (3.25). Aby uzyskać bardziej czytelną informację o rozpatrywanym

zjawisku konieczne jest wyznaczenie co najmniej kilkudziesięciu, a w praktyce nawet kilku tysięcy (w zależności od turbulencji rozpatrywanego przepływu) par zdjęć, na podstawie których w dalszej kolejności można oszacować średni rozkład pola prędkości w układzie (rys. 3.10b). Podstawowa trudność związana ze sposobem dokonywania pomiaru, wynika z prawidłowego oszacowania wielkości dr . Gdyby wektor prędkości wyznaczany był w oparciu o ruch przykładowo jednej cząstki znacznikowej, nie było by to, aż tak trudne zagadnienie, ponieważ cząstka taka byłaby w pełni rozróżnialna. Gdyby cząstek tych było kilka, również można by podjąć próbę ich identyfikacji. W rzeczywistości jednak wielkość przemieszczenia wyznaczana jest w oparciu o ruch znacznej grupy cząstek znacznikowych (chmury cząstek), dlatego niezbędne jest w tym celu wykorzystanie analizy statystycznej, na podstawie której, porównując ze sobą dwa obrazy z poszczególnej pary zdjęć poszukuje się struktur najbardziej do siebie podobnych.

Algorytm obliczający wartość funkcji korelacji krzyżowej dla poszczególnego okna interrogacji przedstawiono poniżej. Zgodnie z nim, biorąc przykładowy element objętości na powierzchni noża świetlnego (rys. 3.11) można utworzyć zbiór Γ zadany wzorem 3.26, w którym poszczególny element X_i określa położenie pojedynczej cząstki znacznikowej w czasie t, opisane składowymi X_i, Y_i, Z_i w kartezjańskim układzie współrzędnych.

1 2 ... n X X X       Γ =_{ }       , gdzie X_i =

(

X Y Z_i, ,_{i i}

)

(3.26) 1 2 ... n x x x γ       =_{ }     , gdzie x_i =

(

x y_i, _i

)

(3.27)

W podobny sposób można utworzyć zbiór

γ

zawierający informację o położeniu tych samych cząstek znacznikowych na płaszczyźnie matrycy CCD (równanie (3.27)), w którym współrzędne położenia danej cząstki wynoszą x_i, y_i. [21]. Ponieważ, w metodzie DPIV pomiaru dokonuje się na płaszczyźnie noża świetlnego, którego grubość ∆Z₀ jest dużo mniejsza od rozmiaru przestrzeni pomiarowej można założyć, iż∆Z₀ ≈ 0.

Ponadto, jeśli zbiór

γ

reprezentuje obraz wszystkich cząstek znacznikowych ze zbioru Γ, można założyć, iż co do rozmiaru Γ = γ. Najprostszą zależność przedstawiającą związek pomiędzy współrzędnymi cząstki na płaszczyźnie noża świetlnego i płaszczyźnie obrazu można przedstawić w postaci wzorów (3.28) i (3.29) [21], gdzie M jest wielkością stałą i oznacza powiększenie obrazu (jest to przypadek wyidealizowany).

i i x X M = (3.28) i i y Y M = (3.29)

W ogólnym przypadku, w którym uwzględnia się chociażby zniekształcenia obrazu wartość ta zadana jest nieco bardziej złożonym równaniem, które przykładowo przedstawione zostało w podrozdziale 3.3.4.

Zakładając, że powyższy zbiór, który możemy oznaczyć jako Γ₀ zawiera informację

o położeniu cząstek znacznikowych w czasie t0 (pierwsze zdjęcie) można utworzyć zbiór Γ₁, który zawierał będzie współrzędne położenia tych samych cząstek znacznikowych na zdjęciu drugim zarejestrowanych w czasie t₁. Zależność pomiędzy położeniem cząstek w zbiorze Γ₀ i Γ₁, przy pominięciu przemieszczenia w kierunku osi z można zapisać w postaci równania (3.30) [21], gdzie kolejne wielkości DXi reprezentują przemieszczenie cząstki w czasie

dt = t1 – t0. Podobnie jak w przypadku równań (3.28) i (3.29) związek pomiędzy przemieszczeniem cząstek na płaszczyźnie noża świetlnego DXi i na płaszczyźnie obrazu

i

dx można powiązać równaniem (3.31) [21].

Podczas pomiaru w metodzie DPIV znane jest położenie cząstek na płaszczyźnie obrazu w czasie t1 oraz t0. Znane więc są (uwzględniając wykonaną kalibrację układu optycznego) wszystkie wielkości w macierzy X1 oraz X0. Celem stosowanego algorytmu obliczeniowego jest wyznaczenie wszystkich wartości w macierzy DX. Aby tego dokonać wykorzystuje się m.in. funkcję korelacji krzyżowej.

Obraz cyfrowy zarejestrowany na matrycy CCD zapisany jest w postaci macierzy o wymiarach NxN, gdzie N jest wielkością zależną od rozdzielczości obrazu. Każdy element takiej macierzy (jeden piksel) zawiera informację o intensywności zarejestrowanego obrazu. Jeżeli założymy, że funkcja I0(x0,i, y0,i) reprezentuje poszczególne elementy macierzy obrazu

zarejestrowanego w czasie t₀, a funkcja I₁(x_1,i, y_1,i) elementy macierzy obrazu zarejestrowanego w czasie t₁, gdzie x_1,i = x_0,i + dx i y_1,i = y_0,i + dy, to wartość funkcji korelacji

krzyżowej w danym punkcie można obliczyć wg wzoru (3.32) [21].

Jak już zostało to zaznaczono powyżej, aby zastosować powyższy wzór do wyznaczenia wielkości przemieszczenia cząstek znacznikowych jakie nastąpiło pomiędzy czasem t0 i t1 w pierwszym kroku oba zdjęcia dzielone są na niewielkie podobszary, czyli okna interrogacji. 1 0 1 1 1 2 2 2 ... ... ... n _t n _t n DX X X X X DX X X _DX                 =  +_{ }→                     _{ } X1 = X0 + DX (3.30) i dx DX M = (3.31)

( )

( ) ( )

0 , 1, 0 , 1, 0 0, 0, 1 1, 1, , ,

, , ,

i i i i L K II i i i i x x L y y K

R dx dy I x y I x dx y dy

=− =−

=

∑ ∑

⋅ + +

(3.32)

Problem polega następnie na odnalezieniu na zdjęciu drugim struktury z przykładowego okna interrogacji (np. O₀, zgodnie z rys 3.12a). zapisanej w chwili t₀ na zdjęciu pierwszym. Najprościej można wyobrazić sobie, iż w celu zlokalizowania struktury oznaczonej symbolem O₀, ze zdjęcia wykonanego w czasie t₀, na zdjęciu wykonanym

w czasie t₁, należałoby wyznaczyć wartość funkcji korelacji krzyżowej dla każdej pary okien

interrogacji składającej się z okna interrogacji O₀ i każdego kolejnego okna interrogacji ze zdjęcia drugiego (rys. 3.12b). Byłby to jednak proces bardzo czasochłonny i wymagający wykonania znacznej liczby obliczeń. W rzeczywistości, cząstki zlokalizowane początkowo w oknie O0 w czasie t1 – t0 pokonują bardzo niewielką odległość. W algorytmach

obliczeniowych zakłada się nawet, iż nie może być ona większa niż połowa wielkości okna interrogacji [21], co w przypadku np. okna o wymiarach 32x32 piksela oznacza, iż przemieszczenie to nie może być większe niż 16 pikseli w każdym z kierunków. Dzięki temu obszar na zdjęciu 2 (rys, 3.12c), w którym wyszukuje się struktury cząstek, zawartych początkowo w oknie O0, znacznie się zawężą (w przypadku okna interrogacji o wymiarach 32x32 piksela wynosi on 64x64 piksela), co pozwala również skrócić czas niezbędny do wykonania obliczeń.

(a) (b) (c)

Rys. 3.12 (a) Przykładowe okno interrogacji na zdjęciu 1, zawierające pewien charakterystyczny rozkład cząstek znacznikowych zarejestrowanych w czasie t₀, (b) schemat wyszukania obszaru okna O0 na zdjęciu 2 polegający na porównaniu ze sobą każdej pary okien interrogacji, (c) schemat wyszukania obszaru okna O0 na zdjęciu 2 zawężający obszar poszukiwań do najbliższego otoczenia danego okna interrogacji

Schemat postępowania przy wyznaczaniu funkcji korelacji krzyżowej dla metody zilustrowanej na rysunku 3.12c można przedstawić w następujący sposób. Załóżmy, że macierz A (rys. 3.13) reprezentuje okno interrogacji o wymiarach 4x4 piksela, w którym poszczególne wielkości (wybrane przykładowo z przedziału 0 – 10) reprezentują wartość intensywności światła zapisaną na poszczególnym pikselu. Z kolei macierz B jest obszarem wyszukiwania o wymiarach 8x8 piksela. Zgodnie z przedstawionym schematem, w każdym kolejnym kroku obliczeniowym, wyznaczana jest wartość funkcji korelacji krzyżowej dla zadanego okna interrogacji i odpowiadającego mu rozmiarem okna, przesuniętego o kolejne wartości dx, dy w obszarze wyszukiwania. Przykładowo, dla przesunięcia dx = -1, dy = -1 wartość macierzy korelacji wynosi 421. Wynikiem kompletnego zbioru operacji jest macierz korelacji RII o wymiarach zadanych ogólnym wzorem (N+1)x(N+1), gdzie N oznacza

długość boku okna interrogacji, co w omawianym przykładzie daje macierz korelacji o wymiarach 5x5 piksela.

W zaprezentowanym przykładzie okno interrogacji O0 przesunięto do prawego rogu

okna wyszukiwania (dx = 2, dy = 2). W wyniku tego wartość korelacji krzyżowej w odpowiednim miejscu macierzy korelacji osiągnęła wartość największą. Obliczając w kolejnym kroku współrzędne środka okna interrogacji na zdjęciu pierwszym oraz współrzędne środka obszaru reprezentującego jego przesunięcie można wyznaczyć wektor przemieszczenia d , co ostatecznie pozwala oszacować prędkość cząstek znacznikowych,

a co za tym idzie, prędkość badanego płynu, w obrębie danego okna interrogacji. Powtarzając powyższy schemat postępowania dla kolejnych okien interrogacji oblicza się rozkład wektorowego pola prędkości w całym obszarze pomiarowym.

Rys. 3.13 Schemat wyznaczania kolejnych wartości funkcji korelacji krzyżowej dla przykładowego okna interrogacji o wymiarach 4x4 piksela i okna wyszukiwania o wymiarach 8x8 piksela

Powyższy przykład przedstawia jedynie ogólny schemat postępowania podczas wyznaczania wektorowego pola prędkości w metodzie DPIV. W rzeczywistości obliczenia są bardziej złożone. Po pierwsze rozmiar okna interrogacji jest większy, przez co sama macierz korelacji ma więcej elementów. Ponadto, podczas analizy danych pochodzących z pomiarów pojawia się szereg aspektów, które komplikują obliczenia. Są nimi m.in. spadek intensywności obrazu cząstek znacznikowych wraz ze wzrostem ich odległości od źródła światła, jak również szumy występujące podczas rejestracji obrazu oraz różna gęstość

posiewu. Pomimo, iż same obliczenia są bardziej złożone, to jednak co do zasady są one zgodne ze schematem przedstawionym w powyższym rozdziale. Przykładowy wykres funkcji korelacji krzyżowej dla okna interrogacji o wielkości 32x32 piksela przedstawiono na rysunku 3.14 (wykres sporządzono bezpośrednio w oparciu o dane pochodzące z badań przeprowadzanych w ramach prezentowanej rozprawy doktorskiej).

Rys. 3.14 Przykładowy wykres funkcji korelacji krzyżowej, dla pojedynczego wektora prędkości otrzymanej w czasie pomiarów, dla okna interrogacji o rozmiarze 32x32 piksela

W dokumencie Index of /rozprawy2/10505 (Stron 48-57)