Производственный и научно-исследовательский институт
по инженерным изысканиям в строительстве (ПНИИИС)

Рекомендации
по решению
на ЭВМ задач
геофильтрации
для трехслойного
безнапорно-напорного потока
подземных вод
при подтоплении
городских
территорий

Москва Стройиздат 1984

Рекомендовано к изданию решением секции гидрогеологии и гидрологии Научно-технического совета ПНИИИС от 1 сентября 1982 г.

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

Для инженерно-технических работников изыскательских и научно-исследовательских организаций.

Ил. 6.

Разработаны ПНИИИС Госстроя СССР (кандидаты геол.-минерал. наук B.C. Зильберг, Г.М. Великина, канд. техн. наук В.М. Лившиц)

СОДЕРЖАНИЕ

1. ОБЩИЕ ПОЛОЖЕНИЯ

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

1.2. Основными гидродинамическими особенностями геофильтрации в трехслойной системе являются:

наличие двух плоскоплановых потоков (безнапорного и напорного), разделенных слабопроницаемым прослоем с вертикальной фильтрацией подземных вод и объединенных в единый безнапорный поток при отсутствии слабопроницаемого прослоя или наличии в нем цитологических «окон»;

гидравлическая разнонаправленная взаимосвязь поверхностных и грунтовых вод;

фильтрационная неоднородность водовмещающих отложений;

наличие переменных во времени и по площади источников и стоков, к которым относятся водозаборные скважины, водонесущие коммуникации, дренажные и мелиоративные системы;

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

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

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

Решения их для области фильтрации с неправильной внешней конфигурацией и вышеприведенными гидродинамическими особенностями наиболее рационально осуществлять численными методами с реализацией на ЭВМ. В связи с этим был составлен на алгоритмическом языке «ФОРТРАН-IV» применительно к операционной системе ДОС ЕС-1022 пакет прикладных программ численного решения систем дифференциальных уравнений эллиптического и параболического типов. Их практическая разработка базировалась на решении двух геофильтрационных задач: оценке запасов подземных вод и прогнозе уроненного режима грунтовых вод в связи с самоподтоплением и подпором от водохранилища. Методика решения последней задачи приводится в настоящих Рекомендациях в качестве иллюстрационного примера.

2. ГИДРОГЕОЛОГИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ

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

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

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

2.3. На искусственных внешних плановых границах задаются отток или приток подземных вод с помощью ГУ III рода или ГУ II рода, включая представление их в виде непроницаемого контура.

2.4. Внутренними источниками возмущения являются водоотбор подземных вод скважинами (ГУ II рода) и дренажными системами (ГУ II рода), разгрузка грунтовых вод во внутреннюю гидрографическую сеть (ГУ III рода), фильтрация из нее (ГУ II или III рода), утечки из водонесущих коммуникаций (ГУ II рода), задаваемые в виде соответствующих граничных условий.

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

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

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

3. МАТЕМАТИЧЕСКОЕ ОПИСАНИЕ ГЕОФИЛЬТРАЦИОННОЙ ЗАДАЧИ

3.1. Для математической записи поставленной задачи необходимо в области фильтрации, расположенной внутри криволинейной границы Г, выделить замкнутыми контурами g внутренние источники и стоки. Оставшаяся после их исключения многосвязная область D должна быть разделена на зоны трехслойного D₁ и однослойного D₂ строения. Тогда неустановившийся процесс фильтрации грунтовых и подземных вод в гидравлической постановке, в соответствии с которой в хорошо проницаемых слоях рассматриваются только горизонтальные составляющие скорости фильтрации, а в разделяющем слое - только вертикальные, запишется в виде системы дифференциальных уравнений для зоны D₁, в которой в верхнем горизонте безнапорный режим фильтрации, в нижнем - напорный режим,

                               (1)

и одним уравнением для зоны D₂

                                 (2)

где

h(х, у, t), H(х, у, t) - уровень и напор в верхнем и нижнем водоносных горизонтах соответственно, м;  - уровень водоносного горизонта в однослойной толще, м; Т_V - водопроводимость верхнего (безнапорного) водоносного горизонта, м²/сут; T_N - водопроводимость нижнего (напорного) водоносного горизонта, м²/сут; m(х, у) и m(х, y) - свободная и упругая водоотдачи соответственно; W(x, y, t) - интенсивность инфильтрации, м/сут; (x, y, h) - интенсивность испарения; м/сут; T_VN - водопроводимость однослойной толщи, м/сут; q_п - единичный переток между горизонтами, вычисляемый по формуле

                                  (3)

где H_кр(x, у) - абсолютная отметка кровли напорного водоносного горизонта, м; k₀(x, у) и т₀(х, у) - коэффициент фильтрации и мощность раздельного слоя, м/сут и м.

3.2. Интенсивность испарения  зависит от глубины стояния грунтовых вод z и вычисляется по зависимости

                                                          (4)

где w₀ - интенсивность испарения у поверхности земли (при z = 0); z_k - критическая глубина стояния грунтовых вод, ниже которой уровень грунтовых вод уже не влияет на величину испарения; z = z₀ - h, где z₀ - абсолютная отметка поверхности земли.

3.3. Установившаяся фильтрация подземных вод для рассматриваемой расчетной схемы описывается также уравнениями (1) - (2) при условии равенства нулю левых частей этих уравнений.

3.4. Однозначность решения уравнений (1) - (2) достигается заданием граничных и начальных условий.

Граничные условия, задаваемые на внешнем контуре Г и внутренних g, имеют вид:

первого рода (ГУ I)

                                            (5)

где

Г₁ + g₁ - часть внешнего к некоторые внутренние контуры с заданными ГУ I (рис. 1);

Рис. 1. Фильтрационная схема

1 - 3 - граничные условия: 1 - первого рода H = H(r, t) на контурах; а - внешних линейных Г₁; б - внутренних линейных (g₁); в - внутренних точечных (g₁); 2 - второго рода на внешних контурах Г2 линейных; а - для Q = Q(p, t); б - для Q = 0; внутренних (g₂); в - линейных для Q = Q(р); г - точечных для Q = Q(p, t); 3 - третьего рода на контурах; а - внешних Г₃; б - внутренних линейных g₃; в - внутренних точечных g₃; 4 - граница перехода трехслойной зоны D₁ в однослойную D₂

второго рода (ГУ II)

                                             (6)

где

Г₂ + g₂ - часть внешнего и некоторые внутренние контуры с заданными ГУ II (см. рис. 1);

третьего рода (ГУ III)

                                            (7)

где  - уровень воды в поверхностном водоеме; Н_пр(s) - абсолютная отметка подошвы подрусловых отложений; c(s) - параметр взаимосвязи подземных и поверхностных вод; s - точка, принадлежащая Г₃ + g₃ - части внешнего и некоторым внутренним контурам с заданными ГУ III (см. рис. 1).

Принято, что Г = Г₁ + Г₂ + Г₃; g = g₁ + g₂ + g₃.

3.5. Начальные условия имеют вид

                                                    (8)

3.6. Для однозначности решения уравнений (1) - (2) в условиях стационарной фильтрации подземных вод достаточно задания только граничных условий (5) - (7); при этом все величины, входящие в эти формулы не зависят от времени t.

4. АЛГОРИТМ ЧИСЛЕННОГО РЕШЕНИЯ НА ЭВМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

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

4.2. Для представления дифференциальных уравнений в конечно-разностном виде на область фильтрации D накладывают сетку Д и непрерывной в области D функции H ставят в соответствие сеточную функцию H_ij, такую, чтобы выполнялось равенство

H_ij = Н(x_i, y_j).                                                                (9)

Выбор сетки не однозначен. Наибольшее применение получили прямоугольные неравномерные сетки, которым отдастся предпочтение.

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

4.4. Для составления разностных схем на неравномерной прямоугольной сетке, шаблон которой изображен на рис. 2, наиболее эффективно применение балансового метода.

4.5. Разностная аппроксимация дифференциального уравнения [например, второго из (1)], записанного для условий стационарной фильтрации подземных вод, имеет вид

S_{i, j-1}(H_{i, j-1} - H_ij) + S_ij(H_{i, j+1} - H_ij) + r_i-1,j(H_{i-1, j} - H_ij) + r_ij(H_i+1,j - H_ij) + Q_{п, ij} = 0,     (10)

где

                                           (11)

Рис. 2. Шаблон неравномерной сетки

                                 (12)

              (13)

