# Abstract

Mathematical statement of direct and inverse problem of Stefan for horizontal layer of rock massif with homogenous and isotropic thermophysical properties is presented. It is assumed as a hypothesis that heat transfer in vertical direction is negligible compared to heat exchange in horizontal plane. At the initial moment, the rock massif has a uniform temperature and the temperature on surfaces of freezing columns was the same for all columns and constant in time. A method proposed allows getting an approximate solution of the direct Stefan problem for a single freezing column with a small consumption of computational resources. Based on a proposed method, a high-speed algorithm for solving inverse Stefan problem for the case of a single freezing column is built. An algorithm is based on the gradient descent method. The effect on the solution of different types of functions used is analyzed. Functions approximate the temperature field in a cooling zone. It is established that time dependence of the radius of a phase transition front essentially depends on the type of an approximation function. The most preferable is an integral exponential function that is a solution to the one-dimensional heat equation in cylindrical coordinates. Then, proposed technique and algorithm are considered for the case of variety of freezing columns that form circle counter and random number of control wells. Results of the calculation of inverse Stefan problem for conditions of the shaft No. 1 of the mine being under construction at the Petrikovsky ore mining and processing enterprise are presented. The calculation included well inclinometry based on geological data. It was studied how measurements of the temperature made at different wells can affect obtaining solution. Options of interpretation of inconsistency of temperatures measured in control wells are offered. Probabilistic analysis of ice wall thickness is carried out.

# Full Text

