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

Математические заметки СВФУ Апрель—июнь, 2017. Том 24, № 2

УДК 519.63

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ПРОМЕРЗАНИЯ ГРУНТА В. В. Попов

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

БОТ 10.25587/SVFU.2017.2.9248

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

Введение

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

1. Формулировка задачи

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

дТ д ( дТ\ Зу

= 0<Ж<^М>0. (1)

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

С = (1-т)Сара +тиСш +т(1 — v)Ci А = (1-т)\а +т(1-V

Работа выполнена при поддержке Российского фонда фундаментальных исследований (код проекта 17.01.00732).

© 2017 Попов В. В.

где т — пористость; Са,ра,Ха соответственно удельная теплоемкость, плотность, теплопроводность скелета пористой среды; те же теплофизи-ческие свойства для льда; С—то же самое для воды; Ст, Хт и С±, соответственно объемная теплоемкость и теплопроводность мерзлого и талого грунтов; дщ — теплота фазового перехода вода — лед; Т — температура.

В области фазовых переходов имеем уравнение диффузии

д^с д ( дс\

(2)

где Б — коэффициент молекулярной диффузии, с — концентрация, а также уравнение термодинамического равновесия фаз:

В талой области имеем

дТ

с = -Т/а, 0 <х<Б(г).

д2т

(3)

<9£ 0,2 дх2

В-

Б <х <1, г> о,

Б < х <1, г > о.

(4)

(5)

дс

дг дх2

На границе раздела талой и двухфазной областей х = Б (г) поддерживается неразрывность потоков тепла и массы и неразрывность искомых функций Т(х, г), с(х, г):

" дТ 1

о, т_ = т+ = т*, (6)

дх

дс дх

= 0, с_ = с+ = с* = — Т * /а. (7)

Система уравнений замыкается граничными и начальными условиями:

Т(о, г) = Тд, г> о, (8)

Т(1,г) = Тп, с(1,г) = со, г> о, (9)

Т(х, о) = Тп, с(х, о) = со, х е [о, г]. (1о)

2. Задача с точным решением температурного поля

Ниже будем сравнивать численное решение задачи (1)—(1о) (при малых временах, так как I — конечное число) с решением следующей задачи:

д2Т

ах

дх2 д2Т

дТ дТ _

<9£ 0,2 дх2 дгус д / дс дг дх дх

с = —Т/а, о <х< Б, г> о,

о <х <б, г> о,

Б < х < г> о, о<х< б, г> о,

(11)

(12)

(13)

(14)

дТ

О, ж = 5, Ь > О, (15)

дх

Т_ = Т+ = —ас*, с_ = с+ = с*, х = 4 > 0, (16)

т(0,г) = Тд, г> о, (17)

т(+<х>,г)= Тп, с(+<х>,г) = со, г> о, (18)

Т (х, 0) = Тп, с(х, 0) = со,, х > 0. (19)

При такой постановке данной задачи нет прямой зависимости температуры от водонасыщенности. Задачу (11)—(19) будем решать, используя автомодельную переменную £ = . Уравнение диффузии с учетом условия термодинамического равновесия выразим через температуру:

^Т) £ + Л / ¿п = 0) с > 0) (20)

¿£ 2 <£ \ ¿£ или в более компактном виде:

(иту^ + (0иту = 0, £>0. (21)

Раскроем скобки:

£ £

г/Г| + г/Г'! + БиТ" + Г>г/Г' = 0, £ > 0. (22)

Перепишем в виде

^(Т£ + 2БТ') = — v(T'£ + 2БТ"), £> 0. (23)

Сгруппируем по искомым функциям:

V' (Т '£ + 2ВТ'')

„ (Т£ + 2ПТГ ^

Левая часть данного уравнения — это производная от логарифма водонасыщен-

ности:

Т '£ + 2ПТ''

= С>0. (25)

Проинтегрируем левую и правую части уравнения от £о до в:

в в

( , ГТ '£ + 2ПТ''

■ +

€о €о

где £о — произвольно взятая точка между 0 и в, определяемая как

р VI'

Левая часть легко интегрируется:

в

( Т'£ + 2ВТ''

Ни{р)) - 1п(К£о)) = - У ТЦ2Е>Т, (27)

Так как на границе £ = в водонасыщенность равна начальному значению (v = vo), получим

в

С T+ 2DT''

Info) - МКСо)) = -j ТЦ2Ш, dt (28)

€о

Приходим к выражению для логарифма водонасыщенности:

в

С T+ 2DT''

ln(KCo))) = + J ТЦ2Ш, dt (29)

€о

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

v(£o) = voe1, (30)

где

в

I = J Tt + 2DT> * (31)

€о

