Математическое моделирование процесса промерзания грунта с учетом движения влаги в талой зоне Текст научной статьи по специальности «Физика»

УДК 532.546

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ПРОЦЕССА ПРОМЕРЗАНИЯ ГРУНТА С УЧЕТОМ ДВИЖЕНИЯ ВЛАГИ В ТАЛОЙ ЗОНЕ

А. Р. Павлов, М. В, Матвеева

Математические модели теиломассоиереноса в промерзающих-протаивающих дисперсных средах в зависимости от степени их влагона-сыщенности имеют различные виды. Это связано с содержанием свободной воды. Как установлено многими исследователями, фазовые превращения свободной воды происходят на границе промерзания, а связанная — кристаллизуется по мере понижения температуры.

1. Постановка задачи. В одномерном случае температурная задача, когда процесс промерзания идет с левого конца стержня 0 < х ^ I, состоит из следующих уравнений:

»<*<««•111

дТ д ( дТ \ с2ж = ^(а2-), т<х<1,Т>Т.. (2)

На границе раздела фаз выполняются условия

Т = Т„ х = ф), (3)

= х = №■ (4)

д X д X ЫС

х

Ньютона:

Т

А 1— = а(Т-Тс), х = 0, (5)

х

© 2009 Павлов А. Р., Матвеева М. В.

а на правом конце х = 1 — условие отсутствия теплового потока:

дТ

- = 0, х = 1. (6)

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

Т(х,0) = ^(х), х < I. (7)

Влажностная задача строится в двух различных видах в зависимости от принятого предположения о миграции влаги в процессе промерзания.

В случае учета миграции влаги как в талой, так и мерзлой зонах области уравнения для влажности пишутся в каждой зоне:

