УДК 519.63:517.958
ЧИСЛЕННОЕ РЕШЕНИЕ НЕСТАЦИОНАРНОЙ ЗАДАЧИ ИСКУССТВЕННОГО ЗАМОРАЖИВАНИЯ ФИЛЬТРУЮЩИХ ГРУНТОВ
М. В, Васильева, Н, В, Павлова
1. Введение
Искусственное замораживание грунтов является универсальным специальным способом при сооружении подземных выработок в водо-насыщенных грунтах. Способ искусственного замораживания грунтов применяется во многих областях строительства: шахтном, промышленном, при разработке полезных ископаемых, коммунальном, гидротехническом, при сооружении метрополитенов и других подземных тоннелей, конструкций. Наиболее приемлемым является способ, заключающийся в создании противофильтрационной завесы из ледопород-ного тела, получаемого в результате замораживания стенки выработки. Для этой цели в северных регионах в зимнее время могут быть использованы низкие температуры атмосферного воздуха, а для интенсификации процесса обычно применяют специальные жидкостные охлаждающие устройства. В настоящее время такие процессы моделируются на основе различных моделей тепломассопереноса в фильтрующих грунтах. Разработанные математические модели некоторых процессов искусственного замораживания фильтрующих грунтов и их решения рассмотрены в работах [1-3] и др.
2. Постановка задачи Рассмотрим процесс искусственного замораживания глубоко залегающих пористых сред (горных пород), насыщенных фильтруюгци-
© 2010 Васильева М. В., Павлова Н. В.
миея подземными водами, с помощью бесконечно длинного кругового цилиндра радиуса Я, на стенке которого поддерживается температура Тс, которая значительно ниже температуры замерзания пластовой воды Т*. Считаем, что температура горных пород в достаточном удалении от замораживающей колонки равна Т > Т*. Допустим также, что течение насыщающей горные породы пластовой воды в начальное время плоско-параллельное и имеет на бесконечности заданную скорость и, перпендикулярную оси рассматриваемого цилиндра. Пренебрегая геотермическим градиентом и влиянием температуры дневной поверхности, будем рассматривать двумерную математическую модель изучаемого процесса в декартовой системе координат.
В силу сделанных выше предположений вокруг замораживающей колонки образуется ледопородное тело. В мерзлой зоне распределение температуры подчиняется нестационарному уравнению теплопроводности
где — область мерзлого грунта (рис. 1). Г0 — стенка замораживающей колонки. — поверхность раздела талого и мерзлого грунтов.
Коэффициент Ах уравнения У Л (1) вычисляется по правилу сме-
со следующими граничными условиями:
Т(х,у,г)=Тс, (х,у) €Г0
(2)
Т(х,у,1) = Т*, (х,у) еГх
(3)
си
Рис. 1. Схема расчетной области.
№
где индекс в соответствует скелету, индекс г — льду, а т — пористости грунта.
В талой зоне требуется определить как распределение температуры, так и распределение давления, задающее поле скоростей. Вследствие движения поровой влаги распространение тепла в этой зоне описывается нестационарным уравнением конвективной теплопроводности
с —- — (х—\ + — (х—\+с * (+ (5)
дЬ дх у " дх ) ду у " ду ) г" ¡л у дх дх ду ду ) ' (х,у) е О2, г > 0,
где П2 — область талого грунта. Здесь предполагается, что движение поровой влаги подчиняется закону Дарси, к — проницаемость горных пород, ^ — коэффициент динамической вязкости воды, Срш — объемная теплоемкость воды, Лш — коэффициент теплопроводности воды, а коэффициент Л2 кондуктивной теплопроводности в талой зоне также удовлетворяет правилу смеси
Л2 = (1 — ш)Л8 + шЛш. (6)
Подлежит определению то решение уравнения (5), которое удовлетворяет граничному условию (3) и
Т(х,у,*)=Т0, (х,у) еГ2,г>0. (7)
р2 — внешняя граница области П2.
Для определения положения поверхности, разделяющей талую и мерзлую зоны, воспользуемся уравнением, выражающим баланс энергии:
дТ 1
дп
= —шЬргУп, (х,у) € Г 1. (8)
Здесь Ь — скрытая теплота фазового превращения поровой влаги в лед5 р1 — плотность ль да, аУп — скорость движения границы фазового перехода по нормали.
Теперь для определения распределения пластового давления фильтрующей воды из уравнения неразрывности и закона Дарси имеем
уравнение фильтрации упругой жидкости в деформируемой пористой среде
^ Ы дх (/л дх) ^ ду (/л ду) ' ^ ' ^ ^ 2' > ' ^ ^
где в — коэффициент совместной упругоемкости воды и пористой среды [4].
На границах выполняются следующие граничные условия:
--Iе = о, (ю)
Мду
хе [о,у, у = о,/2, ¿>о, (11)
М ду
= ж = 0'/ь уе[0,у, ¿>о. (12)
Соотношения (1)—(12) при известном распределении давления р(х, у, ;), (ж, у) € О2 однозначно определяют температурное поле и положение поверхности раздела талой и мерзлой зон.
Начальное приближение распределения давления в пласте задаем из решения задачи об обтекании кругового цилиндра плоскопараллельным потоком, перпендикулярным к оси цилиндра [5]:
к ( Д2 \ Р(х,у,0)=р0 + и-у[1 + -г-—5- , (ж,у)е О, (13)
М у ж + У /
Где р0 — давление в центре замораживающей колонки.
Чтобы сделать возможным применение разностных схем сквозного счета, объединим с помощью методики А. А. Самарского и Б. Д. Мо-исеенко [6] температурные задачи. В этом случае вместо (1), (5) и условия Стефана на неизвестной подвижной границе (10) получим одно уравнение:
Ср— - — {л—) + — (х—
дЬ дх у дх) ду у ду + + („)б„, (>0. ,14,
Здесь
Л
{Л:: Т^Т:: т—т.,.{
1, Т > Т., О, Т < Т.,
Ср = Ср + т
Ьры + СрыТ* 1--
Рш
¿(Т — Т.
Ср
Срь Т > Т., Ср2, Т < Т.,
, Т. — < Т < Т. ,
Ъгр — т ) = ) Д1+Д2:
I О, (Т < Т. — Д1)и (Т>Т. + Д2).
Поскольку коэффициенты при конвективных слагаемых уравнения (14) зависят от производных функции давления по пространствен-ху
евую задачу (7), (14) с краевой задачей (9)—(13), определенной в области со сложной геометрией.
Используя метод фиктивных областей [7], для функции р(х,у,г) ставим краевую задачу
~др = д_ (~ <9р\ + д_ (~ др
дг дх у дх / ду у ду
(х,у) е о2, г > о,
(15)
с граничными условиями др
-к-^- = и, х = о,/ь у е [о,/2], г > о, = о, ж е [о, гх], у = о,/2,
(16)
Здесь
в =
1/е, Т < Т.
е, (х,у) еП\П2,
(17)
е
3. Разностная задача
В области
Л = {х|х = (ж1, х), 0 ^ ха ^ 1а, а = 1, 2}
введем равномерную сетку ^ с шагами Н\ и Л-2 по пространственным
хх
По времени вводим квазиравномерную сетку
^Т = {¿П = ¿п-1 + Тп, гг = 1,2,...; ¿0 = 0},
Г 1,1 • Тп-1 , Т„_1 < Т*,
Тп = 1
[ т* в противном случае.
Используя безусловно устойчивую чисто неявную дискретизацию, запишем конечно-разностный аналог уравнения теплопроводности (14) в операторном виде с аппроксимацией конвективных слагаемых односторонними разностными отношениями «против потока»:
Ау = /, г е 57т, (18)
здесь
А=Л + С, (19)
где
= = - J, а = 1,2,
а=1
т т
к л
Vy = —Cpw— } {'П(Рха)РхаУха + (1 - 'П(Рха))РхаУха) ■ ^ а=1
Начальное условие имеет вид
у(х1,х2,0) = Т0, x£ZJh. (20)
Начально-краевую задачу для давления (15)—(17) аппроксимируем на той же пространственно-временной сетке:
Лр = f, (21)
здесь
ЛР = — (кр8„), а = 1, 2,
(22)
(23)
Переход с (п — 1)-го на п-й временной слой осуществляется в следующем порядке: сначала решается дискретный аналог краевой задачи для распределения давления (21)-(23), а затем дискретный аналог краевой задачи для распределения температуры (18)-(20) с помощью двухслойного неявного метода минимальных невязок [8], который применим для уравнений с несамосопряженной матрицей, полученной вследствие конвективного члена. Для внешних итераций используется метод простых итераций по нелинейности.
4. Результаты численных расчетов
Построенный вычислительный алгоритм реализован в виде комплекса прикладных программ. Расчеты проводились при следующих числовых значениях входных параметров задачи:
Л2 = 1, 77 Вт/(м -К); Л1 = 1, 55 Вт/(м - К); т = 0, 2; Ср! = 2,43 - 10® Дж/(м3- К); Ср2 = 1,94 - 10® Дж/(м3- К);
Ср- = 4,2 - 10е Дж/(м3-К); Ь = 3,34 - 10б Дж/кг; р = 920кг/м3;
Т. — Т
к/р = 1,7 - 1 О-0 м2/(Па - с); в = О-0/Па;
р0 = 10® Па; е = 10"25; к = 100 м; /2 = 50 м; Я = 3 м;
П = 201; п2 = 101; т0 = 1000 с т. = 86400 с.
В рассматриваемой задаче представляет интерес динамика границы фазового перехода во время замораживания при различных начальных данных температуры и скорости фильтрационного потока.
На рис. 2-7 представлены распределения температуры или давления в момент времени I = 200 суток.
Так, на рис. 2, 3 представлены графики распределения температуры при заданной температуре на стенке замораживающей колонки
*
Рис. 3. Распределение температуры при скорости ^ = 0,3 м/сут.
Тс = 243 К при скоростях набегающего фильтрующегося потока на левой границе, принимающего значения 0.1 и 0.3 м/сутки. соответственно. Результаты счета показывают, что фильтрующаяся жидкость за ледопородным телом сносит температуру тем дальше, чем больше ее
Рис. 4. Распределение давления при скорости V = 0,1 м/сут.
Рис. 5. Распределение температуры при скорости V = 0,1 м/сут.
70 -60 -50 -40 -
10 Г
о _I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_
О 25 50 75 100
х
Рис. 6. Распределение температуры при скорости V = 0,2 м/сут.
60 -50 -40 -
10 -
I —\—!—:—|—1—|—&—!—I—!—1—|—|—I—|—:—|—|—.—
% 25 50 75 100
х
Рис. 7. Распределение температуры при скорости V = 0,4 м/сут.
скорость. На рис. 4 приведено распределение давления при скорости
0.1.м/сутки. Из рисунка видно, что при встрече с препятствием и за препятствием в виде ледопородного тела скорость потока выходит на нуль (движение фильтрующейся жидкости останавливается).
Представленные изотермы значений температуры на рис. 5-7 показывают, что при увеличении скорости фильтрационного потока распределение температуры вокруг замораживающей колонки имеет большее отклонение в сторону направления скорости течения насыщающей горные породы пластовой воды.
Численные расчеты показали, что и внутренние и внешние итерации сходятся очень быстро. С ростом температуры ледопородное тело принимает продолговатую форму.
ЛИТЕРАТУРА
1. Бондарев Э. А., Васильев В. И. Искусственное замораживание фильтрующих грунтов // Численные методы решения задач фильтрации многофазной несжимаемой жидкости. Новосибирск, 1987. С. 38-47.
2. Васильев В. И., Максимов А. М., Петров Е. Е., Цыпкин Г. Г. Тепломассоперенос в промерзающих и протаивающих грунтах. М.: Наука, 1996.
3. Сапунов Н. Е. Исследование процесса промерзания полусферической подземной ледопородной емкости, размещенной в фильтрующем пласте // Численные методы решения задач фильтрации многофазной несжимаемой жидкости. Новосибирск, 1980. С. 228-235.
4. Щелкачев В. Н. Упругий режим пластовых водонапорных систем. М.: Госто-птехиздат, 1948.
5. Валландер С. В. Лекции по гидроаэромеханике. Л.: Изд-во Ленингр. ун-та, 1978.
6. Самарский А. А., Моисеенко В. Д. Экономичная схема сквозного счета для многомерной задачи Стефана // Журн. вычислит, математики и мат. физики. 1965. Т. 5. С. 816-827.
7. Вабищевич П. Н. Метод фиктивных областей в задачах математической физики. М.: Изд-во Моск. ун-та, 1991.
8. Самарский А. А., Николаев Е. С. Методы решения сеточных уравнений. М.: Наука, 1978.
г. Якутск
31 мая 2010 г.