Введение Одним из наиболее универсальных и надежных способов проведения вертикальных шахтных стволов в сложных гидрогеологических условиях является искусственное замораживание горных пород [1]. Суть способа заключается в бурении кольца скважин вокруг проектируемого сечения ствола, после чего в скважинах монтируются замораживающие колонки, по которым циркулирует хладоноситель. Таким образом, со временем вокруг будущей горной выработки создается защитное ограждение из мерзлой водонасыщенной породы, предохраняющее выработку от поступления подземных вод при ее строительстве. Эффективность данного защитного ограждения, называемого ледопородным ограждением (ЛПО), зависит в конечном счете от его герметичности и толщины. Производство работ по проходке подземных выработок после замораживания разрешается только после образования замкнутого замороженного контура проектной толщины. При этом в процессе выполнения работ по замораживанию породного массива важно осуществлять непрерывный контроль текущего состояния ЛПО и своевременно определять момент достижения достаточной толщины и сплошности контура из замороженных пород [2]. За последние годы несовершенство существующих способов контроля формирования ледопородного ограждения и недостаточная точность расчетов его параметров являлись причиной проблем с обеспечением герметичности ледопородного ограждения на таких горных предприятиях, разрабатывающих месторождения с повышенной обводненностью, как Гремячинский ГОК ОАО МХК «ЕвроХим», Гарлыкский ГОК ГК «Туркменхимия», а также рудоуправления Верхнекамского месторождения калийно-магниевых солей [3]. В итоге возникшие аварийные случаи при строительстве шахтных стволов привели к снижению безопасности ведения горных работ и появлению дополнительных финансовых затрат. В связи с этим в настоящее время отделом аэрологии и теплофизики ГИ УрО РАН внедряется новая перспективная система определения герметичности и толщины ЛПО, основанная на использовании технологии оптоволоконной термометрии [3, 4]. Данная технология позволяет получать в любой момент времени распределение температур по глубине контрольных скважин, расположенных на некотором удалении от контура замораживающих колонок и служащих для мониторинга текущего теплового режима породного массива. При этом появляется задача восстановления распределения температур во всем породном массиве по данным измерений температуры в контрольных скважинах, количество которых, как правило, невелико (3-5 шт.). Математически это приводит к постановке обратной задачи Стефана [4]. Ранее в существующих исследованиях формирования ЛПО не рассматривались обратные задачи теплопроводности (ОЗТ) и обратные задачи Стефана (ОЗС). Математический аппарат, применявшийся для исследования формирования ЛПО в России и за рубежом, - это исключительно прямые задачи теплопроводности и Стефана [1, 5-12]. При этом в литературе существует ряд исследований, посвященных анализу корректности постановки ОЗТ и ОЗС в общем случае [13, 14], анализу численных методов решения ОЗС и ОЗТ [15-18]. Классически выделяется две группы методов получения устойчивых решений обратных задач: саморегуляризация, когда осуществляется управление мерой близости получаемого решения к точному посредством вариации параметров вычислительных алгоритмов, и регуляризация по методу А.Н. Тихонова, подразумевающая построение регуляризующего оператора, сглаживающего функционала и его последующую минимизацию [14]. В литературе также представлен класс методов дескриптивной регуляризации, когда вводятся дополнительные стабилизирующие ограничения для искомого решения [18]. Общими недостатками существующих методов являются алгоритмическая сложность, существенные затраты вычислительных ресурсов и необходимость дополнительных проверок на сходимость и устойчивость решения. Целью настоящей работы является разработка скоростного метода решения обратной задачи Стефана для анализа процесса формирования ледопородного ограждения строящихся шахтных стволов. Актуальность работы обосновывается необходимостью на практике быстро и с приемлемой точностью определять фактические толщину и сплошность ледопородного ограждения. Идея работы состоит в определении приближенного аналитического вида распределения температур в горизонтальном разрезе массива для случая множества скважин и его использовании для решения обратной задачи Стефана, которое будет таким образом сведено к определению явной зависимости для функционала рассогласования температур в контрольных скважинах и его минимизации методом градиентного спуска. Математическая постановка прямой и обратной задач Стефана Рассматривается горизонтальный слой породного массива с однородными и изотропными теплофизическими свойствами. Принимается, что теплообмен в вертикальном направлении пренебрежимо мал по сравнению с теплообменом в горизонтальной плоскости. В начальный момент времени породный массив имеет однородную температуру T0, а температура Tw на поверхностях замораживающих колонок считается одинаковой для всех колонок и неизменной во времени (рис. 1). Предполагается, что происходит полное превращение воды в лед [19, 20]. Рис. 1. Контур замораживающих колонок в горизонтальном разрезе породного массива: 1 - зона льда; 2 - зона охлаждения, пунктирная линия - текущее положение границы фазового перехода В этом случае математическая постановка двухмерной прямой задачи Стефана в полярных координатах сведется к следующей системе уравнений [5]: , (1) , (2) , (3) , (4) , (5) , (6) (7) где k1 - теплопроводность среды в зоне льда, Вт/(м·°С); k2 - теплопроводность среды в зоне охлаждения, Вт/(м·°С); c1 - удельная теплоемкость среды в зоне льда, Дж/(кг·°С); c2 - удельная теплоемкость среды в зоне охлаждения, Дж/(кг·°С); ρ - плотность среды, кг/м3; Tph - температура фазового перехода, °С; L - удельная теплота плавления льда, Дж/кг; w - влагосодержание массива, кг/кг; Гw - поверхность замораживающих колонок; Гph - поверхность фронта фазового перехода; Rv - радиус внешней границы расчетной области, м; Tw - температура стенок замораживающих колонок, °С; T0 - температура непотревоженного породного массива на удалении, °С. Коэффициентная обратная задача Стефана к прямой задаче (1)-(7) формулируется следующим образом: при условии, что известна зависимость температуры от времени на поверхностях Гсw контрольных скважин , (8) требуется найти значения теплофизических свойств породного массива k1, k2, c1, c2, L при 0 < t ≤ T и найти функцию T(r, φ, t), которые бы удовлетворяли задаче (1)-(7). Дополнительно вводятся ограничения на теплофизические свойства породного массива - соответствующие минимально возможные и максимально возможные значения, которые определяются исходя из инженерно-геологических изысканий [21, 22]. Аппроксимация решения прямой задачи Стефана для случая одиночной замораживающей колонки Прежде всего в данном исследовании был рассмотрен наиболее простой случай одиночной замораживающей колонки. Целью являлось нахождение приближенного аналитического решения прямой задачи Стефана для данного случая. В этом случае система (1)-(7) упрощается: исчезает зависимость от угла φ, а поверхность Гw замораживающих колонок сводится к окружности радиуса rw. Предлагается записать общее решение системы (1)-(7) для случая одиночной замораживающей колонки в виде суперпозиции локальных решений отдельно для зоны льда и зоны охлаждения. Сопряжение двух данных локальных решений будет происходить по границе фазового перехода rph = rph(t). В зоне льда в качестве локального решения принимался логарифм, соответствующий случаю квазистационарного распространения теплоты [9]: (9) Как показали предварительно проведенные численные расчеты прямой задачи Стефана для случая одной скважины, такое приближение привносит ошибку в определение радиуса фронта ледопородного ограждения не более чем на 10 % в исследуемом диапазоне времен и толщин ледопородного ограждения [3]. В зоне охлаждения рассматривалось несколько вариантов функционального вида локальных решений: экспоненциальный вид, тригонометрический вид и вид интегральной показательной функции: , (10) , (11) . (12) Здесь a2 - температуропроводность среды в зоне охлаждения, м2/с. Единственным физически обоснованным вариантом здесь является интегральная показательная функция, которая представляет собой элементарное решение цилиндрического одномерного уравнения теплопроводности [23]. Только данная функция обеспечивает выполнение закона о сохранении количества тепловой энергии. Два дополнительных варианта локальных решений следует трактовать как чисто математический прием подбора аппроксимации решения прямой задачи Стефана для одиночной замораживающей колонки в зоне охлаждения. При этом параметр a2 для формул (11) и (12) уже не является в строгом смысле температуропроводностью. Предложенные базисные функции (10)-(12) представлены на рис. 2. В качестве расчетных параметров задачи приняты: температура непотревоженного массива +10 °С, температура поверхности замораживающих колонок -20 °С, расчетное время - 40 сут. На рис. 2 видно, что экспоненциальная функция достаточно близка к точному решению цилиндрического уравнения тепловой диффузии - интегральной показательной функции. При этом тригонометрическая функция дает заниженные (до 1,5 °С) значения температур в зоне охлаждения, однако это можно исправить нормировкой коэффициента при температуропроводности в формуле (10). Рис. 2. Приближенное решение прямой цилиндрической задачи Стефана с использованием различных аппроксимирующих функций в зоне охлаждения: 1 - экспоненциальный вид; 2 - тригонометрический вид; 3 - вид интегральной показательной функции Приближенное аналитическое решение (9)-(12) должно не только правильно отражать распределения температур в зоне льда и охлаждения, но и давать адекватный прогноз временной динамики фронта фазового перехода r(t). Функция r(t) определяется из решения уравнения (3) после подстановки в него градиентов температур, рассчитанных из функции (9) и одной из функций (10)-(12). На рис. 3 представлены зависимости радиуса фронта фазового перехода от времени, рассчитанные для каждой из функций (10)-(12). Расчет производился численно с использованием явного метода Рунге-Кутты 4-го порядка [24]. Как следует из рис. 3, все три выбранные для зоны охлаждения функции справляются с задачей адекватного прогноза временной динамики фронта фазового перехода (рассогласование менее 6 %). Наиболее жесткую оценку снизу для радиуса фронта фазового перехода дает решение в виде интегрально-показательной функции. Таким образом, при решении только прямой задачи Стефана лучше использовать интегральную показательную функцию для представления решения. Рис. 3. Временная диаграмма фронта фазового перехода с использованием различных аппроксимирующих функций в зоне охлаждения: 1 - экспоненциальный вид; 2 - тригонометрический вид; 3 - вид интегральной показательной функции Алгоритм решения обратной задачи Стефана для случая одиночной замораживающей колонки Ранее в настоящей работе было отмечено, что наиболее жесткую оценку снизу для радиуса фронта фазового перехода при решении прямой задачи Стефана дает интегральная показательная функция. При анализе обратной задачи Стефана данная оценка теряет значение, поскольку расчет толщины ледопородного ограждения производится по замерам в контрольных скважинах. На основании условия (8) можно записать следующий функционал рассогласования температур в контрольных скважинах, который в дальнейшем будет подвержен процедуре минимизации: (13) Суммирование в выражении (13) осуществляется по контрольным скважинам ei. В качестве нормы в формуле (14) для оценки отклонения решения от требуемых значений на временном интервале [t1, t2] принимается следующая скалярная величина, являющаяся нормой в лебеговом пространстве L2 [25]: (14) Требуется определить параметры, по которым осуществляется минимизация функционала (13). Классически [13] выделяются следующие обратные задачи теплообмена: 1) ретроспективная, когда регулируемыми параметрами являются начальные условия; 2) граничная, когда регулирование осуществляется граничными условиями; 3) коэффициентная, для которой варьируются физические свойства теплопроводной среды (производится идентификация оператора теплопроводности). Как показывают исследования [21, 22, 26, 27], при построении математической модели формирования ЛПО в породном массиве наименьшую достоверность имеют физические свойства породного массива и содержащихся в них грунтовых вод. Это связано с неоднородностью реального породного массива, несовершенством методик определения физико-механических свойств породного массива по извлеченным образцам керна, недостаточностью статистических выборок образцов и пр. По этим причинам в данном исследовании обратная задача Стефана считается коэффициентной. В качестве варьируемых параметров может быть принято до четырех независимых теплофизических параметров среды. Это связано с тем, что решение T(r, φ, t) зависит от четырех безразмерных чисел задачи (1)-(7): Фурье в зоне льда Fo1 и зоне охлаждения Fo2, Стефана в зоне льда Ste1 и зоне охлаждения Ste2. В данной работе рассматривается случай двух варьируемых теплофизических параметров - теплопроводностей среды в зоне льда k1 и зоне охлаждения k2. Для каждого из данных параметров устанавливается область допустимых значений и на основании инженерно-геологических изысканий. Таким образом, поставлена задача минимизации функционала (13) по параметрам теплопроводности k1 и k2 породного массива. В данной работе был выбран следующий итерационный алгоритм решения обратной задачи Стефана: 1. Определение начальных приближений для параметров минимизации - теплопроводностей k1 и k2. 2. Численный расчет зависимости радиуса фронта фазового перехода rph(t) от времени с использованием явного метода Рунге-Кутты 4-го порядка. 3. Численное интегрирование рассогласований температур для каждой контрольной скважины для текущих теплопроводностей (k1, k2), а также для теплопроводностей, отклоненных от (k1, k2) на малую величину: (k1 + Δk, k2) и (k1, k2 + Δk). Вычисление интеграла рассогласований I и его первых производных. Сравнение интеграла рассогласований I с требуемой точностью Iacc. В случае, если I < Iacc, выход из процедуры, завершение расчета. 4. Расчет приращений для параметров минимизации методом градиентного спуска [28, 29]. Задание новых значений теплопроводностей k1 и k2. Проверка новых значений теплопроводностей на предельные минимальные и максимальные значения. 5. Возврат к пункту 2 (новая итерация). Пример решения обратной задачи Стефана. Случай одиночной контрольной скважины На рис. 4, 5 представлен пример решения обратной задачи Стефана методом минимизации аппроксимирующей функции для случая одиночной контрольной скважины с заданной зависимостью температуры от времени в виде . (15) Рис. 4. Временная диаграмма температуры массива на поверхности контрольной скважины: 1 - целевая функция; 2 - начальное приближение решения; 3 - решение обратной задачи Стефана методом минимизации аппроксимирующей функции Входные параметры задачи: теплопроводность среды в зоне льда 2,5 Вт/(м·°С); теплопроводность среды в зоне охлаждения 1,5 Вт/(м·°С); удельная теплоемкость среды в зоне льда 800 Дж/(кг·°С); удельная теплоемкость среды в зоне охлаждения 850 Дж/(кг·°С); плотность среды 2000 кг/м3; температура фазового перехода 0 °С; удельная теплота плавления льда 330 кДж/кг; влагосодержание массива 0,3 кг/кг; температура стенок замораживающей скважины -20 °С; температура непотревоженного породного массива на удалении +10 °С. В качестве аппроксимирующих функций рассмотрены интегральная показательная, тригонометрическая и экспоненциальная. Достичь хорошей сходимости решения удалось для всех трех видов функций, однако получаемые при этом распределения температур, результирующие значения коэффициентов теплопроводности и толщин ледопородного ограждения существенно различны между собой (см. рис. 4, 5). Рис. 5. Решение обратной задачи Стефана; временная диаграмма фронта фазового перехода с использованием различных аппроксимирующих функций в зоне охлаждения: 1 - вид интегральной показательной функции; 2 - тригонометрический вид; 3 - экспоненциальный вид Из анализа графиков рис. 4 становится ясно, что каждое решение T = T(r,t) приблизилось к заданной функции Te(t) настолько близко, насколько это было возможно с учетом ограниченности функционального вида этого решения, а также граничных и начальных условий. Различия функционального вида решения обусловили также различия в результирующих коэффициентах теплопроводности среды, к которым сошлись итерационные процедуры алгоритма решения ОЗС (табл. 1), и различия в получаемых толщинах ледопородного ограждения (см. рис. 5). Таблица 1 Коэффициенты теплопроводности, полученные при решении ОЗС Коэффициент теплопро-водности Интегральная показательная функция Экспонен-циальная функция Тригономет-рическая функция В зоне льда 3,5 3,0 2,8 В зоне охлаждения 2,6 0,6 1,5 Если принять кривую rph(t) для интегральной показательной функции за достоверную, то ошибка определения толщины ледопородного ограждения на 15-е сутки для выбранных расчетных параметров при использовании тригонометрической функции составляет 20 %, а при использовании экспоненциальной функции - 50 %. Решение обратной задачи Стефана для случая множества замораживающих колонок Предложенный метод решения обратной задачи Стефана может быть перенесен на случай множества замораживающих колонок и множества контрольно-термических скважин (рис. 6). Рис. 6. Распределение температур в горизонтальном срезе породного массива, полученное в результате решения обратной задачи Стефана В случае множества замораживающих колонок аппроксимационная зависимость температуры от времени, пространственных координат и теплофизических параметров задачи может быть определена следующим образом. Расчетная область разбивается на три подобласти: 1) зона льда около замораживающих колонок, 2) зона охлаждения внутри контура замораживающих колонок, 3) зона охлаждения снаружи контура замораживающих колонок. Функциональный вид поля температур в зоне льда около замораживающих колонок по-прежнему определяется по формуле (9). Поскольку функция (9) является потенциальной, ее суперпозиция для множества скважин не приводит к нарушению равенства в уравнении (1): (16) Здесь Nf - количество замораживающих колонок; rwi - радиус-вектор центра замораживающей колонки № i; rph - полутолщина ледопородного ограждения (расстояние от фронта фазового перехода до контура замораживающих колонок); ccal - калибровочный параметр, вводимый для выполнения граничного условия (5). В рамках данного подхода используется приближение, согласно которому внутренний и внешний фронты фазового перехода движутся с одинаковой скоростью. Данное приближение позволяет оперировать такой величиной, как полутолщина ледопородного ограждения rph, которая равна расстоянию от любого из двух фронтов фазового перехода до контура замораживающих колонок. Функциональный вид поля температур в зоне охлаждения снаружи контура замораживающих колонок определяется функцией (12). Поле температур на удалении от контура замораживающих колонок на величину нескольких радиусов колонки rw становится осесимметричным (см. рис. 5), таким образом, поле температур в данной задаче становится подобно полю температур от некоторой одиночной абстрактной замораживающей колонки с радиусом, приблизительно равным сумме радиуса контура замораживающих колонок rc и текущего значения полутолщины ледопородного ограждения rph: . (17) Функциональный вид поля температур в зоне охлаждения внутри контура замораживающих колонок определяется с помощью ограниченного ряда Фурье-Бесселя [30]: (18) Количество членов ряда Фурье определяется на основании предварительного анализа сходимости ряда для данного момента времени. Для рассматриваемых времен образования ледопородного ограждения t ≈ 10-100 сут сходимость ряда Фурье обеспечивается уже при рассмотрении только 5-10 первых его членов. Состыковка функций (16)-(18) осуществляется на внутреннем и внешнем фронтах фазового перехода. В уравнениях (16)-(18) величина радиуса фазового перехода rph в общем случае является параметром задачи. Она может быть рассмотрена как параметр оптимизации или может быть найдена из численного решения обыкновенного дифференциального уравнения первого порядка (3) после подстановки в него градиентов температур, рассчитанных из функций (16)-(18). Следует отметить, что в рамках предложенного подхода считается, что замораживающие колонки образуют правильный круговой контур (отсутствие инклинометрии). Предложенный подход будет иметь дополнительную погрешность тем большую, чем больше инклинометрия скважин. Далее определяется функционал рассогласования температур в контрольных скважинах и минимизируется методом градиентного спуска. На рис. 6 представлено распределение температур в горизонтальном срезе породного массива, полученное для условий шахтного ствола № 1 строящегося рудника Петриковского горно-обогатительного комбината. На рис. 6 также представлен контур из 41 замораживающей колонки и четыре контрольные скважины: КТ-1 (-1,5 °С), КТ-2 (-1,1 °С), УКЛО (2,5 °С), ГН-2 (5,7 °С). Температуры контрольных скважин соответствуют времени замораживания 60 сут. Входные параметры задачи: температура фазового перехода -2,1 °С; температура стенок замораживающей скважины -20 °С; температура непотревоженного породного массива на удалении +6,3 °С; радиус контура замораживающих колонок 8,2 м; радиус замораживающей колонки 0,17 м; плотность среды 1960 кг/м3; начальная теплопроводность среды в зоне льда 1,5 Вт/(м·°С); начальная теплопроводность среды в зоне охлаждения 1,29 Вт/(м·°С); удельная теплоемкость среды в зоне льда 1421 Дж/(кг·°С); удельная теплоемкость среды в зоне охлаждения 1566 Дж/(кг·°С). В результате минимизации функционала рассогласования температур получено, что теплопроводности в зоне льда и зоне охлаждения равны 1,24 и 0,97 Вт/(м·°С) соответственно. Толщина ЛПО составила 1,94 м. Функционал рассогласования температур в конце итерационной процедуры равен 0,6 °С. Ненулевое значение функционала рассогласования температур обусловлено: - погрешностью численного метода; - несовместностью условий (8) для четырех различных скважин. Случай несовместности измерений температур несколькими контрольными скважинами При поиске решения обратной задачи Стефана для случая наличия нескольких контрольных скважин возникает вопрос о несовместности условий (8) для различных скважин. К примеру, если рассмотреть задачу о формировании симметричного ледопородного ограждения около единично взятой скважины, рассмотреть две контрольные скважины, удаленные с разных сторон на одинаковое расстояние от замораживающей колонки (рис. 7), и задать в них различные значения температур вида (8), то удовлетворить обоим этим условиям одновременно при решении обратной задачи Стефана не удается. 1 Рис. 7. Геометрическая модель породного массива с замораживающей колонкой и двумя контрольными скважинами: 1 - замораживающая колонка; 2 - фронт фазового перехода; 3 - контрольные скважины; 4 - зона льда; 5 - зона охлаждения породного массива Если формально построить функционал вида (13), представляющий сумму рассогласований по всем контрольным скважинам, и минимизировать его каким-либо корректным способом, то скорее всего удастся получить некоторое среднее решение, которое как-то рассогласуется с каждой из контрольных скважин, но в целом по всем скважинам дает минимальное рассогласование. Такой формальный подход, скорее всего, применим в случае, когда несовместность показаний температур по всем контрольным скважинам мала и среднее решение, получаемое из минимизации формально построенного функционала (13), слабо рассогласуется с каждой из контрольных скважин. Однако на практике часто возникают случаи сильно несовместных показаний температур по контрольным скважинам. Это может быть связано как с погрешностью измерительной процедуры в контрольных скважинах, так и с погрешностями измеренных теплофизических параметров породного массива (неоднородность свойств), неточностью технологических показателей (температура и расход закачки хладагента в замораживающие колонки). Зачастую получение более точных сведений об этих параметрах на практике для реальных строящихся стволов невозможно. Данные факты указывают на необходимость привлечения вероятностного подхода для анализа толщины ледопородного ограждения. В рамках данного подхода толщина ледопородного ограждения будет представляться вероятностной функцией от времени t и удаленности r от замораживающей колонки (или контура замораживающих колонок). Для реализации вероятностного подхода требуется задаться какими-либо гипотезами о характере рассогласования модельного распределения температур с фактическим экспериментально-измеренным для каждой скважины. Возможными вариантами являются гипотеза о недостоверных показаниях скважин и гипотеза о природной аномалии. 1. Гипотеза о недостоверных показаниях скважин. В рамках данной гипотезы считаем, что рассогласование связано с погрешностью экспериментальных измерений и вызвано малым количеством контрольных скважин по сравнению с их общим количеством. Поскольку на практике зачастую приходится иметь дело с 3-5 контрольными скважинами, можно принять, что рассогласование вызвано одной контрольной скважиной. Для устранения высокого рассогласования модельного распределения температур с фактическим предлагается убрать наиболее противоречивую контрольную скважину из рассмотрения. В случае, когда ошибка показаний наиболее противоречивой контрольной скважины существенна, данная скважина может быть исключена с помощью ручной обработки. В тех случаях, когда невозможно заранее определить контрольную скважину с недостоверными данными, предлагается произвести поочередное отбрасывание каждой скважины из рассмотрения и N раз минимизировать функционал (13) для каждой комбинации из N-1 скважины. На выходе получится выборка из N радиусов rph(t) ледопородного ограждения и N соответствующих им минимальных достижимых рассогласований (значение функционала (13)) внутри каждой комбинации из N-1 скважины. Исходя из полученных рассогласований Ii рассчитывается функция распределения вероятности толщины ледопородного ограждения в виде точек (ri, Pi) (табл. 2). (19) (20) Таблица 2 Расчет вероятности толщины ледопородного ограждения в рамках гипотезы о недостоверных показаниях скважин Показатель Номер исключенной скважины 1 2 3 4 Радиус ЛПО , м 0,25 0,28 0,29 0,34 Рассогласование , °С 0,582 0,408 0,336 0,142 Вероятность, % 100 84 59 32 Графически функция распределения вероятности толщины ледопородного ограждения представлена на рис. 8. Рис. 8. Функции распределения толщины ледопородного ограждения r: 1 - гипотеза о недостоверных показаниях скважин; 2 - гипотеза о природной аномалии 2. Гипотеза о природной аномалии. Данная гипотеза, так же как и предыдущая, предполагает, что рассогласование вызвано малым количеством контрольных скважин по сравнению с их общим количеством. Однако здесь имеется существенное отличие в том, что наименее достоверные скважины не отбрасываются, а также учитываются в расчете толщины ледопородного ограждения. Сильное отклонение в показании отдельной контрольной скважины по сравнению с остальными может быть вызвано не только погрешностью показаний. Оно также может быть связано с каким-либо геологическим нарушением (трещина, инородное включение в породах), приводящим к существенному нарушению однородности физических свойств породного массива в локальной области. Это, в свою очередь, может привести к замедлению скорости формирования ледопородного ограждения вокруг контура замораживающих скважин в локальной зоне, подверженной влиянию этого геологического нарушения. Тем самым толщина ледопородного ограждения в этой локальной зоне может быть ниже, чем в остальной части около контура замораживающих скважин. Предлагается построить N решений обратной задачи Стефана для каждой контрольной скважины в отдельности. В результате получен массив толщин ледопородного ограждения. На выходе также получаем выборку из N радиусов rph(t) ледопородного ограждения и N соответствующих им минимальных достижимых рассогласований для каждой скважины. Считая наличие геологических нарушений равновероятным, рассчитываем вероятности Pi равными для каждой скважины независимо от полученных рассогласований Ii, которые в данном случае рассмотрения каждой контрольной скважины должны быть достаточно малы. В табл. 3 представлены Таблица 3 Расчет вероятности толщины ледопородного ограждения в рамках гипотезы о природной аномалии Показатель Номер рассмотренной скважины 1 2 3 4 Радиус ЛПО , м 0,41 0,32 0,29 0,14 Вероятность, % 25 50 75 100 результаты расчета радиуса ri и вероятности того, что толщина ледопородного ограждения составляет не менее данного радиуса ri, с учетом гипотезы о природной аномалии. На рис. 8 представлены функции распределения толщины ледопородного ограждения r для каждой из двух рассмотренных гипотез. Заключение В данной работе проведено исследование обратной задачи Стефана применительно к проблеме контроля состояния ледопородного ограждения при проходке шахтных стволов. Основными научными результатами являются: 1. Математическая постановка прямой и обратной задачи Стефана для горизонтального слоя породного массива с однородными и изотропными теплофизическими свойствами. 2. Метод, позволяющий с малыми затратами вычислительных ресурсов получить аппроксимационное решение прямой задачи Стефана для случая одиночной замораживающей колонки и для случая произвольного количества замораживающих колонок, расположенных по круговому контуру. 3. Численная реализация метода решения обратной задачи Стефана для случая одиночной замораживающей колонки и для случая произвольного количества замораживающих колонок, расположенных по круговому контуру. 4. Анализ влияния на решение вида используемых функций, аппроксимирующих поле температур в зоне охлаждения. 5. Расчет обратной задачи Стефана для условий шахтного ствола № 1 строящегося рудника Петриковского горно-обогатительного комбината. 6. Исследование и интерпретация несовместности измеренных в контрольных скважинах температур, вероятностный анализ толщины ледопородного ограждения.

