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

УДК 536.24

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

П.П. Пермяков

Институт физико-технических проблем Севера СО РАН. г. Якутск E-mail: p.p.permaykov@iptpn.ysn.ru

Основываясь на методах решения некорректных задач, восстановлены теплофизические и массообменные характеристики промерзающих-протаивающих дисперсных грунтов, разработаны алгоритмы, которые использованы для решения прямых задач тепломассообмена.

Введение

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

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

В механике мерзлых грунтов широкое практическое распространение имеют математические модели многофазных сред [1]. Для получения математической модели тепломассообмена применяют законы сохранения массы, импульса и энергии для каждой фазы, а затем, суммируя, получают общее уравнение сохранения для всей системы. При этом систему уравнений выписывают в рамках следующих предположений [2]: скелет пористой среды - лед, поровая жидкость и воздух (газ) слабосжимаемы; пренебрегают связанной водой; при образовании льда примесь отторгается в объеме незамерзшей воды и в осадок не выпадает. Система уравнений замыкается равновесным термодинамическим условием существования льда, порового раствора и воздуха. В мерзлых дисперсных грунтах большую роль играет физическое состояние связанной воды, которое усиливается в зависимости от типа техногенного загрязнения. В зоне аэрации (в деятельном слое) обычно пренебрегают воздухом и сжимаемостью пористой среды (скелета), льда и воды. Тогда из модели многофазной среды вытекает общеизвестная модель А.В. Лыкова.

DW„

дТ did T \ dT т d W

с-=-1 X-I- ceV-+ L-л-

дт д z \ 3z I dz дт

дв„

(

дт д z

кл (1 - i(T))

д (P - z)

Yz

л

дв

л

дт

(1)

(2)

дт д z

к-

~д7

д К дЖ

д z дт

(2*)

9WeCe дт

-Ц D дС д z I д z

д VC6 д W„C„ д N

дт

д z

дт

(3)

д-Г = р\ се - N |,

дт 6 kd '

0 < z < H,

тт >т> 0.

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

К = = Жне (Т, С) (4)

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

Сл~кжСв. (5)

где С и Св - объемная теплоемкость грунта и воды, Дж/(м3-К); Т - температура, К; X - теплопроводность, Вт/(м-К); т - время, с; Ь - теплота фазового перехода, Дж/м3; Ж - суммарная влажность, % (масса всей влаги/ масса скелета); в - объемная влажность, % (отношение объема влаги к объему пор); Н=Р-г - напор, м; Р - всасывающее давление влаги, м; г, г (ось г направлена вертикально вниз) -пространственные координаты, м; к - коэффициент диффузии, м2/с; кф - коэффициент фильтрации, м/с; Б - коэффициент конвективной диффузии примеси, м2/с; ¥=(¥„¥) - скорость фильтрации, м/с; Се,Сл - концентрации примеси в воде и льду, %; N - концентрация примеси в твердой фазе, %; в - коэффициент скорости обмена, 1/с; кё (кх) -(межфазный) коэффициент распределения примеси.

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

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

- зависимость от суммарной влажности, которая легко восстанавливается на практике.

д

Ключевыми моментами при построении математических моделей тепломассопереноса в промерзающих и протаивающих грунтах являются: способ локализации области фазового перехода; выбор формы представления уравнения влагопере-носа. Известны два (структурных) подхода: согласно первому из них фазовый переход локализован на поверхности раздела фаз (при определенной температуре); второму - фазовое превращение происходит в протяженной области (так называемая модель фазового перехода в спектре температур). Становление второго подхода связано в первую очередь с экспериментальными работами, которые свидетельствовали о протекании этих процессов в мерзлых грунтах. А.Г. Колесниковым [3] предложена математическая формулировка задачи о температурном режиме грунта, промерзающего в спектре температур. В дальнейшем температурная задача была дополнена уравнением массопереноса.

Вначале рассмотрим, какие трудности встречаются при различных видах (структурных подходах) задания математических моделей теплосолевлаго-переноса. При потенциальной форме записи уравнения влагопереноса возникает проблема представления его в зоне промерзания - протаивания. Влажностная форма математической модели учитывает процесс влагосолепереноса только в талой зоне и поэтому трудно оценить, сколько воды переходит в лед (или, наоборот) при промерзании (про-таивании). А также неизвестно, сколько соли захватывается льдом при промерзании порового раствора грунта. В связи с этим, возникает "проблема задания граничного условия" для влаги и соли на подвижном фронте фазового перехода со стороны талой зоны [4]. Для решения задач теплосолемассооб-мена в насыщенных средах В.М. Ентов, А.М. Максимов, Г.Г. Цыпкин [5] успешно применили второй подход (с зоной промерзания - протаивания). При этом проблема граничных условий первого подхода автоматически снимается, но вместе с тем возникает проблема неопределенности массообменных характеристик в зоне промерзания - протаивания. В связи с этим они с большим удивлением писали "почему все это не было детально изучено ранее?". Причиной этого считаем отсутствие соответствующих приборов и установок для экспериментального изучения и методов идентификации параметров модели. Адекватность вышеуказанных моделей проверена в лабораторных и натурных условиях и имеет хорошую сопоставимость [6]. Параметры выражений (4) и (5) находятся из эксперимента.

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

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

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

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

- Во-первых, эксперимент проводится в цилиндрической ячейке (сосуде) достаточно малого радиуса (Л<0,015 м, высота h>6R), такого, что поток влаги, мигрирующий по радиусу, незначителен.

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

- В-третьих, при проведении эксперимента должны быть соблюдены условия теоремы единственности решения обратной задачи теплопроводности Н.В. Музылева [8]. Температурное поле исследуемого образца в цилиндрическом сосуде описывается следующим уравнением теплопроводности:

ст)рдТ = ()дТ) 0<т<т™■ 0<г<я (6)

