Моделирование фильтрации в насыщенном влагой грунте

Актуальность проблемы в строительстве

Значения скоростей фильтрации воды в грунте во многом определяют конструкцию и материалы, применяемые при возведении фундаментов. Воздействие фильтрации в грунте на себе испытывают также стены подвалов, линий подземных коммуникаций и других заглубленных в грунт частей сооружений. Устойчивость дна и берегов водохранилищ и каналов тоже зависит от фильтрации в прибрежном грунте. Учет фильтрации повышает качество моделирования иных физических процессов в грунте. Например, при распределении тепловых полей, конвективный перенос тепла обусловлен наличием фильтрации.

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

УРАВНЕНИЕ ДАРСИ

Перейдем к детальному рассмотрению уравнения Дарси. Данное уравнение названо в честь французского инженера-гидравлика Анри Филибера Гаспара Дарси и является уравнением в частных производных второго порядка:

(1)   \begin{equation*}\bigtriangledown\ * (-K\bigtriangledown H)=Q_s\end{equation*}

где Q_s - источник жидкости, [Q_s]=c^{-1}, K – коэффициент фильтрации, [K] = м/с; H - гидравлический напор, [H] = м; \bigtriangledown=\frac{\partial} {\partial x}\vec{i}+\frac{\partial}{\partial y}\vec{j}+\frac{\partial}{\partial z}\vec{k} - оператор Гамильтона.

Вектор скорости фильтрации v, исходя из закона Дарси, равен

(2)   \begin{equation*}v= -K\bigtriangledown H\end{equation*}

Для уравнения (1) задают следующие граничные условия:

1) граничное условие первого рода H=H_{0} , здесь H_{0} – заданный гидравлический напор;

2) граничное условие второго рода \vec{n}K\bigtriangledown H =q_{0}, здесь q_{0} – поток жидкости через граничную поверхность, [q_{0}] =м/с, а \vec{n} – единичная внешняя нормаль к поверхности.

Если на всей границе задано граничное условие первого рода, то имеем задачу Дирихле для уравнения Дарси, если же на всей границе задано граничное условие второго рода, то имеем задачу Неймана для уравнения Дарси. В общем случае будем считать, что на границе заданы граничные условия первого и второго рода. Таким образом, далее рассматривается задача со смешенными граничными условиями.

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

ЧИСЛЕННОЕ РЕШЕНИЕ УРАВНЕНИЯ ДАРСИ:

Пусть расчетная область является параллелепипедом и задана расчетная сетка узлов \omega{{(x_{i}, y_{j}, z_{k})}}:i= 1,\ldots, N_{x}; j= 1,\ldots, N_{y}; k= 1,\ldots, N_{z}. Далее под границей будем понимать как границу между материалом (грунтом) и внешней средой, так и границу расчетной области. Каждому узлу p={(x_{i}, y_{j}, z_{k})} расчетной сетки сопоставим вектор d^{p}={(d_{x}^{p}, d_{y}^{p}, d_{z}^{p})}, который определяется следующим образом:

1) p=(0,0,0), Если узел p не принадлежит границе;

2) если же узел p принадлежит границе, то d^{p}=(sign(n_{x}), sign(n_{y}), sign(n_{z}))}, где n_{x}, n_{y}, n_{z} – компоненты вектора нормали \vec{n} в узле p.

Перейдем к построению системы уравнений для численного решения задачи (1). Под H_{i,j,k} обозначим приближенное значение напора H (x,y,z) в узле p=(x_{i},y_{j},z_{k}). Для нахождения неизвестных значений H_{i,j,k} составим систему линейных алгебраических уравнений. Рассмотрим два случая в зависимости от геометрического положения узла p. Если узел p=(x_{i},y_{j},z_{k}) не принадлежит границе, то уравнение (1) аппроксимируется следующим образом:

Если узел p находится на границе, то возможны две ситуации в зависимости от рода граничного условия, заданного в узле p. Пусть в узле p заданно граничное условие первого рода, в этом случае следует использовать уравнение H_{i,j,k}= H_{0}(x_{i},y_{j},z_{k}). Если же в узле p заданно граничное условие второго рода, то в этом случае используют уравнение

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

Ah=b,                                                                       (4)

где матрица коэффициентов A имеет размеры N\times N,N=N_{x}\times N_{y}\times N_{z}.

Отметим, что в общем случае порядок N системы (4) может быть равен нескольким миллионам, таким образом, систему (4) приходится решать итерационными методами. Например, используются метод простых итераций, метод Зейделя, метод сопряженных градиентов, стабилизированный метод бисопряжённых градиентов и др. Для ускорения решения системы уравнений (4) применяется параллельное вычисление на многоядерных CPU, а также на GPU (в частности, используется технология CUDA).

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

