РАСЧЕТ ПРОЦЕССОВ ПЕРЕНОСА ВЛАГИ, ТЕПЛА И СОЛЕЙ В НЕНАСЫЩЕННОЙ И НАСЫЩЕННОЙ ЗОНАХ ПОЧВОГРУНТА

Рекс Л.М., Якиревич А.М.

I. Постановка задачи

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

                                  (1)

                                  (2)

                         (3)

             (4)

где:

Н- обобщенный потенциал почвенной влаги (напор) ,м;

Н=P-Х, Р - капиллярный потенциал (всасывающее давление ),м;

Т-температура, 0с;

С- минерализация порового раствора, г/л;

-коэффициент влагопроводности, м/сут;

Дт()-коэффициент термодиффузии почвенной влаги, м2/сут 0с ;

lэф -эффективный коэффициент теплопроводности, Дж/м. 0с.сут;

Ср - объемная теплоемкость раствора, Дж/мс;

Д - коэффициент конвективной диффузии солей, м2/сут;

ДТтс – коэффициент термодиффузии солей,

W – объемная влажность, м33;

;

W  -объемное содержание незамерзающей влаги;

Л – объемная льдистость;

X -вертикальная координата (положительно ориентированная книзу),м;

-коэффициент свободной емкости,

Сэф - эффективная объемная теплоемкость грунта, кДж/м3.0С;

е - удельный отбор влаги корнями растений, сут-1;

-коэффициент, характеризующий поглощение солей корнями растений (0£dе£1 при dе=0 - соли не поглощаются);

- коэффициент, характеризующий замерзание солей(0£dл£1при dл =0 соли не замерзают);

b -коэффициент скорости растворения, сут-1;

Сн -концентрация предельного насыщения, г/л;

t -время, сут.

В случае использования модели солепереноса в среде с двойной пористостью(“сквозной” и “тупиковой” Coafs K.U. Smith В.Д. I964) уравнение (4) можно заменить системой уравнений:

    (5)

                             (6)

где:

С,N - минерализация порового раствора в "сквозных" и "тупиковых" порах соответственно;

d -коэффициент скорости массообмена.сут.-1;

c -доля “сквозных” пор в общей пористости, принятой зa I.

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

                           (7)

Первые два уравнения описывают передвижение ионов солей Na+ и Са+7 в почвенном растворе. Третье уравнение - изотерма ионного обмена (И.П.Айдаров и ;др.,1978), а последнее уравнение означает  постоянство содержания суммы ионов в почвенно-поглощающем комплексе.

С(х,t) - концентрации ионов в почвенном растворе,

Ni(x,t) - концентрации ионов в почвенно-поглощающем комплексе,

если , то g=1, если, то g - объемная масса почвы г/см3;

Дi -коэффициент конвективной диффузии;

de -коэффициенты, характеризующие поглощение ионов растением

i=1 и i=2 относятся соответственно к характеристикам ионов Na+ и Са++.

2. Начальные и краевые условия.

2.1. Начальные условия.

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

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

                                       (8)

где:

Р0(х), Т0(х), С0(х) - исходное распределение исследуемых компонентов. Значения W(х) можно задавать вместо Pо(х) только в том случае, если рассматриваемая область не содержит первоначально зон полностью насыщенных влагой. Если исходные значения потенциала, температуры и концентрации известны лишь в небольшом числе точек, то исходная информация дополняется путем интерполяции между известными значениями по законам ступенчатой или ломаной функции.

Примечание: Для решения систем (1)-(4) и (1),(2),(5),(6) составлена npoгpaмма на языке Фортран для ЭВМ Минск-32 и серии ЕС “Расчет тепловлагосолепереноса в почвогрунтах (одномерная задача)” (авторы Л.М.Рекс,А.М.Якиревич), а для системы (1),(2),(7) - ”Расчет влагосолепереноса в почвогрунтах) (конвективная диффузия двух-катионного обмена с учетом обменной адсорбции при неустановившемся движении растворителя” (авторы: И.П. Айдаров, И.Г. Глобенко, Лю.М. Рекс, А.М. Якиревич).

2.2.Краевые условия: для уравнения влагопереноса.

При решении задач влагопереноса можно использовать краевые условия 1-го, 2-го и 3-го родов.

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

                        

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

                    (11)

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

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

            Краевые условия второго рода задают потоки воды через границы области:

                           

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

Комбинируя выражения (9), (10), (12), (13),можно записать краевые условие для уравнения влагопереноса в более общем виде:

                   (14)

                   (15)

При  и  получаем условия I-го рода, а при   и  - условие
2-го рода на верхней и нижней границах области. Возможна также их комбинация.

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

На верхней границе, в этом случае, значение испарения можно задать, зависящим от влажности (а следовательно и от потенциала), по какой-либо эмпирической формуле, например, по А.И.Будаговскому (1964):

                   (16)

где:

Е0 - испаряемость (функция влажности, температуры приземного слоя воздуха), средняя влажность корнеобитаемого слоя,

WP - влажность разрыва капилляров,

Wз - влажность завядания.

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

Величина потока в этом случае определяется как:

                          (17)

где: 

h(tio) - напор над дреной на середине междренья в момент времени tio , определяемый решением уравнения Буссинеска; d - водоотдача. В  частности, для почвогрунта однородного сложения, используют – модифицированную зависимость С.Ф.Аверьянова, получим:

                         (18)

где:

hio-напор над дреной на середине междренья, определяемый решением уравнения влагопереноса на каждом шаге по времени:

   

где:

Тd - глубина от дрены до водоупора;

В -междренное расстояние;

d - диаметр дрены.

При наличии напорного питания или в случае слоистой толщи величина  должна определяться по формулам, отражающим  эти условия.

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

а) происходит сток воды с поверхности почвы, условие 2-го рода заменяется на Н(0,t)=0     (19)

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

                       (20)

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

               (21)

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

2.3. Краевые условия для уравнения теплопереноса.

Для решения задачи прогнозирования теплового режима почвогрунта могут быть использованы краевые условия 1-го, 2-го, 3-го родов.

Краевые условия 1-го рода задаются, если известен временной ход изменения температуры почвы на поверхности или на нижней границе области:

              

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

Если задача решается в периодическом варианте и первая гармоника преобладает над остальными, то можно принять:

где:

w - частота колебаний, равная 2П сут-1;

f0  - начальная фаза;

Т1 и Т2 – амплитудные значения температуры.

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

                      (25)

Здесь Т1 представляет собой алгебраическую сумму радиационного, турбулентного, конвективного, транспирационного потоков, а также источников тепла, если они имеются на деятельной поверхности.

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

                                 (26)

Краевое условие третьего рода с учетом конвективного переноса тепла движущейся влагой его можно записать в виде:

                                        (27)

где:

ТЭ(t) – температура воздуха;

 -коэффициент теплообмена.

