Приближенный метод решения задач тепломассопереноса при замерзании влаги в горных породах криолитозоны

УДК 622.45:536.421

В.И. Попов, А.С. Курилко

ПРИБЛИЖЕННЫЙ МЕТОД РЕШЕНИЯ ЗАДАЧ ТЕПЛОМАССОПЕРЕНОСА ПРИ ЗАМЕРЗАНИИ ВЛАГИ В ГОРНЫХ ПОРОДАХ КРИОЛИТОЗОНЫ

Аннотация. Предложен метод решения задач тепломассопереноса с фазовым превращением (замерзание-оттаивание), основанный на процедуре расщепления по физическим процессам исходной задачи тепломассопереноса с фазовым превращением. В качестве таковых выбраны процессы диффузионного переноса тепла, влаги и растворимых солей между узлами сеточного разбиения и процессы перераспределения в узловых областях, рассматриваемые как изолированные системы. Разработанный вычислительный алгоритм позволяет, при интегрировании в пакеты прикладных программ решения нелинейных уравнений теплопроводности, рассчитывать сложные процессы с фазовым превращением. Его независимость от пространственных параметров (на этапе фазовых превращений) позволяет в рамках единого методического подхода применять как для процессов с узким температурным спектром зоны фазового перехода (задача типа Стефана), так и с широкой зоной в соответствии с видом уравнения фазового равновесия в том числе многофронтовых и многомерных. Тестирование на автомодельных и численных решениях показывает удовлетворительное (с точностью до 2%) соответствие. Это дает хорошие перспективы дальнейшего практического использования.

Ключевые слова: тепломассоперенос, фазовые превращения, расщепление по физическим процессам, метод решения .

Введение

Одной из первых математических моделей описания процесса промерзания влажного дисперсного материала, является задача Стефана, отличающаяся тем, что уравнение состояния поровой влаги формулируется как одномоментное изменение фазового состояния при достижении определенной температуры. Данная формулировка нашла широкое распространение при исследовании процессов замерзания — оттаивания в виде различных постановок задач типа Стефана и их решения методом ловли фронта в узел сетки [1], выпрямления фронтов при решении многофронтовых задач [2].

В работе [3] предложена отличная от постановки Стефана математическая

DOI: 10.25018/0236-1493-2018-12-0-57-64

формулировка задачи о температурном режиме грунта, соответствующая промерзанию в спектре температур. Если в задаче Стефана рассматривались две зоны — талая и мерзлая то в постановке [3] зона фазовых переходов помимо льда содержит жидкую фазу — равновесное по отношению к температуре количество незамерзшей воды. В случае существования двухфазной зоны, функция мощности теплового источника фазовых превращений заменяется на распределенную в некотором температурном интервале плотность функции тепловыделения при фазовых превращениях А = А + Lд'^/дг, где IV — содержание незамерзшей воды. Учет миграции влаги при промерзании и таянии грунтов проводится с использованием системы

ISSN 0236-1493. Горный информационно-аналитический бюллетень. 2018. № 12. С. 57-64. © В.И. Попов, А.С. Курилко. 2018.

уравнении совместного тепломассопе-реноса [4]. И хотя в оригинальном виде данная система уравнении содержит ряд трудноопределимых параметров, ее модификации получили широкое распространение [5—7]. При этом основной вычислительной схемой расчетов стал метод эффективной теплоемкости [8], который обобщается и на более сложные функциональные зависимости для источника, но требующий существенных итерационных процедур [5].

Для комплексных задач тепломассо-переноса включающих значительное число переносимых субстанций положение несколько осложняется [5]. Компромисс заключается в развитии подхода, более эффективно учитывающего связи, накладываемые на параметры состояния уравнением фазового равновесия. Этими параметрами, рассматриваемыми в статье, являются температура, содержание незамерзшей воды (уравнение фазового равновесия) и концентрация растворенного компонента.

Интенсивная миграция влаги в промерзающих грунтах объясняется резким падением потенциала влаги при образовании кристаллов льда [9].

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

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

а

дг д'М

д_ дх

дг дх

^-и - ад

дх

дг

где потоки ] равны

1 дт

+ Сш РЛ

дх

^ = -к

дШ дх

^ = -ШЮс дС + С]„ дх

Здесь член /Р = -дт/дЬ представляет мощность источника фазовых превращений (вода-лед); комбинация членов к£/Р определяет захват солей ледовой фазой; параметр кг — коэффициент захвата солей (ионов); X, К, йс — коэффициенты теплопроводности, диффузии влаги, диффузии солей; с = сз + м*^ + w*с. теплоемкость среды; индексы м, ^ s, 0 — относятся соответственно к воде, льду, минеральному скелету и начальным данным; р — плотность вмещающей среды.

