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

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

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

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

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



Рис. 2. Иерархическая структура интегральной трехмерной нестационарной математической модели ПС и управления УГВ и РУКС



Рис. 3. Расчетная схема ПС, состоящая из четырех проводящих каналов

Таблица 2

№ контура	Cоедин. каналов	Длина, м
1	1-2	1200
2	2-3	1800
3	2-4	2000
4	2-5	1600

1-ый контур обхода графа состоит из дуг: 1, 2 и вершин: 1, 2, 3;

2-ой контур обхода графа состоит из дуг: 1, 3 и вершин: 1, 2, 4;

3-ий контур обхода графа состоит из дуг: 1, 4 и вершин: 1, 2, 5.

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

(5)

где h, u -уровень и скорость движения воды в канале, м и м/с; u_q- скорость бокового притока, м/с; q - боковой приток, м ²/с; g - ускорение свободного падения, м/с ²; J₀ - продольный угол наклона дна: J_{f -} уклон трения; Q - расход воды в канале, м ³/с; F - площадь живого сечения, м ².

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

, (6)

где Н - уровень грунтовых вод, м; м - коэффициент водоотдачи; K_f - коэффициент фильтрации, м/с; о - функция источника (стока) влаги, м/с.

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

, (7)

где Q_др - расход воды в дренах, м ³/с; щ, d - площадь сечения и диаметр ДТ, м ² и м; Н ₀ - начальное значение УГВ, м; р - глубина закладки ДТ, м; л - коэффициент гидравлического сопротивления.

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

R = f (F, Qн, Ho, t, a, d, р, L1, L2, К_f, h_O, h_B, Ф), (8)

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

t = Т(F, Qн, Ho, a, d, р, L1, L2, К_f, h_O, h_B, Ф). (9)

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

Таким образом, в рамках интегральной трехмерной нестационарной математической модели ПС и управления РУКС система дифференциальных уравнений в частных производных (5) - (7) описывает пространственно-временные распределения уровня воды в СПОК и УГВ в горизонтальных плоскостях.

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

, (10)

- поток влаги м³/с³; л_x, л_y и л - коэффициенты влагопроводности почвы вдоль горизонтальной координатной плоскости x, y и вертикальной координатной оси z, м²/с; F - капиллярный потенциал почвенной влаги, м²/с ²; z - вертикальная координатная ось с началом координат на поверхности грунтовых вод и направлена вверх;

б = 0.02?G

- коэффициент скорости потери влаги в почве за счет транспирации воды корневой системой растений G.

Алгоритм численного решения трехмерных дифференциальных уравнений параболического типа является весьма трудоемким по затратам машинного времени. Однако, стратегия управления уровнем грунтовых вод, направленная на установление равномерного распределения УГВ в горизонтальной координатной плоскости Х, У, позволяет не учитывать горизонтальные переносы влаги от поверхности грунтовых вод, вызванные практически нулевыми горизонтальные градиентами влажности почвы. Указанное благоприятное обстоятельство РР УГВ, позволяет рассчитывать только вертикальный перенос капиллярного потенциала влаги от поверхности грунтовых вод, описываемый одномерным дифференциальным уравнением:

. (11)

Начальные условия. В начальный момент времени при t₀ =0 для системы нестационарных уравнений в частных производных (5) - (7) и (11) начальные условия задаются в виде:

h(x, t₀) = H(x,y,t₀) = 3; u(x, t₀) = Q (x,y,t₀) = 0. (12)

В качестве начальных условий для функции капиллярного потенциала влаги F, описывающего уравнением (11), задавалась линейная функция таким образом, что на поверхности почвы (z = H) F = 0, а на поверхности грунтовых вод (z = 0), значение функции F соответствовало значению потока влаги U_max.

Граничные условия. Для системы уравнений Сен-Венана (5), описывающих уровневый режим в СПОК на концах контуров ПС, задаются значения расхода воды Q_H(t) = Q_Н, равные производительности насосов Q_Н, если насосы (насос) отсутствуют, то естественно, Q_H(t) = 0.

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

; ; , (13)

