Generalization of the Fourier problem of temperature waves in half-space

Full Text

Введение

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

Если иметь в виду потребности практики, то большое значение для проведения оценочных расчетов имеет теоретическое рассмотрение указанных волновых задач для областей с простой геометрией, а именно для областей в виде шара, цилиндра, пластины и полупространства. В теории электромагнетизма здесь можно указать на задачу о расчете коэффициентов отражения и пропускания при падении плоской гармонической электромагнитной волны на однородную пластину [1, с. 581–585], на задачу о выводе формул для коэффициентов отражения и проникновения (формул Френеля) при падении такой же волны на однородное полупространство [2, с. 162–172], а также на задачу для телеграфного уравнения при исследовании волн в конечных и полубесконечных линиях передачи [3, с. 557–568]; в теории тепломассопереноса – на статьи [4–6], где рассматривается сушка влажной пластины электромагнитными волнами при периодическом изменении интенсивности излучения; в теории теплопроводности – на классическое руководство по теоретической физике [7, с. 542–549], в котором исследуется прохождение тепловых волн сквозь пластину при гармонических граничных условиях 1-го и 3-го рода, и на монографию [8, с. 298–313], где при той же постановке проблемы исследуется еще и задача для полупространства, а также на известную в науке задачу Фурье о колебаниях температуры в поверхностном слое земной коры, вызванных сезонными изменениями температуры воздуха [9, с. 238–247]. Решения всех перечисленных выше задач находят важные применения на практике; в частности, решение задачи Фурье, которая ставится как задача для полупространства, принадлежит к числу основополагающих теоретических фактов в мерзлотоведении (геокриологии) и представляет собой важный инструмент при решении проблем метеорологии, климатологии и охраны окружающей среды, а также при строительстве зданий и развитии сельского хозяйства в области распространения мерзлых пород [10; 11]. Существенным недостатком формул Фурье и вытекающих из них законов Фурье, который ограничивает применение содержащихся в них результатов для исследования проблем геокриологии, является то, что они не учитывают наличия в почве влаги, ее перемещения под действием градиентов температуры и влагосодержания и испарения в толще материала и с его поверхности. Исправление этого недостатка и является целью предлагаемой статьи. Применяя теорию тепломассопереноса А.В. Лыкова [8; 12; 13] к задаче для полупространства с периодически изменяющимися граничными условиями, мы получим формулы для асимптотических по времени полей температуры и влагосодержания, в которых движение влаги и ее превращения будут корректно учтены. Такая задача в теории тепломассопереноса еще никем не рассматривалась.

Математическая модель тепломассообмена полупространства с воздушным потоком

Имея в виду задачу Фурье о колебаниях температуры и влагосодержания в поверхностном слое земной коры, рассмотрим показанное на рисунке однородное, содержащее влагу полупространство x > 0, граница которого x=0 обдувается воздухом, имеющим за пределами пограничного слоя температуру Т_в и влажность φ. Материал полупространства будем считать состоящим из твердой основы (капиллярно-пористое тело) и воды. Примем также, что интенсивность теплообмена Q и интенсивность массообмена J поверхности x=0 с воздушной средой слабо изменяются вдоль этой поверхности, т. е. эти величины зависят только от времени t. В описанной ситуации распределения температуры Т и влагосодержания U будут зависеть только от x и t, т. е. искомыми функциями будут Т(x, t) и U(x, t). Система уравнений и краевых условий для расчета этих функций будет иметь следующий вид [12; 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 };0<x<\infty ;$.     (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};0<x<\infty ;$.    (2)

$Q\left(\tau \right)+r\left(1-\gamma \right)\cdot J\left(\tau \right)=\lambda \frac{\partial T}{\partial x}\left(\mathrm{0,}\tau \right);x=0; \left(3\right)$

$J\left(\tau \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];x=0;$     (4)

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

$\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{в}}\right)\right];\\ P\left(T\right)=\mathrm{6,03}\cdot {10}^{-3}\exp \frac{\mathrm{17,3}T}{T+T_1}.\end{array}$       (6)

