Физика Земли, 2021, № 5, стр. 254-266

1. ВВЕДЕНИЕ

Один из важных процессов добычи полезных ископаемых на горнорудных открытых карьерах связан с взрывной отбойкой горной массы. Вследствие дробления горной породы взрывом этот процесс сопровождается интенсивной эмиссией в атмосферу продуктов детонации и частиц и последующим их переносом ветром на значительные расстояния, что оказывает негативное воздействие на состояние окружающей среды [Patra et al., 2016]. Взрывные работы являются одним из основных источников твердых частиц на карьерах [Бересневич и др., 1990; Roy et al., 2010; Адушкин и др., 2000; Fişne et al., 2011; Monjezi et al., 2009]. Существенный вклад в объем выбрасываемых частиц вносит использование взрывных технологий, связанных с почти одновременным подрывом системы линейно расположенных на борту карьера зарядов ВВ (взрывчатых веществ), в том числе и достаточно мощных. По оценкам, более 50% общего количества твердых частиц, выбрасываемых в приземный слой атмосферы при разработке полезных ископаемых на открытых карьерах, приходится на долю подобного вида множественных взрывов, так называемых массовых взрывов [Адушкин и др., 2010]. В 1959 г. на Высокогорском ГОКе (горно-обогатительном комбинате) в скважины общей длиной 10 км поместили 190 т аммонита [Кутузов, 2008]. Взрывом было отбито более 400 тыс. т руды. В декабре 2020 г. на Лебединском ГОКе в карьер заложили в 1643 скважины 1943 тонны промышленного взрывчатого вещества¹¹. Количество взрывов с суммарной массой более 500 т в период 1997–1999 гг. на Лебединском ГОКе составило в среднем 25 взрывов в год [Адушкин и др., 2000]. В 2015 г. на Лебединском ГОКе был проведен массовый взрыв с расходом ВВ 3155.2 т, а объем отбитой горной массы составил 2855 тыс. м³ [Угаров и др., 2017]. При этом массовые взрывы на карьерах проводятся один раз в месяц. Грубая оценка годового объема отбитой горной массы на Лебединском ГОКе дает значение 34 × 10⁶ м³, а оценка массы взрывчатых веществ, использованных за год, составляет ~37.9 × 10⁶ кг. По данным статистических отчетов общие выбросы твердых частиц в атмосферу от всех регистрируемых источников на Лебединском ГОКе составили в 2017 г. – 6.28 × 10⁶ кг, а в 2018 г. – 7.05 × 10⁶ кг [Звягинцева и др., 2020]. Соответствующим образом увеличилась и доля массовых взрывов.

Несмотря на повсеместное использование взрывов при разработке месторождений открытым способом, эмиссия твердых частиц в атмосферу при взрывном дроблении грунта – процесс еще недостаточно исследованный. Связано это с неполнотой данных и ненадежностью доступных тестов [Sairanen et al., 2017]. В этих условиях детализацию процессов, сопровождающих развитие пылевого облака, можно существенно дополнить методами численного моделирования. Важную информацию о распространении частиц внутри карьера и вне его приносят исследования переноса частиц локальными потоками ветра с учетом возмущений, создаваемых сложным рельефом карьера [Silvester et al., 2009; Joseph et al., 2018; Torno et al., 2010]. Обычно в этих случаях численно решаются уравнения Навье–Стокса с использованием тех или иных моделей турбулентного переноса. Всесторонний обзор моделей переноса пыли, разработанных для предсказания распространения частиц в процессе добычи полезных ископаемых открытым способом, представлен в работе [Reed, 2005]. Пылевое облако взрыва в этих расчетах моделируется экспериментально определенным потоком частиц из заданных областей инжекции, которые в случае массовых взрывов представляют собой протяженные горизонтально расположенные параллелепипеды. В подобных моделях игнорируется образование горячего, заполненного продуктами детонации, воздухом и пылью объема в окрестности эпицентра взрыва и последующее конвективное всплытие горячего газопылевого облака в стратифицированной атмосфере. Однако эти процессы становятся определяющими для переноса пыли, генерируемой мощными взрывами. Верхняя граница газопылевого облака поднимается на высоты, превышающие уровень дневной поверхности карьера. Последующее распространение пылевых частиц определяется господствующими на этих высотах ветрами. Корректное описание переноса ветром частиц после того, как облако достигло верхнего положения, требует определения характеристик облака к этому моменту, которые, в свою очередь, определяются всем комплексом газодинамических процессов, сопровождающих формирование и подъем газопылевого облака.