где L1 - длина открытого канала и осушаемого массива, м; L2 - ширина осушаемого массива, м; L_HK-дополнительное фильтрационное сопротивление на контуре взаимодействия руслового потока с грунтовым потоком, м.

Для дифференциального уравнения (10), определяющего расход воды в дренах, краевые условия задаются только на левой границе для каждого ОМ в виде (см. подпункт 5.2.1), при УГВ выше уровня воды в проводящих каналах (процесс осушения мелиорированых земель)

> 0:

, (14)

< 0

. (15)

На поверхности почвы (z = H) граничное условие для капиллярного потенциала влаги F, задавалось в зависимости от значения потока суммарного испарения E и интенсивности осадков X. В данных расчетах поток Е принимался равным нулю и значение производной капиллярного потенциала влаги F по координате z, через поток влаги U определялся следующем образом:

, (16)

где Е- поток суммарного испарения, м³/с; Х - поток атмосферных осадков, м³/с.

Выражение для вычисления значения капиллярного потенциала влаги на поверхности грунтовых вод (z = 0) через частную производную

определяется как:

. (17)

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

Интерпретатор полученных результатов. Так, как интегрирование моделирующей системы дифференциальных уравнений проводится в локальных системах координат z_l (ЛСК l-го канала), в которых ось х_l направлена вдоль проводящего канала, то при анализе численных результатов расчетов, необходимо значения рассчитанных параметров представлять в общей системе координат. Вектор преобразований координат z_l из ЛСК в общую систему координат Z_i (ОСК) имеет вид:

Z_l = Z_0l + C(б_l)?DL_l, (18)

где Z0_l - вектор координат начала l-го проводящего канала в ОСК; C(б_l)-матрица направляющих косинусов; DL_l-вектор длин l-го проводящего канала в ЛСК; l = 0, …, k; k - число проводящих каналов ПС.

Вектор преобразования координат (19) в матричной форме записывается следующим образом:

, (19)

где d_{l -} длина l-го проводящего канала в ЛСК.

Таким образом, алгоритм вектора преобразований координат z_l ОМ и проводящих каналов из ЛСК в общую систему координат Z_l (ОСК), полностью определен соотношениями (18) - (19), а, следовательно, все искомые функции, вычисленные в узлах разностной сетки осушаемого массива, могут быть интерпретированы в ОСК.

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

Численное решение системы уравнений Сен-Венана, имеющих недивергентные члены, не может быть непосредственно построено из-за двух проблем. Запишем систему уравнений Сен-Венана (8) в следующей форме:

(20)

; ; ;

а = 4 м - ширина канала.

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

(21)

;

- скорость звука в воде;

; ; .

При реалистическом предположении u < c уравнения (21) имеют два семейства характеристик - с положительным тангенсом угла наклона характеристик (первое уравнение

dx/dt = u + c > 0)

и с отрицательным (второе уравнение

dx/dt = u - c < 0)

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

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

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

=; =, (22)

i=1,2, … N_х - 1; j= 0,1, … K - 1;

где N_х и K - число пространственных узлов по оси Х и временных слоев на разностной сетке; ф, h - шаги интегрирования по времени t и оси Х.

Система разностных уравнений (22) записывается в явном трехточечном виде, следующем образом:

; (23)

;

;;

=; .

Вектор искомых функций Т_i^j+1 определяется следующим рекуррентным соотношением:

; i = n - 1, …, 0; j =0, …, k - 1. (24)

Прогоночная матрица E_i и вектор W_i определяются в прямой прогонке (индекс i возрастает) по следующим формулам:

(25)

Обратная матрица в формулах (25) имеет вид:

;

?_i = D11_i D22_i - D12_i D21_i;

D21_i = B21_i; D22_i=B22_i; D11_i =B11_i- (C11_i E11_i+ C12_i E21_i);

D12_i = B12_i - (C11_i E12_i +C12_i E22_i).

Компоненты матрицы прогоночных коэффициентов Е_i+1 определяются как:

E11_i+1 = - A21_i D12_i/?_i; E12_i+1 = - A22_iD12_i/?_i;

E21_i+1 = A21_i D11_i/?_i; E22_i+1= = A22_iD11_i/?_i,

а компоненты вектора прогоночных векторов W_i+1:

W1_i+1 = (D22_i W1_i - D12_iW2_i)/ ?_i;

W2_i+1 = (D11_i W2_i - D12_i W1_i)/ ?_i, i = 0,…, N_x - 1.

Для начала расчета прогоночных коэффициентов (25) необходимо задать их начальные значения в точке i = 0 с учетом наклона характеристик системы уравнений (21) и значений расхода воды на левой границе контура обхода графа ПС (для случая, когда в начале контура обхода установлена насосная станция). Так как на левую границу контура обхода ПС приходят характеристики второго уравнения, что следует из записанных в каноническом виде уравнений (21), то к граничному условию присоединяется это уравнение, и система разрешаемых уравнений на левой границе имеет вид:

- F_0,j+1 u_0,j+1 = - Q_T; (26)

A21₀ F_1,j+1 +A22₀ u_1,j+1 - B21₀ F_0,j+1 - B22₀ u_0,j+1 = F2₀.

Или в векторном виде:

B₀ T_0,j+1 = A₀ T_1,j+1 + F₀. (27)

Рекуррентное соотношение (24) для точки i = 0 в векторной форме записывается в следующем виде:

T_0,j+1 = E₁ T_1,j+1 + W₁. (28)

Умножая слева соотношение (28) на обратную матрицу В ₀^-1 и сравнивая полученное выражение с формулой (26), определяем начальные значения прогоночных коэффициентов E₁ и W₁:

E11₁ = A21₀ / B21₀; E12₁ = A22₀ / B21₀; E21₁ = 0; E22₁ = 0; (29)

W1₁ = - Q_T B22₀ / (F_0,j+1 B21₀) + F2₀ / B21₀; W2₁ = Q_T / F_0,j+1.

Искомые функции F_i,j+1,u_i,j+1 в узлах разностной сетки рассчитываются обратной прогонкой (индекс i - уменьшается) на основе скалярных выражений:

F_i,j+1 = E11_i+1 F_i+1 + E12_i+1 u_i+1,j+1 + W1_i+1; i = N_x - 1,…, 0;

u_i,j+1 = E21_i+1 F_i+1 + E22_i+1 u_i+1,j+1 + W2_i+1. (30)

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

u_i+1,j+1 = (F1_i - D11_i F_i+1,j+1) / D12_i. (31)

Тогда прогоночные соотношения (30), после подстановки формулы (36), преобразуются в рекуррентное выражение:

F_{i, j+1} = (32)

u_{i, j+1} = ,

KF1_i+1 = E11_i+1 - E12_i+1 D11_i /D12_i;

KF2_i+1 = W1_i+1 + E12_i+1 F1_i /D12_i;

KF3_i+1 = F1_i / D12_i;

KF4_i+1 = D11_i /D12_i; s

- количество обходов открытых каналов. Далее по тексту N = N_x.

При наличии откачивающих насосных станций на правых концах контуров обходов графа ПС задаются значения расхода воды Q_H. В результате решения квадратного уравнения определяется выражение для расчета значений площади живого сечения F_N,j+1 и скорости движения воды в русле каналов u_N,j+1, которое образуется при умножении второго уравнения системы уравнений (31) на значения F_N,j+1:

u_N,j+1 = KF3_N / 2 - (KF3_N2 / 4 - KF4_N Q_H)^0.5;

F_N,j+1 = (KF3_N - u_N,j+1) / KF4_N. (33)

Возможность введения новой функции (21), определяющей значения скорости движения воды в руслах каналов, обуславливается семейством характеристик первого канонического уравнения системы уравнений (21) dx / dt = u + c, которые приходят на правую границу при переходе на следующий временной слой. Поэтому, при задании граничных условий на правой границе контура должно использоваться это уравнение и значение расхода воды, аналогично методике задания граничных условий на левой границе (см. систему уравнений (26)). Проведя аналогичные расчеты расчетам граничных параметров на левой границе, можно получить формулу (31). Имея рассчитанные значения F_N-1,j+1 и u_N-1,j+1 и используя их в качестве "очередных" граничных условий, рассчитываются параметры в следующей точке N-2 и т. д.

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

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