Система уравнений (1) дополняется уравнением фазового состояния 7Р = 7Р (№, С) — в виде аналитического или интерполяционного выражения.

Постановка задачи

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

Процедура численного решения

Для иллюстрации метода решения в качестве расщепленных физических процессов [10] выбраны: а) диффузионное выравнивание полей температуры, влажности и концентрации, между расчетными узлами сеточного разбиения и б) локальный фазовый переход и вызванное им перераспределение компонентов между фазами в расчетных уз-

лах. Таким образом, задача (1) расщепляется на внешний тепломассообмен между узлами сеточного разбиения и локальный — в конкретном узле. Поэтому решение исходной системы (1) можно представить в виде последовательного решения задач:

а) диффузии

а(ср7) дг

dW _

~дГ _

d(WC) дг

—j

дх q

-a j

Л W

дх

дх

Jr

б) фазового перехода в области расчетного узла

д

дг

dW ~дг d(WC) дг

= LpIF

= -kCL

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

Решение системы уравнений (3) следует известным способам численного расчета процессов диффузии (параболических систем уравнений). Отличительной особенностью предлагаемого в статье метода является способ решения системы (4). Так как 1Г =-дт / дг мощность стока влаги при переходе (вода-

лед), то подставляя его во второе и третье уравнения системы (4) получим для влажности и концентрации соли на каждом шаге по времени (рис. 1):

W, = Щ + Am; C3 = C2 \W/W

Для определения стока влаги Ат приходящейся на фазовый переход на каждом временном шаге используем уравнение теплового баланса (первое уравнение из системы (4)), а также условие, что температура в процессе фазового перехода определяется диаграммой состояния Т3 = ТР = ТР(W,С):

APsTs - AP2T = -LPAm TF = TF

При этом, переход (Т1, С1 ^ Т2, М2, С2) определяется решением разностных аналогов системы (3), и температура Т2, отражая энергетический баланс процесса теплопроводности носит фиктивный характер (не учитывает фазовый переход). Переход Т2^Тр определяется решением системы (6).

Рассматривая путь Т1^Тр на диаграмме рис. 1, полагаем с точностью до

Рис. 1. Диаграмма расчета расщепленных процессов тепломассообмена. Вектор 1-2 — процессы диффузии тепла и массы; вектор 2-3 — процесс фазового превращения; вектор 1-3 — путь по траектории фазового превращения, определяемый уравнением состояния Fig. 1. Split heat and mass transfer processes chart. Vector 1-2—heat and mass diffusion; 2-3— phase transformation; 1-3—phase transformation path from the equation of state

1 -k

бесконечно малых величин первого порядка:

дУтAm

T = T

После подстановки (7) в систему (6) получим выражение для интенсивности стока влаги в узле при промерзании, ограничиваясь (без существенной потери точности) членами со степенями Ат не выше первой:

с

Am = -

с

/дW (( )Т Определив тем самым интенсивность стока влаги при промерзании на временном шаге, в узле расчетной сетки.

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

1. Решается система уравнений (3) диффузионного переноса тепла, влаги и солей без фазовых превращений.

2. По формуле (8) рассчитывается количество замерзшей влаги.

3. По формулам (5) и (6) рассчитываются равновесные составы фаз и компонентов, а затем температура фазового превращения ТР - ТР С).

4. Осуществляется обновление исходных значений нижнего временного уровня и переход на следующий временной шаг (пункт 1).

На рис. 2—3 приведены результаты расчета перераспределения температуры и влагосодержания по системам уравнений (3)—(4).

Начальные и граничные условия для системы (3) имеют вид:

начальные условия: Т0 = 4 °С; = = 0,2; С0 = 1/58 моль/л;

на границе х = 0 ставится условие конвективного теплообмена

J =*

где

Tcp (t ) = A,

Bcp sin

2%t

V ^per

л Q. >1 H П5 Q. (D

X

0) X CD X

n (D Я

iS

Q. S 0)

m П5 Q.

m

w

LQ

1,4

1,0

0,6

0,2

-0,2

1 It

VI ill UF'

/

2 ' / /

0

0,125

0,250

0,375 Глубина массива, м

Рис. 2. Распределение полей температуры

л Q.

>1

н го

Q. CD 1=

5 :

0)

I-

о;

CD X ГО I п CD

15

О. S

CD

т ГО Q. т CD Ш

1,4

1,0

0,6

0,2

-0,2

1 1

I1 I

u/L' 'Л'1мл А

1 \ 1 \ 3 \ / \:

11 и I

