О способах повышения точности расчета ударных волн1 Текст научной статьи по специальности «Математика»

Математика

УДК 519.63

О СПОСОБАХ ПОВЫШЕНИЯ ТОЧНОСТИ РАСЧЕТА УДАРНЫХ ВОЛН1

Д. А. Варфоломеев2, К. И. Казин3, В.Ф. Куропатенко4

Рассматривается неоднородный разностный метод расчета ударных волн в Лагранжевых координатах. Метод позволяет явно выделять в решении ударные волны в виде разрывов первого рода. Предлагаются способы повышения точности расчета выделенных ударных волн в рамках этого метода. В частности, для определения величин перед фронтом ударной волны наряду с разностным подходом предлагается использовать элементы метода характеристик. Приводится алгоритм модифицированного метода. На примере расчета методических задач показано, что применение модифицированного метода позволило повысить монотонность и точность расчета ударных волн.

Ключевые слова: ударная волна; сильный разрыв; неоднородный разностный метод.

Введение

Для математического моделирования ударных волн в сплошных средах широко применяются однородные разностные методы, в которых ударные волны «размазываются» на несколько интервалов сетки. Неоднородные разностные методы позволяют явно выделять в решении ударные волны в виде разрывов первого рода, которые разграничивают области с гладким решением. Неоднородные методы являются алгоритмически более сложными, но при этом дают возможность на одинаковой сетке получать результаты точнее, чем результаты, полученные с применением однородных методов. Важной особенностью методов с выделением разрывов является отсутствие энтропийных следов в решении.

Одним из первых неоднородных методов, который позволяет выделять в решении сильные и слабые разрывы, был метод характеристик [1, 2]. Метод характеристик использует Эйлеровы координаты.

В 60-х годах создан неоднородный разностный метод в Лагранжевых координатах [3]. Одной из особенностей данного метода является то, что он способен выделять в решении наиболее «важные» разрывы, а остальные разрывы «размазывать». Такая особенность позволяет с успехом применять его для расчета очень сложных течений, акцентируя внимание на определяющих течение областях. Разрывы во второстепенных областях размазываются, что позволяет сократить вычислительные затраты и при этом сохранить точность вычислений.

В данной работе для метода [3] предлагаются способы повышения точности расчета выделенных ударных волн. Для определения величин перед фронтом ударной волны наряду с разностным подходом применяются элементы метода характеристик. Для повышения точности расчета величин за фронтом ударной волны при взаимодействии регулярной сетки и сетки особенностей предлагается алгоритм, позволяющий корректно учитывать возмущения, догоняющие фронт. Для интервалов перед фронтом и за фронтом предлагается подход, который сводит расчет величин в данных интервалах к расчету величин на регулярной сетке.

1 Работа выполнена при поддержке РФФИ ГРАНТ №13-01-00072

2 Варфоломеев Денис Александрович - начальник сектора, Российский Федеральный Ядерный Центр - Всероссийский Научно-Исследовательский Институт Технической Физики им. академика Е.И. Забабахина, г. Снежинск.

E-mail: [email protected]

3 Казин Кирилл Иванович - младший научный сотрудник, Российский Федеральный Ядерный Центр - Всероссийский Научно-Исследовательский Институт Технической Физики им. академика Е.И. Забабахина, г. Снежинск.

E-mail: [email protected]

4 Куропатенко Валентин Федорович - доктор физико-математических наук, профессор, главный научный сотрудник, Российский Федеральный Ядерный Центр - Всероссийский Научно-Исследовательский Институт Технической Физики им. академика Е.И. Заба-бахина, г. Снежинск.

E-mail: [email protected]

Рис. 1. Шаблон расчета особого интервала

Шаблон для определения сеточных величин и основные этапы расчета

В данной работе за основу взят неоднородный разностный метод В.Ф. Куропатенко [3]. В этом методе используется разнесенный шаблон - координаты и скорости определены на границах интервалов, а термодинамические величины в их серединах. Будем придерживаться терминологии, используемой в [3]. Интервалы, в которых отсутствуют разрывы, называются регулярными. Интервалы, в которых присутствуют разрывы, называются особыми. Если разрыв попадает в регулярный интервал, то этот интервал становится особым. Когда разрыв покидает данный интервал, он снова становится регулярным. Таким образом, выделенные сильные разрывы рассчитываются как особенности, перемещающиеся по регулярной сетке.