Моделирование сетей проводящих каналов ПС. Для водо-изолированных и горизонтально расположенных проводящих каналов (приток воды в открытые каналы равен нулю q = 0), имеющих прямоугольное сечение, система дифференциальных уравнений Сен-Венана (8), записывается в следующем виде:

(34)

Начальные условия для системы уравнений (34) задаются следующим образом:

F(x, 0) = 12 м ²; u(x, 0) = 0 м/с. (35)

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

Q(0, t) = 1.0 м ³/с, (36)

а на правой границе контуров обходов графа ПС c учетом непроницаемости торцевых стенок открытых каналов

Q(x_n, t) = 0.0 м ³/с. (37)

Интегрирование системы дифференциальных уравнений Сен-Венана (34) вдоль контуров обходов графа Г польдерных систем осуществляется на основе алгоритма численного решения систем гиперболических уравнений (22) - (33).

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

Проводящие каналы имеют разную длину (см. табл. 3), но одинаковую ширину, равную а = 4 м, при этом шаги интегрирования по координате h = 15 м и по времени ф = 120 с задаются одинаковыми для каждого канала. Это связано с необходимостью, с одинаковой точностью аппроксимировать систему дифференциальных уравнений конечно-разностными уравнениями, описывающими динамику воды в проводящих каналах различной длины. Действительно, точность аппроксимации дифференциальных уравнений полностью определяется шагами интегрирования по координате h и времени ф и выбранной разностной схемой. Сравнительно небольшое значение шага интегрирования по времени ф определяется не условием счетной устойчивости алгоритма, а характерными временами процесса переноса воды в открытых каналах в момент запуска насосной станции.

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

Таблица 3. Состав и структура ветвей направленного графа ПС

№	Cоедин. каналов	Длина, м
1	1-2	495
2	2-3	1000
3	3-4	1000
4	4-5	1000
5	5-6	1000
6	6-7	2500
7	5-8	1000
8	6-9	3900
9	4-10	3000
10	3-11	2500
11	11-13	2000
12	11-12	1500
13	12-14	1000

1-ый контур обхода графа состоит из дуг: 1, 2, 3, 4, 5, 6 и вершин: 1, 2, 3, 4, 5, 6, 7;

2-ой контур обхода графа состоит из дуг: 1, 2, 10, 11 и вершин: 1, 2, 3, 11, 13;

3-ий контур обхода графа состоит из дуг: 1, 2, 3, 9 и вершин: 1, 2, 3, 4, 10;

4-ый контур обхода графа состоит из дуг: 1, 2, 3, 4, 7 и вершин: 1, 2, 3, 4, 5, 8;

5-ый контур обхода графа состоит из дуг: 1, 2, 3, 4, 5, 8 и вершин: 1, 2, 3, 4, 5, 6, 9;

6-ой контур обхода графа состоит из дуг: 1, 2, 10, 12, 13 и вершин: 1, 2, 3, 11, 12, 14.

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

На рис. 5 приведены результаты численных расчетов пространственно-временных распределений площадей живого сечения F и потоков воды Q для открытых проводящих каналов ПС (время t приводится в часах, расстояния - в км). При этом использовались следующие обозначения: F1_i1, F2_i2, F3_i3, F4_i1, F5_i2, F6_i3 и Q1_i1, Q2_i2, Q3_i3, Q4_i1, Q5_i2, Q6_i3 - пространственные распределения площадей живого сечения F и потоков воды Q в каналах вдоль 1-го, 2-го, 3-го 4-го, 5-го и 6-го контуров обходов графа польдерной системы для момента времени t, равного одному часу; (F_j)_1,0, (Fj)_5,67, (Fj)_6,167, (Fj)_7,67, (Fj)_10,167, (Fj)_11,80, (Fj)_13,40 и (Qj)_1,0, (Qj)_5,67, (Qj)_6,80, (Qj)_7,35, (Qj)_4,65, (Qj)_11,80, (Qj)_13,40 - временные распределения площадей живого сечения F и потоков воды Q в каналах, для номеров пространственных узлов i, равных 0, 35, 40, 67, 80, 167, что соответствует значениям координат x_i - 0 км, 0.525 км, 0.6 км, 1.0 км, 1.2 км, 2.5 км для 1-го, 5-го, 6-го, 7-го, 10-го, 11-го и 13-го проводящего канала.

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

