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

К числу классических задач математической физики относятся задачи о построении пространственно одномерных асимптотических решений уравнения диффузии в полупространстве и пластине при гармонических по времени граничных условиях Дирихле, Неймана и смешанного типа. Решения такого рода задач имеют вид бегущих затухающих гармонических волн, распространение которых сопровождается дисперсией, а глубина проникновения быстро уменьшается с ростом частоты. В теории электромагнетизма в качестве примера здесь можно указать на работы [1, с. 457– 460; 2, с. 449–457], в которых исследуется явление скин-эффекта, в теории теплопроводности – на работы [3, с. 542–549; 4, с. 99–103; 5, с. 656–659; 6, с. 238–247; 7, с. 85–87, с. 109–112; 8, с. 131–140;

не учитывают наличия в почве влаги и связанных с ней процессов испарения и конденсации. В работах авторов [11; 12] предпринята попытка этот недостаток устранить. Там разработан общий алгоритм для решения задачи Фурье о температурных колебаниях в однородном полупространстве с учетом содержащейся в материале влаги. Однако полученное авторами общее решение было проанализировано ими лишь для случая простейшей модели тепломассопереноса, в которой не учитываются явления термодиффузии и внутреннего испарения. Другим недостатком указанных работ является то, что они не содержит численных расчетов, позволяющих сравнить предсказанные теорией результаты с опытными данными. В настоящей статье мы проанализируем модель общего вида, проведем и сопоставим с известными экспериментальными данными конкретные расчеты при характерных для практики значениях переменных.

1.    Постановка задачи для уравнений распространения тепла и влаги

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

Известно, что совместные начально-краевые задачи для полей Т и U , в отличие от задач для одного только поля Т , большей частью формулируются как нелинейные и что одной из основных причин нелинейности является формула для интенсивности массообмена J . Здесь мы примем для этой величины выражение в форме закона испарения Дальтона :

J ( т ) = a _m [ P ( T (0, т ) ) -ф P ( т_в ( т ) ) ] ;

,    17,3 T

P ( T ) = 6,03 ■ 10 ^{- 3} exp-------.

т + T 1

Здесь am - коэффициент массообмена по перепаду давления пара; Р(Т) – функция, моделирующая известную из опыта зависимость относительного давления насыщенного водяного пара от его температуры Т при общем нормальном давлении; Т 1 = 238 °С - постоянная. Формула справедлива лишь в ситуации, когда водяной пар вблизи поверхности x = 0 является насыщенным. В настоящей статье мы будем изучать ситуацию, когда температура поверхности T(0, т) и температура воздуха TB( т) совершают малые колебания вблизи одной и той же температуры Т0. Для малых отклонений указанных температур от Т0 зависимость J(т) можно линеаризовать, разложив функцию Р(Т) в ряд Тейлора в окрестности точки Т0. Выполнив линеаризацию и приняв ф = 1, т. е. считая, что водяной пар находится в насыщенном состоянии не только вблизи поверхности материала, но и за пределами пограничного слоя, получим приближенную формулу для закона Дальтона:

J ( т ) = a m [ T (0, т ) - T в ( т ) ] ;

^a m =a m P ‘ ( ^T 0 ) .

Здесь a _m - коэффициент массообмена по перепаду температуры.

Отметим, что закон Дальтона для интенсивности массообмена, записанный в таком виде, совпадает по форме с законом Ньютона для интенсивности теплообмена:

Q(т) = aw [T(0, т) - Tb(т)], который мы и будем использовать при формулировке граничных условий. В последней формуле a w - коэффициент теплообмена поверхности образца с воздушной средой.

Линеаризация закона Дальтона позволяет переписать в линейном приближении всю систему уравнений и краевых условий, описывающих процессы распространения тепла и влаги в нашей задаче. Чтобы записать эту систему, введем следующие обозначения: c, р, X, у, am, 5 - теплофизические характеристики материала, а именно удельная теплоемкость, плотность в сухом состоянии, коэффициент теплопроводности, критерий испарения, коэффициент диффузии влаги, относительный коэффициент термодиффузии влаги; a w = X /(c р) - коэффициент диффузии тепла (коэффициент температуропроводности); r – удельная теплота парообразования воды. Необходимая нам система будет иметь следующий вид [13; 14]:

д T = д ² T r ydU_ дт ^w 3 _x ² c дт ;

д U

= a дт ^m

d ² U     . д ² T

--+ a 5----- ; д x ^{2 m} д x ²

am [T(0, t)- Tb(t)] = amp — (0, t) + 5—(0, t) ; (3) ^ dx dx aw [T(0, t) - Tb(t)] + (4)

+ r (1 - Y ) a _m [ T (0, t ) - T _b( t ) ] = \-~ (0, t ). d x

Здесь (1) и (2) – уравнения для потоков тепла и влаги в области 0 <  x < » , занятой материалом, а (3) и (4) – граничные условия для этих потоков на поверхности x = 0. Условие (3) выражает равенство двух потоков влаги, отводимого от поверхности материала в воздух через пограничный слой и подводимого к этой поверхности изнутри материала, а условие (4) выражает соотношение между тепловыми потоками на поверхности x = 0, а именно между потоком тепла, отводимым по закону Ньютона с поверхности в воздух, потоком тепла, необходимым для испарения с поверхности подступающего к ней из материала потока жидкости, и потоком тепла, подходящим к поверхности изнутри материала.

2.    Гармонические условия на границе и постановка задачи об асимптотике

Пусть при т <  0 температура материала и его влагосодержание имели постоянные по всему полупространству значения Т ₀ и U ₀, температура воздуха Т _в равнялась температуре материала Т ₀, а влажность воздуха ф была равна 1. Очевидно, что в таком состоянии система может находиться неограниченно долго, поскольку все уравнения нашей задачи оказываются удовлетворенными, причем для интенсивностей тепломассообмена мы будем иметь Q = 0 и J = 0.

Примем далее, что с момента т = 0 температура воздуха Т _в начинает совершать вблизи Т ₀ малые гармонические колебания по закону

T _b( t ) = T 0 +^А T _Bsin( roT + y _B), (5) где Л T _B, го , ф _в - заданные величины, а влажность воздуха остается неизменной и сохраняет принятое выше значение ф = 1. Несложно убедиться, что тогда, рассмотрев систему (1)–(5), мы сможем поставить вопрос о нахождении решения этой системы вида

T (x, t) = Tq + T (x )sin (roT + vt( x)), (6) U (x, t) = U о + U (x )sin (гот +фц( x)), в котором температура материала Т и его влагосо-держание U будут совершать при каждом x малые гармонические колебания вблизи своих первона- чальных значений Т0 и U0, а интенсивности тепло- и массообмена Q и J будут совершать такие же колебания вблизи своих первоначальных значений Q = 0 и J = 0. В этом квазистационарном состоянии системы, которое приходит на смену ее исходному стационарному состоянию, нагревание поверхности материала и связанное с ним испарение влаги с поверхности будут периодически замещаться охлаждением поверхности и ее пропиткой, а средние за период колебаний температура и влагосодержание при любом x будут, как отмечалось выше, оставаться неизменными.

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

T ( x ) ^ 0 и U ( x ) ^ 0 при x ^» . (7)

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

Сформулированная задача о нахождении решений вида (6) относится к числу задач без начальных данных [6]; искомые в этой задаче функции дают асимптотику полей Т и U при т ^ » .

3.    Постановка задачи для комплексов

Построение асимптотического решения (6) проведем методом комплексных амплитуд . Вместо функций T ( x , t ), U ( x , t ) введем новые искомые функции T ( x , t ), U ( x , t ):

T (x, t) = T (x, t) - Tq = T (x )sin (гот + vt( x)), (8) U (x, t) = U (x, t) - Uq = U (x )sin (гот +фц( x)), которым, в свою очередь, сопоставим их комплексы T(x) и U(x):

T ( x , t ) о T ( x ) = T ( x )exp ( i v _t( x ) ) , (9) U ( x , t ) о U ( x ) = U ( x )exp ( i ф _ц( x ) ) .

Выражая в системе (1)–(4) функции Т и U через T и U и пользуясь правилами работы с комплексами, вместо уравнений для T и U получим уравнения для комплексов этих функций T и U :

i го T ( x ) = a.	d 2 T ( x ) w dx 2	+ i го —UJ ( x ); c		(10)
•	d 2 U ( x )		s d 2 T ( x ).	(11)
	m dx 2	+ a	m Л 2 ; dx
*m [ T ( 0 ) —	A T B ] = a	m P	du (0) + 5 dT (0) dx      dx	;             (12)
a w [ T (0) —	A T B ] = >	dT dx	(0).	(13)

Здесь обозначено a w = a w + r (1 - Y )a m, A TB = A TB exp( i vB)

– эффективный коэффициент теплообмена и заданное комплексное число соответственно.

В задаче для комплексов условие на бесконечности (7) будет выглядеть так:

|T ( x )p 0 и |lj ( x )|^ 0 при x ^» .             (14)

Решив задачу (10)-(14) для комплексов T и U J , мы затем с помощью (8) и (9) найдем отвечающие этим комплексам асимптотические решения (6).

• 4.    Решение задачи для комплексов

Пользуясь методом Эйлера [15], построим общее решение системы дифференциальных уравнений (10), (11), удовлетворяющее условиям на бесконечности (14), а затем из системы краевых условий (12), (13) найдем входящие в это общее решение произвольные постоянные, чем и завершится процесс нахождения решения для комплексов. В подробном виде эта процедура изложена авторами в статье [11], здесь же мы приведем лишь ее результаты. Для записи решения в указанной работе были введены следующие обозначения:

v = (1 + aw/a m +5 ry/ c) /2;

T1,2 = r-1---------, 1            ;(16)

c 1 ^-v +     ^- a w/ a m

^1,2 = -в1,2(1 + i),

₽ i,2=7^^ ^{2 a} w⁾ ] ^v±7^v 2 - a w / a m .