Из формулы (30) видно, что если интеграл I расходится, то водонасыщенность в этой точке равна нулю, а если сходится, то водонасыщенность больше нуля.

Температура в протяженной области фазовых переходов вычисляется по формуле

T-Tg + (Tz-Tg)e^/2^_y (32)

В талой зоне

ч 1 - erf (£/2^) /

T = Tn + (T,-Tn)l (33)

Первая производная температуры равна

т' = —__Tz ~ Тд e~^2/4ai f34)

y^raT erf(/3/2^/ai) " 1 ;

Вторая производная температуры

у" = £__Tz-T9 c-S2/4ai^ ^gg^

2aiA/7rai erf(/3/2^/aT) '

где

_ (1~Д/01) TZ — Тд

v^I erf(/3/2v/ai)' Отсюда в числителе подынтегральной функции имеем

T + 2DT'' = ki£e-«2/4ai. (37)

Здесь

1 T — T

к2 = 2D 7-.. (38)

д/тгоТ eri(/3/2v/a~r)

Рис. 1. График функции /(£)

Знаменатель подынтегральной функции имеет вид

Т£ + 2БТ' = £Т + к2 е-5'/4а1.

Введем обозначения

Г (£) = Т£ + 2БТ' = £Т + к2 е-5'/4а1. £е-5'/4а1

р(£ ) =

/ (£)

Перепишем интеграл I в следующем виде:

в

I = к I р(£)<£, £ > 0.

5о

(39)

(40)

(41)

(42)

Проведем качественное исследование функции Г (£). Вторая производная

имеет вид

ПО = ^е

2

в

В

___^

а1 2а\)

£2

Заметим, что а1 больше чем Б. Видно, что в интервале

0 < £ < £ *

\

2

И

_1___

2а\ 2а?

(43)

(44)

будет Г''(£) > 0, т. е. кривая, определяемая этой функцией, вогнута вверх (рис. 1).

Допустим, что в интервале (0,£*) существует точка £ = £*, где Г(£) обращается в нуль

Г(£*) = 0, 0 <£* <£*. (45)

Если £* < £ *, то в интервале £* < £ < £ * функция ф (£), определяемая формулой

Ф (£ ) = Г' (£* )(£ — £*),

(46)

будет не больше чем f (£).

На рис. 5-8 приведены аналогичные результаты при Tg = — 4С, Tn = 4С, со = 0.035, vo = 1, T* = — coa — 0.04, a = 1.76273 • 10-7. В этом примере мерзлая область находится в (0.000003789082490396234, 0.000003789082490396235), а граница, разделяющая двухфазную и талую зоны, соответствует в = 0.000328598. Введем еще одну функцию

= У (47)

Поскольку в интервале £* < £ < £* будет ф(£) < f (£), подынтегральная функция р(£) не больше функции р2(£).

в

Если £* < в то интеграл J p(£)d£ можно представить в виде суммы

С»

в с* в

j p(£)d£ = j P(£)d£ + j P(£)d£. (48)

С* С» С*

в С*

Интеграл f p(£)d£ конечный, а интеграл f p(£)d£ равен —ж, так как

С* С*

СГ е-в2/4аг С £

С* С*

(признак сравнения с учетом знака). Если £* > в, то

в в

¡p(£) d£<Jp 2(£)d£ = —ж. (49)

С* С*

Следовательно,

в

J P(£)d£ = —ж. (50)

С*

Таким образом, при выполнении условий (44), (45) или если (более мягкое условие) в интервале (0, в) существует единственная точка, где функция f(£) равна нулю, предел водонасыщенности справа в точке £ = £* равен нулю.

Приведем два примера. Расчеты проведены на сетке с количеством узлов n = 4000 при значениях параметров:

Cs = 920 Дж/кг K; ps = 2000кг/м3; С =2 • 103 Дж/кг K; рг = 900кг/м3; Cw = 4.19 • 103 Дж/кг K; pw = 1000кг/м3; As = 2Вт/м K; Аг = 2,23 Вт/м K; Aw = 0,58 Вт/м K; qw = 3.34 • 105 Дж/кг; a = 66,7; m = 0.1; D = 1,45 • 10-9м2/с; vo = 0.001, 0.035; vo = 1. На всех рисунках линии 1 соответствуют численному расчету, 2 — модели с постоянными коэффициентами в уравнении теплопроводности.

£

Рис. 2. Распределение водонасыщенности

Рис. 5. График функции /(£)

Рис. 6. Распределение водонасыщенности

Рис. 7. Распределение температуры

0.000 0.001 0.002 0.003 0.004 0.005

£

Рис. 8. Распределение температуры