Анализ динамики и значений характерных величин газопылевого облака одиночного тротилового взрыва в диапазоне масс 1–1000 т ТНТ (тринитротолуола) приведен в работе [Khazins et al., 2020]. Численное моделирование газодинамических процессов в этой работе разбито на две стадии, для каждой из которых используются различные приближения и методы. Первая стадия характеризуется быстропротекающими (сотни миллисекунд) процессами, к которым можно отнести генерацию детонационных и ударных волн, инициированных подрывом ВВ, поступление частиц грунта различного размера из формирующегося кратера, расширение продуктов взрыва и их перемешивание с воздухом и пылью в процессе торможения продуктов, понижение давления до атмосферного в окрестности точки подрыва и формирование горячего всплывающего газопылевого объема (термика). Предполагается, что взрыв происходит на поверхности, влияние фактора, связанного с заглублением ВВ в реальных условиях, не учитывается. Вторая стадия, длительность которой составляет несколько минут, характеризуется относительно медленным подъемом термика в стратифицированной атмосфере под действием сил плавучести, его расширением и увлечением окружающего воздуха в процессе турбулентного перемешивания [Белоцерковский и др., 2000]. Первоначальная форма газопылевого облака разрушается, формируется тороидальный вихрь, а само газопылевое облако приобретает грибовидную форму. Подъем облака ограничивается высотой, на которой достигается гидростатическое равновесие с окружающим воздухом.

В настоящей работе рассматривается взаимодействие вихревых структур, сопровождающих подъем термиков, инициированных массовыми взрывами, а именно, системой линейно расположенных зарядов ВВ массой 1 т ТНТ при различной локализации зарядов. Как отмечено выше, для массового взрыва на Высокогорском ГОКе в 1959 г. в скважины общей длиной 10 км поместили 190 т ВВ. Если перевести это соотношение в расчете на 1 т ВВ, то получится, примерно, 1 т ВВ на 50–100 м. Это расстояние мы выбрали в качестве базового в наших расчетах.

В работе [Белоцерковский и др., 2000] представлен ряд работ, в которых рассмотрены пары, ансамбли термиков и вихревых колец, но среди них только результаты для пары приземных термиков можно использовать в какой-то мере для оценки процессов, развивающихся в рассматриваемых нами условиях. Реально число зарядов при массовых взрывах в карьерах существенно больше двух, и мы изучаем изменение характеристик облачных структур при увеличении числа линейно расположенных зарядов от трех, пяти до множества зарядов.

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

В разделе 2 мы описываем постановку задачи и численные модели: модель гидродинамического течения в рамках уравнений Эйлера для расчета ударно-волновых процессов, инициированных детонацией ВВ множественных взрывов, и процессов формирования горячих, заполненных пылью и продуктами детонации объемов (термиков); модель в рамках уравнений Навье–Стокса с вихревой вязкостью, определяемой методом крупных вихрей, для расчетов совместного всплытия термиков. В разделе 3 приводятся результаты расчетов ударно-волновых процессов для взрывов массой 1 т ТНТ. В разделе 4 исследуется эволюция облачных структур, подъем которых определяется линейно расположенными взрывами при разном расстоянии между ними и при различном их числе. Некоторые выводы представлены в последнем разделе.

2. ПОСТАНОВКА ЗАДАЧИ И МЕТОДЫ МОДЕЛИРОВАНИЯ

Пусть на гладкой горизонтальной поверхности полупространства вдоль прямой на одинаковом расстоянии друг от друга расположены полусферические заряды тротила массы М и пусть одновременно в центре каждой полусферы на границе ВВ и грунта инициируется детонация.