Величина v является безразмерной, Т12 имеет размерность °С, а размерностью в12 является 1/м. Первые две величины – постоянные, определяемые свойствами материала, а третья величина кроме свойств материала зависит еще от частоты го. Отметим, что во всех формулах под 4z мы понимаем, как обычно, главное значение корня квадратного из комплексного числа z (Re zjz > 0). Кроме указанных пяти величин, нахождение ко- торых может быть произведено непосредственно по исходным данным, следует еще, согласно [11], определить две безразмерные и зависящие от частоты го величины С12, которые находятся как решение следующей системы линейных алгебраических уравнений:

[ T 1 - a m PH (1 + 5 T 1 )/ « m ] ^C 1 +                     (18)

+ [ T 2 - a m P^ 2 (1 + 8 T 2 )/«m ] C 2 =A T b ;

T 1 (1 -X_H/a w ) C 1 +

+ T 2 (1 -X^ Л w ) C 2 =A T b .

После нахождения С ₁₂ решение задачи (10)–(14) для комплексов T и U J в матричном виде может быть записано как

T ( x ) U ( x )

= C ₁

T 1 exp ( ц 1 x ) exp ( ц 1 x )

)

+ C 2

^T2 ^exP ( ^ 2 ^x exp ( ^ 2 x )

Приняв к сведению это решение из работы [11], мы должны теперь в качестве следующего шага перейти от комплексов ( изображений ) к вещественным функциям ( оригиналам ) и проанализировать зависимость этих функций от координаты, времени, частоты и характеристик материала. Рассмотрим сначала простую ситуацию, когда у = 0 и 5 = 0, что соответствует математической модели тепломассопереноса, в которой пренебрегают явлениями внутреннего парообразования и термодиффузии, вследствие чего функции Т и U оказываются связанными только через граничные уравнения (условия применимости такой модели к задачам практики и анализ получаемых при ее использовании решений можно найти в работе [16]). Решение для комплексов (15)–(19) в случае такой модели становится существенно более простым; переход от этого решения к оригиналам и анализ получаемых при этом функций проведены авторами в статье [12]. Здесь же мы впервые осуществим полный анализ решения для комплексов (15)–(19), которому будет отвечать математическая модель тепломассопереноса общего вида, свободная от указанных выше ограничений. Наше исследование будет выглядеть следующим образом: мы зафиксируем параметры материала и частоту, и, рассчитав при этих условиях коэффициенты в формуле (19), проанализируем зависимости полей Т и U от координаты и времени.

5.    Расчет коэффициентов в формуле для комплексов

Материалом полупространства будем считать глину , ее теплофизические характеристики имеют следующие значения [13]: Х = 0,93 Вт/(м ■ °С); с = 1,9 ■ 10³ Дж/(кг ■ °С); р = 1,5 ■ 10³ кг/м³; y = 0,1; 5 = 1,5 ■ 10^-3 1/°С); a _m = 2,6 ■ 10^-8 м²/с; a _w = 0,32 х х 10^-6 м²/с.

