Содержание к диссертации
Введение
Глава I. Моделирование трехмерных задач переноса консервативных загрязняющих примесей в пористой среде 13
1. Математическая модель 13
2. Применение противотокового весового метода конечных элементов 15
Глава II. Алгоритм решения задачи 39
3. Автоматическая подготовка данных 39
4. Описание численного алгоритма решения 48
5. Расчёт моделей переноса консервативных примесей в пористой среде с неоднородными условиями 51
Глава III. Пространственные задачи переноса реагирующих и нереагирующих веществ в пористой среде 68
6.Трёхмерный перенос загрязнений с учетом физической сорбции. 68
7. Задача распространения многокомпонентного раствора в присутствии реакций сорбции и комплексообразования ,. 78
8. Задача пространственного взаимодействия фильтрационного течения с лучевым дренажем S9
Литература 102
- Применение противотокового весового метода конечных элементов
- Описание численного алгоритма решения
- Расчёт моделей переноса консервативных примесей в пористой среде с неоднородными условиями
- Задача распространения многокомпонентного раствора в присутствии реакций сорбции и комплексообразования
Введение к работе
Актуальность темы. Проблема охраны водных запасов в связи с опасностью загрязнения подземных вод в результати антропогенного и техногенного воздействии на окружающую среду приводит к необходимости совершенствовать методы прогнозировании миграции загрязняющих примесей. Для описания распространения иримесен применяются в основном двумерные модели, построение которых связано с осреднением характеристик течения по одной из координат. Однако, при наличии в области течения участков с резко изменяющейся гидравлической проводимостью, действующих источников и стоков, при сложной конфигурации границ адекватное описание фильтрационного процесса возможно только в рамках трёхмерных моделей. Дополнительным аргументом в пользу применения трёхмерных моделей служит сам пространственный характер явления дисперсии. При исследовании переноса пекопсервативных примесей необходимо также учитывать физико-химическое взаимодействие компонент раствора между собой и с пористой средой.
Прогнозные расчёты переноса примесей базируются на использовании математических моделей, учитывающих движение грунтовых под и конвективно-дисперсионный характер распространения присутствующих в растворе веществ. Если в составе растворов содержатся несколько компонент, то взаимодействие их в виде ионного обмена, сорбции, ком-нлексообразования может существенно влиять на характер распределения концентраций. В этом случае принимаются во внимание кинетические соотношения, отражающие физико химические процессы, математические модели усложняются,
они становится нелинейными, что в свою очередь приводит к дополнительным трудностям при решении соответствующих начально-краевых задач. В основном в литературе но многокомпонентному переносу хорошо исследованы одномерные и двумерные модели с постоянной скоростью фильтрации.
Особую сложность представляет решение задач с преобладающей конвективной составляющей. Развитие вычислительной техники в последнее десятилетие способствовало усилению интереса к пространственным задачам как у нас, так и за рубежом.
Цель работы. Создание трёхмерных математических моделей массопереноса, учитывающих различные неоднородные условия, физико-химические механизмы взаимодействия компонент раствора. Разработка эффективных численных методов реализации этих моделей. Анализ поведения характеристик процесса в зависимости от исходных параметров.
Научная новизна. Разработанные математические модели переноса загрязнений с учётом взаимодействия веществ в системе "растпор-порода"оригинальны. Впервые применяется математическая модель взаимосвязанных потоков грунтовых вод и потоков в лучевом дренаже. Разработаны экономичные вычислительные алгоритмы, позволившие получить свободные от осцилляции решения систем уравнений. Проведены численные исследования поведения решений в зависимости от физических и химических параметров, входящих в основные системы уравнений.
Методика исследований. В работе используются методы конечных элементов, конечно-разностных схем, методы решения систем линейных алгебраических уравнений, метод ав-
тематического разбиения областей.
Практическая ценность. Разработанные математические модели и численные методы их реализации позволяют исследовать распространение примесей на фоне усложнённой фильтрационной картины, обусловленной неоднородностью пористой среды, наличием источников и стоков, а также в присутствии химических реакций.
Апробация работы Результаты работы докладывалась на семинарах и научных сессиях Института гидродинамики СО РАН; на Всесоюзном семинаре «Математическое моделирование гидрогеологических процессов» (Душанбе, 1991 г.); на Первой Всесибирской конференции по математическим проблемам экологии (Новосибирск, 1992 г.); на Второй Всероссийской конференции по математическим проблемам экологии (Новосибирск, 1994 г.); на Международном симпозиуме «Гидрогеологические и экологические процессы в водоемах и их водосборных бассейнах» (Новосибирск, 1995 г.); на Международной конференции «Математические модели и численные методы механики сплошных сред» (Новосибирск, 199G г.), а также были представлены па 22-ой Генеральной Ассамблее Европейского геофизического общества (Вена, 1997 г.) и на 24-й Генеральной Ассамблее EGS (Гаага, 1999 г.).
Публикации. По теме диссертации опубликованы 13 работ.
Объем работы. Диссертация состоит из введения, трех глав и заключения. Объем работы 108 стр., 35 рис., 1 табл. Перечень литературы включает 40 наименований.
Применение противотокового весового метода конечных элементов
После введения соответствующих источников и граничных условий полное поэлементное интегрирование по всей области моделирования дает систему обыкновенных дифференциальных уравнений dC [Л]{С] + [В}{—} = {F}, (49) где [А] — матрица размерности п (п —число вершин конечных элементов), полученная от интегрирования дисперсионного и конвективного членов уравнения (9); [В] — матрица той же размерности, являющаяся результатом интегрирования производной концентрации по времени; {F} — правая часть, объединившая источники и некоторые граничные условия.
Одним из способов решения системы (49) является применение конечно-разностной аппроксимации к производной по времени. Таким образом получаем систему линейных алгебраических уравнений ([А] + Щ){С}1+м = [-{C}t + {F}t+A„ (50) где {C}t+At — искомая неизвестная концентрация на (к+1)-м текущем шаге тго времени, a {C}t известная концентрация тга, предыдущем шаге. Если существует условие первого рода, то его необходимо учесть в системе (50). Один из способов ввода граничного условия является следующее. Запишем (50) в общем виде [Н]{С} = {Щ. (51) Тогда, введя условие первого рода, вектор неизвестных будет содержать часть значений с заданной концентрацией, и систему можно записать в виде (52) I №} [Щи [H]l2 [Н]21 № где {С2} —вектор известных значений. Тогда систему (52) можно записать в сокращенном виде {Н]п{Сі} = {Ri} - [Н]12{С2}. (53) Таким образом можно уменьшить ранг матрицы системы (52) на число узлов, в которых заданы известные значения первого рода.
Весовые противотоковые коэффициенты
Весовые коэффициенты в базисных функциях определяются с помощью так называемых локальных чисел Пекле, то есть чисел Пекле, определенных для каждого элемента в том виде, в котором они были предложены в работах [37, 38]. В трехмерном случае удобно определить эти числа отдельно в горизонтальном и вертикальном направлениях. Рассмотрим произвольный элемент е (рис. 1). Пусть 14 — горизонтальная проекция вектора скорости внутри элемента, отнесенная к основанию призмы треугольнику ijk. Определим Vl3, Уг, VI как проекции VL-Следующее выражение представляет собой локальное число Пекле между узлами і и j и показывает положение этих узлов друг к другу относительно потока. Например, если Tij 0, то узел і находится выше по потоку от узла j. Сумма любых двух чисел Пекле показывает положение внутри треугольника одной вершины против двух. Исходя из выше сказанного, веса определяются следующим образом. В горизонтальном направлении — и _,/,.„ vl M v?\ki\ г l\ \Vh\ 3 vt т - v? ш &i\ Ш yki \кі\ - - vf 1 jk\ Ш\ — 1/6 + Qh и2 = 1/3 + ah ft м \ г "", (54) w3 = l/3 + oft xi\Vh\ — где T4mn — проекция Vh на m — n сторону треугольного основания, m — i,j,k; n — i,j,k, где i,j, к — вершины треугольника Aijk, ij, кг, jk — его стороны (Рис. 1); ah — положительный коэффициент. Весовые — коэффициенты 074, 5: 6 имеют тот же вид, но в качестве Vh берётся проекция скорости внутри элемента на верхнее основание. В вертикальном направлении весовые коэффициенты имеют вид, подобный (54), например, для узла 1(г) (Рис.1): Лі = 1/2 + av -\ (55) щУг,і\ — — где Vzj - вертикальная проекция скорости V на ребро 1-4, примыкающему к узлу ц Az - высота элемента; av — положительный коэффициент. Используя выражение для проекции вектора на отрезок прямой, на-ходим проекции Vh на стороны треугольного основания утп _ VxfijXm Хп) + Уу,к(Ут Уп) \1{Хт - Xnf + (ут - уп)2 Горизонтальные составляющие для Vh ищутся из решения фильтрационной задачи (1)-(2). Имея распределение напора в области решения, по закону Дарси вычисляются пространственные компоненты скорости внутри элемента. Относительно основания каждого элемента принимаем Н(х, у, z, t) = NiH! + N2H2 + АГ3Я3, (56) где Ni — линейные двумерные базисные функции. Тогда составляющие вектора скорости можно представить в виде Vx,h = -КХХ(ЬІЩ + bjHj + bhHk)/2A, (57) Vx, h = -Куу(аЩ + CjHj + ckHk)/2A. Здесь коэффициенты bi, сі (I = i,j,k) удовлетворяют равенствам (19) и (27); Hi — вычисленные значения напора в узлах; А — площадь основания ijk. Вертикальная компонента скорости вычисляется для каждой из вершин i,j,k, например У» = KJ . (58) Следует отметить, что положительные коэффициенты ah, av определяются путем подбора в процессе счета и выбираются такими, чтобы исключить осцилляции в решениях, но численную дисперсию при этом свести к минимуму. В каждом конкретном случае это делается индивидуально и величина коэффициентов зависит от особенностей задачи.
Описание численного алгоритма решения
На верхней границе CQ—І на участке 40 м х 80 м, у—6,5 м,0 z 120 м в течение 5 лет, остальные границы были непроницаемые. Физические параметры были взяты, как в двумерном случае.
Конечно-элементным методом локального баланса было найдено стационарное решение фильтрационной задачи. Область была разделена Рис. 13: Изолинии концентрации при t=5, 8, 12 лет (метод конечных элементов, профильная задача).
Изолинии концентрации при t—5, 8, 12 лет в вертикальном разрезе плоскостью X-Z при у=60 м (участок загрязнения — полоса). на 20 слоев по 105 элементов в каждом слое. На рис. 14 приведены расчётные данные о поведении пятна загрязнения в различные моменты времени в вертикальном сечении плоскостью X-Z при у=60 м. Характер движения подобен двумерному случаю. Шаг по времени варьировался от ОД до 0,5 лет, число итераций составило 4-5 на каждом шаге. Осцилляции в решении исчезают при ah — а„=0,013.
2. Участок загрязнения представляет собой квадратную площадку длиной и шириной 40 м на верхней границе области. На рис. 15 представлено физическое описание этой задачи. Все параметры расчёта, начальные и граничные данные для фильтрационной и конвективно-дисперсионной задач соответствовали предыдущему случаю. Область была разделена на 10 слоев по 243 элемента в каждом слое, в совокупности использовалось сетка из 1573 узлов и 2430 элементов. Шаг по времени брался от 0,1 до 0,2 лет.
На рис. 16 представлено поведение пятна загрязнения в виде изолиний концентрации с течением времени в вертикальном разрезе области плоскостью X-Z по координате у=60 м. По мере попадания загрязнения в область высокой проводимости пятно существенно смещалось вправо. Однако по сравнению с предыдущим примером пятно продвигалось медленнее, что связано с тем, что часть концентрации вещества расходовалось на поперечную дисперсию. Рис. 17 демонстрирует, как весовые коэффициенты подавляют осцилляции в решении. Показаны расчетные распределения концентрации вдоль координаты х при фиксированных z 6.5
Задача 3.Следующий пример представлял собой исследование переноса неконсервативной примеси от объёмного источника загрязнения в виде площадки в неограниченном однородном пласте. Такой источник можно представить системой закачных скважин с общим дебитом. Описание задачи и физические параметры были взяты из работы [35] в которой рассматривается модель расчета линий тока.В области течения размером 150м х 150м х 31м на верхней границе располагался источник в виде прямоугольной площадки размером 12м х 4м х 0, 03м. Дебит за о качных скважин составлял 5зак=1 2 м/сут. На расстоянии 14 м по координате х и 20 м — по у от источника на глубине 3 м от поверхности располагался водозабор в виде откачной скважины с длиной фильтра — =о — -о
На левой и правой границах были заданы значения напора R ± = 100 м и Н\ = 99.1 м, соответственно. Остальные участки границы считались непроницаемыми. Значения компонент гидравлической проводимости постоянны — Кхх = Куу = Kzz = 61 м/сут. При расчетах переноса примеси использовались следующие параметры: пористость среды в — LI—I I IIII llll 11 II II—LLU—l_l l_l LI l_l 1 LI Рис. 19: Конечно-элементная сетка. 0.2, концентрация загрязняющего вещества в источнике С = 100 г/л. С целью исследования влияния дисперсионных процессов на поведение концентрации загрязняющего вещества коэффициенты продольной и поперечной дисперсии варьировались в пределах интервала значений (10 м, 0.1 м). Для иллюстрации ниже приведены наиболее характерные результаты, подтверждающие значительное влияние дисперсии.
Область была покрыта конечно-элементной сеткой,состоящей из 2808 узлов и 4656 элементов. На рис. 19 показана сетка, на которую была разбита область.
При расчетах конвективно-дисперсионной задачи методом подбора было установлено, что при аи = 0.05, av = 0 получаются удовлетворительные результаты, свободные от осцилляции. Расчетное распределение концентрации для случая оц = 10 м, at — 1 м на рис. 20 а), б), где приводится горизонтальное и вертикальное сечения изолиний С, показывает влияние дисперсии на форму пятна загрязнения, которое в этом случае выглядит более расплывчатым по сравнению с распределением С при а; = 1 м, а = 0.1 м (рис. 21 а), б)). В последнем случае след загрязнения имеет более крутой фронт и сосредоточен в верхней части толщи. В этих примерах изолинии С приведены на время t = 100 сут, когда распространение концентрации выходит на стационарный режим.
Расчёт моделей переноса консервативных примесей в пористой среде с неоднородными условиями
Подставляя представление Н в виде (102) в уравнение (104) и интегрируя по всем элементам, получаем систему обыкновенных дифференциальных уравнений ИМ [B]Q} - [А]{Й]. (105) Аппроксимация производной по времени с помощью конечно-разностного приближения, приводит к системе линейных алгебраических уравнений: ([Л] + ){Я}"+1 = {Я}". (106) Конвективно-дисперсионное уравнение (98) решалось также методом локального баланса, но с введением противотоковых весов, чтобы избежать осцилляции в решениях. Для решения задачи движения жидкости в трубах (96)-(97) использовалась конечно-разностная аппроксимация со стабилизирующим членом 6щ. При этом принимались во внимание только стационарные решения. Неявная схема для уравнения движения в трубе имеет вид «(«Г1 - «?) = -{ S f+i p- - Q (107) Полученные системы линейных уравнений для расчета напора и концентрации вместе с конечно-разностной задачей (107) решались совместно итерациями методом Гаусса-Зейделя , при этом учитывались условия сопряжения (101). Итерационный процесс завершался при достижении заданной точности.
Для численного исследования модели (94)-(101) рассматривалась пространственная задача переноса консервативной примеси на фоне течения к лучевому дренажу, размещенному в однородном напорном пласте конечной мощности. Схема модельной задачи приведена на рис. 32. В области размером 300м х 300м х 20 м на левой границе задавался постоянный напор Щ = 100 м, остальные границы считались непроницаемыми. В середине толщи на глубине 10 м от поверхности находился крестообразный водозабор с длиной лучей 60 м и диаметром дренажных труб 0,2 м. Физические параметры задачи задавались следующие; гидравлическая проводимость постоянна во всех направлениях Кх Ку — К2 50 м/сут; 7 — \ — 2,4 106, где п —коэффициент шероховатости; R — , где R— гидравлический радиус, d—диаметр дрены. Напор в области водосборного колодца Нд — 90 м. На правой границе области задавалось постоянное значение концентрации. Начальное распределение концентраций было нулевым; значения коэффициентов дисперсии щ = 0,5 м и at — 0.1 м.
Для сравнения были проведены расчеты задачи, когда напор на дрене задавался постоянным по всей длине. В этом случае решалась только У мм? Рис. 32: Горизонтальная проекция сетки. задача фильтрации, для которой область дренажа представлялась внутренней границей с постоянным значением напора на ней. На рис.33 приведены сравнительные результаты расчётов в виде изолиний напоров в горизонтальном разрезе по координате, соответствующей расположению лучей в плоскости сечения (z = 9,9 м), и в вертикальном сечении плоскостью, проходящей через дрены, параллельные оси X и перпендикулярные оси Z по координате у — 149,9 м (слева — распределение Я, когда течение в трубах отсутствует, справа — с учётом течения).
Учет движения жидкости в дренах меняет фильтрационную картину в области течения. На графике (рис.34) показано распределение напора по длине одной из труб от концевого участка до водосборного колодца. Видно, что распределение и носит неравномерный характер. Дебит скважины действующего лучевого дренажа Qr = 10100 м /сут. При расчётах в предположении постоянного напора в дренах было получено Qr — 11800 м /сут. Все результаты приведены для момента времени, когда режим течения приближался к установившемуся.
Изолинии концентрации в горизонтальном и вертикальном сечениях на различные моменты времени показаны на рис. 35. Форма следа загрязнения определяется гидродинамической картиной течения. Последующие расчёты показали, что при данных условиях задачи лучевой дренаж препятствует дальнейшему распространению примеси.
Задача распространения многокомпонентного раствора в присутствии реакций сорбции и комплексообразования
Подставляя представление Н в виде (102) в уравнение (104) и интегрируя по всем элементам, получаем систему обыкновенных дифференциальных уравнений ИМ [B]Q} - [А]{Й]. (105) Аппроксимация производной по времени с помощью конечно-разностного приближения, приводит к системе линейных алгебраических уравнений: ([Л] + ){Я}"+1 = {Я}". (106) Конвективно-дисперсионное уравнение (98) решалось также методом локального баланса, но с введением противотоковых весов, чтобы избежать осцилляции в решениях.
Для решения задачи движения жидкости в трубах (96)-(97) использовалась конечно-разностная аппроксимация со стабилизирующим членом 6щ. При этом принимались во внимание только стационарные решения. Неявная схема для уравнения движения в трубе имеет вид «(«Г1 - «?) = -{ S f+i p- - Q (107) Полученные системы линейных уравнений для расчета напора и концентрации вместе с конечно-разностной задачей (107) решались совместно итерациями методом Гаусса-Зейделя , при этом учитывались условия сопряжения (101). Итерационный процесс завершался при достижении заданной точности.
Для численного исследования модели (94)-(101) рассматривалась пространственная задача переноса консервативной примеси на фоне течения к лучевому дренажу, размещенному в однородном напорном пласте конечной мощности. Схема модельной задачи приведена на рис. 32. В области размером 300м х 300м х 20 м на левой границе задавался постоянный напор Щ = 100 м, остальные границы считались непроницаемыми. В середине толщи на глубине 10 м от поверхности находился крестообразный водозабор с длиной лучей 60 м и диаметром дренажных труб 0,2 м. Физические параметры задачи задавались следующие; гидравлическая проводимость постоянна во всех направлениях Кх Ку — К2 50 м/сут; 7 — \ — 2,4 106, где п —коэффициент шероховатости; R — , где R— гидравлический радиус, d—диаметр дрены. Напор в области водосборного колодца Нд — 90 м. На правой границе области задавалось постоянное значение концентрации. Начальное распределение концентраций было нулевым; значения коэффициентов дисперсии щ = 0,5 м и at — 0.1 м.
Для сравнения были проведены расчеты задачи, когда напор на дрене задавался постоянным по всей длине. В этом случае решалась только У Горизонтальная проекция сетки. задача фильтрации, для которой область дренажа представлялась внутренней границей с постоянным значением напора на ней. На рис.33 приведены сравнительные результаты расчётов в виде изолиний напоров в горизонтальном разрезе по координате, соответствующей расположению лучей в плоскости сечения (z = 9,9 м), и в вертикальном сечении плоскостью, проходящей через дрены, параллельные оси X и перпендикулярные оси Z по координате у — 149,9 м (слева — распределение Я, когда течение в трубах отсутствует, справа — с учётом течения).
Учет движения жидкости в дренах меняет фильтрационную картину в области течения. На графике (рис.34) показано распределение напора по длине одной из труб от концевого участка до водосборного колодца. Видно, что распределение и носит неравномерный характер. Дебит скважины действующего лучевого дренажа Qr = 10100 м /сут. При расчётах в предположении постоянного напора в дренах было получено Qr — 11800 м /сут. Все результаты приведены для момента времени, когда режим течения приближался к установившемуся.
Изолинии концентрации в горизонтальном и вертикальном сечениях на различные моменты времени показаны на рис. 35. Форма следа загрязнения определяется гидродинамической картиной течения. Последующие расчёты показали, что при данных условиях задачи лучевой дренаж препятствует дальнейшему распространению примеси.
