Метод конечных элементов с использованием ограниченной вспомогательной переменной для решения уравнения Ричардса

Основная информация

• ID статьи: 2510.13012
• Название: A finite element method using a bounded auxiliary variable for solving the Richards equation
• Авторы: Abderrahmane Benfanich (Университет Оттавы), Yves Bourgault (Университет Оттавы), Abdelaziz Beljadid (Университет Оттавы и Политехнический университет Мохаммеда VI)
• Классификация: math.NA cs.NA physics.comp-ph
• Дата публикации: 14 октября 2025 г. (препринт arXiv)
• Ссылка на статью: https://arxiv.org/abs/2510.13012

Аннотация

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

Исследовательский контекст и мотивация

Определение проблемы

Уравнение Ричардса описывает ненасыщенный поток в пористых средах и является фундаментальной моделью в гидрогеологии, инженерной экологии и почвоведении. Уравнение имеет несколько форм:

1. Форма напорного напряжения (Ψ): применима для полностью насыщенного и частично ненасыщенного потока, но численно нестабильна при S→0
2. Форма степени насыщения (S): лучше работает в сухих почвах, но имеет численные трудности в насыщенной области (S=1)
3. Смешанная форма (S,Ψ): улучшает свойства сохранения массы, но увеличивает сложность

Ограничения существующих методов

1. Методы переключения основных переменных: могут давать нефизические решения на границе раздела насыщенной/ненасыщенной области
2. Техники регуляризации: требуют выбора дополнительных параметров моделирования
3. Традиционные методы конечных элементов: имеют численные трудности при работе с полностью сухими (S=0) и полностью насыщенными (S=1) областями

Исследовательская мотивация

Разработка унифицированного численного метода, способного:

• Избежать переключения переменных и регуляризации
• Стабильно обрабатывать сухие и насыщенные области
• Сохранять хорошие свойства сохранения массы
• Применяться к неоднородным пористым средам

Основные вклады

1. Предложена новая u-форма: введена ограниченная вспомогательная переменная u для исключения неограниченных членов в слабой формулировке
2. Разработана стабильная численная схема: сочетание полуявной временной дискретизации и итераций Ньютона
3. Реализован метод декомпозиции области: обработка неоднородных сред, таких как слоистые почвы
4. Проверена эффективность метода: множественные численные эксперименты подтверждают стабильность и точность метода

Описание метода

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

Решение уравнения Ричардса: $\frac{\partial\theta}{\partial t} + \nabla \cdot q = 0$

где θ — объемное содержание воды, q — поток воды согласно закону Дарси-Букингема: $q = -K_s(x)K_r(x,S)\nabla(\Psi(x,S) + z)$

Конструкция новой переменной u

Для всех моделей производная функции Леверетта может быть записана как: $J'(S) = C(x)S^{-a}(1-S^c)^{-b}$

Определяется вспомогательная переменная: $u = \int_0^S (1-s^c)^{-b} ds$

Представление через неполную бета-функцию: $u = \frac{1}{c}B(S^c, \frac{1}{c}, 1-b)$

Уравнение Ричардса в u-форме

В однородной среде u-форма имеет вид: $\phi \frac{\partial S(u)}{\partial t} - \nabla \cdot \left(K_sK_r(S(u))\left(h_{cap}CS(u)^{-a}\nabla u + e_z\right)\right) = 0$

Вариационная формулировка

Найти u ∈ L²(I;H¹(Ω)) такую, что: $\int_\Omega \phi \frac{\partial S(u)}{\partial t} v dx + \int_\Omega K_sK_r(S(u))\left(h_{cap}S(u)^{-a}\nabla u + e_z\right) \cdot \nabla v dx + \text{граничные члены} = 0$

Временная дискретизация

Используется полуявная схема Эйлера: $\int_\Omega \phi \frac{S(u_h^{n+1}) - S(u_h^n)}{\Delta t} v_h dx + \int_\Omega K_sK_r(S(u_h^n))\left(h_{cap}S(u_h^n)^{-a}\nabla u_h^{n+1} + e_z\right) \cdot \nabla v_h dx = 0$