Эти данные позволяют вычислить коэффициенты (15) и (16): v = 6,7432; Т 1 = -10,347 °С; Т ₂ = 12,790 ■ 10 3 °С.

Рассчитаем теперь круговую частоту колебаний го = 2 п / Т , выбрав в качестве периода колебаний Т один год , а затем найдем коэффициенты (17): го = 1,9912 ■ 10^-7 с^-1; 3 1 = 1,9725 1/м; 3 2 = 0,55540 1/м.

При расчете коэффициентов тепло- и массо-обмена будем опираться на фундаментальное руководство по этому вопросу [17], на статью авторов [14], где для этих коэффициентов получены приближенные формулы, а также на имеющиеся в [13] опытные данные. На основании этих работ в качестве характерных значений коэффициентов тепло- и массообмена в опытах, где влажная глина будет обдуваться потоком воздуха, примем a _w = 5 Вт/(м² ■ °С); а _т = 5 ■ 10^-3 кг/(м²<).

В качестве температуры, вблизи которой происходят колебания температуры воздуха, возьмем Т 0 = 20 °С. Тогда расчет коэффициентов a _w и a _m приводит к таким значениям: a _w = 19,5 Вт/(м² ■ °С); a _m = 7,13 ■ 10^-6 кг/(м²^°С).