4.6. Разнесшая аппроксимация остальных уравнений из (1) и (2) осуществляется аналогично, только коэффициенты в (11) и (12) будут зависеть от искомых решений h_ij и  и к левой части (10) добавятся слагаемые Q_w,ij и , определенные как

                                                       (14)

                                        (15)

4.7. Разностная аппроксимация временной производной осуществляется в трех модификациях:

для явных и неявных вычислительных схем

                                                            (16)

для экономичных вычислительных схем (например, для схем метода переменных направлений)

                                                           (17)

для абсолютно устойчивой явной схемы Дюфорта-Франкела

                                                         (18)

Коэффициент при временной производной, умноженный на площадь блока (водоемкость блока), в разностной форме записывается в виде

                                                         (19)

4.8. При построении разностных аналогов граничных условий к разностной системе (10) необходимо добавить уравнения для узлов сеточной границы Г^* + g^*, которая пройдет, в общем случае, только вблизи границы Г + g области фильтрации D, совпадая в ней, быть может, лишь в отдельных точках. Граничные условия задаются на сетке путем «сноса».

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

4.10. Моделирование скважины, работающей с заданным напором, требует обязательного введения дополнительной фильтрационной проводимости T_ф,ij, которая вычисляется по формуле

                                               (20)

где

Dх - размер блока; r_с - радиус скважины (для неоднородной области и неравномерной сетки водопроводимость и размер блока осредняются). Тогда величина дебита скважины, расположенной в блоке (ij), выражается зависимостью

Q_c,ij = T_ф,ij(H_ij - H_c),                                                     (21)

где

Н_с и Н_ij - напор на скважине и в узле модели соответственно.

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

Н_c,ij = H_ij - Q_c,ij/Т_ф,ij,                                                     (22)

где Q_c,ij - дебит скважины.

4.12. Разностный аналог граничного условия III рода, соответствующий выражению (7) на внешнем или внутреннем контуре, имеет вид

                                       (23)

                                                        (24)

где

, H_пр,ij - абсолютные отметки уровня воды в поверхностном водоеме и подошвы подрусловых отложений, соответственно, в узле с координатами (ij); c_ij - параметр взаимосвязи подземных и поверхностных вод, учитывающий несовершенство поверхностного водоема; N_ij - длина водоема в ячейке (ij); DL_ij - эквивалентное несовершенство водоема; l_ij - минимальное расстояние между узлом ячейки и контуром водоема.

4.13. Граничные условия первого рода на моделях задаются обычно, как ГУ III, ибо в природных условиях поверхностные водоемы всегда имеют несовершенство по степени вскрытия, но если водоем расположен в центре блока и его несовершенство пренебрежительно мало, то c придается значение (10 ¸ 30)×10³ м²/сут.

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

                         (25)

4.15. Системы разностных уравнений для верхнего горизонта трехслойной зоны или для водоносного горизонта в однослойной толще - аналогичны. Следует учесть при этом, что первое уравнение из (1) и уравнение (2) - нелинейные.

4.16. Для возможности решения краевой задачи стационарной фильтрации система разностных уравнений (10) приведена к виду

r_i-1,jH_i-1,j + r_ijH_i+1,j - a_ijH_ij + s_i,j-1H_i,j-1 + s_ijH_i,j+1 = -F_ij, i = 1, 2, ..., M; j = 1, 2, ..., N,    (26)

где

                                                     (27)

                                                                            (28)

4.17. Решение системы (26) ведется методом покомпонентной верхней релаксации, который реализуется по следующей вычислительной схеме:

для однослойной зоны

                                            (29)

где  - промежуточное приближение;

для трехслойной зоны

                      (30)

                  (31)

причем вначале просчитывается верхний водоносный горизонт, а затем - нижний.

4.18. Величины ,  вычисляются по формулам

  (32)

   (33)

 и  рассчитываются по формулам аналогичным (33) и (34) только к их правым частям необходимо добавить, соответственно, выражения

 

 и  вычисляются по формулам (27) и (28).

4.19. Итерационный параметр t_опт может принимать значения 1 < t_опт < 2. Для решения практических задач t_опт рекомендуется находить подбором или вычислять различными приближенными методами. Так, В.В. Шаманским в книге «Численное решение задач фильтрации грунтовых вод на ЭЦВМ» (Киев, «Наукова думка», 1969) предложена следующая методика. Счет начинается при t = 1 и после каждой итерации вычисляется

                                                (34)

Когда первые две - три цифры в значении l перестанут изменяться, следует положить l_max = l, вычислить t_опт по формуле

                                                     (35)

и после этого счет продолжить с t = t_опт.

4.20. Для решения краевой задачи нестационарной фильтрации рекомендуется применять: явные схемы (классическую и Дюфорта-Франкела), неявную классическую и явно-неявную схему переменных направлений. Характерные особенности данных схем:

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

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

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

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

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

4.22. Размер временного шага, определяющего точность решения задачи, должен находиться в оптимальных пределах: с одной стороны решение должно удовлетворять заданной точности, а с другой - продолжительность счета на ЭВМ должна быть приемлема. Рекомендуется применять переменный временной шаг, увеличивающийся в логарифмической зависимости. Критерием размера временного шага может служить выполнение баланса на сетке, особенно в тех узлах, где заданы граничные условия II рода.

4.23. Учет нелинейности первого уравнения из (1) и уравнения (2) осуществляется в результате их линеаризации (осреднения мощности потока) не на весь период решения, а на некоторую ее часть, в течение которой положение свободной поверхности подземных вод изменяется в небольших пределах. Затем следует уточнение водопроводимости в соответствии с новым положением уровня грунтовых вод и дальнейшее решение линейного уравнения до нового нарушения пределов линейного участка. Погрешность решения не превосходит 1 - 2 % при изменении мощности водоносного горизонта в пределах одного блока на 10 - 15 %, поэтому было принято следующее условие пересчета фильтрационных проводимостей:

                                                       (36)

где e = 10 ¸ 15 %; b_1,ij - абсолютная отметка подошвы безнапорного водоносного горизонта в блоке (ij).

4.24. Для выполнения непосредственно численного решения краевой задачи, описывающей фильтрацию подземных вод в водоносных горизонтах трех - и однослойного строения в условиях стационарной и нестационарной фильтрации, был разработан пакет прикладных программ для ЭВМ ЕС-1022 под общим названием FILTR12, состоящий из 12 программных единиц типа SUBROUTINE.

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

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

4.27. Подпрограмма BLANK выводит на печать двумерные массивы размерностью М´N в форме, соответствующей области фильтрации, что повышает наглядность и скорость анализа решения.

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

4.29. Подпрограмма BALN12 выполняет расчет баланса для нестационарной задачи по следующим статьям: емкостные запасы для каждого горизонта; питание или разгрузка в поверхностные водоемы; питание или разгрузка, задаваемые в виде расходов на боковых границах, в скважинах и в виде инфильтрации; переток между горизонтами в трехслойной зоне; испарение.

Вызвать подпрограмму можно либо после каждого временного шага, либо в моменты времени, указанные в исходной информации. Отличие подпрограмм BALN12 и SNAВ12 состоит в том, что в последней выполняется расчет бокового питания - расхода для каждого узла.

4.30. Подпрограмма GRAN осуществляет временные изменения граничных условий II и III рода.

4.31. Четыре подпрограммы предназначены для решения нестационарной задачи:

JAWA12 - классическим явным методом или, если не выполняется условие устойчивости, методом Дюфорта-Франкела.

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

MPN1 и PROGON - явно-неявным методом переменных направлений с выполнением монотонной прогонки.

4.32. Подпрограмма PARM12 выполняет пересчет водопроводимости при решении нелинейной задачи.

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

4.34. При объеме памяти ЭВМ ЕС-1022 256 кб максимальное число узловых точек модели составляет примерно 3,2 тысячи. Время решения задачи зависит от применяемой вычислительной схемы. Для максимальной задачи один временной шаг по явной схеме выполняется за 5 - 7 секунд. Примерно столько же времени занимает одна итерация при решении стационарной или нестационарной задачи.

4.35. Для удобства работы с пакетом программ FILTR12, а также ввиду большого количества двумерных массивов по 2500 чисел в каждом, использовался внешний носитель памяти, в частности, личный пакет дисков ЭВМ ЕС-1022. Для организации данных на внешнем носителе (пакете дисков) необходимо произвести следующие операции: инициализацию пакета дисков, создание личной библиотеки оттранслированных программных модулей (комплекта подпрограмм FILTR12), каталогизацию модулей в эту библиотеку, организацию рабочих файлов, куда будут записываться целые и вещественные двумерные массивы. Рассмотрим подробно все эти операции.

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