Для одиночного заряда массой 500 т ТНТ в работе [Khazins et al., 2020] приводятся результаты численного моделирования как формирования горячего квазисферического объема (термика), так и его всплытия. Для расчета всплытия системы термиков данных, полученных на стадии формирования одиночного термика, достаточно, но необходим ответ на вопрос о степени влияния взаимодействующих ударных волн, инициируемых детонацией зарядов ВВ, на начальные характеристики термиков при разных расстояниях между взрывами. С этой целью был выполнена серия расчетов начальной стадии развития взрыва зарядов тротила в трехмерной постановке по модели, аналогичной [Khazins et al., 2020], но без учета формирования кратера и выброса пыли. Модель реализована на основе гидродинамической программы SOVA [Shuvalov, 1999; Shuvalov et al., 1999]. Эта программа позволяет моделировать в рамках уравнений Эйлера сложные гидродинамические течения с сильными разрывами физических параметров и аккуратным описанием границ между областями с разными физико-химическими свойствами, в данном случае границ между грунтом, продуктами детонации и воздухом. Вихревые процессы, возникающие на рассматриваемой стадии развития газопылевого облака, достаточно успешно описываются турбулентным перемешиванием, определяемым схемной вязкостью.

Допускается использование различных видов уравнений состояния, аналитических и заданных в виде таблиц. В приведенных ниже расчетах использовались полуэмпирические уравнения состояния продуктов детонации [Андреев и др., 2004а; 2004б; Васильев, Ждан, 1981], таблицы термодинамических свойств воздуха [Кузнецов, 1965]. В модели учитывается действие силы тяжести и стратификация атмосферы, хотя на рассматриваемой начальной стадии они практически не влияют на результат.

Для описания процесса распространения детонационной волны в рамках эйлеровой схемы применялась следующая модель. После того как внутренняя энергия ε в какой-либо ячейке, через которую проходит ударный фронт, увеличивается до критического значения ε_cr, значение ее увеличивается на величину ${{Q\rho _{{dm}}^{0}} \mathord{\left/ {\vphantom {{Q\rho _{{dm}}^{0}} {{{\rho }_{{{\text{exp}}}}}}}} \right. \kern-0em} {{{\rho }_{{{\text{exp}}}}}}}$, где $\rho _{{dm}}^{0}$ – начальная плотность детонирующей смеси; ${{\rho }_{{{\text{exp}}}}}$ – плотность ВВ в текущий момент времени (которая превышает $\rho _{{dm}}^{0}$ за счет сжатия в ударной волне); Q – теплота реакции. Такая процедура обеспечивает сохранение полной энергии, правильные скорость детонации и параметры за фронтом детонационной волны.

По мере завершения волновых процессов формируется система термиков, схематично представленных на рис. 1.

Согласно результатам работы [Khazins et al., 2020], полученным для взрыва массой 500 т, давление ко времени завершения формирования одиночного термика снижается, практически до атмосферных значений, а температура в области термика уменьшается от значений, соответствующих распространению детонационной волны (~5000 К), до 1000–1500 К. В силу гидродинамического подобия эта температура не должна зависеть от массы заряда. Так, в таком же диапазоне меняется температура термика, определенная с помощью высокоскоростных радиометрических измерений, в экспериментах [Spidell et al., 2010] для зарядов с массой до 50 кг ТНТ. Результаты работы [Khazins et al., 2020] позволяют определить и радиус термика. Для взрыва массой 500 т ТНТ он равен ~140 м. Из гидродинамического подобия следует, что радиус термика должен быть пропорционален корню кубическому из массы заряда. Примерно такая зависимость была обнаружена в экспериментах [Kansa, 1997]:

Здесь $r_{T}^{0}$ – радиус термика в метрах, W − масса ВВ в тоннах. Значение радиуса термика $r_{T}^{0}$, определенного в работе [Khazins et al., 2020] для взрыва массой 500 т ТНТ, хорошо согласуется со значением, следующим из эмпирического соотношения (1).

Таким образом, результаты работы [Khazins et al., 2020], полученные для ВВ с массой 1–1000 т ТНТ, позволяют оценить начальную температуру и размер одиночного термика для любой массы взрыва. Для расчетов мы выбрали значение начальной температуры равным 1250 К. В расчетах предполагалось, что уравнение состояния – идеальный газ с показателем адиабаты γ = 1.4 как в области термика, так и в окружающем воздухе, а невозмущенная атмосфера описывается стандартным вертикальном распределением газодинамических параметров воздуха [ГОСТ…, 2004].