Наконец, в качестве характеристик колебаний температуры воздуха в формуле (5) примем Л Т _в = 5 °С, ^ = 0.

Теперь мы можем обратиться к системе (18). После подсчета коэффициентов она приобретает такой вид:

(0,27501 + i 10,622) С 1 + (12851 + i 61,322) С 2 = 5;

( - 11,320 - i 0,97336) С 1 + (13129 + i 338,78) С 2 = 5.

Здесь скобки и правые части уравнений имеют размерности °С. Решая эту систему, найдем, что С 1 = | С 1| exp( i arg С 1 );                                   (20)

|С 1| = 9,306510 ^{- 3}; arg С 1 =- 0,024990;

^С 2 = | ^С 2I ^exp⁽ i ^arg ^С 2 );

I С ₂| = 0,38875 ■ 10 ^{- 3}; arg С ₂ =- 0,024554.

Таким образом, все коэффициенты в решении для комплексов (19) нам стали известными.

6.    Анализ поля влагосодержания

Обратимся сначала к полю U . Согласно формулам (19), (20) и (17),

U J ( x ) = С 1 exp( p 1 x ) + С 2 exp( ^ 2 x ) = = С 11 exp( i arg С 1 )exp [—3 1 (1 + i ) x ] + + | С 21 exp( i arg С 2 )exp [-3 ₂ (1 + i ) x ] .

Выделяя в каждом из двух слагаемых модуль и аргумент, найдем соответствующие этим двум комплексам оригиналы; добавив к получившейся сумме постоянную U ₀, будем иметь оригинал для поля влагосодержания в следующем виде: U ( x , t ) - U 0 =

= С I exp ( -3 1 x ) sin ( гот + arg С 1 - 3 1 x ) +            (21)

+ | С 21 exp ( -З 2 x ) sin ( ror + arg С 2 -З 2 x ) •

Каждое из двух слагаемых в правой части этой формулы представляет собой бегущую затухающую гармоническую волну . Рассмотрим первое слагаемое. Согласно терминологии, принятой в теории волн, величина 3 1 есть коэффициент затухания этой первой волны влагосодержания, а обратная величина Л 1 = 1 / 3 1 - ее глубина проникновения . Другими характеристиками волны являются ее фазовая скорость V 1 = го / 3 1 и длина волны Л 1 = 2 п / 3 1 . Для гармонической волны общего вида в формулах для фазовой скорости и длины волны вместо коэффициента затухания 3 1 должен стоять коэффициент фазы 3 1 , который является коэффициентом при координате x под знаком синуса, но в нашем случае эти два коэффициента совпадают.

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