где

 - регистрационный номер тома.

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

Личная библиотека - это дополнительная внешняя память. Для трансляции программ создается библиотека объектных модулей (RL). Эта библиотека имеет оглавление (содержание), где хранится информация о программе (модуле), т.е. ее имя, адрес. Наличие оглавления уменьшает затраты времени на поиск нужного модуля, так как оглавление значительно меньше самой библиотеки.

Пакет задания для создания личной библиотеки выглядит так:

// JOB RL

Оператором ASSGN назначается физическое устройство с номером 131, т.е. то устройство, на котором находится пакет дисков с именем ЛПГГУ7. В операторах DLBL и EXTENT указан адрес участка диска, где создается библиотека и ее имя (БОМ), которое будет записано в оглавлении тома. Библиотека БОМ занимает 10 цилиндров, начиная с 35-го. Срок ее хранения - бесконечен. Программа CORGZ выполняет операцию создания библиотеки. В операторе NEWVOL указывается количество цилиндров, выделенных под библиотеку - 10, и в скобках - количество дорожек, выделенных под ее оглавление - 20.

4.38. Единицей хранения в RL является объектный модуль. В библиотеку можно помещать (каталогизировать) модули, полученные в результате трансляции любым из трансляторов, и затем на этапе редактирования объединять любые из них в одну программу, готовую к выполнению. Каталогизацию в RL производит программа MAINT с помощью оператора CATALR. Существует режим трансляции DECK, в результате которого оттранслированная программа в виде объектного модуля помещается на логическое устройство SYSPCH. Этому логическому устройству назначают магнитную ленту.

Пакет задания на каталогизацию программ в RL с промежуточным хранением информации на магнитной ленте:

текст программы

С помощью этого пакета задания каталогизируются в RL все 12 подпрограмм пакета FILTR12. Для этого каждому тексту программы ставится в соответствие оператор CATALP GRAN, CATALR BLANK и т.д.

4.39. В ФОРТРАНЕ ДОС ЕС допускается последовательный тип организации файлов. Последовательная организация файлов возможна на перфокартах, магнитной ленте, дисках. При работе с FILTR12 предлагается использовать диск для организации рабочих файлов. По способу доступа к данным файлы ФОРТРАНА делятся на файлы последовательного и прямого доступа. Каждый файл, созданный на дисках, должен содержать метку, которая служит для идентификации этого файла. Информация о метках задается в операторах DLBL и EXTENT.

Обращение к файлу производится по его номеру, который указывается в операторах ввода-вывода и должен соответствовать логическому устройству ввода-вывода, указанному в операторе DLBL и EXTENT.

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

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

Чтобы организовать файл прямого доступа на диске, надо данный участок диска форматизовать, т.е. произвести разбивку участка диска на записи фиксированной длины. Делается это программой CLRDSK. Она предварительно очищает назначенный участок и затем производит его разметку. Размечаемый файл всегда имеет имя UOUT.

Пакет задания при этом имеет вид

В результате выполнения этого задания на пакете дисков с регистрационным номером ЛПГГУ7 организован файл прямого доступа с идентификатором «FILE6», срок хранения - бесконечен. Начальный адрес файла 82-й цилиндр, нулевая дорожка и занимает он 14 цилиндров. Файл разбит на записи длиной 200 байт. Эта же величина указана и в операторе описания файла DEFINE FILE 9(5000, 200, L, K4), который находится в основной программе. Вся информация о метке файла находится в оглавлении тома.

Для работы с программой FILTR12 были организованы два файла прямого доступа с идентификаторами: «FILE6» - для вещественных массивов и «FILE7» - для целых массивов.

Максимальная длина записи для вещественных массивов - 200 байт, а для целых - 100 байт. Оператор DEFINE FILE в основной программе имеет вид

DEFINE FILE 9 (5000, 200, L, K4), 11 (10000, 100, L, K).

4.40. Для записи вещественных и целых массивов на диск в «FILE6» и «FILE7» написана специальная программа DISK. Пакет задания на редактирование ее и выполнения записи чисел в файлы прямого доступа имеет вид

числа (исходная информация)

/^*

/&

4.41. Исходная информация дли программы DISK:

Простые переменные

N - признаковое число (N = 0 - записывается на диск целый массив. N = 1 - вещественный);

вводится на отдельной перфокарте в формате I5;

K1 - признаковое число (K1 = 0 - записываемый массив не выводится на печать; K1 = 1 - выводится на печать); вводится на отдельной перфокарте в формате I5;

K2 - номер записи, начиная с которой целый или вещественный массив будет записываться на диск в соответствующие файлы прямого доступа; вводится на отдельной перфокарте в формате I5.

Массивы

ITB - массив целых чисел, вводится при N = 0;

формат (13I6), 193 перфокарты;

ТВ - массив вещественных чисел, вводится при N = 1;

формат (13F6.2), 193 перфокарты.

4.42. Исходная информация для программы FILTR12.

Простые переменные

M, N - число сеточных узлов по горизонтали и вертикали прямоугольника, оконтуривающего область фильтрации;

Примечание 1. Начало отсчета - левый нижний угол; индекс I увеличивается по горизонтальной оси слева направо, индекс J - по вертикальной оси снизу вверх.

Примечание 2. С целью удобства программирования сеточная область фильтрации, дополненная до прямоугольника фиктивными блоками, оконтуривается еще одним рядом фиктивных блоков.

LPR - число элементов массива IR (управляющий массив, подробная информация о котором будет ниже);

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

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

NELS - признаковое число решения линейной (NELS = 0) или нелинейной (NELS = 1) прогнозной стационарной задачи.

NY - начальное приближение для стационарной задачи, задаваемое одним числом.

NT - конечный момент времени решения нестационарной задачи.

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

KL - число записей в файле прямого доступа, занимаемых одним массивом. Для данной программы это число всегда равно 50;

KС - номер первой записи первого (целого) массива типа INTEGER´2;

KW - номер первой записи первого вещественного массива;

LL5 - число блоков, в которых заданы граничные условия 1 рода;

Примечание 5. Перечисленные простые переменные, перфорируются на одной перфокарте (первая перфокарта) в формате (13I6).

I5, I6, J5, J6 - минимальные и максимальные координаты прямоугольника по индексу I и J соответственно, для которого будет вычисляться баланс (см. примеч. 16). Перфорируется на отдельной перфокарте (вторая перфокарта) в формате (20I4).

BI, CI - значения свободной и упругой водоотдачи.

Примечание 6. В программе BI и CI умножаются на 10^-4, поэтому в исходных данных необходимо их задавать умноженными на 10⁴.

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

EPS1 - допустимая сработка уровней (в относительных единицах) без корректировки водопроводимости.

EPSN - малое число, когда максимальная разность между двумя последовательными решениями стационарной линейной задачи становится меньше EPSN, итерационный процесс для решения стационарной прогнозной задачи завершается.

ТDV - ускоряющий итерационный параметр.

Примечание 7. Перечисленные простые переменные вещественного типа перфорируются на одной перфокарте (третья перфокарта) в формате (10F8.4).

R - множитель для изменения сопротивления подрусловых отложений при задании ГУ III.

GLKR - критическая глубина стояния грунтовых вод, ниже которой уровень грунтовых вод не влияет на величину испарения.

PS - показатель степени в формуле (4).

W0 - интенсивность испарения у поверхности земли; вводится с минусом.

Примечание 8. Перечисленные простые переменные вещественного типа перфорируются на одной перфокарте (четвертая перфокарта) в формате (10F8.4).

Массивы, необходимые для решения стационарной задачи

Ввод массивов должен производиться в память ЭВМ в той последовательности, в которой они ниже описаны.

IR - управляющий массив, содержащий признаковые числа 0 или 1, посредством которых назначается (IR(I) = 1) или отменяется (IR(I) = 0) тот или иной режим работы программы. Каждый элемент массива перфорируется на отдельной перфокарте в формате (I4) и число перфокарт должно совпадать с числом в простых переменных. Массив всегда вводится с перфокарт и всегда выводится на печать (значения элементов массива будут приведены ниже).

D1(50) - одномерный массив ячеек сетки по горизонтальной оси (в км). Формат (10F6.2), 5 перфокарт. Если простая переменная М меньше 50, то остальные элементы массива должны быть нулями. Массив вводится с перфокарт.

D2(50) - одномерный массив ячеек сетки по вертикальной оси (в км). Формат (10F6.2), 5 перфокарт. Если простая переменная меньше 50, то остальные элементы массива должны быть нулями. Массив вводится с перфокарт.