(Н= F/a)

относительно меньше, чем амплитуды колебания скорости воды в проводящих каналах (здесь и далее - (Q_j)_5,67 - значение потока воды Q на j-временном слое 5-го канала в 67-ом пространственном узле). Сдвиг фаз волновых процессов для различных проводящих каналов польдерной системы обуславливается временем "добегания" волны до соответствующего канала.

Отметим, что в точках ветвления проводящих каналов потоки воды Q делятся пропорционально длинам расходящихся каналов (см. рис. 5 с, Q1_i1- распределение потока воды Q вдоль 1-го контура обхода графа Г ПС), так как в более коротких проводящих каналах быстрее падает уровень воды и, соответственно, значения потоков имеют меньшие значения по абсолютной величине (знак минус указывает, что поток воды Q направлен в противоположную направлению оси Х сторону).



Рис. 5. Пространственно-временные распределения площадей живого сечения открытых каналов F, потоков Q и скоростей воды u для тринадцати проводящих каналов различной длины и шести контуров обходов графа ПС

Расчет параметров СПОК ПС с учетом рельефа местности. Для водо-изолированных проводящих каналов (приток воды в каналы равен нулю, q = 0), имеющих прямоугольное сечение

F = a · h,

где а - ширина проводящего канала), система уравнений Сен-Венана (8) записывается в следующем виде:

(38)

Для двух соединенных вместе проводящих каналов длиной l₁ и l₂, вертикальное сечение с обозначением начального h₀ и установившегося уровня воды h₁, h₂ (ломаная и сплошная прямая, соответственно), приведено на рис. 6. При этом продольный угол наклона дна первого проводящего канала равен J0₁ = 0, а второго -

J0₂ = Hyk₂ / l₂

(Hyk₂ и l₂ - высота уклона и длина второго проводящего канала).Как следует из геометрических построений, приведенных на рис. 6, задачу математического моделирования распределения уровней воды в сети проводящих каналов необходимо решать в два этапа в общей системе координат Н 0Х. На первом этапе, на основе численного интегрирования системы дифференциальных уравнений Сен-Венана (38), вычисляется установившийся уровень воды вдоль контуров обходов графа польдерной системы, исходя из начальных значений уровней воды h₀, заданных в локальных системах координат (ЛСК) с осями х ₁ и х ₂, направленными вдоль проводящих каналов и преобразования этих уровней в общею систему координат (ОСК) Н 0Х (пунктирная кривая на рис. 6). На втором этапе, с учетом углов наклонов J0₁ и J0₁ координатных осей х ₁ и х ₂ ЛСК, связанных с первым и вторым проводящим каналом к координатной оси Х ОСК, значения уровней воды h, вычисленные вдоль контуров обходов графа ПС, преобразуются из ОСК в ЛСК, связанными с проводящими каналами, в которых значения уровней воды равны h₁ и h₂.



Рис. 6. Установившийся уровень воды в двух взаимодействующих проводящих каналах (сплошная кривая), имеющих разные углы наклонов дна J0

Интегрирование системы дифференциальных уравнений Сен-Венана (38) вдоль контуров обходов графа Г польдерных систем осуществляется на основе алгоритма численного решения систем гиперболических уравнений (22) - (33).

Начальные условия для системы уравнений (38) задаются в виде:

h₁(x, 0) = 3 м; h₂(x, 0) = 3 + J0₂·x2 м; u(x, 0) = 0 м/с. (39)

Граничные условия задаются в виде потока воды Q, используя краевые условия 1-го рода. На левой границе контуров обходов графа ПС, при наличии насосной станции, задается ее производительность:

Q(0, t) = 1.0 м 3/с, (40)

а на правой границе контуров обходов графа ПС c учетом непроницаемости торцевых стенок открытых каналов:

Q(x_n, t) = 0.0 м 3/с. (41)