Итерация Ньютона

Для нелинейного члена S(u^{n+1,m+1}) используется разложение Тейлора первого порядка: $S(u^{n+1,m+1}) = S(u^{n+1,m}) + \frac{\partial S}{\partial u}(u^{n+1,m})(u^{n+1,m+1} - u^{n+1,m}) + O((u^{n+1,m+1} - u^{n+1,m})^2)$

Метод декомпозиции области

Для слоистых почв используется неперекрывающийся метод Шварца:

• Независимое решение уравнения Ричардса в каждой подобласти
• Связь через условия на границе раздела: непрерывность напорного напряжения и сохранение потока
• Итеративное решение с использованием условий передачи типа Робина

Экспериментальная установка

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

• Программное обеспечение: FreeFEM++ и FENICSx
• Решатели линейных систем: UMFPACK, PETSc, MUMPS
• Сетка: равномерная треугольная сетка, линейные конечные элементы (k=1)

Критерии оценки

Ошибка L²: $E_{h,\Delta t}(T)^2 = \|u(\cdot,T) - u_{h,\Delta t}(\cdot,T)\|_{L^2(\Omega)}^2$

Порядок сходимости:

• Пространственный порядок сходимости: $p = \log_2(E_{h,\Delta t}(T)/E_{h/2,\Delta t}(T))$
• Временной порядок сходимости: $q = \log_2(E_{h,\Delta t}(T)/E_{h,\Delta t/2}(T))$

Тестовые случаи

1. 1D полностью ненасыщенный волокнистый образец: проверка стабильности в сухих условиях
2. 1D смешанная насыщенно-ненасыщенная область: тестирование способности обработки разрывных начальных условий
3. Производная решение для проверки: вычисление порядков сходимости
4. 2D полностью насыщенный поток: сравнение с программным обеспечением Hydrus
5. Задача слоистой почвы: проверка метода декомпозиции области
6. Неоднородная среда с включением: тестирование обработки сложной геометрии

Результаты экспериментов

Основные результаты

Тест волокнистого образца 1D

• Сетка: 500 равномерных узлов сетки
• Временной шаг: Δt = 1 с
• Результаты: результаты u-формы и S-формы полностью совпадают, хорошо согласуются с экспериментальными данными

Смешанная насыщенно-ненасыщенная область

• Сетка: 5000 равномерных узлов сетки
• Временной шаг: Δt = 10⁻⁵ часов
• Сходимость Ньютона: в среднем 1-5 итераций
• Стабильность: решение сохраняет монотонность, отсутствуют численные колебания

Анализ сходимости производного решения

Пространственная сходимость (фиксированный Δt = 5×10⁻⁵):

Временная сходимость (фиксированный h = 1×10⁻³):

Сравнение 2D насыщенного потока

• Сетка: 95,142 треугольных элемента, 47,972 вершины
• Временной шаг: Δt = 0.001 часов
• Результаты: высокое совпадение с результатами программного обеспечения Hydrus

Тест слоистой почвы

• Сетка: 89,306 треугольных элементов, 45,064 вершины
• Временной шаг: Δt = 0.05 часов
• Сходимость: допуск 10⁻³, максимум 10 итераций, λ = 25
• Результаты: совпадают с результатами из литературы 20

Сравнение методов

В тесте неоднородной среды с включением:

• u-форма: остается стабильной при ε=0, в среднем 4 итерации
• Ψ-форма: требует >1000 итераций при ε=10⁻⁵, остаток допуска все еще 0.04

Связанные работы

Основные направления исследований

1. Методы переключения переменных: Diersch & Perrochet (1999), Forsyth et al. (1995)
2. Техники регуляризации: Schweizer (2007), Pop & Schweizer (2011)
3. Смешанные формы: Celia et al. (1990), Maina & Ackerer (2017)
4. Декомпозиция области: Lions (1990), Manzini & Ferraris (2004)