Здесь функции ТЭ(t) и   представляет собой некоторую эквивалентную температуру среды и некоторый обобщенный коэффициент теплообмена соответственно. Разработана методика определения этих величин (Д,А.Куртенер, ,А.Ф. Чудновский, 1979).

Объединяя различные виды краевых условий, получим:

       

Комбинируя параметрами (0 или i) и  (0, 1 или коэффициент теплообмена), имеем различные краевые условия на верхней и нижней границах области ().

2.4. Краевые условия для уравнения солепереноса.

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

               (30)

                                (31)

Варьируя параметрами bci можем получать различные краевые условия. При этом рекомендуется на верхнее границе использовать условие З-го (bci=1) рода, а  на нижней границе 2-го рода (bc2=1; С2=0) или 1-го рода (bc2=1,С2 -задает минерализацию грунтовых вод).

3.Параметры уравнений переноса влаги, тепла и солей в почвогрунтах.

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

3.1.Параметры уравнения влагопереноса.

3.1.1.3ависимость потенциала почвенной влаги от влажности и льдистости.

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

В зависимости от степени увлажнения почвогрунта разными авторами предлагается ряд методик для определения зависимости ( ) при положительных температурах. Для определения зависимости  (  ) используются пластинчатые или мембранные прессы, а также тензиометры в зависимости от интервала изменения потенциала (A.M.Глобус,1969). Косвенными методами являются измерения криоскопическим методом (И, И.Судницин,1979) или методом центрифугирования (А.А.Роде, 1965). В КГУ разработан экспериментальный метод определения зависимости  (  ) в нестационарных условиях (И.Е.Жернов, Н.В.Дзекунов, Б.А.Файбишенко, 1979). Анализируя кривые зависимости  ( ) можно выделить два вида (рис.1). Кривая I характерна для песчанистых пород. Кривая 2 характерна для глинистых и торфяных почвогрунтов.

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

Так кривая типа I может быть хорошо описана формулами вида (В. В. Ведерников,
Б. Ф.Никитенков, 1976):

                  при                        

где:

hk максимальная высота капиллярного поднятия;

m – пористость;

W0 -максимальная гигроскопичность;

Рм -потенциал влаги при  W=W0 (Р= -500м.вд.ст.).

Аналогичный вид имеет зависимость

      g<1                          (34)

где параметры h*k и g находятся по методу наименьших квадратов при

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

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

        (И.С. Пашковский, 1973),            (35)

 


Рис.1 -Зависимости капиллярного потенциала от влажности:

1.-для песчанистых грунтов;

2.-для глинистых грунтов.

Рис. 2       Аппроксимация зависимости P(w) |

-   1

I           _„:...

•1 - суглинок тяжелей (Кисловская орос, система), 2 - песок тонкозернистый (Энгельская орос. система, эксперим. данные В.М.Яшина.ВНИИГиМ).


или

                    (А.М. Якиревич)                 (36)

W0 – влажность прекращения движения влаги в жидком виде;

h*k  в формуле (35) приведенная высота капиллярного поднятия;

h*k  в формуле(36) при W0 ,соответствующей максимальной молекулярной влагоемкости, соответствует максимальной высоте капиллярного поднятия, а при W0 равной максимальной гигроскопичности h*k   примерно в 1,5-2 раза больше высоты капиллярного  поднятия.

Влияние  защемленного  воздуха на процессы инфильтрации учитывается заменой пористости в формулах (32) - (36) на величину полной влагоемкости. На рис.2 приведены примеры аппроксимации опытных данных формулой (36) для различных почвогрунтов.

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

3.1.2.Зависимость коэффициента влагопроводности от влажности и температуры.

Методики определения коэффициента влагопроводности излагаются в значительном числе работ (А.М.Глобус, 1969, И.И. Судницын,1979,В. Кулик, 1979). Определение этого параметра проводят как в стационарных, так и в нестационарных условиях (экспресс-методы).

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

                                   (37)

где:

Ко- коэффициент фильтрации;

W0 - влагоемкость прекращения движения влаги в жидком виде;

m – пористость, п=3¸6.

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

формула недостаточно хорошо аппроксимирует коэффициент влагопроводности для большинства почв. Для тонкодисперсных и торфяных почв формула (37) дает завышенные (до 1  3порядков) значения коэффициента влагопроводности, особенно при высокой влажности (И.С. Пашковский, 1973; Л.Г.П.пов, 1980).

Более гибкой в этом отношении является зависимость (А.М.Якиревич):

            n*<1;              (38)

По предварительным данным для суглинков 10-4£ n*£ 5´10-3,

для cупесей 5´10-3£ n*£ 5´10-3, для песков 10-2£ n*£1.

Уменьшением n* можно резко уменьшить К(W) при небольшом изменении влажности от полного насыщения. На рис. 3 приведены семейства кривых K(W), вычисленных по
формулам (36) и (37).

Для одновременного определения зависимостей Р(W) и К(W) можно рекомендовать следующий метод (Л.М. Рекс,В.М. Яшин, А.И.Якиревич 1980) с использованием калилляриметра. Исследовались образцы почвогрунта ненарушенной структуры высотой 5см., отобранные на участке поверхностного полива в совхозе “Мелиоратор” (Кисловская оросительная система в Волгоградском Заволжье).

Образцы насыщались капиллярно до постоянного веса и устанавливались в капилляриметре. Под образцами создавалось разрежение и на каждой ступени наблюдалась динамика вытекания веды (рис.4). В результате были получены кривые Р(W ) и кривые  изменения количества отсасываемой воды со временем на каждой ступени разрежения. После этого был поставлен численный эксперимент по вычислению этого количества

воды путем решения уравнения влагопереноса с краевыми условиями вида:

                    (39)

где:

*   -высота исследуемого образца (0,05 м), Р(t);

Рi(t) - давление, задаваемое на данной ступени разрежения. 

Кривые P(W) аппроксимировались зависимостью  (Зб). При  известном Ко,

варьируя показателями n и n* в формулах (37) и (38), и решая поставленную задачу, нами были получены теоретические кривые количества вытекшей влаги на разных ступенях разрежения, которые сравнивались с экспериментальными. Использование формулы С.Ф. Аверьянова не дает хорошего совпадения на всех ступенях разрежения –6 и  -8 м.вод.ст. (рис.4). Дальнейшее увеличение n , приводит к совпадению теоретических и экспериментальных кривых на начальных ступенях  разрежения: -0,5 и  -1 м.вод.ст.,, однако при этом теоретические кривые для разрежений  -2, -6, -8 м.вод.ст. ложатся ниже экспериментальных точек. Напротив, формула (38) дает достаточно хорошее совпадение на всех ступенях разрежения (при n*=0,0002). При положительных температурах коэффициент фильтрации вычисляется но формуле Пуазейля:

                          (40)

В области отрицательных температур коэффициент фильтрации можно определять по формуле В Д. Комарова (1975);

                                         (41)


Рис.3 Кривые зависимости K(W),вычисленные по формулам а) С(37) , б) (38 ).


