Проверка кода для численного моделирования тепловых процессов в пористой среде с учетом фазового перехода "лед – вода"

Цель исследования – проверить доработку кода, разрабатываемого для численного моделирования тепловых процессов в пористой среде при наличии фазового перехода "лёд – вода". Код создаётся для использования при решении следующих проблем:

• 1)    исследование теплового воздействия подземных атомных станций малой мощности (АСММ), размещаемых в труднодоступных районах Восточной Сибири на многолетнемёрзлые горные породы (ММГП);

• 2)    оценка тепловой безопасности ММГП на потенциальном объекте окончательной изоляции ОЯТ Билибинской АЭС.

• 2.    Алгоритм усовершенствования и результаты сравнения

Первая проблема решается в рамках НИР "Научно-техническое обоснование и разработка концепции создания заглублённых и подземных АСММ модульного типа для энергоснабжения горнопромышленных предприятий и населённых пунктов в труднодоступных регионах России" (научные руководители: академик РАН Н.Н. Мельников, профессор В.П. Конухин).

Вторая проблема исследуется в рамках дипломной работы студентки физико-энергетического факультета КФ ПетрГУ Е.В. Резец (научный руководитель: к.т.н., доцент П.В. Амосов).

В принципе, исследования теплового состояния ММГП можно было бы проводить на базе известных программных продуктов, например, PORFLOW (метод конечных разностей) или COMSOL (метод конечных элементов), которые позволяют решать подобные или достаточно близкие к ним задачи. Вместе с тем, иметь собственный "инструмент", который можно модернизировать под новые близкие проблемы (например, замёрзшие отвалы или хвостохранилища в условиях Крайнего Севера), несомненно лучше. Тем более, что собственный проверенный программный продукт можно использовать не только в научных исследованиях.

Модернизации подвергается одна из версий программы FFM, которая более 20 лет используется в Горном институте КНЦ РАН (разработчик: с.н.с. А.В. Наумов) при решении уравнения теплопроводности (метод конечных разностей, алгоритм Патанкара) ( Мельников и др ., 2001). Учёт фазового перехода выполнен по алгоритму, изложенному в русскоязычной литературе, например, в коллективной монографии сотрудников ИГД Севера им. Н.В. Черского СО РАН ( Курилко и др ., 2011).

Проверка программы осуществлялась на следующем модельном примере. Имеем замёрзшую (–20 °С) пористую (10 %) породу размером 1 X 1 м. Граничные условия следующие: на трёх границах (нижней и боковых) условия изоляции, на верхней – фиксированное значение температуры (100 °С).

Амосов П.В. Проверка кода для численного моделирования…

Сравнительному анализу подвергается центральное вертикальное распределение температуры на 24 часа процесса оттаивания.

Естественно, что первоначально описанная задача была решена с помощью указанных выше верифицированных кодов, и было достигнуто, на взгляд автора, хорошее совпадение результатов (см. таблицу). При этом параметры моделей, описывающих фазовый переход, приняты по умолчанию в соответствии с рекомендациями разработчиков.

Коротко опишем подход, используемый д.т.н. Курилко А.С. с коллегами ( Курилко и др. , 2011). В двухмерной постановке процесс распространения тепла в массиве горных пород с учётом фазовых переходов описывается следующей системой уравнений и условий:

[ C ( T) + LWp6 ( T- T * )]( дT / dt ) = д [ X ( T )( дT / дх )]/ дх + д [ X ( T)(дT / ду )]/ ду ,

• Ic 1 p i ,                  T < T ,

с ( T) = I

Ic 2 P 2 ,                 T >  T * ,

|X 1 ,                  T < T ,

X ( T) = I

|X 2 ,                    T >  T * ,

0 < x < H,        0 < y < H„        t > 0, xy где T - температура породы, К; T* - температура фазовых переходов влаги в породе, К; t и х,у -временная (с) и пространственные координаты (м); L - скрытая теплота плавления (замерзания) льда (воды), Дж/кг; W - влажность породы, доли единицы; p - плотность воды, кг/м3; c 1, p 1, X 1 (c2, p2, X2) -удельная теплоёмкость (Дж/(кг-К)), плотность (кг/м3) и коэффициент теплопроводности (Вт/(м-К)), соответственно для мёрзлых (талых) пород; 6(T-T) - 6-функция Дирака; Hx и Ну - границы области моделирования по осям координат, м.

Из множества существующих вариантов аппроксимации эффективной теплоёмкости широкое применение получили методы, в которых влияние 6 -функции Дирака распространяется только на две смежные точки пространственной сетки (в нашем случае попарно вдоль продольной и поперечной осям) и вариант, когда эффективная теплоёмкость непрерывна в точках T -AT и T + AT .

