Численный прогноз динамики криопэгов в мерзлых грунтах Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 9, № 2, 2004

ЧИСЛЕННЫЙ ПРОГНОЗ ДИНАМИКИ КРИОПЭГОВ В МЕРЗЛЫХ ГРУНТАХ

П. П. Пермяков Институт физико-технических проблем Севера СО РАН,

Якутск, Россия e-mail: p.p.permyakov@iptpn.ysn.ru

A mathematical model of salt, heat and moisture transfer in freezing and thawing frozen soils is presented. This model was solved numerically using the directed differences scheme. The results of numerical simulations and the pollution dynamics in the active layer of frozen soil are presented.

Введение

В последние годы в связи с потеплением климата и техногенным загрязнением грунта в районах многолетней мерзлоты усилился процесс образования криопэгов. Криопэги — высокоминерализованные отрицательно-температурные подземные воды. Концентрация в них порового раствора может достигать 200 г/ли более, он не замерзает при -20 °С и ниже. Появление криопэгов вызывает заболачивание территорий и потерю несущей способности грунта. Причинами образования криопэгов являются нарушения планировки местности при строительстве, проникновение в грунт промстоков при аварийных выбросах (высокоминерализованных рассолов карьеров, нефтепродуктов и других экологически опасных загрязнителей).

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

1. Математические модели

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

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2004.

При первом подходе возникает "проблема задания граничного условия" для влагопере-носа на подвижном фронте фазового перехода со стороны талой зоны, подробно освещенная в работах [1, 2]. Эта проблема автоматически снимается при втором подходе [2].

Для простоты изложения рассмотрим одномерный случай. Математические модели конвекции-диффузии при фазовом переходе порового раствора имеют вид [2]

дТ д {дТ\ дТ ЗШЛ дт дг \ дг) дг дт

двв д ( д(Р — г)\ двл

дШв дт

дШвОв дт

__д_(кШ>

дг \ дг

г г

к

Ф

дШл

г т дУСв. дШлСл

г

т

0 < г < Н, тт > т > 0.

Уравнения (1)-(3) замыкаются равновесным уравнением незамерзшей воды

= Жнв = Шнв(Т,Ш,С)

(1)

(2) (2*) (3)

(4)

и условием замерзания соли вместе с почвенной влагой:

Сл кзхСв1

(5)

где с, св — объемные теплоемкости грунта и воды, Дж/(м3-К); Т — температура, К; А — теплопроводность, Вт/(м-К); т — время, с; Ь — удельная теплота фазового перехода, Дж/м3; Ж — суммарная влажность, % (масса всей влаги/масса скелета); 9в, 9л, Wв, Шл — объемные (отношение объема влаги к объему скелета) и весовые влажности воды и льда, %; Н = Р — г — напор, м; Р — всасывающее давление влаги, м; г (ось г направлена вертикально вниз) — пространственная координата, м; к — коэффициент диффузии, /с; кф — коэффициент фильтрации, м/с; О — коэффициент конвективной диффузии

м

примеси, м2/с; У — скорость фильтрации, м/с; Св, Сл — концентрации примеси в воде и во льду, %; кзх — межфазный коэффициент распределения.

Уравнение (1) учитывает процесс промерзания-протаивания порового раствора с учетом фильтрации жидкой фазы. Движение самого порового раствора (воды) с учетом льдовыделения (4) также описывается аналогичными уравнениями параболического типа (2), (2*). Выражения (3), (5) характеризуют процесс солепереноса в промерзающих-протаивающих грунтах. Выражение (2) называется потенциальной формой уравнения влагопереноса, а (2*) — влажностной. Параметры выражений (1)-(5) находятся из экспериментов с использованием методов обратных задач.

2. Идентификация параметров модели

Недостающие параметры математической модели (теплофизические и массообменные характеристики, особенно в зоне промерзания-протаивания, граничные условия) определяют

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

т£ 3(и) = т£ (Л(и) + 32(и)),

«еято и&ято

где

Л(«) = Е /рДт)(Т (и,т) - ТЭ(т))2^т,

г=1 0

пш Ч

3г(и) = ^ / р,(ж)(Ж, (и, ж) - Ж;(ж))2^ж, •?=1 0

Тэ(т) — замеры температуры по времени в г-й точке образца; Ж? (ж) — замеры суммарной влажности по длине образца в ]-й момент времени; пт — число датчиков для замера температуры по длине образца; п^ — число замеров суммарной влажности по времени; рДт), р,(ж) — весовые множители соответственно с размерностями К-2с-1, (%)-2м-1; ТЦм, т), Ж,-(и, ж) — расчетные значения температуры и влажности. Идентификации параметров модели с учетом процесса промерзания-протаивания порового раствора и результаты численного эксперимента более подробно приведены в работе [1].

