Обобщение задачи Фурье о температурных волнах в полупространстве

Полный текст

Введение

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

Если иметь в виду потребности практики, то большое значение для проведения оценочных расчетов имеет теоретическое рассмотрение указанных волновых задач для областей с простой геометрией, а именно для областей в виде шара, цилиндра, пластины и полупространства. В теории электромагнетизма здесь можно указать на задачу о расчете коэффициентов отражения и пропускания при падении плоской гармонической электромагнитной волны на однородную пластину [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 р_а.

Об авторах

Анатолий Михайлович Афанасьев

Email: a.m.afanasiev@yandex.ru

доктор технических наук, доцент, профессор кафедры информационной безопасности Института приоритетных технологий 

Юлия Сагидулловна Бахрачева

Автор, ответственный за переписку.
Email: bakhracheva@yandex.ru

кандидат технических наук, доцент кафедры информационной безопасности Института приоритетных технологий 

Список литературы

Дополнительные файлы

Доп. файлы

Действие

1. JATS XML

Скачать

2. Рис. Обдуваемое воздушным потоком полупространство

Скачать (256KB)

Метаданные