Solving geocryology problems based on generalized Fourier theory for temperature waves in half-space

Full Text

Введение

К числу классических задач математической физики относятся задачи о построении пространственно одномерных асимптотических решений уравнения диффузии в полупространстве и пластине при гармонических по времени граничных условиях Дирихле, Неймана и смешанного типа. Решения такого рода задач имеют вид бегущих затухающих гармонических волн, распространение которых сопровождается дисперсией, а глубина проникновения быстро уменьшается с ростом частоты. В теории электромагнетизма в качестве примера здесь можно указать на работы [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; 9, с. 298–313], объектом исследования в которых являются волны температуры. Одним из наиболее известных результатов, касающихся температурных волн, являются полученные еще Фурье и Пуассоном формулы, моделирующие колебания температуры в поверхностном слое земной коры, вызванные суточными и годовыми колебаниями температуры окружающей среды. В настоящее время формулы Фурье и являющиеся их следствием законы Фурье [6, с. 244–245] нашли широкое применение в мерзлотоведении (геокриологии), где они используются для решения различных проблем, связанных с организацией хозяйственной деятельности в области распространения мерзлых пород [10]. Существенным недостатком формул Фурье, который ограничивает их применение для решения указанных проблем, является то, что они не учитывают наличия в почве влаги и связанных с ней процессов испарения и конденсации. В работах авторов [11; 12] предпринята попытка этот недостаток устранить. Там разработан общий алгоритм для решения задачи Фурье о температурных колебаниях в однородном полупространстве с учетом содержащейся в материале влаги. Однако полученное авторами общее решение было проанализировано ими лишь для случая простейшей модели тепломассопереноса, в которой не учитываются явления термодиффузии и внутреннего испарения. Другим недостатком указанных работ является то, что они не содержит численных расчетов, позволяющих сравнить предсказанные теорией результаты с опытными данными. В настоящей статье мы проанализируем модель общего вида, проведем и сопоставим с известными экспериментальными данными конкретные расчеты при характерных для практики значениях переменных.

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

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

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

$\begin{array}{l}J\left(\tau \right)={\alpha }_{\text{m}}\left[P\left(T\left(\mathrm{0,}\tau \right)\right)-\phi P\left(T_{\text{в}}\left(\tau \right)\right)\right];\\ P\left(T\right)=\mathrm{6,03}\cdot {10}^{-3}\exp \frac{\mathrm{17,3}T}{T+T_1}.\end{array}$

Здесь ${\alpha }_{\text{m}}$ – коэффициент массообмена по перепаду давления пара; Р(Т) – функция, моделирующая известную из опыта зависимость относительного давления насыщенного водяного пара от его температуры Т при общем нормальном давлении; ${Т}_1=238$°С – постоянная. Формула справедлива лишь в ситуации, когда водяной пар вблизи поверхности x = 0 является насыщенным. В настоящей статье мы будем изучать ситуацию, когда температура поверхности $T\left(\mathrm{0,}\tau \right)$ и температура воздуха $T_{\text{в}}\left(\tau \right)$ совершают малые колебания вблизи одной и той же температуры Т₀. Для малых отклонений указанных температур от Т₀ зависимость $J\left(\tau \right)$ можно линеаризовать, разложив функцию Р(Т) в ряд Тейлора в окрестности точки Т₀. Выполнив линеаризацию и приняв $\phi =\mathrm{1,}$ т. е. считая, что водяной пар находится в насыщенном состоянии не только вблизи поверхности материала, но и за пределами пограничного слоя, получим приближенную формулу для закона Дальтона:

$\begin{array}{l}J\left(\tau \right)={\tilde{\alpha }}_{\text{m}}\left[T\left(\mathrm{0,}\tau \right)-T_{\text{в}}\left(\tau \right)\right];\\ {\tilde{\alpha }}_{\text{m}}={\alpha }_{\text{m}}{P}^{\text{'}}\left(T_0\right).\end{array}$

Здесь ${\tilde{\alpha }}_{\text{m}}$ – коэффициент массообмена по перепаду температуры.

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

$Q\left(\tau \right)={\alpha }_{\text{w}}\left[T\left(\mathrm{0,}\tau \right)-T_{\text{в}}\left(\tau \right)\right],$

который мы и будем использовать при формулировке граничных условий. В последней формуле ${\alpha }_{\text{w}}$ – коэффициент теплообмена поверхности образца с воздушной средой.

Линеаризация закона Дальтона позволяет переписать в линейном приближении всю систему уравнений и краевых условий, описывающих процессы распространения тепла и влаги в нашей задаче. Чтобы записать эту систему, введем следующие обозначения: c, ρ, λ, γ,  a_m, δ – теплофизические характеристики материала, а именно удельная теплоемкость, плотность в сухом состоянии, коэффициент теплопроводности, критерий испарения, коэффициент диффузии влаги, относительный коэффициент термодиффузии влаги; $a_{\text{w}}=\lambda /\left(c\rho \right)$ – коэффициент диффузии тепла (коэффициент температуропроводности); r – удельная теплота парообразования воды. Необходимая нам система будет иметь следующий вид [13; 14]:

$\frac{\partial T}{\partial \tau }=a_{\text{w}}\frac{{\partial }^{2}T}{\partial x^2}+\frac{r\gamma }{c}\frac{\partial U}{\partial \tau };$ (1)

$\frac{\partial U}{\partial \tau }=a_{\text{m}}\frac{{\partial }^{2}U}{\partial x^2}+a_{\text{m}}\delta \frac{{\partial }^{2}T}{\partial x^2};$ (2)

${\tilde{\alpha }}_{\text{m}}\left[T\left(\mathrm{0,}\tau \right)-T_{\text{в}}\left(\tau \right)\right]=a_{\text{m}}\rho \left[\frac{\partial U}{\partial x}\left(\mathrm{0,}\tau \right)+\delta \frac{\partial T}{\partial x}\left(\mathrm{0,}\tau \right)\right];$ (3)

${\alpha }_{\text{w}}\left[T\left(\mathrm{0,}\tau \right)-T_{\text{в}}\left(\tau \right)\right]+$ (4)

$+r(1-\gamma ){\tilde{\alpha }}_{\text{m}}\left[T\left(\mathrm{0,}\tau \right)-T_{\text{в}}\left(\tau \right)\right]=\lambda \frac{\partial T}{\partial x}\left(\mathrm{0,}\tau \right).$

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

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

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

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

$T_{\text{в}}\left(\tau \right)=T_0+\Delta T_{\text{в}}\sin (\omega \tau +{\psi }_{\text{в}}),$ (5)

где $\Delta T_{\text{в}},$ $\omega ,$ ${\psi }_{\text{в}}$ – заданные величины, а влажность воздуха остается неизменной и сохраняет принятое выше значение $\phi =1.$ Несложно убедиться, что тогда, рассмотрев систему (1)–(5), мы сможем поставить вопрос о нахождении решения этой системы вида

$T(x,\tau )=T_0+T\left(x\right)\sin \left(\omega \tau +{\psi }_{\text{t}}\left(x\right)\right),$ (6)

$U(x,\tau )=U_0+U\left(x\right)\sin \left(\omega \tau +{\psi }_{\text{u}}\left(x\right)\right),$

в котором температура материала Т и его влагосодержание U будут совершать при каждом x малые гармонические колебания вблизи своих первоначальных значений Т₀ и U₀, а интенсивности тепло- и массообмена Q и J будут совершать такие же колебания вблизи своих первоначальных значений Q = 0 и J = 0. В этом квазистационарном состоянии системы, которое приходит на смену ее исходному стационарному состоянию, нагревание поверхности материала и связанное с ним испарение влаги с поверхности будут периодически замещаться охлаждением поверхности и ее пропиткой, а средние за период колебаний температура и влагосодержание при любом x будут, как отмечалось выше, оставаться неизменными.

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

$T\left(x\right)\to 0\text{\hspace{0.33em}\hspace{0.33em}и\hspace{0.33em}\hspace{0.33em}}U\left(x\right)\to 0\text{\hspace{0.33em}\hspace{0.33em}при\hspace{0.33em}\hspace{0.33em}}x\to \infty .$ (7)

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

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

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

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

$\tilde{T}(x,\tau )=T(x,\tau )-T_0=T\left(x\right)\sin \left(\omega \tau +{\psi }_{\text{t}}\left(x\right)\right),$ (8)

$\tilde{U}(x,\tau )=U(x,\tau )-U_0=U\left(x\right)\sin \left(\omega \tau +{\psi }_{\text{u}}\left(x\right)\right),$

которым, в свою очередь, сопоставим их комплексы $\dot{T}\left(x\right)$ и $\dot{U}\left(x\right)\text{:}$

$\tilde{T}(x,\tau )\leftrightarrow \dot{T}\left(x\right)=T\left(x\right)\exp \left(i{\psi }_{\text{t}}\left(x\right)\right),$ (9)

$\tilde{U}(x,\tau )\leftrightarrow \dot{U}\left(x\right)=U\left(x\right)\exp \left(i{\psi }_{\text{u}}\left(x\right)\right).$

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

$i\omega \dot{T}\left(x\right)=a_{\text{w}}\frac{d^2\dot{T}\left(x\right)}{d{x}^2}+i\omega \frac{r\gamma }{c}\dot{U}\left(x\right);$ (10)

$i\omega \dot{U}\left(x\right)=a_{\text{m}}\frac{d^2\dot{U}\left(x\right)}{d{x}^2}+a_{\text{m}}\delta \frac{d^2\dot{T}\left(x\right)}{d{x}^2};$ (11)

${\tilde{\alpha }}_{\text{m}}\left[\dot{T}\left(0\right)-\Delta {\dot{T}}_{\text{в}}\right]=a_{\text{m}}\rho \left[\frac{d\dot{U}}{dx}\left(0\right)+\delta \frac{d\dot{T}}{dx}\left(0\right)\right];$ (12)

${\tilde{\alpha }}_{\text{w}}\left[\dot{T}\left(0\right)-\Delta {\dot{T}}_{\text{в}}\right]=\lambda \frac{d\dot{T}}{dx}\left(0\right).$ (13)

Здесь обозначено

${\tilde{\alpha }}_{\text{w}}={\alpha }_{\text{w}}+r(1-\gamma ){\tilde{\alpha }}_{\text{m}},$ $\Delta {\dot{T}}_{\text{в}}=\Delta T_{\text{в}}\exp \left(i{\psi }_{\text{в}}\right)$

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

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

$\left|\dot{T}\left(x\right)\right|\to 0$ и $\left|\dot{U}\left(x\right)\right|\to 0$ при $x\to \infty .$ (14)

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

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

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

$\nu =\left(1+a_{\text{w}}/a_{\text{m}}+\delta r\gamma /c\right)/2;$ (15)

$T_{\mathrm{1,2}}=\frac{r\gamma }{c}\frac{1}{1-\nu \mp \sqrt{{\nu }^2-a_{\text{w}}/a_{\text{m}}}};$ (16)

${\mu }_{\mathrm{1,2}}=-{\beta }_{\text{1,2}}(1+i),$ (17)

${\beta }_{\text{1,2}}=\sqrt{\omega /\left(2a_{\text{w}}\right)}\sqrt{\nu \pm \sqrt{{\nu }^2-a_{\text{w}}/a_{\text{m}}}}.$

Величина v является безразмерной, Т_1,2 имеет размерность °С, а размерностью  является 1/м. Первые две величины – постоянные, определяемые свойствами материала, а третья величина кроме свойств материала зависит еще от частоты $\omega .$ Отметим, что во всех формулах под $\sqrt{z}$ мы понимаем, как обычно, главное значение корня квадратного из комплексного числа z $(\mathrm{Re}\sqrt{z}\ge 0).$ Кроме указанных пяти величин, нахождение которых может быть произведено непосредственно по исходным данным, следует еще, согласно [11], определить две безразмерные и зависящие от частоты $\omega$ величины С_1,2, которые находятся как решение следующей системы линейных алгебраических уравнений:

$\left[T_1-a_{\text{m}}\rho {\mu }_1(1+\delta T_1)/{\tilde{\alpha }}_{\text{m}}\right]C_1+$ (18)

$\begin{array}{l}+\left[T_2-a_{\text{m}}\rho {\mu }_2(1+\delta T_2)/{\tilde{\alpha }}_{\text{m}}\right]C_2=\Delta {\dot{T}}_{\text{в}};\\ T_1(1-\lambda {\mu }_1/{\tilde{\alpha }}_{\text{w}})C_1+\\ +T_2(1-\lambda {\mu }_2/{\tilde{\alpha }}_{\text{w}})C_2=\Delta {\dot{T}}_{\text{в}}.\end{array}$

После нахождения С_1,2 решение задачи (10)–(14) для комплексов $\dot{T}$ и $\dot{U}$ в матричном виде может быть записано как

$||\begin{array}{c}\dot{T}\left(x\right)\\ \dot{U}\left(x\right)\end{array}||=C_1||\begin{array}{c}{T}_1\exp \left({\mu }_{1}x\right)\\ \exp \left({\mu }_{1}x\right)\end{array}||+$ (19)

$+C_2||\begin{array}{c}{T}_2\exp \left({\mu }_{2}x\right)\\ \exp \left({\mu }_{2}x\right)\end{array}||.$

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

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

Материалом полупространства будем считать глину, ее теплофизические характеристики имеют следующие значения [13]: $\lambda =\mathrm{0,93}$ Вт/(м ⋅ °С); $с=\mathrm{1,9}\cdot {10}^3$ Дж/(кг ⋅ °С); $\rho =\mathrm{1,5}\cdot {10}^3$ кг/м³; $\gamma =\mathrm{0,1};$ $\delta =\mathrm{1,5}\cdot {10}^{–3}$ 1/°С); $a_{\text{m}}=\mathrm{2,6}\cdot {10}^{–8}$ м²/с; $a_{\text{w}}=\mathrm{0,32}\cdot {10}^{–6}$ × 10^–6 м²/с.

Эти данные позволяют вычислить коэффициенты (15) и (16): $\nu =\mathrm{6,7432};$ ${Т}_1=–\mathrm{10,347}$°С;  ${Т}_2=\mathrm{12,790}\cdot {10}^3$°С.

Рассчитаем теперь круговую частоту колебаний $\omega =2\pi /Т,$ выбрав в качестве периода колебаний Т один год, а затем найдем коэффициенты (17): $\omega =\mathrm{1,9912}\cdot {10}^{–7}$ с^–1; ${\beta }_1=\mathrm{1,9725}$ 1/м; ${\beta }_2=\mathrm{0,55540}$ 1/м.

При расчете коэффициентов тепло- и массообмена будем опираться на фундаментальное руководство по этому вопросу [17], на статью авторов [14], где для этих коэффициентов получены приближенные формулы, а также на имеющиеся в [13] опытные данные. На основании этих работ в качестве характерных значений коэффициентов тепло- и массообмена в опытах, где влажная глина будет обдуваться потоком воздуха, примем ${\alpha }_{\text{w}}=5$ Вт/(м² ⋅ °С); ${\alpha }_{\text{m}}=5\cdot {10}^{–3}$ кг/(м²·с).

В качестве температуры, вблизи которой происходят колебания температуры воздуха, возьмем ${Т}_0=20$°С. Тогда расчет коэффициентов ${\tilde{\alpha }}_{\text{w}}$ и ${\tilde{\alpha }}_{\text{m}}$ приводит к таким значениям: ${\tilde{\alpha }}_{\text{w}}=\mathrm{19,5}$ Вт/(м² ⋅ °С); ${\tilde{\alpha }}_{\text{m}}=\mathrm{7,13}\cdot {10}^{–6}$ кг/(м²·с·°С).

Наконец, в качестве характеристик колебаний температуры воздуха в формуле (5) примем  $\Delta {Т}_{\text{в}}=5$°С, ${\psi }_{\text{в}}=0.$

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

$\begin{array}{l}(\mathrm{0,27501}+i\mathrm{10,622})C_1+(12851+i\mathrm{61,322})C_2=5;\\ (-\mathrm{11,320}-i\mathrm{0,97336})C_1+(13129+i\mathrm{338,78})C_2=5.\end{array}$

Здесь скобки и правые части уравнений имеют размерности °С. Решая эту систему, найдем, что

$C_1=\left|C_1\right|\exp \left(i\arg C_1\right);$ (20)

$\begin{array}{l}\left|C_1\right|={\mathrm{9,306510}}^{-3};\arg C_1=-\mathrm{0,024990};\\ C_2=\left|C_2\right|\exp \left(i\arg C_2\right);\\ \left|C_2\right|=\mathrm{0,38875}\cdot {10}^{-3};\arg C_2=-\mathrm{0,024554.}\end{array}$

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

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

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

$\begin{array}{l}\dot{U}\left(x\right)=C_1\exp \left({\mu }_{1}x\right)+C_2\exp \left({\mu }_{2}x\right)=\\ =\left|C_1\right|\exp \left(i\arg C_1\right)\exp \left[-{\beta }_1(1+i)x\right]+\\ +\left|C_2\right|\exp \left(i\arg C_2\right)\exp \left[-{\beta }_2(1+i)x\right].\end{array}$

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

$\begin{array}{l}U\left(x,\tau \right)-U_0=\\ =\left|C_1\right|\cdot \exp \left(-{\beta }_{1}x\right)\cdot \sin \left(\omega \tau +\arg C_1-{\beta }_{1}x\right)+\\ +\left|C_2\right|\cdot \exp \left(-{\beta }_{2}x\right)\cdot \sin \left(\omega \tau +\arg C_2-{\beta }_{2}x\right).\end{array}$ (21)

Каждое из двух слагаемых в правой части этой формулы представляет собой бегущую затухающую гармоническую волну. Рассмотрим первое слагаемое. Согласно терминологии, принятой в теории волн, величина ${\beta }_1$ есть коэффициент затухания этой первой волны влагосодержания, а обратная величина ${\Delta }_1=1/{\beta }_1$ – ее глубина проникновения. Другими характеристиками волны являются ее фазовая скорость $V_1=\omega /{\beta }_1$ и длина волны ${\Lambda }_1=2\pi /{\beta }_1.$ Для гармонической волны общего вида в формулах для фазовой скорости и длины волны вместо коэффициента затухания ${\beta }_1$ должен стоять коэффициент фазы ${{\beta }^{\text{'}}}_{\text{1}},$ который является коэффициентом при координате x под знаком синуса, но в нашем случае эти два коэффициента совпадают.

Поскольку коэффициент ${\beta }_1$ является функцией частоты $\omega$ (его указанное выше значение мы нашли при конкретной частоте), то фазовая скорость и длина волны также оказываются зависящими от частоты, т. е. распространение волн влагосодержания сопровождается дисперсией.

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

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

$A_1\left(x\right)=\left|C_1\right|\exp (-{\beta }_{1}x);{\phi }_1\left(x\right)=\arg C_1-{\beta }_{\text{1}}x;$ (22)

$A_2\left(x\right)=\left|C_2\right|\exp (-{\beta }_{2}x);{\phi }_2\left(x\right)=\arg C_2-{\beta }_{\text{2}}x.$

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

$U(x,\tau )-U_0=A\left(x\right)\sin \left[\omega \tau +\phi \left(x\right)\right];$ (23)

$A\left(x\right)=\sqrt{A_1^2+2A_1{A}_2\cos ({\phi }_1-{\phi }_2)+A_2^2},$ (24)

$\begin{array}{l}\sin \phi \left(x\right)=\frac{A_1\sin {\phi }_1+A_2\sin {\phi }_2}{A},\\ \cos \phi \left(x\right)=\frac{A_1\cos {\phi }_1+A_2\cos {\phi }_2}{A}.\end{array}$

Здесь величины А₁, А₂, ${\phi }_1,$ ${\phi }_2$ являются функциями координаты x и вычисляются по формулам (22).

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

$U(x,\tau )-U_0=$ (25)

$=\Delta U_{\text{п}}\exp (-x/{\Delta }_1)\sin \left[\omega (\tau -\Delta {\tau }_1)-{\beta }_{1}x\right],$

где $\Delta U_{\text{п}}=\left|C_1\right|=\mathrm{9,31}\cdot {10}^{-3}$ – амплитуда колебаний влагосодержания на поверхности x = 0; ${\Delta }_1=1/{\beta }_1=\mathrm{0,507}$ м – глубина проникновения волны влагосодержания; $\Delta {\tau }_1=-\left(\arg C_1\right)/\omega =\mathrm{1,45}$ сут – время запаздывания колебаний на поверхности x = = 0 относительно колебаний температуры воздуха.

Общее время запаздывания колебаний на произвольной глубине x будет определяться формулой

$\Delta {\tau }_1\left(x\right)=\Delta {\tau }_1+{\beta }_{1}x/\omega .$

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

${\beta }_{1}x/\omega =\mathrm{3,82}\text{\hspace{0.33em}мес}\gg \Delta {\tau }_1.$

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

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

$\dot{T}\left(x\right)=T_1{C}_1\exp \left({\mu }_{1}x\right)+T_2{C}_2\exp \left({\mu }_{2}x\right).$

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

$T_1\left|C_1\right|=-\mathrm{0,0963}°\text{C};T_2\left|C_2\right|=\mathrm{4,97}°\text{C.}$

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

$\begin{array}{l}T(x,\tau )-T_0=\\ =\Delta T_{\text{п}}\exp (-x/{\Delta }_2)\sin \left[\omega (\tau -\Delta {\tau }_2)-{\beta }_{2}x\right].\end{array}$

Здесь $\Delta T_{\text{п}}=T_2\left|C_2\right|=\mathrm{4,97}°\text{C}$ – амплитуда колебаний температуры на поверхности x = 0, которая мало отличается от амплитуды колебаний температуры воздуха $\Delta {Т}_{\text{в}}=5°\text{C;}$ ${\Delta }_2=1/{\beta }_2=\mathrm{1,80}$ м – глубина проникновения волны температуры; $\Delta {\tau }_2=-\left(\arg C_2\right)/\omega =\mathrm{1,43}$ сут – время запаздывания колебаний на поверхности x = 0 относительно колебаний температуры воздуха.

Общее время запаздывания колебаний на произвольной глубине x будет определяться формулой

$\Delta {\tau }_2\left(x\right)=\Delta {\tau }_2+{\beta }_{2}x/\omega .$

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

${\beta }_{2}x/\omega =\mathrm{4,30}\text{\hspace{0.33em}мес}\gg \Delta {\tau }_2.$

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

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

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

В простейшей модели (будем говорить, что это модель А) материал влаги не содержит, исследуется только поле температуры. Температурные волны для такого случая описываются известной в литературе формулой Фурье [6], которая в обозначениях, принятых в настоящей статье, выглядит так:

$T(x,\tau )=T_0+$ (26)

$+\Delta T_{\text{п}}\exp (-{\beta }_{\text{w}}x)\sin \left[\omega \tau +{\psi }_{\text{п}}-{\beta }_{\text{w}}x\right].$

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

${\beta }_{\text{w}}=\sqrt{\omega /\left(2a_{\text{w}}\right)}$ (27)

– коэффициент затухания волны, а $\Delta T_{\text{п}}$ и ${\psi }_{\text{п}}$ есть амплитуда и начальная фаза колебаний температуры на поверхности материала (в отличие от нашего рассмотрения, в задаче Фурье изменяется по гармоническому закону и считается заданной не температура воздуха, а температура поверхности x = 0).

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

Ниже для каждой из моделей мы приводим итоговые формулы для характеристик полей. В случае модели С, принимая сформулированные выше приближения, мы считаем волны затухающими гармоническими. Нижний индекс указывает, о каком поле, Т или U, идет речь в данной формуле, а верхний индекс определяет принадлежность формулы к моделям А, В или С. Формулы статьи [11] приведены в соответствие с обозначениями, принятыми в настоящей статье. Мы обозначили также $Lu=a_{\text{m}}/a_{\text{w}}$ – критерий Лыкова.

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

${\beta }_T^A={\beta }_T^B={\beta }_{\text{w}};$ ${\beta }_T^C={\beta }_{\text{w}}\sqrt{\nu -\sqrt{{\nu }^2-L{u}^{-1}}};$

${\beta }_U^B={\beta }_{\text{w}}\sqrt{L{u}^{-1}};$ ${\beta }_U^C={\beta }_{\text{w}}\sqrt{\nu +\sqrt{{\nu }^2-L{u}^{-1}}}.$

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

$\Delta {\tau }_T^A=\Delta {\tau }_T^B=\frac{{\beta }_{\text{w}}x}{\omega };$ $\Delta {\tau }_T^C=\frac{{\beta }_T^{C}x}{\omega };$

$\Delta {\tau }_U^B=\frac{{\beta }_{\text{w}}x}{\omega }\sqrt{L{u}^{-1}};$ $\Delta {\tau }_U^C=\frac{{\beta }_U^{C}x}{\omega }.$

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

$\Delta t_T^B=\frac{1}{\omega }\text{arctg}\left(\frac{\lambda {\beta }_{\text{w}}}{{\tilde{\alpha }}_{\text{w}}+\lambda {\beta }_{\text{w}}}\right);$ $\Delta t_T^C=-\frac{\arg C_2}{\omega };$

$\Delta t_U^B=\Delta t_T^B;$ $\Delta t_U^C=-\frac{\arg C_1}{\omega }.$

г) Амплитуды колебаний на поверхности $\Delta Т$ и $\Delta U\text{:}$

$\Delta {Т}^B={\tilde{\alpha }}_{\text{w}}\frac{\Delta {Т}_{\text{в}}}{\sqrt{{({\tilde{\alpha }}_{\text{w}}+\lambda {\beta }_{\text{w}})}^2+{\left(\lambda {\beta }_{\text{w}}\right)}^2}};$

$\Delta T^C=\frac{r\gamma }{c}\frac{1}{1-\nu +\sqrt{{\nu }^2-L{u}^{-1}}}\left|{С}_2\right|;$

$\Delta U^B=\frac{{\tilde{\alpha }}_{\text{m}}\lambda }{{\tilde{\alpha }}_{\text{w}}\rho \sqrt{a_{\text{w}}{a}_{\text{m}}}}\Delta T^B;$ $\Delta U^C=\left|C_1\right|.$

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

 

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

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

Модель	$\beta ,$ 1/м		$\Delta \tau ,$ мес		$\Delta t,$ сут		Амплитуды
	Т	U	Т	U	Т	U	$\Delta Т$	$\Delta 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

 

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

1. Как видно из таблицы, модели В и С дают практически одинаковые результаты, поэтому при расчетах можно ограничиться более простой моделью В, в которой не учитываются явления термодиффузии и внутреннего парообразования.
2. Учет влаги с помощью модели В, согласно пунктам а) и б), не изменяет формул для характеристик температурного поля, полученных в рамках модели А, т. е. формул для глубины проникновения ß и времени запаздывания колебаний $\Delta \tau$ на произвольной глубине относительно колебаний на поверхности. Следовательно, даже при наличии в почве влаги теория Фурье для тепловых волн должна приводить к удовлетворительному согласию с экспериментом.
3. Приведенные в таблице значения параметров ß и $\Delta \tau$ для поля температуры находятся в хорошем согласии с данными наблюдений на метеорологических станциях [6].

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

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

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

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

Заключение

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

About the authors

Anatoliy M. Afanasyev

Volgograd State University

Author for correspondence.
Email: a.m.afanasiev@yandex.ru
ORCID iD: 0000-0003-3703-3167
SPIN-code: 3749-8456
Scopus Author ID: 8317-433500
ResearcherId: A-8779-2018
https://volsu.ru/struct/institutes/ipt/itsecurity/employees/emp.php?id=000000889

Doctor of Technical Sciences, professor of the Department of Information Security, Institute of Priority Technologies

Russian Federation, 100, University Avenue, Volgograd, 400062

Yulia S. Bakhracheva

Volgograd State University

Email: bakhracheva@yandex.ru
ORCID iD: 0000-0002-0558-5736
SPIN-code: 5296-8187
Scopus Author ID: 140281
ResearcherId: KHU-6620-2024

Candidate of Technical Sciences, associate professor of the Department of Information Security, Institute of Priority Technologies

Russian Federation, 100, University Avenue, Volgograd, 400062

References

Supplementary files