Файл: Чайлдс Э. Физические основы гидрологии почв.pdf

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

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

X — Xt'/!,

где X, заменившая %из уравнения (12.26), есть функция, зависящая только от с. Ее можно определить и представить в виде графической

очень мало, поскольку в начальной стадии градиент de]dz бесконечно велик, и потому членом К в уравнении (11.10) можно пренебречь по сравнению с Ddcjdz. Это уравнение такое же, как (11.9). Следова­ тельно, на самых ранних стадиях развитие вертикального профиля протекает так же, как и горизонтального профиля, и решение уравне­ ния имеет вид

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

Позднее, когда член Ddcjdz не слишком превышает К , этот по­ следний член из уравнения (1 1 .10) способствует тому, что продвиже­ ние вертикального профиля происходит быстрее, чем горизонтальное впитывание. Поэтому можно ожидать, что решение для этого случая будет состоять из уравнения (12.52) с включением добавочных членов, тем меньших, чем меньше t. Эти члены представляют собой возраста­ ющие степени t, и поскольку первый из них содержит t4 *, логично попробовать использовать возрастающие степени t'!>. Результатом и явится уравнение (12.51).

Уравнение (12.51) имеет очевидные недостатки при очень больших t и потому применимо лишь к ранним стадиям инфильтрации. Напри­ мер, при дифференцировании по t для нахождения скорости опуска­

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