Соотношения (1) и (2) представляют собой уравнения распространения тепла и влаги в области, занятой материалом; уравнениями (3) и (4) задаются краевые условия на границе x = 0; формулами (5) и (6) определяются интенсивности тепло- и массообмена на этой границе (теплообмен по закону Ньютона и массообмен по закону Дальтона). В приведенных соотношениях: c, r, λ, g, a_m, δ – теплофизические характеристики материала, а именно удельная теплоемкость, плотность в сухом состоянии, коэффициент теплопроводности, критерий испарения, коэффициент диффузии влаги, относительный коэффициент термодиффузии влаги; a_w = λ/(cρ) – коэффициент диффузии тепла (коэффициент температуропроводности); r – удельная теплота парообразования воды; α_w и α_m – коэффициенты тепло- и массообмена поверхности образца с воздушной средой; Р(Т) – функция Г.К. Филоненко, моделирующая зависимость относительного парциального давления насыщенного водяного пара от его температуры Т при общем нормальном давлении; Т₁ = 238 °С – постоянная.

По причине, о которой будет сказано ниже, вводить в рассмотрение начальные условия для функций Т и U мы не будем.

Постановка задачи об асимптотике полей тепломассопереноса

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

Пусть теперь, начиная с момента t = 0, температура воздуха Т_в начинает совершать малые колебания вблизи температуры Т₀. При малых отклонениях температуры поверхности Т(0, t) и температуры воздуха Т_в(t) от фиксированной температуры Т₀ зависимость (6) можно линеаризовать, разложив функцию Р(Т) в ряд Тейлора в окрестности точки Т₀. Сделав это, и приняв φ = 1, вместо исходной формулы (6) для интенсивности массообмена получим приближенную формулу

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

где ${\alpha }_{\text{t}}={\alpha }_{\text{m}}\frac{dP}{dT}\left(T_0\right)$ – коэффициент массообмена по перепаду температуры.

Представление функции J(t) в виде (7) превращает введенную нами систему уравнений в линейную систему, на что в дальнейшем мы будем существенным образом опираться.

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

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

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

