Numerical solution for an inverse problem about determination of volumetric heat capacity of rock mass during artificial freezing

Abstract


The paper is devoted to development and implementation of algorithms for a numerical solution of a coefficient inverse Stefan problem. This problem arises in a modeling for a process of ice wall formation around a projected horizontal section of a mine shaft. The ice wall is formed by an artificial ground freezing to convert soil pore water into ice. The temperature field modeling is based on an enthalpy form of a two-dimensional Stefan problem. The aim of the study is to determine the volumetric heat capacity for the rock layer on the base of additional information about temperature evolution in thermal wells. The problem of coefficients’ identification is stated as a variation form of the coefficient inverse Stefan problem. As a result, two algorithms for a numerical solution of the stated inverse Stefan problem have been developed. The first algorithm is based on the conjugate gradient iterative optimization method. The second algorithm is based on the steepest descent iterative optimization method. Under the first algorithm the calculation of the discrepancy functional gradient and determination of parameters for the optimization method are performed by solving a sensitivity problem and an adjoint problem. Forms of these problems have been obtained for the stated direct Stefan problem. For the second algorithm the discrepancy functional gradient and parameters for the descent step method are determined by calculating sensitivity coefficients. Solutions of the direct problem, the sensitivity problem and the adjoint problem are performed by the finite element method. The special feature of the used optimization methods is that these methods have regularizing properties. In order to verify effectiveness of the proposed algorithms, the computational experiments have been performed. The first and second experiments are related to determining only one unknown volumetric heat capacity for an ice domain or a cooling domain. The third experiment is devoted to determining the volumetric heat capacity for the both domains. The results of the experiments show that both algorithms allow to determine the volumetric heat capacity with a sufficiently good accuracy. However, a convergence rate of the second algorithm is higher than the rate of the first algorithm. The presented approach to modeling the process of ice wall formation as the coefficient inverse Stefan problem and the developed algorithms can be used for designing and improving the initial data for building mine shafts with a technology of artificial ground freezing.

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

  1. Трупак Н.Г. Замораживание грунтов в подземном строительстве. - М.: Недра. - 1974. - 280 с.
  2. Трупак Н.Г. Замораживание пород при сооружении вертикальных стволов шахт. - М.: Недра, 1983. - 170 с.
  3. Andersland O.B., Ladanyi B. Frozen ground engineering. - John Wiley & Sons, 2004. - 352 p.
  4. Вабищевич П.Н., Васильева М.В., Павлова Н.В. Численное моделирование термостабилизации фильтрующих грунтов // Математическое моделирование. - 2014. - Т. 26, № 9. - С. 111-125. doi: 10.1134/S2070048215020106
  5. Математическое моделирование искусственного замораживания грунтов / П.Н. Вабищевич [и др.] // Вычислительные технологии. - 2014. - Т. 19, № 4. - С. 19-31.
  6. Амосов П.В., Лукичев С.В., Наговицын О.В. Влияние пористости породного массива и температуры хладоносителя на скорость создания сплошного ледопородного ограждения // Вестн. Кольск. науч. центра РАН. - 2016. - № 4 (27). - С. 43-50.
  7. Левин Л.Ю., Зайцев А.В., Семин М.А. Контроль теплового режима породного массива на основе применения оптоволоконных технологий мониторинга температур в скважинах // Горное эхо. - 2016. - № 1. - С. 35-37.
  8. Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. - 2-е изд. - М: Либроком, 2009. - 782 с.
  9. Lewis R.W., Ravindran K. Finite element simulation of metal casting //International journal for numerical methods in engineering. - 2000. - Vol. 47. - No. 1-3. - P. 29-59. doi: 10.1002/(SICI)1097-0207(20000110/30)47:1/3<29::AID-NME760>3.0.CO;2-X
  10. Voller V.R., Swaminathan C.R., Thomas B.G. Fixed grid techniques for phase change problems: a review // International Journal for Numerical Methods in Engineering. - 1990. - Vol. 30. - No. 4. - P. 875-898. doi: 10.1002/nme.1620300419
  11. Крылов Д.А., Сидняев Н.И., Федотов А.А. Математическое моделирование распределения температурных полей // Математическое моделирование. - 2013. - Т. 25, № 7. - С. 3-27.
  12. Вакуленко И.С., Николаев П.В. Анализ и перспективы развития способа искусственного замораживания горных пород в подземном строительстве // Горный информационно-аналитический бюллетень (научно-технический журнал). - 2015. - № 3. - С. 338-346.
  13. Zabaras N., Ruan Y., Richmond O. Design of two-dimensional Stefan processes with desired freezing front motions // Numerical Heat Transfer, Part B Fundamentals. - 1992. - Vol. 21. - No. 3. - P. 307-325.
  14. Hematiyan M.R., Karami G. A boundary elements pseudo heat source method formulation for inverse analysis of solidification problems // Computational mechanics. - 2003. - Vol. 31. - No. 3. - P. 262-271. doi: 10.1007/s00466-003-0429-0
  15. Okamoto K., Li B.Q. A regularization method for the inverse design of solidification processes with natural convection // International journal of heat and mass transfer. - 2007. - Vol. 50. - No. 21. - P. 4409-4423. doi: 10.1016/j.ijheatmasstransfer.2006.10.019
  16. Voller V.R. Enthalpy method for inverse Stefan problems // Numerical Heat Transfer, Part B. Fundamentals. - 1992. - Vol. 21. - No. 1. - P. 41-55.
  17. Xu R., Naterer G.F. Inverse method with heat and entropy transport in solidification processing of materials // Journal of Materials Processing Technology. - 2001. - Vol. 112. - No. 1. - P. 98-108. doi: 10.1016/S0924-0136(01)00556-8
  18. Khosravifard A., Hematiyan M.R., Wrobel L.C. Simultaneous control of solidus and liquidus lines in alloy solidification // Engineering Analysis with Boundary Elements. - 2013. - Vol. 37. - No. 2. - Р. 211-224. doi: 10.1016/j.enganabound.2012.10.001
  19. Zabaras N., Kang S. On the solution of an ill-posed design solidification problem using minimization techniques in finite-and infinite-dimensional function spaces // International Journal for Numerical Methods in Engineering. - 1993. - Vol. 36. - No. 23. - P. 3973-3990.
  20. Kang S., Zabaras N. Control of the freezing interface motion in two-dimensional solidification processes using the adjoint method // International Journal for Numerical Methods in Engineering. - 1995. - Vol. 38. - No. 1. - P. 63-80. doi: 10.1002/nme.1620380105
  21. Hinze M., Ziegenbalg S. Optimal control of the free boundary in a two-phase Stefan problem // Journal of Computational Physics. - 2007. - Vol. 223. - No. 2. - P. 657-684. doi: 10.1016/j.jcp.2006.09.030
  22. Optimal operation of alloy material in solidification processes with inverse heat transfer / A.A. Nejad [et al.] // International Communications in Heat and Mass Transfer. - 2010. - Vol. 37. - No. 6. - P. 711-716. doi: 10.1016/j.icheatmasstransfer.2010.03.002
  23. Słota D. Identification of the cooling condition in 2-D and 3-D continuous casting processes // Numerical Heat Transfer. Part B: Fundamentals. - 2009. - Vol. 55. - No. 2. - P. 155-176. doi: 10.1080/10407790802605232
  24. Hetmaniok E., Słota D., Zielonka A. Experimental verification of immune recruitment mechanism and clonal selection algorithm applied for solving the inverse problems of pure metal solidification // International Communications in Heat and Mass Transfer. - 2013. - Vol. 47. - P. 7-14. doi: 10.1016/j.icheatmasstransfer.2013.07.009
  25. Hetmaniok E., Słota D., Zielonka A. Using the swarm intelligence algorithms in solution of the two-dimensional inverse Stefan problem // Computers & Mathematics with Applications. - 2015. - Vol. 69. - No. 4. - P. 347-361. doi: 10.1016/j.camwa.2014.12.013
  26. Application of meshfree methods for solving the inverse one-dimensional Stefan problem / K. Rashedi [et al.] // Engineering Analysis with Boundary Elements. - 2014. - Vol. 40. - P. 1-21. doi: 10.1016/j.enganabound.2013.10.013
  27. Johansson B.T., Lesnic D., Reeve T. A meshless method for an inverse two-phase one-dimensional linear Stefan problem // Inverse Problems in Science and Engineering. - 2013. - Vol. 21. - No. 1. - P. 17-33. doi: 10.1016/j.matcom.2014.03.004
  28. Johansson B.T., Lesnic D., Reeve T. A meshless regularization method for a two-dimensional two-phase linear inverse Stefan problem // Advances in Applied Mathematics and Mechanics. - 2013. - Vol. 5. - No. 06. - P. 825-845. doi: 10.1017/S2070073300001259
  29. Johansson B.T., Lesnic D., Reeve T. A meshless method for an inverse two-phase one-dimensional nonlinear Stefan problem // Mathematics and Computers in Simulation. - 2014. - Vol. 101. - P. 61-77. doi: 10.1016/j.matcom.2014.03.004
  30. Gol'dman N.L. Inverse Stefan Problems. - Springer Science & Business Media, 2012. - 250 p.
  31. Гольдман Н.Л. Классы единственности решения двухфазных коэффициентных обратных задач Стефана // Доклады Академии наук. - 2013. - Т. 449, № 5. - С. 507-512. doi: 10.7868/S0869565213110066
  32. Hafid M., Lacroix M. An inverse heat transfer method for predicting the thermal characteristics of a molten material reactor // Applied Thermal Engineering. - 2016. - Vol. 108. - P. 140-149. doi: 10.1016/j.applthermaleng.2016.07.087
  33. Alifanov O.M. Inverse heat transfer problems. - Springer Science & Business Media, 2012. - 348 p.
  34. Ozisik M.N. Inverse heat transfer: fundamentals and applications. - CRC Press, 2000. - 330 p.
  35. Hasanov A., Pektaş B., Erdem A. Comparative analysis of inverse coefficient problems for parabolic equations. Part I: adjoint problem approach // Inverse Problems in Science and Engineering. - 2011. - Vol. 19. - No. 5. - P. 599-615. doi: 10.1080/17415977.2011.579605
  36. Самарский А.А., Вабищевич П. Н. Численные методы решения обратных задач математической физики: учеб. пособие. - 2-е изд. - М., 2007. - 478 с.
  37. Колесник С.А. Метод идентификации нелинейных компонентов тензора теплопроводности анизотропных материалов // Математическое моделирование. - 2014. - Т. 26, № 2. - С. 119-132. doi: 10.1134/S2070048214050044
  38. Mohebbi F., Sellier M. Estimation of thermal conductivity, heat transfer coefficient, and heat flux using a three dimensional inverse analysis // International Journal of Thermal Sciences. - 2016. - Vol. 99. - P. 258-270. doi: 10.1016/j.ijthermalsci.2015.09.002
  39. Gilyazov S.F., Gol'dman N.L. Regularization of ill-posed problems by iteration methods. - Springer Science & Business Media, 2013. - 340 p.

Statistics

Views

Abstract - 82

PDF (Russian) - 66

Cited-By


PlumX


Copyright (c) 2017 Zhelnin M.S., Plekhov O.A., Semin M.A., Levin L.Y.

Creative Commons License
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License.

This website uses cookies

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

About Cookies