от г о г I о г )

11т г*Х(Т ^

г —>о д г

= 0'

- г*х{т

д г

= г'чТ)

Т (г ,0)= Т0 (г )

0 <т<тт

0 <т<тт

0 < г < Я,

где

с(Т) = Сэф = еск + 0 + (св - с,)Шнв(Т) + Ь

дТ

(7)

(8)

(9)

(10)

Ш (Т) - ш

х(т) = х + (х-х)^^—(11)

V/ м V т и) ш0 - ш

ПС

V = 0, 1 - соответствует пластине и цилиндру.

На границе при г = 0 выполняется условие ограниченности решения (7).

Требуется восстановить параметры (с,^,Ст,Х,^,Хт,Шн1Т))Т=(и1,иг,...,ия)т=и по известным замерам температуры ТЭ=ТЭ(г,т) в точках Г; (;=1,ит).

Задачу (6-9) сформулируем в виде задачи оптимального управления. В качестве оптимальности выберем целевой функционал среднеквадратичес-кого отклонения замеренных температур от расчетных значений:

г=0

г=Я

1 (и) = X }Р, Т (и)- Тэ) дт = \\т - Тэ

и е и с Л'

(12)

где д - весовые множители; Т(и), Т? - расчетные и замеренные температуры; пТ - число замеров по радиусу образца; СеХ2Пт[0,тт] - пространство замеров.

Разработана следующая последовательность отыскания неизвестных: из термограммы опыта выбираются области мерзлых и талых зон и в каждой из них определяются коэффициенты (смДм) и (сТДТ); затем с помощью (10) вычисляются сск,^с; на последнем этапе, используя все замеры температуры, восстанавливаем параметры функции незамерзшей воды.

Минимизацию целевого функционала можно осуществить методами скорейшего спуска, сопряженных градиентов или модифицированным методом скорейшего спуска [7, 9].