pиc. 4    Зависимости от времени объема воды, вытекшего из

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

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

-—— - вычисленные на ЭВМ с использованием формул С. Ф. Аверьянова
(I –n=6, 2-
n= 8 ) и (38)(n*=0,0002), К.=-0.32 м/сут. (Данные В.М. Яншина)


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

3.1.3..Зависимость  интенсивности транспирации влаги растением от влажности.

В ряде работ рассмотрены самостоятельные более или менее сложные математические модели поглощения влаги растением. У А.И. Голованова (1975) величина транспирации принимается прямо пропорционально влажности почвы:

                            (42)

                   (43)

где:

Wi - средняя влажность в слое мощностью hi;

Wз —влажность завядания;

 мощность всей корнеобитаемой зоны;

ек -суммарная скорость отбора влаги, м/сут,

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

Такой процесс можно описать функцией вида (Л.М. Рекс):

           (44)

где:

Wз и WППВ – влажности, соответствующие влажности завядания и предельной полевой влагоемкости соответственно.

            Распределение корневой системы по глубине почвы можно задать функцией:

                       (45)

В зависимости от значений параметров распределения (45} Х1, Xл, sл можно получить достаточно широкий спектр кривых (рис.5),аппроксимирующих корневую систему растения.

С учетом формул (44) и (45) получим функцию вида:

                 (46)

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

hmax - максимальная глубина распространения корневой системы на данный момент времени.

 

3.1.4. 3ависимость содержания не замерзшей влаги от температуры.

Согласно принципу равновесного состояния фаз воды в мерзлом грунте Н.А. Цытовича (1973), мерзлый грунт данного состава в широком интервале температур содержит не замерзшую воду, количество которой при постоянном внешнем давлении зависит только от температуры. Экспериментальную кривую W(Т) можно получить различными методами: калориметрическим, который считается стандартным в строительстве /21/, с использованием сублимационных камер (Э.Д. Ершов и др. ,1975) и другими. Д.М.Андерсон и А.Р. Тайс,(1973) предлагают следующую зависимость между содержанием не замерзшей воды и температурой:

                (47)

где  и bw  постоянные для данной почвы коэффициенты.

Авторами выведено уравнение для определения весового содержания не замерзшей воды:

                   (48)

где:

S - удельная поверхность грунта, м2/г;

а=0,2618,  в=0,5519,  с=-1,449,  d=-0,264

Приводя формулу (48) к антилогарифмическому виду, получим для объемного содержания не замерзшей воды:

                   (49)

где g -объемная масса,

W - в объемных %.

С помощью зависимости (47) проведена аппроксимация экспериментальных кривых W(Т), полученных разными авторами. Для представленных на рис.6 кривых получены следующие значения параметров  и bw (табл.1).

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

 

 


 


Рис, 6     Содержание не замерзшей воды в почвогрунтах в  зависимости от температуры:

а) и б) I-бетонит, 2-киевская глина,3-каолин,4-суглинок (Е.Л.Ершов и др.,1975г.)

в) 1-глина,содержащая  монтмориллонит, 2-глина,.З-суглинок, 4-супесь,5-песок (Н.Л.Цытович ,1973г.)

 


Таблица I. Значения параметров  и bw в формуле (47)

для различных почвогрунтов.

№№

Тип почвы.

bw

Рисунок.

1.

Бетонит

0,427

0,389

-0,197

-0,198

6-a

6-б

2.

Киевская глина

0,174

0,168

-0,308

-0,266

6-a

6-б

3.

Суглинок

0,210

0,109

-0,9

-0,847

6-a

6-б

4.

Супесь

0,077

0,076

-0,368

-0,345

6-a

6-б

5.

Глина, содержащая монтмориллонит

0,248

-0,229

6-в

6.

Глина

0,140

-0,166

6-в

7.

Суглинок

0,09

-0,164

6-в

8.

Супесь

0,045

-0,125

6-в

9.

Песок

0,0075

-0,180

6-в

 

3.15.Зависимость коэффициента термодиффузии почвенной влаги от влажности/

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

А.М.Глобусом (1969) установлено, что при влажности почвы большей максимально-гигроскопической коэффициент термодиффузии почвенной влаги варьирует в довольно узких пределах: 10-6 ¸ 10-7 см2/с0. Причем величина его существенно  зависит от влажности. Типичный вид зависимости изображен на рис.7. Кривые такого вида удобно аппроксимировать формулой (А.М.Якиревич):

                           (50)

где: DTO – максимальное значение  при влажности WT.

 


Рис. 7.     Зависимость коэффициента термодиффузии почвенной влаги от влажности:

—— экспериментальная (Д. Б. Циприс, В.Е. Щадилов, 1978)

--- теоретическая   

1-дерновоподзолистая тяжелосуглинистая почва (aТ=I79,56),

2-мощный чернозем (aТ =183,18).


3.2.Теплофизические характеристики почвогрунтов.

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

3.2.1.Объемная теплоемкость почвы.

Объемная теплоемкость Сск сухой почвы равна произведению удельной теплоемкости Суд на объемную массу почвы g :

                    (51)

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

Табл.2 Значение   (под.А.Куртнеру, Д.Ф.Чудновскому, 1979).

№№

Тип почвы.

Суд

1.

Обыкновенный чернозем

1,05

2.

Супесь

0,84

3.

Песок

0,75

4.

Суглинок

0,96

5.

Торф

2,18

6.

Каштановые

0,86

 

Формула для определения эффективной объемной теплоемкости имеет вид ( с учетом мерзлой фазы)

                          (52)

где:

Сл -объемная теплоемкость льда;

СР -объемная теплоемкость раствора;

L -скрытая теплота плавления льда.

Последний член в выражении (52) учитывает выделение-поглощение тепла при фазовом переходе вода-лед.

3.2.2.Коэффициент теплопроводности почвы.

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

                (53)

где: коэффициенты mi определяются из таблиц 3 и 4.

Таблица, 3 (Д.А.Куртенер, А.Ф.Чудновский ,1979).

№№

Тип почвы

m1

m2

m3

m4

Примечания

1.

Обыкновенный чернозем

-0,013

3,1

1,21

20

Cаратов

2.

Темно-каштановая

-0,017

 

2,2

1,90

18

Cаратов

3.

серозем

-0,0062

2,7

-0,20

18

Фрунзе

4.

Южный чернозем

-0,0104

2,4

0,68

20

Одесса

5.

Дерново-глеевая, подзолистая

-0,020

3,1

1,4

20

Рига

Таблица 4 (Д.А.Куртенер, А.Ф.Чудновский ,1979).

№№

Тип почвы

10-3m1

m2

m3

m4

Примечания

1.

Дерново-подзолистая, средне-суглинистая

0,5-0,0005

1,0

1,8

20

 

2.

Темно-серая тяжелосуглинистая

-1,9

1,5

0,6

30

 

3.

