Применение машинных методов в гидрогеологических расчетах

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

где h — уровень подземных вод в точках пласта с координатами х, у на момент времени t; T — водопроводимость пласта; u — водоотдача водоносных пород; W — интенсивность питания или расходования подземных вод.

Наибольшее распространение получили конечно-разностные методы численного решения уравнения (154), основанные на приближенной аппроксимации его различными разностными выражениями. Самыми простыми и легко реализуемыми в вычислительных программах ЭВМ являются явные конечно-разностные схемы. Однако они обладают счетной устойчивостью лишь при весьма небольших промежутках времени. Л. Томлинсон и К. Раштон предлагают модификационную явную схему (явный метод переменных направлений), позволяющую несколько сгладить условие устойчивости. Ho в общем, судя по имеющимся источникам, явные конечно-разностные аппроксимации почти не используются на практике.

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

где Ai,j, Bi,j, Di,j и Ei,j — коэффициенты, определяемые в зависимости от величин водопроводимости и размеров рассматриваемого блока (i, j) и соседних с ним блоков:

где At — временной шаг; i, j и k — соответственно номер строки, столбца рассматриваемого узла и номер временного шага.

Выражения (155) представляют систему линейных алгебраических уравнений относительно искомых уровней в узлах сетки. Число таких уравнений, равное числу узлов сетки разбивки рассматриваемой области, при решении практических задач составляет обычно 0,5—1,5 тыс. Поэтому прямое обращение матрицы (156) требует значительного машинного времени, большой памяти машины и, судя по имеющимся работам, практически не используется при решении двумерных гидрогеологических задач.

Одним из основных способов решения неявного уравнения (155) является метод переменных направлений ADIP или ADI (Alternating Direction Implicit Procedure). Суть этого метода заключается в приближенной замене каждого выражения системы (155) уравнениями для двух полушагов по времени. На первом полушаге расчеты ведутся в одном из направлений (например, по оси х) по неявной схеме, а в другом (ось у) — по явной схеме. На втором полушаге порядок расчетов меняется на обратный. Вычислительная схема может быть записана следующим образом:

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

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

Часто используют также итерационные методы решения конечно-разностных уравнений. Решение системы (1,55) достигается при этом путем последовательных приближений. Применяемые в настоящее время основные итерационные методы получили название: итерационный неявный метод переменных направлений (PIPT), линейная последовательная верхняя релаксация (Line Successive Overrelaxation — LSOR) и чисто неявный итерационный метод (Strongly Implicit Procedure — SIP).

В итеративной схеме метода переменных направлений для получения n-го приближения по известному (n—1) (первое приближение задается произвольно) решаются две системы алгебраических уравнений. Сначала определяются промежуточные значения уровня hn-1/2. При этом, так же как и в методе ADIP, используется явная схема расчетов в одном из направлений и неявная — в другом, например:

где w — итерационный параметр: 1 < w < 2.

Окончательное значение уровня на данной итерации и определяется на втором этапе:

Уравнения (160) и (162) численно решаются методом прогонки.

Примеры использования метода переменных направлений можно найти в работах И. Бредехофта и Дж. Пиндера, Р. Кули, Р. Волкера и И. Мухи.

По методу LSOR расчеты также проводятся в два этапа. Сначала по формуле итераций Гаусса — Зейделя, находятся промежуточные значения уровня h:

Затем определяется окончательная величина уровня на данной итерации

Расчет проводится последовательно вдоль каждой строки (столбца). Значения уровня h* при этом определяются по методу прогонки.

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

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

Метод SIP, основанный на факторизации матрицы А, рассмотрен П. Трескоттом и С. Ларсоном, Т. Уинтером, Р. Кули. В последней работе приведены вычислительные программы, реализующие итерационные методы, рассмотренные выше,

П. Трескотт и С. Ларсон приводят сравнительную оценку итерационных методов. Авторы применяют ADIPIT, LSOR и SIP одновременно для решения четырех тестовых задач. Сопоставляя результаты расчетов, авторы приходят к выводу, что лучшим из рассмотренных методов является SIP. Он характеризуется наибольшей скоростью сходимости итерационного процесса, по этому методу получены хорошие результаты решения всех тестовых задач. Недостаток последнего заключается в том, что он требует большей машинной памяти, чем другие. Следующим по эффективности, по мнению авторов, стоит метод LSOR.

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

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