3. Методы прямой задачи

Рассмотрим численные методы решения прямых задач тепломассообмена с учетом фазового перехода порового раствора.

I. Метод направленных разностей (для влажностной формы). Для решения системы дифференциальных уравнений (1)-(5) в области П = [0, Н] * [0, тт] введем неравномерную сетку. На множестве внутренних узлов систему уравнений аппроксимируем разностной схемой [2]:

СэфТг + АТ + с^Т =

Ж + АЖ = ^2,

$= + + = ^э,

где

АТ = - (АТ*)*, = - )*)*, = - )*).;

С1Т = + V+Т-), С25 = ((V+ (V)-);

= ЬгЖ-, ^2 = ^э = 0, Сэф = с + Ь———' 5 = (ЖвСв + ЖлСл)/100, Ж = + Жв; Жв = Жнв(Т), С = Св + СЛ, Св = V = V+ + VV + = 0, 5 (V + IV|) > 0, V- = 0, 5 (V -IV|) < 0; г = г(Т) = Жл /Ж; щ = ^(Т) = 1 - г(Т) > 0, П2 = П2(Т) = 100/ (к^ + Жв(Т)) .

При этом следует отметить, что разностные операторы С1? С2 построены с учетом знака скорости фильтрации.

пт Т

II. Метод линеаризации Ньютона. Этот метод применяется для решения задач теплосолевлагопереноса в потенциальной форме. В левой части уравнения влагоперено-са (2) применяется метод линеаризации Ньютона для локализации кривой водоудержива-ния = ШВ/<9Р):

s + 1

s V

f, P - p , ег - ег f--1--

т

т

s+1

кф(1 - i) P

- k

ф

(6)

где s — номер итерации.

Численное решение системы уравнений с трехдиагональной матрицей проводится методом прогонки.

III. Метод с двусторонним ограничением решения. Алгоритм с двусторонним ограничением решения применяется для решения задачи влагопереноса (метод Семенова).

В реальных условиях влажность мерзлого грунта меняется в пределах от прочносвя-занной влажности (максимальная гигроскопическая влажность) Wnc = W1 до полного насыщения Wn = W2 (полная влагоемкость). Особенность предлагаемого численного метода — строгое выполнение для решения W(x, т) неравенства W1 < W(x, т) < W2 при любых режимах фильтрационных течений [5]. Основываясь на идеях этой работы, предложим для решения уравнения во влажностной форме (2*) следующий метод:

0т öz

д Л önA dk

öz J + Wl öz \ öz ) öz

ф

0 < z < L, 0 < т;

k

^ + Wik^ - кф = д(т), z = 0, т> 0;

öz

öz

-k^ - W1kdn1 + кф = 0, z = H, т > 0;

dz dz

Я(z, 0) = ?o(z), 0 < z < H, т = 0;

0 < z < L, т > 0;

= A - W A /^д^Л + dki

5т öz 1 öz / дД dW öz

к^ПТ - + кф = -9(т), 2 = 0, т> 0;

+ иу^ - кф = 0, 2 = Я, т> 0;

п(2, 0) = По(2), 0 < 2 < Я, т = 0,

где с = W - п = Ж2 - Ж; д(т) — плотность потока воды, м/с. Искомая переменная находится из выражения

(7)

(8)

(9) (10)

(11) (12)

(13)

(14)