# About the authors

### Lev Yu. Levin

Mining Institute of the Ural Branch of the Russian Academy of Sciences - Branch of the Federal State Budgetary Institution of Science Perm Federal Research Center of the Ural Branch of the Russian Academy of Sciences

Author for correspondence.
Email: aerolog_lev@mail.ru
78 Sibirskaya st., Building A, Perm, 614007, Russian Federation

Doctor of Engineering, Head of the Department of Aerology and Thermophysics

### Mikhail A. Semin

Mining Institute of the Ural Branch of the Russian Academy of Sciences - Branch of the Federal State Budgetary Institution of Science Perm Federal Research Center of the Ural Branch of the Russian Academy of Sciences

Email: mishkasemin@gmail.com
78 Sibirskaya st., Building A, Perm, 614007, Russian Federation

PhD in Engineering, Research Fellow

### Oleg S. Parshakov

Mining Institute of the Ural Branch of the Russian Academy of Sciences - Branch of the Federal State Budgetary Institution of Science Perm Federal Research Center of the Ural Branch of the Russian Academy of Sciences

Email: olegparshakov@gmail.com
78 Sibirskaya st., Building A, Perm, 614007, Russian Federation

Junior Research Fellow

### Evgeniy V. Kolesov

Mining Institute of the Ural Branch of the Russian Academy of Sciences - Branch of the Federal State Budgetary Institution of Science Perm Federal Research Center of the Ural Branch of the Russian Academy of Sciences