При решении фильтрационных задач обычно применяют два способа реализации конечных элементов: метод Рэлея — Ритца и метод Галеркина.

По методу Рэлея — Ритца решение дифференциального уравнения (154) сводится к отысканию минимума интегрального функционала, который строится на основе вариационной формулы Эйлера. Физически это означает, что находится такая поверхность подземных вод, которая обладает наименьшей потенциальной энергией. Например для уравнения (154) указанный функционал можно записать в виде

где w — исследуемая область фильтрации.

Подставив в это равенство приближенное выражение, дающее дискретное представление уровня, находят величину I. Для достижения его минимума используется условие

где i — номера узлов. Из (168) получается система алгебраических уравнений относительно значений уровня подземных вод в узлах модели области фильтрации. При небольшом числе узлов (менее 100—200) эта система может быть решена методами исключения Гаусса, а при большом числе узлов — итерационными методами, изложенными выше. Такова общая схема решения фильтрационных задач по методу Рэлея — Ритца. Теоретическое обоснование его применительно к гидрогеологическим задачам дается И. Яванделем и П. Уитерснуном, Р. Ченгом, а примеры использования в работах.

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

где h — искомая величина; L(h) — дифференциальное уравнение.

После подстановки базисных функций Ni в уравнение (169) и интегрирования получается система алгебраических уравнений, обращение которой дает приближенное решение задачи. Достоинство метода Галеркина — исключение необходимости вариационной формулировки задачи.

Современное состояние техники моделирования методом конечных элементов освещено Р. Ченгом. Применение метода формулируется как с использованием методики Рэлея — Ритца, так и Галеркина. Даются примеры применения метода для изучения различных проблем фильтрации через пористые среды.

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

Во многих работах рассмотрены преимущества методов конечных элементов и конечных разностей. Отмечено, что метод конечных элементов является, в общем, более мощным вычислительным средством, позволяющим получать результаты для самых сложных гидрогеологических схем. Так, например, при рассмотрении анизотропной среды, в которой главные направления анизотропии меняются в зависимости от направления потока, задача легко решается методом конечных элементов, тогда как ее практически невозможно решить методом конечных разностей. Конечными элементами легко аппроксимируются сложные и криволинейные области и отдельные участки в них. Метод конечных элементов позволяет учесть сложные граничные условия.

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

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

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

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

Методика аналогового моделирования с учетом зарубежного опыта в отечественной литературе хорошо освещена Л. Лукнером, В,М. Шестаковым, И. Д. Жерновым, И.Н. Павловцом, И.И. Kpaшиным и др.

Большая часть электрических моделей, предназначенных для исследования фильтрации подземных вод, строится гидрогеологическими организациями. Известны лишь несколько моделирующих устройств, изготовленных промышленностью: анализатор поля АР-600 (Польша), аналоговые устройства типа «Аналог» (ГДР). Разрабатываются также гибридные вычислительные комплексы, включающие электрические сетки и цифровые вычислительные машины. Однако широкого применения они не получили.

В последние годы во многих странах практические фильтрационные расчеты проводят большей частью с использованием ЭВМ, аналоговые устройства применяются реже. Так, в геологической службе США лабораторией аналогового моделирования, действующей с 1959 г., было создано более 100 гидрогеологических моделей. В последние 10 лет для оценки запасов подземных вод, влияния водоотбора на поверхностные воды и т. д. с применением ЭВМ было изучено более 250 объектов.

Таким образом, в последнее время отмечается качественное изменение методики гидрогеологических расчетов. Если раньше основным методом решения сложных практических задач, связанных с расчетом водозаборов, изучением региональных закономерностей формирования и движения подземных вод было аналоговое моделирование, то теперь те же задачи могут быть успешно решены с применением ЭВМ. Как показывает приведенный обзор, за рубежом интенсивно разрабатываются численные методы и вычислительные программы для количественного исследования фильтрации подземных вод в различных гидрогеологических условиях. Эти методы уже использованы на многих конкретных объектах.