W(z, т) = W1 + (W2 - W1K(z, т)/(я(z, т) + n(z, т)). Уравнения (7)—(10), (11)—(14) аппроксимируем следующими разностными схемами:

s

s

s

л

Л

z

z

Ят + ^Я - С3Я = , Пт + + С4П = ^5,

(15)

где

2k

—V ' i = 0,

D = D2, D5n H -(k(mn)-)z, i = - 1,

2kN-1+0

(nin)z, i = n,

к

Разностные

операторы ^5, С4 , С3 во внутренних узлах сетки представим в виде: = -[кг+«((П1П)г+1 - (П1П)г)/К+ - £¿-1+0((П1П)» - (ШП)»-1)/К—]/Й», « = 1, N - 1,

C4n = [bi+a((nin)i+i - (nin)i) + bi-1+e((nin)i - (nin)i-i)]/2^i, i = 1, N - 1, 0, Zi+i - Zi < 0, _ ,. , _ Г 0, Zi - Zi-i < 0,

a = a(i) = { ' i+i в = ^(i - 1) = < i - (16)

w \ 1, Zi+i - Zi > 0, H ' \ 1, Zi - Zi-i > 0; v y

Zi = Zi(Z, T) = ci(z, т) + (1 - ni(Ti))ni(Z, T) + Wi; (17)

C3я = [ai+i/2((niC)i+i - (niС)i) + «i-i/2((ni^)i - (ni^)i-i)/2^i, i = 1, N - 1; 0 < a = a(W, Wi) = [k(W) - k(Wi)]/(W - Wi) < 0 < b = b(W2, W) = [k(W2) - k(W)]/(W2 - W) <

Численная реализация состоит из трех этапов:

ii V s , ,s+1 ,s+1

-S+i г W„- I 7~)S/ r ^ r1s i

«Г - Ci)/T + DS( я) - CS( я) = ^4, i = 0, N

, -I v s+i s+i

(nS+i - ni)/T + DS( n ) + C4( n ) = ^5, i = 0, N;

s+i s+i s+i s+i

W(z,t) = Wi + (W2 - Wi) я (z,t)/( я (z,t) + n (z,t)),

где s — номер итерации; " V" означает нижний слой по времени.

Основная идея заключается в согласованной аппроксимации "переноса" C4n и "диффузионного" D5n членов путем учета знака zx — использования направленных разностей, правых или левых значений k и b относительно границ ячеек разностной сетки. Операторы D4, C3 имеют диагональное преобладание по столбцам и, безусловно, монотонны [4], т. е. выполняется принцип максимума. Из принципа максимума следует неотрицательность решения Z. Разностные операторы D5 , C4 построены по методике С.С. Маханова и А.Ю. Семенова [6] и имеют неотрицательное решение (nin > 0). Отсюда следует, что n > 0, так как ni > 0. Таким образом, доказано свойство, что при любых (z, т) £ (0, H) * (т > 0) выполняются неравенства Z(z, т) > 0, n(z,T) > 0. Кроме того, вышеуказанные алгоритмы обладают свойством консервативности и безусловно устойчивы [6].

Этапы (15)—(17) выполняются до сходимости итераций. В работе [5] доказано, что итерационный процесс сходится со скоростью геометрической прогрессии.

IV. Метод Кабаре. Уравнение (3) расщепляем по физическим процессам: диффузионному и конвективному солепереносам. Согласно методу расщепления по физическим процессам переход с одного временного слоя на следующий осуществляется в два этапа.

На первом этапе решается уравнение конвективного переноса:

д^ д

дТ = - Ж <V'"2S>, (18)

Я>(Я,т„) = £о(тп), т е (тп, Тп+1),

Я(г,т„) = Я(г,тга).

На втором этапе ищем решение краевой задачи для уравнения диффузии:

дт=* к^>)• (19)

Для численного решения задачи (??) воспользуемся разностной схемой "Кабаре" [7]: трехслойная разностная схема с расщепленной временной производной и направленными разностями по пространственной переменной:

£т+1 _ Я?

г+1 + (1 _ а1)Сг + а1Сг+1, ^ (^-1/2 < 0) and (^+1/2 < 0), 1 2'

-^+1/2(П2Я), _ ^-1/2(П2я%- _ ^я»^ + Я»^) + (1 _ («1 + «2)/2)Сг +

+ ^С-1 + ^С+1, if (^-1/2 > 0) and (^+1/2 < 0),

т 2 2

_2Уг-1/2(п2Я), _ Я^-1 + (1 _ «2)Сг + «2^-1, if (^-1/2 > 0) and (Уг+1/2 > 0), (1 _ («1 + «2)/2)Сг + О2С-1 + О- С+ъ if (^ < 0) and (%1/2 > 0),

где аг — весовые множители, г = 1, 2; (п2Я(п25% — левые и правые разностные производные:

(П2Я)г _ (П2Я)¿-1 (П2Я)¿+1 _ (П2Я)г ;

— функция запаса, которая обеспечивает консервативность схемы:

с! = 0, СП+1 = (Я _ я?+1)/т,

где

Ятах), if (Яг >Ятах)),

ЯГ+Я,

Я^, if (Яг <Я(т1п)),

я(тах) = тах(Яп, Яп-1), Я™п) =т1п(Яп,Яп-1).

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

4. Результаты численного моделирования

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

модельной задачи взяты природно-климатические условия Центральной Якутии [2]. Численные решения нелинейных задач осуществляются методом прогонки с использованием итерации. Итерационный процесс заканчивается при выполнении условия

тах|Ук(г)Ь^ < £к,

г

