Задача Стефана о впитывании влаги в почву

ЗАДАЧА СТЕФАНА О ВПИТЫВАНИИ ВЛАГИ В ПОЧВУ

Одним из первых направлений количественного метода изучения движения влаги в почве, как известно, явилось установление зависимости скорости впитывания от времени. Существует много различных формул, выражающих закон инфильтрации, которые основаны на различных моделях и подходах от эмпирических до рассмотрения и учета определенных физических свойств почв, краткий анализ которых проведем по данным [4]. В числе первых законов инфильтрации считается уравнение для скорости впитывания Грина и Эмта (1911). Исходя из предположения о наличии четко выраженной границы между увлажненной и неувлажненной частями почвы, они получили следующее соотношение для скорости u передвижения фронта увлажнения

u = A + B / L ,

где A = к - влагопроводность почвы при полном насыщении, L – глубина проникновения фронта увлажнения к моменту времени t, B = к( Н0 - Нб ), Н0 и Нб - гидростатический напор и давление барботирования, соответственно. Так как константы входящие в это соотношение можно определить только опытным путем, то полученное уравнение следует считать эмпирическим.

К чисто эмпирическим относятся также закон инфильтрации Костякова (1932)

u =uн, и уравнение Хортона (1940)

u = uк - (uн - uк)exp(-a t ),

где параметры uн, uк, a и F должны определяться путем измерений скорости впитывания в различные моменты времени.

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

u =

где константы А и В имеют другой смысл, по сравнению с уравнением Грина-Эмта. Поскольку постоянные А и В выводятся из формулы профиля влажности, то уравнение Филипа справедливо только для начальных фаз впитывания. Однако, как отмечено в [4], если эти константы рассматривать как эмпирические постоянные, то формулу Филипа можно применять как для малых, так и больших значений времени.

Другое направление было связано с теоретическими и экспериментальными исследованиями (Дарси, Навье – Стокс, Пуазейль) факторов, влияющих на скорость потока воды в различных средах. Полученные результаты позволили Букингему и Ричардсу (1931) дать общую формулировку закона движения почвенной влаги в виде уравнения [1], которое представим здесь следующим образом

∂θ/∂t = Ñ(D*× (θ) Ñθ) - ∂k*×(θ)/∂z..

Здесь оператор Ñ = ∂ /∂t + j∂ /∂t + k∂ /∂t (остальные обозначения приведены ниже). Это выражение с различными граничными условиями позволяет получать количественную информацию о движении воды в почве. Однако часто данные можно получить только путем численного интегрирования уравнения Ричардса.

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

Одномерное уравнение переноса влаги во влажностях запишем в виде

∂θ/∂t = ∂(D*× (θ) ∂θ/∂z – k*× (θ))/∂z, (1)

где D*×(θ) = k*×(θ) dh /dθ, z - вертикальная координата, направленная вниз.

Для простоты дальнейшего рассмотрения, по аналогии с [5], предположим, что коэффициент влагопроводности к* и влажность q связаны с потенциалом влаги h (h ≤ 0) следующими зависимостями:

к*(h) = кs eah, q(h) = (qs - qr) eah + qr,

где qs , кs - соответственно влажность и влагопроводность при максимальном насыщении; qr - гигроскопическая влажность; a - константа. С учётом приведённых соотношений уравнение (1) будет выглядеть как

= D - к, (2)

где D = , k = = .

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

Пусть до начала увлажнения (t = 0) почвенная сре - 0

да занимала полупространство z > = 0 с постоянной от q1 z = z(t)

времени влажностью, которую в общем случае можно

представить в следующем виде q2

q2(z) = q2к – (q2к - q20) exp( - bz), z

где b = const > 0, q20, q2к - некоторые постоянные вели - Рис.1. Схема процесса

ны. увлажнения

С момента времени t > 0 на границе z = 0 поддерживается постоянная влажность q1(0,z) = g qs. В этом случае вблизи граничной поверхности возникает увлажненный слой, с влажностью q1(z, t), толщина которого с течением времени увеличивается (рис.1).

Фронт увлажнения z = z(t) в любой момент времени отделяет горизонт с влажностью q1 от соответствующего горизонта q2, двигаясь с некоторой скоростью u = dz/dt в направлении более сухого слоя с влажностью q2(z, t). Считая, что влажности в обеих частях почвенного профиля могут изменяться скачком, запишем уравнения влагопроводности (2) для двух областей в виде

= D - к, t > 0, 0 < z < z(t), (3)

= D, t > 0, z(t) < z < ¥.

Во втором уравнении (3) отсутствует конвективное движение влаги, из-за принятого допущения о первоначально сухой почве. Учитывая, что в начальный момент времени существует только жидкая фаза q2(z), для t = 0 имеем

q2( z,0) = q2(z), z> = 0.

Краевые условия сформулируем следующим образом:

1) на границах области

θ1 (0,t) = g θs, t > 0, z = 0;

θ2 (z, t) = θ2k t > 0, z → ∞

2) на фронте увлажнения z = z(t)

θ1│ z=x-0 = θ2½ z=x+0

Кроме того, предполагается, что на движущейся границе z = z(t)±0 изменение влажности происходит в результате диффузионных процессов. Отсюда можно получить второе граничное условие (так называемое условие Стефана), рассматривая уравнение баланса влаги, которое представим здесь таким образом

= . (3а)

Перейдём в уравнении (3) к безразмерным переменным x, t и введём новую функцию W1(x,t):

x = kz/ D и t = k2t/ D;

