CC BY
38
Желнин М.С., Плехов О.А., Семин М.А., Левин Л.Ю. Численное решение обратной задачи определения объемной теплоемкости породного массива в процессе искусственного замораживания // Вестник Пермского национального исследовательского политехнического университета. Механика. - 2017. - № 4. - С. 56-75. DOI: 10.15593/perm.mech/2017.4.05
Zhelnin M.S., Plekhov O.A., Semin M.A., Levin L.Yu. Numerical solution for an inverse problem about determination of volumetric heat capacity of rock mass during artificial freezing. PNRPU Mechanics Bulletin, 2017, no. 4, pp. 56-75. DOI: 10.15593/perm.mech/2017.4.05
ВЕСТНИК ПНИПУ. МЕХАНИКА
№ 4,2017 PNRPU MECHANICS BULLETIN
http ://vestnik.pstu. ru/mechanics/ab out/inf/
001 10.15593/регш.шесЬ/2017.4.05 УДК 622.253.35:519.6
ЧИСЛЕННОЕ РЕШЕНИЕ ОБРАТНОЙ ЗАДАЧИ ОПРЕДЕЛЕНИЯ ОБЪЕМНОЙ ТЕПЛОЕМКОСТИ ПОРОДНОГО МАССИВА В ПРОЦЕССЕ ИСКУССТВЕННОГО ЗАМОРАЖИВАНИЯ
М.С. Желнин1, О.А. Плехов2, М.А. Семин3, Л.Ю. Левин3
Пермскии государственный национальным исследовательским университет, Пермь, Россия 2Институт механики сплошных сред УрО РАН, Пермь, Россия 3Горный институт УрО РАН, Пермь, Россия
О СТАТЬЕ
АННОТАЦИЯ
Получена: 11 августа 2017 г. Принята: 10 ноября 2017 г. Опубликована: 29 декабря 2017 г.
Ключевые слова:
искусственное замораживание грунтов, обратная задача Стефана, метод сопряженных градиентов, метод наискорейшего спуска, итеративная регуляризация, численное моделирование.
Статья посвящена разработке, реализации и сравнительному анализу эффективности двух алгоритмов численного решения коэффициентной обратной задачи Стефана, возникающей при моделировании процесса формирования ледопородного ограждения вокруг проектируемого сечения вертикального шахтного ствола. Образование ледопородного ограждения происходит в результате искусственного замораживания грунтовых вод, содержащихся в породном массиве. Моделирование распределения температурного поля в породном массиве основывается на двумерной задаче Стефана в энтальпийной постановке. Целью работы является определение объемной теплоемкости слоя породного массива по данным изменения температуры от времени в нескольких термических скважинах. Задача определения неизвестных теплофизических параметров формулируется в виде коэффициентной обратной задачи Стефана, записанной в вариационной постановке.
В результате исследования было разработано два алгоритма численного решения коэффициентной обратной задачи Стефана. Первый алгоритм основан на итеративном методе оптимизации сопряженных градиентов, второй - на итеративном методе оптимизации наискорейшего спуска. В первом алгоритме расчет градиента функционала невязки и определение параметров метода оптимизации происходит путем решения задачи в приращениях температуры и сопряженной задачи. Во втором алгоритме для расчета градиента функционала и параметров метода наискорейшего спуска используются коэффициенты чувствительности. Решение прямой задачи, задачи в приращениях температуры и сопряженной задачи выполняется методом конечных элементов. Особенность ме-
© Желнин Максим Сергеевич - магистр, e-mail: [email protected] Плехов Олег Анатольевич - доктор физико-математических наук, e-mail: [email protected] Семин Михаил Александрович - кандидат технических наук, e-mail: [email protected] Левин Лев Юрьевич - доктор технических наук, e-mail: [email protected]
Maxim S. Zhelnin - MD, e-mail: [email protected]
Oleg A. Plekhov - Doctor of Physical and Mathematical Sciences, e-mail: [email protected] Mikhail A. Semin - CSc in Technical Sciences, e-mail: [email protected] Lev Yu. Levin - Doctor of Technical Sciences, e-mail: [email protected]
тодов оптимизации, используемых в предложенных алгоритмах, заключается в том, что они обладают регуляризующими свойствами.
Сравнение эффективности предложенных алгоритмов проводилось при решении обратной коэффициентной задачи Стефана, заключающейся либо в определении объемной теплоемкости для зоны льда или зоны охлаждения, либо в определении объемных теплоемкостей для обеих зон. Полученные результаты свидетельствуют о том, что оба алгоритма позволяют получить решения поставленной коэффициентной обратной задачи Стефана с достаточно хорошим уровнем точности. Однако скорость сходимости второго алгоритма выше, чем первого.
Представленный подход к моделированию процесса формирования ледопо-родного ограждения как коэффициентной обратной задачи Стефана и разработанные алгоритмы могут быть использованы для проектирования и уточнения исходных данных при строительстве вертикальных шахтных стволов с использованием технологии искусственного замораживания.
© ПНИПУ
NUMERICAL SOLUTION FOR AN INVERSE PROBLEM ABOUT DETERMINATION OF VOLUMETRIC HEAT CAPACITY OF ROCK MASS DURING ARTIFICIAL FREEZING
M.S. Zhelnin1, O.A. Plekhov2, M.A. Semin3, L.Yu. Levin3
Perm State National Research University, Perm, Russian Federation
2
Institute of Continuous Media Mechanics, UB RAS, Perm, Russian Federation
3
Mining Institute of UB RAS, Perm, Russian Federation
ARTICLE INFO 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.
© PNRPU
Received: 11 August 2017 Accepted: 10 November 2017 Published: 29 December 2017
Keywords:
ground artificial freezing, inverse Stefan problem, conjugate gradient method, steepest descent method, iterative regularization, numerical modeling.
Введение
Возведение вертикальных шахтных стволов в слабых, неустойчивых водоносных горных породах требует проведения специальных мероприятий. Универсальным и надежным способом повышения прочности пород и преграждения проникновению грунтовых вод в строящуюся выработку является искусственное замораживание. Цель искусственного замораживания горных пород заключается в создании ледопородного ограждения, представляющего собой замороженный породный массив.
Проектирование ледопородного ограждения требует научного обоснования предлагаемых технических решений. Большинство методов определения технологических параметров замораживающих скважин, прогнозирования геометрии и толщины формируемого ледопородного ограждения основано на упрощенных формулах и эмпирических коэффициентах [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. Геометрия расчетной области Q. Точки Wt (i = 1,2) соответствуют расположению наблюдательных скважин Fig. 1. Geometry of the computational domain Q. The points Wt (i = 1,2) correspond to the location of observation wells
Искомая функция u = u(x, y, t), определенная в замкнутой области Q (см. рис. 1) при t е [0, T], удовлетворяет уравнению
Ht(u) = (X(u)ux)x + (^u)uy)y, (x,y) e Q, t e[0,T] (1)
с начальным условием
u
t=0
= Ф, (x, y) eQ
и граничными условиями
u |г = vw, u |r
= v„
du dn
= 0, t e [0,T], k = 1,2,
(2)
(3)
где и [К] - температурное поле в рассматриваемой области О породного массива в момент времени t, с; ф(х, у) - распределение температуры в начальный момент времени; Гк - граница замораживающей скважины; Гиит - внешняя граница расчетной области; Г, - радиальные границы (к = 1,2); ^) - температура на границе Г№ замораживающей скважины; - температура породного слоя на удалении от замораживающей скважины, остающаяся постоянной во все время охлаждения; п - вектор внешней нормали к указанным границам; X [Вт/(м К)] - теплопроводность среды, нижние индексы обозначают частные производные по соответствующим переменным.
Принимая во внимание, что температурное поле вблизи замораживающей скважины является осесимметричным, исследуется только участок горизонтального сечения, в котором расположены наблюдательные скважины. Предполагается, что наблюдательные скважины оказывают малое влияния на распределение температуры в породном массиве, поэтому их радиус не учитывается.
Функция Н(и) [Дж/(кгК)] выражает удельное теплосодержание, т.е. энтальпию, единицы объемы породного массива:
H (u ) =
Csdu, u < u1 < uph,
I h
u
u1 < u < u2,
(4)
Clq (u " uph ) + Csduph + PA
sd ph
uph < u2 <U,
где Ск (к = ъй, \ц) [Дж/(м3К)] - объемная теплоемкость породного массива в зоне льда или зоне охлаждения; ирк [К] - температура фазового перехода; рг. [кг/м3] - плотность льда; Ь [Дж/кг] - удельная теплота кристаллизации воды.
г
k
1=0
Функция коэффициента теплопроводности двухфазного материала k(u) [Вт/(мК)] имеет вид
k(u) =
k d, и < щ < uph,
3
^ 1Щ, u1 <и < и2, (5)
К , Uph < u2 < u
где (& = sd, lq) [Вт/(мК)] - теплопроводность среды в зоне льда или зоне охлаждения.
Положение границы фазового перехода Гph (t0) в момент времени t0 определяется как множество решений уравнения [8]:
H(u(*, y, to)) = H U + H {Щ) = const. (6)
Для решения прямой задачи Стефана (1)-(3) применяется метод конечных элементов, реализованный в коммерческом пакете численного моделирования Ansys (Academic Research Mechanical license and CFD № 1064623). Подробное изложение алгоритмов численного решения прямой задачи Стефана в энтальпийной постановке может быть найдено в [4, 5, 8-11].
Для обеспечения сходимости метода конечных элементов функции (4) и (5) на малом отрезке [и1з и2] подвергались сглаживанию. При использовании сглаживания третьего порядка коэффициенты h в (4) имеют вид
= u12(3u2uph ~2u2 ~U1Uph) (C _C ) + иКщ _3u2) L
h0 = , 4 3 (Clq Csd) + , 4 3 LPi,
(u1 _ u2) (u1 _ u2)
= u1 (uf + u1u2 + 4u22 - 6u2uph) C u2 (4Uj2 + u1u2 + u22 - 6uuph) c 6u1u2 T
h = 7 73 Clq 7 73 Csd +7 733
(u1 _ u2) (u1 _ u2) (u1 u2)
h = (3u1Uph + 3U2Uph _ 2U' _ 2u1u2 _ 2u22) (C _ C ) _ 3(U1 + u2)
h2 , 43 (Clq Csd) / 43
(u1 _ u2) (u1 _ u2)
(u1 + u2 _ 2uph) 2
h3 =-7-75-(Clq _ Csd ) + 7-73 TPi ,
(u1 _ u2) (u1 _ u2)
коэффициенты ki в (3) записываются как
(7)
/ = u1 (u1 3u2) k _ u2 (u2 3u1) k / = 6u1u2 (k _ - )
'0 " (u _u )3 lq (u _u )3 sd, ^ " (u _u )3 ( lq KsdЛ
m2Z U2) Vm1 U2)
3(u1 + u2)
/ \3 "lq ~~sas? '3 s \3
(u1 _u2) (u1 _u2)
(8)
l2 =_ t^-7t(klq _ksd X l3 =~,-73 (klq _ksd ).
Частная коэффициентная обратная задача к задаче Стефана (1)-(3) заключается в определении объемной теплоемкости породного массива в зоне льда и зоне охлаждения по данным изменения температуры от времени, предоставленным наблюдательными скважинами. При этом все остальные теплофизические параметры и функции, входящие в задачу (1)—(3), предполагаются известными.
1=0
Вариационная постановка обратной задачи формулируется следующим образом. Пусть известно изменение температуры от времени У1 = У1 ^) при t е[0, Т] в т внутренних точках (х1, yi) еП:
и|( ) = У,, i = 1,...,т. (9)
1( х,, у1) 1 4 у
Требуется найти значения (С^, С^ ) коэффициентов (Сй, С1д) в задаче (1)-(3), при которых функционал невязки J достигает наименьшего значения:
т
(С С^ х1 J С , ^), J С , С,) = !}\и(Сй , С, )|(ху) - ¥1 И1[0,Т ] . (10)
Здесь 1Сй х с М2 - декартово произведение отрезков, в котором осуществляется поиск решения; Ь2[0,Т] - пространство функций, суммируемых с квадратом; ||-||Ь2[0Т] - норма в Ь2[0,Т]; и - решение прямой задачи (1)-(3) с коэффициентами (Сй,С1д).
Если один из коэффициентов известен, то обратная задача сводится к определению только неизвестного коэффициента.
Несмотря на то, что прямая задача Стефана (1)-(3) может быть сведена к осесиммет-ричной путем использования полярной системы координат, численное решение прямой задачи и коэффициентной обратной задачи выполняется в декартовой системе координат. Это связано с тем, что дополнительная информация изменения температуры от времени (9), необходимая для решения обратной задачи, задается в отдельных точках исследуемой области П . Вследствие этого обратная задача не является осесимметричной.
Также следует отметить, что представленные ниже алгоритмы могут быть использованы с небольшими изменениями для численного решения аналогичных коэффициентных обратных задач в двумерных областях с более сложной геометрией.
2. Алгоритмы численного решения обратной задачи Стефана
В ходе исследования было разработано два алгоритма численного решения поставленной коэффициентной обратной задачи Стефана. Известно, что обратные задачи [30] являются неустойчивыми по входным данным. Для получения устойчивого решения в качестве численного метода минимизации функционала невязки (10) в первом алгоритме используется итерационный метод сопряженных градиентов, а во втором итерационный метод наискорейшего спуска [33, 34, 39]. В [39] показано, что данные методы оптимизации обладают регуляризующими свойствами.
Важную роль при реализации методов оптимизации играют процедуры вычисления градиентов минимизируемых функционалов. Один из основных подходов к определению вида градиентов функционалов и расчету их значений заключается в построении и решении сопряженных задач.
Первый этап построения сопряженной задачи заключается в определении вида задачи в приращениях температуры.
Для прямой задачи Стефана (1)-(3) было установлено, что задача в приращениях температуры имеет следующий вид:
ни (и) Ди + (и) ДС^ + (и )АС1д = X, (и )(Дихх + Аи^,) +
+2(Хи (и)их Дих + Хи (и)иу Диу ) + (11)
(и)(и„ + и№) + Хии(и)(их2 + и2у) -Ни((и)]Ди, (х, у) 6П, г е [0,Т],
Ди и = 0, (х, у) еО, (12)
дДи
Ди |г = 0, Ди |г = 0,
м пит
дп
= 0, г е [0, Т], к = 1,2. (13)
Отсюда путем применения метода множителей Лагранжа [33] был определен вид сопряженной задачи
Ни + Х(и XV хх +Ууу) +
уу) т
+ 2^ Ц^ - Шх - X, у - у ) = 0, (х, у) еО, г е[0, Т],
(14)
V 1г = 0, V 1г = 0, —
^ м -1- пит дп
Vи = 0, (х,у)еО, (15)
дv
= 0, г е [0, Т], к = 1,2, (16)
где 8(х - х^, у - ) - дельта-функция. Также при построении сопряженной задачи были получены формулы для градиентов функционала (10) по переменным и С1д:
gradсk^ (, ск ) = -£ |пНСгу&Шу, к = , Щ. (П)
Следует отметить, что сопряженная задача (14)-(16) записана относительно обратного времени. Однако заменой т=Т-г ее можно свести к задаче с прямым временем.
Поставленная коэффициентная обратная задача заключается в определении объемной теплоемкости в зоне льда и зоне охлаждения. Путем проведения вычислительных экспериментов было установлено, что при одновременном определении объемных тепло-емкостей не удается восстановить их значения с достаточным уровнем точности. Вследствие этого разработанные алгоритмы предполагают поочередное определение объемных теплоемкостей: сначала для зоны охлаждения, а затем для зоны льда.
Алгоритм численного решения коэффициентной обратной задачи на основе итеративного метода сопряженных градиентов формулируется следующим образом. Пусть на шаге 5 (з = 0,1,...) известно приближенное решение С(щ вариационной задачи (10). Тогда вычисление более точного приближенного решения С1(5+1) осуществляется путем выполнения нижеизложенной процедуры.
1. Находится численное решение и прямой задачи Стефана (1)-(3) с коэффициентами (( с*)).
2. Находится численное решение V сопряженной задачи (14)-(16), из (17) определяется градиент gradCI 3 функционала невязки по переменной С1д.
г
к
г=1
г
к
3. Вычисляется коэффициент сопряжения Ps:
gradQ/ (C? )■( gradc/ (C^, Cj-1))-grad^ (C^, Cj)))
Ps =•
grad^ (( CT)
n(s-1)
s = 1,2,... .
(19)
4. Вычисляется направление спуска r(s):
r(0) = grad Clq J ( C£»), r(s) = grad Clq J (C^, C«
5. Находится численное решение Ли задачи в приращениях (11)—(13) с коэффициен тами C0), C\sq)) и ЛCsd = О, ЛС„ = r(s).
6. Вычисляется длина шага а :
а, =
i =1
/ I
Z(u((cq)) -Y,Ли((с»)
(ч, y)
L2[0,T ]
m .
Z Ли(C«)
i=1
(ч, У)
(20)
L2[0,T ]
Определяется приближенное решение Ct
(s+1) : lq ■
C\r = piciq (C
(s) -а r(s)
), s = 0,1,... .
(21)
Здесь Р1 - оператор проектирования на отрезок 1С ; С^ е 1С - начальное приближение для коэффициента С,,(. Критерием остановки итерационного процесса служит выполнение следующего условия:
С«)<в, 8>О, (22)
где в - заданная точность.
После того как коэффициент С1д определен, аналогичным образом находится коэффициент С( .
Предложенный алгоритм кроме решения прямой задачи также предполагает решение сопряженной задачи и задачи в приращениях температуры. Такой подход является особенно эффективным, когда целью решения обратной задачи является определение неизвестной функции, задающей начальное или граничное условие. Однако, когда требуется определить только постоянные коэффициенты, имеет смысл использовать более простой подход. С этой целью был разработан второй алгоритм численного решения поставленной коэффициентной обратной задачи. Основное отличие второго алгоритма от первого заключается в том, что минимизация функционала J выполняется итеративным методом наискорейшего спуска, а градиент функционала находится с использованием коэффициента чувствительности и фиксированного приращения неизвестного коэффициента.
Градиенты функционала J
по переменным С5( и С^ могут быть записаны в следующем виде:
m 1
gradckJ(Csd, Cq ) = 2X J (u(Csd, Cq ) - Y ) uck ^)dt, k = sd, lq,
(23)
где uC
(ч, У)
- коэффициенты чувствительности. При k = lq данный коэффициент вычис-
ляется по формуле (при k = sd аналогично)
k
uC
_ u (Csd, Cq + ACsd ))) - u (Csd, Clq ) (xi,») ACiq
^, (24)
где ДС1д - фиксированные приращения, которые выбираются исходя из точности входных данных.
Второй алгоритм, основанный на методе наискорейшего спуска, формулируется следующим образом. Пусть на шаге 5 (5 = 0,1,.) известно приближенное решение С) вариационной задачи (10). Тогда приближенное решение С^5+1) определятся путем выполнения нижеизложенной процедуры.
1. Находится численное решение и прямой задачи Стефана (1)-(3) с коэффициентами (С50), сЦ)).
2. Находится численное решение и прямой задачи Стефана (1)-(3) с коэффициентами (С^, с%) +ДС1д).
3. Из (23) вычисляется градиент gradCJ и определяется направление спуска г(5):
г (5) = gradс/(C5(0), С%)). (25)
4. Вычисляется длина шага а5:
' Г (5)
а 5 =-т-2-. (26)
S[u(Cj™C-)|(_) -Y,u,
m
х
i _1
uC
r (s)
( , yi )
L2l0,tq ]
5. Определяется приближенное решение Cl(qs+1):
Cqs+1) _P^C-asr(s)), s _ 0,1,.... (27)
Критерием остановки итерационного процесса является условие (22). После определения коэффициента Clq аналогичным образом выполняется поиск коэффициента Csd .
В (26) скалярное произведение и норма определены в пространстве L2[0, tlq], где tlq — момент времени, при котором значение хотя бы одной функции Yi (i _ 1,..., m) уменьшится до величины u2. Соответственно, при определении значения объемной теплоемкости в зоне льда в (26) рассматривается пространство L2[tsd,T], где tsd — момент времени, при котором все функции Y (i _ 1,.,m) примут значение меньшее либо равное величине u1. Проведенные расчеты показали, что при замене отрезка [0,T] на указанные происходит существенное увеличение скорости сходимости.
Вычисление приближенных значений объемных теплоемкостей и параметров представленных алгоритмов осуществляется в системе компьютерной алгебры Wolfram Ma-thematica 10.3.0. Численное решение прямой задачи (1)—(3) выполняется на четырехугольной сетке. Для численного решения задачи в приращениях температуры (11)—(13) и сопряженной задачи (14)—(16) используется метод конечных элементов, реализованный в Ansys. Решения строятся на той же сетке, что и решение прямой задачи.
Численное интегрирование осуществляется путем применения одномерного или двумерного варианта формулы трапеций. Численное дифференцирование от Нс^ по t
в (17) выполняется с использованием конечных разностей.
Поиск решений задач (11)—(13) и (14)-(16) выполняется исходя из следующих допущений. Из (4) и (5) следует, что в уравнениях (11) и (14) функция X и частная производная Ни принимают постоянные значения при всех и , кроме и е (м1з и2], а частные производные Хи, Хии, Нш не обращаются в ноль только при и е (и1зи2]. Вследствие этого при решении задач (11)—(13) и (14)-(16) значения функции X и частной производной Ни в тех точках области О, в которых и е (и1з и2], заменяются на константы, равные средним интегральным значениям этих функций на отрезке [и1з и2 ], а слагаемые, в которые входят частные производные Хи, Хии, Ни{, не учитываются.
Таким образом, численное решение задач (11)—(13) и (14)-(16) сводится к численному решению задач теплопроводности с разрывными коэффициентами и внутренними источниками (стоками) тепла. При этом в задаче (14)-(16) предварительно выполняется замена т= Т — t.
Проведенные вычислительные эксперименты показывают, что данные допущения позволяют получить достаточно точное решение поставленной коэффициентной обратной задачи.
3. Результаты вычислительных экспериментов
Построенные вычислительные алгоритмы были применены для решения примеров коэффициентных обратных задач Стефана, приближенных к условиям Петриковского участка Старобинского месторождения калийных солей.
Геометрические параметры области О: радиус замораживающей скважины Як =| ОА |= 0,1 м, внешний радиус расчетной области Япит =| ОВ |= 1,26 м, ¿СОВ = % /32.
Температура на границе замораживающей скважины = 253 К, температура на внешней границе расчетной области = 283 К, начальное распределение температуры в породном массиве ф = 283 К.
Теплофизические параметры: объемная теплоемкость породного массива в зоне льда С( = 1,947-106 Дж/(м3К), в зоне охлаждения С1д = 2,896-106 Дж/(м3К); коэффициент теплопроводности породного массива в зоне льда X ( = 1,5 Вт/(мК), в зоне охлаждения Х1д = 1,29 Вт/(м К); температура фазового перехода ирк = 271 К, плотность льда
р.. = 916 кг/м3, удельная теплота плавления Ь = 3,3-105 Дж/кг. Теплофизические параметры были взяты из исходных данных для проходки шахтных стволов рудника Петриков-ского участка Старобинского месторождения калийных солей и соответствуют глиняному слою, расположенному в интервале глубин 14,7 - 23,5 м.
Наблюдательные скважины расположены во внутренних точках области О с координатами:
X Д) = (0,6272; 0,0308), {х2, У2) = (0,1718; 0,0084). (28)
Функции Yi (i = 1,2) изменения температуры от времени были определены из решения прямой задачи (1)-(3) с указанными выше данными, рассмотренного в точках (28).
Сглаживание функций энтальпии (4) и коэффициента теплопроводности (5) выполнялось на отрезке \их,и2], где щ = 270,5 К, и2 = 271,5 К. Длина промежутка [м1,и2]в 1 К
является наименьшей, при которой уравнение баланса энергии для прямой задачи Стефана (1)—(3) на каждом временном шаге выполняется с погрешностью, не превышающей 0,005 Дж. При этом была использована сетка, состоящая из четырехугольных элементов с длиной ребра, не превышающей 0,008 м. Общее количество элементов 2320.
В критерии остановки итерационного процесса (22) точность была задана на уровне 8 = 16.
Первый пример заключался в определении объемной теплоемкости для зоны охлаждения в предположении, что в зоне льда объемная теплоемкость известна. В качестве дополнительной информации использовалась функция YJ, заданная при 0 < t < 2-105 c. Во втором алгоритме было задано фиксированное приращение AC1q = 1-105 Дж/(м3К).
Результаты решения примера представлены в табл. 1. Для коэффициентов C(fq) (s = 0,5), полученных в ходе решения примера первым алгоритмом, на рис. 1, 2 представлены графики численных решений и(s) прямой задачи, рассмотренные в точке (лх, д) и вдоль радиального
направления в момент времени t = 2 • 105 c. Там же изображены график функции Y и график численного решения uext прямой задачи с точно заданным коэффициентом.
Таблица1
Полученные в первом (I) и втором (II) алгоритме значения корня 41 из функционала невязки, коэффициента C^) [Дж/(м3К)] и относительной погрешности 8 в зависимости от числа итераций s
Table 1
Values of square root from the discrepancy functional, the coefficient C(lq) [J/(m3K)]
and the relative error 8 obtained in the first algorithm (I) and the second algorithm (II)
for various number of iterations s
I II
s VJ Cq •ю-6 5 s 41 C^ •Ю-6 5
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
Второй пример заключался в определении объемной теплоемкости для зоны льда в предположении, что объемная теплоемкость для зоны охлаждения известна. В качестве дополнительной информации использовалась только функция У2, определенная при
0 < / < 4-105 с. Фиксированное приращение ЛС^ = 1-105 Дж/(м3 К).
Рис. 2. Графики изменения от времени в точке (xl,л) = (0,6272; 0,0308) функции Yl и численных решений u(s) прямой задачи Стефана (1)-(3) с коэффициентами C^ (s = 0,5) (а); графики распределений вдоль радиального направления r , Rw < r < Rnum , в момент времени t = 2 • 105 c численных решений прямой задачи Стефана (1)—(3): uext для точного значения коэффициента Clq (б);
u( s) для коэффициентов Cllqs) (s = 0,5) Fig. 2. Plots of evolution in the point ( x1, y1) = (0.6272,0.0308) of the function Y1 and numerical solutions u(s) of the direct Stefan problem (1)-(3) with coefficients Cq) (s = 0,5) (a); plots of distribution along radial direction r, Rw < r < Rnum, at time t = 2 -105 s of numerical solutions of the direct Stefan problem (1)-(3): uext for the exact value of the coefficient Clq (b); u(; s) for coefficients C^ (s = 0,5)
Результаты решения примера представлены в табл. 2. На рис. 3 изображены графики численных решений u( s ) прямой задачи для коэффициентов C(d ) (s = 0,7), полученных в ходе решения примера первым алгоритмом.
Таблица 2
Полученные в первом (I) и втором (II) алгоритме значения корня из функционала невязки, коэффициента C/J [Дж/(м3 К)] и относительной погрешности ô в зависимости от числа итераций s
Table 2
Values of square root VJ from the discrepancy functional, coefficient C{ssd} [J/(m3 K)]
and the relative error ô obtained in the first algorithm (I) and the second algorithm (II)
for various number of iterations s
I II
s 41 C« •ю-6 ô s 41 C« •Ю-6 ô
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. Графики изменения от времени в точке (x2,y2) = (0,1718; 0,0084) функции Y2 и численных решений u(s) прямой задачи Стефана (1) с коэффициентами C^J (s = 0,7) (а); графики распределений вдоль радиального направления r , Rw < r < Rnum , в момент времени t = 4 -105 с численных решений прямой задачи Стефана (1): uext для точного значения коэффициента Csd (б);
u(s) для коэффициентов C(sd (s = 0,7) Fig. 3. Plots of evolution in the point (x2, y2) = (0.1718; 0.0084) of the function Y2 and numerical solutions u(s) of the direct Stefan problem (1) with coefficients CS'J (s = 0,7) (a); plots of distribution along radial direction r, Rw < r < Rnum, t = 4 -105 s of numerical solutions of the direct Stefan problem (1): uxt for the exact value of the coefficient Csd (b); u(s) for coefficients C^ (s = 0,7)
Третий пример заключался в определении объемной теплоемкости для зоны льда и для зоны охлаждения. В начале находился коэффициент Clq, в то время как коэффициент Csd оставался равным начальному значению. При этом в качестве дополнительной информации использовалась только функция Y¡, заданная при 0 < t < 2 -105 c. После определения коэффициента Clq был найден коэффициент Csd с использованием только функции Y2, заданной при 0 < t < 4 -105 c. Фиксированные приращения коэффициентов были
такие же, как в первом и втором примерах.
Проведенные расчеты показывают, что разработанные алгоритмы позволяют определить объемную теплоемкость для зоны льда и зоны охлаждения как по отдельности, так и вместе с относительной погрешностью менее 1 %. При этом согласно данным, приведенным в табл. 1-3, в ходе решений обратных задач вторым алгоритмом было совершено меньше итераций до выполнения критерия остановки итерационного процесса, чем при решении задач первым алгоритмом, т.е. скорость сходимости второго алгоритма выше, чем первого. С учетом того, что вычислительные затраты на решение прямой задачи, задачи в приращениях температуры и сопряженной задачи примерно одинаковы, а во втором алгоритме необходимо находить решение только двух прямых задач, поиск решения обратной задачи вторым алгоритмом выполняется быстрее.
Из рис. 2, 3 видно, что значение коэффициента Clq оказывает большее влияние на распределение температурного поля в породном массиве по сравнению со значением коэффициента Csd . В связи с этим одновременное определение коэффициентов приводит к тому, что при уменьшении значения функционала J до заданной точности значение коэффици-
ента Clq становится близко к точному, а значение коэффициента Csd отличается от точного
в несколько раз. При этом увеличение числа итераций не способствует ни значительному уменьшению значения функционала, ни корректировке значений коэффициентов.
Таблица 3
Полученные в первом (I) и втором (II) алгоритме значения корня 41 из функционала невязки, коэффициентов C^) и C^ [Дж/(м3К)], относительной погрешности 8 в зависимости от числа итераций s
Table 3
Values of square root 41 from the discrepancy functional, coefficients C^) and Cd) [J/(m3K)] the relative error 8 obtained in the first algorithm (I) and the second algorithm (II) for various number of iterations s
I II
s 41 C« -10-6 5 s 41 cqs) -ю-6 5
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 1 95,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
s 41 Cd -ю-6 5 s 41 c<d -ю-6 5
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 данные свидетельствуют о том, что при достижении функционалом J требуемого уровня точности приближенные значения коэффициентов С^ , С1д отличаются от точных значений с малыми относительными погрешностями. Как видно из табл. 1, 3, замена точного значения коэффициента С^ на начальное приближение С(0 привело к увеличению величины всего на 2,9, что позволило определить значение коэффициента С1д как первым, так и вторым алгоритмом с малой относительной погрешностью и за такое же число итераций, что в первом примере. Использование найденных значений коэффициента С1д вместо точных, как следует из табл. 2, 3, также привело к незначительному увеличению величины 4J на 2,2 для первого алгоритма и на 2.1 для второго, но при определении значения коэффициента Сы количество итераций, необходимых для достижения требуемого уровня точности, возросло по сравнению со вторым примером.
Следует отметить, что такой подход к определению объемных теплоемкостей возможен благодаря тому, что разработанные алгоритмы основаны на регуляризованных методах оптимизации и являются устойчивыми к возмущениям в исходных данных.
Из приведенных в табл. 3 данных можно заметить, что при уменьшении относительной погрешности между приближенным и точным значениями коэффициента С^ значение функционала J увеличилось. Примечательно, что оба алгоритма продолжили сходиться, несмотря на эту особенность, и значение коэффициента С^ было определено
с малой погрешностью. Однако во втором алгоритме возникновение данной особенности привело к нарушению монотонной сходимости приближенного решения к точному, в то время как в первом алгоритме монотонная сходимость сохранилась.
В третьем примере последовательное определение значений объемных теплоемкостей проводилось на основе данных изменения температуры от времени в двух внутренних точках. В действительности можно обойтись измерениями температуры, проведенными только в одной точке. Главное, чтобы за время наблюдений снижение температуры в наблюдательной скважине при замораживании породного массива было достаточным для определения неизвестных значений. В данной работе при определении объемной теплоемкости для зоны охлаждения рассматривалась точка, в которой произошло снижение температуры примерно на 2 К, а восстановление объемной теплоемкости для зоны льда выполнялось по точке, в которой температура снизилась примерно на 20 К, при этом снижение температуры после фазового перехода составило примерно 8 К. Тем не менее установление зависимости между точностью определения значений объемных теплоемкостей и изменением температуры в наблюдательных скважинах требует дополнительных исследований.
Заключение
Представлен подход к исследованию процесса формирования ледопородного ограждения на основе решения двумерной коэффициентной обратной задачи Стефана. Основной целью работы является определение объемной теплоемкости для зоны льда и зоны охлаждения по данным изменения температуры от времени в ограниченном числе внутренних точек расчетной области. В работе предложены алгоритмы численного решения поставленной задачи. Первый алгоритм разработан на основе метода сопряженных градиентов, второй - на основе метода наискорейшего спуска. Для реализации первого алгоритма были установлены вид задачи в приращениях температуры, вид сопряженной задачи, вид градиента функционала невязки.
Вследствие того, что значение объемной теплоемкости для зоны льда оказывает существенно меньшее влияние на распределение температуры в породном массиве по сравнению со значением объемной теплоемкости для зоны охлаждения, определение неизвестных теплофизических параметров происходит последовательно, начиная с объемной теплоемкости для зоны охлаждения.
Эффективность разработанных алгоритмов была проиллюстрирована путем решения модельных примеров, заключающихся в определении объемной теплоемкости либо для одной зоны, либо для обеих зон. В результате было установлено, что для всех рассмотренных случаев оба алгоритма позволяют определить значения объемных теплоемкостей с относительной погрешностью менее 1 %. При этом скорость сходимости второго алгоритма выше, а вычислительные затраты ниже, чем у первого.
Представленный подход и разработанные алгоритмы могут быть использованы для проектирования и уточнения исходных данных процесса формирования ледопородного ограждения при строительстве вертикальных шахтных стволов с использованием технологии искусственного замораживания породного массива.
Благодарности
Работа выполнена при финансовой поддержке РНФ (грант № 17-11-01204). Acknowledgments
The work has been carried out with the financial support from the Russian Science Foundation, RSF (Grant Nr. 17-11-01204).
Библиографический список
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.C0;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. Slota 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., Slota 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., Slota 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. - С. 119132. 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.
References
1. Trupak N.G. Zamorazhivanie gruntov v podzemnom stroitel'stve [Freezing of soils in underground construction]. Moscow, Nedra, 1974, 280 p.
2. Trupak N.G. Zamorazhivanie porod pri sooruzhenii vertikal'nykh stvolov shakht. [Freezing of rocks in the construction of mineshafts]. Moscow, Nedra, 1983, 170 p.
3. Andersland O.B., Ladanyi B. Frozen ground engineering. John Wiley & Sons, 2004, 352 p.
4. Vabishchevich P.N., Vasil'eva M.V., Pavlova N.V. Numerical modeling of thermal stabilization of filter ground. Mathematical Models and Computer Simulations, 2014, vol. 26, no. 9, pp. 111-125. DOI: 10.1134/S2070048215020106
5. Vabishchevich P.N., et al. Mathematical modeling of the artificial freezing of soils. Computational Technologies, 2014, vol. 19, no. 4, pp. 19-31.
6. Amosov P.V., Lukichev S.V., Nagovitsyn O.V. Vliianie poristosti porodnogo massiva i temperatury khla-donositelia na skorost' sozdaniia sploshnogo ledoporodnogo ograzhdeniia [Influence of rock massif porosity and coolant's temperature on velocity of solid ice wall creation]. Vestnik Kol'skogo nauchnogo tsentra RAN, 2016, vol. 4 (27), pp. 43-50.
7. Levin L.Iu., Zaitsev A., Semin M.A. Kontrol' teplovogo rezhima porodnogo massiva na osnove primene-niia optovolokonnykh tekhnologii monitoringa temperatur v skvazhinakh [Dynamics of ice wall under conditions of damaged freezing pipes when shaft sinking]. Gornoe ekho, 2016, no. 1, p. 35-37.
8. Samarskii A.A., Vabishchevich P.N. Vychislitel'naia teploperedacha [Computational heat transfer]. Moscow, Librokom, 2009, 782 p.
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, pp. 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, pp. 875-898. DOI: 10.1002/nme.1620300419
11. Krylov D.A., Sidnyaev N.I., Fedotov A.A. Mathematical modelling of temperature distribution. Mathematical Models and Computer Simulations, 2013, vol. 25, no. 7, pp. 3-27.
12. Vakulenko I.S., Nikolaev P.V. Analiz i perspektivy razvitiia sposoba iskusstvennogo zamorazhivaniia gornykh porod v podzemnom stroitel'stve [Analysis and outlook for development of artificial freezing of rocks in underground construction]. Gornyi informatsionno-analiticheskii biulleten' (nauchno-tekhnicheskii zhurnal), 2015, no. 3, pp. 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, pp. 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, pp. 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, pp. 4409-4423. DOI: 10.1016/j.ij heatmasstransfer.2006.10.019
16. Voller V.R. Enthalpy method for inverse Stefan problems. Numerical Heat Transfer, Part B Fundamentals, 1992, vol. 21, no. 1, pp. 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, pp. 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, pp. 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, pp. 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, pp. 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, pp. 657-684. DOI: 10.1016/j.jcp.2006.09.030
22. Nejad A.A. et al. Optimal operation of alloy material in solidification processes with inverse heat transfer. International Communications in Heat and Mass Transfer, 2010, vol. 37, no. 6, pp. 711-716. DOI: 10.1016/j.icheatmasstransfer.2010.03.002
23. Slota 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, pp. 155-176. DOI: 10.1080/10407790802605232
24. Hetmaniok E., Slota 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, pp. 7-14. DOI: 10.1016/j.icheatmasstransfer.2013.07.009
25. Hetmaniok E., Slota 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, pp. 347361. DOI: 10.1016/j.camwa.2014.12.013
26. Rashedi K. et al. Application of meshfree methods for solving the inverse one-dimensional Stefan problem. Engineering Analysis with Boundary Elements, 2014, vol. 40, pp. 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, pp 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, pp. 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, pp. 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. Gol'dman N.L. Uniqueness classes of solutions of two-phase coefficient inverse stefan problems. Dokla-dyMathematics, 2013, vol. 87, no. 2., pp. 205-210. DOI: 10.1134/S1064562413020269
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, pp. 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, pp. 599-615. DOI: 10.1080/17415977.2011.579605
36. Samarskii A.A., Vabishchevich P. N. Numerical Methods for Solving Inverse Problems of Mathematical Physics. Walter de Gruyter, 2007, 438 p.
37. Kolesnik S.A. A method for identification of nonlinear components of heat transfer tensor for anisotropic materials. Mathematical Models and Computer Simulations, 2014, vol. 26, no. 2, pp. 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, pp. 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.