При численном моделировании всплытия термиков, как и в работе [Khazins et al., 2020], использовался разностный аналог полной системы уравнений Навье–Стокса для сжимаемой жидкости в гипозвуковом приближении [Затевахин и др., 1994; Хазинс, 2010; 2012], но не в осесимметричной постановке, а с учетом всех трех пространственных переменных (3D). Гипозвуковое приближение приводит к системе уравнений параболического типа, и модель по простоте численной реализации не уступает модели несжимаемой жидкости.

Многочисленные тесты показали применимость такого подхода для расчета всплытия термиков в широком диапазоне перепадов температур в термике и окружающем воздухе. Однако характеристики всплывающего термика существенным образом зависят от вихревой вязкости. Перенос энергии между крупномасштабными и мелкомасштабными вихрями моделировался в рамках численной схемы с помощью турбулентной вязкости, определяемой методом крупных вихрей LES (Large Eddy Simulation) [Lesier, Metais, 1996; Хазинс, 2010]. LES-метод занимает промежуточное положение между простейшими алгебраическими моделями, характеризующимися, например, постоянным числом Рейнольдса, и многокомпонентными системами дифференциальных уравнений в частных производных [Launder, Spalding, 1974; Anderson et al., 1984]. При использовании метода крупных вихрей кинематическая вихревая вязкость ${{\nu }_{T}}$ определяется как:

где: ${{C}_{{sm}}}$ − коэффициент пропорциональности (параметр Смагоринского); l – размер сеточной ячейки; $S$ – тензор деформации:

Здесь все размеры ‒ в метрах, скорость ‒ в м/с, ${{C}_{{sm}}}$ ‒ безразмерный параметр.

Сопоставление с результатами экспериментов показало [Khazins et al., 2020], что высота всплытия термика при одних и тех же начальных условиях зависит от значения параметра Смагоринского, которое, в свою очередь, зависит от величины шага разностной сетки и массы заряда. В качестве критерия для определения параметра Смагоринского было выбрано удовлетворительное согласие высоты верхней кромки H всплывающего термика и экспериментальных данных [Church, 1969], которые на момент времени 2 мин для тротиловых взрывов массой от 1 до 1000 т можно описать соотношением:

где: H ‒ в метрах; W − масса ВВ в тоннах. В работе [Khazins et al., 2020] было показано, что при 10 счетных интервалах на радиус термика зависимость параметра Смагоринского от массы взрыва ведет себя степенным образом:

Однако при 10 счетных интервалах на радиус термика время расчета оказывается очень большим при 3D-моделировании, и мы во всех вариантах использовали 5 счетных интервалов на радиус термика. В этом случае максимальный размер расчетной сетки составлял 280 × 150 × 150 счетных интервалов. Методом проб и ошибок было установлено, что при выбранном шаге разностной сетки для массы взрыва 1 т значение ${{C}_{{sm}}}$ должно быть примерно равным 3, а не 7.36, как это следует из (4). При этом качество 3D-расчетов почти не снижается по сравнению с расчетами на 10 точках по радиусу, а высота термика к двум минутам достигает, примерно, 500 м, как это и следует из соотношения (3) для взрыва массой 1 т.

Для определения геометрических характеристик всплывающего термика и визуализации течения в расчет включено уравнение переноса концентрации пылевой примеси:

где: ρ ‒ плотность газопылевой смеси (кг/м³); ${{c}_{k}}$ ‒ концентрация примеси (${{\rho }_{{ad}}} = \rho {{c}_{k}}$, ${{\rho }_{{ad}}}$ ‒ плотность примеси (кг/м³)); U ‒ вектор скорости (м/с); t ‒ время (с); D ‒ турбулентный коэффициент диффузии (м²/с), определяемый по значениям вихревой вязкости ${{\nu }_{{\text{T}}}}$ из условия, что эти величины связаны между собой турбулентным числом Шмидта $Sc = {{{{\nu }_{T}}} \mathord{\left/ {\vphantom {{{{\nu }_{T}}} D}} \right. \kern-0em} D}$. Турбулентное число Шмидта может меняться в широких пределах (от 0.1 до 2.5) в зависимости от рассматриваемой задачи [Gualtieri et al., 2017]. Мы выбрали значение, равное 1.