Светлосерая тяжелосуглинистая, оподзоленная

-4,6

2,2

0,2

19

 

4.

Дерново-среднеподзолистая легкосуглинистая

-4,7

2,1

0,7

28

 

5.

Дерново-карбонатная, тяжелосуглинистая

-0,9

1,6

0

33

 

6.

Подзолистая

-10

2,0

2,2

16

Почвы Латвии

7.

Окультуренная

-20

2,0

2,3

14

 

8.

Дернрово-карбонатная

-9

1,2

2,7

19

 

9.

Желтозем(поляна)

-0,8g

2,2

0

20

Почвы Азербайджана

10.

Желтозем (лес)

-4,6

2,0

0

20

 

11.

Коричневая

-5,6

2,7

0

27

 

12.

Серозем

-9,4

1,9

0

20

 

 

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

                                (54)

где: lТ -коэффициент теплопроводности талых пород;

lМ -коэффициент теплопроводности мерзлых пород.

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

Для песков:                          (55)

Для глин:                                   (56)

Wисх - выбранное (исходное) значение влажности, при которой определяются значения  и Dlисх. Последние могут быть определены или из опыта, или найдены из номограммы, представленным на рис.8. Н.С.Иванов (1969) рекомендует выбирать Wисх для песков 5%, для глинистых пород 10%.

3.3.Параметры уравнений солепереноса.

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

Величина коэффициента конвективной диффузии принимается пропорциональной влажности и скорости влагопереноса (А.И.Голованов,1975)

                           (57)

где:

DM -коэффициент молекулярной диффузии, м2/сут;

lС -параметр дисперсии пористой среды, м;

V -скорость влагопереноса;  nD 1¸2 ;

W -влажность (формула(57) ) в настоящее время является гипотетической) .


Рис.: 8 { Н.С.Иванов, 1969 )

а) Зависимость критерия для промерзших песчаных горных пород от влажности при различных значениях Wисх.

б) Зависимость величины Dl для промерзших глинистых горных пород от влажности при различных значениях  Wисх.

Примечание для пересчета Dlисх  в кдж/м. сут.град необходимо значение, снятое с графика умножить на 100.


Коэффициенты скорости растворения и обмена также принимается прямо пропорциоиальными влажности:

                 (58)

где:

b*-параметр, определенный при полном насыщении. При наличии градиента температуры имеет место дополнительный поток солей, величина которого пропорциональна этому градиенту (эффект Соре). Коэффициент пропорциональности Дтс –коэффициент термодиффузии солей. А его отношение к коэффициенту молекулярной диффузии растворенного вещества называется коэффициентом Соре. Определение его достаточно трудоемко и для мелиоративной практики не имеет смысла делать. По литературным данным d =10-310-5 для растворов в жидкостях (Г.Д.Рабинович и др.,1975). Если учесть, что коэффициент молекулярной диффузии колеблется в пределах 10¸10 м2/сут, то для коэффициента термодиффузии получим величину порядка

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

4. Численная реализация модели тепловлагосолепереноса на ЭЦВМ.

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

Область (0 £ X £, 0 £ t £То) разобьем сеткой. Шаг сетки по координате обозначим hx .Шаг сетки по времени обозначим hT. Величины hx и hT могут меняться, но как показывают расчеты с достаточной точностью можно пользоваться равномерной сеткой. В этом случае, число узлов cетки по времени обозначим Кт(КТ=То/hT),а по координате- Кх(Кх=+ I, значения hx выбираем кратными значению ). Для аппроксимации уравнении (9)-(14) используем неявные схемы (А.А,Самарский,1977), которые являются абсолютно устойчивыми и позволяют выбирать более крупный шаг по времени, что увеличивает скорость счета на ЭВМ. Система разностных уравнений решается методом немонотонной прогонки (А.Н.Боголюбов, В.И.Телегин, 1974) последовательно для каждой искомой функции (напора, температуры и концентрации).

Предложенные алгоритмы реализованы в виде программ, составленных на языке ФОРТРАН. Краткая типовая блок-схема представлена на рис.9. Программы организованы таким образом, чтобы уравнения, входящие в систему, можно было решать отдельно, вместе и в любой их комбинации. Для этого вводятся переменные целого типа, задающие режим счета:

I, влагоперенос рассчитывается; 0, влагоперенос не рассчитывается;

I, теплоперенос рассчитывается; 0,теплоперенос не рассчитывается;

I,  солеперенос рассчитывается, по уравнению(4), 0,солеперенос не рассчитывается;

2,солеперенос рассчитывается по системе (5)-(6).

Переменная, задающая режим счета MR.

МR=0. Расчет перераспределения компонент по глубине почвогрунта на один заданный период времени ТО. В режиме МR =0 могут подбираться naраметры уравнений путем сопоставления теоретических и экспериментальных эпюр влажности, температуры или солей. При этом для ряда значений параметров уравнений рассчитываются теоретические эпюры (исходные эпюры остаются прежними для каждых новых вводимых данных).

МR=1.Расчет режима. В ЭВМ вводится последовательность временных интервалов ТОi с заданным законом изменения краевых условий. На каждый новый расчетный период времени за исходные принимаются распределения компонент, полученные на конец предыдущего периода. На печать выдаются эпюры.

МR =2. Расчет выходных кривых с шагом ДТ по времени до времени ДР. Служит для подбора параметров по выходным кривым. На печать выдаются изменение компонент на любой заданной глубине ХKF - с заданным шагом. Расчет может повторяться для других параметров.

МR =3.Определение времени и нормы промывки. Расчет выполняется при МС=0 и задании нисходящих потоков влаги. Счет продолжается до тех пор, пока минерализация порового раствора в заданной точке XKF не станет меньше допустимой ДР. На печать выдаются значения компонент в этой точке с шагом ДТ по времени. После выполнения условия СKF<ДР происхо­дит возврат назад на один шаг ДТ, ДТ, уменьшается в 10 раз и расчет повторяется.

МR=4. Расчет режима. Аналогичен случаю МR=1, но на печать выдаются значения компонент только на заданной глубине ХKF. Помимо этого в программе предусмотрена возможность использования различных аппроксимирущих функций Р(W),К(W), e(W).

На рис.9 TV- текущее время, НТ- шаг сетки по времени. Переменная ДК предусматривает возможность автоматического расчета режима орошения (при LR =1)с учетом ограничений, налагаемых на водный и солевой режимы.

                   (59)

                    (60)

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

КХ- Кх: число узлов разностной сетки по координате, Кх£ 41;

КТ-Кт: число узлов разностной сетки по времени; ЕPS-критерии сходимости (рекомендуется принимать 0,001);


 


IRH-рабочий параметр (полагать при расчетах 0);

KIH-максимальное число итераций для уравнения влагопереноса (рекомендуется принимать 5¸10);

КRK-число вариантов счета в режимах МR=0;2;3(в остальных случаях КRК=1). Программа считает данный пример КRK раз, требуя для ввода перфокарты с данными, характеризующими водно-физические, теплофизические и гидрохимические свойства грунта.

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