График восстановленной функции количества незамерзшей воды лежит между кривыми 2, 3 (рис. 1), полученными в результате расчета количества незамерзшей воды, начиная с талой (баланс воды) и мерзлой (баланс льда) зон, соответственно. Некоторые различия между ними наблюдаются в области начала и конца фазовых переходов. Расхождение в начале фазовых переходов можно объяснить следующим образом. При температуре 273 К замерзает значительная часть свободной воды и кривая незамерз-шей воды имеет "скачкообразный" вид. Этот "скачок" при дискретной реализации приближенно аппроксимируется дельта функцией (линейно).

V У /

/2

---2

Т, К 262 266 268 27/ 274

Рис. 1. Сравнение результатов восстановления функции количества незамерзшей воды, полученным численным (1) и квазистационарным (2, 3) методами

В вычислительной теплопередаче широкое практическое применение получил метод сглаживания [10, 11]. Это равносильно предположению, что фазовый переход происходит не при одной определенной температуре Тф, а в некотором интервале температур, длина которого определяется величиной параметра сглаживания А. При численном эксперименте обычно задают различные значения этого параметра, при этом изменяется ширина диапазона температур фазового перехода. С другой стороны, мощность и средняя температура источника тепла фазового перехода остаются постоян-

ными, но меняется вид функции плотности распределения источника тепла.

Мощность источника тепла выражается условием: ~ 1 (Т, С)

\8(Т, А)дТ = 1

или

О

(Щ - Щ„)

О

д Т

-дТ = 1.

В реальных условиях плотность распределения зависит от дисперсности, минерального состава грунтов, типа растворенных ионов и т.д. В результате обработки данных экспериментальных исследований при засолении песка хлоридом натрия получены различные виды графиков плотности распределения источника тепла (рис. 2). С повышением концентрации растворенных солей усиливаются растворение льда и понижение средней температуры фазового перехода. Идентификация параметров функции количества незамерзшей воды не занимает много времени и труда [9].

Аппроксимационные формулы 5-функции соответственно имеют вид (рис. 2): а)

8(Т, А) = <• 2А

|о,

Т <А Т >А

в)

с)

8(Т, А) =

8(Т) =

3 Т 2 1А (1 - ^)' 4А а2

О

Т <А Т >А

дЩш д Т

-Щ - ) =

1

(Щ - Щс )

о,

,|Ь+1

аЬ / |Т - 2731Ь (Щ - А)/(Т3 -Т2), О,

Т < Ть

Т < Т < Т2, Т2 < Т < Т3, Т3 < Т,

д Щ

8(Т) = /(Щ - Жпс) =

1

(Щ - Щпг )

дТ О,

(А1- ЖпС)/(Т2 - Т1), аЬ /\Т - 273 |Ь+\ (Щ - А) /(Т4 - Тз), О,

Т < ТЪ Т < Т < Тъ

Т2 < Т < Т3, Т3 < Т< Т4, Т4 < Т,

где А, а, Ь, А, А1, Т (/=1,2,3,4), Щ, Жпс - параметры функции количества незамерзшей воды.

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

в

0 7 -----------;------

0.6 ""!...........:......

: ! : ! ! : !

! 1 ! ! ! ....;...........!...........1..........;...........1...........[..........

| | | | |

1 ....:...... 1 ; 1 ........:........ 1 ;

; -Ч...... ; ----...... ! .....!. ; ! ! ......и- .... -

! ....]______ -!...... ______ ________!________ -! -- ______[... -- -

; 1 ; ; ; ; ;

; 0.4 -—■1...........:...........-г..........;...........:...........',—

о.з -■■;...........:......

о.1 -—■;...........;.....................;...........;...........---

___________1.1

: :

1 1 """!...........1...........Г..........1........... ; 1 ! ........... ;

11 1

'""": 1 ; : 1 ! ! 1 ! !

'""": ! ; : 11 V" \

...........1...........;..........!........... 1 1 ...........1...... !

I 1 ! !

1 ; ] 1 1 1 1

| 1 ! ! ! 1 1 1

! -ч........... | ! ! ! ...........1..........1...........1........... ! ! !

; ! ; - 1 ----'

260

т. К

С