При одновременном моделировании тепловых полей и фильтрации жидкости в грунте приходится учитывать эффект промерзания грунта. При промерзании грунта поры заполняются льдом, вследствие чего уменьшается коэффициент фильтрации. Полагаем, что коэффициент фильтрации K(x,y,z) для промерзшего грунта равен K(x,y,z)= (1-(w(x,y,z))K_{0}(x,y,z), где w(x,y,z) – доля льда в грунте в точке с координатами (x,y,z), а K_{0}(x,y,z) – коэффициент фильтрации растепленного грунта в той же точке.

ПРИМЕР РАСЧЁТА СКОРОСТЕЙ ФИЛЬТРАЦИИ В ПРОГРАММЕ FROST 3D UNIVERSAL

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

Пусть рассматриваемый участок грунта имеет форму параллелепипеда размеров 15м\times10м\times10м, обозначим его P (см. рисунок 1). В начальный момент времени грунт в параллелепипеде находится в талом состоянии. На левой грани параллелепипеда P (грань x=0м) имеет место гидравлический напор равный H_{x=0}= 0.5м, а на правой грани (x=15м) гидравлический напор равен H_{x=15}=0м. Под действием этого градиента напоров в параллелепипеде будет происходить фильтрация грунтовых вод со скоростью:

v_{x}= K \frac {\Delta H}{\Delta x} = 0.001 \frac {0.5-0}{15-0} = 3.33 \times 10^{-4} м/с

Здесь K – коэффициент фильтрации.

Предположим, что в параллелепипеде P имеется область грунта С цилиндрической формы с радиусом основания (основание параллельно плоскости O_{xy}) r=4м и высотой h=8м, центр цилиндра имеет координаты (5м,5м,5м). В этом цилиндре в начальный момент времени грунт находится в мерзлом состоянии. Поскольку грунт в области C находится в мерзлом состоянии, то скорость фильтрации в этой области будет равна нулю.

Определим с учетом заданного поля скоростей фильтрации распределение тепловых полей в рассматриваемом участке грунта при значении температуры на левой грани параллелепипеда P (грань x=0 м) T_{x=0} =+2 C_{0}, а на правой грани (x=15м) T_{x=15} =-2 C_{0}. Заданные в расчете теплофизические свойства грунта для областей P и С и граничные условия приведены в таблицах 1 и 2. В результате расчета получаем распределение температур, представленное на рисунках 2 и 3.

Модель оценки влияния фильтрации на теплопередачу

Рисунок 1 – Расчетная область

Таблица 1. – Теплофизические свойства

Параллелепипед Цилиндр
Объемная теплоемкость талого грунта djm3xk
2140000 2140000
Объемная теплоемкость мерзлого грунта djm3xk
2140000 2140000
Теплопроводность талого грунта vtmxk
1.8 1.8
Теплопроводность мерзлого грунта vtmxk
1.8 1.8
Суммарная весовая влажность грунта 1 1
Плотность сухого грунта km3
1000 1000
Количество незамерзшей воды 01temperaturefazperehod
Температура фазового перехода c
0 0
Коэффициент фильтрации ms
0.001 0.001
Начальная температура c
+4 -4

Таблица 2. – Граничные условия

Граничное условие тепловой задачи Граничное условие фильтрационной задачи
грань x=0 m Теплообмен по Ньютону:
температура – +20C
коэффициент теплообмена –
Гидравлический напор – 0.5 м
грань x=15 m Теплообмен по Ньютону:
температура – -20C
коэффициент теплообмена –
Гидравлический напор – 0.0 м
остальные грани Тепловой поток – Скорость фильтрации – 0 м/с

 

Температурное распределение при статически заданных значениях фильтрации

Рисунок 2 – Распределение тепловых полей в градусах Цельсия на глубине 5 м через 10 суток для заданной скорости фильтрации

Температурные поля при заданных статических значениях фильтрации

Рисунок 3 – Изолинии температур в градусах Цельсия на глубине 5 м через 10 суток для заданной скорости фильтрации

В действительности скорость фильтрации в параллелепипеде P из-за наличия области C не просто равна 3.33x10-4 м/с, а имеет сложный вид. Найти это распределение скорости фильтрации можно, решив уравнение Дарси. Так, из решения уравнения Дарси в программе Frost 3D Universal получим распределение полей скорости фильтрации (рисунок 4 и 5) и проанализируем, как изменится тепловое поле в случае учета более корректных значений скоростей фильтрации грунтовых вод. Полученные результаты расчетов теплового поля приведены на рисунках 6 и 7.

Фильтрационные потоки по уравнению Дарси

Рисунок 4 – Проекция скорости фильтрации на ось X (в микрометрах в секунду), вычисленная на основании уравнения Дарси

Скорости фильтрации при передачи тепла

Рисунок 5 – Проекция скорости фильтрации на ось Y (в микрометрах в секунду), вычисленная на основании уравнения Дарси

Температурное распределение при правильно рассчитанном значении фильтрации

Рисунок 6 – Распределение тепловых полей в градусах Цельсия на глубине 5 м через 10 суток для вычисленной скорости фильтрации на основании уравнения Дарси

Температурные поля при рассчитанных значениях фильтрации по уравнению Дарси

Рисунок 7 – Изолинии температур в градусах Цельсия на глубине 5 м через 10 суток для вычисленной скорости фильтрации на основании уравнения Дарси

Сравнив результаты вычислений тепловых полей с расчетом скорости фильтрации на основании уравнения Дарси (рисунки 6 и 7) и без расчета (рисунки 2 и 3) видим, что имеется существенное отличие в распределении температур. Особенно критичным является тот факт, что разница температуры в одних и тех же узлах, для этих расчетов достигает по абсолютной величине более 10 0C.

ВЫВОД

В приведенной статье рассмотрен подход вычисления скоростей фильт-рации в грунте, основанный на численном решении уравнения Дарси. Данное уравнение правомерно применять в случае, когда грунт насыщен жидкостью.

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

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *


*

Можно использовать следующие HTML-теги и атрибуты: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>