Для облегчения анализа распределений плотности примеси в качестве начальных данных для уравнения (5) внутри термика мы выбрали равенство плотности примеси ${{\rho }_{{ad}}}$ плотности воздуха у поверхности ${{\rho }_{0}}$, или, в безразмерных переменных, ${{{{\rho }_{{ad}}}} \mathord{\left/ {\vphantom {{{{\rho }_{{ad}}}} {{{\rho }_{0}}}}} \right. \kern-0em} {{{\rho }_{0}}}} = 1$. Вне термика ‒ ${{{{\rho }_{{ad}}}} \mathord{\left/ {\vphantom {{{{\rho }_{{ad}}}} {{{\rho }_{0}}}}} \right. \kern-0em} {{{\rho }_{0}}}} = 0$.

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

3. ВЛИЯНИЕ УДАРНО-ВОЛНОВЫХ ПРОЦЕССОВ НА ФОРМИРОВАНИЕ ТЕРМИКА

В рассматриваемой задаче нас интересует качественное влияние близко расположенных взрывов на форму образующегося термика, а не его точные характеристики. Рассмотрим предельный случай ‒ линейно расположенных зарядов ВВ бесконечно много. В этом случае задача упрощается, т.к. плоскости ABCD и A'B'C 'D' (рис. 1), расположенные посередине между эпицентрами взрывов, становятся плоскостями симметрии. В качестве начальных данных задавался полусферический объем ВВ (тротил с плотностью 1.65 г/см³ и скоростью детонации 7 км/с) радиусом 0.65 м, расположенный на жесткой горизонтальной поверхности. Масса заряда 1 т ТНТ. Детонация инициировалась в центре полусферы на границе ВВ и грунта. Расчетная сетка состояла из 290 × 290 × 290 ячеек в направлениях X (вдоль линии зарядов), Y (в поперечном горизонтальном направлении) и Z (в вертикальном направлении). Начальный размер ячейки составлял 0.8 см, по мере увеличения возмущенной области сетка несколько раз удваивалась.

На рис. 2а показана зависимость размеров облака взрыва по осям Х и Y от времени для случая, когда расстояние между взрывами очень большое, и они не влияют друг на друга. В этом случае размер облака по обеим осям примерно одинаков, сначала быстро растет, затем, достигнув значения около 18 м, почти не меняется (до тех пор, пока облако не начнет подниматься и перемешиваться с окружающим воздухом). У облака нет четкой границы (см. рис. 3), плотность воздуха плавно меняется – от минимальной до фоновой. При этом из-за развития неустойчивостей на начальной стадии разлета продуктов детонации [Khazins et al., 2020] форма облака не идеально сферическая, и наблюдается некоторое отличие размеров облака по осям X и Y. Для количественной оценки мы называли границей облака изолинию ρ/ρ₀ = 0.7, где ρ₀ – фоновая плотность воздуха. Для наглядности и удобства сопоставления результатов с описываемыми выше расчетами подъема термиков мы используем безразмерные переменные $\bar {x} = {x \mathord{\left/ {\vphantom {x {r_{T}^{0}}}} \right. \kern-0em} {r_{T}^{0}}}$, $\bar {y} = {y \mathord{\left/ {\vphantom {y {r_{T}^{0}}}} \right. \kern-0em} {r_{T}^{0}}}$, $\bar {z} = {z \mathord{\left/ {\vphantom {z {r_{T}^{0}}}} \right. \kern-0em} {r_{T}^{0}}}$ и $\bar {t} = {t \mathord{\left/ {\vphantom {t {{{t}_{0}}}}} \right. \kern-0em} {{{t}_{0}}}}$, где $r_{T}^{0}$ – окончательный размер облака для одиночного взрыва (в нашем случае 18 м), ${{t}_{0}} = {{r_{T}^{0}} \mathord{\left/ {\vphantom {{r_{T}^{0}} c}} \right. \kern-0em} c}$, с – скорость звука в воздухе.

При расстоянии между взрывами $d\bar {x} = 4$ и $d\bar {x} = 3$ различие между размерами облака взрыва по осям X и Y тоже невелико, не превышает 10–20%. Максимальное отличие наблюдается в моменты прохождения отраженной волны через границы облака, после многократного прохождения отраженных волн облако становится почти симметричным. На рис. 3 показано течение в некоторые моменты времени для одиночного взрыва и взрывов на расстоянии $d\bar {x} = 4$. Здесь также видно, что форма облака в обоих случаях примерно одинакова.

Таким образом, можно считать, что форма термика при расстояниях между взрывами 3–4$r_{T}^{0}$ и более такая же, как при одиночном взрыве, и в качестве начальных данных можно использовать сферическую форму термика.

4. РЕЗУЛЬТАТЫ РАСЧЕТОВ ВСПЛЫТИЯ СИСТЕМЫ ТЕРМИКОВ

Рассмотрим серию численных экспериментов, в которых расстояние между центрами термиков достаточно мало по сравнению с размерами генерируемых ими облачных структур. Максимальный радиус шапки облака ${{R}_{m}}$ для одиночного взрыва можно оценить по соотношению, установленному в работе [Khazins et al., 2020]:

Здесь $r_{T}^{0}$ и ${{R}_{m}}$ – в метрах, W − масса ВВ в тоннах. Следовательно, для взрывов с массой 1 т при расстояниях между эпицентрами, превышающими 18$r_{T}^{0}$, взаимодействие облачных структур не будет происходить. Пусть расстояние между эпицентрами составляет 4$r_{T}^{0}$.

Сопоставление результатов проведем со случаем всплытия одиночного термика, сформированного взрывом массой 1 т. На рис. 4 представлены характерные поля переменных в плоскости X0Z (рис. 1). Точно такими же полями характеризуется течение и в других направлениях в силу осевой симметрии. Поэтому наблюдающиеся на рис. 4а и 4в вихри представляют собой вихревые торы. Большая часть пассивной примеси попадает в шапку облака и служит удобной характеристикой ее размеров, а максимальная часть примеси со временем сосредотачивается в окрестности оси тороидального вихря. Верхняя кромка облака достигает высоты 500 м к 2 мин (рис. 4а), что соответствует значению, определяемому эмпирическим соотношением (3). Термик к этому моменту все еще характеризуется положительной плавучестью β, где $\beta = {{\left( {\rho - \rho {\text{*}}} \right)} \mathord{\left/ {\vphantom {{\left( {\rho - \rho {\text{*}}} \right)} {\rho {\text{*}}}}} \right. \kern-0em} {\rho {\text{*}}}}$, ρ – плотность воздуха, а ρ*(z) – плотность невозмущенной атмосферы на высоте z (рис. 4б). Основная часть объема термика пока еще легче окружающего воздуха, и наибольшее отклонение плотности воздуха внутри термика от плотности окружающего воздуха достигается в области тороидального вихря. Облако остывает по мере подъема вследствие расширения и подмешивания холодного воздуха, и его плотность становится больше плотности окружающего воздуха. Начало этого процесса наблюдается уже к 2 мин в виде затемненной области в верхней части облака (рис. 4б). Однако облако продолжает подниматься по инерции, сохраняя свою структуру (рис 4в), и к 4 мин плотность внутри облака повсеместно становится больше плотности окружающего воздуха (рис. 4г). Скорость подъема замедляется и, как следует из расчета, примерно к 5 мин облако достигает максимальной высоты.

На рис. 5 показано одновременное всплытие трех термиков, расположенных на одной прямой линии при расстоянии между центрами соседних термиков 4$r_{T}^{0}$. К двум минутам газодинамические потоки, генерируемые одиночными термиками, объединяются и течение схоже с подъемом одиночного термика, но чуть больших размеров (рис. 5а и 5б). Форма облака несколько асимметрична. В сечении, проведенном через область максимальных значений примеси (рис. 5в), асимметрия проявляется довольно слабо. К четырем минутам сохраняется влияние неоднородности начального пространственного энерговыделения, и в горизонтальном сечении шапка облака представляет собой вытянутый в направлении, перпендикулярном к линии расположения термиков в начальный момент времени, эллипс (рис. 5е). Максимальные скорости к 4 мин примерно в 1.5 раза выше, чем в одиночном термике к этому же моменту времени.

В случае одновременного всплытия пяти термиков, расположенных на одной прямой линии, картина течения несколько меняется и к 2 мин течение в плоскостях X = 0 и Y = 0 существенно отличается (рис. 4а и 4б). Если построить поля скорости и плотности примеси не только в представленном здесь сечении X = 0 (рис. 6б), но и в любом сечении $X = {{X}_{0}} < $ 80 м, то отличие будет крайне незначительное. Об этом свидетельствует и распределение плотности примеси в сечении Z = 375 м ‒ на расстояниях ±80 м от оси Y (рис. 6в). Следовательно, в области, ограниченной плоскостями X = ±80 м, развивается течение, которое по характеристикам близко к плоскому осесимметричному. Значение 80 м близко к 4$r_{T}^{0}$, но это совпадение, скорее всего, случайное. Со временем эта качественная характеристика пропадает, и к 4 мин плотность примеси в сечении облака плоскостью Z = 510 м приобретает вид, аналогичный случаю всплытия 3 термиков.

Рассмотрим предельный случай ‒ линейно расположенных термиков бесконечно много. В силу симметрии течение достаточно исследовать только в объеме, заключенном между плоскостями ABCD и A 'B 'C 'D ' (рис. 1), поставив на них граничное условие “жесткой стенки” (нормальные к плоскостям компоненты скорости равны нулю). Если для пяти термиков переход к плоскому осесимметричному течению представлял собой лишь эпизод в развитии трехмерных потоков, то в случае бесконечного числа термиков уже к двум минутам картина газодинамического течения становится одинаковой во всех плоскостях $X = {{X}_{0}},~~ - 2r_{T}^{0} \leqslant {{X}_{0}} \leqslant 2r_{T}^{0}$, и эта тенденция сохраняется на всем протяжении последующего всплытия термиков. Поэтому мы привели на рис. 7 характеристики течения только в сечении X = 0.

Распределения примеси в облаках, формирующихся при различной начальной конфигурации расположения термиков, имеют однотипный характер, и облака различаются, в основном, геометрическими размерами и концентрацией примеси. На рис. 8 представлена зависимость горизонтальных и вертикальных размеров облаков от времени для приведенных выше численных экспериментов и некоторых других. Кривая 1 на рис. 8 соответствует всплытию одиночного термика. К 5 мин облако достигает максимальной высоты (рис. 8а), соответствующей положению гидростатического равновесия, а зависимость от времени боковых размеров одинакова в плоскостях X = 0 и Y = 0 (рис. 8б и 8в) в силу осевой симметрии. Облако, генерируемое тремя термиками (рис. 8, кривые 2), поднимается несколько выше, процентов на 7–8, но горизонтальный размер в плоскости X = 0 увеличивается в полтора раза по сравнению с одиночным термиком. Наиболее любопытно в этом случае изменение горизонтального размера в плоскости Y = 0. Начиная со времени, примерно, 1.5 мин и до трех минут размер в этом направлении начинает уменьшаться, что свидетельствует о пространственной перестройке течения. Затем рост размера возобновляется. Отметим сразу, что подобная тенденция характерна и для случая всплытия 5 термиков (рис. 8, кривые 3). При всплытии 5 термиков максимальная высота подъема несколько уменьшается (примерно на 3–4%) по сравнению со всплытием 3 термиков, но все еще больше, чем для одиночного термика. Здесь важен сам факт понижения высоты с увеличением числа термиков, а не увеличения. Снижение максимальной высоты подъема облака, по-видимому, будет продолжаться по мере увеличения числа всплывающих термиков, и предельное значение, достигаемое в случае бесконечной системы термиков (кривые 4 (сплошная линия), рис. 8), составляет ~15% от высоты всплытия одиночного термика. Однако размер облака в плоскости X = 0 практически не меняется по сравнения со всплытием 5 термиков.

Весь предыдущий анализ относился к случаю, когда расстояние между центрами термиков составляет 4$r_{T}^{0}.$ С увеличением расстояния между термиками взаимодействие вихревых структур будет начинаться позднее, а далее процессы, в общем-то, не будут отличаться от описанных выше. Поэтому мы обратим внимание здесь только на изменение геометрических характеристик облака при увеличении расстояния между термиками в случае бесконечной системы термиков. Как и следовало ожидать, и максимальная высота подъема термиков, и горизонтальный размер в плоскости X = 0 начинают стремиться к значениям, соответствующим всплытию одиночного термика (рис. 8, кривые 5, 6, соответствующие расстояниям 8$r_{T}^{0}$ и 12$r_{T}^{0}$). Причем при расстоянии 12$r_{T}^{0}$ уже можно говорить о всплытии независимой системы термиков.

