Диссертация (1151743), страница 12
Текст из файла (страница 12)
Подтверждением этого являются карты речного стока, построенные поофициальным данным Гидрометслужбы (Атлас…, 1986).В этих природных условиях орошение должно быть незначительным, слегкаувеличивающимувлажненностьиобязательнопроисходитьнафонеестественного или искусственного дренирования. Движение почвенной влагидолжно учитывать неполное насыщение в зоне аэрации, наличие гравитационных,каркасно-капиллярныхсил,атакжезависимостьвлагопроводностиотвлагонасыщения (С.Ф. Аверьянов, 1949).Основанием для модели явилось уравнение двумерного потока влаги,реализованное в виде конечно-разностной численной схемы (3.1).
Исследуемая63толща разбивалась на элементарные слоиh j (1 ≤ j ≤ Nx − 1) переменнойтолщины, от0,1 м вблизи поверхности до 1 м вблизи водоупора, h0 = hNx = 0 . Для учета размераплощади, обслуживаемой одной капельницей, рассматриваемый пласт шириной(нормально к плоскости чертежа) B разбивался вертикальными плоскостями дляобразования столбцов и расчетных блоков. Ширина этих блоков (по длинекатены) bi (1 ≤ i ≤ Ny − 1) принималась различной в зависимости от ее длины, приэтом b0 = bNy = 0.Конечно-разностный аналог дифференциального уравнения передвиженияпочвенной влаги и подземных вод по неявной схеме, исходя из баланса влаги вi, j блоке:Cwn +1i, jH in, +j 1 − H in, j∆t=H in, +j −11 − H in, jhjRвi , j −1−H in, +j 1 − H in, j +1hjRвi, j+H in−+11,j − H in, jbi Rгi −1, j−H in, +j 1 − H in+1, jbi Rгi, j− ein, j ;(3.1)В уравнении (3.1) Hin, +j 1 – напор, м, на расчетный момент времени n + 1; приотсчете напоров от поверхности землиH in, +j 1 = − xi , j + ψ in, +j 1 ;(3.2)ψ in, +j 1 – напор, м, эквивалентный каркасно-капиллярному давлению в зоненеполного насыщения (ψ < 0 ) и эквивалентный гидростатическому давлению взоне полного насыщения;Сwin, +j 1 – коэффициент влагоемкости, м3в/м4:n +1i, jCwω in, +j 1 − ω in, jω in, +j 1 − ω in, j∂ω;===∂H H in, +j 1 − H in, j ψ in, +j 1 − ψ in, j(3.3)где ωin, +j 1 , – объемная влажность почвы, м3в/м3 (м3в – кубический метр почвеннойвлаги).
При полном влагонасыщении Сw = 0. Связь между каркасно-капиллярнымпотенциалом и влажностью почвы принята в виде (3.4). ψω − ωМ= exp− ω= µhkm − ωМn;(3.4)где m – пористость, м3/м3; ωМ – максимальная гигроскопичность, м3/м3; hk –максимальная высота капиллярного поднятия, м; µ и n– коэффициенты,64зависящие от механического состава и структуры почвы, для суглинистых почвпринято µ = 1 , показатель степени n= 1.Коэффициент влагоемкости при n= 1 с учетом (3.3) равен:C (ω ) =∂ω ∂ω ω − ω М==.∂H ∂ψµhk(3.5)Riв, j - вертикальное сопротивление потоку влаги между центрами i, j и i, j + 1блоков, сут;Riв, j = 0,5(h j / Kωi , j + h j +1 / Kωi , j +1 ) ;(3.6)Riг, j - горизонтальное сопротивление потоку влаги между центрами i, j и i + 1, jблоков, сут;Riг, j = 0,5(bi / Kωi , j + bi +1 / Kωi +1, j ) ;(3.7)K (ω ) – коэффициент влагопроводности м в/м /сут, зависящий от коэффициента32фильтрации Kф; пористости почвы m; объемной влажности почвы ω; дляопределения коэффициента влагопроводности K (ω ) в зависимости от влажностипочвы пользуются формулой С.Ф.
Аверьянова:3, 5 ω − ω∗ ,K (ω ) = K ф * m −ω (3.8)где ω*- максимальная молекулярная влагоемкость или влажность разрывакапилляров, для более широкого диапазона влажности эту зависимость можнозаменить на ω − ωМK (ω ) = K ф m − ωМ5 .(3.9)Расходование влаги на испарение принято зависящим от погодных условийи от влажности почвы, оно разделялось на испарение с поверхности почвы,которое учитывалось как граничное условие, и на транспирацию, последняяраспределялась по корнеобитаемому слою пропорционально влажности почвы иплотности корней и входила в уравнение в виде интенсивности влагоотборакорнями растений из единичного объема почвы ei,j, м3в/м3/сут.
С этой целью длякаждой декады теплого периода по известным средней температуре воздуха T, оС65и относительной влажности воздуха a, % подсчитывалось потенциальное (приоптимальной влагообеспеченности) суммарное испарение (эвапотранспирация)Epot по формуле Н.Н. Иванова:E pot = 0,0061Kб (25 + T ) (1 − 0,01a ) ,мм/сут2(3.10)где Kб – биологический коэффициент, учитывающий особенности конкретногоценоза.
Потенциальнаяэвапотранспирация разделялась напотенциальноефиспарение с поверхности почвы E potи потенциальную транспирацию E tpotпропорционально затененности почвы растительным покровом fр, котораяфизменялась по декадам: E pot= (1 – fр) E pot и E tpot = fр E pot . Эти потенциальныевеличины испарения редуцировались на каждом временном шаге:фE ф = εE pot; ε = 2 w 0 − w 02 ; w0 =ωп −ω м;0,8 p − ω м(3.11)при влажности поверхностного 2…5 см слоя почвы, если ωп>0,8p,тоε' = 1; этизависимости согласуются, например, с исследованиями А.И. Будаговского.Фактическая транспирация редуцируется в зависимости от неоптимальностисредней влажности корнеобитаемого слоя почвы:E t = ε w E tpot , где ε w = 2 w k − w k2 , wk = ωoptk − ВЗ ;ω k − ВЗ(3.12)εw– коэффициент, учитывающий уменьшение транспирации при отклонениивлажности почвы от оптимальной, вид этой зависимости соответствуетисследованиямА.Р.Константинова(1968);ωk–средняявлажностькорнеобитаемого слоя почвы, переменная во времени; ω kopt - то же, оптимальная вданную декаду; ВЗ – влажность завядания.Скорректированная величина транспирации E t распределялась по глубинекаждого столбца в заданном корнеобитаемом слое пропорционально влажностипочвы и массе корней в виде интенсивности влагоотбора корнями растений изединичного объема почвы ei,j, м3в/м3/сут.Определение напоров почвенной влагиHin, +j 1с помощью системыалгебраических уравнений (3.1) представляет собой громоздкую вычислительную66задачу, т.к.
сводится к нахождению порядка 700 неизвестных (при принятойразбивке на блоки) с шагом около 1 суток на протяжении нескольких десятковлет. Следует также отметить существенную нелинейность этой системыуравнений, в которой емкостной коэффициент и проводимость существеннозависят от напоров почвенной влаги, следовательно, и от влажности почвы, чтотребует 3…7 итераций на каждом временном шаге.
Поэтому алгоритм решенияэтой системы должен быть наиболее эффективным. В настоящее время таковымявляется метод матричной прогонки, который введением вектора напоров по всемi – тым столбцам для каждого слоя j позволяет понизить размерность задачи доодномерной:n+1Uj = | H1n, +j 1 ; H 2n,+j1 ; H3n,+j1 ; . . . H Nr−1, j | при этом j=0,1,2,3...Nx.(3.13)С помощью этого вектора система уравнений (3.1) запишется в матричномвиде:AA jU j −1 − CC jU j + BB jU j +1 = F jгде -AA jиBB j(3.14)- квадратные диагональные матрицы размером ( Nr − 1) ⋅ ( Nr − 1) ,учитывающие вертикальные потоки влаги между i, j − 1 и i, j блоками и между i, j иi, j + 1CCW jблоками:A11AAj = 000A220010 и т.д., где Ai ,i =;h j Riв, j −1A33B11BB j = 000B220010 и т.д., где Bi ,i =;h j Riв, jB33- трехдиагональная матрица размером ( Nr − 1) ⋅ ( Nr − 1) , учитывающаягоризонтальные потоки влаги между i − 1, j и i, j блоками и между i, j и i + 1, jблоками, а также емкостной член:CC j = AA j + BB j + CCW j − DD jCw1100Cwin, +1jCCWj = 0Cw220 и т.д., где Cwi ,i =;∆t00Cw33(3.15)67DD j =D11D21D12D220D23000D32D33D3400D43D44и т.д., гдеD ii = − D i ,i −1 − D i ,i +1 ; Di ,i −1 =11; Di ,i +1 = ггbi Ri , jbi Ri −1, j(3.16)Левое граничное условие, т.е.
отсутствие потока в центре рассматриваемогопласта учитывается особыми правилами вычисления элементов этой матрицы D111.b1 R1г, jи D12; D11 = –D12; D12 =Аналогичное правое граничное условие учитывается при вычислениипоследних элементов этой матрицыD Nr −1 , Nr −1 = − D Nr −1 , Nr − 2;D Nr −1 , Nr − 2 =1гb Nr −1 R Nr−2, j.Вслучае, если в каком-то блокеi,j имеется источник или сток (дрена или канал), ониучитываются при вычислении соответствующих элементов матрицыDD j .Входящий в систему уравнений (3.1) вектор Fj объединяет все свободныечлены:Fj=| e1n, j −Сw1n, +j 1∆tH 1n, j ; e2n, j −Сw 2n,+j1∆tH 2n, j ; e3n, j −Сw 3n,+j1∆tnH 3n, j ; e Nr−1, j −n +1Сw Nr−1, j∆tnH Nr−1, j |.(3.17)При наличии источников или стоков на вертикальных границах или внутриобласти фильтрации, они учитываются при вычислении этого вектора.Решение системы матричных уравнений (3.14) ищется в виде рекуррентнойформулы:U j = PPj ⋅ U j+1 + QQj ;(3.18)для этого при прямой прогонке вычисляют матрицы прогоночных коэффициентовPPj и прогоночные векторы-столбцы QQj:PPj =|| CC j − AAj ⋅ PPj −1 || −1 ⋅BBj ;QQj =|| CC j − AAj ⋅ PPj −1 || −1 ⋅ || AAj ⋅ QQ j −1 − Fj || .В этих формулах|| a ik || −1(3.19)(3.20)обозначает обращенную матрицу.Матрицы PPj квадратные, размером ( Nr − 1) ⋅ ( Nr − 1) , их общее количестворавно Nx.