КS -число слоев с различными свойствами в рассматриваемой области, KS <5;

IP-признак эпюры. При IР=0 исходные эпюры аппроксимируются на сетке ступенчатой функцией, а при IР=1-ломаной.

MR,MH,MT,MC-переменные, задающие режим счета (см.выше),

МД- переменная, задающая расчет оттока в дренаж.

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

НI -координаты, в которых определены исходные эпюры, М; НI (1)=0;

HI (К)-задает нижнюю границу области расчета.

НS-верхняя граница каждого слоя, м;

HS=0;

РО- m :пористость каждого слоя , доли един.;

GА- <g.-плотность каждого слоя (объемный вес грунта),кг/м3;

КДХ, КДР - целые числа, управляющие выводы информации на печать при MR=0 или I. Ha печать выводятся значения компонент в координатах

Х(I), Х(I+KДХ),.. ,Х(1+2КДХ).... Х(I+NКДХ), I+NКДХ £ КДР £ I¸(N+I)КДХ.

При КДХ=1 происходит вывод в каждом узле сетки от I-го узла с номером КДР.

ДТ - шаг по времени, с которым выводятся на печать выходные кривые при М=2 или 3, сут. (при  ДТ¹0);

ДР - конечное время, до которого считается выходная кривая при =2,-сут.

ДР - предельная концентрация при МR=3, г/л(или % если МН=0).

ДР=0 при МR=4.

КF- узел сетки, для которого считается выходная кривая при МR=2;3 или 4.

MHW=1; в качестве исходных вводятся значения P, м;

          =2; в качестве исходных вводятся значения W33;

MPS=1, зависимость P(W) по формулам (32),(33);

         =2, зависимость P(W) по формуле (34);

 

 

          =3, зависимость P(W) по формуле  (35);

          =4, зависимость P(W) по формулам  (36);

МКW=I, зависимость К( W) по формуле     (37);

         =2, зависимость K(W) по формуле  (38);

ME =l, зависимость e(W) по формулам    (42), (43);

       =2, зависимость е(W) по формулам   (44),(46);

Р - значения капиллярного потенциала в исходной эпюре, м;

W- значения объемной влажности в исходной эпюре, м33.

ТД - глубина от дрены до водоупора, м;

В -междренное расстояние, м;

ДДР  диаметр дрены, м;

НДР -глубина заложения дренажа, м;

Ко -коэффициент фильтрации, м/сут;

НК – hK или h*K в формулах (32)-(36), м;

WMMW* в  формулах (32)-(33) м33.,

WDW0 - h/o в формулах(33)-(38), м33 ;

PV – полная влагоемкость, м33;

PPV – предельно-полевая влагоемкость, м33;

WZ -влажность завядания; м33;

РМ - Рм в формуле (33) (м) или в формуле (34);

NKW - n или N* в формуле (37),(38);

ДTW –ДТО в формуле (50), м3/сут град;

АТ, WT - aT,WT в формуле (50),

Т - значения температуры в исходной эпюре, С0;

CS - -СOK объемная теплоемкость скелета, кдж/м3, град;

М1, М2, М3, М4 – m1, m2, m3, m4  в формуле (53)

NGP=1,если слой почвы глинистый,

         =0,если слой почвы песчанистый;

TKI= Dlисх, кдж/м.сут.град, если NGP=0;

       = ,если NGP=1; ; (формулы (55),(56))

WISWисх в формулах (55),(56),м33.

AW, BW,-aw,-bw в формуле (47) параметры кривой W (T))

С - значения минерализации порового раствора в "сквозных" порах , г/л;

N - значения минерализации порового раствора в "тупиковых" порах, г/л;


 


 

СА -c: доля "сквозных" пор в общей пористости, принятой за I (активная пористость

вычисляется в программе как ). Если МН=1 и МС=1,то СА=1.

NDnD в формуле (57);

CLD - lC: параметр дисперсии, м;

ДМ-ДМ : коэффициент молекулярной диффузии, м3/сут;

BE - b*: параметр скорости растворения или обмена, сут-1;

СН – СН :концентрация предельного насыщения, г/л.

ДТС-ДТС: :концентрация термодиффузии солей; м3 /сут.град.г/л; (в практических расчетах полагается равным нулю).

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

ТО -расчетный интервал времени, сут(одинаково набивается в перфокартах 21,25,26).

Н1-Н1, Н2-Н2, ВН1- bН1, ВН2- bН2, MV, МN -параметры, определяющие краевые условия для уравнения влагопереноса. Здесь могут быть следующие комбинации.

На верхней границе:

1)ВН1=0-условия I рода, HI-задает изменение обобщенного потенциала (Н) на верхней границе,

При MV=0 происходит скачкообразная смена предыдущего значения Нх =0 на значение HI(полив по полосам), которое держится в течение всего интервала ТО.

При М V =1 происходит линейное изменение потенциала на верхней границе от исходного до значения HI на момент времени TV =ТО.

При МV=5 происходит скачкообразная смена краевого(как и при MV=0), но в данном случае ТО задает поливную норму в м /га, а время полива TV вычисляется.

2) ВН1=1-условие 2-го рода, HI задает поток через верхнюю границу. Н1<0  задает интенсивность суммарного испарения,| м/сут, (при этом MV =0).

HI 0 задает интенсивность впитывания дождя, м/сут.

При MV=0 поток HI не меняется в течение времени ТО. Если поток на какой-то момент времени превысит впитывающую способность почвы, то при MV=2,происходит смена на краевое условие 1-го рода Нх=0=0;

При МV =3 условие второго рода меняется на (20) и расчет продолжается до момента ТО, т.е. моделируется накопление слоя воды на поверхности с учетом впитывания; при MV=4 происходит тоже, что и при МV=3, но расчет продолжается до тех пор, пока весь накопившийся за время ТО слой воды не впитается в почву.

На нижней границе:

1).ВН2=0-условме I рода, Н2 задает изменения обобщенного потенциала на нижней границе. При этом можно задавать МN=0 и МN=1 аналогично верхней границе.

2).ВН2=1-услсвие 2-го рода.

Н2<0 задает интенсивность оттока через нижнюю границу;

Н2>0 задает интенсивность напорного питания, м/сут;

При MN=0 значения Н2 остаются постоянными в течение интервала  ТО.

При MN=2 отток вычисляется по формуле (18) с учетом параметров дренажа; а Н2, в этом случае, задает интенсивность фильтрационных потерь( из каналов и т.д.) в м/сут с положительным знаком.

При MN=4  задается условие (21), а Н2=0.

КЕ=I, учитывается отбор влаги корнями растений,

      =0, отбор влаги корнями не учитывается.

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

Если КК=1,то программа выполняет решение задачи и выдачу результатов на печать (для MR=I или 4).