1 1 1 1

! -1----------- 1 -1........... 1 ....;........... : : 1 ! ! !

: : : ! ! ! ...........}..........1...........1........... 1 1 1 ...........■..........■...........!...........

! ...:________Г ! ! ! -1________:__________!___________!___________

: 1 ! : —^

е f

Рис. 2. Плотность распределения источника тепла: а) прямоугольник; Ь) парабола; с) Св = 5; в) 15; е) 20; f) 30 %

Сама граница фазового перехода явно не выделяется и не участвует в построении разностной схемы, а при необходимости свободная граница идентифицируется как изотерма средней температуры фазового перехода Тф;:

^ 1

' ф ш - ш

'' п '' м.

7 дШнв (Т, С) , ■ | т—нвК ' 'аг,

дТ

3. Численная реализация

Для численного решения системы дифференциальных уравнений (1-3) в области 0=[03]х[0,тт] введем неравномерную сетку:

тЬх = ШИ хо>т, ШИ = {хо = X = Х-1 + Н, ' = 1, ат= { т] = ] т, ] = 1,2,3,...}, к, = (Н++ И~)/2, ' =1, N -1.

На множестве внутренних узлов ю4={х;|, /=1,Л—1} систему уравнений аппроксимируем разностной схемой [12]:

СэфТг + °1Т + СВС1Т = Фl,

* + 1 * * Р—Р_ 1

¡л —-'- + -

+ А Р + (кф )Л=Ф2,

+ 035 + С2 5 = ф3

Ь

а

где

Б? = -(ХТ-)-, БЩ = -(кфп(Р\)-, В'^Щ = -(ЩЩ)-)Б3Б = -(Б{щ Б)-);, С? = (V Т; + V+Т-), С25 = ((V-щБ); + (У+щБ)-),

Ф!= МЩ-, Ф2 =Фз = 0 Б = (ЩеСе +ЩлСл )/100, Щ = Щл + Щ, в=вл +вв, Ще = Щне(Т), Св =п2Б, V = V ++V-, V+ = 0,5 (V + V) > 0, V" = 0,5 (V - ^|)< 0,

I =

Щ,

П = 1 -и Ъ = 100/(кЩ + Щв).

Разностную схему температурной задачи получили следующим образом. Используя функцию незамерз-шей воды We=Wue(T,W,S) и учитывая независимость потенциалов Т^Б, множитель теплоты фазового перехода дWJдт ур. (1) можно записать в виде:

дЩ

дЩ, дТ + дЩ.Ж +дЩ, дБ

д Т д Т

д Щл д 1(Т )Щ д Щ = д Щ д Жл д (Щ - Жв)

д Т '

= КТ),

д Щ

д Б

д Б

д Б

получим: д Щ

дт

д Щ д Т . ,„, д Ж д Щ д Б -+г(Т )-

дТ дт

дт дБ дт

зии и солепереноса имеют диагональное преобладание по столбцам, и достаточные условия устойчивости выполняются в гильбертовом пространстве ¿„(м), ¿1(м) [13].

Динамика температуры

5,0 10,0 г 15,0 20,0 25,0

Динамика влажности 0,0 5,0 10,0 г 15,0 20,0 25,0

дт дТ т дЩ дт дБ дт'

Отсюда при помощи частных производных функции W;

д Щл д (Щ - Ще) д Жв

0,0

2,6

5,2 г

7,8

10,4

13,0

17 20 17

\ А

17,20,25,

Рис. 3. Распределения температуры и суммарной влажности при Св=0

Динамика температуры 5,0 10,0 г 15,0 20,0 25,0

Первое слагаемое данного выражения учитывается в эффективной теплоемкости (10). Последним членом выражения можно пренебречь, т.к. его вклад в реальных условиях достаточно мал (дWuJдS»0).

В левой части уравнения фильтрации (2) применяется метод Ньютона для локализации кривой водоудерживания (р=дв/дР). Введение 5-функции для приближенного решения задачи Стефана позволяет строить разностные схемы со сглаженными коэффициентами, т.е. совершается переход к обычной задаче параболического типа. Преимуществом данного подхода является то, что не нужно выбирать параметр сглаживания Д. Численная реализация нелинейной задачи (1-3) осуществляется с помощью итерационных методов сквозного счета, что намного упрощает процесс решения многомерных задач тепломассообмена в мерзлых грунтах при использовании экономичных аддитивных локально-одномерных разностных схем. Через функцию количества незамерзшей воды выражаются доли порового раствора п(Т) и засоленности ц2(Т), которые участвуют в процессе фильтрации (диффузии). Разностные операторы уравнений диффу-

Динамика влажности 5,0 10,0 Г 15,0 20,0 25,0

0,0

2,6

5,2

1.

7,8

10,4

13,0

Динамика засоленности 0,0 5,0 10,0 Г 15,0 20,0 25,0

■-М I I г

Рис. 4. Распределения температуры, суммарной влажности и засоленности при Св=100 г/.л

Для проверки возможностей математической модели приведем пример. Основная цель - методическая: показать какого рода численные исследования подземных фильтрационных процессов можно проводить в промерзающих - протаивающих грунтах и насколько эффективно работает при этом предложенный алгоритм. Алгоритм решения при промерзании-протаивании мерзлого грунта построен с учетом знака скорости фильтрации с помощью направленных разностей. Численные расчеты проведены с учетом природно-климатических условий, радиационно-теплового баланса поверхности и массообменных характеристик грунтов Центральной Якутии [6].

Рассмотрим случай, когда в грунт поступает (в летний период) надмерзлотная загрязненная грунтовая вода при различных концентрациях. Техногенная вода, загрязненная промстоками, проникает по полосе 0,5...0,75 м с напором - 0,5 м и фильт-

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

1. Нигматулин Р.И. Основы механики гетерогенных сред. — М.: Наука, 1978. —336 с.

2. Васильев В.И., Максимов А.М., Петров Е.Е., Цыпкин Г.Г. Теп-ломассоперенос в промерзающих и протаивающих грунтах. — М.: Наука, 1997. —224 с.

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

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

5. Ентов В.М., Максимов А.М., Цыпкин Г.Г. Образование двухфазной зоны при промерзании пористой среды. — М.: ИПМ АН СССР. Препринт № 269, 1986. — 56 с.

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

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

рирует вниз в сторону реки. На рис. 4 приведены распределения температуры, суммарной влажности и засоленности грунта (кг/м3) в конце июля при коэффициенте фильтрации kф = 0,48 м/сут. С повышением концентрации промстока процесс растворения льда усиливается, что сопровождается понижением температуры, а также увеличением подвижности жидкой фазы. Численный эксперимент проведен с учетом сезонной динамики уровня реки. Из рисунка видно, что поступление техноген-но-загрязненной воды действует на динамику температурного и влажностного режимов массива.

Вывод

Предложенные алгоритмы решения и программные средства можно применять для прогноза динамики техногенного загрязнения в мерзлых грунтах.

Работа выполнена при частичной поддержке РФФИ, проект № 03-05-65408.

8. Музылев Н.В. О единственности одновременного определения коэффициентов теплопроводности и объемной теплоемкости // Ж. вычисл. матем. и матем. физ. — 1983. — Т. 23. — № 1. — С. 102—108.

9. Будак Б.М., Соловьева Е.Н., Успенский А.Б. Разностный метод со сглаживанием коэффициентов для решения задач Стефана // Ж. вычисл. матем. и матем. физ. — 1965. — Т. 5. — № 5. — С. 828—840.

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

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

12. Пермяков П.П. Численное моделирование задач конвективно-диффузионного переноса при фазовом переходе // Проблемы прочности материалов и конструкций для регионов холодного климата: Тр. 1-го Евразийского симпозиума. 16—20 июля 2002 г. В 4 ч. — Якутск, 2002. — Ч. IV. — С. 149—156.

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