Аналогичный комментарий можно дать и для второго слагаемого в формуле (21).

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

A1 (x) = |C1| exp(-P1 x); ф1(x) = arg C1 -01 x;(22)

A ₂ ( x ) = C ₂1 exp( -p ₂ x ); ф ₂ ( x ) = arg C ₂ - P ₂ x .

Тогда формулу (21) для поля влагосодержания можно будет переписать так:

U (x, t) - U о = A( x) sin [мт + ф( x)];(23)

A( x) = A^ + 2 Ai A 2 cos(ф1 - Ф2) + A 2 ,(24)

A., sinф + A, sin ф9 sin ф( x) = —-----1----2

A

, , A. cosф + A9 cosф9 cos ф( x) = —-----1----2

A

Здесь величины А 1 , А 2 , ф 1 , ф 2 являются функциями координаты x и вычисляются по формулам (22).

Обратим теперь внимание на то, что второе слагаемое в формуле (21), согласно (20), имеет незначительную амплитуду по сравнению с первым. Воспользовавшись этим, получим, что поле U приближенно можно считать «чистой» затухающей гармонической волной

U ( x , т ) - U 0 =                                           (25)

= ЛUп exp(-x/A1 )sin[m(t-Ati)-P1 x], где AUп = |C1| = 9,31-10-3 - амплитуда колебаний влагосодержания на поверхности х = 0; Ai = = у Pi = 0,507 м - глубина проникновения волны влагосодержания; Ati =-(argCi)/го = 1,45 сут -время запаздывания колебаний на поверхности x = = 0 относительно колебаний температуры воздуха.

Общее время запаздывания колебаний на произвольной глубине x будет определяться формулой AT 1 ( x ) = AT 1 + P 1 x / ro .

Здесь с ростом глубины второе слагаемое в правой части быстро становится много больше первого; например, уже на глубине x = 1 м будем иметь

P 1 x/ ro = 3,82 мес » AT 1 .

7.    Анализ поля температуры

Аналогичным образом может быть проведен анализ и поля температуры. Согласно (19),

T ( x ) = T 1 C 1 exp( ^ 1 x ) + T 2 C 2 exp( ^ 2 x ).

Поскольку множители Т1 и Т2 – вещественные, то оригинал для поля температуры будет иметь тот же вид (21), но только теперь слева будет стоять выражение T (x, t) - To, а справа к первому и ко второму слагаемым добавятся множители Т1 и Т2 соответственно. Анализ получившегося выраже- ния не отличается от того анализа, который мы провели для поля влагосодержания. Завершая его, мы в конце должны будем заметить, что

T 1 | C 1| = - 0,0963 ° C; T ₂ C ₂| = 4,97 ° C.

Значит, для приближенного описания температурного поля мы можем ограничиться лишь вторым слагаемым в полученной формуле, и тогда эта приближенная формула, представляющая «чистую» затухающую гармоническую волну, будет выглядеть как

T ( x , t ) - T 0 =