KPV(2, 80) - двумерный массив координат узлов верхнего горизонта, значения уровней в которых необходимо выводить в частичном выводе. Первая координата I, вторая - J. Формат (20I4), 8 перфокарт. Если переменная ISV имеет значение не 80, а 20, 40 или 60 (другие значения не рекомендуются), то 6,4 или 2 перфокарты остаются пустыми. Массив вводится с перфокарт.

KPN(2, 60) - то же самое для нижнего горизонта и простой переменной.

Примечание 9. В идентификаторах простых переменных и массивов последние символы V и N обозначают, как правило, верхний и нижний горизонты соответственно.

KTV (50, 50) - двумерный массив порядковых номеров зон водопроводимости для верхнего слоя верхнего горизонта. Записывается при обходе сеточной области фильтрации, начиная с нижнего левого угла направо и вверх (так записываются все большие двумерные массивы). Формат (13I6), 193 перфокарты. В фиктивных узлах надо записать 0. В память ЭВМ вводится с пакета дисков.

STV (65) - одномерный массив значений коэффициентов, на которые умножается массив ТSV в соответствии с порядковыми номерами зон. Формат (13F6.2), 5 перфокарт. Вводится с перфокарт.

KTV 1(50, 50) - двумерный массив порядковых номеров зон водопроводимости нижнего слоя верхнего горизонта. Остальное аналогично массиву KTV.

STV1 (65) - одномерный массив значений коэффициентов, на которые умножается массив TSV1 в соответствии с порядковыми номерами зон. Остальное аналогично массиву STV.

KTN (50, 50) и STN (65) - то же самое для нижнего горизонта.

Примечание 10. Массивы KТV, KТV1, KTN и STV, STV1, STN предназначены для более удобного и менее трудоемкого изменения водопроводимости при решении обратных задач. Они относятся к числу необязательных массивов и если необходимость в зонной корректировке водопроводимости отсутствует, их можно не вводить. При этом должно быть IR(9) = 0.

Подготовка этих массивов состоит в следующем. На карте литологического состава пород соответствующего водоносного слоя (горизонта) выполняется районирование по величине коэффициента фильтрации K_ф. Далее все зоны пронумеровываются. Если в разных частях области фильтрации имеются зоны с одинаковыми K_ф, то номера зон должны быть разными. Затем на карту накладывается сетка, контуры зон аппроксимируются прямолинейными границами, и каждому узлу в зоне присваивается один и тот же номер - номер зоны. Эти номера зон и являются элементами массивов KТV, KТV1 и KTN.

В массивах STV, STV1 и STN содержится по одному коэффициенту для каждой зоны и порядковый номер элемента этих массивов соответствует номеру зоны. Перед вычислением фильтрационного сопротивления между двумя узлами первоначальная водопроводимость TSV, ТSV1 или TSN умножается на соответствующий коэффициент из STV, STV1 или SТN. Так, например, в третьей зоне верхнего слоя верхнего горизонта водопроводимость надо уменьшить в 4 раза, а в пятой зоне увеличить в 2 раза. Остальные массивы водопроводимости должны быть исходными. Тогда массив STV должен быть: 1, 1, 0,25, 1, 2 ..., остальные единицы, а массивы STV1 и STN состоят из одних единиц.

TSV (50, 50) - водопроводимость в ячейках сетки для верхнего слоя верхнего горизонта (м²/сут). Формат (13I6), 193 перфокарты, в фиктивных узлах записываются нули. В ЭВМ вводится с пакета дисков.

TSV1 (50, 50) - то же самое для нижнего слоя верхнего горизонта.

TSN (50, 50) - то же самое для нижнего горизонта.

Примечание 11. Для однослойной толщи водопроводимость делится пополам и значение каждой половины записывается в массивы TSV и TSN.

При решении обратной задачи по уточнению фильтрационных свойств пород может оказаться, что необходимо изменить границы зон указанных в массиве KTV, KTV1 или KTN. Для того чтобы не перебивать данные массивы составлена специальная подпрограмма REDTS, с помощью которой можно корректировать значения водопроводимости внутри выделенных зон, а также значения водопроводимости в отдельных блоках. С этой целью после массива STN. или, если IR/9/ = 0, то после массива KPN вводится редактирующий массив REDT. Перфорируется он следующим образом: на отдельной перфокарте в формате 4I3 набиваются минимальные и максимальные координаты прямоугольника по индексу I и J соответственно, для которого будет корректироваться водопроводимость, затем в 14-й позиции поля перфокарты - 1 (если вводится корректирующий коэффициент) или 2 (если будет вводиться измененное значение водопроводимости) и далее в формате F8.4 пробивается либо значение коэффициента, либо само новое значение водопроводимости. Число таких перфокарт в принципе неограничено. Пустая перфокарта служит признаком завершения редакции массива ТSV. Аналогично производится корректировка массивов TSV1, TSN. Если необходимость в дополнительной корректировке отсутствует, то при этом должны быть IR/72/ = IR(73) = 0. В противном случае ore; равны единице.

Z (50, 50) - массив абсолютных отметок поверхности земли. В местах отсутствия испарения записывается 100. В фиктивных узлах записывается ноль. Формат (13F6.2), 193 перфокарты. Вводится при IR(79) = 1.

KWV (50, 50) - двумерный массив порядковых номеров зон инфильтрации и ГУ II для верхнего горизонта. Аналогичен массиву KТV.

SWV (195) - одномерный массив значений коэффициентов, на которые умножается массив W2V в соответствии с порядковыми номерами зон. Перфорируется в формате (13F6.2), 15 перфокарт. Вводится с перфокарт.

KWN (50, 50) - аналогичен массиву KWV, только для нижнего горизонта. Для однослойного потока в соответствующее для него место в файле прямого доступа на пакете дисков записываются единицы в формате (13I6), 193 перфокарты.

SWN (39) - аналогичен массиву SWV, только для нижнего горизонта. Вводится с перфокарт после массива SWV в формате (13F6.2), 3 перфокарты. Для однослойного потока его значения всегда равны единице.

Перечисленные массивы KWV, SWV, KWN, SWN относятся к числу необязательных массивов и вводятся при задании IR/10/ = 1.

W2V (50, 50) - инфильтрация и ГУ II в ячейках сетки верхнего горизонта. Формат (13F6.2), 193 перфокарты, в фиктивных узлах записываются нули. Ввод в память ЭВМ осуществляется с пакета дисков.

Примечание 12. Для определения инфильтрации (ГУ II) следует поступить аналогично определению водопроводимости (см. примеч. 10). Но необходимо это помнить и при задании дебитов скважин и их изменений в процессе решения задачи. Так, если в третьей зоне инфильтрация уменьшена, например, в четыре раза относительно первоначальной, задаваемой массивом W2V, то дебит скважины в этой зоне следует задавать в 4 раза больше необходимого. Ибо этот дебит, будучи умноженным на 0,25 дает нужный.

W2N (50, 50) - то же самое для нижнего горизонта.

Примечание 13. В однослойной толще и инфильтрация и дебит скважин задается только в массив W2V, а массив W2N остается пустым.

Решение обратных задач, связанных с уточнением величин утечек из во донесу щих коммуникаций, с помощью массивов KWV и SWV удобно лишь на первом этапе. Поэтому, в дальнейшем корректировка массива W2V должна осуществляться с помощью подпрограммы REDTS. Делается это с помощью массива REDW, который вводится сразу же после массива STN (если IR/10/ = 0), либо после массива REDT (если IR/73/ = 1). Перфорируется этот массив аналогично массиву REDT, после последней редактирующей перфокарты ставится пустая карта. Если необходимость в дополнительной корректировке отсутствует, то следует положить IR/74/ = 0, в противном случае - единице.

H2 (50, 50) - абсолютные отметки кровли нижнего горизонта (в однослойной толще в массив следует записать нули). Формат (13I6), 193 перфокарты, в фиктивных узлах записываются нули. В ЭВМ вводится с пакета дисков.

TR (50, 50) - параметр раздельного слоя: K₀/m₀×10⁶. В однослойной толще следует записать нуль. Формат (13I6), 193 перфокарты, в фиктивных узлах записываются нули. Вводится с пакета дисков.

KR (50, 50) - двумерный массив порядковых номеров зон коэффициентов фильтрации раздельного слоя. Аналогичен KTV.

SR (65) - одномерный массив значений коэффициентов, на которые умножается массив TR в соответствии с порядковыми номерами зон: Формат (13F6.2), 5 перфокарт, вводится с перфокарт.

