Математическое моделирование солевого режима почв с фрактальной структурой Текст научной статьи по специальности «Математика»

УДК 517.958:531.34 С.Ю. Беданокова

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ СОЛЕВОГО РЕЖИМА ПОЧВ С ФРАКТАЛЬНОЙ СТРУКТУРОЙ

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

1. Основные уравнения модели и определение начально-краевых условий. Значительный интерес представляет разработка математических моделей, учитывающих влияние фрактальной структуры почвы на их водный и солевой режимы. На важность установления функциональной связи между водным и солевым режимами почв в процессе совместного движения воды и солей при полном насыщении почвогрунтов с растворимыми солями обратили внимание многие исследователи [1, стр. 25].

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

ди д2и ди

— = - а— + Ь(ит - и), (1)

дЬ дх2 дх

где и(х, у) — концентрация с(х, Ь) [г/л] почвенного раствора в точке х почвенного слоя 0 < х < г мощности г в момент времени Ь > 0 [сут]; х — расстояние [м] от поверхности; а = со/т\ — фактическая скорость движения воды в порах грунта; с0 — постоянная скорость фильтрации [м/сут]; т1 — порозность; ит — предельная концентрация насыщения; Ь — коэффициент растворимости [1/сут]; Бк — коэффициент конвективной (фильтрационной) диффузии [м2/сут].

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

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

/ (и) = Ь

I

и-

ит — I и(х, Ь) йх о

(2)

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

Уравнение движения (1) не учитывает, что почвогрунт, как правило, имеет фрактальную структуру [3]. Учёт этого фактора принципиально меняет уравнение движения солей, превращая его в дифференциальное уравнение движения солей дробного порядка следующего вида:

щ = Б]?дахи($, Ь) - аО0а-п ди^ Ь) + ¥ [и], (3)

где — коэффициент фрактальной диффузии;

- п дпЩ, Ь)

дЪхи($’ Ь) = Бо- п ’ п - 1 < а ^ п’ п е 14

Бвх и д%х — операторы дробного интегро-дифференцирования порядка |в| по Риману—Лиувил-лю и по М. Капуто порядка а [4], а ¥ [и] = Ь (ит - и) или ¥ [и] = / (и).

При а = 2 и ¥ [и] = Ь (ит - и) уравнение (3) совпадает с моделью (1) движения почвенного раствора.

Как установлено в [2], «почвенные коллоидные структуры чаще всего представляют собой массовые фракталы». Показатель Порода, т. е. число х в соотношении ^I(к) ~ х^к, характеризующем зависимость интенсивности рассеяния нейтронов ^I(к) от передаваемого импульса к, при реализации метода малоуглового рассеяния нейтронов используют для нахождения фрактальной размерности объектов. Для массовых фракталов число х совпадает со значением фрактальной размерности ц = Б.

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

Анализ данных фрактальных размерностей зональных почв, полученных в работе [2] методом малоуглового рассеяния нейтронов, показывает, что Б удовлетворяет неравенству 2 < Б < 3,9.

Основная цель этой работы — качественный анализ уравнения в частных производных дробного порядка а как математической модели движения солей для всех а е ]п —1, п] (п е N1) при этом особый акцент сделан на случае, когда 2 < Б < 3,9.

При хорошо растворимых солях и малом их содержании в твёрдой фазе уравнения (1) и (3) можно заменить следующими уравнениями:

ди д2и ди

д = Бк дх2 - а а? (4)

^ = Бгаа,и((, () - а^. (5)

В дальнейшем будем предполагать, что реализована процедура обезразмеривания зависимых и независимых переменных, а также всех параметров исследуемых задач. Предполагается также известным распределение концентрации почвенного раствора в начальный момент времени Ь = 0:

и(х, 0) = т(х), 0 ^ х ^ г, (6)

и начальная эпюра солесодержания такова, что т(х) е С [0, г] п С2}0, г] в случае уравнения (4)

и т(х) е Сп]0, г], т(п\х) е Ь [0, г] в случае уравнения (5).