= A T _n exp( - x/ A 2 )sin [ro ( T-AT 2 ) -P 2 x J .

Здесь A T _n = T 2 | c ₂| = 4,97 ° C - амплитуда колебаний температуры на поверхности x = 0, которая мало отличается от амплитуды колебаний температуры воздуха A Т _в = 5 ° C; A 2 = 1/ Р 2 = 1,80 м -глубина проникновения волны температуры; AT 2 = =- (arg C 2 )/ го = 1,43 сут - время запаздывания колебаний на поверхности x = 0 относительно колебаний температуры воздуха.

Общее время запаздывания колебаний на произвольной глубине x будет определяться формулой AT 2 ( x ) = AT 2 + Р 2 x / ГО .

Здесь с ростом глубины второе слагаемое в правой части быстро становится много больше первого; например, уже на глубине x = 4 м будем иметь P 2 x /ro = 4,30 мес » AT 2 .

Поскольку коэффициент P 2 является функцией частоты го , то распространение волн температуры, как и волн влагосодержания, сопровождается дисперсией. Явление дисперсии тепловых волн в полупространстве, при условии, что материал не содержит влаги, исследовано А.В. Лыковым методом преобразования Лапласа в работе [6, с. 306– 308]. Согласно полученным там результатам, фазовая скорость волны оказывается пропорциональной корню квадратному из произведения коэффициента диффузии тепла а _w на частоту. Можно убедиться, что в нашем случае, когда наличие влаги при расчете было учтено, это утверждение останется в силе, если только мы ограничимся приближением «чистой» затухающей гармонической волны.

8.    Обсуждение результатов

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

В простейшей модели (будем говорить, что это модель А ) материал влаги не содержит, исследуется только поле температуры. Температурные волны для такого случая описываются известной в литературе формулой Фурье [6], которая в обозначениях, принятых в настоящей статье, выглядит так: T ( x , т ) = T 0 + (26) ^{+ A} T n ^exP⁽ -^P w x ^)sin [^{ит +} ^ п — ^P w x ] •

В этой формуле

P w = V^ю/ (2 a w ) (27) - коэффициент затухания волны, а Л T _n и ^ _п есть амплитуда и начальная фаза колебаний температуры на поверхности материала (в отличие от нашего рассмотрения, в задаче Фурье изменяется по гармоническому закону и считается заданной не температура воздуха, а температура поверхности x = 0).

В двух других следующих по сложности моделях наличие в материале влаги учитывается в рамках теории тепломассопереноса А.В. Лыкова, но в одной из них полагают у = 0 и 5 = 0, что соответствует пренебрежению явлениями внутреннего парообразования и термодиффузии ( модель В ), а в другой модели этими явлениями не пренебрегают и такого ограничения не вводят ( модель С ). Модель В исследовалась авторами в статьях [11; 12], ее формулы получаются как частный случай более общих формул настоящей статьи, в которой исследовалась модель С . В моделях А и В волны являются затухающими гармоническими, а в модели С они могут считаться такими лишь в рамках принятых в предыдущих двух пунктах приближений. Можно говорить, что полученные в рамках моделей В и С результаты, содержащие описание взаимосвязанных процессов распространения полей температуры и влагосодержания, представляют собой обобщение теории Фурье для температурных волн в полупространстве.

Ниже для каждой из моделей мы приводим итоговые формулы для характеристик полей. В случае модели С, принимая сформулированные выше приближения, мы считаем волны затухающими гармоническими. Нижний индекс указывает, о каком поле, Т или U, идет речь в данной формуле, а верхний индекс определяет принадлежность формулы к моделям А, В или С. Формулы статьи [11] приведены в соответствие с обозначениями, принятыми в настоящей статье. Мы обозначили также Lu = a ml a w - критерий Лыкова.

• а)    Коэффициент затухания в :

P T = P B = P w > P B = P w v- \ 2 - и - ' ;

P U = P w LL-^ ; P U = P w J v+Vv ² - Lu ^{- 1} •

• б)    Время запаздывания колебаний Лт на глубине x относительно колебаний на поверхности:

Лт A =Лт B ^{P w} x ; Лт B = ^P T x ; ю        ю

Лт B ^P- x 'LF ¹ ; Лт B =^P U x •

U ю           U ю

• в)    Время запаздывания колебаний Л t на поверхности относительно колебаний температуры воздуха:

Л t B = 1 arctg (—1 ; Л t B =-^ ;

T ю   Ba w ⁺^P w J    T ю ’

Л t B =Л t B ; Л t B =- arg B •

ю

• г)    Амплитуды колебаний на поверхности Л Т и Л U :

  Л T B

  
  

  ;

  V ( a w +^P w )² + ( ^P w )²

  
  

л tb = r¹----^=— B 2I ;

c 1 -v + yv ² - Lu ¹

Л UB =--- ^a m ^   Л TB ; Л UB = B l

6 w ^p w m

Ниже в таблице приведены результаты расчетов по этим формулам для разных моделей и в условиях, принятых в настоящей статье: материалом полупространства мы считаем влажную глину, в качестве периода колебаний выбираем один год, характеристиками воздуха являются T о = 20 °С и Л Т _в = 5 °С, расчет времени запаздывания Лт производится на глубине x = 4 м.

