CC BY
27
ГИАБ. Горный информационно-аналитический бюллетень / MIAB. Mining Informational and Analytical Bulletin, 2020;(11-1):118-134 ОРИГИНАЛЬНАЯ СТАТЬЯ / ORIGINAL PAPER
УДК 911.2, 504.05, 551.582 DOI: 10.25018/0236-1493-2020-111-0-118-134
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ РАСПРЕДЕЛЕНИЯ ПЫЛИ ПО СКЛОНАМ ГОРНОГО УЩЕЛЬЯ ОТ ХВОСТОХРАНИЛИЩА, РАСПОЛОЖЕННОГО В АЛАГИРСКОМ УЩЕЛЬЕ
Е.С. Каменецкий12, А.А. Радионов2, В.Ю. Тимченко2, О.С. Панаэтова2, Г.И. Свердлик1
1 Северо-Кавказский горно-металлургический институт (государственный технологический университет), Владикавказ, Россия 2 Южный математический институт Владикавказского научного центра РАН, Владикавказ, Россия, esk@smath.ru
Аннотация: Добыча полезных ископаемых является источником переносимых по воздуху частиц загрязняющих веществ (ЗВ) в результате прямых выбросов и ветровой эрозии хвостохранилищ. Летние и сухие погодные условия способствуют насыщению атмосферы пылью и аэрозолями, содержащими тяжелые металлы, что является определяющим фактором для миграции элементов в природной среде. Рассеивание пыли и аэрозолей от Унальского хвостохранилища Садонского свинцово-цинкового комбината (ССЦК), расположенного вблизи селения Унал, Алагирский район, РСО-Алания, РФ, происходит в горном ущелье. Эти хвосты сильно загрязнены свинцом, цинком, мышьяком и другими элементами. С использованием вычислительной гидродинамической модели атмосферы моделируется перенос пыли от модельного источника в районе хвостохранилища и исследуется распределение пыли на склонах Алагирского ущелья. В модель включены региональные особенности ландшафта, характерные для местности розы ветров, а также уравнение переноса пыли. Расчеты концентрации ЗВ при различных направлениях внешнего ветра усреднялись с коэффициентами, равными вероятности соответствующего ветра в розе ветров. Установлено, что особенности ландшафта местности и характерное направление скорости ветра являются важными факторами, которые определяют пространственное распределение осажденной пыли. Сопоставление расчётов с результатами полевых измерений показывает удовлетворительное согласие.
Ключевые слова: горное ущелье, пыль, хвостохранилище, математическая модель, OpenFOAM, перенос загрязняющих веществ, ландшафтно-геохимическая миграция элементов, спутниковые метеорологические измерения, атмосфера в горных ущельях.
Для цитирования: Каменецкий Е.С., Радионов А.А., Тимченко В.Ю., Панаэтова О.С., Свердлик Г.И. Математическое моделирование распределения пыли по склонам горного ущелья от хвостохранилища, расположенного в алагирском ущелье // Горный информационно-аналитический бюллетень. - 2020. - № 11-1. - С. 118-134. DOI: 10.25018/0236-1493-2020111-0-118-134.
© Е.С. Каменецкий, А.А. Радионов, В.Ю. Тимченко, О.С. Панаэтова, Г.И. Свердлик. 2020.
Mathematical modeling of dust propagation on mountain ravine walls from tailings storage facilities: A case-study of the Alagir Ravine
E.S. Kamenetsky1,2, A.A. Radionoff2, V.U. Timchenko2, O.S. Panaetova2, G.I. Sverdlik1
1 North Caucasian Institute of Mining and Metallurgy (State Technological University), Vladikavkaz, RNO-Alania, Russia; 2 Southern Mathematical Institute — the Affiliate of Vladikavkaz Scientific Center of the Russian Academy of Sciences, Vladikavkaz, RNO-Alania, Russia, esk@smath.ru
Abstract: Mining is a source of airborne metal and metalloid contaminants, both due to direct emissions and due to wind erosion of mine tailings. Summer dry weather contributes to air pollution with contaminated atmospheric dust is a major condition for migration of elements in the mountain environment. Dispersion of dust occurs in a mountain valley from the Unal's tailings, located near the village Unal, Alagir district, North Ossetia-Alania, Russia. These tailings are heavily contaminated with lead, zinc, arsenic and other elements. Using a computational hydrodynamic model of the atmosphere, dust transfer from a tailings source is simulated and dust distribution on the slopes of the Alagir's valley is investigated. The model incorporates regional topographic features, the local weather patterns (wind rose) and dust transfer equation. Calculations of pollution concentration in different directions of the external wind were averaged with coefficients equal to the probability of the corresponding direction in the wind rose. It is shown as the result of calculations that the regional topographical features of terrain and the wind rose are important factors that determine the spatial deposition of dust on the slopes of valley. Comparison of calculations with the results of field measurements is in satisfactory agreement.
Key words: mountain valley, dust, tailings, mathematical model, OpenFOAM, transport of contaminants, landscape-geochemical migration of elements, satellite meteorological measurements, atmosphere in mountain valleys.
For citation: Kamenetsky E.S., Radionoff A.A., Timchenko V.U., Panaetova O.S., Sverdlik G.I. Mathematical modeling of dust propagation on mountain ravine walls from tailings storage facilities: A case-study of the Alagir Ravine. MIAB. Mining Inf. Anal. Bull. 2020;(11-1):118-134. [In Russ]. DOI: 10.25018/0236-1493-2020-111-0-118-134.
Введение
Ландшафтно-геохимические условия миграции элементов в природной среде помимо других факторов обусловливаются и хозяйственно-промышленной деятельностью человека, активно ведущейся в горных территориях. Пыль от эродирующих отходов деятельности горнодобывающих комбинатов, называемых хвостами, попадает в атмосферу и приводит к перераспределению химических элементов, регистрируемых в почвах склонов горных ущелий. Потенциально вредное воздействие этих химических
элементов на экологию района нуждается в изучении. В летних условиях и при слабом внешнем ветре (ветер, дующий выше пиков горных вершин) пылевой аэрозоль, поднимаясь от хвостохранилищ, распространяется в основном вдоль горного ущелья вблизи его склонов в сторону гор или равнин. При заметном внешнем ветре ситуация меняется: во внутренних участках атмосферы горного ущелья возникают области значительной турбулентности, при попадании в которые пылевой аэрозоль может интенсивно подниматься вверх и выветриваться.
На процесс распространения ЗВ в горных ущельях влияет множество факторов, детальный учет которых возможен только в рамках полных математических моделей. Однако для некоторых практических приложений достаточно знать оценки влияния основных факторов на процессы распространения ЗВ в сложных горных ландшафтах.
Унальское хвостохранилище Садон-ского свинцово-цинкового комбината (ССЦК) расположено вблизи селения Унал, Алагирский район, РСО-Алания, РФ, в горном ущелье Северной части Кавказского хребта, на 42,862° северной широты и 44,145° восточной долготы, на высоте около 1700 метров над уровнем моря. Ширина ущелья вблизи хвостохранилища варьируется от 1000 м до 3000 м на высотах до 1000 м над дном ущелья. Высота склонов горного ущелья в месте расположения хвостохранилища достигает 2571 м. Чаша Унальского хвостохранилища ССЦК расположена на дне ущелья на высоте около 1677 м над уровнем моря в излучине реки Ардон и занимает площадь более 150 000 кв. м. Хвостохранилище эксплуатируется многие десятки лет. В нем хранится порядка 2,6 млн тонн промышленных отходов, в которых 0,21 % масс. содержания свинца, 0,32 % масс. содержания цинка, множество других элементов.
Для предотвращения выбросов пыли Унальское хвостохранилище покрыто слоем воды. Это защитное озеро в северной части периодически мелеет и просыхает, тогда сухие хвосты подвергаются эрозии и становятся источником пыли. Такая ситуация может возникать не каждый год, в летнее время. Количество пыли, вынесенной в атмосферу за время существования хвостохранилища, не поддается оценке. Дополнительным источником пыли, вероятно, является расположен-
ная в ущелье Транскавказская магистраль: осевшие на поверхность дороги частицы пыли могут дробиться и вновь подниматься в атмосферу большегрузным автотранспортом, движущимся по трассе в непосредственной близости от хвостохранилища. Таким образом, мощность источника пыли является неизвестной величиной при рассмотрении задачи переноса пыли.
Хвостохранилище ССЦК расположено в 700 метрах от ближайшего к нему селения Унал, немногим дальше расположено селение Зинцар. Приблизительно в двадцати километрах вниз по ущелью при выходе на равнину расположен густонаселенный район: город Алагир и прилегающие селения; выше по ущелью расположены базы отдыха, курорты. Распределение химических элементов эродирующего хвостохрани-лища по склонам Алагирского ущелья и при выходе ущелья на равнинные территории нуждается в изучении.
Полевым наблюдениям по осаждению пыли от эродирующих хво-стохранилищ и обработке результатов измерений посвящено множество работ. Экспериментально проведенные исследования, где также определяются условия отрыва пылинок от поверхности, представлены в [1 — 9]. Основными факторами, обусловливающими отрыв пылинки от поверхности, является сила ветра, который возникает над сухими эродирующими хвостами, и турбулентность воздуха, которая также связана с температурой и солнечным нагревом. Закономерности распределения осажденной пыли по поверхности со сложным ландшафтом обусловливаются топографией и направлениями преобладающих ветров. Работы [6, 7] показывают успешное применение математического моделирования аэродинамики и распространения пыли в реальных условиях.
В работах [10 — 13] изучаются особенности миграции химических элементов в горных ландшафтах Уналь-ского хвостохранилища. В работе [13] представлены результаты анализа образцов почв и грунтов, взятых в различных точках Алагирского ущелья на предмет содержания химических элементов, и также показаны фоновые значения их содержания. Найдено, что содержание некоторых химических элементов в почвах склонов ущелья может существенно превышать как фоновые значения, так и значения, регламентируемые предельно допустимыми концентрациями (ПДК). В [14] рассматриваются некоторые источники ЗВ, возникающие при горнодобывающей деятельности, и методы их учета.
Использование статистических методов математического моделирования для задачи миграции химических элементов в горных территориях представлено в работах [15, 16].
В обзорах [17, 18] описаны современные математические модели, используемые численные методы и закономерности аэродинамики горных ущелий. Для описания распространения пылевого аэрозоля могут использоваться стандартизованные хорошо отлаженные алгоритмы решения гидродинамических уравнений, содержащихся в вычислительных пакетах, таких как ANSYS для метода конечных элементов, OpenFOAM, реализующий метод конечных объемов, и других пакетов. Также используются различные модификации метода частиц, фокусирующиеся на механике движения отдельной частицы пыли в турбулентном потоке воздуха. Распределение пыли по поверхности описывается ансамблем реализаций множества траекторий таких пылинок от источника до осаждения на поверхность. Из-за сложности и многофакторности течений в горных
ущельях [17—20] внимание исследователей уделяется не только конкретным задачам, но и горным ущельям идеализированной формы при упрощенных синоптических условиях. Использование упрощенной формулировки задачи математического моделирования (например, без учета микрофизических процессов взаимодействия пыли с образующимися капельками влаги, что важно при конденсации влаги прохождении циклонов над хвостохрани-лищем) позволяет выделить влияние на перенос пыли основных факторов, таких как направление и сила ветра над хвостохранилищем и влияние особенностей ландшафта.
Для оценки вредного воздействия хвостохранилища на здоровье человека и экологию, а также ландшафтно-гео-химических условий миграции элементов, содержащихся в пыли Унальского хвостохранилища СЦКК, необходимо знать распределение ЗВ на склонах алагирского ущелья по всей его протяженности и предгорных равнинных территорий. Создание модели, точно учитывающей большое разнообразие атмосферных условий в горной атмосфере (таких как изменчивость солнечного излучения, смену сезона, влияние редко дующих сильных фёнов и осадков различных видов), за десятки лет существования Унальского хво-стохранилища — задача весьма трудоемкая и ресурсоемкая. По этой причине в настоящей работе предпринята попытка оценить распределение пыли по склонам Алагирского ущелья в упрощенной постановке с учетом только основных факторов.
Методы исследования
В работе используется ряд упрощений для решения сформулированной задачи. Климатическими «средними» условиями для каждого расчета выби-
рались условия сплошной облачности, в которых не принимаются во внимание сезоны, не учитываются осадки и различие температурных режимов склонов ущелья. Такой подход оправдывает себя при расчетах широкого класса задач, связанных с атмосферными ЗВ [9]. В модель включались региональные топографические особенности и характерные для местности розы ветров, которые получены из общедоступного источника (https:// power.larc.nasa.gov/data-access-viewer/). При моделировании в месте нахождения хвостохранилища располагался модельный источник ЗВ в виде пыли с безразмерной нормированной мощностью, равной единице. Этот источник действовал в течение всего времени расчета. Проводилось не менее 16-ти расчетов с различными направлениями ветра, которые потом усреднялись с учетом розы ветров над хвостохранилищем. Полученное поле концентрации ЗВ сопоставлялось с результатами полевых измерений, проведенных в нескольких точках ущелья [13], и нормировалось на эти значения.
Например, результаты расчета показывают, что концентрация ЗВ в районе
селения Унал от единичного источника, расположенного в районе хвостохранилища, составляет Cu = 1,8 • 10-6. Для свинца измеренное вблизи селения Унал в полевых условиях значение есть C^ = 438 мг/кг [13]. Используя эту информацию, можно получить нормировку для рассчитанного поля ЗВ по всей протяженности ущелья. После нормировки можно оценивать концентрацию ЗВ по всей протяженности ущелья в любой его точке.
Исходные данные
В табл. 1 приводятся данные измерений содержания в грунтах свинца и цинка, полученные в полевых исследованиях в селениях Архон, Унал, Зинцар и опубликованные в работе [13]. В качестве маркеров при моделировании выбраны свинец и цинк, поскольку их содержание в выбросах пыли хвостохранилища наибольшее и на поверхности нет заметных фоновых значений концентрации этих элементов. В табл. 1 ПДК указываются согласно Федеральному закону Российской Федерации «О санитарно-эпидемиологическом благополучии населения» № 52-ФЗ от 30 марта 1999 г.
Таблица 1
Результаты полевых измерений [13], мг/кг
Table 1
The results of field measurements [13], mg/kg
Элемент Сельхозугодия Пастбища Равнины
Архон (по три пробы)
Свинец (ПДК 32) 424, 423, 441 511, 5, 55 326, 315, 3117
Цинк (ПДК 23) 424, 423, 441 425, 447, 725 653, 642, 2979
Зинцар (по три пробы)
Свинец 564, 609, 576 317, 318, 321 818, 38
Цинк 695, 72, 699 368, 362, 381 471, 324
Унал (по три пробы)
Свинец 341, 341, 344 328, 534, 453 533, 541, 532
Цинк 485, 468, 473 568, 476, 473 623, 64, 606
С помощью данных НАСА, полученных с использованием спутниковых измерений близповерхностных характеристик земной атмосферы (таких как температура, скорость ветра, облачность, давление) с 1.1.1995 года по сентябрь 2019 года (https://power.larc.nasa. gov/data-access-viewer/), построена роза ветров над унальским хвостохранили-щем ССЦК (ри^ 1). Эта роза ветров построена для высоты 50 м для летних месяцев года с 16 апреля по 15 октября каждого года измерений. Зимние месяцы мало сказываются на распространении пыли.
Из рис 1 видно, что над хвостохра-нилищем преобладает восточно-юго-восточный ветер, а более сильные ветра дуют с юго-востока, что может приводить к переносу ЗВ в сторону равнины и расположенного там города
Алагира. Роза ветров является осреднением измерений по довольно большой площади (порядка 18х18 км) и характеризует дующие в этом районе ветра.
Описание модели. Для моделирования использовался вычислительный пакет OpenFOAM (https://openfoam. org/release/7/), свободно распространяемый по лицензии GNU General Public Licence (https://openfoam.org/licence/). OpenFOAM — Open Field Operation and Manipulation — это набор средств для задач механики сплошной среды, включая вычислительную гидродинамику (CFD), реализованный на объектно-ориентированном языке программирования высокого уровня С++. Этот пакет находит широкое применение как в инженерных вычислениях, так и для геофизических приложений. Пакет распространяется в исходных
Рuc. 1. Данные метеорологического спутника над районом Унальского хвостохранилища. а — роза ветров на высоте 50 м для летних месяцев периода измерений с 01.01.1995 года по сентябрь 2019 года. Кривая показывает количество наблюдений ветра данного направления за период измерений, вверх — северный ветер, вправо — восточный; б — величина скорости ветра на высоте 50 м (м/с), сплошной кривой отмечена средняя скорость ветра, пунктирной — максимальная скорость ветра
Fig 1. Meteorological satellite data over the region of the Unal tailings: a — the wind rose at altitude of 50 m for the summer months from 1.1.1995 to September 2019. The curve indicates the number of wind observations in a given direction, upward — the north wind, the right — east wind; b — magnitude of wind velocity at height of 50 m (m/s), the solid curve is depicted the average wind velocity, the dashed curve is depicted the maximum wind velocity
программных кодах, что позволяет легко модифицировать стандартные решатели (солверы) для учета действующих факторов, характерных для рассматриваемой задачи. Этот пакет использовался для описания распространения пыли [6] от промышленных предприятий; для описания распространения пыли в атмосфере над сложной подстилающей поверхностью [7], для решения других подобных задач (https://openfoam.Org/release/7).
Для нашей модели был выбран солвер pisoFoam [21], предназначенный для описания нестационарных течений несжимаемой жидкости и входящий в стандартную поставку OpenFOAM. Для решения поставленной задачи было сделано две модификации солвера: 1) введено слагаемое f х V, учитывающее силу Кориолиса, где f = (0, 20сов(ф), 20в1п(ф)) - вектор параметра Кориолиса, частота вращения Земли О = 7,3 • 10-5 1/с, Ф - географическая широта; 2) введено уравнение для концентрации пыли [8], имеющее вид:
дС дС дС , ч дС
-+ и-+ V— + м - wn)-=
т дх ду п' дг
=¿С УС)+
где С - концентрация пыли; V = (и^^) -вектор скорости ветра; (х, у, z) - координаты, направленные на восток, на север и вертикально вверх, соответственно; t - время; V — оператор «набла»; 5с = 1,3 - число Шмидта; 5 = 1 — мощность нормированного источника пыли; wn - скорость осаждения пыли; Мт — турбулентная вязкость воздуха, рассчитываемая по (к-в)-модели.
Использовалось RAS-моделирование турбулентности. Настроечные параметры (к-в)-модели, используемые при описании расчета в OpenFOAM, задавались равными: Сти = 0,03; С1 = 1.21; С2 =
= 1,92; C3 = 0; sigmak = 1; sigmaEps = 1,3; kappa = 0,41, согласно [22]. На поверхности земли в качестве граничного условия для энергии турбулентности использовалось kLowReWallFunction, а для скорости диссипации энергии турбулентности значение epsilonWallFunction.
Граничные условия для скорости ветра на входных боковых границах расчетной области выбирались из геострофического приближения с зависимостью от высоты по формуле слоя Экмана [23, 24]:
u = Ug - e"z (Ug cos (az) - Vg sin (az)), v = Vg - eaz (Vg cos (az)- Ug sin (az)), w = 0, a2 = 2Q sin (ф)/ц5, |s » 50 Па-с.
Здесь Ug = (dpdy)/f - геострофический перепад давления по оси Ох, Vg = (dp/dx )/ f - геострофический перепад давления по оси Оу. В некоторых расчетах, когда ветер дует с гор, принимался логарифмический профиль скорости ветра. На верхней границе для скорости задавался вектор скорости, для давления граничным условием является нулевой градиент. Для выходных боковых границ расчетной области задавались нулевыми градиенты скорости. На нижней границе — нулевая скорость ветра и нулевой градиент давления
Пыль моделировалась одной фракцией, для пылинок которой характерна скорость осаждения wn ~ 10-4 м/с. Это соответствует размеру частиц порядка от 2 до 10 микрон, способных перемещаться в потоке воздуха на заметные расстояния (для частички 2 мкм скорость осаждения порядка wn ~ 10-4 м/с). Большего размера пылинки далеко не улетают, частица размером 15 микрон имеет скорость осаждения wn ~ = 10-2 м/с.
Генерация сетки. Для выполнения расчетов с помощью пакета OpenFOAM создавалась трехмерная вычислительная сетка, расположенная в атмосфере над регионом Алагир-ского ущелья, в центре которого находится унальское хвостохранилище. Для достоверного воспроизведения высот склонов ущелья и горных пиков использовались топографические данные из открытых источников НАСА (http://terrain.party). Выбирались топографические данные для области 40 х 40 км, которые имеют горизонтальное разрешение около 30 м. Данные о высотах получены в виде файла изображения PNG в градациях серого цвета. Используя свободно распространяемую программу heightmap2stl (https://sourceforge.net/projects/ heightmap2stl/), изображение конвертируется в STL-Binary формат, а затем при помощи программы Meshmixer (http://www.meshmixer.com/) конвертируется в формат STL-ASCII, который является стандартом задания 3D поверхности в пакете OpenFOAM. Полученная топография для этой области сглаживалась на линейных масштабах порядка 100 м для уменьшения влияния мелких подробностей топографии. Для сглаживания использовались возможности программы Meshmixer.
При помощи программы BlockMesh генерировалась базовая сетка в виде прямоугольного параллелепипеда с квадратным основанием, для нее задавались размер и количество базовых точек. Топографические размеры области принимались 36 х 36 км с высотой 7 км, количество точек для грубой сетки было 60 х 60 х 21, для промежуточной сетки (на которой и выполнялись все расчеты) — 80 х 80 х 28, для подробной сетки — 100 х 100 х 35.
Для ограничения базовой сетки снизу высотами поверхности из STL-файла высот использовалась программа snappyHexMesh, входящая в пакет OpenFOAM. Стандартные настройки этой программы незначительно модифицировались, чтобы обеспечить правильное отражение в сетке значительных углов поверхности, и задавался минимальный объем ячейки генерируемой сетки, чтобы полученная сетка не содержала ячеек слишком маленького объема, что сказывается на шагах по времени при проведении расчета. После этого проводился расчет.
Проверка качества сетки. Вычислительная сетка является важным параметром проводимых вычислений, поэтому ее проверка — это необходимая часть исследования, подтверждающая реалистичность результатов и их повторяемость. Для выполнения быстрых расчетов можно рекомендовать грубую сетку, но с числом узлов не меньше чем 60 х 60 х 21. Расчет на такой сетке до 14400 с занимает порядка 5 часов времени работы одного ядра персонального компьютера Intel(R) Core(TM) i7—3770 CPU 3.4GHz и может быть ускорен использованием параллельного расчета. Работа такого компьютера в течение 5 часов является вполне доступным ресурсом для множества исследователей. Расчет на более подробных сетках занимает большее время, но продолжает оставаться доступным с точки зрения стоимости вычислений. Предпочтительнее анализировать результаты, полученные на более подробных сетках, которые разрешают большее количество деталей течения. Алгоритм проверки здесь не приводится, он близок к алгоритму проверки в работе [6] с помощью трех использованных сеток и выполнен в соответствии с рекомендациями, описанными в работах [4, 5].
Результаты исследований
Для скорости ветра было проведено 16 расчетов с различными направлениями скорости ветра, задаваемой на верхней границе расчетной области: северное, северо-северо-западное, северо-западное и так далее с шагом 22,5 градусов. В каждом варианте вычислялось поле ЗВ от модельного источника, расположенного в районе Унальского хвостохранилища. Источник попадает в одну ячейку расчетной области, непосредственно примыкающую к поверхности, высотой около 100 метров. Модуль скорости на верхней границе расчетной области во всех расчетах задавался равным 10 м/с. Заметим, что все измеренные скорости на высоте 50 м ниже 10 м/с (рис. 1). Для каждого варианта использовалась сетка размером 80 х 80 х 28 (размер, задаваемый до включения рельефа поверхности). Все рассчитанные варианты складывались друг с другом с различным весовым множителем. Для каждого направления ветра этот множитель выявлялся на основании розы ветров,
из которой определялось количество повторений данного направления внешнего ветра относительно общего числа наблюдаемых ветров.
Для оценки распределения ЗВ по склонам Алагирского ущелья была выбрана роза ветров, изображенная на рис. 1. Восточный и юго-восточный ветра имели максимальный весовой коэффициент, равный 0,1422, весовой коэффициент северо-западного ветра равнялся 0,0577. В результате осреднения всех расчетов с весовыми множителями, полученными из розы ветров, получается среднее поле концентрации ЗВ, изображенное на рис. 2. Для удобства анализа на рис. 2 показана карта высот поверхности расчетной области, экспортированная из расчетной сетки OpenFOAM. Показаны значения 0 < с < 0,1 С^^ где Стах - конЦентрация ЗВ вблизи источника. Расстояния показаны в метрах.
Данные на рис. 2 позволяют получить представление о приблизительном значении концентрации ЗВ в каждой точке склонов ущелья. Это возможно
Рис. 2. Результаты расчетов. Белым кружком отмечена чаша хвостохранилища: а — карта высот поверхности, экспортированная из OpenFOAM; б — вычисленная концентрация ЗВ Fig. 2. Result of simulations. The tailings impoundment is marked with a white circle: a — the surface elevation map exported from OpenFOAM; b — the simulated field of pollutant concentration.
а
после нормировки на наблюдаемые в полевых условиях значения ЗВ. Для нормировки можно выбрать измерения, выполненные в одной из точек поверхности ущелья или в нескольких точках.
Вычисленная осредненная концентрация ЗВ в поселке Унал равнялась Си = 1,8 • 10-6, соответствующая ей концентрация свинца, полученная при измерениях ЗВ: С^ = 438 мг/кг. Эта информация позволяет нормировать полученное среднее поле ЗВ. Если из поля уже нормированной концентрации ЗВ вычислить концентрацию свинца в поселке Архон, получится Сд0™ = = 202 мг/кг (среднее по измерениям есть С^ = 356 мг/кг), а в поселке Зин-цар О0™ = 414 мг/кг (среднее по измерениям есть С^ = 445 мг/кг). Сопоставление с данными полевых измерений для с. Унал и с. Зинцар показывает удовлетворительное согласие. Вычисления сведены в табл. 2.
В измерениях табл. 2 для с. Архон указано два значения для свинца: 356 и в скобках 624, первое значение получено осреднением измерений (табл. 1)
без учета одиночного очень большого значения 3117 для свинца. Аналогично для цинка 522 (795), первое значение получено осреднением измерений (табл. 1) без учета значения 2979. Эти одиночные большие значения измеренной концентрации ЗВ соответствуют некоторым локальным особенностям местности, например, наносам дождевыми потоками в небольшую низину. Авторы считают такие значения не вполне характеризующими среднюю концентрацию ЗВ, поэтому они исключались при проведении вычислений.
Селение Архон расположено в одном из небольших ущелий, примыкающих к Алагирскому ущелью. Вычисленная концентрация ЗВ в Ала-гирском ущелье, напротив выхода из ущелья селения Архон, составляет С'А = 1,45-10-7, что в 5,7 раза меньше, чем в селении Архон. Вентиляция этой небольшой долины довольно низкая и ЗВ накапливаются в ней, что видно на рис. 2, как локальный максимум юго-западнее хвостохранилища. Таких локальных максимумов просматрива-
Таблица 2
Сопоставление вычислений по модели с полевыми измерениями
Table 2
Comparison of simulations with field measurements
Показатель Унал Архон Зинцар Мизур (прогноз) Алагир (прогноз)
Ненормированные вычисления 1,8-10-6 8,3-10-7 1,7-10-6 5,7-10-8 3,0-10-7
По свинцу (Pb), мг/кг
Полевые измерения, ПДК 32 438 356 (624) 445 - -
Прогноз по измерениям в с. Унал 438 202 414 14 73
Прогноз по измерениям в с. Архон 772 356 729 24 129
Прогноз по измерениям в с. Зинцар 471 217 445 15 78
По цинку (Zn)
Полевые измерения, ПДК 23 470 522 (795) 421 - -
Прогноз по измерениям в с. Унал 470 217 444 15 78
Прогноз по измерениям в с. Архон 1132 522 1069 35 189
Прогноз по измерениям в с. Зинцар 446 205 421 29 152
ется несколько в разных точках рис. 2. Все они соответствуют особенностям локальной топографии, в которых пыль плохо выветривается и ЗВ накапливаются.
Количество пыли, оседающей в таких областях, зависит от времени расчета. При расчете более чем на 4 часа концентрация ЗВ в таких местах может иметь большее значение. В использованной методике все расчеты ограничены 4 часами. Приблизительно через 2 часа расчета поле скоростей ветра и турбулентности перестает меняться во времени, еще через час перестает существенно меняться и поле концентрации ЗВ. Однако в некоторых расчетах поле концентрации ЗВ медленно растет в нескольких приземных низинах или на подветренной стороне гор - там пыль плохо выветривается и накапливается.
Используемая методика может быть применима для приблизительной оценки поля ЗВ от источника пыли, имеющего неопределенную мощность, действующего продолжительное время и расположенного в горной или на холмистой территории, где существенно влияние топографических особенностей местности. В результате серии расчетов получается осредненное поле концентрации ЗВ, которое нужно нормировать на результаты полевых измерений, полученных хотя бы для одной точки поверхности.
Обсуждение результатов
исследования
Топография поверхности оказывает существенное влияние на пространственное распределение и распространение мелкой пыли. В горных ущельях и вблизи горных склонов закономерности распространения ЗВ и мелкодисперсной пыли существенно отличаются от закономерностей, известных
для равнинных территорий. В горных ущельях ЗВ в основном распространяются вдоль ущелья, ориентация которого и определяет направление движения ЗВ, и мало зависит от дующего над ущельем ветра.
Скорость ветра над вершинами гор меняется мало и сохраняет значение, определяемое ветром, заданным на верхней границе. Скорость ветра в ущельях ниже вершин гор имеет сложную структуру. Если его скорость над ущельем имеет перпендикулярную оси ущелья компоненту, то в ущелье возникает крупное вихревое течение с горизонтальной осью, и крупные вихри могут аккумулировать в себе ЗВ, которые слабо распространяются за пределы этих вихрей. В ущелье за одиночно стоящей горой также возникает вихревое течение, которое обладает сложной трехмерной структурой, ЗВ в нем накапливаются меньше, чем в заполняющих ущелья вихрях с горизонтальной осью вращения.
Отметим существенное влияние еще одного фактора. Наличие горных вершин существенно меняет пространственное распределение турбулентной вязкости воздуха. Пространственное распределение турбулентной вязкости уникально для каждого из шестнадцати рассчитанных вариантов, но в целом области максимумов этой величины вытянуты вдоль направления ветра с подветренной стороны наибольших горных вершин. Ниже по ветру от горной вершины максимальные значения турбулентной вязкости образуют длинный шлейф от вершины и до выходной границы расчетной области, высота этого шлейфа может несколько превышать высоту горной вершины, ширина его сравнима с шириной вершины. Зачастую таких турбулентных следов в расчетной области наблюдается три или более от разных вершин, и рас-
положены они на разных расстояниях друг от друга. В пространственных пределах этого турбулентного следа от горной вершины турбулентность значительна, и эти области существенно способствуют переносу ЗВ. Источник ЗВ в нескольких расчетах находится вблизи одной из таких областей повышенной турбулентности, которая возникает от расположенной в тысяче метрах вершины (отметка 2571 м). В области повышенной турбулентности ЗВ могут подниматься на высоту, сопоставимую с высотой горной вершины, и далее переноситься на значительные расстояния на большой высоте, при этом концентрация ЗВ вблизи поверхности под этим шлейфом остается малой.
Например, наиболее повторяющийся ветер над районом хвостохра-нилища — это юго-восточный ветер. В случае такого направления ветра на верхней границе расчетной области скорости ветра в ущелье непосредственно над хвостохранилищем являются юго-западными, однако ЗВ переносится в северо-восточном направлении только до начала вихревого шлейфа расположенной рядом горы — приблизительно в 3-х км ниже по ущелью. В этом турбулентном шлейфе они поднимаются вверх на наветренный склон ущелья, выше гряды гор высотой около 1800 м, и далее переносятся юго-восточным ветром. При этом ветер внутри ущелья дует вдоль оси ущелья и переносит часть ЗВ в северо-восточном направлении. В этой точке поток ЗВ разделяется на два, один из которых перемещает ЗВ над поверхностью в виде приподнятого факела на северо-запад, другой — вдоль по ущелью вблизи поверхности на северо-восток в сторону города Алагир.
Не часто повторяющиеся северные ветра и эффективное выветривание
турбулентными потоками приводит к тому, что в поселке Мизур, расположенном вверх по ущелью от хвостох-ранилища, прогнозируемые моделью концентрации ЗВ не превышают ПДК (табл. 2). Часто дующие над хвостох-ранилищем юго-восточные и юго-юго-восточные ветра переносят ЗВ вдоль ущелья и достигают г. Алагира, где прогнозируются значения концентрации ЗВ, превышающие ПДК (табл. 2).
Модельное среднее распределение концентрации ЗВ (рис. 2) показывает несколько мест, где вблизи поверхности локально очень высокая концентрация ЗВ. Эти места видны на рис. 2 как пятна повышенной концентрации ЗВ. Они являются результатом топографических особенностей подстилающей поверхности. Это низины и небольшие горные долины. Пыль осаждается в этих местах и слабо рассеивается как ветром, так и турбулентностью. Одно из таких мест расположено вблизи с. Архон, где измерения показывают большее значение концентрации ЗВ, чем расчетные значения. Хотя модель также показывает локальный максимум концентрации ЗВ в этой долине, но не настолько большой.
Необходимо рассмотреть некоторые допущения использованной модели с точки зрения физики атмосферы. Первое, что нужно отметить, это использование в модели предположения о несжимаемости воздуха [21]. Такое предположение может быть справедливым только при сплошной облачности. Влияние сжимаемости атмосферного воздуха сказывается при учете изменений температуры на высотах от 2 км и выше [23]. Однако это требует задания поля температуры, что вводит в систему решаемых уравнений еще несколько существенных параметров, например альбедо поверхности. Кроме того, температуру над склонами гор
можно вычислить только для конкретного времени суток с учетом сезона, времени года и высоты Солнца над горизонтом, переизлучения солнечной радиации склонами ущелья, а также с учетом затенения облаками поверхности земли. Турбулентность также зависит от стратификации атмосферы и от времени суток. Учет изменения температуры воздуха потребует более длительных вычислений, моделирующих все летнее время в ущелье.
Еще одно использованное предположение состоит в том, что не удаляются из расчета пылинки, достигшие поверхности земли. Это означает, что все достигшие поверхности пылинки участвуют во вторичном переносе: достигнув поверхности, они сдуваются ветром и остаются в атмосфере. На самом деле только часть пылинок после осаждения вновь отрывается от поверхности. Способность отрыва пылинки от поверхности зависит от свойств поверхности, которая в наших расчетах задавалась идеализированной, одинаковой для поросших лесом склонов, для скал и для пастбищ.
В работе не используются какие-либо критерии выноса частиц пыли с поверхности хвостохранилища. Экспериментально определенные условия
отрыва пыли от поверхности [1 — 6] могут быть использованы на относительно коротких временных интервалах, когда известны свойства вещества в хвостохранилище, его площадь и мощность источника ЗВ, для этих случаев имеет смысл и более подробное моделирование физических процессов, связанных с учетом температуры.
Предложенная модель является достаточно простой и может применяться для оценки многолетних выбросов ЗВ и для оценки мощности источника по результатам полевых измерений, на которые нормируются расчетные концентрации ЗВ.
Заключение
Представленная модель может оказаться полезной для оценки концентрации пыли и описания картины ее распределения на склонах горных ущелий. В иерархии математических моделей различной сложности представленная модель может рассматриваться среди простых и работоспособных инструментов. Сопоставление результатов расчетов с данными полевых измерений показывает приемлемую точность, достаточную для многих практических приложений.
СПИСОК ЛИТЕРАТУРЫ
1. EPA. Iron King Mine and Humboldt Smelter. 2010. URL: https://www.epa.gov/ air-emissions-factors-and-quantification/ap-42-compilation-air-emissions-factors (дата обращения 09.10.2019).
2. Gillies J.A., 2013. Fundamentals of aeolian sediment transport: dust emissions and transport — near surface. In: Shroder, J. (Editor in Chief), Lancaster, N., Sherman, D.J., Baas, A.C.W. (Eds.). Treatise on Geomorphology. Academic Press, San Diego, CA, vol. 11, Aeolian Geomorphology, pp. 43—63.
3. Kok J.F., Mahowald N.M., Fratini G., Gillies J.A., Ishizuka M., Leys J.F., Mikami M., Park M.-S., Park S.-U., Van Pelt R.S., and Zobeck T.M. An improved dust emission model — Part 1: Model description and comparison against measurements. Atmos. Chem. Phys., 2014, 14, pp. 13023—13041, DOI: 10.5194/acp-14—13023—2014.
4. Roache P.J. Verification and Validation in Computational Science and Engineering. Albuquerque, NM: Hermosa Publishers; 1998.
5. Stern F, WiLson RV, CoLeman HW, Paterson EG. Comprehensive approach to verification and validation of CFD simulations Part 1: methodology procedures // J. Fluids Eng., 2001, 123, pp. 793-802.
6. Stovern M. at al. Simulation of windblown dust transport from a mine tailings impoundment using a computational fluid dynamics model // AeoLian Res. 2014, 14, pp. 75-83. doi:10.1016/j.aeoLia.2014.02.008.
7. Turpin C., Harion J.L. Effect of the topography of an industrial site on dust emissions from open storage yards // Environ Fluid Mech. 2010, 10, pp. 677. https://doi.org/10.1007/ s10652 — 010—9170—3
8. Алоян А.Е. Динамика и кинетика газовых примесей и аэрозолей в атмосфере / Курс лекций. — М.: ИВМ РАН, 2002. — 201 с. — ISBN 5 — 901854—05—5.
9. Берлянд М.Е. Современные проблемы атмосферной диффузии и загрязнения атмосферы. — Л., Гидрометеоиздат, 1975. 448 с.
10. Дегтярева Т.В., Лиховид А.А., Лысенко А.В., Караев Ю.И. Региональные структуры миграции химических элементов в ландшафтах Северного Кавказа // Устойчивое развитие горных территорий. — 2018. - Т. 10. -№ 4. - С. 481 — 492.
11. Лолаев А.Б., Бадтиев Б.П., Бутюгин В.В., Бадоев А.С. Определение консолида-ционных характеристик хвостов намывных техногенных месторождений // Устойчивое развитие горных территорий. — 2017. - Т. 9. -№ 4. - С. 355 — 361.
12. Голик В.И., Соболев А.А., Дзапаров В.Х., Харебов Г.З. Перспективы разработки месторождений Садона // Устойчивое развитие горных территорий. — 2018. - Т. 10. -№ 3. - С. 420—426.
13. Лолаев А.Б., Гурбанов А.Г., Дзебоев С.О., Илаев В.Э. Загрязнение прилегающих территорий в районе деятельности Садонского свинцово-цинкового комбината (Республика Северная Осетия-Алания) // Вестник ВНЦ РАН. - 2017. Том 6. №2. -С. 177—180.
14. Атаева А.Ю., Свердлик Г.И. Разработка экспериментальной установки для исследования процессов очистки пылегазовых выбросов предприятий предгорной зоны // Устойчивое развитие горных территорий. — 2017. - Т. 9. -№ 1. - С. 92 — 97.
15. Петров Ю.С., Соколов А.А., Раус Е.В. Математическая модель оценки техногенного ущерба от функционирования горных предприятий // Устойчивое развитие горных территорий. — 2019. - Т. 11. -№ 4. - С. 554—560.
16. Петров Ю.С., Хадиков М.К. Математическая модель и алгоритм анализа экологической устойчивости транспортной системы горного региона // Устойчивое развитие горных территорий. — 2018. - Т. 10. -№ 3. - С. 427 — 435.
17. Teixeira M.A.C., Kirshbaum D.J., Olafsson H., Sheridan P.F. and Stiperski I., eds. The atmosphere over mountainous regions. Frontiers in Earch Science. — Frontiers Media SA, Lausanne, Switzerland, 2016, pp. 162. ISBN 9782889450169.
18. Chow F.K., De Wkkker S.F.J., Snyder B.J. Mountain Weather Research and Forecasting, Recent Progress and Current Challenges. — Springer-VerLag Berlin Heidelberg, 2013, pp. 750. ISBN 978—94—007 — 4097 — 6.
19. Pathirana A., Herath S., and Yamada T. SimuLating orographic rainfaLL with a Limited-area, non-hydrostatic atmospheric modeL under ideaLized forcing // Atmos. Chem. Phys. 2005, 5, pp. 215—226
20. Lehner M., Whiteman C.D., Dorninger M. Inversion BuiLd-Up and CoLd-Air Outflow in a SmaLL ALpine SinkhoLe //Boundary-Layer MeteoroLogy, 2017, 163, pp. 497—522, DOI 10.1007/s10546—017—0232—7
21. Issa R.I. SoLution of the impLicitLy discretised fluid flow equations by operatorspLitting // JournaL of ComputationaL Physics, 1986, V. 62, 1, pp. 40—65. https://doi. org/L0.1016/0021—9991(86)90099—9
22. Hargreaves D.M. and Wright N.G. n the use of the k-EpsiLon model in commercial CFD software to model the neutral atmospheric boundary Layer // J. of wind engineering and industrial aerodymanics, 2007, 95, pp. 355 — 269
23. Shapiro A., Fedorovich E. Unsteady convectively driven flow along a vertical plate immersed in a stably stratified fluid // J. Fluid Mech. 2004, vol. 498, pp. 333—352. DOI 10.1017/S0022112003006803
24. Tan Z.M., Farahani M.M. An analytical study of the diurnal variations of wind in a semigeostrophic ekman boundary layer model // Boundary-Layer Meteorology, 1998, 86, pp. 313—332. https://doi.org/L0.1023/A:1000694732459 EES
REFERENCES
1. EPA. Iron King Mine and Humboldt Smelter. 2010. URL: https://www.epa.gov/ air-emissions-factors-and-quantification/ap-42-compilation-air-emissions-factors (data obrashcheniya 09.10.2019).
2. Gillies J.A., 2013. Fundamentals of aeolian sediment transport: dust emissions and transport near surface. In: Shroder, J. (Editor in Chief), Lancaster, N., Sherman, D.J., Baas, A.C.W. (Eds.). Treatise on Geomorphology. Academic Press, San Diego, CA, vol. 11, Aeolian Geomorphology, pp. 43—63.
3. Kok J.F., Mahowald N.M., Fratini G., Gillies J.A., Ishizuka M., Leys J.F., Mikami M., Park M.-S., Park S.-U., Van Pelt R.S., and Zobeck T.M. An improved dust emission model Part 1: Model description and comparison against measurements. Atmos. Chem. Phys., 2014, 14, pp. 13023—13041, DOI: 10.5194/acp-14—13023—2014.
4. Roache P.J. Verification and Validation in Computational Science and Engineering. Albuquerque, NM: Hermosa Publishers; 1998.
5. Stern F, Wilson RV, Coleman HW, Paterson EG. Comprehensive approach to verification and validation of CFD simulations Part 1: methodology procedures. J. Fluids Eng., 2001, 123, pp. 793 — 802.
6. Stovern M. at al. Simulation of windblown dust transport from a mine tailings impoundment using a computational fluid dynamics model. Aeolian Res. 2014, 14, pp. 75 — 83. doi:10.1016/j.aeolia.2014.02.008.
7. Turpin C., Harion J.L. Effect of the topography of an industrial site on dust emissions from open storage yards. Environ Fluid Mech. 2010, 10, pp. 677. https://doi.org/10.1007/ s10652 — 010—9170—3.
8. Aloyan A.E. Dinamika i kinetika gazovyh primesej i aerozolej v atmosfere [Dynamics and kinetics of gas contaminants and aerosols in the atmosphere]. Kurs lekcij. Moscow: IVM RAN, 2002. 201 p. ISBN 5-901854-05-5. [In Russ]
9. Berlyand M.E. Sovremennye problemy atmosfernoj diffuzii i zagryazneniya atmosfery [Modern problems of atmospheric diffusion and air pollution]. Leningrad: Gidrometeoizdat, 1975. 448 p. [In Russ]
10. Degtyareva T.V., Lihovid A.A., Lysenko A.V., Karaev Yu.I. Regional patterns of chemical elements migration in the landscapes of the North Caucasus. Sustainable Development of Mountain Territories. 2018. T. 10. no. 4. pp. 481 — 492. [In Russ]
11. Lolaev A.B., Badtiev B.P., Butyugin V.V., Badoev A.S. Determination of consolidation characteristics of tailings of anthropogenic alluvial deposits. Sustainable Development of Mountain Territories. 2017. T. 9. no. 4. p. 355 — 361. [In Russ]
12. Golik V.I., Sobolev A.A., Dzaparov V.H., Harebov G.Z. Prospects for the development of Sadon deposits. Sustainable Development of Mountain Territories. 2018. T. 10. no. 3. pp. 420 — 426. [In Russ]
13. Lolaev A.B., Gurbanov A.G., Dzeboev S.O., Ilaev V.E. Contamination of adjacent territories in the area of the Sadon Lead-Zinc Plant (Republic of North Ossetia-Alania). Vestnik VNC RAN. 2017. Tom 6. no. 2. pp. 177—180. [In Russ]
14. Ataeva A.Yu., SverdLik G.I. DeveLopment of experimentaL instaLLation for research of processes of cLeaning dust and gas emissions of enterprises in the foothiLL zone Sustainable Development of Mountain Territories. 2017. T. 9. no. 1. pp. 92—97. [In Russ]
15. Petrov Yu.S., SokoLov A.A., Raus E.V. MathematicaL modeL for estimating technogenic Losses from the operation of mining enterprises. Sustainable Development of Mountain Territories. 2019. T. 11. no. 4. pp. 554—560. [In Russ]
16. Petrov Yu.S., Hadikov M.K. MathematicaL modeL and anaLysis aLgorithm of the transport system environmentaL sustainabiLity in mountain region. Sustainable Development of Mountain Territories. 2018. T. 10. no. 3. pp. 427 — 435. [In Russ]
17. Teixeira M.A.C., Kirshbaum D.J., OLafsson H., Sheridan P.F. and Stiperski I., eds. The atmosphere over mountainous regions. Frontiers in Earch Science. Frontiers Media SA, Lausanne, SwitzerLand, 2016, pp. 162. ISBN 9782889450169.
18. Chow F.K., De Wekker S.F.J., Snyder B.J. Mountain Weather Research and Forecasting, Recent Progress and Current ChaLLenges. Springer-VerLag BerLin HeideLberg, 2013, pp. 750. ISBN 978-94-007-4097-6.
19. Pathirana A., Herath S., and Yamada T. SimuLating orographic rainfaLL with a Limited-area, non-hydrostatic atmospheric modeL under ideaLized forcing. Atmos. Chem. Phys. 2005, 5, pp. 215 — 226.
20. Lehner M., Whiteman C.D., Dorninger M. Inversion BuiLd-Up and CoLd-Air Outflow in a SmaLL ALpine SinkhoLe. Boundary-Layer MeteoroLogy, 2017, 163, pp. 497—522, DOI 10.1007/s10546-017-0232-7.
21. Issa R.I. SoLution of the impLicitLy discretised fluid flow equations by operator-spLitting. JournaL of ComputationaL Physics, 1986, V. 62, 1, pp. 40—65. https://doi.org/10.1016/0021-9991(86)90099-9.
22. Hargreaves D.M. and Wright N.G. n the use of the k-EpsiLon modeL in commerciaL CFD software to modeL the neutraL atmospheric boundary Layer. J. of wind engineering and industriaL aerodymanics, 2007, 95, pp. 355—269.
23. Shapiro A., Fedorovich E. Unsteady convectiveLy driven flow aLong a verticaL pLate immersed in a stabLy stratified fluid. J. FLuid Mech. 2004, voL. 498, pp. 333—352. DOI 10.1017/S0022112003006803.
24. Tan Z.M., Farahani M.M. An anaLyticaL study of the diurnaL variations of wind in a semigeostrophic ekman boundary Layer modeL. Boundary-Layer MeteoroLogy, 1998, 86, pp. 313—332. https://doi.org/10.1023/A:1000694732459.
ИНФОРМАЦИЯ ОБ АВТОРАХ
Каменецкий Евгений Самойлович1,2 — докт. физ-мат. наук, профессор, главный научный сотрудник лаборатории математического моделирования, e-maiL: esk@smath.ru ORCID: 0000-0002-7105-3578;
Радионов Анатолий Анатольевич2 — канд. техн. наук, научный сотрудник лаборатории математического моделирования, ORCID: 0000-0002-6934-6873; Тимченко Василий Юрьевич2 — соискатель лаборатории математического моделирования;
Панаэтова Ольга Софокловна2 — аспирант лаборатории математического моделирования, ORCID: 0000-0003-2718-5108;
Свердлик Григорий Иосифович1 — докт. техн. наук, профессор;
1 Северо-Кавказский горно-металлургический институт (государственный технологический университет), Республика Северная Осетия-Алания, г. Владикавказ, Россия;
2 Южный математический институт Владикавказского научного центра РАН, Республика Северная Осетия-Алания, г. Владикавказ, esk@smath.ru.
Для контактов: Каменецкий Е.С., e-maiL: esk@smath.ru
INFORMATION ABOUT THE AUTHORS
Kamenetsky E.SЛ2, Professor, Dr. Sci. (Phys. Mathem.), Chief Researcher, esk@smath.ru; ORCID: 0000-0002-7105-3578;
Radionoff A.A.2, Cand. Sci. (Eng.), and Researcher, e-mail: aar200772@mail.ru ORCID: 0000-0002-6934-6873;
Timchenko V.Y.2 graduate student, Department of Mathematical Modeling, timchenko.vasily@ mail.ru;
Panaetova O.S.2, graduate student, ORCID: 0000-0003-2718-5108; Sverdlik G.I.1, Dr. Sci. (Eng.), Professor;
1 North Caucasian Institute of Mining and Metallurgy (State Technological University), Republic of North Ossetia-Alania, Vladikavkaz, Russia;
2 Southern Mathematical Institute of Vladikavkaz Scientific Center of Russian Academy of Sciences, Republic of North Ossetia-Alania, Vladikavkaz, Russia, esk@smath.ru; Corresponding author: Kamenetsky E.S., e-mail: esk@smath.ru.
Получена редакцией 26.05.2020; получена после рецензии 03.07.2020; принята к печати 10.10.2020. Received by the editors 26.05.2020; received after the review 03.07.2020; accepted for printing 10.10.2020.