На рис. 2-4 приведены результаты расчетов, проведенных при граничных и начальных условиях Тд = —2С, Тп = 4С, с0 = 0.001, Т* = — с0а — 0.03 и коэффициенте а = 5.06266 • 10-8. По результатам численного расчета точка, где водона-сыщенность равна нулю, лежит между точками = 00.000009197372406967679 и £ = 0.000009197372406967680, а граница раздела талой и двухфазной областей находится в точке в = 0.000385802. На рис. 1 приведены графики функции /(£), где для численного метода ее производная заменена разностным аналогом.

На рис. 2 дано сравнение распределения водонасыщенности. На рис.3, 4 показаны распределения температуры.

Заключение

Результаты вычислительных расчетов показали существование в протяженной области фазовых переходов точки, где выражение {Т-щ + обращается в нуль как для исходной задачи (задача (1)-(10)), так и для задачи с точным решением температуры в полубесконечной области. Также следует отметить хорошую согласованность полученных распределений температуры и водонасыщенности в большей части двухфазной области при высоких концентрациях воды, насыщающей пористую среду. Разница вблизи границы мерзлой области обусловлена большим градиентом водонасыщенности.

ЛИТЕРАТУРА

1. Ентов В. М., Максимов А. М., Цыпкин Г. Г. Образование двухфазной зоны при промерзании пористой среды. Препринт № 269. М., 1986.

2. Васильев В. И., Максимов А. М., Петров Е. Е., Цыпкин Г. Г. Тепломассоперенос в промерзающих и протаивающих грунтах. М.: Наука, 1996.

3. Васильев В. И., Максимов А. М., Цыпкин Г. Г. Математическая модель замерзания-таяния засоленного мерзлого грунта // Прикл. механика техн. физика. 1995. Т. 36, № 5. С. 57-66.

4. Васильев В. И., Попов В. В. Численное решение задачи промерзания грунта // Мат. моделирование. 2008. Т. 20, № 7. С. 119-128.

5. Попов В. В. Двумерная математическая модель промерзания засоленного влажного грунта // Вестн. Северо-Вост. федерал. ун-та. 2014. Т.11, № 5. С 9—23.

6. Попов В. В. Автомодельное решение задачи о промерзании пористой среды, насыщенной раствором соли // Наука и образование. 1998. № 1. С. 92—95.

Статья поступила 28 марта 2017 г. Попов Василий Васильевич

Северо-Восточный федеральный университет имени М. К. Аммосова,

ул. Белинского, 58, Якутск 677000;

Институт проблем нефти и газа СО РАН,

ул. Октябрьская, 1, Якутск 677000

1ш1. рт. р¥¥1@шаИ . ги

Математические заметки СВФУ Апрель—июнь, 2017. Том 24, № 2

UDC 519.63

MATHEMATICAL MODEL OF SOIL FREEZING V. V. Popov

Abstract. We compare two mathematical models of the freezing process of a moist soil saturated with an aqueous solute of salt. For the mathematical model with constant coefficients in the heat equation in a two-phase domain, the possibility of a continuous transition of water saturation to zero at the interface between two-phase and frozen domains is shown.

DOI 10.25587/SVFU.2017.2.9248

Keywords: phase transition, thermodynamic equilibrium, self-similar solution, diffusion, thermal diffusivity, water saturation, moisture.

REFERENCES

1. Entov V. M., Maximov A. M., and Tsypkin G. G., "Formation of two-phase domain during the freezing process in a porous medium," Preprint No. 269, Moscow (1986).

2. Vasil'ev V. I., Maximov A. M., Petrov E. E., and Tsypkin G. G., Teplomassoperenos v Promerzauschih i Protaivaushih Gruntah [in Russian], Nauka, Moscow (1996).

3. Vasil'ev V. I., Maximov A. M., Petrov E. E., and Tsypkin G. G., "A mathematical model for the freezing—thawing of a saline frozen soil," J. Appl. Mech. Tech. Phys., 36, No. 5, 57—66 (1995).

4. Vasil'ev V. I. and Popov V. V., "Numerical solution of the soil freezing problem," Mat. Modelir., 20, No. 7, 119-128 (2008).

5. Popov V. V., "Two-dimensional mathematical model of the freezing process of a moist soil saturated with a solute of salt," Vestn. Severo-Vostoch. Federal. Univ., 11, No. 5 19-23 (2014).

6. Popov V. V., "Automodel solution to the problem of the freezing process of a potous medium saturated with an aqueous solute of salt," Nauka Obrazov., No. 1, 92-95 (1998).

Submitted March 28, 2017 Vasily V. Popov

M. K. Ammosov North-Eastern Federal University, Institute of Oil and Gas problem SB RAS, 58, Belinsky St., Yakutsk 677000, Russia imi.pm.pvvl@mail.ru

The work was supported by the grant from the Russian Foundation for Basic Research (project code 17.01.00732).

© 2017 V. V. Popov