Для однослойной толщи массивы KR и SR не вводятся, следовательно, IR(55) всегда равно нулю.

HR (50, 50) - массив абсолютных отметок воды в поверхностном водоеме. Формат (13F6.2), 193 перфокарты; в узлах, где реки отсутствуют и в фиктивных узлах записывается ноль. Вводится с пакета дисков.

Примечание 14. Если отметка воды в водоеме действительно нуль, то в исходных данных следует записать 0.01.

DL (50, 50) - массив параметров взаимосвязи подземных и поверхностных вод (м/сут). В фиктивных узлах и в узлах, где рек нет, в массив следует записать ноль. Формат (13I6), 193 перфокарты, вводится с пакета дисков.

Примечание 15. Максимально допустимое значение параметра взаимосвязи не должно превосходить 32 тыс.м²/сут.

H4 (50, 50) - массив абсолютных отметок подошвы подрусловых слабопроницаемых отложений. Когда уровень подземных вод становится ниже отметок подошвы подрусловых отложений, происходит отрыв уровней подземных и поверхностных вод и ГУ III переходит в ГУ II. В фиктивных узлах, где рек нет, следует записать нуль. Формат (13F6.2), 193 перфокарты, вводится с пакета дисков.

KRB (2, LL5) - двумерный массив координат узлов, где заданы ГУ I в верхнем водоносном горизонте. Первая координата I, вторая - J . Формат (20I4), число перфокарт зависит от значения LL5, которое должно быть не более 200. Массив вводится с перфокарт, при этом должно быть IR(8) = 1.

HRB (LL5) - одномерный массив значений ГУ I (уровни воды в поверхностном водоеме). Формат (13F6.2), число перфокарт зависит от значения LL5, которое должно быть не более 200. Последовательность заполнения массива HRB должна совпадать с последовательностью заполнения массива KRB. Массив вводится с перфокарт.

KBAL (2, 80) - двумерный массив координат узлов, в которых должен вычисляться баланс по подпрограммам BALN12 и SHAB12 при IR(66) = 1. Формат (20I4), 8 перфокарт, вводится с перфокарт. Первая координата I, вторая - J.

Примечание 16. При вычислении баланса предусмотрено два режима. В первом режиме баланс считается для всех узлов прямоугольника, координаты углов которого заданы простыми переменными I5, I6, J5, J6. Во втором режиме баланс считается только в узлах, координаты которого указаны в массиве KBAL.

UNV (50, 50) - массив начальных условий для решения стационарной или нестационарной задачи. Задание его не обязательно: для решения стационарной задачи начальное приближение можно задать одним числом (простая переменная NY), а для решения нестационарной задачи берется результат решения стационарной задачи. Однако этот массив может использоваться в других целях. Если записать в него значения уровней, которые должны получаться в результате решения стационарной обратной задачи, то на печать можно выводить решение стационарной задачи в разностях: между тем, что должно быть и тем, что получилось. При этом надо задать IR(4) = 0; IR(50) = IR(52) = 1. Такой вывод значительно упрощает анализ результатов решения обратной стационарной задачи.

В фиктивных узлах в массив следует записать ноль. Формат (13F6.2), 193 перфокарты, вводится с пакета дисков.

UNN (50, 50) - то же самое для нижнего горизонта. Для однослойной толщи этот массив совпадает с массивом UNV. Вводится с пакета дисков.

При необходимости решения только стационарной задачи (например, при решении обратной стационарной задачи требуется многовариантный счет) необходимо задать IR(58) = 1.

Для решения нестационарной задачи следует добавить: Н1(50, 50) - массив абсолютных отметок подошвы нижнего горизонта Формат (13F6.2), 193 перфокарты, в фиктивных узлах следует записать ноль. Вводится с пакета дисков. Для однослойного горизонта этот массив выступает в роли абсолютных отметок его подошвы.

H3 (50, 50) - массив абсолютных отметок подошвы верхнего слоя верхнего горизонта. Формат (13I6), 193 перфокарты, в фиктивных узлах в массив следует записать ноль. Вводится с пакета дисков. Для однослойного горизонта в этот массив следует записать нули.

Информация, необходимая для изменения ГУ II во времени и в процессе решения нестационарной задачи

Для верхнего горизонта:

Т1 (30) - моменты времени изменения ГУ II, сут. Формат (10F8.2), 3 перфокарты. Последним (отличным от нуля) элементом массива T1 должен быть конечный момент решения нестационарной задачи;

W1N (2847) - новое значение ГУ II, м³/сут. Формат (13F6.2), 219 перфокарт;

KD1 (2, 2840) - координаты узлов, где должны быть изменены ГУ II. Формат (20I4), 284 перфокарты;

NW1 (40) - число одновременно изменяемых граничных условий (в возрастающей последовательности). Формат (20I4), 2 перфокарты. Все перечисленные выше массивы вводятся с перфокарт.

Примечание 17. Пример. В начальный момент времени (t = 0) в узлы с координатами 8/14 и 12/15 необходимо задать утечки соответственно равными 11,8 и 3 м³/сут. Через 0,5 года величина утечек в первом блоке должна увеличиться вдвое, а еще через 0,5 года утечки надо ликвидировать. Прогноз выполнить на 10 лет. Тогда массивы будут иметь вид:

T1 = /0, 182.5, 365, 3650, 0, 0, .../

W1N = /11.8, 3, 23.6, 11.8, 0, 0, …/

KD1 = /8, 14, 12, 15, 8, 14, 8, 14, 12, 15, 0, 0, .../

NW1 = /2, 3, 5, 0, 0, 0, .../

Для нижнего горизонта:

Т2 (30) - вводится аналогично массиву T1;

W2P (65) - вводится аналогично массиву W1N. Формат (13F6.2), пять перфокарт;

KD2 (2, 70) - вводится аналогично массиву KD1. Формат (20I4), 7 перфокарт;

NW2 (40) - вводится аналогично массиву NW1. Формат (20I4), две перфокарты.

Информация, необходимая для изменения ГУ III во времени и в процессе решения нестационарной задачи

T3 (40) - аналогичен массиву T1, четыре перфокарты;

HRN (2262) - новые значения абсолютных отметок уровней поверхностных вод. Формат (13F6.2), 174 перфокарты;

KD3 (2, 2260) - аналогичен массиву KD1. Формат (20I4), 226 перфокарт;

NHR (40) - аналогичен массиву NW1. Формат (20I4), 2 перфокарты.

Все перечисленные массивы вводятся с перфокарт.

MIV (50, 50) - массив коэффициентов свободной водоотдачи для верхнего слоя верхнего горизонта. Используется, если водоотдача не задана одним числом (BI). Величину водоотдачи следует записывать в массив умноженной на 10⁴. В фиктивных узлах записывается ноль. Формат (13I6), 193 перфокарты, вводится с пакета дисков.

MIN (50, 50) - то же самое для нижнего горизонта, если считать его безнапорным (до тех пор, пока горизонт напорный, упругая водоотдача задается одним числом CI).

Примечание 18. При решении обратной нестационарной задачи часто возникает необходимость в корректировке свободной водоотдачи.

REDM - корректирующие массив величин свободной водоотдачи. Аналогичен массиву REDT. Признаком того, что массив MIV отредактирован служит пустая перфокарта. Вводится при IR(81) = 1.

D(130) - одномерный массив размеров временных шагов, сут. Формат (10F8.2), 13 перфокарт, вводится с перфокарт.

ТВ (130) - одномерный массив моментов времени изменения временного шага, сут. Первый элемент - нуль, последний (отличный от нуля) - конечный момент времени решения нестационарной задачи. Формат (10F8.2), 13 перфокарт, вводится с перфокарт.

ТВV (70) - массив моментов времени вывода на печать частичного решения для верхнего горизонта, сут. Формат (10F8.2), 7 перфокарт, вводится с перфокарт.

TVN (70) - то же самое для нижнего горизонта. Для однослойного горизонта вводятся семь пустых перфокарт.

TBL (70) - массив моментов времени расчета баланса, сут. Формат (10F8.2), 7 перфокарт, вводится с перфокарт. На печать выводится суммарный баланс.

TAB (40) - массив моментов времени вывода на печать таблицы локального (поточечного) баланса. Аналогичен ТВV, 4 перфокарты.

РР(40) - массив моментов времени вывода на печать массивов решения НV и HN. Формат (10F8.2), 4 перфокарты, вводится с перфокарт.

На этом ввод исходной информации заканчивается.

4.43. Выбор режима работы программ осуществляется числовым массивом IR. Ниже приводится перечень всех режимов работы.

