Численное решение нестационарной задачи искусственного замораживания фильтрующих грунтов Текст научной статьи по специальности «Физика»

УДК 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 г.