сказываемая уравнением скорость опускания профиля не стабили­ зируется ни при каких влажностях, в отличие от того, что имеет место в действительности, а продолжает возрастать со временем и влажностью, причем форма профиля непрерывно изменяется. Однако в определенных пределах времени это уравнение дает хоро­ шее согласие с опытом, что было показано Янгсом [1761для несколь­ ких материалов. Его кривые воспроизведены на рис. 12.15.

Влажность,% от объема

Рис. 12.15. Стадии раз­ вития нисходящего про­ филя влажности в слан­ цевой пыли (а) и стек­ лянных шариках (б).

ки — экспериментальные точки.

12.8.Горизонтальное впитывание

внеоднородный материал

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

Если вариации ограничены линейными размерами геометрически подобных материалов, как в параграфе 1 0 .1 , а относительно зависи­ мостей между влагопроводностью и сосущей силой и между сосущей силой и влажностью сделаны определенные допущения, то, согласно Филипу [125], решение существует.

В параграфе 10.1 указывалось, что если К есть влагопроводность материала, обладающего единичным характеристическим размером, то влагопроводность материала с характеристическим размером N равна

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

При равных влажностях соответствующие поры заполняются, а радиусы кривизны менисков становятся пропорциональными ха­ рактеристическому размеру N. Поскольку, в соответствии с пара­ графом 8.2 и, в частности, е уравнениями (8.1а), (8.2а) и (8.3а), со­ сущая сила обратно пропорциональна радиусу кривизны, она в дан­ ном случае будет обратно пропорциональна и характеристическому размеру. Следовательно,

где Н — отрицательное давление (сосущая сила с обратным знаком) в материале с единичным характеристическим размером. Отсюда, сделав подстановки из уравнений (12.54) и (12.55), можно получить уравнение неразрывности, аналогичное уравнению (11.28):

где N — неопределенная функция х.

Гипотетическими зависимостями между К, Н и с, позволяющими продолжить анализ, служат следующие функции:

Подставляя уравнения (12.57) и (12.59) в уравнение (12.56), получаем

Поскольку N зависит только от х и не зависит от времени, диф­ ференцирование уравнения (12.55) по t приводит к результату

dH/dt^NdH/dt.

Подставляя это соотношение и само уравнение (12.55) снова в урав­ нение (12.60), получаем:

В качестве граничных условий примем, что до впитывания ко­ лонка находилась в статическом равновесии при постоянном отри­ цательном напоре Н 0 и что впитывание началось при изменении напора на поверхности до Н1, каким он здесь остается и далее. Иначе говоря,

Н = Н0 для всех значений х при f = 0,

(12.62)

# = # х для всех t приж = 0 .

Это уравнение имеет ту же форму, что и уравнение (12.25), причем коэффициент (Яд/ЯД11 играет роль зависящего от влажности коэф­ фициента диффузии D из уравнения (12.11), так как если указано значение П, то коэффициент тем самым определен. Следовательно, уравнение (12.68) можно, как и уравнение (12.25), решить с помощью подстановки Больцмана и вычислить зависимость %от П посредством метода итераций, описанного в параграфе 12.6. Ввиду очень жестких ограничений, накладываемых на условия применимости этого ре­ шения, мы не будем здесь рассматривать его более детально.

12.9. Профили влажности при нестационарной фильтрации в грунтовые воды

В параграфе 12.6 мы успешно применяли подстановку Больцмана, поскольку не накладывалось физических ограничений на предель­ ную глубину проникновения фронта увлажнения в почвенную толщу. Эквивалентность влияния на математическую модель бесконечно малого времени и бесконечно большой глубины проникновения в поч­ ву данной изоплеты (данного с) позволяло выразить как граничные условия, так и сами переменные через единую обобщенную перемен­ ную xt~4 K Когда почвенный профиль ограничен не только сверху— дневной поверхностью, но и снизу — например, уровнем грунтовых вод, такая эквивалентность больше не соблюдается, и подстановка Больцмана уже неприменима. Подобный случай имеет место, когда избыточно увлажненному профилю, возникающему, например, при поливе напуском, удается стать стационарным за счет дренажного стока в грунтовые воды. При осушении поверхности стационарность

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

Метод решения такой задачи предложен Дэем и Лэсеном [55] и состоит в непосредственном численном интегрировании уравнения неразрывности (11.30):

произвольное число п элементов (каждый толщиной А) и заменив уравнения (12.69) соответствующим конечно-разностным уравнением для отдельного элемента.

Следуя этим авторам, разметим границы между элементами в ин­ тервале от 0 до п, где 0 — уровень грунтовых вод, п — поверхность почвы. Индекс г относится к верхней поверхности r-го элемента, счи­ тая от уровня грунтовых вод. Так же обозначаются влагопровод­ ность и потенциал на этой поверхности, т. е. Кг — влагопроводность на разделе г, а Ф,. — потенциал на том же разделе.

Принимая, что элементы достаточно тонки, можно допустить, что изменения К и Ф в пределах каждого из них линейны, так что влаго­ проводность А г_і / 2 в среднем сечении элемента г есть среднее арифме­ тическое влагопроводностей на границах элемента, а градиент по­ тенциала есть разность граничных потенциалов, деленная на А — толщину элемента. Тогда

Kr-ч, == (Кг-i + K r)j2,

- # г - і) /А

и точно так же

*,♦*/, = {Кг+ К г+1)/2,

(дЩдг)„чг = (Нг+1- Н г)ІА.

Кроме того,

{д [К (dH/dz)]/dz}r= [Яг+Ѵі (дН/дг)г+Чг - к т.Чг (дН/дг)г_Чшу д

(дК/dz), = (K„4t - iCr-./,)/А = (Kr+l - кг.,)/2Д.

Подставляя эти выражения вместо частных производных в урав­ нение (12.69), получим после небольших преобразований

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

на поверхности Нп имеет заданную величину. Затем делают пробную

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

àQ/6t == К0(grad Ф)0.

Рис. 12.16. Ста­ дии развития про­ филя сосущей силы в вертикаль­ ной колонке, ко­ торая дренируется уровнем грунто­ вых вод от исход­ ного насыщен­ ного состояния.

а — эксперименталь­ ные кривые, б — тео­

ретические.

Д а в л ен и е, см вод. cm.

При этом определяется интервал времени ôt между двумя последо­ вательными состояниями профиля, после чего можно по величине бі и разности влагозапасов в конце и начале интервала времени вычислить dcjdt для любой части профиля.

На следующем этапе вычислений для двух верхних элементов с помощью уравнения (12.70) рассчитывают величину Нп_ 1на основе предположений о величинах Нп, Нп_2 и соответствующих им зна­ чений Кп, Кп_2, сп и с„_2. Полученную таким образом величину Нп_ г используют затем вместо предполагаемой ранее величины для вычисления Нп_ г с помощью второй пары элементов и уравнения (12.70) так же, как и ранее. Этот прием повторяют, пока не дойдут до уровня грунтовых вод, после чего возвращаются к поверхности

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

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

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