В соответствии с методом сглаживания 6 -функция Дирака приближённо заменяется 6 -образной функцией 6 ( T-T , AT) , которая отлична от нуля на интервале ( T -AT , T + AT) и удовлетворяет условию нормировки

T * +AT

/6 ( T- T * , AT) = 1.

T* –∆T

Записывается эффективная теплоёмкость [ C ( T)"] = C ( T) + LWp6 ( T- T * , AT) , которая далее интегрируется в пределах ( T-AT , T + AT) . В результате получаем следующее условие постоянства энтальпии на интервале сглаживания

T * + AT

/ [ C ( T )] dT = LWp + ( c 1 p 1 + c 2 p ₂) AT .

T*–∆T

Предлагается эффективную теплоёмкость записывать в следующем виде

Ic i p 1 ,                                                                              T < T * - AT ;

[ C ( T )] = I ( c 1 p 1 + c ₂ p ₂)/2 + ( c 1 p 1 - c₂p ₂)( T- T * )/2/ AT + LWp (1 -|T- T * | / AT )/ AT ,          |T- T * | <  AT ;

• Ic ₂ p ₂,                                                                            T> T * + AT .

Разрывность коэффициента теплопроводности устраняется посредством соединения точек ( X ₁, T^* - AT ) и ( X ₂, T^* + AT ) прямой линией

• IX 1 ,                                   T < T * - AT ;

[ X ( T ) ] = I ( X 1 + X 2 )/2 + ( X 2 - X 1 )( T - T * )/2/ AT , IT - T * I <  AT ; IX ₂,                                  T> T * + AT .

Таким образом, начальное уравнение теплопроводности при численной реализации заменяется уравнением

[ c ( T )]( дтi дt ) = д [[ x ( T ) ]( дт / дх )]/ дх + д [[ x ( T )]( дт / ду )]/ ду ,

0 < х < H x ,              0 <у < H y ,               t > 0.

В начальный момент времени задаётся распределение температуры:

T ( x , y ,0) = T 0 ,    0 ≤ x ≤ H x ,       0 ≤ y ≤ H y .

В задаче, выбранной для проверки, использованы следующие граничные условия:

на верхней границе –                T ( x , H_y ,0) = T_up , 0 ≤ x ≤ H_x (условие Дирихле);

на нижней и боковых границах –    λ ( ∂T / ∂y ) |_{y =0} = 0, 0 ≤ x ≤ H_x и

λ(∂T/∂x)|x=0,x=Hx = 0, 0 ≤ y ≤ Hy (условие Неймана, в частности, нулевой поток тепла).

Результаты численных экспериментов с использованием всех программных продуктов сведены в таблицу.

По мнению автора, полученные результаты позволяют с определённым оптимизмом смотреть на возможное использование в будущих расчётах собственного кода. Действительно, согласие неплохое, хотя некоторое расхождение в цифрах наблюдается. Например, собственный код выдаёт заниженные положительные температуры (относительная ошибка по модулю 7,8 %), но в области фазового перехода собственные расчётные значения попадают в середину интервала верифицированных кодов.

В качестве примера потенциального применения программы может быть задача по определению скорости протаивания замёрзшей горной породы при наступлении оттепели. В плоской геометрии на модели размером 30X20 м прослежено продвижение "фронта" нулевой температуры вглубь массива при задании на поверхности условия 3-го рода. Выполненная оценка позволила получить значение скорости протаивания на уровне 3-4 см/сутки за время порядка 60 суток. Результат достаточно физичен и отвечает набору исходных параметров.

Таблица. Пространственное распределение температуры на 24 часа процесса оттаивания, °С

Координата, м	Код COMSOL	Код PORFLOW	Собственный код
0,0	–17,2	–17,1	–17,6
0,1	–16,7	–16,6	–17,1
0,2	–14,7	–15,0	–15,6
0,3	–12,3	–12,1	–12,7
0,4	–7,9	–7,7	–8,1
0,5	–1,8	–0,9	–1,3
0,6	10,7	10,9	9,7
0,7	28,2	28,2	26,0
0,8	49,6	49,6	47,2
0,9	74,1	74,0	72,5
1,0	100,0	100,0	100,0

3. Заключение

Достигнутая точность расчётных значений температуры в сравнении с результатами, полученными с использованием верифицированных программ POFRLOW и COMSOL, позволяет автору рекомендовать созданную программу для численного моделирования процессов переноса тепла при наличии фазового перехода "лёд – вода".