Автомодельные течения вязкого газа в каналах с теплои массообменом на стенке

Тип работы:
Реферат
Предмет:
Механика


Узнать стоимость

Детальная информация о работе

Выдержка из работы

_________УЧЕНЫЕ ЗАПИСКИ Ц, А Г И
Т о м VII 19 7 6
№ 2
УДК 532. 542:532. 526
АВТОМОДЕЛЬНЫЕ ТЕЧЕНИЯ ВЯЗКОГО ГАЗА В КАНАЛАХ С ТЕПЛО- И МАССООБМЕНОМ НА СТЕНКЕ
А. П. Бир кин
Рассмотрены автомодельные течения вязкого газа в плоских и осесимметричных каналах с прямолинейными стенками, в которых параметры потока по длине изменяются по степенному или экспоненциальному закону. Такие течения реализуются при одном и том же профиле чисел М в поперечных сечениях канала и при наличии тепло- и массообмена на стенке. Построено автомодельное решение полной системы уравнений Навье-Стокса для течений газа в цилиндрических каналах, когда параметры потока по длине изменяются по экспоненциальному закону. Предельным состоянием рассматриваемых течений газа при М -" 0 является течение Пуазейля.
1. Будем рассматривать ламинарное стационарное течение газа в плоском канале постоянного сечения или цилиндрической трубе при наличии отсоса (вдува) на стенке при таких значениях чисел Ие, когда влияние вязкости распространяется на все поперечное сечение канала.
В предположении, что канал длинный, (у ш — высота или
радиус канала, — длина), и скорость отсоса (вдува) является величиной порядка ы^/Ие (и*- скорость газа на входе в канал), для описания течения вязкого газа будем использовать приближенные уравнения, совпадающие по форме с уравнениями пограничного слоя [1, 2]
-- ди р и — дх
, — - ди
1- ю — =
ду
-= (ри) 4-дх
дх
+ - у ду
і- ди
ду
4 г 4= (У" ?*) = о,
У дУ
ри Л= (А-?-
дх
-М* и2
+а±(н+
д
ду
-ч — дк
ду
+ (х-і)м-±?=
у дУ
-V- - ди
У
ду
рхМ1 = ріг, ?а = /г& quot-.
(1)
Система (1) выписана применительно к совершенному газу и степенной зависимости коэффициента вязкости от температуры (энтальпии). Здесь все величины безразмерные и связаны с размерными соотношениями
— X — у — и ^ п
х =- я-, у = -, и- -, v= - ке,
уш Ие ' ут ' ы* ' и*
Т Н — р — р — и.
— А 1 Р 2 & gt- Р р & gt- ^ & gt-
«* р* И Р*
где х, у — расстояния, измеряемые вдоль оси и по нормали к оси канала- и, V — продольная и поперечная составляющие скорости- Л, р, [л, М — соответственно энтальпия, давление, плотность, коэффициент динамической вязкости, число М- и*, Л*, р*, р*, М*-
значения соответствующих величин на оси канала в начальном
сечении х = 0- х — показатель адиабаты- Рг — число Прандтля- Яе = р* и* Уо,/р* - число Рейнольдса- индекс ги отвечает условиям на стенке канала, величина V принимает значение 0 в плоском и
1 — в осесимметричном случае. В рассматриваемом приближении давление постоянно по сечению канала, т. е. р = р (х).
Система уравнений (1) должна решаться при следующих граничных условиях
ди дк п «- А
-= = -===(), -Ц = 0 при «у — 0, оу ду
и = 0, v = vw (x), к = Ът (х) при у — ,
(2)
где Ут (х), /ги, (х) — заданные функции. В качестве начальных условий должны быть заданы профили величин и = и (у), Н = Ъ (у) и значение М* при х = 0..
Будем далее искать такие решения системы уравнений (1) с граничными условиями (2), которые имеют вид (черту в обозначениях опускаем)
и = и1(х)и2(у), v — v1(x)vг (y), к = к1(х)^(у),
______/.ч. /& lt-л (3)
%М2
Подстановка соотношений (3) в систему (1) показывает, что искомые решения будут иметь место при условиях
/ Г
= -55а = а — С0Пз1' -42 =* = СОП8^ -Ёр- = с = СОП81,
'-У 1
где а, Ь, с — произвольные константы, причем величина Ь без ограничения общности может быть положена равной единице- штрих означает дифференцирование по х.
Из первого условия при этом следует, что для рассматриваемых течений профили чисел М в поперечных сечениях канала, соответствующие продольной составляющей скорости, не изменяются по длине канала. Определив связь между р1 и ии указанные условия можно представить в следующем виде:
Л, = «2, р, = ис/а, = ~2п+2_с/а = а.
и
При интегрировании последнего соотношения получаем два следующих закона изменения величины их (и соответственно величин
Л» ри от
при 2 и-|-2 — са ф 1 имеем степенной закон
и, = (1 + кх) а' - при 2п + 2 — са =1 имеем экспоненциальный закон
и1 = е^'х,
где к, а» р* - новые произвольные константы.
Рассмотрим в отдельности оба этих случая,
а) Степенной класс. В соответствии с полученными выше результатами решения (3) для этого класса нужно искать в следующем виде:
и = (1 + и2(у), V = (1 + кх) а& gt- v2 (у),
А -(1 + /гху*И2(у), р = -^ (1 + кхУ' «(4)
р = (1 + ?*)"Е РгЫ.
где постоянные а2, а3, а4, а5 определяются подстановкой выписанных соотношений в исходную систему уравнений (1). В результате получаем
а2 = а, -1, а3 = 2а., а4 = 1 + (2я + 1), а5 = 1 + {2п — 1),
а а, и 4 считаются параметрами. При этом система обыкновенных дифференциальных уравнений для определения функций щ{у), (У), й2 (У), ?Лу) имеет вид '-
с*! Ар, и* + р, V-, и' = - к -1 + °д + ^ (-У'- ?ь
у '- & quot-
(1 + 2 а, я) ?р2 ы2 + (_у'-- р2 ъ2)'- = 0-
2 а 1А р2 Мг --2- Мо и21 -(- р2Х'2 -|-2- Мо Нг1 *=
= -^гу (Г К)'- + (х — ]) моу (У& quot- н-2 щ «?У-
р2= 1/Л2-2 =2.
Решение системы (5) должно удовлетворять следующим условиям на оси канала:
«2 (0) = /?2 (0) = 0, и2(0) = Л2(0) = 1, ®2(0)~0, (6)
которых достаточно, чтобы сформулировать задачу Коши для этой системы.
При этом величины а, и М0 (число М на оси) для заданных х, Рг, п будем считать параметрами. Величина, А должна определиться в результате интегрирования уравнений (5) из условия равенства нулю величины и2 на стенке канала (условие прилипания). Приведенные энтальпия газа на стенке Л2(1) и поперечная скорость г& gt-2(1) будут искомыми величинами.
Это объясняется тем, что в данном случае законы изменения всех газодинамических величин по длине определены, поэтому для реализации рассматриваемых законов при заданных а, и М0 величины Л2(1) и г& gt-2(1) должны быть найдены из решения всей задачи в целом.
(5)
В общем случае в рассматриваемых течениях должен иметь место как тепло-, так и массообмен. Требованию как сле-
дует из второго уравнения (5), отвечает условие
1+204» = 0, ocj = - ½ и.
б) Экспоненциальный класс. Для данного класса решения (3) следует искать в виде
и = е^х и2 (у), v = е$-х Vo (у), h = e^h2{y), |
Р = ^'-Х& gt- Р = е-**Р2(У) — j
Постоянные р2, Рз& gt- ?4& gt- Рб) как и в предыдущем случае, найдем в результате подстановки условий (7) в исходную систему уравнений (1). В результате получим
Pi = Pi, Рз = 2р1(р4 = р1(2л+1)& gt- Рв — Рг (2 я — 1),
где р, является параметром.
В отличие от степенного класса здесь продольная и поперечная составляющие скорости по длине изменяются одинаковым образом.
Система обыкновенных дифференциальных уравнений для определения функций и2, v2, /г2, р2 в данном случае будет
р, р2 и + р2 v2 в'- = - Pl (-^а& quot- !) + ~ (У' V-2 «2)'-:
о У
2 р!» Рз «2 Н-V (vv Рз 'Ог)'- = 0-
У
2 Pi Р2 м2 (^2 & quot-Ь ¦-9- Мо + р2 '-Оз [h2 -}-2- Мо u2J =
= -?Г (Г !^2 AJ)'- + (* - 1) М2 (_у& gt-2 и2 и-)'--
р, = 1/Л2- и2 = Л& quot-.
Начальные условия для уравнений (8) даются соотношениями (6).
При решении системы (8) с заданными *, Рг, «параметром будет число М0. Величина р, должна определиться в результате интегрирования уравнений (8) из условия равенства нулю продольной составляющей скорости на стенке. Как и для степенного класса, при получении решения в данном случае энтальпия газа на стенке и скорость отсоса (вдува) будут искомыми величинами. Условию для данного класса будет отвечать п = 0, т. е.
независимость коэффициента вязкости от температуры. Отметим, что автомодельные решения указанных классов при отсутствии массообмена на стенке были рассмотрены ранее в [3].
В качестве примера для х = 1,4, Рг = 0,71, п = 0,76 и числа М0 = 1 проведено численное интегрирование систем (5), (8) с начальными условиями (6).
Получение численного решения с заданной точностью проводилось методом проб (подбирались такие значения k или р, при которых величина yw, соответствующая и2 = 0, мало отличается от единицы) в сочетании с методом последовательных приближений.
В табл. 1 для решений степенного класса приведена сводка отдельных численных результатов.
Величина, А здесь определяет предельное значение координаты хпр (положительное или отрицательное), где газодинамические величины стремятся к нулю или к бесконечности
хпр = -/к.
В окрестности точки х = хпр используемые уравнения, вообще говоря, будут несправедливыми.
Из таблицы видно, что с увеличением & amp- величины /г2 (1) и г& gt-2 (1) монотонно убывают, причем может даже менять знак. Значе-
Таблица № 1
№ V к а1 А"0) «*(1)
1 0 -2,157 0 1,124 1,630
1 -4,729 0 1. 130 1. 315
0 -0,451 1 0,727 0,605
2 1 -0,991 1 0,759 0. 535
0 0,753 -1 0,409 0. 128
3 1 1,651 -1 0,453 0,128
0 1,353 -½ п 0,222 «?еО
4 1 2,959 -½ и 0,272 и2 = 0
ние отвечающее /г2(1) = 0, по физическому смыслу является максимально ВОЗМОЖНЫМ (?шах).
Зависимость параметра от к является немонотонной. На фиг. 1 для М0=1 эта зависимость приведена для плоского и осесимметричного случаев. Здесь соответствует значению ^ = 0,
II 3? 1 і
1
1 1 ! И
I
1 с 1

I- А/ / ! К- 1/1 —
и і
& gt-
У
— * и I
I й щ
2_ 1/ «т, ,

/
-/ / *
I /

1
. 1
Фиг. 1
к2- значению а, = - ½ п. Область I на графике отвечает течению с нарастанием продольной скорости (и температуры) и отводом массы по длине канала, область II — течению с уменьшением продольной скорости (и температуры) и отводом массы, область III — течению с уменьшением продольной скорости и подводом массы. С увеличением числа М0 область III может исчезать, а область II может быть ограничена значением? = & amp-шах ввиду достижения на стенке температуры абсолютного нуля.
Для значений М0, при которых ктах& gt-0, при к -& gt--+ 0 величина Я] -& gt-- + оо по закону & lt-з.х = сопэ^'-й. Последнее означает, что степенные решения при к — 0 вырождаются в экспоненциальные. Отсюда следует, что течения, соответствующие экспоненциальным решениям, возможны при отводе тепла от газа к стенке, т. е. при ^& lt-0,
Фиг. 2
'-
когда скорость и температура газа по длине убывают. В качестве примера на фиг. 2 представлены профили величин и2, V2, Л, в поперечном сечении канала для степенного класса при а^О. Видно, что приведенные профили газодинамических величин по сечению для плоского канала и цилиндрической трубы сравнительно мало отличаются друг от друга.
Случай (*! = () характерен тем, что продольная составляющая скорости, температура и величина №)т по длине постоянны, а давление убывает по линейному закону (в плоском и осесимметричном случае с разной интенсивностью, определяемой значением & amp-). При М0 = 1 этот режим течения имеет место практически при теплоизолированной стенке.
Отметим, что при М0 -* 0 рассматриваемые автомодельные течения газа переходят в течение Пуазейля с тождественно равной нулю поперечной составляющей скорости и постоянной температурой как по длине, так и по сечению. При этом для степенных решений & amp- ~ Мо, для экспоненциальных ^ - Мб.
2. Рассмотрим течения вязкого газа в плоских клиновых или конических каналах с малыми углами раствора (0да& lt-С1) — Для описания этих течений газа также будем использовать приближенную систему уравнений, которая в полярной системе координат г, 9 имеет вид
да
dp, г dr, dF
дг
(Р»)
Р"^гh +
дг
1
ди
гдт — 1
Re Ь2 г* ч» дг.
. ди
(р®) —
JJL + v-f- a+f) = °-
м’и2) +
[jV
h +
Рг Ref
1 д
гЦ'-1
dh
д^
(х — 1) М2
X — 1. «2 ,
-2- М. и2 1 й
Re в*
ди
рх М& quot- = рА, ^ = /г& quot-.
Здесь все величины безразмерные и связаны с размерными соотношениями, приведенными в п. 1, т| = 6? вт. Отличие состоит лишь в том, что радиальная координата отнесена к г%, поперечная составляющая скорости отнесена к Ие = р* и* а индекс «* и
соответствует величинам при 0 = 0 на расстоянии г* от начала координат. Кроме того, в этих соотношениях величина и* рассматривается как алгебраическая. Таким образом, Ие& gt-0 будет соответствовать течению в расширяющихся каналах, Ре& lt-0 — в сужающихся каналах.
Как и в предыдущем пункте, будем искать решения системы уравнений (9) вида
и = и1(г)и2(1Г1) — V =Vl®vt (r^)¦, А = А1(г)А2(-«}) —
Р=-
Р i® — р = р1(г)р,('-ч).
(10)
Подстановка выражений (10) в систему (9) дает следующие условия существования таких решений:
r2Plul rPlVl и j.
~--а — const- - 2п~ == b = const-
г/, = и, — hl = u\
w
2п+2
г2 Р1
и
2л+1
1
¦с = const,
где а, Ь, с — произвольные постоянные.
С учетом полученных условий и их совместности выражения (10) можно представить в виде
и — гт-ИзМ- V = гг"®2(•"]) — Л=г^Л2(т () — р гт" — р = ГЪр2(т)) — (11)
1)-1,
где Т2 = Т1, Тз = 2т» Т4 = Т1 (2л + 1) — 1, ?5 = Т1(2га
7! — свободный параметр.
Система обыкновенных дифференциальных уравнений для нахождения функций «2.2. А2, р2 при этом имеет вид
71 Ь «2 + Р2 ^"2=-
7, (2я + 1) — 1
(г1' Ъи'-оУ]
(V + 2 7, л) р2 «2 + - (V р2 г-,)'- = 0- ч
2 Т] р2 и2 (4-------------------------5----- Мо #2
ро г"2 (Л2
1, , … («- 1) М2
РгКев* V (71 2) + Неб?,
¦Мо Мг) =
1
¦ (V [*, иг и'-)'--,
р2 --- 1 /Но р*2 — Л2 •
Начальными данными для решения полученной системы уравнений (12) будут условия (6). При заданных теплофизических свойствах газа параметрами для данной системы будем считать ?! и М0-Выбор величины Иеб^ подчиним условию и2(1) = 0. Энтальпия газа на стенке и скорость отсоса (вдува), как и в предыдущих рассмотрениях, будут определены в результате решения задачи. Заметим, что при у 4~ 2-(1п = 0 имеем v~0.
Следует сказать, что данный класс решений для течения газа в каналах был рассмотрен ранее в работах [4, 5], причем для описания течения в них использовались полные уравнения Навье- Стокса. Однако результаты численных расчетов, приведенные в [5], отвечают только отдельным значениям 7,.
Для различных значений 75 и М0 = 1 (х= 1, 4, Рг = 0,71, я = 0,76) проведено численное интегрирование системы (12) при условиях (6). В табл. 2 приводится сводка некоторых численных результатов.
! Таблица 2
№ V Ке02№ Ь ?2(1) о. (1)
1 0 -0,664 1 0,425 -0. 524
1 — 1,350 1 0,470 -0,406
2 0 0,485 — 1 0,725 0,797
1 1,120 — 1 0,759 0,108
3 0 0,684 -½ п 0,766 0,546
] 1,606 -½ п 0,799 1*2 = 0
4 0 2,994 0 1,141 О III N & amp-
1 8,058 0 1,154 -0,242
Как было сказано выше, Ке6^& lt-0 отвечает течению в сужающихся каналах, Ке0^& gt-О- в расширяющихся каналах. Из табл. 2 видно, что с увеличением Ие6^ приведенная энтальпия газа на стенке Л2(1) возрастает, а отсос газа сменяется вдувом. Значение Яе6|, отвечающее Л2(1) = 0, является минимально — возможным. Функция от Иеб^ является немонотонной, причем 71-* + оо при + ° (фиг. 3).
В предельном случае при Яе 6^, = 0(0ц, = 0) автомодельное течение будет иметь место в плоском канале постоянного сечения или
Фиг. 3
цилиндрической трубе и, как показывает анализ, с экспоненциальным распределением параметров потока по длине. Такие течения рассмотрены выше в п. 1.
На фиг. 4 для 71 = 0 приведены графики величин и2(т), •г'-гС7]). Л,(т|). Случай 7! = 0 соответствует такому режиму течения в рас-
7
1,0
0,8
О, в
ол
0,2
Фиг. 4
ГгО

ч V 1
ч
к 1& gt-
г N II
г 1
1 ч ч 4
/ * А
V,
/ Ч
/ N
1 ч
/ к Ч к I1
/ N N 1
1 Л
/ 1
/ \
/
/ 1 1 1 V: н
/
/
О 0,2 0,4 0,6 0,8 1,0 агТи1уИг
3 Ученые записки ЦАГИ № 2 33
ширяющихся каналах, когда составляющие скорости и температура по длине постоянны, а давление убывает по гиперболическому закону. При этом для 7 = 0 стенка должна быть теплоизолированной и ъ = 0, для у=1 должен иметь место вдув газа практически при теплоизолированной стенке.
Предельным состоянием данных течений газа при М0 -+ 0 также является течение Пуазейля.
3. В предыдущем разделе показано, что автомодельные течения газа в каналах постоянного сечения с экспоненциальным распределением параметров потока по длине являются промежуточными для автомодельных течений в сужающихся и расширяющихся каналах с прямолинейными стенками. Для последних существуют автомодельные решения полных уравнений Навье — Стокса. Это указывает на существование для полных уравнений Навье-Стокса также экспоненциальных автомодельных решений, что находится в соответствии с результатами работы [6].
Система уравнений Навье-Стокса для случая течения газа в рассматриваемых цилиндрических каналах имеет вид
ди, ди др, 1 (4 д / ди, д / ди
дх ^ ду & lt-Заг Ре | 3 дх дх ] ду ду
. д / да & quot-2 д (
& quot-И"- ду (^ дх) 3 дх I ^
да
р идГ
ди
да дУ др
1 /ди, ди 2 д / V
у дх I 3 дх у) | '
1(4 д / ді' д (ди ,
йе 3 ду (^ ду) дх дх)
д / ди_ 2 д / да
& lt-3лг ^ у 3 ду дх
2а-
(да і'
?7» V) —
2 _д_ 3 ду
?-(Р»)+^-(Р*) + *Ру = 0.
ду
1
+
, 1 дк + ^-37
+
дх
(*-1)м I
Яе
ду I
Ие Рг дх
дА_
дх
9 V /-'-і2 (- 4- _ 2- (д-±. дл.
«^ у) ду) 3 ^ ду
рхМ^ - рН, [х = Лп.
(13)
Как и выше, здесь все величины безразмерные и связаны с размерными зависимостями, приведенными в п. 1. Однако в данном случае поперечная составляющая скорости отнесена к и*, а размерная координата х — к половине высоты или к радиусу канала.
Исходя из полученных выше результатов, интересующие нас решения будем искать в виде
и = и{ (х)и2 (у) = е^хи2(у), v — vl {х)ъ2 {у) — е'-^ху*(у), А = Аі {х)Н2 (у) = Л2 (у), р == -?у рх (х) р2 (у) =
1 X М& quot-
1
е (2п +)*р2(у), р = Рі (Ж) р2(У) = Є!
5, (2л- !)дг
92 (У),
(14)
где 8, — параметр, значение которого находится в результате численного решения задачи (см. ниже).
Из этих соотношений, в частности, следует, что число Ие, определенное по параметрам потока на оси и высоте канала, будет постоянным по длине канала..
Подставляя (14) в уравнения (13), получим систему обыкновенных дифференциальных уравнений для определения функций и2, Щ, К, р2, Р'2
51 92 А + р2 «2 ---°і(^2 Рі + (4 81 (2л + 1) ^ 11'-'- ~
+ (^2 и2У «Ь (н*2 УзУ з~ 81 (2 Л 4- 1) |12 V,-, +
Р-2 — («о + 81®г)~
--±- 8Д2"4- 1)^-3-
)
і Рг иг '-и2
Р, «г2 гг'----+ ?{_3_^зг, 2), + 5?(2п + +
-(2 /г -{- 1) ?а2 и2 3- Зі (ц-г иіУ + 4 2(*2у & quot-у"-) 3& quot-(^2]у) } '
2 §! про м, + (р2 г'-з)'- 4- г -у- = О,
2§! р2 // | /ї2 & quot-І- Р2 ^2
йе Рг
4 §і (п + 1) [а2 /г2 + (ц, А2)'-+& gt-р-2 ¦
(*-1)Мр
Ке
¦Р*2
+ 2 К)* + 2 V (Г + (8, «2 + «2)г — -І-І 81 и2 + гг2 +
^2
(15)
Р2 — Рг ^2& gt- 12 —2
Полученные обыкновенные дифференциальные уравнения (15) ¦будем решать для условий на оси канала
и2 (0) = /г, (0) = р, (0) = 1, к2 (0) = 0,
и2 (0) — К (0) = о, ®-(0)-?^т-
(16)
Второе условие для у2 получаем из третьего уравнения (15) с учетом симметрии на оси [р2(0) = 0]. Значение параметра определим в результате интегрирования системы (15) с начальными данными (16) при заданных х, Рг, п, числе Ие и М0 на оси канала из условия равенства нулю продольной скорости на стенке [и2(1) = 0]. Приведенная энтальпия газа на стенке канала к"{) и приведенная скорость отсоса ^2(1) будут при этом искомыми величинами. Как и в случае использования приближенных уравнений, результаты расчетов показывают, что имеет место ограничение режимов таких течений по числам М0 из-за достижения на стенке температуры абсолютного нуля. Предельно возможное число М0 отвечает Йе оо и нулевой температуре стенки.
Отметим асимптотические свойства полученных решений при больших числах Ие. Анализ и расчеты показывают, что при больших числах Ие (практически уже при Яе& gt-10) давление постоянно по сечению канала с точностью до членов 1/Не2 и имеют место
одни и те же профили величин uju0, (v/u0) Re, А/А0 по у. При этом Sj Re = pl5 где — параметр, введенный при использовании приближенной системы уравнений (1). Это означает, что экспоненциальные решения уравнений (1) являются предельными для данных решений уравнений Навье-Стокса при Re -* оо.
В качестве примера на фиг. 5 приведены результаты расчета для значений х=1,4, Рг = 0,71, и = 0,76, M0=l, Re = 5. Видно, что даже при столь малом числе Re давление по сечению практически постоянно.
При М0-«0 рассматриваемые течения газа переходят в течение Пуазейля.
ЛИТЕРАТУРА
1. Williams J. С. Viscous compressible and incompressible flow in slender channels. A1AA Journal, vol. 1, N J, 1963.
2. Б ы p к и h A. П., Межиров И. И. О расчете течения вязкого газа в канале. «Изв. АН СССР, МЖГ 1967, № 6.
3. Быркин А. П. Об автомодельных течениях вязкого газа в канале при наличии теплообмена. «Изв. АН СССР, МЖГ», 1969, № 5.
4. Щенников В. В. Об одном классе точных решений уравнений Навье-Стокса для случая сжимаемого теплопроводного газа.
ПММ, т. 33, № 3, 1969.
5. Б ы р к и н А. П. О точных решениях уравнений Навье -Стокса для течения сжимаемого газа в каналах.. Ученые записки ЦАГИ», т. I, № 6, 1970.
6. Shidlovski V. P. Special case of viscous gas motion in cylindrical tube in slip flowregime. Rarefied gas dynamics sixth symposium,
July, 1968. Academic Press, vol. 1. 1969.
Рукопись поступила 2/VII 1974 г¦

ПоказатьСвернуть
Заполнить форму текущей работой