FILTR12. Число слева является порядковым номером элемента массива IR; режим описанный справа, выполняется при равенстве единицы данного элемента.

1 - вывод на печать массива ТSV;

2 - вывод на печать массива TSN;

3 - вывод на печать массивов D1 и D2, а также вызов подпрограмм BALN12 вместе с подпрограммой SNAB12 для расчета баланса;

4 - вывод на печать массивов KРV и KPN, а также запись результатов решения стационарной задачи на диск в массивы UNV и UNN;

5 - вывод на печать массива KTV;

6 - вывод на печать массива KTN;

7 - вывод на печать разности массивов (HV - HN) для установления направленности потоков между горизонтами;

8 - ввод массивов KRB и HRB для задания ГУ I;

9 - ввод массивов KTV, KТV1, KTN, STV, STV1, STN;

10 - ввод массивов KWV, KWN, SWV, SWN;

11 - вывод на печать массива KWV;

12 - вывод на печать массива KWH;

13 - вывод на печать массива W2V;

14 - вывод на печать массива W2N;

15 - всегда единица;

16 - всегда ноль;

17 - вывод на печать массива H1;

18 - вывод на печать массива Н2;

19 - вывод на печать массива H3;

20 - вывод на печать массива TR;

21 - вывод на печать массива HR;

22 - вывод на печать массива DL;

23 - вывод на печать массива H4;

24 - ввод массивов UNV, UNN;

25 - вывод на печать массивов UNV, UNN

26 - вывод на печать массива МIV;

27 - вывод на печать массива MIN;

28 - вывод на печать массива TFV - массив коэффициентов фильтрации верхнего слоя верхнего горизонта;

29 - вывод на печать TFN - массив коэффициентов фильтрации нижнего горизонта;

30 - вывод на печать массива НV (уровни верхнего горизонта) после решения стационарной задачи;

31 - вывод на печать массива HN после решения стационарной задачи. Для однослойного потока всегда ноль;

32 - ввод массивов T1, W1N, KD1, NW1 и реализации ГУ II для верхнего горизонта;

33 - вывод на печать массивов T1, W1N, KD1, NW1;

34 - ввод массивов Т2, W2Р, KD2, NW2 и реализация ГУ II для нижнего горизонта. Для однослойного потока всегда ноль;

35 - вывод на печать массивов T2, W2P, KD2, NW2;

36 - ввод массивов T3, HRN, KD3, NHR и реализация ГУ III;

37 - вывод массивов T3, HRN, KD3, NHR;

38 - вызов подпрограммы REL12 из JAWА12;

39 - вывод на печать массивов D, ТВ, ТВV, TBN, TBL, TAB, РР;

40 - вызов подпрограммы JAWA12;

41 - вызов подпрограммы STAL12;

42 - вызов подпрограммы REL12 из основной программы;

43 - вызов подпрограммы PARM12;

44 - вывод части решения для верхнего горизонта;

45 - вызов подпрограммы BALN12;

46 -

47 - всегда ноль;

48 - вывод таблицы баланса;

49 - ввод массивов MIV и MIN;

50 - подготовка вывода решения стационарной задачи в разностях;

51 - вывод решения в разностях (НУ - UNV) для верхнего горизонта;

52 - вывод решения в разностях (HN - UNN) для нижнего горизонта. Для однослойного горизонта всегда ноль;

53 - всегда ноль;

54 - всегда единица;

55 - ввод массива KR;

56 - вывод на печать массива KR;

57 - ввод массивов H1 и H3;

58 - решение только стационарной задачи;

59 - вывод на печать массива решения НV в соответствии с данными массива РР;

60 - вывод на печать массива решения HN в соответствии с данными массива PP. Для однослойного горизонта всегда ноль;

61 - всегда ноль;

62 - вызов подпрограммы BALN12 после решения стационарной задачи;

63 - вызов подпрограммы SHAВ12 после решения стационарной задачи;

64 - вывод части решения для нижнего горизонта. Для однослойного потока всегда ноль;

65 -

66 -

67 - вызов подпрограмм MPN1 и PROGON;

68 - всегда ноль;

69 - ввод массива ТSV1;

70 - вывод массива TSV1;

71 - вывод массива KТV1;

72 - редактирование массива TSV;

73 - редактирование массива TSV1;

74 - редактирование массива W2V;

75 - запись на диск отредактированного массива TSV;

76 - запись на диск отредактированного массива TSV1;

77 - запись на диск отредактированного массива W2V;

78 - вывод на печать всего решения нестационарной задачи для верхнего слоя верхнего водоносного горизонта в разностях между HV на данный момент времени (в соответствии с массивом РР) и НV на начальный момент времени;

79 - ввод массива Z;

80 - вывод массива Z;

81 - редактирование массива MIV;

82 - запись на диск отредактированного массива МIV.

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

При выдаче части решения или всего массива решения печатается надпись: «Время с начала отсчета T =» с текущим временем и именем массива. При выводе части решения сначала выводится строка координат узлов, а под каждой парой координат печатается значение уровня или величина понижения.

Вывод результатов расчета баланса по BALN12 и SHAB12 предшествует выводу уровней. Если вводится таблица баланса, то в соответствующих колонках проставлены координаты узлов, для которых рассчитывается баланс. Суммарный баланс включает в себя составляющие либо тех узлов, координаты которых выведены в таблице, либо узлов, попавших в заданный прямоугольник.

Если последним управляющим оператором является «конец задания», то счет заканчивается сообщением о продолжительности работы программы (в часах, минутах, секундах) с момента ввода первого управляющего оператора.

4.45. Пакет задания для редактирования, вызываемых из «БОМ» подпрограмм, и выполнения счета имеет следующий вид:

исходная информация, вводимая с перфокарт.

5. ВЫВОДЫ

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

Приложение

ПРИМЕР
РЕШЕНИЯ НА ЭВМ ГЕОФИЛЬТРАЦИОННОЙ ЗАДАЧИ
ПРОГНОЗА УРОВЕННОГО РЕЖИМА ГРУНТОВЫХ ВОД,
ФОРМИРУЮЩЕГОСЯ В УСЛОВИЯХ ТРЕХСЛОЙНОГО
БЕЗНАПОРНО-НАПОРНОГО ПОТОКА
ПРИ ПОДТОПЛЕНИИ И ПОДПОРЕ ОТ ВОДОХРАНИЛИЩА

1. Основные методологические положения решения на ЭВМ рассматриваемой геофильтрационной задачи

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

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

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

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

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

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

2. Характеристика природной обстановки исследуемой области фильтрации

Исследуемая область фильтрации представляет зону промышленной и селитебной застройки площадью ~ 40 км², ограниченную на юге рекой, на севере - озерами и болотами. Восточная и западная границы проведены по нормали к реке на расстоянии 2,2 - 2,5 км от боковых контуров промышленной зоны. Гидрографическая сеть представлена, кроме того, мелкими реками, озерами и болотами, а также осложнена системой каналов (рис. 3).

В геоморфологическом отношении территория находится в пределах поймы, первой и второй надпойменных террас.

В геологическом строении принимают участие отложения пермской и четвертичной систем (рис. 4). Отложения пермской системы представлены сакмарским ярусом нижнего отдела системы, нижнеказанским подъярусом и татарским ярусом верхнего отдела системы. Образование четвертичной системы сложены несколькими генетическими типами пород. Наиболее широко развиты современные и верхнечетвертичные аллювиальные отложения, представленные песками, супесями и суглинками. Более ограниченное распространение имеют озерно-болотистые и техногенные отложения.

Рис. 3. Схематическая карта техногенных условий моделируемой области фильтрации

1 - территория промзоны и микрорайонов (I, II, III); 2 - водонесущие коммуникации; 3 - участки с наличием инфильтрации атмосферных осадков и испарения с уровня грунтовых вод

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

Современные аллювиальные отложения слагают пойму реки.

Верхнечетвертичные аллювиальные отложения в пределах исследуемой территории развиты на первой и второй надпойменной террасах.

Верхнечетвертичный и современный аллювиальные водоносные горизонты, плановая граница между которыми проходит по тыловому шву поймы, объединены в один двухслойный аллювиальный водоносный горизонт. В верхнем слое коэффициент фильтрации мелкозернистых песков изменяется от 2,4 до 18 м/сут, коэффициент водопроводимости от 26 до 310 м/сут, максимальное значение коэффициента водоотдачи от 0,02 до 0,04, коэффициент фильтрации пылеватых песков - 1 м/сут, коэффициент водоотдачи - 0,04.