д~Ш д ( д\¥л \

-9Г = ^ ' °<*<т,Т<Т., (8)

дШ д ( дШ \

= т<*<1,т>т.. (9)

Граничные условия формулируются в зависимости от их фазового состояния. В талой области задаются условия дШ

к2-— = а*(№-ф), х = 0 ,Т>Т*, (10)

дх

дШ

— = 0, х = 1, Т>п. (11)

дх

Для мерзлого состояния указанные условия принимают вид , дШ1

дх дШ1

= а»(Ш - ф), х = 0, Т < Т», (12)

= 0, х = 1,Т < Т». (13)

дх

Начальное распределение влаги во всей области задано так:

Ш(х,0) = ^»(х), х < I. (14)

Условие на границе фазового перехода следует из подсчета баланса массы аналогично условию Стефана:

к1д^-к2™ = ш2(%)§, х = т- (15)

д х д х

Здесь введены следующие обозначения: Т(х, 1),Т* — температура точки области с координатой х в момент времени С и температура фазового перехода соответственно; с\, с2 (А1, А2) — объемные теплоемкости (коэффициенты теплопроводности) мерзлой и талой зон; Тс — температура среды; Ь — теплота кристаллизации 1м3 воды; Ж(х,-Ь) — весовая влажность талой зоны; Ж — влажность по жидкой фазе в мерзлой зоне; — весовая влажность по твердой фазе (льду); ]¥ = —

суммарная влажность в мерзлой зоне; к\, к2 — коэффициенты диффузии влаги соответственно в мерзлой и талой зонах; а* — коэффициент влагоотдачи; ф(х) — равновесная влажность, ниже которой нет влаго-обмена с окружающей средой.

Анализ приведенной математической модели (1)^(5) показывает, что температурная задача может быть определена полностью, если известно Ж включая ее значение при Т = Т*, которое входит в условие (4). Указанная величина ^г(Т) определяется с помощью решения влажностной задачи, которая, в свою очередь, требует знания температурного поля. Такая связь между искомыми величинами показывает, что задача тепломассопереноса уравнениями (1)-(15) не определена полностью. Для получения недостающего условия исходим из следующих известных экспериментальных результатов. Во-первых, температурное поле в окрестности границы фазового перехода изменяется медленно под влиянием выделяющейся в этой области скрытой теплоты кристаллизации воды. Во-вторых, температура как функция про-

х

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

х

ному закону. Это позволяет воспользоваться линейной интерполяцией значений Т(х, С) в окрестности границы промерзания, выбрав два значения, между которыми находится температура фазового перехода Т*. Этот факт дает искомое дополнительное условие. Пусть температура фазового перехода Т* в момент времени С находится между ее значе-

ниями T и T2, принимаемыми соответственно в точках x = a и x = Ь. Тогда функция T(x,t) в этом промежутке представляется в виде

Т = 1\ + ^^(Т2-Т1), a < х < Ь. (16)

Ь — a

В случае, когда миграция влаги происходит в одной талой зоне, из предыдущей задачи сохраняются уравнения (9)—(11). Когда талая зона занимает область £(t) ^ x ^ l, возникает необходимость задания граничного условия на подвижной границе промерзания. По этой проблеме нет единого мнения среди исследователей [1-6]. Необходимо заметить, что величина W* зависит от наличия и отсутствия свободной воды. В указанных работах в основном рассмотрен случай отсутствия свободной воды. Следовательно, определение этой величины должно быть различным в зависимости от величины влажности материала.

При отсутствии миграции влаги в мерзлой зоне (нет потока влаги) из уравнения (8) следует равенство

dt '

что означает постоянство суммарной влажности в процессе промерзания. На границе фазового перехода эта величина представляет сумму W(Tt) = Wi(T*) + W2(Tt). С другой стороны, эта сумма равняется значению влажности на границе промерзания со стороны талой зоны Wi(T*) + W2(Tt) = W„ т. е. W* = const. Таким образом, из предположения об отсутствии миграции влаги в мерзлой зоне следует постоянство величины W*.

Из всего сказанного следует, что влажностная задача в случае отсутствия миграции влаги в мерзлой зоне состоит из следующих соотношений:

dW d ( <9W\

-Ж = 1Гх(^), ф)<х<1, Т>П. (17)

В талом состоянии всей области 0 ^ x ^ l функция W удовлетворяет условиям

dW

k2— = a*(W-il>), х = 0 ,Т>Т*, (18)

dx

дШ

— = 0, х = 1, Т>п. (19)

Когда талая зона занимает область ^ х ^ I, граничные условия задаются в виде

Ш = Ш* = Шг(Т*) + Ш2(Т*), х = £(г), (20)

дШ

— =0, х = 1. (21)

Начальное условие задается во всей области:

Ш(х,0) = р*(х), х < I. (22)

2. Разностная схема для задач тепломассопереноса. Сформулированные выше задачи тепломассопереноса: первая, представленная уравнениями (1)—(16), и вторая, состоящая из уравнений (1)-(7), (16)—(22), решаются численным методом. Особенность первой задачи состоит в разрывности ее решения на границе фазового перехода. Ее численная реализация рассмотрена нами в предыдущих работах (см., например, [7]). В настоящей работе изложен метод численного решения второй задачи. Особенность данной задачи заключается в том, что, во-первых, температурная и влажностная задачи решаются в различных областях; во-вторых, влажностная задача задана в области, изменяющейся в процессе промерзания.

Для построения однородной разностной схемы температурную задачу обычно записывают в виде задачи теплопроводности для одного уравнения, заданного во всей области, т. е. уравнения (1), (2), (4) заменяются одним уравнением

где

с _ I с1Р1, Т<Т*, а-!А1' Т<Т*' I С202, Т > Т*, \ А2, Т > Т*,

а функция Ш2 в отличие от соответствующей функции уравнения (1)

Т Т*

правую часть (4).

Граничные и начальные условия (5)—(7) сохраняются и для этой задачи.

Доказательство справедливости такой замены проводится с использованием методики, развитой в работах [8-10].

Для численного решения задачи тепломассопереноса, состоящего из уравнений (23), (5)-(7), (17)-(22), строится на плоскости хоЬ (для простоты записи) равномерная разностная сетка по каждой переменной с шагами Н по оси х и т по Р х^ = гН, г = 0,1, 2,..., Щ ЫН = I;

Ь = Зт, 3 = 0,1,2,....

Во внутренних узлах сетки уравнения (17), (23) аппроксимируются следующими разностными уравнениями:

ОТг= (АТЯ)Х + ЬШ2 (24)

Ш = (к2Шх)х, (25)

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

То^Т+^ь (26)

Тм = к2Тм — + щ, (27)

коэффициенты которых имеют следующие выражения:

А0.5 ^(С0Т0 + Ь(Ш2-Ш2)о) + каТс

«1 = ТТТТ2-:-ГТ; = —

(С0£ + Ао.5 + аК)' (Со^ + Ао.й + а/»)

Адг-0.5 £(СкТк + Ь(№-2-^2)лг)

2т

К2 = Т-—1У?-:-г, щ =

( \ ' & / \ * (Сдг 27 + Адг-о.б) (С*ДГ27 + Ап-О.б)

Граничные условия (18), (19) аппроксимируются аналогично, когда обе границы х = 0 и х = I находятся в талом состоянии:

Ш0 = к3Ш1+"з, (28)

Шм = к4Шм-1+ Щ, (29)

_ №2)0.5 _ +

К3~ + (к2)0.5 + а*Ъ) ' 1/3 ~ + (к2)0.5 + а*Ъ) '