Наряду с условием Коши в задачах оптимального управления водно-солевыми режимами важную роль играет нелокальное условие вида

г

— ^ и(х, Ь) йх = 8(Ь), 0 ^ Ь ^ Т, (7)

0

где Т — расчётное время.

Функция 8 = 8(Ь) — «математическое ожидание» содержания почвенного раствора (солей) в слое мощности г. Предполагается, что функция 8(Ь) принадлежит классу С1 [0, Т] функций, непрерывно дифференцируемых на временном сегменте [0, Т].

В классе достаточно гладких решений и = и(х, Ь) из уравнения (3) в силу (7) имеем

, Бг г а а\ а п и(0, ґ) п „1 1 Г

= даги($, ґ) сіх-- Ба;пи(ї, ґ)------------гп-а + - ¥(и) сіх

г J г [ (п - а +1) г J

(8)

0

Из (8) в случае, когда а = 2 ¥ [и] = /(и) = Ъ[ит — 8(Ь)], непосредственным вычислением получаем

8’(Ь) = -у [их(г, Ь) — их(0, Ь)] — а [и(г, Ь) — и(0, Ь)] + Ъ[ит — 8(Ь)]. (9)

Определяющее уравнение (9) порождает нелокальные краевые условия:

и(г, Ь) — и(0, Ь) = Р0(Ь), 0 ^ Ь ^ Т,

их (г, Ь) — их (0, Ь) = Р1(Ь), 0 ^ Ь ^ Т. ( )

Пусть заданы функции Р0(Ь), Р1(Ь). Тогда среднее солесодержание 8(Ь) меняется со временем по закону

8'(Ь) + Ъ8(Ь) = Ъит + - [Б-^Р1(Ь) — аР0(Ь)]. (11)

г

Положим в уравнении (3) Р [и] = Ъ[ыт — 8(г)] и применим к обеим его частям оператор Б^—а частного дифференцирования дробного порядка п — а, в результате, принимая во внимание, что БП~—аВа—п = Б°х, где Б°х — единичный оператор, получим другую форму записи уравнения

движения

—Бгп—аи($, г) = Бг „ дг °х 1 дхп

дпи(х, г) ди(х, г)

Нетрудно показать, что

— а-

дхп

+ Ъ [ит — 8(г)] Б°;а1

(12)

р,п—а-, _ Б°х 1 =

Г (1 + а — п)

I

1 [ Б п-а1 йх =-------- ----------,

гЛ °х Г (2 + а — п)

(13)

где Г(г) — гамма—функция Эйлера.

Так как ° ^ п — а < 1, то по определению

°п—аи(^’ г) = ~оХ-Боха—1и($> г)

и, стало быть,

/

/

Бп:аи(г)йх = Бп;аи(г) =

— Г

-а — п) J

Г(1 + а — пи (Г — £,)п—а' о

(14)

Проинтегрируем обе части равенства (12) по х от 0 до г .В результате на основании (10), (13) и (14) будем иметь

д

-Бп;а—1и($, г) = Бї

дп 1и(г, г) дп 1и(°, г)

дх п—1

дх п—1

— аР°(г) + Ъ[ит ш] га—п+1. Г (2 + а — п)

(15)

Уравнение (15) порождает нелокальное условие

дп 1и

дх

п-1

д -1 и

дх

1

х=°

(16)

Если допустить, что водно-солевой режим почвенного слоя таков, что

Би({, г) = \а8(г),

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

Л

и(г, г) =-а— гп-1-а8(г), о ^ г ^ Т. (17)

Г ( п - а)

При соблюдении условий (16) и (17) уравнение движения (15) переписывается следующим образом:

Кб'(г) +

Ъ

Г (2 + а — п)

Г а—п+18(г ) = БҐР п—1(г) — аР° (г) +

Ъи

Г (2 + а — п)

т а—п+1

(18)

2. Установившийся модельный вариант распределения солей в почвенном слое. Стационарные варианты уравнения (4) и (5) соответственно имеют вид

а

и (х) — ши (х) = °, ш = —,

Бк

а . а

д0хи(О — Ша и (х) = °, Ша = —

Б

г

(19)

(20)

Из (19) следует, что градиент V(х) = и'(х) концентрации солей в почвенном слое толщины г меняется по следующему экспоненциальному закону:

v (х) = v (°) ехр(шх).

(21)

а-п

а—п

х

х=г

Уравнение (20) эквивалентно системе:

йп 1 V (Г)

Б <^—п----------Г — ша V (х) = °,

йГп—1 а

и'(х) = V(х), ° < х ^ г.

К обеим частям равенства (22) применим оператор Б. В результате получим

йп—1 V

йхп

~ї = Ша Б 0х V (Г), ° ^ п — а < 1.

Любое решение V = V(х) уравнения (24) при п 5 2 является решением уравнения

й

йх

йп—2 V

йхп—2

— ШаБ,

п-1-а

аБ°х

V (Г)

= °.

Рассмотрим алгоритм решения уравнения (25) при

2 < а ^ 3 (п = 3).

(22)

(23)

(24)

(25)

(26)

1. Условие (26) означает, что градиент концентрации солей должен удовлетворять уравне-

нию

й

— [V '(х) — ШаБ2—а V (Г)] =°.

йх

2. Уравнение (27) редуцируется к уравнению

где

V'(х) — ШаБ2—а V(Г) = С2,

С2 = ІІШ V'(х).

х°

3. Из равенства (28), после применения к обеим его частям оператора Б-1, имеем

причем

V(х) — ШаБ°—а V(Г) = С2х + С1,

С1 = V (°).

(27)

(28)

(29)

(30)

(31)

4. Через Уа V обозначим левую часть уравнения (30). На основании формулы Хилле—Тамар-кина [4], которая обращает интегральный оператор Вольтерра Уа, любое решение V(х) уравнения (30) имеет вид

й’ (* V(х) = йх (С1 + С2г)Еа—1 [Ша(х — г)а—1\ йг.

(32)

Здесь

ер №> А = £

° Г (А+І/р)

, Ецр [г]=Ер [г;1], р > °,

функция типа Миттаг—Леффлера и функция Миттаг—Леффлера. Прямые расчёты показывают, что

й

йх

х

/'

]

— Еа—1 [Ша (х — г) йг =£ ^——— — [ (х — г)( ) 1 йг =

Л 1=° Г 1 + (а — 1) 1)йхЗ

° 1 =° °

±(й* V

йх ° 1=° Г {1 + (а — 1) 1)

= Еа- 1 [Шах

а- 1і

х

г

х

х оо ]

й [ гЕа-1 ‘(х - г)а-1\йг = £ _±х2На-1)1 [О (1 - 0(а-1)^г =

йх] 1=0 Г 1 + (а- 1) п йх і

о 1 = v ' о

1

= Ё ‘Г’/,2 Vа -“:)) X 1+(а-1) 1 в(2,1 + а 1 -1) =

,• _о Г (1 + (а - 1) ])

]

= X£ х--/ = X £ = хВ№_и Кх°-1;21

]=о Г (,3 + ( а - 1)7) 7=0 Г [2 + ( а - 1)7)

Поэтому формулу (32) можно переписать в виде

V(х) = Ва-1 [ШаXа-1] С1+ хВ1/(а-1) [ыаXа-1;2] С2. (33)

Как замечено в работе [3], «в целом ряде случаев значение фрактальной размерности достигает 3, что является предельной величиной для фрактальных объектов», и «такое значение

может характеризовать переход от массовых к поверхностным фракталам» [6]. Величина 3

является предельной и для уравнений (20) и (22). При а = 3 эти уравнения записываются следующим образом:

и”'(х) -ш3и'(х) = 0; (34)

V"(х) - ш3 V(х) = 0. (35)

Единственное решение V(х) уравнения (35), удовлетворяющее условию Коши (29), (31), задаётся формулой

V(х) = С1сЬ[у1/Ш3х) + С2ЬЬ(-\/Шзх). (36)

В силу (23) и (36) решение и(х) задачи Коши:

и(0) = с0, и'(0) = с1, и''(0) = с2 (37)

для стационарного уравнения движения почвенного раствора (33) однозначно определяется формулой

и(х) = С0 + С 8Ь (уШЗ х) + С2 х) - С2) / л/ШЗ. (38)

В случае же уравнения (19) решение и(х) задачи Коши: и(0) = с0, и’(0) = С1 имеет вид

и(х) = с0 + — [ехр(ых) - 1] (39)

ш

Из (37)—(39) следует, что эти решения, будучи близкими на поверхности х = 0 почвы, существенно отличаются вдали от неё.

Приведём теперь схему решения задачи для уравнения (20) в общем случае, когда п > 3 и п- 1<а<п.

Пусть т = п - 1, в = п - а. Очевидно, 0 < в < 1 и уравнение (22) эквивалентно уравнению

°0х Игт -‘аV(х)=0, т > 2. (40)

^в dmv(О а о

Известно, что

т Л хв-7

Б-в Vт (О = К- V(О - £ и(т-1) (0) -

і=0 Г (1 + в - .

Поэтому уравнение (40) можно переписать в виде

V(О) -‘аV(х) = Н(х), (41)

где

Н(х) = £

1

Уравнение (41) было объектом исследования работы [6], и его решение, удовлетворяющее условию (42), можно выписать в явном виде через функции типа Миттаг—Леффлера.

3. Нестационарная математическая модель солепереноса. Здесь мы ограничимся рассмотрением модели, в основе которой лежит уравнение

8’(^ = Б]д%хи(О, ?) - аих, 1 < а < 2, (43)

с граничным условием

Б/их(0, () = ф(^, 0 ^ ^ Т. (44)

Решение уравнения (43) принимается за приближенное решение уравнения (5).

Пусть

ш(х, V)= Dfux. (45)

Тогда уравнение (43) можно переписать в виде

/у 1 8'и) а

Баа- ш (О, ^ - Лш (х, Ъ = -^, Л , (46)

БГ БГ

а граничное условие (44), в силу (45), можно записать в форме

ш(0, ^ = ф(А, 0 ^ t ^ Т, (47)

где ф(^ — непрерывная на сегменте [0, Т] функция.

Схема построения решения уравнения (46) с граничным условием (47) выглядит следующим образом.

1. Поскольку в = а - 1 е ]0,1],

хв-1

ш (°1, $=ш (х «- Гв х^Ш (О °

и предел в правой части равен нулю, то уравнение (46) можно записать в виде

(48)

о 8^ (t) хв

ш (х, « - ш (О я = Г ^ •

2. Из уравнения (48) по формуле Хиле—Тамаркина получим

! х

ш(х, V) = 8 () . — [ ^Вя Х(х - ^в dt. (49)

Г 1 + в в1

о

1

К обеим частям равенства (49) применим оператор Б01. Тогда согласно (45) будем иметь

ю^ы(х, г) = —г-—^[ гв Ев Л(х - г )в йг = —гвЕв Л(х - г)в * Г (1 + в і ^ Г (1 + р)}=0 г (1 + рі)J ^

8'(г) ~ л і

X

0 Г (1 + вр0 Г\1 + ,

0 ’’

йг.

Отсюда, принимая во внимание, что

X 1

I гв&-г)в1йг = xв(1+J)+1j' Ов(1 - О)в1йО = xр(1+і)+1B (в+1, + 1),

0

где В^, у) = IГ(XX+У)) — бета—функция, заключаем, что

. “ Л1 xв+1+в 1 . в+1 “ 1

Dfы(x, г) = 8 (г)^ —(--------) = 8 (г)Xе V —(----------).

* к Г (2+р+вї) к Г (2+в+вї)

X

Стало быть,

Dfu(x, t) = б'(t)xв+1E1/в Лxв; в + 2

3. Равенство (5G) приводит к уравнению

Df б(t ) = б'(t )afir

(5G)

(51)

л r

с коэффициентами ar = 1 f xв+1E1/в ^xj; в + 2\ dx.

Число

=

Л j

=o rГ [2 + в + в/')

r

/ОО

x в+1+в jdx = r в+1 У j =0

(Лг в) J

0 Г (3+в+в J)

- = rв+1E1/в ІЛгв; в + 3

r в+1

отлично от нуля при а 5 0. При а = 0 оно равно Г(3+в)'

Согласно (6), (7) условие Коши для уравнения (51) задаётся равенством

8 (0) = т,

(52)

где т = -/t(x) dx среднее значение концентрации почвенного раствора.

С учетом (52) из (51) получаем

б(^=тexp

J

rDf

, б (t) = —j- exp

afy

J

(53)

4. Подставляя найденное значение 8'(^ из (53) в равенство (50), приходим к эффективной, хорошо реализуемой на компьютере, формуле, определяющей распределение солей в почвенном слое мощности г:

тEl/вlЛxв;2 + в] <xв1

u(x t)=^^лJ3+J\ Ы exp-

Dft

гг+іEl/в \Лгв;3 + в}\

(54)

4. Нестационарный вариант модели (З). За приближенное решение u(x, t) уравнения (3) в случае, когда F[u] совпадает с функционалом f (u), определённым по формуле (2), примем решение уравнения

б'(t) + Ь6(t) = Df dOjxud, t) - aD0aX—nЩ(і, t) + bun, (55)

удовлетворяющее соответствующим начальному и краевому условиям типа (6), (16) и (17).

Исследование начально-краевых задач для уравнения (55) приводится по указанным ранее схемам, выражая функцию w = DfUx через б'(t) + Ь6(t).

В случае, когда задано нелокальное краевое условие

б(0 = тexp(-bt), 0 ^ t ^ T,

уравнение (55) переходит в уравнение

Df d^xUd, t) - aD-a0nщ(і, t) = -bum,

которое можно свести к виду уравнений, исследованных в работе [7].

Автор выражает благодарность профессору А.М. Нахушеву за постановку задач и ценные советы.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Аверьянов, C. Ф. Борьба с засолением орошаемых земель [Текст] | C. Ф. Аверьянов. — М.: Колос, 197B. — 2BB с.

2. Сербина, Л. И. Об одной математической модели переноса субстанции во фрактальных средах [Текст] |Л. И. Сер-бина || Мат. моделирование. — 2GG3. — Т. 15, № 9. — С. 17-28.

3. Фрактальные коллоидные структуры в почвах различной зональности [Текст] |Г.Н. Федотов, Ю. Д. Третьяков, В. К. Иванов и др. || Докл. РАН. — 2GG5. — Т. 4G5, № 3. — С. 351-354. — ISSN GB69-5652.

a

r

j

r

r

r

4. Нахушев, А. М. Дробное исчисление и его применение [Текст] |А. М. Нахушев. — М.: Физматлит, 2GG3. — 272 с.

5. Кольцова, Э. М. Нелинейная динамика термодинамики необратимых процессов в химии и химической технологии [Текст] |Ю.Д. Третьяков, Л.С. Гордеев, А. А. Вертегел.—М.: Химия, 2GG1. —4GB с.

6. Baret, J.N. Differential equation of non-integer [Text] |J.N. Baret || Canad. J. Math. — 1954. — Vol. б, No. 4.— P. 529-541.

7. Беданокова, С. Ю. Задача Коши и нелокальная краевая задача для обобщённых дробных осцилляционных уравнений [Текст] |С. Ю. Беданокова || Докл. Адыг. (Черкес.) Междунар. акад. наук. — 2GG5. — Т. B, № 1. — С. 9-15.

НИИ прикладной математики и автоматизации КБНЦ РАН, г. Нальчик рпИртаЗЗЗЭтаИ . сот

Поступила G7.11.2GG6