14 (18)


Ć w i c z e n i e 13
Pomiar pola prędkości w przepływie turbulentnym
metodÄ… termoanemometrii
1. Wprowadzenie
Pomiary pola prędkości w przepływie turbulentnym są zwykle dokonywane za
pomocą techniki termoanemometrycznej. Stałotemperaturowy układ pomiarowy
termoanemometru (rys.1) pozwala zmierzyć napięcie E, które jest związane
zale\nością funkcyjną z prędkością czynnika U oraz jego temperaturą Śa [1].
Zale\ność ta wynika z bilansu cieplnego włókna sondy włączonej do zastosowanego
układu pomiarowego i zapisana mo\e być w postaci:
n
E2 = A + BUef (1)
w której Uef oznacza składową wektora prędkości, odpowiedzialną za efekt chłodzenia
gorącego włókna.
Wielkości A i B z równania (1) dla danego medium o temperaturze Śa = idem oraz dla
konkretnego czujnika o ustalonej temperaturze włókna Św = idem, mogą być
traktowane jako stałe, je\eli prędkość czynnika U nie przekracza 30 m/s. W praktyce
pomiarowej przyjmuje się najczęściej, \e wykładnik potęgowy n nie zale\y od
temperatury płynu i przy umiarkowanej prędkości przepływu (U < 30 m/s) jego
wartość mo\na uznać za stałą, wynoszącą n = 0.5 [2].
Je\eli dany przepływ turbulentny ma charakter stacjonarny1 wówczas zgodnie z
hipotezą Reynoldsa, chwilowe wartości parametrów fizycznych poruszającego się
czynnika traktować mo\na jako wynik superpozycji niezale\nych od czasu wartości
średnich oraz wielkości fluktuacyjnych. Składowe wektora prędkości chwilowej
elementów płynu, w dowolnie wybranym układzie współrzędnych prostokątnych,
mo\na więc przedstawić w postaci:
U1 = U1 + u1; U2 = U + u2; U3 = U + u3; (2)
2 3
lub
Ui = U + ui; i = 1,2,3
i
gdzie:
U , m/s - składowa prędkości średniej,
i
ui , m/s - składowa fluktuacji prędkości.
Dla uproszczenia rozwa\ań przyjmuje się zwykle, \e oś x1 układu współrzędnych jest
równoległa do kierunku wektora prędkości średniej, wówczas:
U1 = U; U = 0; U = 0 . (3)
2 3
W tak zorientowanym przepływie mo\na umieścić sondę pomiarową w taki sposób,
aby jej włókno było równoległe do płaszczyzny (x1, x2), tworząc z kierunkiem
1
Przepływ turbulentny ma charakter stacjonarny wówczas, gdy momenty statystyczne
wszystkich parametrów fizycznych tego przepływu nie zale\ą od czasu.
110
prędkości średniej kąt Ć (rys.1). Prędkość efektywna odpowiedzialna za chłodzenie
U
gorącego włókna określona jest wówczas wzorem:
2
Uef = [(U + u1)sinÕ - u2 cosÕ]2 + u3 (4)
Tak określona prędkość efektywna Uef jest składową prędkości chwilowej U w
kierunku prostopadłym do włókna sondy (rys.1). Składowa prędkości chwilowej
wzdłu\ włókna nie odgrywa istotnej roli w procesie chłodzenia o ile kąt między
kierunkiem prędkości chwilowej a włóknem jest większy od 20o [2].
Napięcie E mostka stałotemperaturowego związane jest z prędkością czynnika
zale\nością:
2
E = E + e = A + B4 [(U + u1)sinÕ - u2 cosÕ]2 + u3 (5)
Aby wydzielić z równania (5) składniki stałe, niezale\ne od czasu oraz składniki
fluktuacyjne, mo\na funkcję występującą po prawej stronie zale\ności (5) rozwinąć w
szereg potęgowy, otrzymując:
2
u1 u2 u3 u1
E + e = A + B U sinÕ + a1 + a2 + a3 + a4 2 +
U U U
U
(6)
2 2
u2 u3 u1u2 u1u3 u2u3
+ a5 2 + a6 2 + a7 2 + a8 2 + a9 2 + ...
U U U U U
111
Rys. 1. Sonda w układzie mostkowym termoanemometru
112
Współczynniki szeregu potęgowego (6) nie zale\ą od czasu i wyra\ają się podanymi
ni\ej zwiÄ…zkami:
B U sinÕ
a1 =
4 A + B U sinÕ
B U sinÕ
a2 = - ctgÕ
4 A + B U sinÕ
a3 = 0
3B U sinÕ + 2A
B2 U sinÕ
a4 = -
32(A + B U sinÕ)
B A + B U sinÕ U sinÕ
3B U sinÕ
B2U sinÕ
a5 = - ctg2Õ
32(A + B U sinÕ)
B A + B U sinÕ U sinÕ
2B U sinÕ + 2A
B2U sinÕ 1
a6 =
sin2 Õ
16(A + B U sinÕ)
B A + B U sinÕ U sinÕ
3B U sinÕ + 2A
B2U sinÕ
a7 = ctgÕ
16(A + B U sinÕ)
B A + B U sinÕ U sinÕ
(7)
a8 = 0
a9 = 0
...........
Po przeprowadzeniu uśrednienia w czasie szeregu (6) otrzymamy:
2 2 2
u1 u2 u3 u1u2
E = A + B U sinÕ + a4 2 + a5 2 + a6 2 + a7 2 + .... (8)
U U U U
lub
E + "E = A + B U sinÕ (8a)
przy czym:
2 2 2
ëÅ‚
u1 u2 u3 u1u2 öÅ‚
"E = -ìÅ‚ a4 2 + a5 2 + a6 2 + a7 2 + ....÷Å‚ (9)
ìÅ‚ ÷Å‚
U U U U
íÅ‚ Å‚Å‚
Z powy\szych związków wynika, \e wartość napięcia średniego E w ogólnym
przypadku zale\y nie tylko od prędkości średniej U (rys. 2), ale równie\ od wariancji
i kowariancji składowych fluktuacyjnych oraz uśrednionych iloczynów wy\szych
rzędów (iloczynów wy\szych rzędów nie uwzględniono w zapisie szeregu (6)).
113
Rys. 2. Odpowiedz termoanemometru na Rys. 3. Odpowiedz termoanemometru przy
zmiany prędkości czynnika niskim poziomie turbulencji przepływu
W przypadku, gdy rozwa\amy przepływ o umiarkowanym poziomie turbulencji
(składowe fluktuacyjne prędkości to w zale\nościach (8) i (9) mo\na
(ui << U)
pominąć wyra\enia rzędu drugiego oraz rzędów wy\szych jako wielkości
öÅ‚
nieskoÅ„czenie maÅ‚e ëÅ‚ uiu j H" 0; i, j = 1,2,3÷Å‚. Zale\noÅ›ci (8) i (9) przyjmÄ… wiÄ™c postać:
ìÅ‚
2
ìÅ‚ ÷Å‚
U
íÅ‚ Å‚Å‚
E = A + B U sinÕ
"E = 0
Takie przybli\enie oznacza, \e odpowiedz napięciowa mostka stałotemperaturowego
E(U) na zmianę prędkości chwilowej w zakresie od Umin do Umax
(dla "U = Umax  Umin << U ), przy ustalonej prędkości średniej (U = idem) w
wybranym punkcie pomiarowym mo\e być aproksymowana odcinkiem stycznej
wyznaczonej dla U = U (rys. 3).
Z przedstawionych rozwa\ań wynika równie\, \e w przypadku ui << U [1, 2]
słuszna jest zale\ność:
B U sinÕ
E + e = A + B U sinÕ + (u1 - u2ctgÕ) (10)
4U A + B U sinÕ
w której wyró\nić mo\na składnik niezale\ny od czasu:
E = A + B U sinÕ (11)
oraz składnik fluktuacyjny:
B U sinÕ
e = (u1 - u2ctgÕ) = s(u1 - u2ctgÕ) (12)
4U A + B U sinÕ
Występująca w zale\ności (12) czułość s jest wielkością charakteryzującą układ
pomiarowy termoanemometru i dana jest wzorem:
"E B U sinÕ B2
s = = = sinÕ (13)
2
"U
4U A + B U sinÕ 4EëÅ‚ E - AöÅ‚
ìÅ‚ ÷Å‚
íÅ‚ Å‚Å‚
114
Analiza związków (10), (12) i (13) wykazuje, \e sonda o włóknie prostopadłym do
Ä„
ëÅ‚Õ öÅ‚
kierunku przepływu średniego = rejestruje, poza prędkością średnią, jedynie
ìÅ‚ ÷Å‚
2
íÅ‚ Å‚Å‚
składową wzdłu\ną fluktuacji prędkości:
E + e = A + B U + su1 (14)
Nale\y jednak wyraznie zaznaczyć, \e związek (14) zachowuje swą wa\ność tylko w
przypadku pomiarów w przepływie turbulentnym o ustalonej temperaturze
Śa = const , w którym wartość skuteczna fluktuacji temperatury RMS(Ń) = Ń2 H" 0 .
Badanie pola prędkości w przepływie turbulentnym, przy u\yciu termoanemometru
współpracującego z sondą o pojedynczym włóknie prostym, polega na określeniu w
wybranych punktach pomiarowych, prędkości średniej U oraz wartości skutecznej
2
składowej wzdłu\nej fluktuacji prędkości
RMS(u1 ) = u1
Przed wykonaniem badań zasadniczych nale\y wyznaczyć parametry stałe A i B
występujące w związku (14) poprzez wzorcowanie układu pomiarowego.
115
Rys. 4. Schemat stanowiska badawczego
116
2. Stanowisko badawcze
Stanowisko badawcze i zastosowana aparatura, pokazana w sposób schematyczny
na rys. 4, pozwalają określić zarówno rozkłady prędkości średniej i wartości skuteczne
fluktuacji prędkości w kołowym przepływie swobodnym, jak i współczynniki stałe A i
B równania równowagi układu pomiarowego (vide równ. (1)).
Z dyszy o średnicy d = 40 mm wypływa powietrze o ustalonej temperaturze
. Dysza ta jest jednocześnie zwę\ką pomiarową, poniewa\ ró\nica ciśnień
Åša = const
statycznych w dwóch przekrojach kontrolnych 1 i 2 jest funkcją prędkości średniej .
U
Sonda podłączona jest do mostka stałotemperaturowego CTA 55M10, którego
napięcie wyjściowe jest zale\ne od prędkości chwilowej czynnika. Włączony w tor
pomiarowy woltomierz cyfrowy 55D31 posiada układ całkujący z regulowanym
czasem uśredniania.
Układ taki pozwala zmierzyć składową stałą sygnału napięciowego która zgodnie z
E
zale\nością:
2
E = A + B U (15)
Jest funkcją prędkości średniej U . Woltomierz RMS 55D35 umo\liwia pomiar
wartości skutecznej składowej zmiennej sygnału napięciowego, która jest miarą
fluktuacji prędkości:
2
e2 = s u1 (16)
Podłączony do mostka CTA 55M10 oscyloskop umo\liwia obserwację przebiegów
czasowych napięcia e(t).
3. Metodyka pomiarów i obliczeń
Ćwiczenie składa się z dwóch części. W części pierwszej nale\y wykonać
wzorcowanie układu pomiarowego, natomiast część druga obejmuje badanie pola
prędkości w kołowej strudze swobodnej.
3.1. Wzorcowanie układu pomiarowego
Wzorcowanie układu pomiarowego termoanemometru polega na wyznaczeniu jego
odpowiedzi napięciowej na zmianę prędkości poruszającego się czynnika i kreśleniu
wartości parametrów stałych równania (15). Przepływ, w którym wykonuje się
wzorcowanie winien charakteryzować się jednorodnym rozkładem prędkości średniej
2
U , ustalonÄ… temperaturÄ… Åša i niskim poziomem turbulencji ëÅ‚ u1 << U i Ń 2 H" 0öÅ‚.
ìÅ‚ ÷Å‚
íÅ‚ Å‚Å‚
Warunki te spełnia przepływ w jądrze potencjalnym kołowej strugi swobodnej (rys.4).
Z postaci związku (15) wynika, \e w układzie ( U ; E2) zale\ność E(U) ma
charakter liniowy. Obserwacje doświadczalne wskazują jednak, \e liniowość ta jest
zachowana w ograniczonym zakresie prędkości czynnika, nie większej od 30 m/s.
Dlatego cechowanie układu pomiarowego nale\y wykonać w zakresie niskich i
umiarkowanych prędkości przepływu, przy ustalonej temperaturze włókna sondy
Åšw = const i ustalonej temperaturze czynnika .
Åša = const
117
Przygotowanie termoanemometru do pomiarów
Przygotowanie termoanemometru do pomiarów wymaga właściwego dostrojenia
mostka CTA 55M10 i jednostki podstawowej 55M01 do podłączonej sondy. Nale\y
równie\ zmierzyć rezystancję włókna sondy Ra w temperaturze przepływu Śa = const
i
ustalić jego temperaturę Św = const poprzez zastosowanie przegrzewu m = Rw/Ra =
1.8.
Czynności przygotowawcze mo\na wykonać według ni\ej podanej instrukcji
opracowanej według [3] (patrz rys. 5).
Rys. 5. Płyta czołowa termoanemometru stałotemperaturowego
1. Włączyć zasilanie wszystkich jednostek układu pomiarowego.
2. Nastawy wstępne:
SQUARE WAVE : OFF
HF FILTER : 3
VOLTS : 10
FUNCTION : STD.BY
PROBE TYPE : WIRE
GAIN : 4
Oporność dekady : 00,00 OHMS.
3. Styki uchwytu sondy zewrzeć końcówką zerującą.
4. Kompensacja oporności przewodów:
a) przełącznik FUNCTION ustawić w poło\eniu RES.MEAS.,
b) wyregulować ZERO OHMS tak, aby wskazówka miernika wychyliła się do
czerwonej kreski,
c) przełącznik FUNCTION ustawić w poło\eniu STD.BY, zdjąć końcówkę
zerującą i podłączyć sondę pomiarową.
5. Pomiar rezystancji włókna sondy:
a) sondę umieścić w przepływie o ustalonej temperaturze Ś (w jądrze
a
potencjalnym strugi),
b) ustawić przełącznik FUNCTION w poło\eniu RES.MEAS.,
c) pokrętłami dekady doprowadzić wskazówkę miernika do czerwonej kreski
skali,
118
d) nastawa dekady wskazuje rezystancję włókna sondy Ra w temperaturze Ś ,
a
e) ustawić przełącznik FUNCTION w poło\enie STD.BY,
f) rezystancję Ra pomno\yć przez stopień przegrzewu m = 1.8; uzyskana wartość
Rw = 1.8 Ra nastawić na dekadzie MAIN UNIT 55M01.
6. Dostrojenie MAIN UNIT i CTA STANDARD BRIDGE do sondy pomiarowej:
a) przełącznik FUNCTION ustawić w poło\enie OPERATE,
b) SQUARE WAVE nastawić na 30 kHz,
c) Pokrętłami Q i L CABLE COMPENSATION oraz HF FILTER i GAIN
wyregulować sygnał wyjściowy mostka CTA tak, aby przebieg czasowy
napięcia uzyskany na ekranie oscyloskopu miał kształt pokazany na rys. 6.
d) SQUARE WAVE przestawić w poło\enie OFF.
7. Woltomierz 55D31 wskazuje napięcie średnie E , woltomierze 55D35 i V531
(Rys. 4) wartość skuteczną napięcia zmiennego e2 , oscyloskop pokazuje
przebieg czasowy napięcia e(t).
Rys. 6. Odpowiedz układu pomiarowego termoanemometru na prostokątną falę testującą
Regulacja prędkości średniej i pomiar napięcia średniego
Prędkość powietrza wypływającego z dyszy (rys. 4) mo\na określić ze wskazań
mikromanometru podłączonego do króćców pomiarowych, wykorzystując wzór:
2Ám Å" g Å"l Å"i
U = , m/s
(17)
i
Ä… Å" Áa
w którym:
Ám , kg/m3 - gÄ™stość cieczy manometrycznej,
i - przeło\enie mikromanometru,
g, m/s2 - przyspieszenie ziemskie,
Áa , kg/m3 - gÄ™stość powietrza,
ą - stała dyszy pomiarowej,
l, m - długość słupa cieczy manometrycznej równowa\ącego ró\nicę
ciśnień statycznych, istniejącą między przekrojami pomiarowymi
dyszy.
119
Prędkość Ui mo\na zmieniać za pomocą pokrętła potencjometru podłączonego do
układu sterującego prędkością obrotową wentylatorów. Potencjometr regulacyjny
umieszczony jest w dogodnym miejscu przy stanowisku pomiarowym.
Występującą we wzorze (17) gęstość powietrza mo\na wyliczyć z zale\ności:
pa
Áa = , kg/m3 (18)
R Å" Åš
a
w której:
pa, N/m2 - ciśnienie statyczne w strudze,
Śa , K - temperatura przepływu wyra\ona w skali bezwzględnej Kelvina,
R = 287,4 J/kg ·K - staÅ‚a gazowa powietrza.
Ciśnienie statyczne w przepływie swobodnym pa mo\na uznać za równe ciśnieniu
atmosferycznemu. Wyliczona ze wzoru (18) gęstość dotyczy powietrza traktowanego
jako gaz doskonały. Przyjęcie takiego uproszczenia jest dopuszczalne, poniewa\
zakres zmian parametrów przepływu w przypadku prowadzonego eksperymentu jest
stosunkowo mały.
Napięcie Ei termoanemometru odpowiadające nastawionej prędkości U
i
wskazuje woltomierz cyfrowy 55D31 po nastawieniu odpowiedniego czasu
uśredniania. Wyniki pomiarów nale\y wpisać do tabeli 1.
3.2. Pomiar pola prędkości
W wybranym przekroju pomiarowym x1 = idem turbulentnej, swobodnej strugi
kołowej (rys. 4) nale\y wyznaczyć rozkład prędkości średniej oraz rozkład wartości
skutecznej składowej wzdłu\nej fluktuacji prędkości. Wielkości te mo\na określić w
oparciu o pomiary E i e2 w wybranych punktach pomiarowych ustalonego
przekroju.
Odpowiednie wzory obliczeniowe mają postać:
- prędkość średnia:
2
2
ëÅ‚
E - AöÅ‚
ìÅ‚ ÷Å‚
U = , m/s (19)
ìÅ‚ ÷Å‚
B
íÅ‚ Å‚Å‚
- wartość skuteczna RMS(u1):
1 1
2
u1 = e2 = KRMS W , m/s (20)
s s
gdzie KRMS oznacza wybrany zakres woltomierza RMS 55D35, a W [V]  wskazanie
woltomierza V531,
- czułość układu pomiarowego:
B2 V
s = , . (21)
2
4EëÅ‚ E - AöÅ‚ m/s
ìÅ‚ ÷Å‚
íÅ‚ Å‚Å‚
Występujące w powy\szych wzorach wielkości A i B nale\y wyznaczyć w sposób
opisany w rozdziale 3.1.
120
4. Szczegółowy program badań
Kolejność postępowania podczas realizacji ćwiczenia jest następująca:
1. Po ustaleniu wartości Śa , Ra i Rw = 1.8 Ra i dostrojeniu układu pomiarowego
termoanemometru do zastosowanej sondy, przeprowadza siÄ™ jego wzorcowanie
zmieniając prędkość przepływu powietrza poprzez zmianę prędkości obrotowej
wentylatorów. Uzyskane dane pomiarowe zamieszcza się w kolumnach 2 i 3
tabeli 1.
2. PrÄ™dkość powietrza U oblicza siÄ™ wg wzoru (17), a gÄ™stość powietrza Áa wedÅ‚ug
i
(18). Po wypełnieniu kolumn 4, 5, 6 tabeli 1 nanosi się punkty (xi, yi) w układzie
( U , E2) .
3. Wartość parametrów stałych A i B równania równowagi termoanemometru oblicza
się metodą najmniejszych kwadratów według wzorów [4]:
n xi yi - (" xi)(" yi) (22)
"
B =
n xi2 - (" xi )2
"
("xi2)(" yi)- ("xi)("xi yi) (23)
A =
n - ("xi)2
"x2
i
Ocenę dokładności wyznaczonych wartości A i B mo\na dokonać poprzez
obliczenie ich odchyłek standardowych:
yi2 - B xi yi - A yi
n
" " "
sB = (24)
n - 2
n - ("xi)2
"x2
i
1
sA = sB xi2 (25)
"
n
i określenie przedziałów ufności dla przyjętego poziomu istotności ą przy n
niezale\nych pomiarach:
P { A - tn,Ä… sA d" Arzeczyw. d" A + tn,Ä… sA } = 1  Ä… (26)
P { B - tn,Ä… sB d" Brzeczyw. d" B + tn,Ä… sB } = 1  Ä… (27)
WartoÅ›ci sum wystÄ™pujÄ…cych we wzorach (22 ÷ 25) sÄ… zawarte w tabeli 1.
Wielkość tn,ą jest zmienną losową rozkładu Studenta, której wartość dla znanego n
i przyjętego ą mo\na odczytać z odpowiednich tablic statystycznych.
4. Uzyskaną funkcję (liniowa w układzie ( U , E2) ):
E2 = A + B U (28)
przedstawić nale\y w postaci graficznej. Charakterystykę taką dla konkretnej
sondy w postaci przykładu pokazano na rys. 7.
5. Po wycechowaniu termoanemometru przeprowadza się pomiar rozkładu wielkości
2
U i u1 w płaszczyznie x1 = idem, przy ustalonej prędkości wypływu powietrza z
dyszy U .
w
4. Wyniki zamieszcza siÄ™ w kolumnach 1, 2, 3, 4 tabeli 2.
2
Wielkości U , u1 , s oblicza się z zale\ności (19), (20) i (21). Uzyskane wyniki w
postaci zredukowanej przedstawia się w formie graficznej. Przykładowe rozkłady
121
prędkości średniej i wartości skutecznej wzdłu\nej składowej fluktuacji prędkości
przy Re = 40 000 pokazano na rys. 8.
Rys. 7. Charakterystyka termoanenometru dla konkretnej sondy pomiarowej
Rys. 8. Rozkład prędkości średniej i fluktuacji prędkości w przekroju pomiarowym
strugi kołowej
Uzyskane wyniki badań pozwalają odpowiedzieć na następujące pytania:
- Jak zmienia się prędkość średnia przepływu w kierunku normalnym do osi
swobodnej strugi kołowej?
- Jaki jest rozkład poziomu turbulencji w przekroju pomiarowym?
Literatura
1. Elsner J.W.: Turbulencja przepływów, PWN, Warszawa 1987
2. Hinze J.O.: Turbulence, 2nd ed. New York, Mc Graw-Hill 1975
122
3. Instruction Manual DISA 55M System with 55M10 CTA Standard Bridge
4. Szydłowski H.: Teoria pomiarów, PWN, Warszawa 1981
Tabele pomiarowo-obliczeniowe
Tabela 1
Dane ogólne:
temperatura medium Åša= K
ciśnienie statyczne pa = N/m2
gÄ™stość powietrza Á = kg/m3
gÄ™stość cieczy manometrycznej Ám = kg/m3
stała gazowa powietrza R = 287,04 J/kg K
przeło\enie manometru i = -
stała dyszy pomiarowej ą = -
rezystancja sondy w temp. Åša Ra = &!
rezystancja sondy w temp. Åšw Rw = &!
li Ei Ui
xi = Ui yi = Ei2
L.p. xi yi xi2 yi2
mm V m/s m0.5s-0.5 V2
1 2 3 4 5 6 7 8 9
1
2
3
n
n
"
i=1
123
Tabela 2
Parametry stałe układu pomiarowego:
A = V2; B =
V2/ m/s
temperatura otoczenia Åšot = K
ciśnienie otoczenia pot = N/m2
Parametry fizyczne czynnika w płaszczyznie wylotowej dyszy:
temperatura czynnika Åša = K
ciśnienie statyczne pa = N/m2
średnica dyszy d = m
stała dyszy pomiarowej ą = -
gÄ™stość cieczy manometrycznej Ám = kg/m3
przeło\enie manometru i = -
stała gazowa czynnika R = J/kg K
gÄ™stość czynnika Áa = kg/m3
kinematyczny współczynnik lepkoÅ›ci ½ = m2/s
prędkość czynnika U = m/s
w
liczba Reynoldsa Re = -
współrzędne płaszczyzny pomiarowej x1 = m
zredukowana odległość płaszczyzny od wylotu dyszy x1/d = -
2
2
x2 KRMS W s x2/d U/Uw u1
E
u1
Uw
mm V - V m/s V·s/m m/s - - -
1 2 3 4 5 6 7 8 9 10
124


Wyszukiwarka

Podobne podstrony:
04 Rozdzial 14 18
Farmacja 18 12 14 II termin
18 (14)
18 14 Luty 2000 Ognie pod żelazną bramą
18 (14)
18 03 14 Pietrzyk
2012 01 16 14 04 18
dictionary 18 14

więcej podobnych podstron