Full Text

Введение Возведение вертикальных шахтных стволов в слабых, неустойчивых водоносных горных породах требует проведения специальных мероприятий. Универсальным и надежным способом повышения прочности пород и преграждения проникновению грунтовых вод в строящуюся выработку является искусственное замораживание. Цель искусственного замораживания горных пород заключается в создании ледопородного ограждения, представляющего собой замороженный породный массив. Проектирование ледопородного ограждения требует научного обоснования предлагаемых технических решений. Большинство методов определения технологических параметров замораживающих скважин, прогнозирования геометрии и толщины формируемого ледопородного ограждения основано на упрощенных формулах и эмпирических коэффициентах [1-3], что не позволяет описать в полной мере процесс образования ледопородного ограждения и может приводить к значительному отклонению реальных значений от проектных. В связи с этим в настоящее время для исследования искусственного замораживания породного массива широко применяются методы численного моделирования [4-7]. Базовой моделью для процесса формирования ледопородного ограждения является двухфазная задача Стефана, которая описывает изменение температуры в сплошной среде с учетом фазового перехода первого рода. Среди наиболее эффективных и распространенных методов численного решения прямой задачи Стефана можно выделить методы с явным определением положения границы фазового перехода и методы с фиксированной сеткой, реализуемые с использованием конечно-разностного или конечно-элементного подходов [8, 9]. В первой группе методов положение фронта фазового перехода отслеживается на каждом временном шаге. В свою очередь, методы с фиксированной сеткой основаны на энтальпийной формулировке задачи Стефана и позволяют строить решение без явного выделения этой границы, что дает возможность сравнительно просто выполнять решение задач, поставленных в двух- и трехмерных областях со сложной геометрией [10, 11]. Для решения прямой задачи Стефана необходимо, чтобы в уравнении и условии Стефана были заданы все теплофизические параметры и были известны все функции, определяющие начальное и граничные условия. Однако в силу технологических особенностей искусственного замораживания горных пород определение всех функций и параметров, входящих в задачу Стефана, затруднительно или требует дорогостоящих и длительных экспериментальных исследований [12]. Ввиду этого при моделировании процесса формирования ледопородного ограждения возникает необходимость в постановке и решении обратных задач. Дополнительной информацией, необходимой для решения обратных задач, могут служить измерения температуры, проводимые в термических наблюдательных скважинах, которые бурятся на некотором расстоянии от контура замораживающих скважин для контроля за формированием ледопородного ограждения. Большое количество работ посвящено разработке методов численного решения граничных обратных двухфазных задач Стефана [13-29], связанных с определением граничных условий первого или второго рода или коэффициента теплоотдачи. В [13-21, 26-28] в качестве дополнительной информации задается изменение температуры от времени во внутренних точках исследуемой области и сведения о положении или движении фронта фазового перехода. В [22-25, 29] предполагается известной только информация об изменении температуры. В большинстве из указанных работ для постановки обратной задачи Стефана используется вариационный подход [13-25], а решение сводится к минимизации функционала невязки. Для минимизации функционала в [13-18] применяются алгоритмы, использующие коэффициенты чувствительности, в [19-22] представлены алгоритмы, в которых выполняется вычисление градиента функционала невязки при помощи решения сопряженной задачи. Эволюционные алгоритмы минимизации применяются в [23-25]. Численное решение сопутствующих прямых задач Стефана в [13-15, 19-21] выполняется численными методами с явным определением положения фронта фазового перехода, методы с фиксированной сеткой используются в [16-18, 22-25]. В [26-29] численное решение поставленных обратных задач выполняется бессеточными методами. При этом обратная задача сводится к системе алгебраических уравнений [26-28] или к минимизации целевой функции [29]. Исследованию и разработке алгоритмов численного решения коэффициентных обратных двухфазных задач Стефана посвящено значительно меньшее количество работ [30-32]. В работе [30] представлено математическое исследование граничных и коэффициентных обратных квазилинейных одномерных задач Стефана в вариационной постановке. Для решения обратных задач был разработан метод дескриптивной регуляризации, который основан на методе проекций сопряженных градиентов с вычислением градиента функционала невязки путем решения сопряженной задачи. Данный метод был применен для определения мощности теплового источника по известному распределению температуры в конечный момент времени. В работе [31] доказывается единственность решения этой задачи. В [32] рассматривается одномерная обратная задача Стефана, которая заключается в одновременном определении тепловых потоков на концах отрезка и постоянных коэффициентов теплопроводности в твердой и жидкой фазе по заданному изменению температуры от времени в нескольких внутренних точках. Алгоритм решения задачи основан на методе оптимизации Левенберга-Марквардта с использованием коэффициентов чувствительности. Настоящая статья посвящена разработке и сравнению эффективности двух алгоритмов численного решения коэффициентной обратной задачи Стефана в двумерной области. Цель работы состоит в определении объемной теплоемкости породного массива по изменению температуры от времени во внутренних точках исследуемой области. Решение обратной задачи выполняется на основе вариационного подхода с использованием энтальпийной формулировки для прямой задачи Стефана. Первый алгоритм основан на методе сопряженных градиентов, при этом для вычисления градиента функционала невязки решается сопряженная задача. Второй алгоритм основан на методе градиентного спуска с использованием коэффициентов чувствительности. Аналогичные методы успешно применяются для решения обратных коэффициентных задач теплопроводности в работах [33-38]. 1. Постановка задачи Процесс теплопереноса в горизонтальном сечении слоя породного массива, окружающего замораживающую скважину, описывается прямой задачей Стефана в энтальпийной постановке. Геометрия задачи представлена на рис. 1. Рис. 1. Геометрия расчетной области . Точки соответствуют расположению наблюдательных скважин Fig. 1. Geometry of the computational domain . The points correspond to the location of observation wells Искомая функция , определенная в замкнутой области (см. рис. 1) при , удовлетворяет уравнению (1) с начальным условием (2) и граничными условиями (3) где [К] - температурное поле в рассматриваемой области породного массива в момент времени , c; - распределение температуры в начальный момент времени; - граница замораживающей скважины; - внешняя граница расчетной области; - радиальные границы - температура на границе замораживающей скважины; - температура породного слоя на удалении от замораживающей скважины, остающаяся постоянной во все время охлаждения; - вектор внешней нормали к указанным границам; l [Вт/(м·К)] - теплопроводность среды, нижние индексы обозначают частные производные по соответствующим переменным. Принимая во внимание, что температурное поле вблизи замораживающей скважины является осесимметричным, исследуется только участок горизонтального сечения, в котором расположены наблюдательные скважины. Предполагается, что наблюдательные скважины оказывают малое влияния на распределение температуры в породном массиве, поэтому их радиус не учитывается. Функция [Дж/(кг·К)] выражает удельное теплосодержание, т.е. энтальпию, единицы объемы породного массива: (4) где [Дж/(м3·К)] - объемная теплоемкость породного массива в зоне льда или зоне охлаждения; [К] - температура фазового перехода; [кг/м3] - плотность льда; [Дж/кг] - удельная теплота кристаллизации воды. Функция коэффициента теплопроводности двухфазного материала [Вт/(м·К)] имеет вид (5) где [Вт/(м·К)] - теплопроводность среды в зоне льда или зоне охлаждения. Положение границы фазового перехода в момент времени определяется как множество решений уравнения [8]: (6) Для решения прямой задачи Стефана (1)-(3) применяется метод конечных элементов, реализованный в коммерческом пакете численного моделирования Ansys (Academic Research Mechanical license and CFD № 1064623). Подробное изложение алгоритмов численного решения прямой задачи Стефана в энтальпийной постановке может быть найдено в [4, 5, 8-11]. Для обеспечения сходимости метода конечных элементов функции (4) и (5) на малом отрезке подвергались сглаживанию. При использовании сглаживания третьего порядка коэффициенты в (4) имеют вид (7) коэффициенты в (3) записываются как (8) Частная коэффициентная обратная задача к задаче Стефана (1)-(3) заключается в определении объемной теплоемкости породного массива в зоне льда и зоне охлаждения по данным изменения температуры от времени, предоставленным наблюдательными скважинами. При этом все остальные теплофизические параметры и функции, входящие в задачу (1)-(3), предполагаются известными. Вариационная постановка обратной задачи формулируется следующим образом. Пусть известно изменение температуры от времени при в m внутренних точках : (9) Требуется найти значения коэффициентов в задаче (1)-(3), при которых функционал невязки достигает наименьшего значения: (10) Здесь - декартово произведение отрезков, в котором осуществляется поиск решения; - пространство функций, суммируемых с квадратом; - норма в ; - решение прямой задачи (1)-(3) с коэффициентами . Если один из коэффициентов известен, то обратная задача сводится к определению только неизвестного коэффициента. Несмотря на то, что прямая задача Стефана (1)-(3) может быть сведена к осесимметричной путем использования полярной системы координат, численное решение прямой задачи и коэффициентной обратной задачи выполняется в декартовой системе координат. Это связано с тем, что дополнительная информация изменения температуры от времени (9), необходимая для решения обратной задачи, задается в отдельных точках исследуемой области . Вследствие этого обратная задача не является осесимметричной. Также следует отметить, что представленные ниже алгоритмы могут быть использованы с небольшими изменениями для численного решения аналогичных коэффициентных обратных задач в двумерных областях с более сложной геометрией. 2. Алгоритмы численного решения обратной задачи Стефана В ходе исследования было разработано два алгоритма численного решения поставленной коэффициентной обратной задачи Стефана. Известно, что обратные задачи [30] являются неустойчивыми по входным данным. Для получения устойчивого решения в качестве численного метода минимизации функционала невязки (10) в первом алгоритме используется итерационный метод сопряженных градиентов, а во втором итерационный метод наискорейшего спуска [33, 34, 39]. В [39] показано, что данные методы оптимизации обладают регуляризующими свойствами. Важную роль при реализации методов оптимизации играют процедуры вычисления градиентов минимизируемых функционалов. Один из основных подходов к определению вида градиентов функционалов и расчету их значений заключается в построении и решении сопряженных задач. Первый этап построения сопряженной задачи заключается в определении вида задачи в приращениях температуры. Для прямой задачи Стефана (1)-(3) было установлено, что задача в приращениях температуры имеет следующий вид: (11) (12) (13) Отсюда путем применения метода множителей Лагранжа [33] был определен вид сопряженной задачи (14) (15) (16) где - дельта-функция. Также при построении сопряженной задачи были получены формулы для градиентов функционала (10) по переменным и : (17) Следует отметить, что сопряженная задача (14)-(16) записана относительно обратного времени. Однако заменой ее можно свести к задаче с прямым временем. Поставленная коэффициентная обратная задача заключается в определении объемной теплоемкости в зоне льда и зоне охлаждения. Путем проведения вычислительных экспериментов было установлено, что при одновременном определении объемных теплоемкостей не удается восстановить их значения с достаточным уровнем точности. Вследствие этого разработанные алгоритмы предполагают поочередное определение объемных теплоемкостей: сначала для зоны охлаждения, а затем для зоны льда. Алгоритм численного решения коэффициентной обратной задачи на основе итеративного метода сопряженных градиентов формулируется следующим образом. Пусть на шаге известно приближенное решение вариационной задачи (10). Тогда вычисление более точного приближенного решения осуществляется путем выполнения нижеизложенной процедуры. 1. Находится численное решение прямой задачи Стефана (1)-(3) с коэффициентами . 2. Находится численное решение сопряженной задачи (14)-(16), из (17) определяется градиент функционала невязки по переменной . 3. Вычисляется коэффициент сопряжения : (18) 4. Вычисляется направление спуска : (19) 5. Находится численное решение задачи в приращениях (11)-(13) с коэффициентами и , . 6. Вычисляется длина шага : (20) Определяется приближенное решение : (21) Здесь - оператор проектирования на отрезок ; - начальное приближение для коэффициента . Критерием остановки итерационного процесса служит выполнение следующего условия: (22) где - заданная точность. После того как коэффициент определен, аналогичным образом находится коэффициент . Предложенный алгоритм кроме решения прямой задачи также предполагает решение сопряженной задачи и задачи в приращениях температуры. Такой подход является особенно эффективным, когда целью решения обратной задачи является определение неизвестной функции, задающей начальное или граничное условие. Однако, когда требуется определить только постоянные коэффициенты, имеет смысл использовать более простой подход. С этой целью был разработан второй алгоритм численного решения поставленной коэффициентной обратной задачи. Основное отличие второго алгоритма от первого заключается в том, что минимизация функционала выполняется итеративным методом наискорейшего спуска, а градиент функционала находится с использованием коэффициента чувствительности и фиксированного приращения неизвестного коэффициента. Градиенты функционала по переменным и могут быть записаны в следующем виде: (23) где - коэффициенты чувствительности. При данный коэффициент вычисляется по формуле (при аналогично) , (24) где - фиксированные приращения, которые выбираются исходя из точности входных данных. Второй алгоритм, основанный на методе наискорейшего спуска, формулируется следующим образом. Пусть на шаге известно приближенное решение вариационной задачи (10). Тогда приближенное решение определятся путем выполнения нижеизложенной процедуры. 1. Находится численное решение прямой задачи Стефана (1)-(3) с коэффициентами . 2. Находится численное решение прямой задачи Стефана (1)-(3) с коэффициентами . 3. Из (23) вычисляется градиент и определяется направление спуска : (25) 4. Вычисляется длина шага : (26) 5. Определяется приближенное решение : (27) Критерием остановки итерационного процесса является условие (22). После определения коэффициента аналогичным образом выполняется поиск коэффициента . В (26) скалярное произведение и норма определены в пространстве , где - момент времени, при котором значение хотя бы одной функции уменьшится до величины . Соответственно, при определении значения объемной теплоемкости в зоне льда в (26) рассматривается пространство , где - момент времени, при котором все функции примут значение меньшее либо равное величине . Проведенные расчеты показали, что при замене отрезка на указанные происходит существенное увеличение скорости сходимости. Вычисление приближенных значений объемных теплоемкостей и параметров представленных алгоритмов осуществляется в системе компьютерной алгебры Wolfram Mathematica 10.3.0. Численное решение прямой задачи (1)-(3) выполняется на четырехугольной сетке. Для численного решения задачи в приращениях температуры (11)-(13) и сопряженной задачи (14)-(16) используется метод конечных элементов, реализованный в Ansys. Решения строятся на той же сетке, что и решение прямой задачи. Численное интегрирование осуществляется путем применения одномерного или двумерного варианта формулы трапеций. Численное дифференцирование от по в (17) выполняется с использованием конечных разностей. Поиск решений задач (11)-(13) и (14)-(16) выполняется исходя из следующих допущений. Из (4) и (5) следует, что в уравнениях (11) и (14) функция и частная производная принимают постоянные значения при всех , кроме , а частные производные , , не обращаются в ноль только при . Вследствие этого при решении задач (11)-(13) и (14)-(16) значения функции и частной производной в тех точках области , в которых , заменяются на константы, равные средним интегральным значениям этих функций на отрезке , а слагаемые, в которые входят частные производные , , , не учитываются. Таким образом, численное решение задач (11)-(13) и (14)-(16) сводится к численному решению задач теплопроводности с разрывными коэффициентами и внутренними источниками (стоками) тепла. При этом в задаче (14)-(16) предварительно выполняется замена . Проведенные вычислительные эксперименты показывают, что данные допущения позволяют получить достаточно точное решение поставленной коэффициентной обратной задачи. 3. Результаты вычислительных экспериментов Построенные вычислительные алгоритмы были применены для решения примеров коэффициентных обратных задач Стефана, приближенных к условиям Петриковского участка Старобинского месторождения калийных солей. Геометрические параметры области : радиус замораживающей скважины м, внешний радиус расчетной области м, . Температура на границе замораживающей скважины К, температура на внешней границе расчетной области К, начальное распределение температуры в породном массиве К. Теплофизические параметры: объемная теплоемкость породного массива в зоне льда Дж/(м3·К), в зоне охлаждения Дж/(м3·К); коэффициент теплопроводности породного массива в зоне льда Вт/(м·К), в зоне охлаждения Вт/(м·К); температура фазового перехода К, плотность льда кг/м3, удельная теплота плавления Дж/кг. Теплофизические параметры были взяты из исходных данных для проходки шахтных стволов рудника Петриковского участка Старобинского месторождения калийных солей и соответствуют глиняному слою, расположенному в интервале глубин 14,7 - 23,5 м. Наблюдательные скважины расположены во внутренних точках области с координатами: (28) Функции изменения температуры от времени были определены из решения прямой задачи (1)-(3) с указанными выше данными, рассмотренного в точках (28). Сглаживание функций энтальпии (4) и коэффициента теплопроводности (5) выполнялось на отрезке , где К, К. Длина промежутка в 1 К является наименьшей, при которой уравнение баланса энергии для прямой задачи Стефана (1)-(3) на каждом временном шаге выполняется с погрешностью, не превышающей 0,005 Дж. При этом была использована сетка, состоящая из четырехугольных элементов с длиной ребра, не превышающей 0,008 м. Общее количество элементов 2320. В критерии остановки итерационного процесса (22) точность была задана на уровне . Первый пример заключался в определении объемной теплоемкости для зоны охлаждения в предположении, что в зоне льда объемная теплоемкость известна. В качестве дополнительной информации использовалась функция , заданная при c. Во втором алгоритме было задано фиксированное приращение Дж/(м3·К). Результаты решения примера представлены в табл. 1. Для коэффициентов полученных в ходе решения примера первым алгоритмом, на рис. 1, 2 представлены графики численных решений прямой задачи, рассмотренные в точке и вдоль радиального направления в момент времени c. Там же изображены график функции и график численного решения прямой задачи с точно заданным коэффициентом. Таблица1 Полученные в первом (I) и втором (II) алгоритме значения корня из функционала невязки, коэффициента [Дж/(м3·К)] и относительной погрешности в зависимости от числа итераций Table 1 Values of square root from the discrepancy functional, the coefficient [J/(m3·K)] and the relative error obtained in the first algorithm (I) and the second algorithm (II) for various number of iterations I II 0 902,8 0,745 0,743 0 902,8 0,745 0,743 1 470,0 1,405 0,515 1 306,7 1,778 0,386 2 198,1 2,095 0,277 2 59,6 2,616 0,097 3 62,5 2,601 0,102 3 3,00 2,882 0,005 4 15,4 2,822 0,026 5 2,3 2,885 0,004 Второй пример заключался в определении объемной теплоемкости для зоны льда в предположении, что объемная теплоемкость для зоны охлаждения известна. В качестве дополнительной информации использовалась только функция , определенная при c. Фиксированное приращение Дж/(м3·К). Рис. 2. Графики изменения от времени в точке функции и численных решений прямой задачи Стефана (1)-(3) с коэффициентами (а); графики распределений вдоль радиального направления , , в момент времени c численных решений прямой задачи Стефана (1)-(3): для точного значения коэффициента (б); для коэффициентов Fig. 2. Plots of evolution in the point of the function and numerical solutions of the direct Stefan problem (1)-(3) with coefficients (a); plots of distribution along radial direction , , at time s of numerical solutions of the direct Stefan problem (1)-(3): for the exact value of the coefficient (b); for coefficients Результаты решения примера представлены в табл. 2. На рис. 3 изображены графики численных решений прямой задачи для коэффициентов , полученных в ходе решения примера первым алгоритмом. Таблица 2 Полученные в первом (I) и втором (II) алгоритме значения корня из функционала невязки, коэффициента [Дж/(м3·К)] и относительной погрешности в зависимости от числа итераций Table 2 Values of square root from the discrepancy functional, coefficient [J/(m3·K)] and the relative error obtained in the first algorithm (I) and the second algorithm (II) for various number of iterations I II 0 50,9 0,600 0,692 0 50,9 0,600 0,692 1 43,1 0,861 0,558 1 15,2 2,323 0,193 2 29,0 1,264 0,351 2 5,4 1,988 0,021 3 18,9 1,539 0,210 3 4,0 1,933 0,007 4 14,0 1,719 0,117 5 9,8 1,843 0,053 6 9,6 1,893 0,028 7 4,0 1,933 0,007 Рис. 3. Графики изменения от времени в точке функции и численных решений прямой задачи Стефана (1) с коэффициентами (а); графики распределений вдоль радиального направления , , в момент времени с численных решений прямой задачи Стефана (1): для точного значения коэффициента (б); для коэффициентов Fig. 3. Plots of evolution in the point of the function and numerical solutions of the direct Stefan problem (1) with coefficients (a); plots of distribution along radial direction , , s of numerical solutions of the direct Stefan problem (1): for the exact value of the coefficient (b); for coefficients Третий пример заключался в определении объемной теплоемкости для зоны льда и для зоны охлаждения. В начале находился коэффициент , в то время как коэффициент оставался равным начальному значению. При этом в качестве дополнительной информации использовалась только функция , заданная при c. После определения коэффициента был найден коэффициент с использованием только функции , заданной при c. Фиксированные приращения коэффициентов были такие же, как в первом и втором примерах. Проведенные расчеты показывают, что разработанные алгоритмы позволяют определить объемную теплоемкость для зоны льда и зоны охлаждения как по отдельности, так и вместе с относительной погрешностью менее 1 %. При этом согласно данным, приведенным в табл. 1-3, в ходе решений обратных задач вторым алгоритмом было совершено меньше итераций до выполнения критерия остановки итерационного процесса, чем при решении задач первым алгоритмом, т.е. скорость сходимости второго алгоритма выше, чем первого. С учетом того, что вычислительные затраты на решение прямой задачи, задачи в приращениях температуры и сопряженной задачи примерно одинаковы, а во втором алгоритме необходимо находить решение только двух прямых задач, поиск решения обратной задачи вторым алгоритмом выполняется быстрее. Из рис. 2, 3 видно, что значение коэффициента оказывает большее влияние на распределение температурного поля в породном массиве по сравнению со значением коэффициента . В связи с этим одновременное определение коэффициентов приводит к тому, что при уменьшении значения функционала до заданной точности значение коэффициента становится близко к точному, а значение коэффициента отличается от точного в несколько раз. При этом увеличение числа итераций не способствует ни значительному уменьшению значения функционала, ни корректировке значений коэффициентов. Таблица 3 Полученные в первом (I) и втором (II) алгоритме значения корня из функционала невязки, коэффициентов и [Дж/(м3·К)], относительной погрешности в зависимости от числа итераций Table 3 Values of square root from the discrepancy functional, coefficients and [J/(m3·K)] the relative error obtained in the first algorithm (I) and the second algorithm (II) for various number of iterations I II 0 905,7 0,745 0,743 0 905,7 0,745 0,743 1 460,8 1,429 0,507 1 307,0 1,779 0,386 2 195,8 2,107 0,272 2 52,7 2,653 0,084 3 61,5 2,613 0,098 3 2,5 2,893 0,001 4 15,1 2,828 0,023 5 3,0 2,890 0,002 0 53,1 0,600 0,692 0 53,0 0,600 0,692 1 37,6 1,017 0,478 1 11,2 2,071 0,064 2 30,7 1,267 0,349 2 6,5 1,994 0,024 3 19,4 1,568 0,195 3 7,2 1,967 0,010 4 15,0 1,694 0,130 4 7,5 1,958 0,006 7 9,2 1,775 0,088 5 7,1 1,971 0,012 8 7,4 1,835 0,058 6 3,7 1,944 0,002 9 10,0 1,881 0,034 10 3,9 1,937 0,005 В случае последовательного определения коэффициентов представленные в табл. 3 данные свидетельствуют о том, что при достижении функционалом требуемого уровня точности приближенные значения коэффициентов , отличаются от точных значений с малыми относительными погрешностями. Как видно из табл. 1, 3, замена точного значения коэффициента на начальное приближение привело к увеличению величины всего на 2,9, что позволило определить значение коэффициента как первым, так и вторым алгоритмом с малой относительной погрешностью и за такое же число итераций, что в первом примере. Использование найденных значений коэффициента вместо точных, как следует из табл. 2, 3, также привело к незначительному увеличению величины на 2,2 для первого алгоритма и на 2.1 для второго, но при определении значения коэффициента количество итераций, необходимых для достижения требуемого уровня точности, возросло по сравнению со вторым примером. Следует отметить, что такой подход к определению объемных теплоемкостей возможен благодаря тому, что разработанные алгоритмы основаны на регуляризованных методах оптимизации и являются устойчивыми к возмущениям в исходных данных. Из приведенных в табл. 3 данных можно заметить, что при уменьшении относительной погрешности между приближенным и точным значениями коэффициента значение функционала увеличилось. Примечательно, что оба алгоритма продолжили сходиться, несмотря на эту особенность, и значение коэффициента было определено с малой погрешностью. Однако во втором алгоритме возникновение данной особенности привело к нарушению монотонной сходимости приближенного решения к точному, в то время как в первом алгоритме монотонная сходимость сохранилась. В третьем примере последовательное определение значений объемных теплоемкостей проводилось на основе данных изменения температуры от времени в двух внутренних точках. В действительности можно обойтись измерениями температуры, проведенными только в одной точке. Главное, чтобы за время наблюдений снижение температуры в наблюдательной скважине при замораживании породного массива было достаточным для определения неизвестных значений. В данной работе при определении объемной теплоемкости для зоны охлаждения рассматривалась точка, в которой произошло снижение температуры примерно на 2 К, а восстановление объемной теплоемкости для зоны льда выполнялось по точке, в которой температура снизилась примерно на 20 К, при этом снижение температуры после фазового перехода составило примерно 8 К. Тем не менее установление зависимости между точностью определения значений объемных теплоемкостей и изменением температуры в наблюдательных скважинах требует дополнительных исследований. Заключение Представлен подход к исследованию процесса формирования ледопородного ограждения на основе решения двумерной коэффициентной обратной задачи Стефана. Основной целью работы является определение объемной теплоемкости для зоны льда и зоны охлаждения по данным изменения температуры от времени в ограниченном числе внутренних точек расчетной области. В работе предложены алгоритмы численного решения поставленной задачи. Первый алгоритм разработан на основе метода сопряженных градиентов, второй - на основе метода наискорейшего спуска. Для реализации первого алгоритма были установлены вид задачи в приращениях температуры, вид сопряженной задачи, вид градиента функционала невязки. Вследствие того, что значение объемной теплоемкости для зоны льда оказывает существенно меньшее влияние на распределение температуры в породном массиве по сравнению со значением объемной теплоемкости для зоны охлаждения, определение неизвестных теплофизических параметров происходит последовательно, начиная с объемной теплоемкости для зоны охлаждения. Эффективность разработанных алгоритмов была проиллюстрирована путем решения модельных примеров, заключающихся в определении объемной теплоемкости либо для одной зоны, либо для обеих зон. В результате было установлено, что для всех рассмотренных случаев оба алгоритма позволяют определить значения объемных теплоемкостей с относительной погрешностью менее 1 %. При этом скорость сходимости второго алгоритма выше, а вычислительные затраты ниже, чем у первого. Представленный подход и разработанные алгоритмы могут быть использованы для проектирования и уточнения исходных данных процесса формирования ледопородного ограждения при строительстве вертикальных шахтных стволов с использованием технологии искусственного замораживания породного массива.

About the authors

M S Zhelnin

Perm State National Research University

O A Plekhov

Institute of Continuous Media Mechanics

M A Semin

Mining Institute of UB RAS

L Yu Levin

Mining Institute of UB RAS

References

Statistics

Views

Abstract - 222

PDF (Russian) - 136