On the numerical solution to the problem of air shock wave propagation in mine workings.pdf В соответствии с «Правилами безопасности в угольных шахтах^» [1] угольные шахты являются опасными производственными объектами. В [2] перечислены виды аварий, характерных для шахт, где взрывы метана и угольной пыли являются наиболее тяжелыми. Опасными факторами взрыва являются избыточное давление, высокая температура, высокая концентрация токсичных продуктов взрыва и высокая скорость воздушных потоков, увлекающая за собой тяжелые предметы. Вследствие этих факторов взрывы часто приводят к человеческим жертвам и почти всегда к многомиллионным убыткам. Согласно статистике [3], с 1991 по 2017 год произошло 208 аварий с общим числом жертв 1600 человек, поэтому вопросы локализации взрыва, гашения воздушных ударных волн остаются до сих пор актуальными. Практическое определение границ зон поражения по разным факторам осуществляется с помощью специальных отраслевых методик. В России применялись или применяются три из них [4-6], при этом, действующая методика [6] основана на газодинамическом подходе с использованием системы классических уравнений газовой динамики, состоящей из уравнений сохранения массы, импульса и энергии. Такой подход обладает значительным преимуществом перед методом с использованием чисто эмпирических соотношений - он позволяет учитывать волновые эффекты и вводить в математическую модель новые физические особенности, оставаясь физически и математически корректным. Таковыми могут быть учет водяных и сланцевых заслонов [7-10], автоматических установок взрыволокали-зации и др. Наряду с преимуществами использование газодинамической модели привносит и свои сложности: скорость расчетов низкая из-за использования при расчетах численных методов по сравнению с расчетами с использованием эмпирических методик; модернизация математической модели влечет за собой существенное изменение программного кода, используемого для решения задач. Тем не менее модернизация математической модели является актуальной в связи с требованиями учета инженерных взрывозащитных сооружений. О численном решении задачи распространения воздушных ударных волн 109 В настоящей статье предлагается методика организации расчетных структур и декомпозиции расчетов при решении задач распространения ВУВ в сети горных выработок, которая упрощает модернизацию численных расчетов при изменении математической модели распространения ВУВ, ускоряет расчет за счет распараллеливания вычислений на современных многоядерных или многопроцессорных компьютерах. Математическая модель Математическая модель, описывающая распространение ВУВ в разветвленной сети горных выработок, описывается системой нестационарных газодинамических уравнений. Подробное описание модели рассматривается в [9]. Уравнения сохранения массы, импульса и энергии газа для прямолинейных участков выработок: ∂ρS + ∂ρuS = 0 ; (1) (2) (3) ∂t ∂x ’ ∂P fS∂P fuS + ^f- = 0; ∂t∂x ∂ρuS + ∂(ρu^2 + p)S = π + ∂S ; + - Z> X + p ; ∂t∂x fr ∂x 5pes + ∂(ρEu + pu) s = q∏ , E = c: T + u2∙ ∂t∂x v 2 ∂x (4) ∂x p=PRT. (5) В зонах изменения конфигурации и сечения горных выработок (сопряжения горных выработок, примыкания вертикальных выработок, повороты) будем использовать трехмерный подход. Течение газа в этой области описывается трехмерными уравнениями газовой динамики. Будем называть эту область выработок «узел». Каждый узел представляется в виде одной трехмерной ячейки, к сторонам которой могут примыкать смежные ветви различного сечения под разными углами. В декартовой системе координат x, y, z уравнения газовой динамики имеют вид ∂ρ + ∂ρ u + ∂ρv + ∂ρw 0 ∂t ∂x ∂y ∂z ’ (6) ∂ρ f∂ρ fu∂ρ fv∂ρ fw f+f+f+f=0 ; ∂t ∂x ∂y ∂z (7) ∂ρu+∂(ρu2 + p)+∂ρuv+∂ρuw=0 ; ∂t ∂x ∂y ∂z (8) ∂ρv+∂ρuv+∂(ρv2 + p)+∂ρvw=0 ; ∂t∂x ∂y ∂z (9) ∂ρw + ∂ρuw + ∂ρvw + ∂(ρw2 + p)=0; ∂t ∂x ∂y ∂z (10) 110 Е.Е. Мазепа, П.И. Кусаинов, О.Ю. Лукашов, А.Ю. Крайнов ∂ρE ∂(ρEu + pu) ∂(ρEv + pv) ∂(ρEw + pw) u2+v2+w2 -----t-L = 0, E = cvT+-----------; (11) ∂t ∂χ ∂y ∂z p=ρRT. (12) Предполагая, что известны координаты, протяженность, значения состояния среды после взрыва в зоне взрыва, а также параметры рудничной атмосферы, ставим следующие начальные условия для прямолинейных участков выработок: p'χ,°'={p° XїВМ,: t'χ,°'={T° X7,Мі. Γρ f X ∈ ВМ, P f(Х,0) = Lf т, u(X,0') = 0, (13) X0 X ∉ВМ, и для сопряжений: p(χ,y,z,0)= Pb в ВМ, T(χ,y,z,0)= Tb в ВМ, р(χ,у,",0' = {p0 вне ВМ, t'χ,у,χ,0> = St0 вне ВМ, fρ f в ВМ, (14) ρf(х,У,z,0) = ^ f f [0 вне ВМ, u(χ,y,z,0)= 0, v(χ, y, z,0) = 0, w(χ, y, z,0) = 0. Граничные условия для уравнений (1) - (4) и (6) - (11) задаются в зависимости от конфигурации и сечения прилегающей выработки, с которой она граничит. Если граница выработки тупик или изолирующая перемычка, то используется условие непротекания u I гр = °. (15) Если выработка выходит на поверхность (16) p I гр ратм , ρ I гр ρ атм . Принятые обозначения: t - время, χ, y, z - оси декартовой системы координат, ρ - плотность газа, ρf - парциальная плотность продуктов сгорания, u, v, w - компоненты вектора скорости, p - давления, E - полная энергия газа, T - температура, S - площадь проходного сечения канала, Π - периметр сечения канала, τfr - сила трения о стенки, q - поток тепла в стенки, R - газовая постоянная, cv - удельная теплоемкость при постоянном объеме. Индексы: 0 - атмосферные условия, b - параметры в зоне взрыва, f - продукты сгорания, ВМ - зона взрыва метана. Для решения системы уравнений используется метод С.К. Годунова [11]. Реализация алгоритма численного расчета В обобщенном виде алгоритм расчета можно описать следующим образом: 1. Прямолинейные участки горных выработок (ветви) разбиваются на серии одномерных расчетных ячеек (рис. 1); 2. В сопряжениях, поворотах и точках изменения геометрических параметров выработок (узлах) размещаются трехмерные ячейки; 3. Во всех ячейках устанавливаются начальные условия; 4. Выполняется расчетный цикл: О численном решении задачи распространения воздушных ударных волн 111 4.1) на текущем временном слое определяется шаг по времени с учетом критерия Куранта; 4.2) на границе каждой пары смежных ячеек решается задача о распаде произвольного разрыва; 4.3) в ячейках производится решение уравнений (1) - (5) или (6) - (12), и нахождение параметров среды на новом временном слое; 5. Процедура итеративно повторяется до факта достижения критерия остановки. Выходы на поверхность Вертикальные выработки Горизонтальная выработка Тупик Рис. 1. Простейшая модельная сеть выработок (а) и соответствующая ей дискретная расчетная область (б) Fig. 1. (a) The simplest model network of workings and (b) the corresponding discrete computational domain Первичная реализация алгоритма [9] с выделением объектов типа «узел» и «ветвь» и объединением этих объектов в массивы (рис. 2), показала, что этот вариант является рабочим, но низкоэффективным в плане адаптации к последующему развитию математической модели и метода решения: - использование массивов усложняет связывание ячеек на границах областей с разными размерностями ячеек; - при распараллеливании расчетов путем геометрической декомпозиции расчетной области сложно выполнять балансировку нагрузки на вычислительные потоки, так как основными единицами являются не отдельные расчетные ячейки, а ветви с существенно различающимися длинами; г ю S Q -Q 03 S - массив 2 - λ m ю S . Q Q 03 S л Iiiiiiiiiiii π□ массив 4 Рис. 2. Демонстрация подхода, когда в качестве основной расчетной единицы принята ветвь с массивом ячеек внутри Fig. 2. Demonstration of the approach when a branch with an array of cells inside is accepted as a main calculation unit 112 Е.Е. Мазепа, П.И. Кусаинов, О.Ю. Лукашов, А.Ю. Крайнов - отсутствует универсальный механизм изменения граничных условий. Такой механизм, например, нужен при наличии глухих взрывоустойчивых перемычек, когда до разрушения перемычки соответствующая ей граница должна рассматриваться как глухая, а после разрушения стать проходной; - отсутствует универсальный механизм замены математической модели. Такой подход может быть полезен, когда в сети выработок имеются, например, водоналивные перемычки - до их срабатывания имеет смысл использовать газодинамическую модель, а после срабатывания - модель с учетом капель воды в газе. Опыт использования и модернизации первичного варианта алгоритма [9] показал, что возможен более эффективный подход к реализации расчетного алгоритма. Он основан на идее, заложенной в методе контрольных объемов - выделении элементарного объема, имеющего связи со смежными ячейками, в котором реализуется алгоритм решения системы уравнений математической модели. Связи между смежными ячейками устанавливаются на каждой границе объема, а количество этих границ соответствует размерности данной ячейки (рис. 3). Это делает ячейку элементарной (описывающую поведение системы в данном элементарном объеме), изолированной (остаются лишь сама ячейка и ее связи со смежными объемами) и самодостаточной (для расчета параметров газа есть все необходимое -уравнения математической модели, граничные и начальные условия). Рис. 3. Контрольные объемы в случаях одномерной (а) и трехмерной (b) ячеек, и границы, через которые устанавливаются связи между смежными ячейками Fig. 3. Control volumes for (a) one-dimensional and (b) three-dimensional cells, and the connecting boundaries between adjacent cells Предлагается универсальный механизм связывания ячеек: достаточно создать экземпляры этих ячеек, связать их и указать используемую математическую модель, описывающую поведение системы именно в этом объеме. Причем таких моделей может быть одна или несколько, они могут динамически меняться в зависимости от разных физических факторов в процессе выполнения численных расчетов (с соблюдением соответствующих правил трансформации между моделями и с соблюдением законов сохранения). В результате экспериментов по реализации алгоритма расчета ВУВ «граница» была выделена в виде самостоятельного элемента алгоритма с наличием данных о двух смежных ячейках, которые она соединяет (рис. 4, а). Основными исходными данными для решения задачи о распаде произвольного разрыва выступают плотность, давление и скорости потоков в смежных ячейках. Эти данные не зависят от размерности смежных ячеек, поэтому «граница» позволяет единообразно связывать ячейки любых размерностей (рис. 4, b). О численном решении задачи распространения воздушных ударных волн 113 Граница p1, ρ1, U1 ...⅛ p2, ρ2, U2 .......‰. ....... Рис. 4. Связывание двух смежных ячеек через границу (а), пример организации связей в сопряжении нескольких ветвей (b) Fig. 4. (a) Binding of two adjacent cells through a boundary; (b) an example of linking in a conjugation of several branches «Границу» также удобно использовать для моделирования особых граничных случаев - выходы поверхность, тупик. Это может достигаться за счет применения как специальных видов ячеек, находящихся зеркально по отношению к соответствующим пограничным ячейкам расчетной области, так и использованием специальных видов «границ», имитирующих наличие таких ячеек. Во втором случае удобно моделировать наличие вентиляционных сооружений (дверей) и взрывоустойчивых перемычек с последующей заменой типа границы в случае разрушения перемычки (рис. 5). Зона взрыва Целая перемычка (глухая граница) Рис. 5. Пример замены типа границы: а - изолированный взрыв; b - распространение ВУВ после взрыва перед моментом разрушения перемычки; c -распространение ВУВ после разрушения перемычки. Разрушение перемычки сопровождается заменой граничного условия в месте установки перемычки с условия не протекания на расчет в ветви Fig. 5. An example of the boundary type alternation: (a) an isolated dead end in which an explosion is expected; (b) the shock wave propagation induced by the explosion before the stopping is destructed; and (c) the shock wave propagation following the destruction of the stopping. The destruction is accompanied by a replacement of the boundary condition on the stopping installation site from the impermeability condition to the calculation in the branch Зона распространения ВУВ после разрешния Зона распространения Целая перемычка (глухая граница) Перемычка отсутствует (проходная граница) тупик, в котором ожидается Необходимо отметить еще один важный момент - необходимость учета углов, под которыми сопрягаются между собой реальные горные выработки. Это выражается в наличии таких углов между крайними одномерными ячейками ветвей с трехмерными ячейками смежных узлов и необходимостью расчета соответствующих проекций - нормальных и тангенциальных составляющих скоростей. Для 114 Е.Е. Мазепа, П.И. Кусаинов, О.Ю. Лукашов, А.Ю. Крайнов автоматизации расчета этих проекций была создан дополнительный элемент алгоритма - «грань» ячейки, которая выполняет эти действия на основе информации о номере стороны расчетной ячейки, к которой она прилегает, и углов сопряжения смежных ветвей (рис. 6). Рис. 6. Дискретная модель поворота горной выработки (а) и соответствующая ему связка ячеек, граней и границ (b) Fig. 6. (a) A discrete model of the bend in a mine working and (b) the corresponding binding of the cells, faces, and boundaries С учетом определения понятий «ячейка», «грань» и «граница» в алгоритме порядок расчета можно детализировать следующим образом: 1. Прямолинейные участки горных выработок (ветви) разбиваются на серии одномерных расчетных «ячеек»; 2. В сопряжениях, поворотах и точках изменения геометрических параметров выработок (узлах) размещаются трехмерные «ячейки»; 3. Для каждой одномерной «ячейки» создается по две «грани», которые присоединяются к участвующим в вычислениях сторонам «ячейки», для каждой трехмерной «ячейки» создается по шесть «граней» - для всех ее сторон; 4. Создаются «границы», каждая из которых связывается с двумя «гранями» смежных «ячеек»; 5. Во всех «ячейках» устанавливаются начальные условия; 6. Выполняется расчетный цикл: 6.1) на текущем временном слое определяется шаг по времени с учетом критерия Куранта; 6.2) проверяются факторы, определяющие замену границ, например непроте-каемых на проходные; 6.3) на границе каждой пары смежных ячеек решается задача о распаде произвольного разрыва, при этом нормальные составляющие векторов скоростей определяются через смежные грани; 6.4) в ячейках производится решение уравнений (1) - (5) или (6) - (12) и нахождение параметров среды на новом временном слое. На этом этапе требуется лишь решение уравнений, непосредственно составляющих исходную математическую модель. При необходимости, модель может быть динамически заменена; 7. Процедура итеративно повторяется до факта достижения критерия остановки. Проверка алгоритма расчета Для проверки корректности расчетов с использованием описанного алгоритма был проведен ряд модельных и тестовых расчетов. Модельные схемы каналов (выработок) показаны на рис. 7. О численном решении задачи распространения воздушных ударных волн 115 о Рис. 7. Модельные схемы каналов: а - окружность, b - замкнутый квадрат, c - «разрезная печь»: 1 - 3D ячейка, 2 - 1D ячейка, 3 - ячейка с открытой границей, 4 - зона повышенного давления, L, n, m - длины ветвей Fig. 7. (a) A circle, (b) a closed square, and (c) a “crosscut”: 1, a 3D cell; 2, a 1D cell; 3, a cell with an open boundary; and 4, a pressure zone; L, n, m are the branch lengths I - I I *∙-»-J n 9 • m ~EΞΓ Все схемы являются модельными, однако позволяют оценить корректность нового алгоритма расчета параметров ВУВ: проверить выполнимость основных законов сохранения, симметричность распространения ударных волн при симметричности геометрии каналов и параметров взрыва. Проведен расчет распространения ВУВ в канале, представляющем собой кольцо (рис. 7, а). Зона взрыва задавалась в нижней части кольца, давление в зоне взрыва составляло 570 кПа. Результаты расчета показаны на рис. 8. Видно, что ВУВ распространяются симметрично в разные стороны. В момент времени t = 0.6 c образуется новый пик давления, обусловленный взаимодействием волн давления, распространяющихся навстречу друг другу. Расстояние, м Рис. 8. Распределения давления (а) и температуры (b) в разные моменты времени в кольцевом канале Fig. 8. Distribution of the (а) pressure and (b) temperature in the annular channel at various time instants Расстояние, м 116 Е.Е. Мазепа, П.И. Кусаинов, О.Ю. Лукашов, А.Ю. Крайнов Результаты расчета распространения ВУВ в замкнутой квадратной конфигурации каналов (горных выработок, рис. 7, b) представлены на рис. 9. В момент времени t = 0.1 c образуется два пика на одинаковом расстоянии, связанные с поворотом канала на 90°. В противоположенной от места взрыва ветви ударные волны сходятся и наблюдается всплеск давления в момент времени 0.6 с. Рис. 9. Распределения давления в разные моменты времени в канале квадратной конфигурации Fig. 9. Pressure distributions in the square channel at various time instants Следующий расчет был проведен для конфигурации каналов «разрезная печь» (рис. 7, с). Результаты расчетов в виде распределений давления по путям распространения ударных волн, изображенным на рис. 10 штриховой и штрихпунктир-ной линиями, представлены на рис. 11. Рис. 10. Схемы отображения топологии рис. 7, c: 1 - путь № 1, 2 - путь № 2 Fig. 10. Topology diagrams for Fig. 7c: 1, way No. 1 and 2, way No. 2 О численном решении задачи распространения воздушных ударных волн 117 Расстояние, м Расстояние, м Рис. 11. Распределения давления по путям № 1 и № 2 «разрезной печи» в разные моменты времени Fig. 11. Pressure distributions along tracks No. 1 and No. 2 in the “crosscut” at various time instants Из полученных результатов, представленных на рис. 11, видно, что ударная волна, распространяющаяся по пути № 1, при прохождении мест сопряжения с перемыкающей выработкой немного уменьшает свою интенсивность. На пути № 2 в примыкающей выработке в момент времени 0.9 с наблюдается всплеск давления, изображенный на рис. 11, b, обусловленный встречным взаимодействием ударных волн, идущих по ней от двух сопряжений с выработками пути № 1. В процессе всех расчетов, результаты которых представлены на рис. 8, 9, 11, контролировалась выполнимость законов сохранения массы и полной энергии газа в выработках. Полная масса газа и полная энергия газа в выработках при численном решении сохраняется с точностью 99.95 %. При симметричной геометрии каналов и места взрыва расчетные распределения давления в каналах при распространении ударных волн также симметричны. Вариант алгоритма [9] используется при составлении плана ликвидации аварий (ПЛА) на горнорудном предприятии. Оперативная часть ПЛА разрабатывается для руководства действиями обслуживающего персонала шахты при возникновении аварийной ситуации или аварии [7]. Оперативной частью ПЛА должны охватываться все виды возможных аварий и аварийных ситуаций в шахте, в связи с чем скорость получения результатов влияет на убытки предприятия и человеческие жертвы. Изложенный в данной статье алгоритм может быть доработан с целью распараллеливания. Это позволит ускорить проведение массовых расчетов для составления ПЛА. Заключение Разработан подход к реализации метода решения задач о распространении ВУВ в разветвленной сети горных выработок с учетом произвольных углов их сопряжения. Подход базируется на использовании метода С.К. Годунова и выделении элементов алгоритма «ячейка», «грань» и «граница». Каждый из указанных элементов алгоритма отвечает за этап вычислений: «граница» - за решение задачи 118 Е.Е. Мазепа, П.И. Кусаинов, О.Ю. Лукашов, А.Ю. Крайнов о распаде произвольного разрыва в параметрах газа, «грань» - за определение нормальных и касательных к «грани» компонентов векторов скоростей в ячейках, «ячейка» - за численное решение уравнений, составляющих математическую модель распространения ВУВ. Предлагаемая декомпозиция расчетов направлена на достижение сразу нескольких целей: - существенное снижение трудозатрат при развитии математической модели за счет введения в нее новых физических процессов, например учета водяных и сланцевых заслонов; - упрощение распараллеливания вычислений; - упрощение установления связей между ячейками с разными размерностями; - существенное упрощение реализации расчетов распространения ВУВ в геометрически сложных областях. Указанный подход реализован в виде компьютерной программы, выполнены тестовые расчеты, решены модельные задачи с проверкой выполнения законов сохранения массы, импульса, энергии. Разработанный подход планируется применять при решении задач промышленной безопасности на угольных шахтах России.

Скачать электронную версию публикации