CC BY
103
УДК 532.5326
ПРОЦЕССЫ ТОРМОЖЕНИЯ СВЕРХЗВУКОВЫХ ТЕЧЕНИЙ В КАНАЛАХ
И.И. Липатов
Центральный аэрогидродинамический институт им. проф. Н.Е. Жуковского, кафедра аэрогидромеханики E-mail: [email protected]
В работе исследованы нестационарные процессы взаимодействия и отрыва ламинарного пограничного слоя в сверхзвуковом потоке. Выведены уравнения, описывающие подобное течение, получены численные решения линеаризованной и нелинейной системы уравнений. Предполагается, что эти результаты позволяют объяснить эффекты, наблюдающиеся при торможении сверхзвуковых течений в каналах. Для полного описания процессов торможения необходимо рассмотрение течения в области присоединения, где, собственно, и возникают автоколебательные процессы.
ВВЕДЕНИЕ
Процессы торможения сверхзвуковых течений в каналах, несмотря на многолетние предшествующие исследования, все еще остаются во многих аспектах «terra incognita» газовой динамики. Причина этого, видимо, лежит в области практики. Пока для конструкторов эти явления не были важны или мало проявлялись в работе двигательных установок или других газодинамических устройств, достаточно было иметь некоторую эмпирическую информацию или располагать набором приближенных методов, предсказывающих интегральные характеристики процессов.
Ситуация определенно изменилась за последние 10-15 лет, когда на повестку дня встала разработка новых поколений авиационных и ракетных воздушно-реактивных двигателей, рассчитанных на работу при больших сверхзвуковых скоростях. Теперь потребовалось иметь более детальную информацию, не только интегральную, но и локальную о газодинамических процессах, сопровождающих торможение высокоскоростных потоков в каналах.
Эта информация была необходима также и потому, что потребовалось управлять процессами торможения. Одной из ключевых проблем стало создание модели — физической и/или математической, на основе которой можно было бы предлагать методы управления, определять границы устойчивой работы газодинамических устройств, да и вообще конструировать такие устройства.
Как писал Л.И. Седов — «исследовать явление — значит создать его модель... ». Можно было бы надеяться на современный уровень развития вычислительной гидродинамики и считать, что решение уравнений Навье - Стокса и будет служить такой моделью. Но и здесь остается достаточно вопросов, связанных с моделированием турбулентности, с нахождением критериев подобия, да и вообще с обоснованием (на физическом и математическом уровнях) постановок задач.
Представляется поэтому, что продвижение в понимании сложного явления может быть обеспечено за счет сочетания экспериментальных, вычислительных и теоретических методов исследования процессов торможения сверхзвуковых течений в каналах.
При исследовании дозвуковых течений обычно не возникает проблем. Входящий в канал поток содержит расположенные на стенках пограничные слои, толщина которых растет в продольном направлении, при достаточной длине канала они сливаются, формируя течение Пуазейля (ламинарное или турбулентное). При этом торможение течения у стенок требует разгона течения в невязком ядре (для выполнения условия сохранения расхода). При этом разгоне поперечный размер невязкого ядра уменьшается, что в сумме с увеличивающимися толщинами пограничных слоев обеспечивает выполнение сохранения толщины канала (если это канал с параллельными стенками. В любом случае течение такого типа, по определению Крокко, является докритическим, т.е. способным передавать возмущения вверх по потоку.
Processes of Supersonic Flows Retardation in Channels I.I. Lipatov
Investigated are unsteady laminar boundary layer interaction and separation processes in supersonic flow. Deduced are equations describing such flow and obtained are numerical solutions of linearysed and nonlinear system of equations. It is supposed that these results will allow explain effects accompanying supersonic flow retardation in channels. For the full description it is needed to investigate flow in the reattachment region where selfoscillating processes arise.
© И.И. Липатов, 2008
49
Совсем иная ситуация возникает при входе в канал сверхзвукового потока. Здесь также на стенках образуются пограничные слои. Но прежнего баланса между развитием пограничных слоев и невязкого сверхзвукового ядра уже нет. Режим реализующегося течения, по определению Крокко, уже может быть закритичским (сверхзвуковым в среднем) и переход в дозвуковое течение под влиянием каких либо причин на выходе (дросселирования, энерговыделения, выдува массы и т.д.), как показывают эксперименты, происходит через сложную структуру скачков и областей отрывного преисеночного течения, называемого в целом псевдоскачком. Интересно, что когда пограничный слой на стенках канала отсасывался, формирующийся в канале скачок по форме был близок к прямому.
Соответствующий пример можно видеть на рис. 1, приведенном в [1].
Рис. 1. Квазипрямой скачок уплотнения в канале после удаления пограничных слоев
Если же отсос пограничного слоя отсутствует реализуется совершенно другой режим течения с системой скачков, с областями отрыва и присоединения пограничного слоя.
То есть образование псевдоскачка целиком связано с проявлением сил вязкости при одновременном проявлении газодинамических эффектов невязких течений.
С пограничными слоями связано явление отрыва. В условиях внешнего обтекания известны соотношения, определяющие условия образования отрыва и параметры зон отрывных течений. При малых интенсивностях, возникающих в канале скачков, возможно безотрывное обтекание (для ламинарных и турбулентных режимов течения это разные интенсивности или разные допустимые числа Маха).
При увеличении числа Маха возникает ситуация, при которой возникающий отрыв (если его параметры такие же как и для отрыва во внешнем течении) должен был бы обязательно запирать течение в канале и формировать выбитую ударную волну. Вместе с тем эксперименты показывают, что этого не происходит, по крайней мере, в определенном диапазоне изменения параметров. Ясно, что при значительном повышении давления на выходе из канала режим с выбитой волной обязательно реализуется. Но возникает вопрос: что происходит в промежуточном диапазоне изменения параметров?
Фотографии структуры течения показывают, что, во-первых, трансформируется система скачков
— они приобретают лямбда-образную или икс-образную структуру (рис.2). При этом возникающие градиенты могут быть такими, при которых отрыв пограничного слоя еще не возникает. Понятно, что
Рис. 2. Торможение сверхзвукового потока в канале при М = 2 [1]
могут быть и другие режимы, при которых все-таки образуется отрыв. При этом, если бы он был таким, какой возникает во внешних течениях, это снова приводило бы к запиранию потока. Следовательно, в канале возникают (могут возникать) другие режимы отрывных течений. Главное их отличие состоит в меньших поперечных размерах (меньших углах отхода слоя смешения). Каким образом может реализоваться такая ситуация? Ответ, видимо, единственный — только в условиях возникновения нестационарного течения в пограничных слоях. Нестационарное течение возникает и в невязком ядре, но там влияние организованных пульсаций относительно мало, в отличие от течения в пограничных слоях.
Как модель можно было бы предположить возникновение Рэлеевского режима течения в пограничном слое. Для этого нужно было бы обеспечить колебательное движение стенок канала в продольном направлении и обеспечить определенные значения амплитуды и частоты колебаний. В силу отсутствия нелинейных конвективных членов в уравнениях пограничного слоя никакого отрыва с изменением поперечных масштабов пограничного слоя не возникнет, хотя и будут наблюдаться в течение определенных временных интервалов возвратные течения. Воздействие на такой пограничный слой ударных волн не приведет к возникновению отрыва (если выполнены определенные соотношения между интенсивностью волны и параметрами колебаний [3]).
Другой возможный вариант связан с возникновением пульсаций давления (при отсутствии колебаний стенки), накладываемых на некоторое среднее движение. Причиной появления таких пульсаций могут быть нестационарные процессы в области присоединения пограничного слоя, которые будут существенно менять структуру области возвратных (или колебательных) движений. Косвенно в пользу такого предположения свидетельствуют как результаты экспериментов, так и найденные в вычислительных экспериментах граничные условия, которые соответствуют не условиям прилипания на стенке, а условиям развития струйных течений, при которых на стенке ставятся условия равенства нулю нормальных производных.
Дальнейшее развитие происходит как последовательное развитие пограничных слоев до области смыкания. Ниже представлены результаты исследования течения в одной из областей псевдоскачка, связанной с началом области взаимодействия вязкого и невязкого течений. Цель этого исследования
— выяснить свойства нестационарного отрыва, проявляющегося в свободном потоке. Для течения в канале такое приближение применимо в начальном участке области взаимодействия.
Анализ экспериментальных данных показывает, что возникновение отрыва в канале отличается от соответствующего процесса в открытом течении. Среди различных гипотез наиболее вероятной представляется гипотеза о проявлении эффектов нестационарности, приводящих к повышению осред-ненного давления в точке отрыва, равно как и давления в области развитого отрыва.
1. ПОСТАНОВКА ЗАДАЧИ
Экспериментальные исследования, а в дальнейшем теоретический анализ, основанный на изучении асимптотических решениях уравнений Навье - Стокса [4-5], позволили выявить и описать такой режим течения, при котором в результате взаимодействия течения в пограничном слое с внешним невязким потоком индуцируется самосогласованное распределение давления. Для подобного режима течения в работе [6] предложен термин — «свободное взаимодействие».
К течениям со свободным взаимодействием относятся течения в локальных областях вблизи точек отрыва пограничного слоя, если отрыв начинается с гладкого участка контура обтекаемого тела, а взаимодействие пограничного слоя с внешним потоком до точки отрыва является слабым. Режим свободного взаимодействия реализуется также при взаимодействии слабых ударных волн с пограничным слоем, обтекании угловых точек контура тела, при сверхзвуковом и дозвуковом обтекании задней кромки пластины и во многих других случаях. Обзор некоторых задач, решенных в рамках теории «свободного взаимодействия», приведен в работе [7].
Исследования нестационарных течений в пограничном слое с самоиндуцированным давлением дали возможность отметить связь между теорией устойчивости и теорией «свободного взаимодействия» [8]. Результаты изучения пространственного течения с самоиндуцированным распределением давления представлены в работах [9-10], в которых рассмотрено обтекание малых неровностей, находящихся на дне ламинарного пограничного слоя в предположении, что возмущенное течение описы-
вается линейной системой уравнений. В настоящем разделе представлены результаты исследований нестационарных пространственных течений со свободным взаимодействием.
2. ОЦЕНКИ МАСШТАБОВ КООРДИНАТ И ФУНКЦИЙ ТЕЧЕНИЯ, УРАВНЕНИЯ И КРАЕВЫЕ УСЛОВИЯ
Рассматривается обтекание плоской поверхности потоком вязкого сжимаемого газа. Пусть на расстоянии I от передней кромки расположена область возмущенного ламинарного течения с самоин-дуцированным распределением давления. Декартова система координат выбрана таким образом, что направление оси ОХ совпадает с направлением скорости невозмущенного набегающего потока, величина которой равна ите, а ось О2 параллельна плоской поверхности. Предполагается, что число Рейнольдса достаточно велико, т.е. выполнен предельный переход
где р^ — плотность невозмущенного набегающего потока, д^ — коэффициент вязкости, соответствующий невозмущенному набегающему потоку.
Ниже для декартовых координат, компонентов вектора скорости, плотности, давления, энтальпии, коэффициента вязкости и времени приняты следующие обозначения: 1х, 1у, и^и, и^-и, р^р, р^и^р, и^Л, д<^д, 1£/и^ соответственно. Согласно теории «свободного взаимодействия», в возмущенной области течения можно выделить три характерные области. Размеры всех областей в направлениях х и г совпадают и равны по порядку величины е3/4(е = Б,е-1/2); следует отметить, что эти размеры больше, чем толщина пограничного слоя перед областью взаимодействия (<5д ~ е).
Размер области I в направлении, перпендикулярном поверхности, при конечных значениях числа Маха М совпадает по порядку величины с размерами в направлении осей ОХ и О2. В области I содержатся струйки тока сверхзвукового или дозвукового невязкого безвихревого слабовозмущенного течения. Пограничный слой на дне области I представляет собой характерную область II. Возмущенное течение в областях I и II оказывается в первом приближении невязким, поэтому в рассмотрение вводится область III с толщиной, равной е5/4 , лежащая на дне области II, где влияние сил вязкости становится сравнимым с силами инерции. Продольная скорость в пристеночной области меняется в главном члене, соответственно в основном порядке меняется толщина вытеснения, определяющая изменение возмущения давления.
В области I независимые координаты и функции можно представить в виде
х = 1+ е3/4 х1, у = е3/4у1, £ = е1/2 ¿4, и =1 + е1/2и1 (х1,у1) +
V = е1/2^1 (Х1 ,У1) + р = 1/тМ^ + е1/2р1 (Х1 ,У1) + р =1 + е1/2р1 (Х1,У1) + (2.2)
Л = 1/(7 - 1)М^ + е1/2Л1 (Х1, У1) + ...,
где М^ — число Маха в невозмущенном набегающем потоке, 7 — отношение удельных теплоемкостей.
Подстановка (2.2) в систему уравнений Навье - Стокса и предельный переход (2.1) приводят к уравнениям теории малых возмущений для сжимаемого потока:
Система уравнений (2.3) описывает квазистационарное возмущенное невязкое течение. Масштаб изменения функций по времени определяется по величине скорости и длине возмущенной области. Из системы уравнений (2.3) можно получить следующее уравнение для возмущения давления:
(2.1)
ди1 + др1 + о ди1 + др1 о + др1
дх1 дх1 ду1 дх1 дх1 дх1 ду1
(2.3)
(2.4)
Для дальнейшего анализа необходимо получить связь между распределением возмущения давления р1 и распределением вертикальной скорости V! при у1 ^ 0, т.е. на внешней границе пограничного слоя. Решение системы уравнений (2.4), устанавливающее подобную зависимость приведено в [11].
Решение уравнения (2.4), удовлетворяющее условию отсутствия приходящих и на поверхность извне возмущений имеет вид
Р1 = с/(х1 - в0У1), (2.5)
где во = (М2 - 1)1/2.
Следует отметить, что для течений в каналах в определенных условиях необходимо учитывать
приходящие (отраженные от противоположной стенки)возмущения. Исследование таких процессов
представлено в работах [12-14]. Наконец, искомое соотношение между возмущениями давления и вертикальной скорости для решения (2.5) имеет вид
Р1(х1, 0, ¿1) = Vl(xl, 0, ¿1 )в-1. (2.6)
В области II координаты и функции течения можно представить в следующей форме:
х = 1 + е3/4 х1? У = еУ2, £ = е1/2£1, и = и0 + е1/4и2 + ..., V = е1/2 v2 + ...,
р = 1/тМ2 + е1/2р2(х1, У1) + ..., р = ро + е1/4р2 + ... , Л = 1/(7 - 1)М2 + е1/4Л2 + ...
(2.7)
Предельный переход (2.1) после подстановки (2.7) в уравнения Навье - Стокса дает следующую систему уравнений:
ди2 , ди2
Ро и0—----+ Ро Ь2^— — 0,
дР2 п ди2 . дь2 . дР2 . дР2 п
— —0, Ро-------+ Ро--------------------------+ ио--+ Ь2— — 0, (2.8)
ду2 дхі ду2 дхі ду2
дР2 . дР2 п
ио ~----------------------+ — 0.
дх1 ду2
Ее решение можно записать в виде
( ^ ^и0 ^ А ^ дА(х1 >У2,¿1) ,от
и2(Х1 ,У2,¿1) — А(Ж1,У2,І1); Ь2(Х1, у2, ¿1) — -ио----------
Наконец, для масштабов координат и функций течения в области III верны формулы
х — 1+ є3/4 х15 у — є5/4 у2, Ь — є1/2£15 и — є1/4 из + ..., ь — єз/4ьз + ...,
р — 1/тД2 + є1/2рз (Х1 ,у1) + ..., Р — Рад + є1/2Рз + ..., Н — + є1/4 Нз + ...
(2.10)
Для главных членов разложений функций течения после подстановки (2.10) в уравнения Навье -Стокса получается следующая система уравнений:
ди3 ди3 ди3 др3 д2 и3
—--------+ рто и3 —----+ рто vз т:-+ т:— = Дш п 2 ,
д£1 дх1 дУ3 дх1 дУз (2 11)
др3 = 0; ди3 + = 0 .
^ 0; ^ ^ 0. дУ3 дх1 дУ3
Краевые условия определяются при сращивании решений в области III с решениями в областях II и I и с решением для невозмущенного пограничного слоя:
и3 ^ а(У3 + А), У3 ^ го,
и3 ^ аУ3, х3 ^ -го,
и3 = vз = 0, У3 = 0, (2.12)
Vl (х1, 0, ¿1) = - дАдх1,£1) = воР3, дх1
¿и0(0)
где а = —--------.
«У2
Известия Саратовского университета. 2008. Т.8. Сер. Математика. Механика. Информатика, вып.3 Введем замены переменных
Х3 = X(р„а5/2^/2)-1/2\И1 - 1|-3/8 , уз = У(р-2а-3^)1/4 \М£ - 1\"1/8 ,
рз = Р(а^)1/2 \М2 - 1\1/2 , ¿з = Т(р„а3/2^/2)-1 \М£ - 1\-1/4
^ = V(р-1 а3/2^/2)1/2 \М2 - 1\1/8 , и3 = и(р-2)1/4 \М^ - 1\-1/8 ,
^3 = W(р-2)1/4 \М2 - 1\-1/8 , А = А(р-2а-3^)1/4 \М£ - 1\"1/8 .
Преобразуем систему уравнений (2.11) и краевых условий (2.12) к виду
ди и ди уди дР _ д2 и
дт + дх + дУ + дХ _ дУ2 ’
дР _ ди дУ _
дУ _ 0 дХ + дУ _ 0
и ^ У + А, У ^го, и ^ У, X ^ -го,
и _ V _ 0, У _ 0,
(2.13)
(2.14)
V (X, 0,Т)_ - № _ Р.
Будем искать решение задачи (2.14) в виде малых возмущений к исходному решению, описывающему течение в пристеночной области пограничного слоя:
и _ У + С0/2(У)ехр(аХ + ¿^Т), Р _ СоР2 ехр(аХ + ¿^Т), (2.15)
где Со << 1
Замена переменных У _ а-1/3У1 приводит систему уравнений (2.14) для решений (2.15) к следующей системе уравнений:
(У1 + Т/ _ /2",
1 (2 16) /(0) _ /2(го) _ 0, /2'(0) _ а2/3Р2,
7 Р
где Т1 _ гша-1. Решение (2.16) выражается через функции Эйри: _ САг(У).
Наконец, дополнительное краевое условие, полученное при сращивании с решениями в областях I и II, позволяет получить дисперсионное соотношение:
Аг' (Т1) _ -«4/^ Аг(У1) ¿У.
Т1
3. РЕЗУЛЬТАТЫ РЕШЕНИЯ
На рис. 3 представлены зависимости действительной и мнимой части параметра а от действительного параметра <^. Полученные в таком диапазоне изменения параметров решения соответствуют бегущим внутри пограничного слоя волнам с изменяющейся по пространству амплитудой. Эти волны распространяются как вверх, так и вниз по потоку. В данном случае представлены решения, описывающие волны, распространяющиеся вверх по потоку. Физически эти результаты соответствуют распространению возмущений от источника, расположенного в некоторой точке поверхности.
ОО
Рис. 3. Зависимость действительной и мнимой части волнового числа от частоты колебаний
Кроме того, были получены решения и нелинейной краевой задачи, описывающей возмущенное течение выше по потоку от источника нестационарных возмущений давления. Это решение в определенном смысле моделирует решение полной задачи, содержащей и область нестационарного течения вблизи точки присоединения пограничного слоя. В расчетах в области возвратных токов использовалась аппроксимация Флюгге - Лотц, для которой конвективные члены при возникновении возвратных токов полагаются равными нулю. Как было ранее показано, такая аппроксимация удовлетворительно работает в течениях с локальными зонами отрыва.
На рис. 4 представлены зависимости от времени возмущения давления и поверхностное трение вблизи точки, где расположен источник нестационарных возмущений.
Рис. 4. Зависимость возмущения давления и поверхностного трения от времени
Библиографический список
1. Газовая динамика. Сер. Аэродинамика больших скоростей / Под ред. Эммонса. М.: Физматлит, 1957.
2. Гуськов О.Н., Копченов В.И., Липатов И.И. и др. Процессы торможения сверхзвуковых течений в каналах. М.: Физматлит. 2008.
3. Липатов И.И., Нейланд В.Я. К теории нестационарного отрыва и взаимодействия ламинарного пограничного слоя // Учен. записки ЦАГИ. 1987. Т. 18, № 1. С. 36-49.
4. Нейланд В.Я. К теории отрыва ламинарного пограничного слоя в сверхзвуковом потоке // Изв. АН СССР. Механика жидкости и газа. 1969. № 4. С. 53-57.
5. Stewartson K., Williams P.G. Self-induced separation // Proc. Roy. Soc. A. 1969. V. 312, № 1509. P. 181-206.
6. Chapman D.R., Kuehn D., Larson H. Investigation of separated flows in supersonic and subsonic streams with emphasis on the effect of transition // NACA Rep. 1958. № 1356. 40 p.
7. Гогиш Л.В., Нейланд В.Я., Степанов Г.Ю. Теория двумерных отрывных течений // Итоги науки. Гидромеханика. 1975. Т. 8. С.5-73.
8. Рыжов О.С., Терентьев Е.Д. О нестационарном по-
граничном слое с самоиндуцированным давлением // ПММ. 1977. Т. 41, вып. 6.
9. Smith F.T., Sykes R.I., Brighton P.W.M. A twodimensional boundary layer encountering a threedimensional hump // J. Fluid Mech. 1977. V. 83. Pt. 1. P.163-176
10. Липатов И.И. Пространственное обтекание малой неровности в режиме слабого гиперзвукового взаимодействия // Учен. записки ЦАГИ. 1980. Т. 11, № 2.
11. Черный Г.Г. Течение газа с большой сверхзвуковой скоростью. М.: Физматгиз. 1959.
12. Smith F.T. Upstream interactions in channel flows // J. Fluid Mech. 1977. V. 79. Pt. 4. P. 631-655.
13. Николаева Е.М., Тригуб В.Н. О самоиндуцирован-ном взаимодействии пограничных слоёв в плоском канале // Изв. РАН. Механика жидкости и газа. 1995. № 4. С. 131-141.
14. Дубинский С.В., Липатов И.И. О распространении возмущений в пограничных слоях во внутренних течениях // Изв. РАН. Механика жидкости и газа. 2004. № 2. С. 65-74.
УДК 539.376
НЕЛИНЕЙНОЕ ВЯЗКОУПРУГОЕ ПОВЕДЕНИЕ НАПОЛНЕННЫХ ЭЛАСТОМЕРНЫХ МАТЕРИАЛОВ
Е.В. Ломакин, Т.А. Белякова, Ю.П. Зезин*
Московский государственный университет им. М.В. Ломоносова, механико-математический факультет кафедра теории пластичности,
^Институт механики МГУ им. М.В. Ломоносова E-mail: [email protected]
Проведены экспериментальные исследования процесса релаксации напряжений для резины плотностью 1200 кг/м3, включающие испытания на релаксацию при нормальной температуре в условиях одноосного растяжения и сжатия при различном уровне деформаций. Для описания вязкогиперупругого деформирования наполненных эластомерных материалов используются определяющие соотношения, являющиеся обобщением нелинейной теории упругости и линейной теории вязкоупругости Больцмана - Вольтерра. Разработан метод определения параметров взаимосвязанных гиперупругих и реологических характеристик деформирования. На основании экспериментальных данных отмечена зависимость релаксационных свойств исследуемых эластомеров от уровня деформации и предложен подход к описанию этой зависимости, основанный на применении определяющих соотношений нелинейной эндохронной теории стареющих вязкоупругих материалов. Наблюдается вполне удовлетворительное соответствие теоретических и экспериментальных зависимостей.
Non-Linear Viscoelastic Behavior of Filled Elastomers
E.V. Lomakin, T.A. Beliakova, Yu.P. Zezin
The analysis of the results of experimental studies of relaxation properties of rubber at room temperature under the conditions of tension and compression at different strain levels is carried out. For the description of viscohyperelastic deformation of filled elastomers, the constitutive relations are used that are the generalization of non-linear elasticity and linear viscoelasticity Boltzmann - Volterra. The method for the determination of interrelated hyperelastic and reological characteristics of deformation is proposed. On the base of experimental data, the dependence of relaxation properties of elastomer on the strain level is discovered and the approach to the description of this dependence based on the constitutive relations of non-linear endochronic theory of aging viscoelastic materials is proposed. The finite element numerical simulation of an experiment under conditions of the confined compression of flat cylindrical specimen and subsequent relaxation is performed. Sufficient correspondence between the calculated values of axial load and the results of experimental study is demonstrated.
ВВЕДЕНИЕ
Возрастающее использование резиноподобных материалов во многих областях современной техники приводит к необходимости описания с высокой точностью кратковременных и длительных характеристик деформирования и разрушения эластомерных элементов конструкций. Повсеместное использование эластомеров в нефтегазовой промышленности, машиностроении, гражданском строительстве,
© Е.В. Ломакин, Т.А. Белякова, Ю.П. Зезин, 2008