Результаты численных расчетов. Рассмотрим польдерную систему, состоящую из четырех открытых проводящих каналов одинаковой длины и ширины (а = 4 м), но с разными продольными углами наклона дна J0, имеющую три контура обходов графа (см. рис. 3 и табл. 2). Значения шагов интегрирования по координате х (h = 12 м) и времени t (ф = 120 с) задавались равными для каждого проводящего канала. По нелинейности и связанности системы дифференциальных уравнений Сен-Венана (38), проводились три итерации по всем контурам обхода графа ПС.

На рис. 7 приведены результаты численных расчетов пространственно-временных распределений площадей живого сечения F и потоков воды Q для открытых проводящих каналов ПС (время t приводится в часах, расстояния - в км). При этом использовались следующие обозначения: h1_i1, h2_i2, h3_i3 и Q1_i1, Q2_i2, Q3_i3 - пространственные распределения уровней h и потоков воды Q в каналах вдоль 1-го, 2-го и 3-го контуров обходов графа польдерной системы для момента времени t, равного двум часам; (Fj)_1,10, (Fj)_2,250, (Fj)_3,250, (Fj)_4,250 и (Qj)_1,10, (Qj)_2,10, (Qj)_3,10, (Qj)_4,10 - временные распределения площадей живого сечения F и потоков воды Q в каналах, для номеров пространственных узлов i, равных 10, 250, что соответствует значениям координат x_i - 0.12 км, 3 км для 1-го, 2-го, 3-го и 4-го проводящего канала.

Рис. 7. Пространственно-временные распределения площадей живого сечения F, уровня h и потоков воды Q для открытых проводящих каналов ПС, для случая а)

В результате процесса установления уровней воды h в двух соединенных вместе проводящих каналах (время установления около двух часов, см. рис. 7 b-d), имеющих разные продольные углы наклонов дна J0, установившиеся уровни воды h в ОСК (рис. 7 а), имеют разные линейные распределения в зависимости от значений продольных углов наклонов дна каналов J0. Наиболее интенсивные перераспределения уровней воды h в проводящих каналах происходит в течение первого часа процесса и носит волнообразный характер, о чем и свидетельствуют временные распределения площадей живого сечения F и потоков воды Q, приведенные на рис. 7 b-d. Установившиеся значения площадей живого сечения проводящих каналов F, выбранные в точках - х ₂₅₀ = 3 км, обратно пропорциональны продольным углам наклонов дна J0 (см. рис. 10 а). Абсолютные значения потоков воды Q для открытых проводящих каналов ПС при времени t ? 2 часов стремятся к нулю, что следует из рис. 10 с-d. Пространственные распределения уровней h и потоков воды Q имеют линейную зависимость, что характерно и для горизонтально расположенных проводящих каналов.

В пятой главе сформулирована и реализована задача двухмерного моделирования динамики уровневого режима грунтовых вод польдерных систем в горизонтальной плоскости, являющая составной частью ИТНММ ПС. Особое внимание при этом обращалось на разработку экономичных численных методов решения двумерного дифференциального уравнения Буссинеска и переноса воды в ДТ, увязки граничных условий на контуре взаимодействия руслов сети проводящих каналов и осушаемых массивов с учетом стока воды по ДТ, а также построению инвариантного алгоритма, синтезирующего результаты расчетов УГВ, полученные в ЛСК и преобразования их в ОСК.

Расчет УГВ с учетом дренажа и заданием граничных условий в области сопряжения проводящей сети с осушаемым массивом.

Граничные условия. Для вывода формулы, позволяющей рассчитывать значения потока воды в ДТ на левой границе у = 0, рассмотрим построения, приведенные на рис. 8.

Рис. 8. Схема системы проводящего канала и прилегающего к нему ОМ для вывода формулы значения потока воды в ДТ на левой границе y= 0 ОМ

Запишем уравнение Бернулли в вертикальных плоскостях Н ₀, d и h, пренебрегая членами уравнения с квадратами скоростей подъема (спуска) воды в канале и ОМ по сравнению с квадратом скорости воды в ДТ, следующим образом:

, (42)

где с - плотность воды, кг/м ³; V - скорость воды при выходе из дренажной трубы, м/с; р - глубина закладки ДТ, м;

h_p= h - z_p;

H_p= H - z_p; Н_р