Если КК=0, то программа выходит на ввод перфокарт для следующего примера, начиная с
1-ой. Если за картой, где КК=0 лежит перфокарта, где в первых 6 позициях набиты нули, то программа заканчивает работу. Для
MR=0, 2или З программа выходит на решение следующего примера после того как просчитает данный пример КRK раз (вводя для каждого варианта карты с некоторыми из общих исходных данных (см. схему ввода с перфокарт).

KL=|0, если мерзлая фаза отсутствует,

KL=  I, если мерзлая фаза присутствует,

RL -признак расчета режима орошения,

Если RL=0,то расчета режима не производится.

Если RL=1, то режим рассчитывается. При этом программой может автоматически назначаться полив, если текущее время ТV<ТО и средняя влажность заданного слоя (ХМАХ, м) меньше WMIN или концентрация солей больше СМАХ. Краевые условия для уравнения влагопереноса при поливе задаются в 23пк (HI,BHI,MV ), ТР, СР -температура и минерализация поливной воды. Полив прекращается, когда в слое ХMAХ выполняется условие W>WMAX и С < СМIN . Затем счет продолжается с прежним краевым условием до момента TV= TO. После каждой смены краевых условий происходит выдача результатов счета на печать(с указанием величины поливной нормы при назначении полива).

ЕК - cуммарная скорость отбора влаги растением из корнеобитаемой  зоны, м/сут;

НМАХ - мощность корнеобитаемой зоны, м/сутГ

Sa - sл, XI--X1, XL-XA - параметры распределения корневой системы в формуле (45) (задаются при МЕ=2, при ME=1 полагаются равными 0). Поток через верхнюю границу определяется при задании суммарного потока влаги HI ,как разница (Н1-EK).

LЕ - переменная, определяющая расчет испарения.

LЕ=0 - испарение задается переменной HI.

LЕ>0может принимать значения при BH1=1. Возможны следующие варианты:

1) LЕ=1, испарение рассчитывается теплобалансовым методом:

Е=(R-РТ-ВТ)/25I,4,

задаются:

К- тепловой баланс, кдж/м2ж

РТ - турбулентный теплообмен, кДж/м2;

ВТ-поток тепла в почву, кДж/м2 (за период ТО)

2) Е=2,испарение рассчитывается биоклиматическим методом:

Е=ВК.SДЕF ,

задаются:

ВК - биоклиматический коэффициент,

SДЕF -сумма среднесуточных дефицитов влажноcти (мб) за период ТО.

3) E=З, испарение расcчитывается по формуле А.И. Будаговского (16),

задаются:

ЕО - интенсивность испаряемости; м/сут;

W1W3  -влажность завядания, м33;

W2Wр- влажность разрыва капиллярной связи, м33;

ХЕ - средняя мощность испаряющего слоя, м(-0,5¸1,25м).

Испаряемость может рассчитываться по формуле Н.Н .Иванова:

Ео=0,0006(ТВ32+25)(I00-A)

В этом случае в 24пк вместо Е0 набивается число +1.000+4 и задаются:

ТВЗ - температура воздуха, С0;

А - относительная влажность воздуха, %,

При LE>0 HI может задавать среднюю интенсивность дождя за период ТО.

Если МН=0, задается скорость фильтрации 0, м/сут.

TI-Т, ТА-Та, F0 – f0, Т2-Т2 AI-a1, A2-a2, BTI - bT1, BT2- bT2, MTV1HTV -параметры, определяющие краевые условия для уравнения теплопереноса в формулах (28),(29).

При ВТ1=0, ВТ2=0 ,А1=1, А2=1 имеем краевые условия 1-го рода. Если MTV =0 и MTN=0 то происходит скачкообразная смена температур на границах области. Если суточных колебаний температуры на верхней границе не задается, то ТА=0, 0=0. При МТ =1 и МТ =1 задается линейное изменение температуры от предыдущего значения до Т1иТ2 на верхней и нижней границах соответственно. Эти условия используются наиболее часто.

При BT2=1 и А2=0 имеем условие  на нижней границе.

СР-Сп, С2-С2, BCI - bci, BC2- bc2,, параметры, определяющие краевые условия для уравнения солепереноса (30),(31).

При ВС1=0 и ВС2=0 имеем краевые условия 2-го рода СР иС2 задают минерализацию поливных и грунтовых вод соответственно.

При ВС1=1 имеем условия 3-го рода, используемое наиболее часто,

СР - минерализация поливной воды, СР=0 в перерывах между поливами.

При ВС2=1 и С2=0 получаем условие  на нижней границе.

ВСЕ - dе, ВLл .

Все вводимые с перфокарт данные за исключением К, КS, KXДХ, ДТ, ДР, КДР, КF, КЕ выводятся на печать.

Результатом счета являются следующие массивы, выводимые на печать:

Для МR=0 или I:

Границы зон полного и неполного насыщения (где Р=0)

TV -текущее время, сут;

Транспирация, м3/га;

Эпюры:

Х- координата, м;

Н-  обобщенный потенциал, м;

Р- капиллярный потенциал, м;

W -влажность;

WN -содержание незамерзающей влаги,

WL-льдистость;

V -скорость влагопереноса, м/сут;

Q -поток влаги через сечение X за расчетный период времени, м3 /га;

С- концентрация солеи в "сквозных" порах;

N-концентрация солей в "тупиковых" порах;

S =cC+(1-c)N - средняя концентрация (для MC=I S=C).

SW-запасы влаги послойно в соответствии с координатой X, м3 /га;

SWN-запасы незамерзающей влаги;

SWN -зaпacы льда;

SC -запасы солей в "сквозных" порах, т/га (если концентрация задана в г/л);

SN -запасы солей в "тупиковых" порах;

SS-общие запаcы солей;

ДW, ДWN, ДWL, ДC ДN ДS   изменение соответствующих запасов по сравнению с предыдущими;

SSW, SSWN, SSWL, SSC, SSN, SSS    -суммарные запасы компонент в зоне от Х(1)до Х(KДР);

SДW, SДWN, SДWL, SДC, SДN, SДS-изменение суммарных запасов по сравнению с предыдущим периодом.

Для МR=2,3,4 на печать выводятся:

TV -текущее время;

Х(К )-координата точки, в которой считается выходная кривая; значения Н, Р, W, WN,WL,Q,C,N,S  в этой точке.

5.0. Рекомендации по проведению расчетов процессов тепловлагосолепереноса  на ЭВМ.

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

Расчетная глубина  почвогрунта выбирается из соображения возможности задания нижнего краевого условия и должна включать наиболее интересующую нас область, но не слишком превосходить ее, чтобы уменьшить время расчета. При близком залегании УГВ и наличии дренажа расчетная глубина принимается ниже уровня заложения дренажа. После этого область расчета разбивается равномерной сеткой с шагом hx путем задания числа узлов сетки Kx:hx=/(Kx-I) (первый узел располагается на верхней, а последний- на нижней границах области).