_ №2) N—0.5 _ ^ N

^ ~ (# + (*2)лг-0.б) ' ^ _ (# + (ЫлГ-О.б) ' В том случае, когда граница промерзания находится внутри области, имеем условие (20). Но подвижная граница х = £(£) может совпасть с одним из узлов сетки, может и находиться между двумя узлами. Пусть граница х = £(£) находится между узлами хда и хда+1- Тогда один из ближайших узлов будем считать граничным узлом и значения функции Ш на этом узле принимаем равными Ш*, т. е. производим простой снос значения функции в точке х = £(£) на ближайший узел сетки. Для повышения точности задания граничного условия можно воспользоваться способом Коллатца [11]. Таким образом, граничное условие согласно (20) можно написать в виде

Шт = Ш* = Шг(П) + Ш2(П). (30)

Начальные условия для разностной задачи получаются из условий (7) и (14):

Т(х<,0) = фн), Ш{х^) = <р* (х<). (31)

Присоединяя к уравнениям (24)-(31) разностный аналог уравнения (16), получаем полную систему уравнений, решение которой дает значения неизвестных функций в узлах сетки.

3. Примеры численных расчетов. Для проверки работоспособности построенного алгоритма выполнены расчеты процесса промерзания грунта при следующих входных данных задачи (для суглинка):

Сск = 1.127 кДж/кг • К, Св = 4.2 кДж/кг • К, Сл = 2.1 кДж/кг • К,

кг кг кг

рск = 1522—, рв = 1000—, рл = 916—,

мл мл мл

Ь = 51748 кДж/м3, а = 83.5 кДж/м2^К, а* = 840 м/ч, Ь2м.

Коэффициенты теплопроводности мерзлого и талого грунтов (суглинка) приняты согласно работе [12]:

Ам = 4.19(0.141 + 2.169 - 1.0492) кДж/м • ч • К, К = 4.19(0.183 + 1.О69 - 0.20992) кДж/м • ч • К, где 9 — степень влагонасыщенностн. Коэффициент диффузии влаги К = 0.378 6 • 1 О"5 ехр( 16.46Ш) м2/ч,

Тс= -20°С, Т* = 0, ф = 0.015, у=\0°С, у* = 0.16951.

В результате расчетов получены распределения температуры во всей области, влажности в талой зоне и суммарной влажности в мерзлой области.

Для примера распределения указанных величин при степени влагонасыщенностн 9 = 0.6 в момент времени £ = 200 ч приведены в табл. 1.

Таблица 1

Распределение температуры и влажности в процессе промерзания по глубине грунта

ж, м т°, С IV,%

0.0 -18.680 15.54

0.2 -11.756 15.54

0.4 -5.472 15.54

0.6 -0.236 15.54