- расстояние между УГВ и горизонтом закладки дрен.

Выражая явно скорость воды при выходе из дренажной трубы V и, учитывая, что поток воды в ДТ, равен Q = щV, получим соотношение для вычисления краевого условия на левой границе для дифференциального уравнения (6) в следующей форме, при УГВ выше уровня воды в проводящих каналах (процесс осушения мелиорированых земель)

. (43)

Уравнение движение сплошной среды по трубопроводу относительно потока жидкости Q записывается в виде:

, (44)

где Р - давление движущей среды; F - внешние массовые силы.

С учетом построений и обозначений, приведенных на рис. 11, дифференциальное уравнение (44), для переноса воды в ДТ, запишется в следующей форме:

. (45)

Алгоритм численного решения уравнения переноса воды в дренах. Запишем дифференциальное уравнение (45) в следующем виде:

, (46)

; .

C учетом знака Q > 0 разностная схема бегущего счета для дифференциального уравнения (46) запишется следующим образом:

, (47)

где m = 1, …, M; j = 0,…, k-1.

Умножая левую и правую часть разностного уравнения (47) на ф, приводя подобные члены, явно выразим искомую функцию Q_mj+1

,

; (48)

m = 1, …, M; j = 0,…, k-1,

Такая методика построения численного решения двухмерного уравнения Буссинеска подразумевает проведение итерационного процесса до получения необходимой точности.

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

Алгоритм численного решения двумерного дифференциального уравнения Буссинеска (6) реализовался на основе двухшаговой консервативной разностной схемы, с введением полуцелых временных слоев, обеспечивающей второй порядок точности по обеим переменным, т.е. 0(ф²+h²) и счетную устойчивость. Следуя приведенному алгоритму в гл. 2, дифференциальное уравнение (6) запишется в виде следующих двух разностных уравнений:

(49)

где h_y - шаг интегрирования по оси у; i = 1, 2, …, n - 1 - номер пространственного узла по оси х; m = 1,2,…,M-1 - номер пространственного узла по оси у; j = 0,…,k-1 - номер временного слоя; h_y - шаг по оси у.

Алгоритм численного решения двумерного уравнения (6) на разностной сетке строился на основе Т-образных разностных шаблонов (см. гл. 2, рис. 2.7-2.8, шаблоны выделены жирными линиями).

Моделирование УГВ реальных польдерных систем. Особенности расчета УГВ реальных польдерных систем связаны с необходимостью моделировать уровневый режим грунтовых вод ОМ в общей системе координат, математически увязывая модели расчетов УГВ для одиночных ОМ, выполненных в локальных системах координат с осью Х, направленной вдоль проводящего канала (см. рис. 3). На рис. 9 приведен фрагмент польдерной системы, состоящий из l-1, l и l+1 проводящих каналов и прилежащих к ним ОМ с указанием их длины: d_l-1, d_l и d_l+1 и ширины ОМ: L2_l-1, L2_l и L2_l+1.

Рис. 9. Фрагмент ПС, состоящий из l-1-го, l-го и l+1-го проводящих каналов, с указанием ЛСК и дренажных труб (выделенные прямые)

Для фрагмента польдерной системы, изображенного на рис. 9 синтезируем общую матрицу преобразований координат С(б_L) из ЛСК_L (L= l - 1, l, l +1) в ОСК по формулам (18) и (19), исходя из преобразований координат в следующем виде:

,

С(б_L)= (50)

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

Необходимо отметить, что приведенный выше алгоритм синтеза общей матрицы преобразований координат С(б_L) из ЛСК_L (L= l - 1, l, l +1) в ОСК, может легко быть распространен на ПС с произвольным числом проводящих каналов и ОМ, что позволяет говорить об инвариантности данного алгоритма.

