Чисельне моделювання аварійних вибухів рудникової атмосфери
Мета. Розробка ефективної схеми чисельного розрахунку спільного розв'язку задачі газової динаміки й хімічної кінетики горіння газоповітряного середовища на основі методу великих часток. Методика. Математичне моделювання, чисельний експеримент, аналіз і узагальнення й результатів. Результати...
Збережено в:
| Дата: | 2020 |
|---|---|
| Автори: | , |
| Формат: | Стаття |
| Мова: | Ukrainian |
| Опубліковано: |
Інститут фізики гірничих процесів НАН України
2020
|
| Назва видання: | Физико-технические проблемы горного производства |
| Теми: | |
| Онлайн доступ: | http://dspace.nbuv.gov.ua/handle/123456789/166520 |
| Теги: |
Додати тег
Немає тегів, Будьте першим, хто поставить тег для цього запису!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Цитувати: | Чисельне моделювання аварійних вибухів рудникової атмосфери / М.М. Налисько, Л.І. Барташевська // Физико-технические проблемы горного производства: Сб. научн. тр. — 2020. — Вип. 22. — С. 85-102. — Бібліогр.: 12 назв. — укр. |
Репозитарії
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
irk-123456789-166520 |
|---|---|
| record_format |
dspace |
| spelling |
irk-123456789-1665202020-02-26T01:25:55Z Чисельне моделювання аварійних вибухів рудникової атмосфери Налисько, М.М. Барташевська, Л.І. Физика горных процессов на больших глубинах Мета. Розробка ефективної схеми чисельного розрахунку спільного розв'язку задачі газової динаміки й хімічної кінетики горіння газоповітряного середовища на основі методу великих часток. Методика. Математичне моделювання, чисельний експеримент, аналіз і узагальнення й результатів. Результати. Для спільного рішення задачі газової динаміки і хімічної кінетики горіння газоповітряної суміші пропонується ввести в чисельну схему методу великих часток концентраційну функцію, яка дозволяє враховувати багатокомпонентний склад газового середовища. Концентраційна функція дає можливість вводити в чисельну схему рівняння хімічної кінетики в вигляді рівняння Арреніуса і розрізняти компоненти хімічної реакції і продукти горіння. У задачі розрахунку детонаційних вибухів виникають сильні градієнти тисків, які, при виході фронту ударної хвилі на границю вільний вихід генерують нефізичні флуктуації параметра. Для виключення їх впливу проводиться аналіз різних видів апроксимації параметрів в фіктивний шар розрахункової схеми. З аналізу фізичних процесів знайдений ефективний вид граничних умов вільний вихід для задачі поширення ударної хвилі в каналі. Наукова новизна. Модифікація чисельного методу великих часток за рахунок введення концентраційної функції дозволяє отримувати спільне рішення задачі газової динаміки і хімічної кінетики вибухового горіння газоповітряної суміші. Для коректної роботи граничних умов вільний вихід в умови розривних течій розроблена схема апроксимації параметра в фіктивний шар на основі ударної адіабати конкретного газу. Практична значимість. Виконана модифікація методу великих часток дозволяє проводити чисельний експеримент з розрахунку безпечних відстаней при аварійних газових вибухах в умовах вугільних шахт, а також на основі розрахунку поширення ударної повітряної хвилі по каналу визначати динамічні навантаження на вибухозахисні споруди. Цель. Разработка эффективной схемы численного счета совместного решения задачи газовой динамики и химической кинетики горения газовоздушной среды на основе метода крупных частиц. Методы. Математическое моделирование, численный эксперимент, анализ и обобщение и результатов. Результаты. Для совместного решения задачи газовой динамики и химической кинетики горения газовоздушной среды предлагается ввести в численную схему метода крупных частиц концентрационную функцию, которая позволяет учитывать многокомпонентный состав газовой среды. Концентрационная функция дает возможность вводить в численную схему уравнения химической кинетики в виде уравнения Аррениуса и различать компоненты химической реакции и продукты горения. В задаче расчета детонационных взрывов возникают сильные градиенты давлений, которые, при выходе фронта ударной волны на границу «свободный выход» генерируют нефизические флуктуации параметра. Для исключения их влияния проводится анализ различных видов аппроксимации параметров в фиктивный слой расчетной схемы. Из анализа физических процессов найден эффективный вид граничных условий «свободный выход» для задачи распространения ударной волны в канале. Научная новизна. Модификация численного метода крупных частиц за счет введения концентрационной функции позволяет производить совместное решения задачи газовой динамики и химической кинетики взрывного горения газовоздушной среды. Для корректной работы граничных условий свободный выход в условия разрывных течений разработана схема аппроксимации параметра в фиктивный слой на основе ударной адиабаты конкретного газа. Практическая значимость. Проведенная модификация метода крупных частиц позволяет проводить численный эксперимент по расчету безопасных расстояний при аварийных газовых взрывах в условиях угольных шахт, а также на основе расчета распространения ударной воздушной волны по каналу определять динамические нагрузки на взрывозащитные сооружения. Purpose. Development of an effective scheme for numerical calculation of the joint solution of the problem of gas dynamics and chemical kinetics of combustion of a gas-air medium on the basis of the large-particle method. Methods. Mathematical modeling, numerical experiment, analysis and generalization and results. Findings. For joint solution of problems of gas dynamics and chemical kinetics of combustion gas environments proposed in the numerical scheme of the method of large particles concentration function, which allows to take into account the multicomponent composition of the gas medium. This function is determined at the stage of formation of the calculation area and in each cell of the calculation scheme it determines the mole fraction of each substance. The function is involved in the calculation of mass flows across the boundaries of the calculated cells, determining the mass of the overflow for each substance. The concentration function makes it possible to introduce into the numerical scheme the equations of chemical kinetics in the form of the Arrhenius equation and to distinguish the chemical reaction components and combustion products. In the problem of calculating detonation explosions, strong pressure gradients arise which, when the front of the shock wave reaches the free exit boundary, nonphysical fluctuations of the parameter are generated. To exclude their influence on the process under consideration, various types of approximation of parameters in the fictitious layer of the design scheme are analyzed. From the analysis of physical processes an effective form of the boundary conditions is found for a free yield for the problem of propagation of a shock wave in a channel. Originality. Modification of the numerical method of large particles due to the introduction of a concentration function allows the joint solution of the problem of gas dynamics and chemical kinetics of explosive combustion of a gas-air medium. For correct operation of the boundary conditions, a free exit into the conditions of discontinuous flows is developed for the scheme of approximation of the parameter in a fictitious layer on the basis of the shock adiabat of a particular gas. Practical implications. The modification of the large-particle method makes it possible to conduct a numerical experiment on the calculation of safe distances in emergency gas explosions in coal mine conditions, and also on the basis of calculating the propagation of a shock air wave through a channel to determine the dynamic loads on explosion-proof structures. 2020 Article Чисельне моделювання аварійних вибухів рудникової атмосфери / М.М. Налисько, Л.І. Барташевська // Физико-технические проблемы горного производства: Сб. научн. тр. — 2020. — Вип. 22. — С. 85-102. — Бібліогр.: 12 назв. — укр. 2664-17716 DOI: https://doi.org/10.37101/ftpgp22.01.007 http://dspace.nbuv.gov.ua/handle/123456789/166520 622.812 uk Физико-технические проблемы горного производства Інститут фізики гірничих процесів НАН України |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| language |
Ukrainian |
| topic |
Физика горных процессов на больших глубинах Физика горных процессов на больших глубинах |
| spellingShingle |
Физика горных процессов на больших глубинах Физика горных процессов на больших глубинах Налисько, М.М. Барташевська, Л.І. Чисельне моделювання аварійних вибухів рудникової атмосфери Физико-технические проблемы горного производства |
| description |
Мета. Розробка ефективної схеми чисельного розрахунку спільного розв'язку задачі газової динаміки й хімічної кінетики горіння газоповітряного середовища на основі методу великих часток.
Методика. Математичне моделювання, чисельний експеримент, аналіз і узагальнення й результатів.
Результати. Для спільного рішення задачі газової динаміки і хімічної кінетики горіння газоповітряної суміші пропонується ввести в чисельну схему методу великих часток концентраційну функцію, яка дозволяє враховувати багатокомпонентний склад газового середовища. Концентраційна функція дає можливість вводити в чисельну схему рівняння хімічної кінетики в вигляді рівняння Арреніуса і розрізняти компоненти хімічної реакції і продукти горіння. У задачі розрахунку детонаційних вибухів виникають сильні градієнти тисків, які, при виході фронту ударної хвилі на границю вільний вихід генерують нефізичні флуктуації параметра. Для виключення їх впливу проводиться аналіз різних видів апроксимації параметрів в фіктивний шар розрахункової схеми. З аналізу фізичних процесів знайдений ефективний вид граничних умов вільний вихід для задачі поширення ударної хвилі в каналі.
Наукова новизна. Модифікація чисельного методу великих часток за рахунок введення концентраційної функції дозволяє отримувати спільне рішення задачі газової динаміки і хімічної кінетики вибухового горіння газоповітряної суміші. Для коректної роботи граничних умов вільний вихід в умови розривних течій розроблена схема апроксимації параметра в фіктивний шар на основі ударної адіабати конкретного газу.
Практична значимість. Виконана модифікація методу великих часток дозволяє проводити чисельний експеримент з розрахунку безпечних відстаней при аварійних газових вибухах в умовах вугільних шахт, а також на основі розрахунку поширення ударної повітряної хвилі по каналу визначати динамічні навантаження на вибухозахисні споруди. |
| format |
Article |
| author |
Налисько, М.М. Барташевська, Л.І. |
| author_facet |
Налисько, М.М. Барташевська, Л.І. |
| author_sort |
Налисько, М.М. |
| title |
Чисельне моделювання аварійних вибухів рудникової атмосфери |
| title_short |
Чисельне моделювання аварійних вибухів рудникової атмосфери |
| title_full |
Чисельне моделювання аварійних вибухів рудникової атмосфери |
| title_fullStr |
Чисельне моделювання аварійних вибухів рудникової атмосфери |
| title_full_unstemmed |
Чисельне моделювання аварійних вибухів рудникової атмосфери |
| title_sort |
чисельне моделювання аварійних вибухів рудникової атмосфери |
| publisher |
Інститут фізики гірничих процесів НАН України |
| publishDate |
2020 |
| topic_facet |
Физика горных процессов на больших глубинах |
| url |
http://dspace.nbuv.gov.ua/handle/123456789/166520 |
| citation_txt |
Чисельне моделювання аварійних вибухів рудникової атмосфери / М.М. Налисько, Л.І. Барташевська // Физико-технические проблемы горного производства: Сб. научн. тр. — 2020. — Вип. 22. — С. 85-102. — Бібліогр.: 12 назв. — укр. |
| series |
Физико-технические проблемы горного производства |
| work_keys_str_mv |
AT nalisʹkomm čiselʹnemodelûvannâavaríjnihvibuhívrudnikovoíatmosferi AT bartaševsʹkalí čiselʹnemodelûvannâavaríjnihvibuhívrudnikovoíatmosferi |
| first_indexed |
2025-07-14T22:02:08Z |
| last_indexed |
2025-07-14T22:02:08Z |
| _version_ |
1837661475925131264 |
| fulltext |
Физико-технические проблемы горного производства 2020, вып. 22
85
УДК 622.812 https://doi.org/10.37101/ftpgp22.01.007
ЧИСЕЛЬНЕ МОДЕЛЮВАННЯ АВАРІЙНИХ ВИБУХІВ
РУДНИКОВОЇ АТМОСФЕРИ
М.М. Налисько
1*
, Л.І. Барташевська
2
1
Державний вищий навчальний заклад «Придніпровська державна академія
будівництва та архітектури», м. Дніпро, Україна
2
Національний технічний університет «Дніпровська політехніка», м. Дніпро,
Україна
*
Відповідальний автор: e-mail: 59568@i.ua
NUMERICAL SIMULATION OF AN EMERGENCY EXPLOSIONS
OF THE MINING ATMOSPHERE
М.М. Nalisko
1*
, L.I. Bartashevska
2
1
State Higher Education Establishment «Pridneprovsk State Academy of Civil
Engineering and Architecture», Dnipro, Ukraine
2
Dnipro University of Technology, Dnipro, Ukraine
*
Corresponding author: e-mail: 59568@i.ua
ABSTRACT
Purpose. Development of an effective scheme for numerical calculation of the
joint solution of the problem of gas dynamics and chemical kinetics of combus-
tion of a gas-air medium on the basis of the large-particle method.
Methods. Mathematical modeling, numerical experiment, analysis and generali-
zation and results.
Findings. For joint solution of problems of gas dynamics and chemical kinetics of
combustion gas environments proposed in the numerical scheme of the method of
large particles concentration function, which allows to take into account the mul-
ticomponent composition of the gas medium. This function is determined at the
stage of formation of the calculation area and in each cell of the calculation
scheme it determines the mole fraction of each substance. The function is in-
volved in the calculation of mass flows across the boundaries of the calculated
cells, determining the mass of the overflow for each substance. The concentration
function makes it possible to introduce into the numerical scheme the equations of
chemical kinetics in the form of the Arrhenius equation and to distinguish the
chemical reaction components and combustion products. In the problem of calcu-
lating detonation explosions, strong pressure gradients arise which, when the front
of the shock wave reaches the free exit boundary, nonphysical fluctuations of the
parameter are generated. To exclude their influence on the process under consid-
eration, various types of approximation of parameters in the fictitious layer of the
design scheme are analyzed. From the analysis of physical processes an effective
https://doi.org/10.37101/ftpgp22.01.007
mailto:59568@i.
mailto:59568@i.
Физико-технические проблемы горного производства 2020, вып. 22
86
form of the boundary conditions is found for a free yield for the problem of prop-
agation of a shock wave in a channel.
Originality. Modification of the numerical method of large particles due to the
introduction of a concentration function allows the joint solution of the problem of
gas dynamics and chemical kinetics of explosive combustion of a gas-air medium.
For correct operation of the boundary conditions, a free exit into the conditions of
discontinuous flows is developed for the scheme of approximation of the parame-
ter in a fictitious layer on the basis of the shock adiabat of a particular gas.
Practical implications. The modification of the large-particle method makes it
possible to conduct a numerical experiment on the calculation of safe distances in
emergency gas explosions in coal mine conditions, and also on the basis of calcu-
lating the propagation of a shock air wave through a channel to determine the dy-
namic loads on explosion-proof structures.
Keywords: air-gas mixture, emergency explosion, numerical calculation, large
particle method, concentration function, non-reflecting boundary
1. ВСТУП
Сучасний рівень проектування техніки і технологій неможливий без за-
стосування математичного моделювання процесів і, зокрема, використання
чисельних методів розв’язання диференціальних рівнянь. Цьому сприяє зна-
чний прогрес у збільшенні обчислювальних потужностей персональних
комп’ютерів. Особливо це актуально для швидкоплинних аварійних проце-
сів, наслідки яких необхідно з достатнім ступенем точності прогнозувати і
враховувати в заходах щодо захисту персоналу і зменшення негативних нас-
лідків надзвичайних ситуацій.
В умовах гірничих виробок вугільних шахт найбільш небезпечні, руйнівні
з недостатньо прогнозованими наслідками види аварій є поземні пожежі та
вибухи рудникової атмосфери. Причому, при підземних пожежах виникає
небезпека неодноразових вибухів атмосфери яка збагачується продуктами
піролізу вугільної речовини. Для захисту гірничорятувальників від ударних
повітряних хвиль використовується захист відстанню та зведення вибухоза-
хисних споруд. Прогнозування параметрів ударних повітряних хвиль в цих
умовах є основним завданням для достовірного визначення безпечних відс-
таней при проведенні аварійних робіт і визначення динамічних навантажень
для розрахунку стійкості вибухозахисних споруд.
2. ПОСТАНОВКА ЗАДАЧІ
Розрахунок параметрів ударних повітряних хвиль цікавив дослідників
дуже давно. Спочатку для оцінки наслідків вибухів на промислових об’єктах
та у військовій справі використовувалися емпіричні методи, які узагальню-
вали дані вибухів. Пізніше стали використовуватися експериментальні ме-
тоди. Так при моделюванні вибухових процесів, на підставі закономірностей
подібності, був сформульований принцип «кубічного кореня» Хопкінса-
Кранца, який покладено в основу прогнозування наслідків зосереджених ви-
бухів:
Физико-технические проблемы горного производства 2020, вып. 22
87
√
⁄ ,
де R – відстань від центру заряду, Е – повна енергія вибуху.
Цей принцип закладений у формулу А. Садовського, яка до теперішнього
часу використовується у нормативних документах з розрахунку вибухових
навантажень на інженерні споруди (атомні електростанції, укриття) від дії
ударних повітряних хвиль:
√
(
√
)
,
де ΔPф – надлишковий тиск у фронті УПХ, G – маса ВР у тротиловому екві-
валенті, R – відстань від осередку вибуху.
В умовах гірничих виробок так само використовувалися експерименталь-
ні методи. Чинні нормативні методики з розрахунку параметрів ударних
хвиль в гірничих виробках, засновані на даних натурних вимірювань, прове-
дених А.М. Чеховських та В.І. Гудковим в штольні шахти «Карагайлінська»
ВО «Кіселевсьвугілля» наприкінці 60-х років минулого століття.
3. МЕТОДИ ДОСЛІДЖЕННЯ
В останнє десятиліття для розрахунку параметрів ударних повітряних
хвиль, застосовується чисельне моделювання розривних газодинамічних те-
чій [1, 2]. Так в роботах [3, 4] використовується метод кінцевих різниць для
розв’язання систем рівнянь газової динаміки, який дозволяє отримувати зна-
чення багатьох параметрів поширення ударних повітряних хвиль по мережі
виробок. В роботі [5] показано застосування цього методу у програмному
комплексі FIRE для розрахунку параметрів детонації. Як показано в роботах
А.А. Самарського, складність застосування цього методу полягає в побудові
однорідних різницевих схем для розрахунку розривних течій. Останні пере-
міщаються по масі, а параметри течії по обидві сторони від розриву
пов’язані умовою Гюгоніо, що викликає різкі коливання сіткової функції за
фронтом ударної хвилі (рис. 1).
Для вирішення цього питання у чисельній схемі необхідно використову-
вати метод «розтягування» фронту за рахунок введення в систему різнице-
вих рівнянь дисипативних членів (псевдов’язкості).
Більш стійким у розрахунках розривних течій є метод великих часток за
рахунок наявності в різницевих схемах схемної в’язкості. В роботі [6] даний
метод застосовується для розрахунку аварійних вибухів в умовах топок га-
зових котлів. Однак, в розрахунку не враховується кінетика процесу вибухо-
вого горіння газових сумішей, яка обумовлює параметри ударної хвилі.
Физико-технические проблемы горного производства 2020, вып. 22
88
Рисунок 1. Розрахунок руху ударної хвилі
При розрахунку розривних течій як правило виникає проблема постанов-
ки коректних умов на штучних межах розрахункової області. У задачі роз-
рахунку детонаційних вибухів виникають сильні градієнти тисків, які, при
виході фронту ударної хвилі на кордон вільний вихід генерують нефізичні
флуктуації параметра. Для виключення їх впливу на процес необхідно вста-
новлювати межі на значних відстанях, що сильно збільшує обсяги обчис-
лень.
4. РЕЗУЛЬТАТИ ДОСЛІДЖЕННЯ
Для розрахунку процесу вибуху газоповітряних сумішей і поширення
ударних повітряних хвиль в гірських виробках пропонується використову-
вати метод великих часток (метод Давидова) [7]. Основні положення цього
методу в розглянутих умовах наступні.
Рух середовища в циліндричній системі координат описується рівняння-
ми Ейлера (у дивергентному вигляді) нерозривності, руху, енергії:
0div W
t
нерозривності,
0
0
u P
div uW
t z
v P
div vW
t r
руху, (1)
0
E
div EW div PW
t
енергії,
де ρ – щільність; Р – тиск; – вектор швидкості; u, v – компоненти швидко-
сті W по осі z та r; z, r – цилиндрические координаты;
( ) –
повна енергія; J – внутрішня енергія газу; τтр – напруження поверхневих сил
тертя газового потоку об стінку; q – щільність теплового потоку в стінку ка-
налу; S, П – поперечний переріз та периметр виробки; qх – тепловий ефект
хімічної реакції горіння вуглеводнів; α – молярна частка вуглеводнів у шах-
тній атмосфері; t – час.
Для замикання цієї системи використовується рівняння стану ідеального
газу:
Физико-технические проблемы горного производства 2020, вып. 22
89
P = (γ – 1) ρ∙J,
де γ – показник адіабати.
Задача вирішена в циліндричній системі координат, в якій розрахункова
область представлена у вигляді циліндричного каналу. По суті, така схема є
ударною трубою з ділянкою, що заповнена газоповітряної сумішшю. Швид-
ке горіння суміші (дефлаграція або детонація) викликає формування і поши-
рення ударних повітряних хвиль у циліндричному каналі (рис. 2).
Рисунок 2. Загальна структура розрахункової сітки в плоскому поданні: М1,
М4 – тип границі розрахункової області «вільний вихід», М2, М3 – тип границі
розрахункової області «непротікання»; v, u – компоненти вектора швидкості
Схема методу великих часток, питання його стійкості і порядку апрокси-
мації досить докладно розглянуті в літературі [12]. Різницева схема рішення
нестаціонарних системи диференціальних рівнянь (1) ґрунтується на ідеї ро-
зщепленні цієї системи за фізичними процесами. Завдання вирішується в три
етапи: Ейлером, Лагранжа і заключний. На ейлеровому етапі, проміжні зна-
чення швидкості ũ, ṽ і енергії потоку Ẽ визначаються з умови «заморожу-
вання» поля щільності (∂ρ/∂t)=0, тому чисельна апроксимація рівнянь руху і
енергії (1), в момент часу tn, в циліндричній системі координат r, z буде
представлена наступними явними кінцево-різницевими алгебраїчними рів-
няннями першого порядку точності за часом і другого порядку за просто-
ром:
̃
̃
̃
[
( )
( )
]
, (2)
де i, j – адреса комірки прямокутної розрахункової сітки.
Физико-технические проблемы горного производства 2020, вып. 22
90
На лагранжевому етапі обчислюємо перетоки речовини за час Δt між осе-
редками за наступними алгебраїчним рівнянням, отриманим з чисельної ап-
роксимації рівняння нерозривності (1).
Потоки мас по осі z (в залежності від напрямку):
{
( )
̃
̃
якщо ̃
̃
( )
̃
̃
якщо ̃
̃
.
Потік мас по осі r (в залежності від напрямку):
{
( )
̃
̃
якщо ̃
̃
( )
̃
̃
якщо ̃
̃
.
На заключному етапі розраховуємо нові значення всіх параметрів потоку
з урахуванням перетікання мас:
, 0,5 0,5,
, 0,5
1 2 1
, , , 1, 0,5, , , , 1,
2
, , 1 , ,
( 0,5) (1) (2) (3)
(4) {( 0,5)
i j i j
i j
n n n n n n n n n n n
i j i j i j i j i j i j i j i j i j
n n n n n
i j i j i j i j
X j r z D X M D X M D X M
D X M X j r z
0,5, , 0,5 0,5,, , ,
, , 0,5
[1 (1)] [1 (2)] [1 (3)]
[1 (4)] } .
i j i j i j
n n n n n n
i j i j i j
n n
i j i j
D M D M D M
D M
де
– функція-ознака напряму перетікання мас через кордон розрахунко-
вих осередків;
, ̃
– значення параметра на новому часовому шарі і
проміжне його значення.
В результаті цих розрахунків нам стають відомими величини ρ
n+1
, u
n+1
,
v
n+1
, E
n+1
на новому часовому шарі.
Знаючи ці величини, можна визначити внутрішню енергію
де
W
2
= u
2
+ v
2
, а отже, визначити значення тиску за формулами рівнянь стану.
Спочатку, чисельна схема методу великих часток дозволяє розраховувати
течії в однокомпонентних газових системах. В даному випадку необхідно
враховувати горіння газоповітряних сумішей, причому горючих компонен-
том може бути кілька, тобто в розрахунковій осередку можуть бути присут-
німи відразу кілька речовин, тому для розрахунку тиску в ній необхідно зна-
ти концентрації цих речовин. Для цієї мети вводиться концентрационная
функція
( ), для якої, наприклад, k = 1 відповідає вуглеводню (метан),
k = 2 – кисень, k = 3 – азот.
Перед тим як розрахувати потоки маси через поверхню розрахункової
комірки (3), (4), необхідно знати масиви
( ), тобто знати концентрації
речовин. Це задається початковими умовами. При розрахунках потоків маси
Физико-технические проблемы горного производства 2020, вып. 22
91
через кордони рахункових осередків, після введення
( ), враховується
той факт, що загальний потік маси дорівнює сумі потоків мас окремих ком-
понентів. Для цього спочатку розраховуються загальні потоки маси, потім
потоки окремих компонентів (передбачається одношвидкісна модель). Для
цього проводиться заміна . Таким чином, визначаються
( ),
α = 1…3.
Після цього визначаємо концентрації:
1
,
1
,1
,
)(
n
ji
n
jin
jiC
.
Розрахунок тиску необхідно вести за формулою:
4
1
1
,
1
,
1
, )()(
n
ji
n
ji
n
ji CPP ,
де
( ) – парціальний тиск визначається за рівняннями стану для кожної
компоненти газового середовища.
Застосування концентраційної функції дозволяє вирішити кінетичні рів-
няння хімічної реакції. Хімічна реакція представляється у вигляді однієї
брутто-стадії: «вихідні речовини продукти», а швидкість реакції предста-
вляється в формі Арреніуса:
exp ivi a
i
dc E
Z Пc
dt RT
, (3)
де Z, Еa, vi – ефективні предекспоненціальний множник, енергія активації і
порядок реакції по i-му компоненту.
У нашому випадку протікають бімолекулярні реакції, тому рівняння (3)
по компоненту «метан» має вигляд:
1 21
1 2exp v vaEdc
Z c c
dt RT
. (4)
З макрокінетичного рівняння горіння метану слід, що швидкість витра-
чання кисню в 2 рази більше швидкості витрачання метану:
1 21
2
dc dc
dt dt
. (5)
Проинтегрировав за часом рівняння (5), отримаємо , де А –
постійна інтегрування, яку можна визначити з початкових умов:
2 12H HA c c .
Індекс "Н" вказує на початкові значення концентрацій.
Физико-технические проблемы горного производства 2020, вып. 22
92
В результаті отримаємо:
( ) ( ). (6)
Таким чином, знаючи залежність зміни концентрації від часу c1(t), визна-
чаємо за формулою (6) функцію c2(t).
Нехай в початковий момент часу c1=c1H, тому c2=c2H. При повному виго-
рянні метану c1 = 0, концентрація кисню буде дорівнює: .
При стехіометричному складі . Тому метан і кисень повністю
вигорають одночасно с1 = 0, с2 = 0. Якщо, , то , тобто спо-
стерігається залишок О2. При кисень вигорає раніше, а залишок
метану дорівнює
.
Для опису реакції горіння метан-кисень слід особливо відзначити, що в
рівнянні (4) в якості концентрації «c», використовується кількість речовини
(в молях), що припадає на одиницю об’єму в см
3
, тобто розмірність
[ ]
.
У рівняння газової динаміки входить щільність речовини. Хімічна реакція
не змінює цю величину, так як вихідні речовини просто замінюються проду-
ктами реакції. Тому, гомогенні щільності речовин, що входять до складу ат-
мосфери нам знадобляться тільки для завдання початкових умов.
Чисельне рішення рівняння хімічної кінетики спільно з рівняннями газо-
вої динаміки і рівняннями стану виглядає наступним чином. Різницевий
аналог рівняння (4) представлявся у вигляді:
1 2
1
1 1, ,
1 2, ,
n n
v vn ni j i j
i j i j
c c
k c c
t
,
де ( ⁄ ) – константа швидкості хімічної реакції горіння; i, j
– цілочисельні координати розрахункової осередки (великої частки); n – но-
мер тимчасового шару.
З рівнянь (3), (5) випливає:
1 21
1 1 1 1, , , ,
v vn n n n
i j i j i j i j
с с t k с с
,
1 1
2 2 1 1, , , ,
2
n n n n
i j i j i j i j
с с с с
.
При розрахунку рівняння збереження енергії в праву частину додається
доданок:
1 2
1 2, ,
exp /
v vn n
a i j i j
t Q Z E RT c c ,
де Q – теплотворна здатність горіння метану в повітрі при стехіометричному
складі метан-кисень, Дж/м
3
.
На границях розрахункової області j = jm, i = im необхідно ставити умови,
які не перекручують параметри у розрахунковій області. Для цього вводять
Физико-технические проблемы горного производства 2020, вып. 22
93
фіктивні розрахункові осередки, в яких значення параметрів визначаються
екстраполяцією (рис. 3а). Відповідно до алгоритму МКЧ на невідбиваючих
границях можна використовувати екстраполяції нульового, першого і друго-
го порядку. Для границі М2, М3 використовуються умова непротікання: вво-
диться нульовий фіктивний осередок, у якому vi, 0 = – vi, 1, тому
.
2
,1,
.5,0
n
ji
n
jin
ji
vv
v
Вибір порядку екстраполяції для невідбиваючих границях М1, М4 вико-
наний на основі дослідження тестового завдання – вибух газової хмари у не-
обмеженому просторі (рис. 3б).
а) б)
Рисунок 3. Структура розрахункової сітки (фіктивні осередки заштриховані) –
а, і схема тестового завдання вибору порядку екстраполяції – б
Нулева екстраполяція yim+1,j = yim,j, yi,jm+1 = yi,jm дає «зрізані» значення па-
раметра (рис. 4). Лінина екстраполяція параметрів може бути отримана з на-
ступних міркувань (рис. 5):
1
1
1 2
imim
imim
im
im
imim yy
x
yy
yx
x
y
yy
,
21
12
1
1
10 2 yy
x
yy
yx
x
y
yy
.
Праві (або ліві) похідні лінійної екстраполяції дають завищені значення
параметрів в фіктивної осередку. Найбільш прийнятні значення виходять
при застосуванні центральної похідної, яка більшою мірою враховує граді-
єнт графіка в околиці кордону розрахункової області (рис. 5):
2
22
1
2
1
2
3
22
imim
imim
im
imim
im
im
imim yy
yy
y
x
yy
yx
x
y
yy .
Физико-технические проблемы горного производства 2020, вып. 22
94
Рисунок 4. Апроксимація параметра у фіктивному шарі
Рівняння екстраполяції параметрів ударної хвилі в фіктивної осередку
2-го порядку:
2
2
2
1
!2
1
!1
1
x
x
y
x
x
y
yy imim
,
де
,1
x
yy
x
y ii
x
x
y
x
y
x
y imim
1
2
2
.
Рисунок 5. Вихід ударної хвилі на межі розрахункової області при лінійній
екстраполяції
Після підстановки отримаємо:
)2(
2
1
2 2111 imimimimimim yyyyyy .
Аналіз результатів розрахунку значень параметрів, визначених за форму-
лою (4), показав, що вони є досить завищеними (рис. 5). Таким чином, най-
більш прийнятним є використання лінійної екстраполяції параметрів.
Физико-технические проблемы горного производства 2020, вып. 22
95
Рисунок 6. Апроксимація параметра в негативну область
Однак такий вид екстраполяції можливий при малих градієнтах функції.
При розривних характеристиках, тобто у фронті ударної хвилі, лінійна апро-
ксимація призводить до виникнення протилежних за знаком значень параме-
тра, що суперечить фізичному змісту задачі (рис. 6). У розрахунку це прояв-
ляється виникненням генерації коливань тисків після проходження піку уда-
рної хвилі через границю (рис. 7). Подібний ефект спостерігається і при ну-
льовій апроксимації: у момент виходу ударної хвилі на границю розрахун-
кової області, у фіктивному шарі, виникає тиск, рівні амплітуді в ударній
хвилі. Утворюється ефект гальмування потоку на границі «вільний вихід».
Рисунок 7. Генерація нефізичних флуктуацій при апроксимації параметра
у від’ємну область: 1 – газоповітряна суміш; 2 – повітря; 3 – границя розділу
середовищ; 4 – границя розрахункової області; 5 – фактичне положення фрон-
ту ударної хвилі; 6 – відбивання від границі розрахункової області
Завдання знаходження коректного рішення виконується проведенням спі-
льного дослідження граничних умов і фізичних процесів. Так в роботі [8]
екстраполяція нульового порядку для швидкості потоку визначається спів-
відношенням, отриманим з умов Гюгоніо на стрибку ущільнень:
))((1 pp
v ,
де a∞, ρ∞, p∞ – швидкість звуку, щільність і тиск у незбуреному потоці газу.
Физико-технические проблемы горного производства 2020, вып. 22
96
У розглянутій задачі ударна хвиля поширюється по суміші газів і значен-
ня її амплітуди в наступний момент часу залежить від термодинамічних вла-
стивостей цієї суміші. Тиск в газі має теплове походження, воно пов’язане з
перенесенням імпульсу частинками, які беруть участь в тепловому русі і
завжди визначається видом ударної адіабати, тобто пропорційно обсягу і те-
мпературі, а в ударної повітряної хвилі і швидкості руху потоку [11]. Тому
для визначення надлишкового тиску в фіктивної осередку можна отримати
рівняння ударної адіабати в системі Р–U. Згідно робіт Я.Б. Зельдович,
Ю.П. Райзер (1966) залежність між швидкістю фронту УВВ і швидкістю ре-
човини за фронтом хвилі, в широкому діапазоні амплітуд, є лінійної:
D = A + Bu, (7)
де A, В – коефіцієнти.
Так само відома швидкість руху фронту УПХ:
11
a
PD , (8)
де ρа, ρ – щільність газу, атмосферна і в ударній хвилі.
Вирішуючи спільно рівняння (5), (6) отримаємо умову «м’якої» стінки:
imima
imimim
im
VV
uBA
P
,
2
1
1
)(
,
де ujm+1 – швидкість речовини в фіктивної осередку; Va, V – питома щіль-
ність газу.
У формулі (5) значення А близько до швидкості звуку в речовині і зада-
ється таблично. Значення коефіцієнта В визначаємо по граничній комірці,
вирішуючи рівняння (7), (8) щодо В:
im
imimimaim
im
u
AVVP
B
)( ,
, (9)
Використання залежності (9) також має свою особливість. Величину
швидкості речовини в ній для граничної осередки логічно визначати як пов-
ну швидкість: 22
jmjmjm vuU . Однак це не дає правильних результатів обчис-
лення тисків в фіктивних осередках з двох причин. Перше: при розрахунку
jmU на нижній межі відбувається втрата знака в напрямку швидкості, так як
швидкість має від’ємне значення при русі речовини вниз (рис. 8).
Физико-технические проблемы горного производства 2020, вып. 22
97
Рисунок 8. Напрямок руху речовини в фіктивних комірках і розрахунковий на-
прямок потоку
Друге: розрахунок перетікання речовини згідно МКЧ визначається через
кордони осередки, тобто строго по вертикалі або по горизонталі. Реальна
швидкість речовини в осередках на периферії від осей симетрії буде спрямо-
вана під кутом до них (максимально 45°). Тому тиск в крайніх фіктивних
осередках буде визначено невірно. Для виключення таких явищ у виражен-
ні (7) у якості швидкості повинна використовуватися відповідна компонента
швидкості потоку. У розрахунку це реалізується встановленням окремих
умов для лівої, правої, верхньої та нижньої меж розрахункової області.
Рисунок 9. Перехід фронту УПХ через штучну межу «вільний вихід» за умови
«м’яка» стінка: 1 – газоповітряна суміш; 2 – повітря; 3 – лінія розділу середо-
вищ; 4 – межа розрахункової області; 5 – фронт УПХ після переходу межи «ві-
льний вихід»; 6 – фронт УПХ перед границею «вільний вихід»
МКЧ передбачає виконання трьох етапів розрахунків. Відповідно і грани-
чні умови виконуються в кінці кожного етапу. Для стійкості рахунку необ-
хідною і достатньою умовою є застосування умови «м’якої» границі в за-
ключній частині лагранжевого етапу для розрахунку тисків, а граничні умо-
ви в ейлеровому етапі повинні визначаться за нульовою апроксимації через
проміжного характеру результатів. В результаті, розрахункова схема дозво-
ляє коректно виконати перехід фронту ударної хвилі через кордон розрахун-
кової області (рис. 9).
Физико-технические проблемы горного производства 2020, вып. 22
98
5. ОБГОВОРЕННЯ РЕЗУЛЬТАТІВ
Проблеми коректної роботи невідбиваючих граничних умов (НГУ) існу-
ють в практично у всіх численних схемах при розрахунку розривних те-
чій [8]. Для забезпечення НГУ на відкритій границі, користуються різними
прийомами: приєднують буферні обсяги з розрахунковими комірками, за-
стосовують евристичні залежності тиску від часу релаксаційоного процесу
або використовують одномірні характеристичні НГУ. Такі прийоми орієнто-
вані на вирішення конкретнї задачі [9, 10]. Загальний підхід до постановки
НГУ, на сьогоднішній день, відсутній. Приклад похибки НГУ можна проде-
монструвати в системі FlowVision. У моделі повністю стискаємого середо-
вища, граничні умови по типу інваріанта Рімана, призначених для розрахун-
ку транс- і надзвукових течій, дають падіння тиску в незбуреному потоці
(рис. 10).
Рисунок 10. Похибка роботи граничних умов вільний вихід за типом інваріанта
Рімана
Отримане рішення крайових умов для розривних функцій (ударної хвилі)
підвищує стійкість явної чисельної схеми методу великих часток, при цьому
похибка
6. ВИСНОВКИ
Модифікація методу великих часток дозволяє розширити область його за-
стосування на спільне вирішення задачі газової динаміки і хімічної кінетики
горіння газоповітряних сумішей. Умова «м’якої» стінки дозволяє не перек-
ручувати параметри розривних течій при переході границі розрахункової
області і, таким чином, дає можливість зменшити обсяги обчислень за раху-
нок зменшення розмірів розрахункової області.
СПИСОК ЛІТЕРАТУРИ
1. Бабкин, А.В., Колпаков, В.И., Охитин, В.Н. & Селиванов, В.В. (2006). При-
кладная механика сплошных сред. Т.3 Численные методы в задачах физики быст-
ропротекающих процессов. Москва: Изд-во МГТУ им. Н.Э. Баумана, 520 с.
2. Кукуджанов, В.Н. Численные методы в механике сплошных сред. (2006).
Москва: «МАТИ»-РГТУ, 158 с.
Физико-технические проблемы горного производства 2020, вып. 22
99
3. Васенин, И.М., Шрагер, Э.Р., Крайнов, А.Ю. & Палеев Д.Ю. (2011). Матема-
тическое моделирование нестационарных процессов вентиляции сети выработок
угольной шахты. Компьютерные исследования и моделирование, (2), 155–163.
4. Агеев, В.Г., Греков, С.П., Зинченко, И.Н. & Салахутдинов Т.Г. (2013). Ком-
пьютерное моделирование развития, распространения и локализации взрывов ме-
тановоздушных смесей в горных выработках. Вісник Харківського національного
університету, (1058), 5–12.
5. Скоб, Ю.А. & Угрюмов, М.Л. (2013). Численное моделирование детонации в
газовых смесях. Вісник Харківського національного університету, (1058), 149–157.
6. Поландов, Ю.Х. & Бабанков, В.А. (2014). Влияние места расположения ис-
точника воспламенения в помещении на развитие взрыва газа. Пожаровзрывобез-
опасность, (3), 68–74.
7. Егоров, М. Ю. Метод Давыдова – современный метод постановки вычисли-
тельного эксперимента в ракетном твердотопливном двигателестроении. (2014).
Вестник Пермского национального исследовательского политехнического универ-
ситета. Серия: Аэрокосмическая техника, (37), 6–70.
8. Ильгамов, М.А. & Гильманов, А.Н. (2003). Неотражающие условия на грани-
цах расчетной области. Москва: Физматлит, 240 с.
9. Лидский, Б.В., Посвянский В.С. & Фролов С.М. (2009). Неотражающие гра-
ничные условия на открытых границах для сжимаемых и несжимаемых многомер-
ных течений. Горение и взрыв, (2), Москва: Торус Пресс, 31–35.
10. Поздеев, С.В., Некора, О.В., Демешок, В.В. & Медведь, Б.Ю. (2016). Иссле-
дование поведения деревянной плиты при пожаре с помощью метода конечных
элементов. Строительство, материаловедение, машиностроение: сб. научн. тру-
дов, (93), 25–31.
11. Чернай, А.В., Соболев, В.В., Илюшин, М.А. & Житник, Н.Е. (1994). О методе
получения механических импульсов нагружения, основанном на лазерном подрыве
покрытий из взрывчатых составов. Физика горения и взрыва, (2), 106–111.
12. Белоцерковский, О.М. & Давыдов, Ю.М. (1982). Метод крупных частиц в
газовой динамике. Москва: Наука, 392 с.
REFERENCES
1. Babkin, A.V., Kolpakov, V.I. Ohitin, V.N. and Selivanov, V.V. (2006), Prikladna-
ya mehanika sploshnyih sred. T.3 Chislennyie metodyi v zadachah fiziki byistroprotekay-
uschih protsessov [Applied continuum mechanics. Vol. 3 Numerical methods in problems
of physics of fast processes], MGTU im. N.E. Baumana, Moscow, Russia, 520 s
2. Kukudzhanov, V.N. (2006), Chislennyie metodyi v mehanike sploshnyih sred [Nu-
merical methods in continuum mechanics], «MATI»-RGTU, Moscow, Russia, 158 s
3. Vasenin, I.M., Shrager, E.R., Kraynov, A.Y. and Paleev, D.Y. (2011), “The mathe-
matical modelling of nonsteady ventilation processes of coal mine working net”, Com-
puter researches and modelling, (2), 155–163.
4. Ageev, V.G., Grekov, S.P., Zinchenko, I.N. and Salahutdinov, T.G. (2013), “Com-
puter simulation development, spread and localization of explosions of methane-air mix-
tures in mines”, The Journal of V.N.Karazin Kharkiv National University, (1058), 5–12.
5. Skob, Yu.A. (2013), “Numerical modeling of detonation in gas mixtures”, The
Journal of V.N.Karazin Kharkiv National University, (1058), 149–157.
6. Polandov, Yu.H. and Babankov V.A. (2014), “The effect of the location of the
source of ignition in the premises for the development of gas explosion”, Fire and Explo-
sion Safety, (3), 68–74.
Физико-технические проблемы горного производства 2020, вып. 22
100
7. Egorov, M. Yu. (2014), “Davydov’s method is a modern method of placing the
computational experiment in solid propellant engine”, PNRPU Aerospace Engineering
Bulletin, (37), 6–70.
8. Ilgamov, M.A. and Gilmanov, A.N. (2003), Neotrazhayuschie usloviya na
granitsah raschetnoy oblasti [Non-reflecting conditions on boundaries of computational
domain], Fizmatlit, Moscow, Russia, 240 s
9. Lidskiy, B.V., Posvyanskiy, V.S. and Frolov S.M. (2009). Nonreflecting boundary
conditions on open boundaries for compressible and incompressible multidimensional
flows, Combustion and explosion, Editor-in-Chief and Chair of Editorial Council Profes-
sor S.M. Frolov, (2). 31–35.
10. Pozdeev, S.V., Nekora, O.V., Demeshok, V.V. and Medved, B.Yu. (2016). Inves-
tigation of the behavior in fire timber frame with finite element method, Construction,
Materials science, Mechanical Engineering: scientific works collection, Series: Life ac-
tivity safety, SHEE “Prindeprovs’ka State Academy of Civil Engineering and Architec-
ture”, (93), 25–31.
11. Chernay, A.V., Sobolev, V.V., Ilyushin, M.A. and Zhitnik, N.E. (1994), The
method of producing the mechanical pulse loading based on the laser-blasting explosive
compositions of coatings, Combustion, Explosion, and Shock Waves, (2), 106–111.
12. Belotserkovskij O.M. and Davydov J M. (1982) Metod krupnyih chastits v
gazovoy dinamike [The large particles method in the gas dynamics]. Moscow, Nauka
Publ., 391 s
ABSTRACT (IN UKRAINIAN)
Мета. Розробка ефективної схеми чисельного розрахунку спільного розв'яз-
ку задачі газової динаміки й хімічної кінетики горіння газоповітряного сере-
довища на основі методу великих часток.
Методика. Математичне моделювання, чисельний експеримент, аналіз і
узагальнення й результатів.
Результати. Для спільного рішення задачі газової динаміки і хімічної кіне-
тики горіння газоповітряної суміші пропонується ввести в чисельну схему
методу великих часток концентраційну функцію, яка дозволяє враховувати
багатокомпонентний склад газового середовища. Концентраційна функція
дає можливість вводити в чисельну схему рівняння хімічної кінетики в ви-
гляді рівняння Арреніуса і розрізняти компоненти хімічної реакції і продук-
ти горіння. У задачі розрахунку детонаційних вибухів виникають сильні
градієнти тисків, які, при виході фронту ударної хвилі на границю вільний
вихід генерують нефізичні флуктуації параметра. Для виключення їх впливу
проводиться аналіз різних видів апроксимації параметрів в фіктивний шар
розрахункової схеми. З аналізу фізичних процесів знайдений ефективний
вид граничних умов вільний вихід для задачі поширення ударної хвилі в ка-
налі.
Наукова новизна. Модифікація чисельного методу великих часток за раху-
нок введення концентраційної функції дозволяє отримувати спільне рішення
задачі газової динаміки і хімічної кінетики вибухового горіння газоповітря-
ної суміші. Для коректної роботи граничних умов вільний вихід в умови ро-
зривних течій розроблена схема апроксимації параметра в фіктивний шар на
основі ударної адіабати конкретного газу.
Физико-технические проблемы горного производства 2020, вып. 22
101
Практична значимість. Виконана модифікація методу великих часток до-
зволяє проводити чисельний експеримент з розрахунку безпечних відстаней
при аварійних газових вибухах в умовах вугільних шахт, а також на основі
розрахунку поширення ударної повітряної хвилі по каналу визначати дина-
мічні навантаження на вибухозахисні споруди.
Ключові слова: газоповітряна суміш; аварійний вибух; чисельний розраху-
нок; метод великих часток; концентраційна функція; невідбиваюча границя
ABSTRACT (IN RUSSIAN)
Цель. Разработка эффективной схемы численного счета совместного реше-
ния задачи газовой динамики и химической кинетики горения газовоздуш-
ной среды на основе метода крупных частиц.
Методы. Математическое моделирование, численный эксперимент, анализ и
обобщение и результатов.
Результаты. Для совместного решения задачи газовой динамики и химиче-
ской кинетики горения газовоздушной среды предлагается ввести в числен-
ную схему метода крупных частиц концентрационную функцию, которая
позволяет учитывать многокомпонентный состав газовой среды. Концен-
трационная функция дает возможность вводить в численную схему уравне-
ния химической кинетики в виде уравнения Аррениуса и различать компо-
ненты химической реакции и продукты горения. В задаче расчета детонаци-
онных взрывов возникают сильные градиенты давлений, которые, при вы-
ходе фронта ударной волны на границу «свободный выход» генерируют не-
физические флуктуации параметра. Для исключения их влияния проводится
анализ различных видов аппроксимации параметров в фиктивный слой рас-
четной схемы. Из анализа физических процессов найден эффективный вид
граничных условий «свободный выход» для задачи распространения удар-
ной волны в канале.
Научная новизна. Модификация численного метода крупных частиц за
счет введения концентрационной функции позволяет производить совмест-
ное решения задачи газовой динамики и химической кинетики взрывного
горения газовоздушной среды. Для корректной работы граничных условий
свободный выход в условия разрывных течений разработана схема аппрок-
симации параметра в фиктивный слой на основе ударной адиабаты конкрет-
ного газа.
Практическая значимость. Проведенная модификация метода крупных
частиц позволяет проводить численный эксперимент по расчету безопасных
расстояний при аварийных газовых взрывах в условиях угольных шахт, а
также на основе расчета распространения ударной воздушной волны по ка-
налу определять динамические нагрузки на взрывозащитные сооружения.
Ключевые слова: газовоздушная смесь; аварийный взрыв; численный рас-
чет; метод крупных частиц; концентрационная функция; неотражающая
граница
Физико-технические проблемы горного производства 2020, вып. 22
102
ABOUT AUTHORS
Nalisko Nikolay, Candidate of Technical Science, Associate Professor, Department of
Life Safety, State Higher Education Establishment “Pridneprovsk State Academy of Civil
Engineering and Architecture”, 24-A Chernishevskogo Str., Dnipro, Ukraine, 49600.
E-mail: 59568@i.ua.
Bartashevskaya Lyudmila, Candidate of physical and mathematical sciences, Associ-
ate Professor, Department of Physics, Dnipro University of Technology, 19 Dmytra Ya-
vornytskoho Ave., Dnipro, Ukraine, 49005. E-mail: bartashevska.l.i@gmail.com
mailto:59568@i.ua
https://context.reverso.net/%D0%BF%D0%B5%D1%80%D0%B5%D0%B2%D0%BE%D0%B4/%D0%B0%D0%BD%D0%B3%D0%BB%D0%B8%D0%B9%D1%81%D0%BA%D0%B8%D0%B9-%D1%80%D1%83%D1%81%D1%81%D0%BA%D0%B8%D0%B9/physical
http://physics.nmu.org.ua/en/
|