0.8 3.681 16.956

1.0 6.346 16.951

1.2 8.015 16.951

1.4 8.981 16.951

1.6 9.496 16.951

1.8 9.737 16.951

2.0 9.806 16.951

Таблица 2

Распределение температуры и влажности, полученное решениями первой и второй задач тепломассопереноса

ж, м Т - 1,°,С Т — 2,°,С IV-1,% \¥ - 2, %

0.0 -18.385 -18.432 23.53 16.71

0.2 -11.862 -12.087 25.60 16.72

0.4 -5.878 -6.214 24.20 16.71

0.6 -0.772 -1.117 24.39 16.71

0.8 3.172 2.924 25.38 24.11

1.0 5.924 5.77 25.44 25.23

1.2 7.702 7.618 25.43 25.40

1.4 8.769 8.745 25.43 25.42

1.6 9.359 9.336 25.43 25.43

1.8 9.645 9.633 25.43 25.43

2.0 9.729 9.720 25.43 25.43

Часть данных, относящихся к области мерзлого грунта (0 ^ х ^ 0.6) третьего столбца, представляют значения суммарной влажности. На границе промерзания влажность со стороны талой зоны меньше .

Аналогичные распределения температуры и влажности получены и при величине влагонасыщенностп ц = 0.9. В этом случае влажность на границе промерзания меньше примерно па 1.7%, чем ее значения в талой зоне. По сравнению со случаем ц = 0.6 процесс промерзания идет более интенсивно и соответственно получается большее значение глубины промерзания. Такой результат связан с выбором коэффициента теплопроводности в виде функции от влагонасыщенности, которая для значений ц является возрастающей функцией.

Полученные результаты решения второй задачи тепломассопереноса сопоставлены с решениями первой задачи. В момент времени £ = 200 ч с начала процесса промерзания полученные распределения

температуры и влажности представлены в табл. 2.

Из представленных в табл. 2 результатов следует: во-первых, процесс промерзания идет более интенсивно по второй модели; во-вторых, суммарная влажность в мерзлой зоне, найденная по первой модели, значительно превышает ее значения по второй модели.

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

ЛИТЕРАТУРА

1. Чистотипов Л. В. Миграция влаги в промерзающих неводонасыщенных грунтах. М.: Наука, 1973.

2. Фельдман Г. М. Прогноз температурного режима грунтов и развитие криогенных процессов. Новосибирск: Наука, 1977.

3. Гречищев С. Е., Чистотипов Л. В., Шур Ю. Л. Криогенное фпзпко-геологпчес-кпе процессы и их прогнозы. М.: Недра, 1980.

4. Harlan N. L. Analysis of coupled heat — fluid transport in partially frozen soil // Water Resource Res. 1973. V. 9, N 5. P. 1314-1323.

5. Золотарь If. А. Расчет промерзания и величина пучения грунта с учетом миграции влаги // Процессы тепло- и массопереноса в горных породах. М.: Наука, 1965. С. 19-25.

6. Хаин В. Я. Расчет промерзания грунта с учетом миграции влаги в талой и мерзлой зонах // Тр. Днепропетр. ин-та инж. ж/д. транспорта. 1969. Т. 15. С. 65-72.

7. Павлов А. Р., Матвеева М. В. Итерационная разностная схема для задачи тепломассопереноса при промерзании грунтов // Вестн. Самарского гос. ун-та. Сер. естественнонаучная. 2007. № 6. С. 242-252.

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

9. Годунов С. К., Рябенький В. С. Разностные схемы. М.: Наука, 1973.

10. Коновалов А. Н. Задачи фильтрации многофазной несжимаемой жидкости. Новосибирск: Наука, 1988.

11. Верезин И. С., Жидков Н. П. Методы вычислений. М.: Изд-во физ-мат. лит., 1959. Т. 2.

12. Основы мерзлотного прогноза при инженерно-геологических исследованиях / Под ред. проф. В. А. Кудрявцева. М.: Изд-во Моск. гос. ун-та, 1974.

г. Якутск

25 декабря 2008 г.