Результаты численных расчетов. В качестве примера расчета основных параметров реальных ПС, включая расход воды в дренах, задавалась польдерная система, состоящая из трех последовательно соединенных проводящих каналов: длиной d_l-1 = 300 м (Nx_l-1 = 60), d_l = 100 м (Nx_l = 50), d_l+1 = 300 м (Nx_l+1 = 60) и шириной прилегающих к ним ОМ: L2_l-1 = 150 м (Ny_l-1 = 15), L2_l = 150 м (Ny_l = 15), L2_l-1 = 100 м (Ny_l-1 = 15,). Шаг интегрирования по времени составлял 3600 с и в момент времени, когда уровень воды в канале достигал 2.2 м насос выключался. Далее повышался до 2.4 м и насосная станция включалась (решались совместно система дифференциальных уравнений Сен-Венана для проводящего канала, двухмерное дифференциальное уравнение Буссинеска и дифференциальное уравнение переноса воды в дренах (5.10)), моделировалась закладка дрен диаметром d = 0.15 м на глубине р = 2.5 м, перпендикулярно руслу канала. При этом междренное расстояние составляло 25 м, и дрены, закладывались по всей длине осушаемого массива (канала).

На рис. 10 приведены результаты численных расчетов пространственных зависимостей УГВ в общей системе координат Х 0Y вдоль осей Х - (а), (с) и Y - (b), (d).

Рис. 10. Пространственные распределения УГВ вычисленные в ОСК вдоль осей Х - (а), (с) и Y - (b), (d)

Верхних два рисунка (см. рис. 10 а - b) подразумевают нижнюю часть рис. 9, до 3-го канала включительно и означают следующее: распределение УГВ вдоль оси координат х_l-1 на расстояние у ₂ = 20 м (а) и вдоль - у_l на расстоянии х ₅ = 25 м. Для верхней части рис. 9 приведено распределения УГВ вдоль оси Х ОСК, учитывающее влияние только 1-го и 3-го канала, отстоящее от нее на расстоянии Y = 110 м (с). На рис. 10 d указано распределения УГВ вдоль оси У ОСК, учитывающее влияние только 2-го и 3-го канала, отстоящее на расстоянии Х = 310 м. Необходимо отметить, что численные расчеты УГВ в ОСК дают адекватное представление о пространственно-временном распределении УГВ для всего ОМ, позволяющее оценивать качества ПС.

В шестой главе завершено решение задачи формирования и реализации ИТНММ ПС и управления РУКС почвы мелиорированных почв, подключением к горизонтальным процессам переноса влаги вертикального переноса влаги, моделирование которого основано на численном решении дифференциального уравнения капиллярного потенциала влаги (11) в зоне аэрации почвы. Заменяем исходное дифференциальное уравнение в частных производных параболического типа (11) конечно-разностным во внутренних узлах по консервативной схеме.

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

A_i F_{i+1,j+1 -} B_i F_i,j+1 + C_i F_i-1,j+1 = - D_i, (52)

; ; ; .

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

, (53)

а значит полученная консервативная дивергентная разностная схема (51), аппроксимирующая исходное дифференциальное уравнение (11), безусловно - устойчива и аппроксимирует его с первым порядком точности по времени и со вторым по координате 0(, h²). Решение уравнения для капиллярного потенциала (11), после его приведения к трехдиагональному виду (52), реализуется на основе метода прогонки.

Управление режимом увлажнения корнеобитаемого слоя. Рассмотрим алгоритм расчета оптимального уровня грунтовых вод, например, для условий глубокозалежных болот (многолетние травы, U₀ =86.5 мм). Оптимизируя целевую функцию потока влаги U(H, h) (см. раздел 6.3, (6.6)), при управление потоком влаги U_max > max и U_max ? U₀ и функциональных ограничениях 5 h 50 и 20 Н 60, вычисляем оптимальное значение уровня грунтовых вод Н_оп = 20 см, высоту корнеобитаемого слоя почвы h_оп = 14.9 см и максимально возможный поток влаги в зону аэрации U_max= 86.44 мм (см. рис. 11). По вычисленному значению максимально возможного потока влаги U_max от поверхности грунтовых вод, определяется краевое условие для капиллярного потенциала влаги F, при z = 0. На рис. 12 приведены значения капиллярного потенциала влаги F и потоков влаги U для трех моментов времени 1.5сут, 5сут и 10сут, рассчитанные на основе численного решение дифференциального уравнения (6.4) с начальными и граничными условиями (6.5) и (6.7) и оптимальном уровне грунтовых вод Н_оп.

...