В нижнем слое значения коэффициентов фильтрации относительно однородны, их средняя величина составляет 25 м/сут, а коэффициент водоотдачи 0,2.

Отложения татарского яруса представлены глинами трещиноватыми с линзами и маломощными прослоями мергеля и известняка, алевролитами, мергелями, гипсами. Мощность татарских отложений увеличивается от 0 до 20,3 м в восточном и юго-восточном направлениях от контура их размыва в северо-западной части участка, кровля татарских отложений сильно эродирована, ее абсолютные отметки колеблются от 43 до 55 м.

Рис. 4. Фильтрационная схематизация гидрогеологического разреза исследуемой территории

Р₂k₁ - индекс водовмещающих отложений; k - коэффициент фильтрации, м/сут; m - свободная водоотдача; m^* - упругая водоотдача

Коэффициент фильтрации глин принят условно равным 0,1 м/сут. Практически этот водоносный комплекс является относительным водоупором, внутри которого имеются литологические «окна» с коэффициентом фильтрации 8 м/сут.

Отложения нижнеказанского подъяруса имеют повсеместное распространение. Они представлены известняками, неравномерно трещиноватыми, иногда кавернозными. Мощность известняков 17,2 - 30,1 м.

Известняки на западе и северо-западе участка перекрываются верхнечетвертичными аллювиальными отложениями, на остальной территории - несогласно, с размывом, отложениями татарского яруса.

Значения коэффициента водопроводимости известняков в зоне выклинивания татарских отложений (мощность татарских глин 0,8 - 2,5 м) составляют 2440 (коэффициент фильтрации 142 м/сут) и 2070 м²/сут (коэффициент фильтрации 112 м/сут), величина водоотдачи соответственно равна 0,06 и 0,001. В зоне, где известняки перекрыты татарскими отложениями мощностью от 3 - 5 до 20 м, значения коэффициента водопроводимости изменяются от 225 до 660 м²/сут (коэффициент фильтрации изменяется от 9,1 до 33,2 м/сут, среднее значение 18 м/сут).

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

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

В пределах распространения нижнеказанского водоносного горизонта выделено две зоны с различными фильтрационными параметрами: 1 зона - K = 120 м/сут, m^* = 0,03; 2 зона - K = 20 м/сут, m^* = 0,001.

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

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

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

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

3. Характеристика техногенной обстановки исследуемой территории

Исследуемая территория (см. рис. 3) характеризуется интенсивной промышленной и селитебной застройкой. Под зданиями и сооружениями находится 2800 га (61 % всей площади).

Селитебная территория площадью 1300 га схематизируется пятью типами застройки: индивидуальной, одноэтажной с приусадебными участками; 2 - 3-этажной старой; 5-этажной новой; многоэтажной (9, 10 и 14-этажной) новой.

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

В пределах селитебной застройки также выделены спортивные площадки, зеленые насаждения общего пользования, заасфальтированные площадки, бани, котельные.

На исследуемой территории расположена 21 промышленная площадка, основная из которых размером в 750 га. На ней расположено (в соответствии с генеральным планом) 189 цехов и различных сооружений. Размеры отдельных цехов: 400´210 м, 640´300 м. Основные цехи завода расположены в несколько рядов. Расстояния между цехами 50 - 200 м. Между цехами и сооружениями уложена густая сеть водонесущих коммуникаций: противопожарный и производственный водопроводы, фекальная канализация, производственно-ливневая канализация. теплотрассы (наземного и подземного заложения), по которым имеются: нормативный, фактический и плановый (ориентировочный) расход воды.

Общая сеть только водопроводных коммуникаций достигает 335 км (см. рис. 4). Суммарное водопотребление оценивается 340 - 350 м/сут.

Эти данные позволяют рассчитать величину нормативных утечек и ориентировочно оценить максимальные возможные утечки (суммарные потери воды ориентировочно могут составлять 10 - 20 % от общего водопотребления завода).

4. Характеристика источников, формирующих уровенный режим грунтовых вод

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

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

режиме атмосферных осадков, коэффициенте их инфильтрации и наличии площадей, на которых возможна инфильтрация;

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

фактическом расходе воды по водонесущим коммуникациям.

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

Величины гидравлического сопротивления подрусловых отложений водотоков и водоемов DL, рассчитанные по данным режимных наблюдений за уровнем грунтовых и поверхностных вод по общепринятой методике, изменяются от 435 до 2650 м.

В схематизации фильтрационных водопроводимостей при решении эпигнозных и прогнозных задач в нестационарной постановке используются уточненные величины DL.

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

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

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

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

На участках каналов, внутренних рек и водоемов, расположенных за пределами зоны влияния подпора, уровенный режим, учитывая незначительную величину его годовой амплитуды (0,2 - 0,5 м), принимался в течение эпигнозного и прогнозных периодов постоянным в соответствии с результатами схематизации для условий стационарного режима фильтрации грунтовых вод.

Схематизация абсолютных отметок дна каналов и водоемов выполнена с учетом их глубины и результатов схематизации в них уровня воды.

Схематизация инфильтрационного питания грунтовых вод выполнена для условий их нестационарного режима фильтрации в соответствии с:

режимом атмосферных осадков при решении инверсной и обратной задач в эпигнозной постановке, а также при прогнозе уровенных режимов грунтовых вод;

коэффициентом инфильтрации атмосферных осадков;

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

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

Коэффициент инфильтрации атмосферных осадков вычислен из отношения объема воды, поступившей на уровень грунтовых вод в пределах выделенной площадки, к количеству выпавших осадков в зимний и летний периоды. Инфильтрация задавалась на площадях с отсутствием застройки и асфальтовых покрытий.

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

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

Для условий нестационарного режима фильтрации грунтовых вод при решении инверсных и обратных задач величина испарения с уровня грунтовых вод на незастроенных и незаасфальтированных площадях схематизировалась в течение эпигнозного периода в виде двух значений интенсивности испарения с уровня грунтовых вод у поверхности земли w₀ = -0,001 и -0,002 м/сут при п = 1 и 2 и критической глубине грунтовых вод 2,5 м.

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

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

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

5. Расчетная схема моделируемой области

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

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

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

Естественными внешними плановыми границами являются река, протекающая на юге исследуемой территории, и болота с озерами на севере. Для условий стационарного режима фильтрации грунтовых вод на них задавались попеременно ГУ I и II родов, а для условий нестационарного режима только ГУ III.

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

На водонесущих коммуникациях и котловане для обоих видов режима фильтрации грунтовых вод задавались ТУ II рода, постоянные во времени.

На каналах, внутренних реках и озерах для условий стационарного режима фильтрации грунтовых вод задавались попеременно либо ГУ I рода, либо ГУ III рода, а для условий нестационарного режима только ГУ III рода.

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

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

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

Вторые являются исходными при решении прогнозных задач с различными обеспеченностями режимов уровня реки и атмосферных осадков. Они соответствуют прогнозируемым установившимся поверхностям грунтовых и подземных вод при НПГ.

6. Система уравнений, описывающих моделируемый фильтрационный процесс

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

T_v = k₁(h - b₁) + T₂,                                                    (37)

где Т_v - суммарная водопроводимость верхнего водоносного горизонта; k₁ - коэффициент фильтрации верхнего слоя аллювиальных отложений; h и b₁ - абсолютные отметки уровня грунтовых вод и подошвы верхнего слоя аллювиальных отложений; T₂ - водопроводимость нижнего слоя аллювиальных отложений.

Водопроводимость Т_VN для зоны D₂ будет включать сумму водопроводимостей T_V и T_N

T_VN = Т_V + Т_N,                                                          (38)

где T_N - водопроводимость нижнего водоносного горизонта, приуроченного к отложениям казанского яруса.

Краевые условия (граничные и начальные) идентичны соотношениям (5) - (8).

7. Построение сеточной модели исследуемой области фильтрации

Построение модели исследуемой области фильтрации сводилось к ее дискретизации путем наложения на нее прямоугольного шаблона. При этом соблюдались следующие требования:

Рис. 5. Модель области фильтрации

1 - модельная аппроксимация внешних плановых границ; 2 - аппроксимация внешних плановых границ в виде непроницаемого контура; 3 - аппроксимация внешних плановых границ в виде ГУ III рода; 4 - аппроксимация внутренних граничных контуров (каналов) в виде ГУ III рода; 5 - аппроксимация внутренних граничных контуров (озер) в виде ГУ III рода; 6 - аппроксимация внутренних граничных контуров (водонесущих коммуникаций) в виде ГУ II рода; 7 - блоки, в которых задавалась инфильтрация атмосферных осадков и испарение с УГВ; 8 - контур распространения отложений татарского яруса