Так как для бесконечной системы термиков с расстояниями между центрами 4$r_{T}^{0}$ достаточно быстро развивается плоское осесимметричное течение, то можно сопоставить результаты со случаем, когда термик с самого начала задается в виде бесконечного цилиндра, заполненного горячим воздухом. Представляется естественным задать радиус сечения этого цилиндра из соображений сохранения энергии, для чего достаточно, чтобы объем такого цилиндра между плоскостями ABCD и A'B 'C 'D ' (рис. 1) был бы равен объему, занимаемому сферическим термиком, т.е. $r_{c}^{0} = {{~r_{T}^{0}} \mathord{\left/ {\vphantom {{~r_{T}^{0}} {\sqrt 3 }}} \right. \kern-0em} {\sqrt 3 }}$, где $r_{c}^{0}$ ‒ радиус цилиндра. Результаты подобного 2D-расчета приведены на рис. 6 точечной кривой. Результаты несколько отличаются от расчета с межцентровыми расстояниями между термиками 4$r_{T}^{0}$ (кривые 4), но для оценок вполне приемлемы. Качественно картина течения не отличается от рис. 7. Можно достичь и лучшей близости кривых, уменьшив $r_{c}^{0}$ в полтора раза. Результаты этого случая отображены на рис. 8 пунктирной кривой.

ЗАКЛЮЧЕНИЕ

На горнорудных открытых карьерах при взрывной отбойке горной массы довольно часто пылегазовое облако образуется в результате подрыва системы линейно расположенных на борту карьера достаточно мощных зарядов ВВ. Для оценки воздействия взрывов на окружающую среду вне пределов горнорудных карьеров необходимы данные об источнике загрязнения, который в рассматриваемом случае представляет собой вытянутую вдоль линии подрывов шапку грибовидного в сечении облака пыли. Основными параметрами являются высота подъема пылевого облака, размер шапки в направлении, перпендикулярном линии подрывов, масса содержащейся в шапке облака пыли и распределение массы пылевых частиц по размерам. Моделирование подъема пыли сопряжено со значительными трудностями, связанными с необходимостью расчета образования кратера и преобразования вылетающего из кратера вещества в дискретный набор частиц с заданным распределением по массе, которое известно лишь в редких случаях [Khazins et al., 2020]. В частности, для рассматриваемых в настоящей работе взрывов с массой 1 т ТНТ каждый таких данных нет. Более того, как показали расчеты [Khazins et al., 2020], масса воздуха в термике значительно больше массы содержащейся в нем пыли, поэтому наличие последней не должно, по крайней мере качественно, влиять на характер взаимодействия между термиками. Поэтому задача рассмотрена в упрощенной постановке без учета процессов, связанных с образованием и переносом пыли в области начального термика. Однако предполагалось, что в области начального термика существуют равномерно распределенные взвешенные в воздухе частицы пыли, двигающиеся со скоростями гидродинамического потока. Расчет распределения этих частиц в пространстве позволяет лучше визуализировать течение и достаточно просто определять геометрические характеристики шапки облака.

Численное моделирование показало, что при расстояниях между взрывами в несколько начальных радиусов термика ударные волны от соседних взрывов практически не влияют на форму и размер термика. Поэтому допустима постановка задачи, в которой начальные термики задаются таким же образом, как и для одиночного термика. В результате численного моделирования было показано, что с увеличением числа линейно расположенных взрывов при относительно небольшом расстоянии между эпицентрами (4$r_{T}^{0}$) высота верхней кромки пылевого облака сначала вырастает до значений, процентов на 5–10 превосходящих аналогичное значение для одиночного взрыва, а затем, начиная с 5 взрывов, начинает снижаться, достигая в предельном случае бесконечного числа зарядов значений, процентов на 15–20 меньших, чем высота кромки для одиночного взрыва. В случае бесконечного числа зарядов и расстояниями между ними 4$r_{T}^{0}$ довольно быстро трехмерное течение трансформируется в плоское осесимметричное течение в плоскостях, параллельных Y0X (рис. 1). С увеличением расстояния между термиками высота верхней кромки облака в случае бесконечного числа зарядов растет, достигая значений, соответствующих одиночному термику.

Размер облака в направлении, перпендикулярном линии подрыва зарядов, с увеличением числа зарядов растет, достигая двукратного значения при бесконечном числе зарядов по сравнению с аналогичным размером облака одиночного взрыва.