«   | 2 / !!

* \ V

г

О

0,175

0,350

0,525 Глубина массива, м

Рис. 3. Распределение полей температуры

Fig. 3. Distribution of temperature field,

причем Ас = -8

Ср2

'С, Вср = 16 °С; а =

15 Вт/(м2град).

На границе х = I Т(хД) = Т0 и Т(х^) = 0 для второго варианта.

Для концентрации и влагосодержания на границах области принято отсутствие потоков (^ = 0, JC = 0). Расчет приведен для граничных условий воздействия знакопеременной температуры на поверхности х = 0 с периодом {ег = 2 ч, I = 0,7 м (рис. 2) и {рег = 24 ч, I =еГ1,0 м (рис. 3). Длительность процесса воздействия составляет 96 и 200 ч соответственно.

В качестве оценки эффективности того или иного метода допускается сравнение полученных результатов с точным аналитическим решением или с данными решений других авторов. Для одномерного случая существует автомодельное решение задачи Стефана, которое представляет собой трансцендентное уравнение. Результаты численных расчетов были сравнены с точным решением задачи Стефана при следующих исходных данных: первоначальная тем-

пература массива горных пород равнялась -5 °С. Коэффициенты теплопроводности соответственно для мерзлых и талых пород составили 2,45 и 1,83 Вт/(м • К). Влажность пород равнялась 34%, плотность пород в талом и мерзлом состояниях одинакова (1745 кг/м3); удельная теплоемкость мерзлых горных пород составляла 1224 Дж/(кг • К), талых — 1800 Дж/(кг • К). На границе области граничное условие первого рода — Т = = +10 °С. Результаты сравнительных расчетов приведены в таблице.

Как видно из таблицы, точность решения полученная по данному методу вполне удовлетворительная, относительная ошибка не превышает 3,1% за период времени продолжительностью 2000 ч.

Заключение

Проведенные контрольные расчеты показывают вполне удовлетворительное описание процессов распределения тепла, влагосодержания. Имеется возможность непосредственно рассчитывать

Сравнение расчетов по данной модели с результатами автомодельного решения задачи Стефана

Comparison of the proposed model data and the results of the similarity solution of a Stefan problem (initial temperature is non-zero). The data for the comparison are borrowed from [13] p. 41

Время, ч Глубина протаивания (точное решение), м Алгоритм А. Самарского и Б.Д. Моисеенко [8] с автоматическим выбором шага Алгоритм Дж. Цяо [12] Предложенный алгоритм

1* 2* 3* 1* 2* 3* 1* 2* 3*

100 0,2957 0,3263 0,0306 10,3 0,3133 0,0176 6,0 0,305 0,0093 3,1

200 0,4181 0,4417 0,0236 5,6 0,4476 0,0295 7,1 0,422 0,0039 0,9

300 0,5121 0,5303 0,0182 3,6 0,5481 0,036 7,0 0,512 1E-04 0,0

400 0,5913 0,6128 0,0215 3,6 0,631 0,0397 6,7 0,587 0,0043 0,7

500 0,6611 0,6885 0,0274 4,1 0,7024 0,0413 6,2 0,654 0,0071 1,1

600 0,7242 0,7506 0,0264 3,6 0,7717 0,0475 6,6 0,714 0,0102 1,4

700 0,7822 0,8031 0,0209 2,7 0,8334 0,0512 6,5 0,77 0,0122 1,6

800 0,8362 0,863 0,0268 3,2 0,889 0,0528 6,3 0,822 0,0142 1,7

900 0,887 0,9079 0,0209 2,4 0,9446 0,0576 6,5 0,871 0,016 1,8

1000 0,9349 0,9613 0,0264 2,8 0,9934 0,0585 6,3 0,917 0,0179 1,9

1100 0,9806 0,9999 0,0193 2,0 1,0433 0,0627 6,4 0,962 0,0186 1,9

1200 1,0242 1,0477 0,0235 2,3 1,0886 0,0644 6,3 1 0,0242 2,4

1300 1,066 1,092 0,026 2,4 1,1335 0,0675 6,3 1,05 0,016 1,5

1400 1,1062 1,1293 0,0231 2,1 1,1757 0,0695 6,3 1,09 0,0162 1,5

1500 1,1451 1,1735 0,0284 2,5 1,2167 0,0716 6,3 1,12 0,0251 2,2

1600 1,1826 1,2049 0,0223 1,9 1,2572 0,0746 6,3 1,16 0,0226 1,9

1700 1,219 1,2456 0,0266 2,2 1,2941 0,0751 6,2 1,2 0,019 1,6

1800 1,2544 1,2842 0,0298 2,4 1,3327 0,0783 6,2 1,23 0,0244 1,9

1900 1,2887 1,3148 0,0261 2,0 1,3693 0,0806 6,3 1,27 0,0187 1,5

2000 1,3222 1,3526 0,0304 2,3 1,4033 0,0811 6,1 1,33 0,0078 0,6

* 1 — глубина протаивания, м; 2 — абсолютная ошибка, м; 3 — относительная ошибка, %.

многозонные задачи, так как фазовое равновесие определяется в каждом расчетном узле. Налицо существенное снижение трудоемкости вычислений, так как задача сводится к решению трех нелинейных уравнений параболического типа без фазового перехода и расчету состава по формулам типа (5)—(6).

Разработанный вычислительный алгоритм позволяет, при интегрировании в

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

так и с широкой зоной в соответствии с видом уравнения фазового равновесия в том числе многофронтовых и многомерных. Тестирование на автомодельных и численных решениях показывает

СПИСОК ЛИТЕРАТУРЫ

удовлетворительное (с точностью до 2%) соответствие.

Это может обеспечить хорошие перспективы дальнейшего практического использования.

1. Будак Б. М., Васильев Ф. П. Егорова А. Т. Об одном варианте неявной разностной схемы с ловлей фронта в узел сетки для решения задач типа Стефана / Вычислительные методы и программирование. Вып. 4. — М.: Изд-во МГУ, 1967. — С. 231—241.

2. Будак Б. М., Гольдман Н.Л., Успенский А. Б. Разностные схемы с выпрямлением фронтов для решения многофронтовых задач типа Стефана // Доклады АН СССР. — 1966. — Т. 167. — № 4. — С. 735—738.

3. Колесников А. Г. К изменению математической формулировки задачи о промерзании грунта // Доклады АН СССР. — 1952. — Т. 82. — № 6. — С. 889—892.

4. Лыков А. В. Явления переноса в капиллярно-пористых телах. — М.: Изд-во техн.-теорет. лит-ры, 1954.

5. Пермяков П. П., Романов П. Г. Тепло- и солеперенос в мерзлых ненасыщенных грунтах. — Якутск: Изд-во ЯФ СО РАН, 2000.

6. Taylor G.S., Luthin J. N. A model for coupled heat moisture transfer during soil freezing // Canad. Geotechnical J. 1978. V. 15. P. 548—555.

7. Jame Y. W., Norum D. J. Heat and mass transfer in freezing unsaturated porous medium // Water Resour. Rec. 1980. V. 16. No 4. P. 811—819.

8. Самарский А. А., Моисеенко Б. Д. Экономичная схема сквозного счета для многомерной задачи Стефана // Журнал вычислительной математики и математической физики. — 1965. — Т. 5. — № 5. — С. 816—827.

9. Комаров И.А. Термодинамика и тепломассообмен в дисперсных мерзлых породах. — М.: Научный мир, 2003.

10. Марчук Г. И. Методы вычислительной математики. — М.: Наука, 1977.

11. Карпов И. К. Физико-химическое моделирование на ЭВМ в геохимии. — Новосибирск: Наука, 1981.

12. Hsiao J. S. An efficient algoritm for finite-difference analyses of heat transfer with melting and solidification // Numer. Heat Transfer. 1985. V.8, No 6. P. 653—666.

13. Хохолов Ю.А., Соловьев Д. Е. Математическое моделирование тепловых процессов в горных выработках шахт и рудников Севера. — Новосибирск: Академическое изд-во «Гео», 2013. — 185 с. ЕШ

КОРОТКО ОБ АВТОРАХ

Попов Владимир Иванович1 — кандидат технических наук, старший научный сотрудник, E-mail: popov.gtf@mail.ru, Курилко Александр Сардокович1 — доктор технических наук, профессор, заместитель директора по научной работе, 1 Институт горного дела Севера им. Н.В. Черского СО РАН.

ISSN 0236-1493. Gornyy informatsionno-analiticheskiy byulleten'. 2018. No. 12, pp. 57-64.

Approximate method to solve problems of heat and mass transfer under water freezing in permafrost rocks

Popov V.I.1, Candidate of Technical Sciences, Senior Researcher, e-mail: popov.gtf@mail.ru, Kurilko A.S.1, Doctor of Technical Sciences, Professor, Deputy Director for Science, 1 Chersky Mining Institute of the North, Siberian Branch, Russian Academy of Sciences, 677000, Yakutsk, Russia.

Abstract. The problems of predicting heat and humidity conditions in various mine structures, as well as taking into account induced pollution of frozen soil and pollutant migration in permafrost strata need new mathematical models. It is required that these models take into account to a fuller extent features of physical freezing and thawing processes in pore water. It is also of the current concern to develop new software tools and methods for the numerical solution of model problems on heat-moisture-salt transfer in multi-dimensional domains with regard to phase ice-solution transition. This article proposes a method to solve the problems of heat and and mass transfer with phase transformation (freezing-thawing). The method involves splitting of the initial problem of heat and mass transfer with phase transformation with respect to physical processes. Such processes are assumed to be diffusion transfer of heat, moisture and soluble salt between points of meshing, as well as the processes of redistribution in the point neighborhoods which are assumed to be isolated systems. The developed algorithm allows calculating complex processes including phase transformations while integrating solutions of nonlinear equations of thermal conduction in application packages. Irrespective of space parameters (at the stage of phase transformations), the algorithm makes it possible to use equations of phase equilibrium both for the processes with narrow and wide temperature range of the phase transfer zone (Stefan problem type), including multifront and multidimensional processes, in the framework of the common methodical approach. Tests of the similarity and numerical solutions show a satisfactory agreement (accurate to 2%). This offers good prospects for further applications.

Key words: Heat and mass transfer, phase transformations, splitting with respect to physical processes, solution method.

DOI: 10.25018/0236-1493-2018-12-0-57-64

REFERENCES

1. Budak B. M., Vasil'ev F. P. Egorova A. T. Ob odnom variante neyavnoy raznostnoy skhemy s lovley fronta v uzel setki dlya resheniya zadach tipa Stefana [A variant of implicit difference scheme with front capturing at mesh point for solving Stefan-type problems], Vychislitel'nye metody i programmirovanie. Issue 4], Moscow, Izd-vo MGU, 1967, pp. 231-241.

2. Budak B. M., Gol'dman N. L., Uspenskiy A. B. Raznostnye skhemy s vypryamleniem frontov dlya resheniya mnogofrontovykh zadach tipa Stefana [Difference schemes with strengthened fronts for solving multifront Stefan problems], Doklady AN SSSR. 1966. Vol. 167, no 4, pp. 735-738. [In Russ].

3. Kolesnikov A. G. K izmeneniyu matematicheskoy formulirovki zadachi o promerzanii grunta [Changing mathematical formulation of problems on soil freezing], Doklady AN SSSR. 1952. Vol. 82, no 6, pp. 889— 892. [In Russ].

4. Lykov A. V. Yavleniya perenosa vkapillyarno-poristykh telakh [Transfer phenomena in capillary-porous bodies], Moscow, Izd-vo tekhn.-teoret. lit-ry, 1954.

5. Permyakov P. P., Romanov P. G. Teplo- i soleperenos v merzlykh nenasyshchennykh gruntakh [Heat and salt transfer in unsaturated frozen soil], Yakutsk, Izd-vo YAF SO RAN, 2000.

6. Taylor G. S., Luthin J. N. A model for coupled heat moisture transfer during soil freezing. Canad. Geo-technical J. 1978. V. 15. P. 548—555.

7. Jame Y. W., Norum D. J. Heat and mass transfer in freezing unsaturated porous medium. Water Resour. Rec. 1980. V. 16. No 4. P. 811—819.

8. Samarskiy A. A., Moiseenko B. D. Ekonomichnaya skhema skvoznogo scheta dlya mnogomernoy zadachi Stefana [Economic end-to-end calculation scheme for multidimensional Stefan problem]. Zhurnal vychislitel'noy matematiki i matematicheskoy fiziki. 1965. vol. 5, no 5, pp. 816—827. [In Russ].

9. Komarov I. A. Termodinamika i teplomassoobmen v dispersnykh merzlykh porodakh [Thermodynamics and heat and mass exchange in dispersed frozen rocks], Moscow, Nauchnyy mir, 2003.

10. Marchuk G. I. Metody vychislitel'noy matematiki [Methods of computational mathematics], Moscow, Nauka, 1977.

11. Karpov I. K. Fiziko-khimicheskoe modelirovanie na EVM vgeokhimii [Computer-aided physicochemi-cal modeling in geochemistry], Novosibirsk, Nauka, 1981.

12. Hsiao J. S. An efficient algoritm for finite-difference analyses of heat transfer with melting and solidification. Numer. Heat Transfer. 1985. V.8, No 6. P. 653—666.

13. Khokholov Yu. A., Solov'ev D. E. Matematicheskoe modelirovanie teplovykh protsessov vgornykh vy-rabotkakh shakht i rudnikov Severa [Mathematical modeling of thermal processes in mines in the North], Novosibirsk, Akademicheskoe izd-vo «Geo», 2013, 185 p.