конечно-разностная аппроксимация исследуемой области фильтрации должна обеспечить наибольшую адекватность модели природным и техногенным условиям;

участки с идентичными, геоморфолого-литологическими признаками и районы, выделенные при схематизации городской застройки, следует вписать в соответствующую систему блоков;

водонесущие коммуникации необходимо разместить в блоках с наименьшими размерами;

как можно точнее выполнить ступенчатую аппроксимацию внешних и внутренних границ с отнесением их к центрам блоков сетки;

разумно минимизировать количество блоков сетки, чтобы не допустить больших затрат машинного времени при моделировании на ЭВМ.

В результате выполненной дискретизации составлена сеточная модель исследуемой области фильтрации (рис. 5), состоящая из 676 (26´26) блоков для одного водоносного горизонта. Из них действительных 536. К зоне D₁ приурочено 506 блоков, к зоне D₂ - 30. Число блоков с ГУ II рода - 318, из них: 34 моделируют внешние плановые границы в виде непроницаемого контура и 284 моделируют утечки из водонесущих коммуникаций. Число блоков с ГУ III рода - 108, из них: 48 моделируют внешние плановые границы, включая основную реку и озера; 60 моделируют внутренние водотоки в виде каналов, озер, мелких рек. Инфильтрация атмосферных осадков задавалась на территории не подверженной экранированию. К ней приурочено 202 блока. На этой же территории учитывалось испарение с уровня грунтовых вод.

8. Перенесение модели области фильтрации на ЭВМ

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

Как отмечалось в п. 4.39, наборы данных в виде двумерных массивов (50´50) организовывались на внешнем носителе (пакете дисков) с помощью рабочих файлов (для целых массивов ¢FILE¢ Т, вещественных - ¢FILE6¢) прямого доступа. Для этого на первом этапе осуществлялась перфорация двумерных массивов, на втором - запись этих массивов в соответствующие файлы на пакете дисков с помощью программы DISK.

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

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

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

9. Установление функциональной надежности расчетной схемы моделируемой области

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

Идентификационные исследования ведутся в двух постановках: в условиях воспроизведения стационарного и нестационарного режима фильтрации грунтовых вод.

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

выявлены величины утечек из водонесущих коммуникаций;

оценены коэффициенты перетока;

обоснован род граничных условий, задаваемый на внешних и внутренних линейных границах;

определена фильтрационная проводимость подрусловых отложений рек, озер, болот и каналов;

уточнены поля распределений водопроводимостей верхнего и нижнего слоев водоносного горизонта в аллювиальных отложениях;

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

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

уточнены значения свободной водоотдачи верхнего слоя аллювиальных отложений;

оценены величины инфильтрации атмосферных осадков. Полученные при этом для каждого из шести периодов коэффициенты инфильтрации были использованы при решении прогнозных задач;

определены параметры задания испарения с уровня грунтовых вод;

обоснована возможность применения при решении прогнозных задач явной вычислительной схемы с уточнением по неявной с переменным временным шагом и соответствующими значениями e.

Таким образом, результаты идентификационных исследований позволили перейти к прогнозу уровенного режима грунтовых вод, обусловленного самоподтоплением исследуемой территории и подпором от водохранилища.

10. Прогноз уровенного режима грунтовых вод в пределах исследуемой территории в связи с самоподтоплением и подпором от водохранилища

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

На I этапе прогноза уровенного режима грунтовых вод осуществлена схематизация режимообразующих факторов, в результате которой:

режим основной реки представлен абсолютными отметками уровня воды при НПГ, зимних попусках и для условий естественных колебаний соответствующими 1 %-ной обеспеченности;

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

режим атмосферных осадков принят 50 %-ной обеспеченности;

испарение принято с параметрами w₀ = -0,0015 м/сут при n = 1,5;

боковые границы представлены в виде непроницаемого контура.

На II этапе прогноза уровенных режимов грунтовых вод осуществляется последовательное воспроизведение на ЭВМ процесса подтопления исследуемой территории путем решения серии нестационарных задач при параметрах испарения w₀ = -0,0015 м/сут и n = 1,5 и непроницаемых боковых границах в условиях:

наличия НПГ в основной реке и исключения атмосферных осадков;

1 %-ной обеспеченности уровенного режима основной реки и исключении атмосферных осадков;

1 %-ной обеспеченности уровенного режима основной реки и 50 %-ной обеспеченности режима атмосферных осадков.

Подготовка дополнительной информации, необходимой для задания ГУ II и III рода при решении прогнозных нестационарных задач с целью воспроизведения вышеприведенных обеспеченностей уровенных режимов рек и атмосферных осадков, выполнена в соответствии с требованиями методики, приведенной в п. 4.42.

Воспроизведение процесса формирования уровенного режима грунтовых вод при НПГ от начальных условий до наступления стационарного положения сопровождалось выдачами на АЦПУ прогнозируемых в каждом блоке модели абсолютных отметок и величин подъема уровня грунтовых вод через каждые 90 сут. Формирование процесса подтопления практически стабилизировалось на 540-е сутки. Однако для получения гарантированной стационарной поверхности, являющейся начальной для воспроизведения режимов грунтовых вод, которые формируются при комплексном влиянии уровня рек, атмосферных осадков и испарения, период прогноза был продлен еще на полгода с выдачей результатов на 630 и 720-е сутки. Сопоставление двух последних выдач подтвердило стабилизацию процесса, а поверхность уровней грунтовых вод, сформировавшаяся на 720-е сутки, принята в качестве начальных условий для остальных прогнозных вариантов. Их воспроизведение охватывало годовой цикл формирования соответствующего уровенного режима, динамика которого представлена семью выдачами на АЦПУ прогнозируемых в каждом блоке модели абсолютных отметок и величин подъема уровня грунтовых вод на 140, 190, 200, 210, 240, 270 и 360-е сутки от начала зимних попусков в реке. Первая выдача (на 140 сутки) характеризует уровенную поверхность грунтовых вод, формирующуюся под влиянием зимних попусков; вторая (на 190-е сутки) - под влиянием восходящей ветви графика уровенного режима половодья реки; третья, четвертая, пятая и шестая (на 200, 210, 240 и 270-е сутки) под влиянием нисходящей ветки графика уровенного режима половодья рек; седьмая (на 360-е сутки) - под влиянием уровенного режима реки, предшествующего зимним попускам в ней.

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

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

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

Каждый атлас содержит карты абсолютных отметок и величин подъемов уровней грунтовых вод в каждом блоке модели для варианта в условиях НПГ - 68 м на 90, 180, 270, 360, 450, 540, 630 и 720-е сутки, для остальных вариантов на 140, 190, 200, 210, 240, 270 и 360-е сутки. Кроме того, приведенные в модели данные трансформированы в виде карты максимальных прогнозных уровней грунтовых вод в аллювиальных отложениях (рис. 6) и сопровождаются обобщенным представлением о динамике балансовых составляющих прогнозируемых уровенных режимов. Для каждого временного периода воспроизведенного на ЭВМ варианта процесса подтопления обобщения содержат суммарную величину утечек водонесущих коммуникаций, потерь из поверхностных водотоков и водоемов, разгрузки подземных вод в поверхностные водотоки, перетока из верхнего водоносного горизонта в нижний и наоборот, сработки или пополнения емкостных запасов верхнего и нижнего водоносного горизонта, испарения, выраженных в м³/сут.

Рис. 6. Схематическая карта максимальных прогнозных уровней грунтовых вод в аллювиальных отложениях

1 - изолиния максимальных прогнозных уровней грунтовых вод в абсолютных отметках, м; 2 - модельная аппроксимация внешних плановых границ области фильтрации

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

уровенный режим грунтовых и подземных вод, формирующихся при НПГ, зимних попусках, 1 %-ной обеспеченности уровенного режима реки, 50 %-ной обеспеченности режима атмосферных осадков и соответствующих параметрах испарения в условиях из раздельного и совместного влияния с учетом непроницаемых боковых границ;

изменение максимальных положений уровня грунтовых вод во времени и по площади с отражением волнового характера его режима и продолжительности каждого периода;

наиболее неблагоприятный уровенный режим грунтовых вод;

участки затопления;

степень раздельного влияния каждого режимообразующего фактора на формирование уровенного режима грунтовых вод;

динамику изменения балансовых составляющих для всех вышеприведенных прогнозируемых уровенных режимов грунтовых вод.