где Ук — невязка уравнения, аппроксимирующего к-е уравнение модели (1)-(3) (к — 1, 2, 3); £к — заданная точность для к-го уравнения; Кг — потоковый шаг сетки по пространству. Для всех уравнений задается одинаковое значение £к = 0.001.

Пример 1. Результаты численного эксперимента, представленые на рис. 1, отличаются друг от друга формой для уравнения влагопереноса и методами его решения. Годовой водооборот на поверхности с учетом атмосферных осадков и испарения равен нулю. Промсток с концентрацией 10 г/л поступает в грунт со снеговой водой в мае.

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

Дисбаланс воды и засоленности рассчитывается по формулам

н н

/Ш(г,гт)(!г — / Ш(г, 0)(г 0 0

т

н

• 100 %,

/Ш(г, 0)(г

=

н н

/ Б (г, Тт)(г — / Бщ (г, Тт)(г 00

н

100 %,

/ Бщ(г,Тт)(г

где Бщ(г,тт) — распределение засоленности, полученное по методу I (тт = 5 лет). На рис. 1 приведены распределения температуры Т, суммарной влажности Ш (незамерзшей воды Шнв), засоленности Б и порового давления Р по глубине массива (Н = 12 м). Области изменений на оси абсцисс подобраны с учетом минимальных и максимальных значений в годичном цикле и таким образом, чтобы не было наложения графиков. Наблюдается расхождение распределения суммарной влажности по глубине исследуемого массива. Это вызвано тем, что коэффициент влагопереноса в случае фильтрации зависит от порово-го давления, а во влажностной форме — от весовой влажности. Влияние влажности на температурное поле незаметно. Баланс воды и соли лучше соблюдается в двух первых методах. Перенос соли осуществляется в основном потоком влаги, поэтому возмущение влажности действует и на распределение соли.

Таким образом, для долгосрочного геокриологического прогноза более целесообразно использовать два первых метода (I, II), а два последних (III, IV) занимают больше машинного времени, так как работают при малых шагах дискретизации по времени.

Дисбаланс воды и соли

Методы I II III IV

, % 0.1717 0.1739 0.0670 0.6545

, % 0.0000 1.4217 1.1481 1.6953

Рис. 1. Распределения температуры Т, суммарной влажности Ш, засоленности порового давления Р и незамерзшей воды Шнв. Июнь: а — влажностная; б — потенциальная формы; в — метод Семенова; г — метод Кабаре.

а

в

а

П р и м е р 2. При аварийном выбросе промстока в деятельном слое мерзлого грунта формируется криопэг. Концентрация промстока равна 20 г/л. На правой границе области учитывается влияние водной среды (озера).

Численная реализация двумерной задачи осуществлена методом покомпонентного расщепления. С течением времени размер области криопэга меняется в зависимости от температуры и дебита промстока (рис. 2).

0 7.0 14.0 21.0 28.0 35.0 Х,М

СО

"10" 1

0 7.0 14.0 21.0 28.0 35.0 X, М

а

Рис. 2. Начальные распределения температуры (слева), незамерзшей воды (справа) и засоленности (внизу) в криопэге (а); состояние криопэга на февраль через три года (б), через десять лет (в).

Приведенный численный эксперимент подтверждает работоспособность и надежность предложенных методов расчета на основе направленных разностей для решения задач конвекции-диффузии с фазовым переходом.

Список литературы

[1] Пермяков П.П. Идентификация параметров математической модели тепловлагопе-реноса в мерзлых грунтах. Новосибирск: Наука, 1989.

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

[3] Алифанов О.М., Артюхин Е.А., Румянцев С.В. Экстремальные методы решения некорректных задач. М.: Наука, 1988.

[4] Вабищевич П.Н., Самарский А.А. Об устойчивости разностных схем для задач конвекции-диффузии // Журн. вычисл. математики и мат. физики. 1997. Т. 37, № 2. C. 188-192.

[5] Семенов А.Ю. Метод расчета нелинейного уравнения фильтрации, обеспечивающий выполнение двусторонних оценок для решения // Журн. вычисл. математики и мат. физики. 1996. Т. 36, № 11. С. 172-175.

[6] Маханов С.С., Семенов А.Ю. Двумерный неотрицательный алгоритм расчета течений жидкости в открытых руслах // Журн. вычисл. математики и мат. физики. 1996. Т. 36, № 4. С. 97-105.

[7] Головизнин В.М., Карабасов С.А. Нелинейная коррекция схемы Кабаре // Мат. моделирование. 1998. Т. 10, № 12. С. 107-123.

Поступила в редакцию 19 ноября 2003 г.