Материал: Hydrogeodynamics101

3. Исходные представления о схемах численного

моделирования нестационарной фильтрации на ЭВМ

Рассмотрим простейшее уравнение одномерной филь­трации (4.1) в конечно-разностной форме:

.Д?₊,-2Я*+Я*_

A t “ (Ал)² ' (4.75)

Пусть для исходного уравнения заданы граничные ус­ловия: Я (0, t) = Я = const, Н (L , f) ^т Н = const и начальное условие И (х, 0) = Й(х) (рис. 4.8). На выбран­ной конечно-разностной сетке краевые условия запишут­ся в виде

Н^к_я=Н_к Н°=П, (для ( > 0 ) ,

где нижний индекс i * 0 отвечает левой границе, a i-M — правой границе (М - Ы А х); верхний индекс к * 0 отве­чает начальному моменту времени.

Перепишем уравнение (4.75) в виде

Я?^+| = Я,*+ ——,(я*,, -2Я,*+ Я* .') ,

¹ (Дх)²1 ^{1+1 1} '7 (4.75а)

где в правой части стоят лишь значения напоров на £-ом временном слое. Положим к=0 и / = 1, тогда из уравнения

75. получаем формулу для расчета напора Н\ в первой узловой точке (т.е. для х * 1 -Ах) на первом временном слое (т.е. на момент t - 1 -A t):

Давая далее индексу i последующие значения (i = 2, 3, ..., М — 1), аналогично определяем все остальные значения Н. на первом временном слое.

Теперь положим в формуле (4.75а) к = 1 и переходим к расчету значений яДт втором временном слое (т.е. для t = 2-Л/), подставляя для этого в правую часть равенства (4.75а) уже известные нам значения Я.¹, определенные для первого временного слоя, и т.д. Так последовательно вы­полняется расчет до последнего N-го временного слоя, Отвечающего конечному расчетному моменту t_N(N'At =

Рассмотренная конечно-разностная схема, описывае­мая уравнениями (4.75), называется явной. Она позво­ляет выразить в явном виде неизвестное значение напора на расчетном временном слое через известные его значе­ния, рассчитанные на предыдущих временных слоях. Это

д^Н

оказывается возможным потому, что производную ₂в

дх

правой части исходного уравнения (4.75) мы выразили через значения напоров, отвечающие началу расчетного временного интервала (т.е. через напоры, отвечающие верхнему положению пьезометрической кривой на рас­смотренной на рис. 4.7а схеме).

На самом же деле, пьезометрическая кривая принима­ет за расчетный интервал времени A t множество последо­вательных положений от Я*до Я^⁺¹, и поэтому логиче­ские соображения никак не препятствуют и другому ва­рианту записи исходного уравнения:

ггЛ+1 тгк гтЛ+1 Л гг^+1 I гг^+1

где пространственная производная в правой части выра­жена через значения напоров, отвечающие концу расчет­ного временного интервала (т.е. через напоры, отвечаю­щие нижнему положению пьезометрической кривой на рис. 4.7а).

ЗАМЕЧАНИЕ. Обратите внимание, что аналогичное уравнение использовано при обосновании схемы Либмана (см. раздел 4.3.2).

Однако, в отличие от уравнения (4.75), уравнение

75. не позволяет явно выразить искомую величину

tff⁺¹ через значения предыдущего временного слоя:

определение #*⁺¹ становится возможным лишь после того, как мы запишем выражения вира (4.76) для всех узлов / = 1,2,..., М — 1 на к +1 -ом слое и решим получен­ную систему уравнений. С этой точки зрения конечно­разностная схема, описываемая уравнениями вида (4.76), получила название неявной.

^li

циалы i/*^—1 в предшествующий момент (£-1). Решение задачи на такой модели проводится от шага к шагу:

[Т ] на концы сопротивлений R_t подаются потенциа­лы t/P, отвечающие заданным начальным значениям на­пора Я.⁰, и в узловых точках снимаются потенциалы Я/, отвечающие напорам Я/для первого расчетного момента времени t_x = 1 A t\

_2 J значения Uj подаются на концы временных со­противлений и в узловых точках снимаются потенциалы Uf и т.д.

Аналогично могут решаться и двумерные плановые задачи, тогда каждый узел сетки будет содержать пять сопротивлений (в уравнении (4.68а) /1 = 4).

Т77УТТТ71^\У1^^Т7Т77%-ПМ /7~/7х

L

Рис. 4.8. Схема пространственной разбивки области фильтрации

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

Наряду с этим используются и ЛС-сетки, в которых время остается непрерывным: емкость водоносной систе­мы реализуется разрядкой во времени электрических ем­костей-конденсаторов (С), подключаемых в узловые точ­ки (вместо временных сопротивлений). В этом случае уравнение электрического тока имеет вид:

" ^иг^и _С»У

l = i ^Ri ^st_M' (4.72)

где, в отличие от формулы (4.70а), в правой части сохра­няется непрерывная форма записи производной. Так как этому случаю соответствует представление правой части

уравнения (4.68) в виде/i а)~~, то для подобия упомяну-

д t

тых уравнений фильтрации и нестационарного электри­ческого тока необходимо потребовать, чтобы

ц* а) С, t — a_ft_M, (4.73)

где t_M — модельное время (замеряемое время разрядки конденсатора).

Масштабные коэффициенты a_t и а_ф должны быть, очевидно, связаны дополнительным критерием подобия:

^at~^a[i ^аФ" (4.74)

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

Простота, доступность и физическая осязаемость се­точных электрических моделей наряду с вполне удовлет­ворительной надежностью и точностью привели в свое

время к исключительно широкому и эффективному их внедрению в практику гидрогеологических исследований (примеры такого рода приведены в разделе 8.3); методи­ческие аспекты моделирования более подробно рассмот­рены, например, в работах [7,14]. Вместе с тем, в послед­ние десятилетия в гидрогеологии наиболее широко стало использоваться математическое моделирование на ЭВМ [19,20,40,48].

Понятно, что и для человека, и для ЭВМ гораздо проще последовательно провести М однотипных расчетов по формуле (4.75а), чем решать систему из М уравнений. Поэтому, казалось бы, логически ясно, что всегда разум­нее считать по явной схеме, нежели по неявной. Однако с конечно-разностными схемами дело обстоит отнюдь не так просто. Они обладают специфическими свойствами, из-за недоучета которых подобные, внешне логичные рассуждения могут оказаться неприемлемыми. Для при­мера на рис. 4.9 приведены в схематизированном виде результаты расчетов одной и той же фундаментальной (см. раздел 4.1.1) задачи по явной и неявной схемам. По мере роста времени (числа временных шагов) неявное решение все теснее сближается с точным (построенным по аналитической зависимости (4.12)). Между тем, явное решение при том же, достаточно большом, временном шаге At ведет себя весьма неестественно: оно начинает постепенно «раскачиваться» и приводит к физически не-

* к.

При этом известными, кроме величин Н- с предыдущего временного слоя, яв­ляются граничные значения реальным результатам (например, оно дает положение уровня воды ниже подошвы водоносного горизонта). Что­бы понять такое поведение решения, нужно вспомнить, что разностные представления аппроксимируют произ­водные, входящие в исходное уравнение, приближенно и, следовательно, на каждом шаге вычислений в значения искомой функции (напора) вносятся какие-то погрешно­сти. Если в процессе вычислений по мере роста числа операций (в нашем случае — числа временных шагов) эти погрешности постепенно подавляют, гасят друг друга, то конечно-разностная схема является устойчивой и не может приводить к результатам, подобным кривой 2 на рис. 4.9. В противном же случае, когда идет накопление погрешностей в процессе счета, схема называется не­устойчивой. Ясно, что вести расчет можно только по устойчивым схемам.

Так вот, оказывается, что неявные схемы (в частности

76. ) устойчивы всегда (отсюда и отмеченная выше на­дежность схемы Либмана), в то время как для устойчивости явных схем необходимо вводить ограничения на шаг по времени A t. Например, схема (4.75) устойчива, если

a At 1

(4.77)

Рис. 4.9. Схематизированное представление результатов расчета фундаментальной фильтрационной задачи:

1 - данные расчета по неявной схеме; 2-то же, по явной схеме; 3 - точное аналитическое решение

причем при соблюдении этого условия результаты счета по явной схеме (4.75) сходятся к точному решению даже быстрее, чем при использовании неявной схемы (4.76) при тех же параметрах сетки Аде и A t.

Обратите внимание, что в условии (4.77) в знаменателе стоит

величина Аде?. Следовательно, увеличение дробности пространст­венной разбивки (уменьшение Ах) не обязательно приводит к росту точности вычислительной схемы (это еще раз показывает, что логику численных методов нельзя уяснить, отталкиваясь от одного лишь здравого смысла).

Выполнение условия (4.77) часто требует, однако, резкого увеличения числа временных шагов и объема расчетных операций в целом. Именно поэтому в практике моделирования на ЭВМ основное развитие получили не­явные схемы, а также смешанные - явно-неявные.

Простейшим примером явно-неявной схемы может служить следующее конечно-разностное представление уравнения (4.1), обобщающее схемы (4.75) и (4.76):

^г к I и к

н^к+1-н¹ *

- = а

— 2 Н?+ 7//+J

(1-е)

+

A t

(Ах)²

т,к+1 ^Hi+1

1. Особенности задач, связанных

с интерпретацией опытно-фильтрационных исследований

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

водности). Для некоторых задач оказывается необходи­мой также оценка инфильтрационного питания и пара­метров перетекания. В качестве специальной характери­стики фильтрационного сопротивления ложа реки или водоема выступает параметр AL.

ЗАДАНИЕ. На основе ранее изученных разделов курса, уточни­те для себя следующие вопросы: 1) от каких факторов зависят воз­можные изменения коэффициента фильтрации (проницаемости), в том числе в трещиноватых и глинистых породах (см. раздел 1.5); 2) как связаны между собой коэффициенты фильтрации и водопрово- димости (см. раздел 2.3); 3) от каких факторов зависят возможные изменения коэффициентов гравитационной и упругой водоотдачи (см. раздел 1.4); 4) как связаны эти величины соответственно с коэф­фициентами уровнепроводности и пьезопроводности (см. раздел 2.3); 5) каким параметром характеризуется интенсивность стацио­нарного перетекания через относительный водоупор (см. раздел 2.3); б) что представляет собой параметр AL (см. раздел 3.4); 7) почему для решения широкого круга прогнозных задач нет необходимости в данных о параметре инфильтрационного питания (см. раздел 3.3).

До сих пор мы предполагали, что эти параметры нам заданы. Теперь, однако, настало время выяснить, как же получают числовые значения фильтрационных парамет­ров для конкретного водоносного горизонта, как его иден­тифицируют (т.е. узнают его фильтрационную сущность). Так как сам горизонт не может быть изучен досконально «изнутри», то его приходится рассматривать как систему типа «черный ящик», о которой можно судить лишь по ее реак­ции на какие-то внешние воздействия. Коль скоро мы говорим об оценке фильтрационных параметров, то эти воз­действия должны, естественно, сводиться к некоторым филь­трационным возмущениям, изменяющим скорости и напоры в изучаемом пласте. Если такие возмущения специально со­здаются искусственно — путем откачки (выпуска) воды из пласта или ее закачки (нагнетания, налива), то их назы­вают опытно-фильтрационными опробованиями (ОФО). Если же для определения фильтрационных параметров используют данные возмущений, сопутствующих раооте какого-либо водозабора, дренажа и т.п., не имеющих своей основной целью оценку параметров, то говорят об опытно-фильтрационных наблюдениях (ОФН).

Общим в обоих указанных случаях, объединяемых общим термином «опытно-фильтрационные работы» (ОФР), является то, что неизвестными в решаемой задаче служат фильтрационные параметры, а заданными — на­поры в отдельных точках пласта и расходы потока на некоторых участках области фильтрации. В математиче­ском плане мы отыскиваем коэффициенты или свободные члены дифференциальных уравнений по заданным част­ным значениям функции и (или) ее производных. Такого рода задачи относятся к классу так называемых обрат­ных задач.В решавшихся нами ранее прямых зада­чах коэффициенты и свободные члены считались извест­ными. Заметим, что к обратным принято относить также задачи, в которых неизвестны те или иные краевые усло­вия.

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

- 2 tff⁺¹ + Н^к+\

+ а

(Д xf

(4.78)

где а — весовой коэффициент (0 < а <1).

Схема (4.78) при а = 1 оказывается неявной, при а * О — явной, а при промежуточных значениях а —смешан­ной, явно-неявной. Теория и численные эксперименты показывают, что такие смешанные схемы могут вобрать в себя достоинства как явных, так и неявных схем. В част­ности, при а > 0,5 явно-неявная схема (4.78) всегда ус­тойчива (подобно неявной схеме (4.76)), но за счет нали­чия явного члена (при а & 1) она сходится к точному решению быстрее, чем схема (4.76) (иначе говоря, для достижения той же точности она требует меньшего числа временных шагов, а в конечном счете — наименьшего машинного времени).

ЗАМЕЧАНИЕ. Пространственная разбивка для рассмотренных здесь конечно-разностных уравнений назначается исходя из выбора элементов длины Ах, что при расчетах планово-неоднородных пла­стов может приводить к большим численным погрешностям. Поэто­му на практике проводят разбивку области фильтрации по фильтра­ционным сопротивлениям, аналогично схеме Либмана (см. раздел 4.3.2). Соответствующие расчетные схемы (см., нарример, уравне­ние (4.68)), получившие название консервативных , позволяют ве­сти счет фильтрационных задач при грубых пространственных сет­ках, т.е. при сравнительно небольшом числе пространственных узлов М. На практике большинство задач решается при числе М, не пре- вышающм 2500-3000.

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

Контрольные вопросы

[Т] Опишите физическую картину осушения напорного пла­ста при быстром снижении уровней на одной из его границ до кровли пласта. Что принципиально меняется при снижении уровня до подо­швы пласта?

|~2~| В чем смысл и значение краевых условий фильтрации? Записать и объяснить краевые условия для фундаментальной задачи (о мгновенном снижении напоров подземных вод на границе пласта).

[~3] Для какой физической ситуации получено фундаменталь­ное решение нестационарной фильтрации? Можно ли им воспользо­ваться для описания безнапорной фильтрации? Как расширяются возможности применения этого решения благодаря принципу супер­позиции и методу недеформируемых лент тока?