W1(x,t) = (gqs - q1(x,t))exp(- x/2 + t/4) / (gqs).

В результате вместо системы уравнений (3) получим

W1(x,t)/∂t = ∂2W1(x,t)/∂x2, t > 0, 0 < x < ζ(t), (4)

∂W2(x,t)/∂t = ∂2W2(x,t)/∂x2, t > 0, ζ(t) < x < ∞,

где W2(x,t) ≡ θ2(x,t).

Начальное и краевые условия задачи (3) для полученных функций Wi(x,t) перепишем следующим образом:

W2(x,0) = q20, t = 0, x > 0;

W1(0,t) = 0, t > 0, x = 0;

W2(x,t) = θ2k, t > 0, x ® ∞.

Из первого соотношения на границе увлажнения x = z(t) имеем

W2(z(t),t) = g θs- g θs W1(ζ(t), t) exp (ζ(t)/2 - t /

С помощью автомодельной переменной μ = x /t (преобразование Больцмана), аналогично работе [2], приведем уравнения (4) к обыкновенным

дифференциальным уравнениям для функций Wί (μ):

d2Wί / dμ2 + (1/2) μ dWί /dμ = 0, ί = 1,2;

общее решение, которого записывается как

Wi (μ) = Aί+Bί Φ (μ / 2), (6)

где Ф(m/2) = (2 /)- интеграл ошибок; A и В - const.

Возвращаясь в полученном соотношении (6) к переменным x t и функциям W(x, t) и W2 (x, t), решение системы уравнений (4) запишем в виде

W(x, t) = A+BФ(x/(2)), t > 0, 0 < x < z(t) (7)

W2(x,t) = A2+B2 Ф(x/(2)), t >0, z(t)< x < ∞.

Из приведенных выше граничных условий для функций Wi (x, t), находим соотношения для постоянных Аi и Вi

А1 = 0 , А2 + В2 = q2к.

Из первого условия на фронте увлажнения (5) с учетом соотношений для констант имеем

A2+B2 Ф(z (t)/(2Ö t)) = gqs–B1gqsФ( z(t)/2Ö t)exp(z(t) / 2 -t /4).

Анализ этого выражения показывает, что для любого t > 0 условие на границе увлажнения может быть выполнено только в том случае, если В1 = 0, а аргумент функции ошибок не зависит от времени. Следовательно, необходимо принять, что соотношение z(t) / Ö t = const = σ.

Таким образом, с точностью до некоторой постоянной σ, закон движения фронта увлажнения описывается уравнением z(t) = σ Ö t , а его скорость

u (t) = dz / dt = σ / (2 Ö t ).

В первоначальных координатах эти соотношения имеют вид

z( t) = σ√D √t , u (t) = (σ√D t -1/2 )/ 2

или, учитывая выражения для k и D в уравнении (2),

z(t) = 2 β σ√t, u (t) = β σ t -½,

где β = (ks / 4α (θs-θr))1/2.

Вернемся к нахождению постоянных коэффициентов А2 и В2 . С учётом вышеизложенного анализа и полученных ранее соотношений следует, что

A2 + B2 = .θ2к,

A2 + B2 Ф(σ/2) = gqs;

откуда имеем

A2 = gqs (1 - θ2кФ(σ/2)/ gqs) / ω,

B2 = ( θ2к - gqs) / ω,

где введено обозначение ω = 1 – Ф(σ/2).

В результате, с учетом полученных соотношений, решение первоначальной задачи (3) запишем следующим образом:

θ1﴾z‚t﴿ = g θs‚ t > 0, 0 < z < z( t); (8)

θ2 ( z,t ) = θ2к (Ф (z/ 2√Dt) – Ф (σ/2) ) / ω + g θs (1- Ф ( z/ 2√Dt )) / ω, t > 0,

z( t) < z < ∞.

Точное выражение для скорости движения увлажненной границы можно получить, если определить значение постоянной σ. Для этого воспользуемся условием Стефана на фронте увлажнения для задачи (3). Так как dФ/dz = 2 exp(-z2)/√π, то после дифференцирования интеграла по параметру t в соотношении (3а) и последующих алгебраических преобразований, условие на границе будет выглядеть следующим образом:

-D (θ2к – gθs) exp(-σ2/4) / (ω√πD√t) = βσ(gθs - θ2к + (θ2к - θ20 ) exp( -b z(t)) ) ,

или после несложных преобразований

() exp (-σ ²/4) = σ(1 – Ф(σ/2)) , (9)

где = (θs – θ2к) ⁄ ma , ma = θs – θ2к + (θ2к - θ20) exp(- b z( t)).

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

Выводы

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

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

u(t) = ks + β σ t -½. (10)

2.Отметим, что полученные из решения задачи Стефана динамика изменения фронта увлажнения z(t) = 2β σ√t и скорость его движения в виде соотношения (10) зависят только от начального состояния и гидрофизических характеристик почвенной среды.

Библиографический список

1. Н, В, С, Шержуков и физико – химические свойства горных пород. М.: Недра. 1997. С.271.

2 , Малов уравнения математической физики. МГТУ им. . М.. 2002.С.368.

3. В, и др. Краевые задачи со свободной границей. Ташкент. 1988.

4. Под ред. Глобуса основы гидрогеологии почв. Л.: Гидрометеоиздат. 1973.С. 428.

5. J.-M. Chen, Y.-C. Tan, C.-H. Chen, and J.-Y. Parlange. Analitical solutions for linearized Richard’s equation with arbitrary time-dependent surface fluxes. Water Resour. Res., 37(4), , 2001.