Рекомендуется число узлов выбирать таким образом, чтобы шаг сетки hx был порядка 10-20 см при  м и порядка 25-50 см при > Зм, т.е. Кх=П¸21. С увеличением Кх пропорционально растет время расчета. При этом желательно разбивку производить таким образом, чтобы границы слоев с разными водно-физическими свойствами попадали в точки сетки.

Затем готовятся сведения о воднофизических, теплофизичееких и гидрохимических свойствах почвогрунта.

Экспериментальные кривые Р(W) и К(W) аппроксимируются с помощью какой-либо из предположенных выше аналитических зависимостей, и по методу наименьших квадратов находятся необходимые параметры (hK, hK*, n, n*).

Если экспериментальных кривых P(W) и K(W) не имеется, то зависимость Р (W) принимается априори по одной из формул, с учетом имеющихся данных (высота капиллярного поднятия, пористость, полная влагоемкостъ, наименьшая влагоемкостъ, максимальная гигроскопичность и т.д.),  а параметры коэффициента влагопроводности K(W) (Ко, n или n* ) подбираются путем сопоставления расчетных и экспериментальных данных по перераспределению влаги в почве при поливе или испарении.

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

На рис. 11а и б проделан ряд , как влияет показатель степени n в формуле С.Ф. Аверьянова на процессы переноса. С увеличением n (уменьшением n* в формуле (38)) фронт становится более крутым. При задании краевого условия второго рода на поверхности почвогрунта, а также при процессах испарения наблюдается та же закономерность.

Установлено, что в условиях переменных по направлению потоков влаги термодиффузионный член мало влияет на влагоперенос (доля в общем потоке влаги не более 5%) и им можно пренебрегать в практических расчетах (ДТW=0).

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

После выбора параметров уравнений весь расчетный период времени разбивается на отдельные промежутки Тoi в течение которых процесс влагопероноса на верхней границе имеет одно направление (осадки, полив или испарение). При расчете каждый интервал Тoi, разбивается временной сеткой. Число узлов сетки по времени Кт. Рекомендуемое число Кт=20¸25. При этом его следует сопоставлять с величиной интервалов То и учесть при выборе последних следующее: желательно, чтобы максимальная величина шага сeтки по времени hTi=ToiТ  для процессов инфильтрации не превышала 0,025¸0,05сут и 0,5¸2 сут для процессов испарения. На каждом интервале Toi задаются краевые условия для уравнений (текущие исходные данные).

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

                 м/сут

где Ос- осадки, м;

И- суммарное испарение, м;

Поливы при расчете следует задавать отдельно. Краевое условие на верхней границе можно задавать 1-го и 2-го рода. Используя условие 2-го рода (ВН1=I, МV =2), величину потока желательно задавать не превышающим коэффициент фильтрации. Для этого искусственно увеличивается время полива. При поливе по полосам задается напор на поверхности почвы. Норма полива (TO,BHI=0,MV =5).При поливе по бороздам напор на верхней границе принимается Hj = -0,5м. Если испарение задается в течение длительного интервала времени (более 25сут) и интенсивность задаваемого потока более 0.0025 м/сут, то в расчете следует задавать его через отбор влаги корнями растений с интенсивностью порядка 80-90% от суммарного потока.


 


Для уравнения теплопереноса предпочтительно задавать краевые условия первого рода на обоих границах (bTi=0; aI=1) МТV = MTN=1нижней границе при значительном удалении от поверхности почвы (более Зм) можно задавать условие второго рода  . Для уравнения солепереноса на верхней границе обычно задается условие 3-го рода (Данквертса -Бренкера), причем для периодов испарения минерализации поливной воды полагается равной нулю.

На нижней границе ставится условие 2-го рода  или задается минерализация грунтовых вод.

Пример1. Определение коэффициента влагопроводности на основе данных по перераспределению влажности при поливе (Бардинская опытно-мелиоративная станция, данные АзНИИГиМ за 1977год). .

Литологическое строение зоны аэрации следующее: до глубины 1м -суглинок легкий, от 1до Зм - глина легкая светло-серая, 3-5м-песок, ниже глина. Уровень грунтовых вод находится на глубине 2,6 м. Данные по изменению влажности при поливе имеются на глубинах от 0 до 3 м через 25см. Исходя из этого, задаемся областью  =3м и шагом сетки hx =0,125 м, число узлов сетки Кх=25. В качестве исходной эпюры задается влажность. Область делим на два слоя с границей на 1 м. Поскольку опытных данных по зависимости потенциала от влажности не имеется, мы эту характеристику принимаем априори по формуле (36) ( с учетом осреднения по глубине воднофизических свойств имеем следующие параметры.

1 слой: hK* =2,2м; Ко=1,5м/сут; m=0,503;. ); Wп=0,450; W0=0,07

2 слой: hK* =2,6м; Ко=1,0м/сут;/; m =0,460; Wп =0,445; W0=0,07

Зависимость K(W) принималась по формуле (38).

На рис.12а показано перераспределение влаги после 1-го полива нормой Q=1760м3 /га и продолжительностью t =1,1сут.

На поверхности почвы задаем краевое условие 2-го рода задающее поток q =0,16м/сут  (bHi=1; =0,I6) со сменой на краевое условие 1-го рода, задающее напор Н=0 в случае полного насыщения верхнего слоя почвы. На нижней границе задаем условие 1-го рода уровень грунтовых вод (bHi=0; H2=-2,6), MV=2.

Для разных значений n получаем различные расчетные эпюры и сопоставляем их с экспериментальными. При n*=0,01 расчетная эпюра ложится левеe экспериментальной (расчетный объем впитавшейся воды меньше экспериментального), следовательно K(W) следует увеличить путем увеличения n*.

При n* =0,05 расчетная эпюра значительно положе экспериментальной, а глубина промачивания больше натурной, следовательно, 0,01 <n*<0,05.

При n* =0,03 имеем неплохое совпадение зпюр, но расчетная глубина промачивания больше натурной.  Изменяя теперь для второго слоя n* подбираем его, уменьшая таким образом, чтобы совпадение улучшилось (n*2 =0,009) .Полученные параметры используем при расчете второго

полива (нормой 1920м3/га, продолжительностью 1,2сут). Результаты расчета дают хорошее совпадение (Рис.126).

Пример 2. Расчет водно-солевого и теплового режима почвогрунтов (режим орошения задан).

Расчет выполнен для темно-каштановых почв Поволжья. Исходный УГВ на глубине 10м. Высота капиллярного поднятия 2,0-2,5м. Известно, что на глубине 4м влажность в течение года мало меняется (21-22% объема). Задаемся этой глубиной :Кх=17, hx=0,25см, Кт=20. Зависимость P(W) принимается по формуле (36), К(W)-(37). Толща делится на 2 слоя с параметрами:

1)пахотный: m =0,49; g= 1600 кг/м3; Ко=0,5м/сут; h*K=2M; W0=0,11; WП=0,47; WППВ=0,35; W3=0,18; Рм—500м; n =5; Дто=0; aT=0; WT=0; СCK=1376 кдж/м3 град; m1 =-0,017; m2=2,2; m3=1,9; m4=18; Dl=8,4кДж/м.сут.град; Wисх=0,1; аw=0,21; вw =-0,9; lc=2,1м;  nд=1; Дм=I0-4м2/cyт; b =0; Сн=0; Дтс=0 ;c=1.

2)нижний : m =0,470; Ко=0,05м/сут; h*K=3м; WП=0,45

остальные параметра те же.

Исходные эпюры представлены в таблице 5.

Таблица 5.

ХМ

0

0,5

1,0

1,5

2,0

4,0

W0

0,174

0,190

0,188

0,180

0,188

0,210

T0

18

16

15

14

13,5

13,0

C

1

3

7

10

15

16

 

Для уравнения влагопереноса на верхней границе задавался поток влаги ; MV =0, если Q< 0; MV =2 если Q>0;

На нижней границе- условие 1-го рода Н =-10м; bH2=0; МN=0 (W=0,210).

Для уравнения теплопереноса на верхней границе задавалась температура (bТi =0, а1=1, МТV=1), на нижней границе -  (bT2=1, a2=0}. Для уравнения солепереноса на верхней границе - условие Данквертса (bCi=1), на нижней границе - . Исходные для расчета данные по периодам в течение года представлены в таблице (культура - яровые зерновые ).

В результате на конец каждого периода получали распределение компонент по глубине. На рис.13 представлены кривые изменения компонент на глубине Х=0,5м. Время счета на ЭВМ МИНСК-32 50мин. Если нет необходимости в расчете какой-либо компонента (например, температуры) то переменная, задающая режим расчета этой компоненты

(МТ - для температуры), полагается равной нулю, и соответствующие исходные данные не вводятся.

            Если в вегетационный период требуется расчет режима орошения, то на этот период переменная aR=1 (пример дан в приложении 7).


Рис. 12 Перераспределение влажности после полива:

а) 1-й   полив,  Q =1760 м3/гa,

б) 2-й   полив, Q=1920 м3/га.

(Экспериментальные данные АзНИИГиМ. Бардинская опытно-мелиоративная станция).

 


 

Таблица 6.

Период

Дата

Т0, сут.

NП3/га

Ос, м3/га

И, м3/га

Ос+NП-И=S

м/сут.

СП,г/л

Тпов, С0

Н2, м

(нижнее краевое)

1. Осенний

1/IX

2/IX-30/X

0,4

60

1200

-

 

780

-

990

+1200

-210

+0,3

-0,00035

0,01

0

+19

+1

-10

-10

2. Осенне-зимний-весенний Снеготаяние

1/XI-31/I

 

1/II-31/III

1/IV-7/IV

92

 

60

7

-

 

-

-

-

-

1020

-

 

-

400

0

 

0

+620

0

 

0

+0,0014

0

 

0

0,01

-7

 

+1

+7

-10

 

-10

-10

3. Весенний

8/IV-18/IV

19/IV

19/IV-7/V

11

0,4

18,6

-

1200

-

-

-

340

158

-

400

-158

+1200

-60

-0,0014

+0,3

-0,00032

0

0,01

0

+14

+12

+18

-10

-10

-10

4.Вегетацион.

 

 

 

 

 

 

 

 

 

 

1.полив

  перерыв

8/V

8/V-19/V

0,167

11,833

500

-

-

170

-

423,8

+500

-253,8

+0,3

-0,00214

0,01

0

+15

+19

-10

-10

2.полив

  перерыв

20/V

20/V-1/VI

0,167

12,833

500

-

-

-

-

733,4

+500

-733,4

+0,3

-0,0057

0,01

0

+17

+21

-10

-10

3.полив

  перерыв

2/VI

2/VI-13/VI

0,233

11,767

700

-

-

60

-

786

+700

-726

+0,3

-0,00032

0,01

0

+20

+25

-10

-10

4.полив

  перерыв

  перерыв

14/V

14/VI-31/VII

1/VII-13/VIII

0,233

47,767

31

700

-

-

-

280

250

-

460

300

+700

-180

-50

+0,3

-0,00038

-0,00016

0,01

0

0

+22

+28

+18

-10

-10

-10

Примечание: Тпов.- температура поверхности почвы на конец периода ТО.

 

 


Cписок  литеpaтуры.

1.Аверьянов С.Ф. Зависимость водопроницаемости почвогрунтов от содержания в них: воздуха. Докл. АН СССР, 1949т. Ст. 69., №2,141-144

2.Айдаров И.П. .КлыковВ.Е. Пестов Л.Ф. ,Шулъгин Д.Ф. Математическая модель динамики ионов натрия и кальция в почвах. Почвоведение. 18,1978г.

3.Андерсон Д.М. , Taйc A.P. Незамерзшая вода в мерзлых грунтах. II междунар. конф. по мерзлотоведению. М., Наука, 1973., с 102

4.Богушевский А.А. Оценка типичности при выборе экспериментальных объектов осушения. В сб. "Режим осушения и методика полевых научных исследований. М., "Колос"., 1971г.

5.Боголюбов А.Н. Телегин В.И. Об одном численном методе решения линейных систем уравнений с трехдиагональной матрицей. Журн. вычисл.  матем. и матем. физики.,I974.,т.14. ,№8., 768-771.

6.Будаговскжй А.И. Испарение почвенной влаги М. .Наука, 1964., 244с.

7.Ведерников В.В. Никитенков Б.Ф. Исследование совместного движения воды и воздуха в почвогрунтах. Почвоведение, 1976. , №12, 73-78.

8.Глобус А,М Экспериментальная гидрофизика почв. Л. ,Гидрометеоиздат,1969., 355 с.

9.Голованов А.И. Прогноз водносолевого режима и расчет дренажа на орошаемых землях. Автореф. диссер. на соиск. уч. степ. докт. техн. наук М. ,1975, 32 с.

10.Долгов С.И. Исследование подвижности почвенной влаги и ее доступности для растений. II.-Л., АН СССР. ,1948., 205с.

11.П.Ершов Э.Д., Кучуков Э.З., Комаров И.А. Сублимация льда в дисперсных породах. М., Изд. МГУ. ,1975.

12.Жернов И.Е., Дзекунов Н.Е., Файбишенко Б.А. .Новое в исследовании влаго- и солепереноса в зоне аэрации. Сб. Инженерные изыскания в строительстве., Киев., 1979., 23 с.

1З.Иванов Н.С. Тепло- и массперенос в мерзлых горных породах. М. .Наука., 1969.,240 с.

14.Комаров В.Д. Исследование водопроводимости мерзлой почвы. Метеорология и гидрология.,1975.,Л.2.,10-13.

15.Кулик В.Я. Прикладные расчеты на ЭЦВМ влагопереноса в зоне аэрации М., Недра, 1979. 160с.