Рассмотрим одиночный сильный разрыв, перемещающийся вправо (рис. 1). Для наглядности здесь и далее на подобных рисунках будем использовать Лагранжевы координаты, в которых частицы не перемещаются, а, следовательно, траектории границ интервалов являются вертикальными линиями. Точка I является фронтом ударной волны и принадлежит сетке особенностей. Пусть за время г разрыв переходит с п на (п + 1) слой, при этом остается внутри интервала и не взаимодействует с другими разрывами и контактными границами.

Особый интервал состоит из двух интервалов: интервал перед фронтом - «а-интервал», интервал за фронтом «^-интервал». Для величин на разрыве используются обозначения: «-» - перед фронтом, «+» - за фронтом.

Перед началом расчета интервала с разрывом известны все величины в регулярных интервалах на момент времени 1п и В интервале с разрывом на момент известны:

• (р, Е, Р,и, С)-, (р, Е, Р,и, С)+ - величины на разрыве; , Вп - положение и скорости разрыва;

• (р,Е,Р,М)"а - величины в интервале перед фронтом (а-величины);

• (р,Е,Р,М)в - величины в интервале за фронтом (^-величины).

Здесь р - плотность, Е - удельная внутренняя энергия, Р - давление, и - скорость, С - скорость звука, М - масса интервала, В - скорость перемещения разрыва в Эйлеровых координатах, Ж - скорость перемещения разрыва в Лагранжевых координатах.

Задача состоит в том, чтобы на момент 1п+х определить все величины в интервале с разрывом и величины на разрыве. Алгоритм расчета особого интервала можно разбить на несколько этапов. Не вдаваясь в конкретные операции каждого этапа, рассмотрим общую структуру. Последовательно определяем:

Величины перед фронтом:

• новое положение фронта и величины перед разрывом (р,Е,Р,и, С)-+1;

• плотность р^1 в а-интервале и массу АМ, которая «перетекла» сквозь фронт за время г;

• термодинамические величины в (Е, Р, М )^+1 а-интервале. Величины за фронтом

• величины за разрывом (р, Е, Р, и, С )++1;

• величины в ^-интервале и его массу Мр

После проведения данных операций можно говорить, что на момент времени ^ с разрывом определены все величины.

,п+1

,п+1

в интервале

Расчет минус величин

Определим скорость фронта ударной волны ""+1 = Ю" +т

СЮ

Л

"+1

где

СЮ

Л

"+1

ускоре-

ние фронта, полученное с применением дополнительного к условиям Гюгонио уравнения для расчета величин вдоль траектории фронта ударной волны [3]. Новое положение фронта

я"1 = Я" +тЮ"+1 .