Выводы из представленных данных.

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

• 2.    Учет влаги с помощью модели В , согласно пунктам а) и б), не изменяет формул для характеристик температурного поля, полученных в рамках модели А , т. е. формул для глубины проникновения β и времени запаздывания колебаний ∆τ на произвольной глубине относительно колебаний на поверхности. Следовательно, даже при наличии в почве влаги теория Фурье для тепловых волн должна приводить к удовлетворительному согласию с экспериментом.

• 3.    Приведенные в таблице значения параметров β и ∆τ для поля температуры находятся в хорошем согласии с данными наблюдений на метеорологических станциях [6].

Таблица. Результаты расчетов в рамках моделей А , В , С

Table. Results of calculations within the framework of models A , B , C

Модель	β , 1/м		∆τ , мес		∆ t , сут		Амплитуды
	Т	U	Т	U	Т	U	∆ Ò	∆ U
А	0,558	–	4,30	–	–	–	Задается	–
В	0,558	1,96	4,30	15,2	1,39	1,39	4,88 °С	0,0112
С	0,555	1,97	4,31	15,3	1,43	1,45	4,97 °С	0,00931

Кроме геокриологии математическое моделирование процессов распространения тепла и влаги в полупространстве является актуальным также и в метеорологии и климатологии. Опытное отслеживание годовых колебаний температуры в поверхностном слое Земли позволяет сделать вывод, что температурные колебания суши затухают на глубине порядка 15–20 метров, в то время как в океанах глубина проникновеня тепловых волн оказывается намного больше (это различие возникает благодаря эффекту перемешивания [18]). За счет этого в летнее время океаны накапливают тепла значительно больше, чем суша, а это влияет на формирование климата, поскольку зимой океаны начинают согревать проходящие над ними массы воздуха.

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

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

В настоящей статье, расчетной основой которой является метод комплексных амплитуд, авторы продолжили исследование возможностей аналитических методов для решения совместных начально-краевых задач для уравнений распространения тепла и влаги. В предыдущих статьях, посвященных этому направлению [21; 22], ими был разработан расчетный алгоритм на основе метода разделения переменных. К исходной задаче для уравнений тепломассопереноса, содержащих не одну, а две искомые функции, этот классический алгоритм применен быть не может; данное затруднение преодолевается авторами с помощью идеи расщепления процесса по физическим факторам.

Заключение

Разработана система уравнений и краевых условий, моделирующих процессы распространения тепла и влаги в однородном, содержащем влагу полупространстве, граница которого обдувается воздухом с изменяющейся по гармоническому закону температурой. Использованы уравнения теории тепломассопереноса А.В. Лыкова, в которых учитываются явления термодиффузии и внутреннего парообразования, а краевые условия для потоков тепла и влаги формулируются на основе закона теплообмена Ньютона и закона испарения Дальтона соответственно. Получены асимптотические по времени распределения температуры и влагосодержания как функции координаты, времени, частоты и характеристик материала. Эти распределения представляют собой суперпозиции бегущих гармонических волн, имеющих различные коэффициенты затухания и фазовые скорости. Распространение таких волн сопровождается дисперсией. В рамках математических моделей тепломассопереноса различного уровня сложности проанализирована зависимость от условий опыта глубины проникновения волны, амплитуд колебаний температуры и влагосодержания на поверхности материала и времени запаздывания колебаний на произвольной глубине относительно колебаний на поверхности. Результаты расчетов в условиях, приблизительно соответствующих ус- ловиям наблюдений на метеорологических станциях, хорошо соотносятся с опытными данными. Разработанный в статье расчетный алгоритм может считаться обобщением теории Фурье для температурных колебаний в полупространстве при отсутствии влаги и при граничных условиях теплообмена первого рода. Материалы работы могут быть использованы в геокриологии в качестве теоретического инструмента при моделировании сезонных колебаний теплофизического состояния мерзлых пород и грунтов. С их помощью моделирование состояния грунтов можно будет производить с учетом содержащейся в них влаги, что ранее можно было делать лишь в рамках простейших моделей.