Преимущества данной работы

1. Унифицированная обработка: без необходимости переключения переменных или регуляризации
2. Численная стабильность: стабильна на всем диапазоне S∈0,1
3. Физическая согласованность: сохраняет непрерывность напорного напряжения
4. Вычислительная эффективность: быстрая сходимость Ньютона

Выводы и обсуждение

Основные выводы

1. Эффективность u-формы: успешно исключены неограниченные члены традиционных форм
2. Численная стабильность: сохраняется стабильность в полностью сухих и насыщенных областях
3. Производительность сходимости: пространственная сходимость второго порядка, временная сходимость первого порядка
4. Практичность: применима к сложным задачам в неоднородных средах

Ограничения

1. Зависимость от модели: требует выполнения условия $\lim_{S\to 0} K(S)S^{-a} < \infty$
2. Нелинейное решение: требует итераций Ньютона в насыщенной области
3. Обработка границ раздела: неоднородные среды требуют техники декомпозиции области
4. Чувствительность параметров: параметр Робина λ в декомпозиции области требует настройки

Направления будущих исследований

1. Теоретический анализ: строгий численный анализ и доказательства сходимости
2. Многофизическое связывание: расширение на перенос растворенных веществ и тепла
3. Адаптивные техники: разработка адаптивной сетки и адаптивного временного шага
4. Параллельные вычисления: крупномасштабная параллельная реализация

Глубокая оценка

Преимущества

1. Высокая инновационность: конструкция u-переменной остроумна, имеет прочную теоретическую основу
2. Численная стабильность: избегаются численные трудности традиционных методов
3. Полные эксперименты: охватывают случаи от 1D до 2D, от однородных до неоднородных сред
4. Инженерная практичность: сравнение с коммерческим программным обеспечением Hydrus повышает достоверность

Недостатки

1. Недостаточный теоретический анализ: отсутствуют строгие доказательства сходимости и стабильности
2. Недостаточный анализ вычислительной сложности: детальный анализ вычислительных затрат итераций Ньютона недостаточен
3. Выбор параметров: выбор параметра Робина λ лишен теоретического руководства
4. Отсутствие трехмерной проверки: отсутствует проверка на трехмерных задачах

Влияние

1. Академическая ценность: предоставляет новый подход к численному решению уравнения Ричардса
2. Практическая ценность: имеет широкие перспективы применения в гидрогеологии и почвоведении
3. Воспроизводимость: предоставляет подробные детали реализации, облегчающие воспроизведение

Сценарии применения

1. Поток в пористых средах: подземный поток, перемещение почвенной влаги
2. Инженерные приложения: фильтрация через плотины, анализ устойчивости склонов
3. Экологические задачи: миграция загрязняющих веществ, восстановление подземных вод
4. Материаловедение: фильтрация жидкости в волокнистых материалах

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

Ключевые работы

1. Richards, L.A. (1931). Capillary conduction of liquids through porous mediums. Physics 1(5), 318-333.
2. Celia, M.A., et al. (1990). A general mass-conservative numerical solution for the unsaturated flow equation. Water Resour. Res. 26(7), 1483-1496.
3. van Genuchten, M.T. (1980). A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44(5), 892-898.
4. Manzini, G., Ferraris, S. (2004). Mass-conservative finite volume methods on 2-d unstructured grids for the Richards' equation. Adv. Water Resour. 27, 1199-1215.

Общая оценка: Это статья с важным инновационным значением в области численного решения уравнения Ричардса. Путем введения ограниченной вспомогательной переменной u остроумно решаются численные трудности традиционных методов при работе с сухими и насыщенными областями. Метод имеет хорошую теоретическую основу и практическую ценность, численные эксперименты полностью подтверждают эффективность метода. Хотя теоретический анализ требует дальнейшего укрепления, в целом это высококачественная исследовательская работа.