$\left\{\begin{array}{l}T\left(x,\tau \right)=T_0+T\left(x\right)\sin \left(\omega \tau +{\psi }_{\text{t}}\left(x\right)\right),\\ U\left(x,\tau \right)=U_0+U\left(x\right)\sin \left(\omega \tau +{\psi }_{\text{u}}\left(x\right)\right),\end{array}\right.$    (9)

в котором и температура, и влагосодержание материала, так же как и температура воздуха, совершают при каждом фиксированном x малые гармонические колебания вблизи равновесных значений Т₀ и U₀. Действительно, вычислив с помощью (9) и (8) разность $\left[T\left(\mathrm{0,}\tau \right)-T_{\text{в}}\left(\tau \right)\right],$ и подставив ее в формулу (7), мы получим

$\begin{array}{l}J\left(\tau \right)=\\ ={\alpha }_{\text{t}}\left[T\left(0\right)\sin \left(\omega \tau +{\psi }_{\text{t}}\left(0\right)\right)-\Delta T_{\text{в}}\sin \left(\omega \tau +{\psi }_{\text{в}}\right)\right].\end{array}$    (10)

Таким образом, функция J(t) оказывается гармонической. Аналогичным образом, на основании формулы (5), доказывается и гармоничность функции Q(t). Но ведь и все другие члены в уравнениях (1)–(4), в соответствии с (9), также будут гармоническими, и потому этим уравнениям можно будет попытаться удовлетворить, подобрав должным образом зависимости

$T\left(x\right),U\left(x\right),{\psi }_{\text{t}}\left(x\right),{\psi }_{\text{u}}\left(x\right)$

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

$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 .$.      (11)

Сформулированная задача о подборе функций (9) относится к числу задач без начальных данных; ее решение дает асимптотику полей Т и U при t → ∞.

Постановка задачи для комплексов гармонических полей

Искомые функции нашей задачи, Т(x, t) и U(x, t), определяются формулами (9). Введем вместо них новые искомые функции, Т*(x, t) и U*(x, t), и, обращаясь к методу комплексных амплитуд, сопоставим этим новым функциям их комплексы  по следующим правилам:

$\left\{\begin{array}{l}T*\left(x,\tau \right)=T\left(x,\tau \right)-T_0=\\ =T\left(x\right)\sin \left(\omega \tau +{\psi }_{\text{t}}\left(x\right)\right)\leftrightarrow \dot{T}\left(x\right)=\\ =T\left(x\right)\exp \left(i{\psi }_{\text{t}}\left(x\right)\right);\\ U*\left(x,\tau \right)=U\left(x,\tau \right)-U_0=\\ =U\left(x\right)\sin \left(\omega \tau +{\psi }_{\text{u}}\left(x\right)\right)\leftrightarrow \dot{U}\left(x\right)=\\ =U\left(x\right)\exp \left(i{\psi }_{\text{u}}\left(x\right)\right).\end{array}\right.$       (12)

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

 $\left\{\begin{array}{l}i\omega \dot{T}\left(x\right)=a_{\text{w}}\frac{d^2\dot{T}\left(x\right)}{d{x}^2}+\frac{r\gamma }{c}i\omega \dot{U}\left(x\right);\\ 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}.\end{array}\right.$       (13)

Здесь мы имеем систему двух обыкновенных дифференциальных уравнений второго порядка относительно двух неизвестных комплекснозначных функций действительной переменной. Система линейная, однородная, с постоянными коэффициентами, неопределенная. Для выделения единственного решения этой системы обратимся к оставшимся уравнениям задачи. В области комплексов уравнения (3) и (4) будут выглядеть так:

$\begin{array}{l}\dot{Q}+r\left(1-\gamma \right)\dot{J}=\lambda \frac{d\dot{T}}{dx}\left(0\right);\\ \dot{J}=a_{\text{m}}\rho \left[\frac{d\dot{U}}{dx}\left(0\right)+\delta \frac{d\dot{T}}{dx}\left(0\right)\right].\end{array}$      (14)

Образуя комплексы от обеих частей (10), и приравнивая их, с учетом первой из формул (12) получим

$\dot{J}={\alpha }_{\text{t}}\left[\dot{T}\left(0\right)-\Delta {\dot{T}}_{\text{в}}\right]$,

где $\Delta {\dot{T}}_{\text{в}}=\Delta T_{\text{в}}\exp \left(i{\psi }_{\text{в}}\right)$ – заданное комплексное число. Аналогичным образом, обращаясь к формуле (5), для комплекса интенсивности теплообмена будем иметь

$\dot{Q}={\alpha }_{\text{w}}\left[\dot{T}\left(0\right)-\Delta {\dot{T}}_{\text{в}}\right].$

Подставляя полученные выражения для $\dot{J}$ и $\dot{Q}$ в (14), после преобразований получим систему двух уравнений следующего вида:

$\left\{\begin{array}{l}\dot{T}\left(0\right)-\frac{\lambda }{{\alpha }_{\text{e}}}\frac{d\dot{T}}{dx}\left(0\right)=\Delta {\dot{T}}_{\text{в}};\\ \dot{T}\left(0\right)-\frac{a_{\text{m}}\rho }{{\alpha }_{\text{t}}}\left[\frac{d\dot{U}}{dx}\left(0\right)+\delta \frac{d\dot{T}}{dx}\left(0\right)\right]=\Delta {\dot{T}}_{\text{в}}.\end{array}\right.$    (15)

Здесь ${\alpha }_{\text{e}}={\alpha }_{\text{w}}+r\left(1-\gamma \right){\alpha }_{\text{t}}$ – эффективная интенсивность теплообмена.

Построив общее решение системы (13), т. е. записав общие выражения для функций $\dot{T}\left(x\right)$ и $\dot{U}\left(x\right)$, мы затем из системы (15), которая играет роль краевых условий, найдем входящие в эти общие выражения произвольные постоянные.

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

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

$\dot{T}\left(x\right)=T\exp \left(\mu x\right),\dot{U}\left(x\right)=U\exp \left(\mu x\right),$       (16)

где Т, U, μ – комплексные постоянные (метод Эйлера, [15]). Подставляя эти выражения для комплексов в (13), после сокращений получим квадратную систему линейных алгебраических уравнений для нахождения введенных постоянных:

$\left\{\begin{array}{l}\left(a_{\text{w}}{\mu }^2-i\omega \right)T+\frac{r\gamma }{c}i\omega U=0;\\ a_{\text{m}}\delta {\mu }^{2}T+\left(a_{\text{m}}{\mu }^2-i\omega \right)U=0.\end{array}\right.$    (17)

Чтобы решение Т, U этой однородной системы было нетривиальным, ее определитель должен равняться нулю, что приводит к характеристическому уравнению

$\left(a_{\text{w}}{\mu }^2-i\omega \right)\left(a_{\text{m}}{\mu }^2-i\omega \right)-\frac{a_{\text{m}}\delta r\gamma }{c}i\omega {\mu }^2=0.$      (18)

Это уравнений четвертой степени относительно μ, и оно является биквадратным. Раскрывая скобки, и вводя новую безразмерную постоянную материала

$\nu =\frac{1}{2}\left(1+\frac{a_{\text{w}}}{a_{\text{m}}}+\frac{\delta r\gamma }{c}\right),$    (19)

приведем уравнение (18) к виду

${\mu }^4-\frac{2\nu }{a_{\text{w}}}i\omega {\mu }^2-\frac{{\omega }^2}{a_{\text{w}}{a}_{\text{m}}}=0.$

Вводя новую переменную $\xi ={\mu }^2$, перепишем это уравнение так:

${\xi }^2-\frac{2\nu }{a_{\text{w}}}i\omega \xi -\frac{{\omega }^2}{a_{\text{w}}{a}_{\text{m}}}=0.$

Решая это уравнение, найдем, что

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

Вспоминая, что $\xi ={\mu }^2,$ получим в итоге для характеристического уравнения (18) четыре корня

${\mu }_{\mathrm{1,}\mathrm{2,}\mathrm{3,}4}=\pm \sqrt{{\xi }_{\mathrm{1,}2}}.$

Отметим, что в последних двух формулах под $\sqrt{z}$ мы понимаем, как обычно, главное значение корня квадратного из комплексного числа z $(\mathrm{Re}\sqrt{z}\ge 0).$

Выполним теперь условие на бесконечности (11). Согласно формулам (16), для их выполнения необходимо иметь Re μ < 0. Следовательно, с учетом свойства главного значения корня, из четырех корней характеристического уравнения нам подойдут только два корня ${\mu }_{\mathrm{1,}2}=-\sqrt{{\xi }_{\mathrm{1,}2}},$ или

$\begin{array}{l}{\mu }_1=-\sqrt{\frac{\omega \nu }{a_{\text{w}}}\left(i+\sqrt{\frac{a_{\text{w}}}{a_{\text{m}}{\nu }^2}-1}\right)},\\ {\mu }_2=-\sqrt{\frac{\omega \nu }{a_{\text{w}}}\left(i-\sqrt{\frac{a_{\text{w}}}{a_{\text{m}}{\nu }^2}-1}\right)}.\end{array}$     (20)

Используя найденные корни характеристического уравнения (18), вернемся к системе (17) и найдем коэффициенты Т и U. Поскольку определитель системы равен нулю, уравнения линейно зависимы, и значит, одно из них можно отбросить; оставшееся уравнение будет иметь бесконечно много решений, и потому одно из двух неизвестных можно положить равным любому числу. Имея это в виду, отбросим второе уравнение, а в первом уравнении положим U = 1. Тогда, приняв в этом уравнении сначала μ = μ₁, а затем μ = μ₂, найдем сначала коэффициент Т₁, а затем коэффициент Т₂. Результат будет таким:

$U_{\mathrm{1,}2}=1;T_{\mathrm{1,}2}=\frac{r\gamma }{c\left(i\omega -a_{\text{w}}{{\mu }_{\mathrm{1,}2}}^2\right)}i\omega .$

Подставляя сюда μ_1, 2 из (20), перепишем, после преобразований, эти формулы так:

$U_{\mathrm{1,}2}=1;T_{\mathrm{1,}2}=\frac{r\gamma }{c}\frac{1}{\left(1-\nu \right)\pm i\nu \sqrt{\frac{a_{\text{w}}}{a_{\text{m}}{\nu }^2}-1}}.$       (21)

Определив корни характеристического уравнения μ_{1, 2} и коэффициенты U_{1, 2}, Т_{1, 2}, мы тем самым построили фундаментальную систему решений

$\left(T_n\exp \left({\mu }_{n}x\right),U_n\exp \left({\mu }_{n}x\right)\right),n=\mathrm{1,}2$

для системы дифференциальных уравнений (13). Тогда общее решение этой системы в матричном виде будет выглядеть так:

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

Здесь μ_{1, 2} вычисляются по формулам (20), Т_{1, 2} – по формулам (21), а С_{1, 2} – произвольные постоянные.

Для нахождения постоянных С_{1, 2} обратимся к краевым условиям (15). Используя (22), вычислим входящие в (15) величины:

$\begin{array}{l}\dot{T}\left(0\right)=C_1{T}_1+C_2{T}_2;\\ \frac{d\dot{T}}{dx}\left(0\right)=C_1{T}_1{\mu }_1+C_2{T}_2{\mu }_2;\\ \frac{d\dot{U}}{dx}\left(0\right)=C_1{\mu }_1+C_2{\mu }_2.\end{array}$

Подставляя эти значения в (15), получим систему линейных алгебраических уравнений для нахождения С_{1, 2}:

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

Определив отсюда С_{1, 2} и подставив эти коэффициенты в (22), будем иметь решение исходной задачи для комплексов, т. е. решение системы дифференциальных уравнений (13) с краевыми условиями (15). После этого останется только перейти от комплексов к оригиналам, т. е. представить полученные решения в виде (9).

Пример конкретного расчета по указанному алгоритму будет произведен в следующем пункте.

Асимптотическое решение для математической модели тепломассопереноса с δ = 0 и γ = 0

Реализация представленного выше расчетного алгоритма представляет собой в общем случае непростую вычислительную работу, сложность которой зависит, прежде всего, от принятой математической модели тепломассопереноса. Здесь мы ограничимся расчетами в рамках одной из самых простых математических моделей, в которой полагают δ = 0 (пренебрегают термодиффузией, т. е. движение влаги к поверхности происходит только за счет перепада влагосодержания) и γ = 0 (пренебрегают внутренним парообразованием, т. е. превращение воды в пар происходит только на поверхности). Условия применимости такой упрощенной модели к задачам теории тепломассопереноса и анализ получаемых при сделанных приближениях решений можно найти в работах [4–6; 16]. Как видно из (1) и (2), уравнения для температуры и влагосодержания в этом частном случае оказываются независимыми (связь между функциями Т и U осуществляется в этой модели только через граничные условия), что вносит существенные упрощения в алгоритм исследования процесса. Именно, отпадает необходимость в использовании теории систем дифференциальных уравнений, на которую мы опирались в предыдущем пункте. Уравнения для комплексов (13) в принявших более простой вид условиях нашей задачи станут такими:

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

Пользуясь методом Эйлера, найдем общие решения этих уравнений, удовлетворяющие условию на бесконечности (11). Они будут выглядеть так:

$\dot{T}\left(x\right)=C_1\exp \left({\mu }_{1}x\right),\dot{U}\left(x\right)=C_2\exp \left({\mu }_{2}x\right),$     (24)

где

$\begin{array}{l}{\mu }_1=-\sqrt{\frac{\omega }{a_{\text{w}}}}\left(\frac{1}{\sqrt{2}}+i\frac{1}{\sqrt{2}}\right),\\ {\mu }_2=-\sqrt{\frac{\omega }{a_{\text{m}}}}\left(\frac{1}{\sqrt{2}}+i\frac{1}{\sqrt{2}}\right),\end{array}$         (25)

а С_{1, 2} – произвольные постоянные. Для их нахождения воспользуемся краевыми условиями (15). Положив в них δ = 0, и подставив в получившиеся уравнения найденные выше общие решения  и  будем иметь:

$\left\{\begin{array}{l}{C}_1-\lambda C_1{\mu }_1/{\alpha }_{\text{e}}=\Delta {\dot{T}}_{\text{в}};\\ C_1-a_{\text{m}}\rho C_2{\mu }_2/{\alpha }_{\text{t}}=\Delta {\dot{T}}_{\text{в}}.\end{array}\right.$

Решая эту систему, после преобразований найдем:

 $C_1=\frac{\Delta {\dot{T}}_{\text{в}}}{1-\lambda {\mu }_1/{\alpha }_{\text{e}}},C_2=k{C}_1,$    (26)

где $k=\frac{{\alpha }_{\text{t}}\lambda }{{\alpha }_{\text{e}}\rho \sqrt{a_{\text{w}}{a}_{\text{m}}}}.$

Итак, решение нашей задачи для комплексов дается формулами (24), в которых коэффициенты определяются формулами (25) и (26). Отметим, что обращение в нуль знаменателя в первой из формул (26) исключено, потому что коэффициент μ₁ является комплексным.

Теперь мы должны вернуться от комплексов (изображений) к исходным гармоническим функциям времени (оригиналам). В методе комплексных амплитуд этот этап расчетов соответствует нахождению обратного преобразования Фурье при решении задач спектральным методом. Найдем сначала оригинал для комплекса $\dot{T}\left(x\right).$ Введем предварительно новое обозначение. Постоянную μ₁ в формулах (25) запишем в виде

${\mu }_1=-{\beta }_{\text{w}}\left(1+i\right),$          (27)

где ${\beta }_{\text{w}}=\sqrt{\omega /\left(2a_{\text{w}}\right)}$ – коэффициент затухания тепловых волн. Обратную величину Δ_w = 1/β_w назовем глубиной проникновения тепловых волн. Эти названия заимствованы из теории электромагнитных волн. Правомерность обращения к таким терминам будет обоснована ниже. Переход к оригиналу основывается на правиле образования комплексов (12) и осуществляется по следующей схеме:

$\begin{array}{l}\dot{T}\left(x\right)=C_1\exp \left({\mu }_{1}x\right)=\\ =\left|C_1\right|\exp \left(i\arg C_1\right)\exp \left[-{\beta }_{\text{w}}\left(1+i\right)x\right]\leftrightarrow \\ \leftrightarrow \left|C_1\right|\exp \left(-{\beta }_{\text{w}}x\right)\sin \left(\omega \tau +\arg C_1-{\beta }_{\text{w}}x\right).\end{array}$

Основную трудность, как это видно из формул (25) и (26), здесь представляет приведение к показательному виду комплексной постоянной С₁, т. е. расчет модуля этой постоянной  и ее аргумента arg C₁. Проделав такие расчеты, и прибавив, в соответствии с (12), к получившейся гармонической функции постоянную Т₀, получим искомое асимптотическое решение для поля температуры:

$\begin{array}{l}T\left(x,\tau \right)=\\ =T_0+\frac{\Delta T_{\text{в}}{\alpha }_{\text{e}}}{\sqrt{{\left({\alpha }_{\text{e}}+\lambda {\beta }_{\text{w}}\right)}^2+{\left(\lambda {\beta }_{\text{w}}\right)}^2}}\exp \left(-{\beta }_{\text{w}}x\right)\times \\ \times \sin \left[\omega \tau +{\psi }_{\text{в}}-\text{arctg}\left(\frac{\lambda {\beta }_{\text{w}}}{{\alpha }_{\text{e}}+\lambda {\beta }_{\text{w}}}\right)-{\beta }_{\text{w}}x\right].\end{array}$      (28)

Если принять во внимание, как изменяется амплитуда этой синусоиды с изменением координаты x, то станет понятным, почему для величины β_w мы применили указанное выше название.

Проделав аналогичные выкладки, для асимптотического поля влагосодержания получим:

$\begin{array}{l}U\left(x,\tau \right)=U_0+\frac{k\Delta T_{\text{в}}{\alpha }_{\text{e}}}{\sqrt{{\left({\alpha }_{\text{e}}+\lambda {\beta }_{\text{w}}\right)}^2+{\left(\lambda {\beta }_{\text{w}}\right)}^2}}\exp \left(-{\beta }_{\text{m}}x\right)\times \\ \times \sin \left[\omega \tau +{\psi }_{\text{в}}-\text{arctg}\left(\frac{\lambda {\beta }_{\text{w}}}{{\alpha }_{\text{e}}+\lambda {\beta }_{\text{w}}}\right)-{\beta }_{\text{m}}x\right].\end{array}$       (29)

Здесь ${\beta }_{\text{m}}=\sqrt{\omega /\left(2a_{\text{m}}\right)}$ – коэффициент затухания для волн влагосодержания.

Формулы (28) и (29) и представляют собой решение поставленной в этом пункте частной задачи. Входящие в них коэффициенты затухания β_w и β_m являются функциями частоты ω, а все другие коэффициенты являются определенными выше постоянными.

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

Обсуждение проведем на примере поля температуры. Согласно формуле (28), если температура воздуха, обдувающего влажное полупространство, совершает вблизи фиксированной температуры Т₀ малые гармонические колебания по закону (8), который характеризуется, кроме значения Т₀, еще и параметрами  $\Delta T_{\text{в}},$ $\omega$ и ${\psi }_{\text{в}},$ то в установившемся режиме колебания температуры вблизи равновесного значения Т₀ в любом месте полупространства оказываются также гармоническими функциями времени, причем имеют место следующие закономерности:

1) зависимость амплитуды этих колебаний ΔТ от частоты ω и глубины x определяется формулой

 $\begin{array}{l}\Delta T=\frac{\Delta T_{\text{в}}{\alpha }_{\text{e}}}{\sqrt{{\left({\alpha }_{\text{e}}+\lambda {\beta }_{\text{w}}\right)}^2+{\left(\lambda {\beta }_{\text{w}}\right)}^2}}\exp \left(-{\beta }_{\text{w}}x\right),\\ {\beta }_{\text{w}}=\sqrt{\omega /\left(2a_{\text{w}}\right)};\end{array}$        (30)

2) время запаздывания Dt максимумов (минимумов) температуры в материале от соответствующих моментов для температуры воздуха изменяется в зависимости от глубины x следующим образом:

$\Delta \tau =\frac{1}{\omega }\text{arctg}\left(\frac{\lambda {\beta }_{\text{w}}}{{\alpha }_{\text{e}}+\lambda {\beta }_{\text{w}}}\right)+\frac{1}{\omega }{\beta }_{\text{w}}x.$         (31)

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

Сформулированные нами законы изменения во времени полей температуры и влагосодержания в полупространстве с обдуваемой воздухом границей можно рассматривать как обобщение известных в литературе законов Фурье [9; 10], которые справедливы лишь в ситуации, когда, во-первых, материал полупространства не содержит влаги, и, значит, рассмотрение ведется только для поля температуры, и, во-вторых, изменяется по гармоническому закону не температура воздуха, а температура границы полупространства (принимаются граничные условия не 3-го рода, а имеющие ограниченное применение при моделировании процессов тепломассопереноса граничные условия 1-го рода). Этими обстоятельствами существенно ограничивается область использования рассматриваемых в литературе законов Фурье для решения проблем геокриологии.

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

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

Развитием предложенных здесь идей может быть построение аналогичных решений не для полупространства, а для пластины. Прохождение гармонических тепловых волн сквозь пластину (важная задача при исследовании теплоизолирующих свойств строительных конструкций) рассмотрена в монографии [7]. Мы планируем обобщить это решение, дополнив тепловые волны волнами влагосодержания. Задача для пластины является актуальной также и в теории электромагнитной сушки, где режимы с периодически изменяющейся интенсивностью излучения используют для организации «щадящих» режимов удаления влаги из термолабильных материалов [4–6]. Исследование таких режимов в теории сушки традиционно производят методом преобразования Лапласа, и такие исследования представляют собой в общем случае незначительный интерес, так как приводят к крайне сложным расчетам [4; 13]. Расчет можно существенным образом упростить, если исключить из таких решений входящий в них переходной процесс, т. е. начальный период сушки. Правомерность пренебрежения переходными процессами в задачах сушки обосновывалась авторами ранее [18]. Несложно понять, что исследование режима установившегося представляет собой намного более простую задачу, чем исследование процесса на всем его протяжении, и такое исследование может быть проведено по алгоритму настоящей статьи.

Заключение

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

Финансирование

Исследование выполнено при финансовой поддержке РФФИ и Администрации Волгоградской области в рамках научного проекта № 19-48-340015 р_а.

The reported study was funded by RFBR and the government of Volgograd region according to the research project no. 19-48-340015 р_а.

About the authors

Anatoly M. Afanasyev

Email: a.m.afanasiev@yandex.ru

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

Yulia S. Bakhracheva

Author for correspondence.
Email: bakhracheva@yandex.ru

Candidate of Technical Sciences, Associate Professor of the Department of Information Security, Institute of Priority Technologies

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. Airflow half space

Download (256KB)

Indexing metadata