Рассмотрим точку Б (рис. 2). Величины в данной точке характеризуют состояние на разрыве перед фронтом в момент ?и+1 и соответствуют «минус» величинам. Координата точки Б совпадает с положением фронта в момент {'+х,

Для нахождения решения в точке Б воспользуемся характеристическими свойствами уравнений газовой динамики [1]. Построим характеристический треугольник. Для этого из точки А проведем линию с наклоном

как у а характеристики СЯ = (и + с )"А, а из точки В с наклоном как у /? характеристики = (и - С)",. Для

Ся

точки Б построим траекторию — = и^ . Координаты и скорости точек А, В и Б неизвестны и их

Ж

необходимо определить. Для этого составим систему из трех уравнений. Первое уравнение - это уравнение а характеристики, второе и третье - выражения для и и С, полученные путем линейной интерполяции на п слое. Запишем систему уравнений для точки А

Я"1 = ЯПА +т(и + с)" ,

Рис. 2. Определение «минус» величин

иА = ип +

сА = с"" +

ип п им р" р"

ЯЮ " Ям

СП у—

- " см ЯЮ - Ям

(ЯА" - Я"), (ЯА - Я").

Решив данную систему, получим выражение для определения координаты точки А

Я

Яа =■

п+1

-т(и-" + с- - Я" (21 + 22))

1 + т

(21 + 22)

где 21 =

и- - им

и" и"

ЯЮ - Ям

2 2 =

с" -с

м

ЯЮ - Ям

Для определения Яв решим аналогичную систему при этом первое уравнение примет вид

Я"+1 = ЯВ +т(и - с)" . Получим

Яв =-

Я-1 -т(у- - с- - Я" (21 - 22 )) 1 + т(2, - 2

(21 - 22 )

Далее линейной интерполяцией на " слое рассчитываем величины (р,Е)АВ и из уравнения

состояния находим РАВ и сАВ .

Запишем уравнения, описывающие изменения величин вдоль а и Р характеристик с учетом симметрии (а = 1, 2, 3)

1 (а- 1)исЖ

сСи ±-ЖР -'--= 0.

рс Я

"

Для каждой из характеристик заменим дифференциальное уравнение разностным и после несложных преобразований получим выражения для определения и-+1 и Р_п+1

U-+1 =.

фл - ФВ

рлСл +pbcb

рп+1 = РлслФв -РвсвФл Paca +Pbcb

Здесь фл = рл + PaCaUл

(a- 1)PaUaCA _ тт (a- 1)PbUbCB

-т, фв = рл + PbCbUв

А = РА + РАСАиА 1 , ^В = РА

ЯА ЯВ

Из точки Б опустим на п слой траекторию и найдем положение точки Б

ЯР = Я1+1 -тип*1.

т.

Величины (р,Е,Р, С)р определяются аналогично точкам А и В. Вдоль траектории энтропия постоянна, следовательно С2 = —Р. Для определения плотности перед фронтом запишем разност-

d р

ныи аналог этого уравнения

рИ+1 р

и+1 . р- - PF

Р- =Pf +-

C 2 CF

+1

Для определения Е- проинтегрируем с заданной точностью уравнение изэнтропы —Е = -Р—У вдоль траектории частицы в пределах (уП, У-+1). Из уравнения состояния найдем

давление Р-+1 и скорость звука С-+1. Таким образом, определено состояние перед фронтом на

,п+1

момент t .

В целях повышения точности проведем пересчет. Для этого из точки Б «опустим» характеристики и определим ЯА, ЯВ следующим образом

ЯА = Яд -т(и + С)_ 1, ЯВ = Яд - т(и - С)-Затем пересчитаем «минус» величины для данных координат.

Расчет плотности в а-интервале и массы, перешедшей через фронт

Для учета изменения плотности в а-интервале за шаг по времени представим, что на левой

границе интервала (япв , Я!"1) с массой М^ отсутствует разрыв, т.е. он является регулярным, его

масса не меняется, левая граница движется со скоростью вещества перед фронтом (рис. 3). Обозначим состояние данного интервала на момент /' 1 индексом С.

где

ч П+1

Координата левоИ границы RAC = RD +tUACc ■

r-rn+1 rrn . _ UAC = U- + T

dU | n

—— J . Правая граница совпадает с границей соседнего регулярного интервала. На момент времени tn+1 определим объем и плотность

(RAC1, RnB+C Не1, РС+1 = Mil &t

n +1 C

Рис. 3. Расчет «альфа» величин

Пусть интервал К характеризует ближайший к разрыву регулярный интервал перед фронтом на момент {п+х. Тогда плотность в а-интервале на мо-

,и+1

мент времени t определим линеинои интерполяциеи

Р1+1

=рС+1+Р+1 -РС+1)

n +1 n +1 rd - rac

2R

n +1

к

-(RBC + Rn+1

AC

)

Далее определим объем а-интервала (RD+1, RBC его массу MI+1 = р

?na+1Bna+1 и массу, пе-

решедшую за время т сквозь фронт ударноИ волны AM = M И-Ma

n+1

Расчет величин в а-интервале

Предположим теперь, что интервал (RF,Щ+1) является регулярным (см. рис. 3). Точка Ь -

середина интервала. Линейной интерполяцией рассчитываем значения термодинамических величин в этом интервале на момент времени На момент времени состояние в данном интервале соответствует значениям величин в а-интервале. Таким образом, задача определения альфа

величин сводится к расчету регулярного интервала (Щр, Щ+1) на момент времени {г+х. Знак вели-

чины Аи = и-г+1 - иВС определяет сжатие или растяжение в интервале. Пользуясь методом расчета регулярных интервалов [3], определяем величины в интервале перед фронтом ударной волны. Таким образом, при расчете регулярных интервалов и интервала перед фронтом используется один метод, что позволяет получить одинаковый порядок аппроксимации.

тр+1

Расчет плюс величин

На фронте ударной волны значения кинематических и термодинамических величин справа и слева в идеальной среде связаны между собой условиями Гюгонио

Ж (Ц+- и-)-( Р+- Р-) = 0, Ж (¥+- V-) + (U+- и-) = 0, (Е+- Е-) + 0,5 (Р+ + Р„){У+- У-) = 0.

Величина скорости фронта ударной волны Вп+ и состояние перед фронтом на момент времени определены. Найдем значение Жп+1

Жп+1 = Р-+1 (Вп+1 - Ц-+1).

Тогда условия Гюгонио в совокупности с уравнением состояния Р = Р(V,Е) содержат четыре

уравнения и четыре неизвестных. Разрешая полученную систему уравнений, определим параметры на фронте ударной волны в момент времени .

Расчет величин в ^-интервале

Положение правой границы ^-интервала соответствует положению сильного разрыва. Левая граница является границей соседнего регулярного интервала (рис. 4). Масса в интервале за фронтом равна сумме массы на момент времени {г и массы, перешедшей через фронт ударной волны за шаг г.

МП+1 = М "р+АМ. Рассчитаем объем

(Щ/-+1, Щр+1 и определим плотность

— П+1 д *П+1 I /}П+1

Рр =Мр / вр .

Для определения величин в интервале за фронтом поступим аналогично определению величин в а-интервале. Представим, что интервал (Щ^1, ) регулярный, т.е. разрыв на правой границе отсутствует. Определим состояние в данном интервале на момент времени 1;п. Из точки с = 0,5 (ЩР+1 + Щ^1), «опустим» со скоростью ир+1 = 0,5 (ип+ + Ц++1) траекто-

Рис. 4. Расчет величин в интервале за фронтом

координатой Щр+1

рию на п слой и получим точку Ь: Щь = Щр -ти" . Зная распределение величин за фронтом на

момент Г, рассчитаем величины в точке Ь. Теперь ^-интервал на момент времени можно рассчитать как регулярный. Заметим, что подобные операции проводятся потому, что масса интервала за шаг по времени изменилась, а, следовательно, середина ^-интервала на моменты времени Г и не принадлежат одной траектории. По этой причине нельзя использовать Ла-гранжев метод напрямую.

После определения величин в интервале за фронтом можно говорить, что особый интервал

,п+1

на момент времени t рассчитан полностью.

Особый случай расчета величин в ^-интервале

В процессе счета сетка особенностей перемещается относительно регулярной сетки. При этом возникают ситуации, когда разрыв переходит через границу регулярного интервала. В таком

случае на момент масса в интервале за фронтом равна нулю, т.е. М^ = 0 (рис. 5) и применение

подхода, описанного выше, приведет к тому, что точка Ь будет лежать на траектории фронта ударной волны. На практике это приводит к тому, что в профиле величин за фронтом возникают немонотонности. Они являются следствием того, что при расчете ^-интервала некорректно учитываются возмущения, «пришедшие» из области за фронтом. Решить данную проблему можно путем измельчения шага по времени, что негативно отразится на расчете всей системы. В рамках рассматриваемого метода был разработан алгоритм расчета величин в интервале за фронтом при переходе разрыва через границу регулярного интервала.

Временно исключим из рассмотрения регулярную точку (/-1). Таким образом, эту ситуацию сводим к расчету величин в интервале за фронтом в общем случае. Рассчитаем величины в интервале К, т.е. интервале

). Масса этого интервала равна:

(К+\ Я

Мк = МП_15 + АМ . Определяем объем интервала К и его плотность

(Г) п+1 Г) п+1 * Яi_2, Я

Ю

)^вк , Рк =

М

Рис. 5. Расчет Д-интервала при выходе разрыва из целой точки

к

в

к

Я

к

Из точки К опустим траекторию на п слой и находим: = 0,5 (Я^1 + Яп+1

Яь = Як _тик,

где

П+2 ), ик = 0,5 (и++1 + иП_+2). Вт. Ь найдем термодинамические величины, за-

чП+1

тем, аналогично расчету величин в ^-интервале, описанному выше, определим (Е, Р)к .

На момент времени масса в ^-интервале соответствует массе перешедшей через фронт за шаг х: Мв+1 =АМ. Предположим, что в области (Як, Я) распределение плотности и энергии

линейные. Найдем такое значение ЯП+, при котором выполнятся условия

' „п+1 АМ Рв ~

в

п+1 в

в =Рк +(Рк -Р+П+1 )

О п+1

р++П Япп

Я

к

Як _ Я

п+1 £

Для этого итерационно вычислим ноль функции

АМ

*(ЯД1 ) = АПг_Рк _(Рк-Р++1)

Определив Я^1, рассчитаем иП_+11 и

вв ' сЮ_ у А

+1

Яв _ Як

„ „п+1

Як _ Яй

+1

/ _1

Г) п+1 пп

ттП+1 _ Я/_1 _ КР и1 _1 =

Т

жи

А )■

+1

ип_+1 _ и+п

/_1

Энергию вычислим линейной интерполяцией на (п+1) слое между «плюс» величинами и величинами в точке К

Е}++ = Ек + (Ек

(Ек _ Е++1)

+1

Яв _ Як

„ „п+1

Як _ Яо

Давление определяется из уравнения состояния.

Т

Применение данного подхода позволяет корректно учитывать возмущения, догоняющие фронт, тем самым, повысить точность и монотонность профилей величин за фронтом.

Тестовые расчеты

Задача о движении поршня с постоянным давлением на границе по линейному профилю плотности.

Задана область 0 < Щ < 1. На левой границе давление = 1,5, на правой границе и в области и = 0. Количество расчетных точек равно 10. В области: энергия Е = 0, распределение плотности р(Щ) = 1,7 - 1,4Щ . Вещество - идеальный газ, у = 5/3. Симметрия плоская.

На левой границе системы после распада разрыва образуется ударная волна, которая распространяется вправо. На рис. 6 представлено распределение плотности в системе на момент времени t = 0,49, полученное с использованием модифицированного метода. До прихода ударной волны профиль плотности перед фронтом не меняется, т.к. давление в области равно нулю. За фронтом наблюдаем сжатие вещества.

7.ОС «.СО

3.00

1.00

Э.00

2.00 1.00

4.М 1.00

Рис. 6. Профиль плотности на момент времени (= 0,49

На рис. 7 представлены зависимости плотности р_ перед фронтом ударной волны от времени, полученные с помощью метода [3] и модифицированного метода, а также аналитическое решение (начальное распределение).

Рис. 7. Зависимости плотности перед фронтом от времени. Метод [3] (—), модифицированный метод (—), точное решение (•)

Рис. 8. Зависимости плотности за фронтом от времени.

Метод [3] (—), модифицированный метод (—)

Из рис. 7 видно, что использование метода, основанного на применении уравнений в характеристической форме для определения величин перед разрывом, позволило повысить монотонность и точность расчета термодинамических величин и скорости. Все неточности, возникающие при определении состояния перед фронтом, переносятся сквозь фронт и отражаются на величинах за фронтом, что видно из рис. 8.

Задача о движении поршня с переменной скоростью по постоянному фону.

Задана область 0 < Я < 1. На правой границе и в области и = 0, на левой границе скорость

изменяется со временем по закону иь ^) = 1 -12 . Количество расчетных точек равно 10. В области: энергия Е = 0,5, плотность р = 1,7. Вещество - идеальный газ, у = 5/3. Симметрия плоская.

От левой границы в систему пойдет ударная волна, за фронтом которой вещество разгружается. На рис. 9 представлено распределение давления в системе на момент времени I = 0,2557.

3.3

з.о

2.5

ь 2.0

1.5

1.0 0.5

0.00 0.20 0.40 0.60 0.80 1.00

Рис. 9. Распределение давления в системе на момент времени (= 0,2557

3. 35 3. 30

3.25

з.го

0.33 0.35 0.37 0.3® 0.41 0.03

Рис. 10. Распределение давления в области за фронтом на момент времени (= 0,2557.

Метод [3] (-•-), модифицированный метод (—), «точное» решение (•)

На рис. 10 представлено распределение давления в области за фронтом на момент времени I = 0,2557. Этот момент времени соответствует моменту перехода разрыва через границу регулярного интервала. Здесь под «точным» решением подразумевается решение, полученное на 1000 точках. Дальнейшее увеличение числа точек при рассмотрении в данном масштабе бессмысленно. Из рис. 10 видно, что профили величин за фронтом ударной волны, полученные модифицированным методом, более монотонные, значения величин на ударной волне ближе к «точному» решению.

Заключение

Предлагаемые способы модификации метода [3] повысили точность и монотонность расчета выделенных ударных волн. Идеология и порядок расчета величин остались прежними, изменились методы их определения. Авторами показана принципиальная возможность использования элементов метода характеристик внутри разностного метода. Данный подход применительно к методу [3] показал свою работоспособность и эффективность. Для любого другого разностного метода элементы метода характеристик можно использовать в качестве интерполяции для определения величин в любой точке пространства и времени. Преимуществом такой интерполяции является учет особенностей уравнения состояния и изменения течения во времени.

Литература

1. Жуков, А.И. Применение метода характеристик к численному решению одномерных задач газовой динамики / А.И. Жуков // Труды Математического института АН СССР, 1960.

2. Хоскин, Н.Э. Метод характеристик для решения уравнений одномерного неустановившегося течения / Н.Э. Хоскин // Вычислительные методы в гидродинамике: сб. науч. тр. - М.: Изд-во Мир, 1967. - С. 264-292.

Варфоломеев Д.А., Казин К.И., О способах повышения точности расчета ударных волн

Куропатенко В.Ф.

3. Комплекс программ «ВОЛНА» и неоднородный разностный метод для расчета неустановившихся движений сжимаемых сплошных сред. Ч. 1. Неоднородный разностный метод / В.Ф. Куропатенко, Г.В. Коваленко, В.И. Кузнецова и др. // ВАНТ. Серия «Математическое моделирование физических процессов». - 1989. - Вып. 2. - С. 9-17.

Поступила е редакцию 3 июня 2014 г.

Bulletin of the South Ural State University Series "Mathematics. Mechanics. Physics" _2014, vol. 6, no. 3, pp. 5-13

ON METHODS FOR INCREASING THE ACCURACY OF SHOCK WAVE CALCULATIONS

D.A. Farfolomeyev*, K.I. Kazin2, V.F. Kuropanenko3

The paper presents a non-homogeneous difference method for shock wave simulation in Lagrangian coordinates. The method explicitly determines shocks as discontinuities of the first kind. The ways of increasing the accuracy of calculated shocks within the scope of the proposed method are suggested. In particular, the author suggests using the difference approach together with some elements of the method of characteristics to determine quantities before the shock front. The algorithm of the modified method is described. A number of methodical calculations are provided. On their basis it's shown that the new method improves the monotonicity and accuracy of shock calculations.

Keywords: shock wave; strong discontinuity; non-homogeneous difference method.

References

1. Zhukov A.I. Trudy Matematicheskogo Instituta ANSSSR.. 1960. (in Russ.).

2. Khoshin N.E. Metod kharakteristik dlya resheniya uravneniy odnomernogo neustanovivshegosya techeniya (A method of characteristics for solving the equations of 1D transient flows). Vychislitel'nye metody v gidrodinamike (Fundamental Methods in Hydrodynamics). Moscow, MIR Publ., 1967. pp. 264-292. (in Russ.). [Hoskin N.E. Solution by Characteristics of the Equation of One-Dimentional Unsteady Flow. Fundamental Methods in Hydrodynamics. New York, London, Academic Press, 1964. pp. 265-293. (in Eng.).]

3. Kuropatenko V.F., Kovalenko G.V., Kuznetsova V.I. et al. VANT. Seriya "Matematicheskoe modelirovanie fizicheskikh protsessov". 1989. Issue 2. pp. 9-17. (in Russ.).

Received 3 June 2014

1 Varfolomeev Denis Aleksandrovich is Head of the sector, Russian Federal Nuclear Center - Zababakhin All-Russian Scientific Research Institute of Technical Physics, Snezhinsk.

E-mail: [email protected]

2 Kazin Kirill Ivanovich is Junior Research Fellow, Russian Federal Nuclear Center - Zababakhin All-Russian Scientific Research Institute of Technical Physics, Snezhinsk.

E-mail: [email protected]

3 Kuropatenko Valentin Fedorovich is Dr. Sc. (Physics and Mathematics), Professor, Senior Researcher, Russian Federal Nuclear Center -Zababakhin All-Russian Scientific Research Institute of Technical Physics, Snezhinsk.

E-mail: [email protected]