Email: kolesovev@gmail.com
78 Sibirskaya st., Building A, Perm, 614007, Russian Federation

Engineer

# References

1. Romenskii A.A. Obosnovanie parametrov prokhodnicheskogo tsikla i ledoporodnogo ograzhdeniia pri stroitel'stve vertikal'nykh stvolov [Substantiation of the parameters of the tunnel and the ice rock wall in the construction of vertical trunks]. Ph. D. thesis. Moscow, 1983, 227 p.
2. Pravila bezopasnosti pri stroitel'stve podzemnykh sooruzhenii PB 03-428-02 [Safety rules for the construction of underground structures PB 03-428-02]. Utverzhdeny Postanovleniem Gosgortekhnadzora RF ot 02.11.2001 no. 49, 167 p.
3. Levin L.Iu., Semin M.A., Zaitsev A.V. Kontrol' i prognoz formirovaniia ledoporodnogo ograzhdeniia s ispol'zovaniem optovolokonnykh tekhnologii [Monitoring and forecasting the formation of the ice rock wall using fiber-optic technology]. Strategiia i protsessy osvoeniia georesursov. Sbornik nauchnykh trudov. Perm', 2016, pp.236-238.
4. Levin L.Iu., Zaitsev A.V., Semin M.A. Kontrol' teplovogo rezhima porodnogo massiva na osnove primeneniia optovolokonnykh tekhnologii monitoringa temperatur v skvazhinakh [Control of the thermal regime of the rock massif based on the application of fiber-optic technologies for monitoring temperatures in wells]. Gornoe ekho, 2016, no.1(62), pp.35-37.
5. Trupak N.G. Zamorazhivanie gornykh porod pri prokhodke stvolov [Freezing of rocks when drilling trunks]. Moscow, Ugletekhizdat, 1954, 896 p.
6. Trupak N.G. Zamorazhivanie gruntov v podzemnom stroitel'stve [Freezing of soils in underground construction]. Moscow, Nedra, 1974, 280 p.
7. Man'kovskii G.I. Spetsial'nye sposoby prokhodki gornykh vyrabotok [Special methods for excavating mine workings]. Moscow, Ugletekhizdat, 1958, 454 p.
8. Dolgov O.A. Metodika rascheta protsessa zamorazhivaniia gornykh porod pri prokhodke stvolov shakht sposobom zamorazhivaniia na bol'shuiu glubinu [Method for calculating the process of freezing rocks when drilling shafts through a method of freezing to great depths]. Zamorazhivanie gornykh porod pri prokhodke stvolov shakht. Moscow, Izdatel'stvo Akademii nauk SSSR, 1961, pp.9-64.
9. Dmitriev A.P., Goncharov S.A. Termodinamicheskie protsessy v gornykh porodakh [Thermodynamic processes in rocks]. Moscow, Nedra, 1990, 360 p.
10. Bel'ferman M.U. Temperaturnoe pole ledoporodnogo ograzhdeniia shakhtnykh stvolov pri dvukhriadnom raspolozhenii zamorazhivaiushchikh kolonok [The temperature field of the ice rock wall of shaft shafts with a two-row arrangement of freezing columns]. Voprosy organizatsii i mekhanizatsii gornoprokhodcheskikh rabot. Moscow, Institut gornogo dela imeni A.A. Skochinskogo, 1976, p.109-116.
11. Harris J.S. Ground freezing in Practice. Thomas Telford Limited, 1995, 290 p.
12. Andersland O.B., Ladanyi B. An introduction to frozen ground engineering. Springer US, 1994, 352 p.
13. Alifanov O.M. Obratnye zadachi teploobmena [Inverse heat transfer problems]. Moscow, Mashinostroenie, 1988, 280 p.
14. Tikhonov A.N., Arsenin V.Y. Solution of ill-posed problems. Washington, DC, Winston & Sons, 1977, 258 p.
15. Vabishchevich P.N., Vasil'eva M.V., Pavlova N.V. Chislennoe modelirovanie termostabilizatsii fil'truiushchikh gruntov [Numerical modeling of thermostabilization of filtering soils]. Matematicheskoe modelirovanie, 2014, vol.26, no.9, pp.111-125.
16. Musakaev N.G., Romaniuk S.N., Borodin S.L. Chislennoe issledovanie zakonomernostei dvizheniia fronta fazovogo perekhoda v mnogoletnemerzlykh porodakh [Numerical study of the laws governing the motion of the phase transition front in permafrost rocks]. Izvestiia vuzov. Neft' i gaz, 2011, no.6, pp.122-128.
17. Borodin S.L. Chislennye metody resheniia zadachi Stefana [Numerical methods for solving the Stefan problem] Tyumen State University Herald. Physical and Mathematical Modeling. Oil, Gas, Energy. 2015, vol.1, no.3, pp.164-175.
18. Goldman N. Inverse Stefan problems. Springer Science & Business Media, 2012, 412 p.
19. Mikkola M., Hartikainen J. Mathematical model of soil freezing and its numerical application. International Journal for Numerical Methods in Engeneering, 2001, vol.52, pp.543-557. doi: 10.1002/nme.300
20. Kruschwitz J., Bluhm J. Modeling of ice formation in porous solids with regard to the description of frost damage. Computational Material Science, 2005, vol.3-4, pp.407-417. doi: 10.1016/j.commatsci.2004.09.025
21. Razrabotka iskhodnykh dannykh dlia proekta prokhodki stvolov, v tom chisle: iskhodnye dannye po skipovomu stvolu [Development of initial data for the project of sinking of trunks, including initial data on the skip shaft]. Otchet o NIR № 58-12, book 1, etap 10.1. Belgorkhimprom, 2013, 192 p.
22. Analiz i obobshchenie rezul'tatov. Vyiavlenie zakonomernostei variatsii teplofizicheskikh i prochnostno-deformatsionnykh kharakteristik gornykh porod v vertikal'nom i gorizontal'nom napravleniiakh na uchastke petrikovskogo gorno-obogatitel'nogo kompleksa. Formirovanie bazy dannykh [Analysis and synthesis of the results. Identification of the regularities in the variation of the thermophysical and strength-deformation characteristics of rocks in the vertical and horizontal directions in the area of the Petrikov Ore Mining and Processing Complex. Formation of the database]. Otchet o NIR №58-12, etap 30.2.6. Belgorkhimprom, 2013, 230 p.
23. Karslou G., Eger D. Teploprovodnost' tverdykh tel [Thermal conductivity of solids]. Moscow, Nauka, 1964, 488 p.
24. Butcher J.C. Numerical methods for ordinary differential equations. New York, John Wiley & Sons, 2008, 440 p.
25. Maddox I.J. Elements of functional analysis. 2nd ed. Cambridge, 1988, 256 p.
26. Alexiades V., Solomon A.D. Mathematical modeling of melting and freezing processes. Washington DC, Hemisphere, 1993, 336 p.
27. Hu H., Argyropoulos S.A. Mathematical modelling of solidification and melting: a review. Modelling and Simulation in Materials Science and Engineering, 1996, vol. 4, pp.371-396. doi: 10.1088/0965-0393/4/4/004
28. Vasil'ev F.P. Metody optimizatsii [Optimization methods]. Moscow, Faktorial Press, 2002, 824 p.
29. Jarny Y., Ozisik M.N., Bardon J.P. A general optimization method using adjoint equation for solving multidimensional inverse heat conduction. International Journal of Heat and Mass Transfer, 1991, vol.34, no.11, pp.2911-2919. doi: 10.1016/0017-9310(91)90251-9
30. Ozisik M.N. Inverse heat transfer: fundamentals and applications. CRC Press, 2000, 352 p.

# Statistics

#### Views

Abstract - 259

PDF (Russian) - 42

PDF (English) - 32

### Refbacks

• There are currently no refbacks.

Copyright (c) 2017 Levin L.Y., Semin M.A., Parshakov O.S., Kolesov E.V.

## This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies