Теплофизика высоких температур, 2022, T. 60, № 2, стр. 249-259
Точное аналитическое решение нестационарной задачи расчета полей температур и механических напряжений в двухслойном цилиндре и его применение для верификации численных моделей
Ю. Н. Токарев 1, Е. В. Моисеенко 1, *, Н. И. Дробышевский 1, Р. А. Бутов 1
1 ФГБУН Институт проблем безопасного развития атомной энергетики РАН (ИБРАЭ РАН)
Москва, Россия
* E-mail: moi@ibrae.ac.ru
Поступила в редакцию 21.12.2020
После доработки 20.05.2021
Принята к публикации 28.09.2021
- EDN: CTTJEH
- DOI: 10.31857/S0040364422020302
Аннотация
Доказана самосопряженность оператора теплопроводности в пространстве кусочно-дифференцируемых функций, заданных в неоднородных областях с разрывным коэффициентом теплопроводности, для которых на линии раздела областей формально выполняется условие идеального теплового контакта. Это позволило применить метод разделения переменных и получить в виде функционального ряда точное аналитическое решение нестационарной задачи расчета полей температур в двухслойном цилиндре с источником тепла. Полученное решение использовано для расчета термомеханических напряжений. Решение может применяться для оценки температур и напряжений в объектах, которые могут быть приближенно представлены в виде двухслойного цилиндра, а также для верификации численных моделей, применяемых для расчета термомеханических процессов в объектах с разрывными физическими свойствами. Представленное решение явилось следствием необходимости верификации кода FENIA на задаче остывания бидона с остеклованными радиоактивными отходами. Показано полное совпадение численного и аналитического решений, в том числе вблизи точки разрыва на границе сред.
ВВЕДЕНИЕ
Как правило, процессы, связанные с эволюцией полей температур и механических напряжений, моделируются численно, с применением соответствующих программных средств [1, 2]. В первую очередь это связано с отсутствием точных аналитических моделей, описывающих сложные объекты – неоднородные среды, в которых проходят нестационарные процессы. Даже при наличии таких моделей (см., например, [3]) они не могут с такой же точностью, как численные модели, учесть все детали и особенности описываемого объекта.
Тем не менее аналитические модели нельзя сбрасывать со счетов. Прежде всего это связано с их применением для получения предварительных приближенных оценок моделируемых процессов. Они позволяют получить качественное и приближенное количественное в${{{\text{и}}}^{'}}$дение протекания процесса без построения численной модели, которое само по себе может быть достаточно трудоемким. Еще одним важным применением является верификация программных средств численного моделирования. Даже при наличии современных подходов к верификации, основанных, например, на методе сконструированных решений [4], аналитические модели, описывающие объекты, сходные с реальными, сохраняют свою ценность. Они позволяют решать не только задачу верификации (проверки корректности решения уравнения в расчетном коде), но отчасти и валидации (проверки применимости кода для решения постулируемых задач [5]) с учетом, в частности, реалистичных начальных и граничных условий и источниковых членов. Кроме того, точные аналитические решения позволяют в рамках верификации исследовать отдельные особенности численных решений (пространственную и временную сходимость, учет разрывов в решениях, отсутствие влияния преобразований координат на точность решения и т.д.) [6].
При моделировании термомеханических процессов многие объекты реального мира могут быть приближенно описаны как двухслойные цилиндры с внутренним тепловыделением и внешним охлаждением. Примерами таких объектов могут служить тепловыделяющий элемент в активной зоне ядерного реактора, металлический бидон с находящимися внутри остеклованными радиоактивными отходами, изолированный проводник с протекающим по нему электрическим током и т.д. Вследствие различия теплофизических свойств материалов в таких объектах в результате теплового расширения могут возникать значительные механические напряжения. Эти напряжения будут терпеть разрыв на границе соприкосновения материалов, который невозможно измерить экспериментально и сложно оценить численно. В зависимости от свойств материалов такие напряжения могут приводить к различным повреждениям: образованию зазоров между слоями, растрескиванию внешнего слоя. Т.е. при численном моделировании критически важно корректно оценивать напряжения на границе материалов. Именно аналитические решения, будучи своего рода идеализированным предельным случаем с точным выделением разрыва на границе материалов, являются наиболее надежной основой для верификации численных моделей, позволяя, в частности, отслеживать возможные проблемы, связанные с пространственной дискретизацией.
В настоящее время имеются решения стационарной и нестационарной задач для однородного цилиндра [7, 8]. Также получено решение нестационарной задачи для двухслойного цилиндра без внутренних источников тепла [9]. Задача нахождения нестационарного температурного поля в многослойной среде в приближении быстро, но непрерывно меняющихся в пространстве физических свойств решена аналитически в [10]. В данной работе предлагается математический аппарат, который позволит рассчитывать поля температур и вызванных тепловым расширением механических напряжений для двухслойного цилиндра с учетом как граничных условий, так и источников тепла в цилиндре. При этом источники могут изменяться в пространстве и во времени.
НЕСТАЦИОНАРНОЕ ТЕМПЕРАТУРНОЕ ПОЛЕ В ДВУХСЛОЙНОМ ЦИЛИНДРЕ
Рассмотрим нестационарное температурное поле в двухслойном цилиндре с внешним радиусом ${{r}_{2}}$, схематично представленном на рис. 1. Радиус внутреннего цилиндра обозначен ${{r}_{1}}$. При наличии внутреннего источника тепла справедливо уравнение
(1)
$\rho c\frac{{\partial T}}{{\partial t}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {kr\frac{{\partial T}}{{\partial r}}} \right) + q\left( {r,t} \right),$
Функция $T\left( {r,t} \right)$ принимается кусочно-непрерывно дифференцируемой, и на поверхности соприкосновения внутреннего и внешнего цилиндров ($r = {{r}_{1}}$) выполняется условие идеального теплового контакта:
(2)
$ - {{\left. {{{k}_{1}}\frac{{\partial T}}{{\partial r}}} \right|}_{{r = {{r}_{1}} - 0}}} = - {{\left. {{{k}_{2}}\frac{{\partial T}}{{\partial r}}} \right|}_{{r = {{r}_{1}} + 0}}}.$На границе внешнего цилиндра задается граничное условие третьего рода:
(3)
$ - {{k}_{2}}{{\left. {\frac{{\partial T}}{{\partial r}}} \right|}_{{r = {{r}_{2}} - 0}}} = \alpha \left( {{{{\left. T \right|}}_{{r = {{r}_{2}}}}} - {{T}_{L}}} \right),$(4)
${{\left. T \right|}_{{t = 0}}} = {{T}_{0}}\left( r \right),\,\,\,\,0 \leqslant r \leqslant {{r}_{2}}.$Докажем самосопряженность оператора теплопроводности для разрывного коэффициента теплопроводности в общем трехмерном случае. Из его самосопряженности будет вытекать [7] взаимная ортогональность собственных функций, соответствующих разным собственным значениям, что позволит разложить решение уравнения (1) по этим функциям в ряд Фурье. Оператор теплопроводности обозначен
(5)
$L = - \nabla \cdot k\nabla ,\,\,\,\,\nabla = i\frac{\partial }{{\partial x}} + j\frac{\partial }{{\partial x}} + k\frac{\partial }{{\partial x}}.$Пусть коэффициент теплопроводности $k$ в формуле (5) является кусочно-непрерывной функцией, терпящей разрыв на границе материалов. В [7] самосопряженность $L$ доказывается для случая, когда коэффициент теплопроводности является непрерывной функцией координат. Докажем, что оператор теплопроводности остается самосопряженным и в случае, когда коэффициент теплопроводности терпит разрыв на границе $S$ двух областей G1 и G2, $G = {{G}_{1}} \cup {{G}_{2}}$ (рис. 2), но при этом на S выполняется условие идеального теплового контакта.
Рис. 2.
Область $G = {{G}_{1}} \cup {{G}_{2}}$, на внутренней границе $S$ выполняются условия сопряжения.
Пусть в области $G$ задана кусочно-непрерывная вместе со своей производной функция $k\left( x \right)$, предствляющая собой коэффициент теплопроводности и терпящая разрыв на поверхности $S$, а на границе ${\text{Г}}$ непрерывная функция $\alpha \left( x \right)$ – коэффициент теплоотдачи. Через ${{U}_{S}}$ обозначено пространство непрерывных функций, имеющих кусочно-непрерывные первые и вторые производные, терпящие разрывы первого рода на границе $S$, для которых выполняется условие
(6)
$k(x)\nabla {{\left. {f{\text{(}}x{\text{)}}} \right|}_{{{{S}^{ - }}}}} = k(x)\nabla {{\left. {f{\text{(}}x{\text{)}}} \right|}_{{{{S}^{ + }}}}}{\text{,}}$(7)
$k(x)\nabla {{\left. f \right|}_{{{{{\text{Г}}}^{ - }}}}} = \alpha {\text{(}}x{\text{)}}{{\left. {f(x)} \right|}_{{\text{Г}}}}{{N}^{ - }}.$Требуется доказать, что оператор теплопроводности $L = \nabla \cdot k\left( x \right)\nabla $ является самосопряженным в пространстве ${{U}_{S}}$. Пусть функции $f,g \in {{U}_{S}}$. Используя вторую формулу Грина [7], получим
(8)
$\begin{gathered} \left( {Lf,g} \right) - \left( {f,Lg} \right) = \\ = \,\,\mathop \smallint \limits_{{{G}_{1}}} \left[ {g\nabla \cdot \left( {k\nabla f} \right) - f\nabla \cdot \left( {k\nabla g} \right)} \right]dx - \\ - \mathop \smallint \limits_{{{G}_{2}}} \left[ {g\nabla \cdot \left( {k\nabla f} \right) - f\nabla \cdot \left( {k\nabla g} \right)} \right]dx = \\ = \,\,\mathop \smallint \limits_{{{S}^{ - }}} k\left( {g\nabla f - g\nabla f} \right) \cdot {{{\mathbf{n}}}^{ - }}dS + \\ + \,\,\mathop \smallint \limits_{{{S}^{ + }}} k\left( {g\nabla f - g\nabla f} \right) \cdot {{{\mathbf{n}}}^{ + }}dS + \\ + \,\,\mathop \smallint \limits_{{{Г}^{ - }}} k\left( {g\nabla f - f\nabla g} \right) \cdot {{{\mathbf{N}}}^{ - }}dS = \\ = \,\,\mathop \smallint \limits_{{{S}^{ + }}} k\left( {g\nabla f - g\nabla f} \right) \cdot {{{\mathbf{n}}}^{ - }}dS + \\ + \,\,\mathop \smallint \limits_{{{S}^{ + }}} k\left( {g\nabla f - g\nabla f} \right) \cdot {{{\mathbf{n}}}^{ + }}dS + \\ + \,\,\mathop \smallint \limits_Г \left( {g\alpha f - f\alpha g} \right) \cdot {{{\mathbf{N}}}^{ - }}dS = 0. \\ \end{gathered} $Первые два слагаемых в (8) дают в сумме ноль, поскольку n+ = –n– и выполняется условие (6). Третье слагаемое было получено с использованием (7) и также равно нулю. Следовательно, $\left( {Lf,g} \right) = \left( {f,Lg} \right)$, и утверждение доказано.
Обезразмеривается уравнение (1), его граничное (3) и начальное (4) условия, а также условие сопряжения (2) и вводится безразмерная координата $r{\kern 1pt} ' = \frac{r}{{{{r}_{2}}}}$. Кроме того, разделив уравнение (1) на коэффициент теплопроводности внутреннего цилиндра ${{k}_{1}}$, можно его привести, опуская штрихи, к виду
(9)
$\begin{gathered} \frac{{\rho c}}{{{{\rho }_{1}}{{c}_{1}}}}\frac{{\partial T}}{{\partial \left( {\frac{{{{k}_{1}}t}}{{r_{2}^{2}{{\rho }_{1}}{{c}_{1}}}}} \right)}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {\frac{k}{{{{k}_{1}}}}r\frac{{\partial T}}{{\partial r}}} \right) + \\ + \,\,q\left( {r,t} \right)\frac{{r_{2}^{2}}}{{{{k}_{1}}}},\,\,\,\,0 \leqslant r \leqslant 1. \\ \end{gathered} $С учетом (9) вводятся обозначения
(10)
$\begin{gathered} \xi = \frac{{{{r}_{1}}}}{{{{r}_{2}}}},\,\,\,\,p = \frac{{{{\rho }_{2}}{{c}_{2}}}}{{{{\rho }_{1}}{{c}_{1}}}},\,\,\,\,k = \frac{{{{k}_{2}}}}{{{{k}_{1}}}},\,\,\,\,a_{1}^{2} = \frac{{{{k}_{1}}}}{{{{\rho }_{1}}{{c}_{1}}}}, \\ a_{2}^{2} = \frac{{{{k}_{2}}}}{{{{\rho }_{2}}{{c}_{2}}}},\,\,\,\,a = \frac{{{{a}_{1}}}}{{{{a}_{2}}}},\,\,\,\,Q\left( {r,t} \right) = q\left( {r,t} \right)\frac{{r_{2}^{2}}}{{{{k}_{1}}}}. \\ \end{gathered} $Параметры $k,p,a$ не являются независимыми, поскольку ${{a}^{2}} = \frac{p}{k}$. Кусочно-постоянные функции определяются как
(11)
$P = \left\{ {\begin{array}{*{20}{c}} {1,r < \xi ,} \\ {p,\xi \leqslant r \leqslant 1,} \end{array}} \right.\,\,\,\,K = \left\{ {\begin{array}{*{20}{c}} {1,r < \xi ,} \\ {k,\xi \leqslant r \leqslant 1,} \end{array}} \right.$Для однородной среды значение ${\text{Fo}}$ позволяет оценить достижение процессом нагрева или охлаждения стационарного состояния или регулярного режима. В рассматриваемом случае параметр ${\text{Fo}}$, определенный формулой (12), позволяет приближенно судить о достижении этих состояний, если ${{r}_{2}} - {{r}_{1}} \ll {{r}_{2}}$ или $\left| {{{a}_{2}} - {{a}_{1}}} \right| \ll {{a}_{2}}$. Такое предположение выполняется, например, для тепловыделяющих элементов ядерного реактора или для бидонов с остеклованными радиоактивными отходами.
Безразмерный коэффициент теплоотдачи – число Био определяется по формуле
Для удобства записи граничного условия вводится величина отклонения от температуры окружающей среды
Уравнение теплопроводности, граничное условие и условие сопряжения с учетом формул (9)–(14) в новых обозначениях записываются в виде
(15)
$P\frac{{\partial U}}{{\partial {\text{Fo}}}} = \frac{1}{r}\frac{\partial }{{\partial r}}\left( {Kr\frac{{\partial U}}{{\partial r}}} \right) + Q\left( {r,t} \right),$(16)
$ - {{\left. {\frac{{\partial U}}{{\partial r}}} \right|}_{{r = 1 - 0}}} = {\text{Bi}}{{\left. U \right|}_{{r = 1}}},$(17)
$ - {{\left. {\frac{{\partial U}}{{\partial r}}} \right|}_{{r = \xi - 0}}} = - k{{\left. {\frac{{\partial U}}{{\partial r}}} \right|}_{{r = \xi + 0}}}.$Начальное условие (4) принимает вид
Выполнение условий (16) и (17) ограничивает поиск решений уравнения (15) элементами ранее введенного пространства US.
Определим собственные функции с весом P оператора теплопроводности в уравнении (15)
(19)
$ - \frac{1}{r}\frac{\partial }{{\partial r}}\left( {Kr\frac{{\partial {{U}_{n}}}}{{\partial r}}} \right) = \lambda _{n}^{2}P{{U}_{n}}.$Формальным решением уравнения (19) в случае, если K и P – константы, является функция
(20)
${{U}_{n}}\left( r \right) = b{{J}_{0}}\left( {\sqrt {\frac{P}{K}} r{{\lambda }_{n}}} \right) + c{{Y}_{0}}\left( {\sqrt {\frac{P}{K}} r{{\lambda }_{n}}} \right),$Значения коэффициентов $b$ и $c$ различны на отрезках $\left( {0,\xi } \right)$ и $\left( {\xi ,1} \right)$. Подстановка в (20) значений $P = 1$ и $K = 1$, соответствующих отрезку $r < \xi $, с учетом ограниченности решения в центре приводит к необходимости положить $c = 0$. Кроме того, следует положить $b = 1$, поскольку полученные собственные функции будут в дальнейшем нормироваться. Тогда
(21)
${{U}_{n}}\left( r \right) = {{J}_{0}}\left( {r{{\lambda }_{n}}} \right),\,\,\,\,r \leqslant \xi .$На отрезке $\xi \leqslant r \leqslant 1$ из (10) имеем $\sqrt {\frac{P}{K}} = a$. Поэтому
(22)
${{U}_{n}}\left( r \right) = {{b}_{n}}{{J}_{0}}\left( {ar{{\lambda }_{n}}} \right) + {{c}_{n}}{{Y}_{0}}\left( {ar{{\lambda }_{n}}} \right),\,\,\,\,\xi \leqslant r \leqslant 1.$Собственные функции должны удовлетворять условию сопряжения (17) и быть непрерывными в точке $\xi $, поэтому коэффициенты ${{b}_{n}}$ и ${{с}_{n}}$ в (22) определятся из системы
(23)
$\begin{gathered} {{b}_{n}}{{J}_{0}}\left( {a{{\lambda }_{n}}\xi } \right) + {{c}_{n}}{{Y}_{0}}\left( {a{{\lambda }_{n}}\xi } \right) = {{J}_{0}}\left( {{{\lambda }_{n}}\xi } \right), \\ ak\left[ {{{b}_{n}}{{J}_{1}}\left( {a{{\lambda }_{n}}\xi } \right) + {{c}_{n}}{{Y}_{1}}\left( {a{{\lambda }_{n}}\xi } \right)} \right] = {{J}_{1}}\left( {{{\lambda }_{n}}\xi } \right). \\ \end{gathered} $Решая систему (23), находим
(24)
$\begin{gathered} {{b}_{n}} = \frac{1}{{2k}}\left( {\pi {{\lambda }_{n}}\xi \left[ {{{J}_{1}}\left( {{{\lambda }_{n}}\xi } \right){{Y}_{0}}\left( {a{{\lambda }_{n}}\xi } \right) - } \right.} \right. \\ - \,\,\left. {\left. {ak{{J}_{0}}\left( {{{\lambda }_{n}}\xi } \right){{Y}_{1}}\left( {a{{\lambda }_{n}}\xi } \right)} \right]} \right), \\ {{c}_{n}} = \frac{1}{{2k}}\left( {\pi {{\lambda }_{n}}\xi \left[ {ak{{J}_{0}}\left( {{{\lambda }_{n}}\xi } \right){{J}_{1}}\left( {a{{\lambda }_{n}}\xi } \right) - } \right.} \right. \\ - \,\,\left. {\left. {{{J}_{1}}\left( {{{\lambda }_{n}}\xi } \right){{J}_{0}}\left( {a{{\lambda }_{n}}\lambda \xi } \right)} \right]} \right). \\ \end{gathered} $Далее в пространстве US вводится скалярное произведение с весом P
(25)
${{\left( {f,g} \right)}_{P}} = \mathop \smallint \limits_0^1 Pfgdr = \mathop \smallint \limits_0^\xi rfgdr + \mathop \smallint \limits_\xi ^1 rpfgdr$(26)
$\begin{gathered} \left\| {{{U}_{n}}} \right\|_{P}^{2} = \mathop \smallint \limits_0^\xi rU_{n}^{2}\left( r \right)dr + \mathop \smallint \limits_\xi ^1 prU_{n}^{2}\left( r \right)dr = \\ = \,\,\frac{1}{2}{{\xi }^{2}}\left[ {J_{0}^{2}\left( {{{\lambda }_{n}}\xi } \right) + J_{1}^{2}\left( {{{\lambda }_{n}}\xi } \right)} \right] + \frac{p}{2}\left\{ {b_{n}^{2}J_{0}^{2}\left( {a{{\lambda }_{n}}} \right) + } \right. \\ + \,\,b_{n}^{2}J_{1}^{2}\left( {a{{\lambda }_{n}}} \right) + 2{{b}_{n}}{{c}_{n}}{{J}_{0}}\left( {a{{\lambda }_{n}}} \right){{Y}_{0}}\left( {a{{\lambda }_{n}}} \right) + \\ + \,\,2{{b}_{n}}{{c}_{n}}{{J}_{1}}\left( {a{{\lambda }_{n}}} \right){{Y}_{1}}\left( {a{{\lambda }_{n}}} \right) + c_{n}^{2}\left[ {Y_{0}^{2}\left( {a{{\lambda }_{n}}} \right) + Y_{1}^{2}\left( {a{{\lambda }_{n}}} \right)} \right] - \\ - \,\,{{\xi }^{2}}\left[ {b_{n}^{2}J_{0}^{2}\left( {a{{\lambda }_{n}}\xi } \right) + 2{{b}_{n}}{{c}_{n}}{{J}_{0}}\left( {a{{\lambda }_{n}}\xi } \right){{Y}_{0}}\left( {a{{\lambda }_{n}}\xi } \right) + } \right. \\ \left. {\left. { + {{{\left[ {{{b}_{n}}{{J}_{1}}\left( {a{{\lambda }_{n}}\xi } \right) + {{c}_{n}}{{Y}_{1}}\left( {a{{\lambda }_{n}}\xi } \right)} \right]}}^{2}} + c_{n}^{2}Y_{0}^{2}\left( {a{{\lambda }_{n}}\xi } \right)} \right]} \right\}. \\ \end{gathered} $Параметры$p,a,{{b}_{n}},{{c}_{n}},\xi $ определены ранее. Собственные значения ${{\lambda }_{n}}$ определены ниже. Ортонормированные собственные функции определяются по выражению
(27)
${{u}_{n}}\left( r \right) = \frac{{{{U}_{n}}\left( r \right)}}{{{\text{||}}{{U}_{n}}{\text{|}}{{{\text{|}}}_{P}}}}.$Числитель в формуле (27) задается равенствами (21), (22), а знаменатель формулой (26).
Очевидно, справедливо используемое ниже равенство
(28)
$\left\| {{{u}_{n}}} \right\|_{P}^{2} = \mathop \smallint \limits_0^1 Pru_{n}^{2}\left( r \right)dr = 1.$Вычислим собственные значения оператора теплопроводности. Собственные функции должны удовлетворять граничному условию (16). После подстановки выражения для собственной функции (22) с коэффициентами ${{a}_{n}},{{b}_{n}}$, определенными в (24), в уравнение (16), получается характеристическое уравнение для нахождения собственных значений ${{\lambda }_{n}}$:
(29)
$\begin{gathered} X\left( \lambda \right) = a{{J}_{0}}\left( {\lambda \xi } \right)\left\{ {{{J}_{1}}\left( {a\lambda \xi } \right)\left[ {{\text{Bi}}{{Y}_{0}}\left( {a\lambda } \right) - a\lambda {{Y}_{1}}\left( {a\lambda } \right)} \right] + } \right. \\ + \,\,\left. {{{Y}_{1}}\left( {a\lambda \xi } \right)\left[ {a\lambda {{J}_{1}}\left( {a\lambda } \right) - {\text{Bi}}{{J}_{0}}\left( {a\lambda } \right)} \right]} \right\} + \\ + {{J}_{1}}\left( {\lambda \xi } \right){{Y}_{0}}\left( {a\lambda \xi } \right)\left[ {{\text{Bi}}{{J}_{0}}\left( {a\lambda } \right) - a\lambda {{J}_{1}}\left( {a\lambda } \right)} \right] + \\ + \,\,{{J}_{1}}\left( {\lambda \xi } \right){{J}_{0}}\left( {a\lambda \xi } \right)\left[ {a\lambda {{Y}_{1}}\left( {a\lambda } \right) - {\text{Bi}}{{Y}_{0}}\left( {a\lambda } \right)} \right] = 0. \\ \end{gathered} $Функция $X\left( \lambda \right)$, определенная в (29), зависит от радиусов слоев, их теплофизических свойств и условия теплообмена на границе. Корни уравнения (29) могут быть найдены численно.
Далее уравнение (15) решается методом разделения переменных. После нахождения ортонормированных собственных функций (27) дальнейший ход решения совпадает со стандартной процедурой метода разделения переменных (см. [7]).
Решение уравнения (15) представляется в виде разложения в ряд Фурье по системе ортонормированных собственных функций с коэффициентами, зависящими от числа Фурье
(30)
$U\left( {r,{\text{Fo}}} \right) = \mathop \sum \limits_{n = 1}^\infty {{c}_{n}}\left( {{\text{Fo}}} \right){{u}_{n}}\left( r \right).$Подставляя (30) в (15) и учитывая соотношение (19), получим
(31)
$P\left( {\mathop \sum \limits_{n = 1}^\infty \frac{{\partial {{c}_{n}}}}{{\partial {\text{Fo}}}} + \mathop \sum \limits_{n = 1}^\infty \lambda _{n}^{2}{{c}_{n}}} \right){{u}_{n}}\left( r \right) = Q.$Умножая уравнение (31) на $r{{u}_{n}}\left( r \right)$ и интегрируя на отрезке (0, 1), получаем
(32)
$\frac{{\partial {{c}_{n}}}}{{\partial {\text{Fo}}}} + \lambda _{n}^{2}{{c}_{n}} = {{q}_{n}},\,\,\,\,{{q}_{n}} = \left( {Q,{{u}_{n}}} \right) = \mathop \smallint \limits_0^1 rQ{{u}_{n}}dr.$При получении уравнения (32) использовалoсь соотношение (28).
Для получения начального условия для уравнения (32) следует умножить равенство (18) на $r{{u}_{n}}\left( r \right)$ и проинтегрировать на отрезке (0, 1):
(33)
${{c}_{n}}\left( 0 \right) = \mathop \smallint \limits_0^1 r{{U}_{0}}\left( r \right){{u}_{n}}dr = c_{n}^{0}.$Следует отметить, что при вычислении интегралов (32) и (33) необходимо разбить отрезок интегрирования на два – $\left( {0,\xi } \right)$ и $\left( {\xi ,1} \right)$, поскольку собственные функции ${{u}_{n}}\left( r \right)$ имеют на каждом из них свое представление.
Пусть $Q$ не зависит от времени. Тогда решением уравнения (32) с начальным условием (33) будет функция
(34)
${{c}_{n}}\left( {{\text{Fo}}} \right) = \left( {c_{n}^{0} - \frac{{{{q}_{n}}}}{{\lambda _{n}^{2}}}} \right){\text{exp}}\left( { - \lambda _{n}^{2}{\text{Fo}}} \right) + \frac{{{{q}_{n}}}}{{\lambda _{n}^{2}}}.$Для часто встречающегося на практике случая, когда начальная температура и мощность источника тепла не зависят от координаты на внутреннем цилиндре, внешнем или сразу на обоих, полезно вычислить явно интегралы $\int_0^1 {r{{u}_{n}}dr} $, содержащиеся в (32) и (33), разбив область интегрирования на два интервала:
(35)
$\begin{gathered} \mathop \smallint \limits_0^1 r{{u}_{n}}dr = {{s}_{{1n}}}\left( \xi \right) + {{s}_{{2n}}}\left( 1 \right), \\ {{s}_{{1n}}}\left( r \right) = \mathop \smallint \limits_0^r r{{u}_{n}}dr = \frac{{r{{J}_{1}}\left( {r\lambda } \right)}}{{{{\lambda }_{n}}{{{\left\| {{{U}_{n}}} \right\|}}_{p}}}},\,\,\,\,0 \leqslant r \leqslant \xi , \\ {{s}_{{2n}}}\left( r \right) = \mathop \smallint \limits_\xi ^r r{{u}_{n}}dr = \frac{\xi }{{2ak{{{\left\| {{{U}_{n}}} \right\|}}_{p}}}} \times \\ \begin{array}{*{20}{c}} \begin{gathered} \times \,\,\,\left\{ {^{{^{{^{\,}}}}}\pi akr{{J}_{0}}\left( {{{\lambda }_{n}}\xi } \right)\left[ {{{J}_{1}}\left( {a{{\lambda }_{n}}\xi } \right){{Y}_{1}}\left( {ar{{\lambda }_{n}}} \right) - } \right.} \right. \\ \left. { - \,\,{{Y}_{1}}\left( {a{{\lambda }_{n}}\xi } \right){{J}_{1}}\left( {ar{{\lambda }_{n}}} \right)} \right] + \\ \end{gathered} \\ \begin{gathered} + \,\,{{J}_{1}}\left( {{{\lambda }_{n}}\xi } \right)\left[ { - \frac{2}{{a{{\lambda }_{n}}}} + \pi r{{Y}_{0}}\left( {a{{\lambda }_{n}}\xi } \right){{J}_{1}}\left( {ar{{\lambda }_{n}}} \right) - } \right. \hfill \\ \left. {\left. { - {{\,}^{{^{{^{\,}}}}}}\pi r{{J}_{0}}\left( {a{{\lambda }_{n}}\xi } \right){{Y}_{1}}\left( {ar{{\lambda }_{n}}} \right)} \right]} \right\}, \hfill \\ \end{gathered} \end{array} \\ \xi \leqslant r \leqslant 1. \\ \end{gathered} $Величины ${{\left\| {{{U}_{n}}} \right\|}_{p}}$ вычислены в (26).
Описанный аппарат может быть расширен и на случай многослойных сред, отличающийся от рассмотренного бóльшим числом уравнений вида (15). Но сложность выражения для коэффициентов (24), а также характеристического уравнения при этом будет возрастать многократно с каждым новым слоем, что неудобно в практическом использовании.
ТЕРМОМЕХАНИЧЕСКИЕ НАПРЯЖЕНИЯ В ДВУХСЛОЙНОМ ЦИЛИНДРЕ
На основании полученного решения уравнения переноса тепла далее рассмотрены механические напряжения, вызванные тепловым расширением в двухслойном цилиндре. При изучении термомеханических напряжений в однородном цилиндре традиционно рассматриваются два случая в зависимости от граничных условий на его торцах. Первый случай – это приближение бесконечно тонкого диска, соответствующее реальной задаче с незакрепленными торцами цилиндра. В другом варианте цилиндр считается бесконечно длинным, что соответствует задаче о цилиндре с закрепленными в осевом направлении торцами. Поскольку два этих разных случая важны для технических приложений, для двухслойного цилиндра результаты получены для каждого отдельно с использованием полученного ранее температурного поля (30).
Первый случай – случай тонкого диска. При изменяющейся только по радиусу температуре тонкого диска только две компоненты тензора деформаций и две тензора напряжений являются ненулевыми (см. [11]). Приняв в качестве температурного фактора величину ${{U}_{\sigma }} = T - {{T}_{0}} = U - {{U}_{0}}$, записать их можно так
(36)
$\begin{gathered} {{\varepsilon }_{r}} = \frac{{\partial {{u}_{r}}}}{{\partial r}},\,\,\,\,{{\varepsilon }_{\phi }} = \frac{{{{u}_{r}}}}{r},\,\,\,\,{{\varepsilon }_{r}} - {{\alpha }_{T}}{{U}_{\sigma }} = \frac{1}{E}\left( {{{\sigma }_{r}} - \nu {{\sigma }_{\phi }}} \right), \\ {{\varepsilon }_{\phi }} - {{\alpha }_{T}}{{U}_{\sigma }} = \frac{1}{E}\left( {{{\sigma }_{\phi }} - \nu {{\sigma }_{r}}} \right), \\ \end{gathered} $(37)
${{U}_{\sigma }}\left( {r,t} \right) = U\left( {\frac{r}{{{{r}_{2}}}},\frac{{a_{1}^{2}t}}{{r_{2}^{2}}}} \right) - {{U}_{0}}\left( {\frac{r}{{{{r}_{2}}}}} \right).$Коэффициенты $E,\nu ,\alpha $ для внутреннего цилиндра указываются ниже с индексом “1”, для внешнего с индексом “2”. Из системы уравнений (36) находим компоненты тензора напряжений
(38)
$\begin{gathered} {{\sigma }_{{rr}}} = \frac{E}{{1 - {{\nu }^{2}}}}\left[ {{{\epsilon }_{r}} + \nu {{\sigma }_{\phi }} - \left( {1 + \nu } \right){{\alpha }_{T}}{{U}_{\sigma }}} \right] = \\ = \,\,\frac{E}{{1 - {{\nu }^{2}}}}\left[ {\frac{{\partial {{u}_{r}}}}{{\partial r}} + \nu \frac{{{{u}_{r}}}}{r} - \left( {1 + \nu } \right){{\alpha }_{T}}{{U}_{\sigma }}} \right], \\ {{\sigma }_{{\phi \phi }}} = \frac{E}{{1 - {{\nu }^{2}}}}\left[ {{{\epsilon }_{\phi }} + \nu {{\sigma }_{r}} - \left( {1 + \nu } \right){{\alpha }_{T}}{{U}_{\sigma }}} \right] = \\ = \,\,\frac{E}{{1 - {{\nu }^{2}}}}\left[ {\frac{{{{u}_{r}}}}{r} + \nu \frac{{\partial {{u}_{r}}}}{{\partial r}} - \left( {1 + \nu } \right){{\alpha }_{T}}{{U}_{\sigma }}} \right]. \\ \end{gathered} $Подставляем выражения (38) в уравнение равновесия
(39)
$\frac{{\partial {{\sigma }_{{rr}}}}}{{\partial r}} + \frac{{{{\sigma }_{{rr}}} - {{\sigma }_{{\phi \phi }}}}}{r} = 0,$(40)
$\frac{\partial }{{\partial r}}\left[ {\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r{{u}_{r}}} \right)} \right] = \left( {1 + \nu } \right){{\alpha }_{T}}\frac{{\partial {{U}_{\sigma }}}}{{\partial r}}.$Тогда с учетом (40) для каждого из слоев цилиндра получаем
(41)
$\frac{\partial }{{\partial r}}\left[ {\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r{{u}_{r}}} \right)} \right] = \left( {1 + {{\nu }_{1}}} \right){{\alpha }_{1}}\frac{{\partial {{U}_{\sigma }}}}{{\partial r}},\,\,\,\,r \leqslant {{r}_{1}};$(42)
$\frac{\partial }{{\partial r}}\left[ {\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r{{u}_{r}}} \right)} \right] = \left( {1 + {{\nu }_{2}}} \right){{\alpha }_{2}}\frac{{\partial {{U}_{\sigma }}}}{{\partial r}},\,\,\,\,{{r}_{1}} \leqslant r \leqslant {{r}_{2}}.$Параметры ${{\alpha }_{1}}$ и ${{\alpha }_{2}}$ являются значениями ${{\alpha }_{T}}$ в каждом из слоев. Интегрируя (41), получаем
(43)
$\begin{gathered} {{u}_{r}} = \frac{{\left( {1 + {{\nu }_{1}}} \right){{\alpha }_{1}}}}{r}\mathop \smallint \limits_0^r r{{U}_{\sigma }}\left( {\frac{r}{{{{r}_{2}}}},t} \right)dr + \frac{{{{C}_{1}}}}{2}r = \\ = \,\,\frac{{\left( {1 + {{\nu }_{1}}} \right){{\alpha }_{1}}}}{r}{{I}_{1}}\left( {r,t} \right) + \frac{{{{C}_{1}}}}{2}r + \frac{{{{C}_{2}}}}{r}, \\ \end{gathered} $(44)
${{I}_{1}}\left( {r,t} \right) = \mathop \smallint \limits_0^r r{{U}_{\sigma }}\left( {r,t} \right)dr.$Константы интегрирования в (43) зависят, очевидно, от времени ${\text{Fo}}$ (12), поскольку от него зависит ${{U}_{\sigma }}$. Учитывая формулы (30) и (35), а также полагая для простоты начальную температуру константой, получаем
(45)
$\begin{gathered} {{I}_{1}}\left( {r,t} \right) = r_{2}^{2}\mathop \smallint \limits_0^{\frac{r}{{{{r}_{2}}}}} \frac{r}{{{{r}_{2}}}}\left[ {{{U}_{\sigma }}\left( {r,t} \right) - {{U}_{0}}} \right]d\frac{r}{{{{r}_{2}}}} = \\ = \,\,r_{2}^{2}\mathop \sum \limits_{n = 1}^\infty {{c}_{n}}\left( {{\text{Fo}}} \right){{s}_{{1n}}}\left( {\frac{r}{{{{r}_{2}}}}} \right) - \frac{{{{r}^{2}}}}{2}{{U}_{0}}. \\ \end{gathered} $Поскольку в центре цилиндра ${{\left. u \right|}_{{r = 0}}} = 0$, то ${{C}_{2}} = 0$, и при $0 \leqslant r \leqslant {{r}_{1}}$ получаем
(46)
$\begin{gathered} {{u}_{r}} = \frac{{\left( {1 + {{\nu }_{1}}} \right){{\alpha }_{1}}}}{r}{{I}_{1}}\left( {r,t} \right) + \frac{{{{C}_{1}}}}{2}r, \\ \frac{{\partial {{u}_{r}}}}{{\partial r}} = \left( {1 + {{\nu }_{1}}} \right){{\alpha }_{1}}{{U}_{\sigma }}\left( {r,t} \right) - \frac{{\left( {1 + {{\nu }_{1}}} \right){{\alpha }_{1}}}}{{{{r}^{2}}}}{{I}_{1}}\left( {r,t} \right) + \frac{{{{C}_{1}}}}{2}. \\ \end{gathered} $Аналогично при ${{r}_{1}} \leqslant r \leqslant {{r}_{2}}$ из (42) получаем
(47)
$\begin{gathered} {{u}_{r}} = \frac{{\left( {1 + {{\nu }_{2}}} \right){{\alpha }_{2}}}}{r}{{I}_{2}}\left( {r,t} \right) + \frac{{{{C}_{2}}}}{2}r + \frac{{{{C}_{3}}}}{r}, \\ \frac{{\partial {{u}_{r}}}}{{\partial r}} = \left( {1 + {{\nu }_{2}}} \right){{\alpha }_{2}}{{U}_{\sigma }}\left( {r,t} \right) - \\ - \,\,\frac{{\left( {1 + {{\nu }_{2}}} \right){{\alpha }_{2}}}}{{{{r}^{2}}}}{{I}_{2}}\left( {r,t} \right) + \frac{{{{C}_{2}}}}{2} - \frac{{{{C}_{3}}}}{{{{r}^{2}}}}. \\ \end{gathered} $Соответственно, для ${{I}_{2}}\left( {r,t} \right)$ в (47) получаем
(48)
$\begin{gathered} {{I}_{2}}\left( {r,t} \right) = \mathop \smallint \limits_{{{r}_{1}}}^r r{{U}_{\sigma }}\left( {r,t} \right)dr = \\ = \,\,r_{2}^{2}\mathop \smallint \limits_0^{\frac{r}{{{{r}_{2}}}}} \frac{r}{{{{r}_{2}}}}\left[ {U\left( {\frac{r}{{{{r}_{2}}}},t} \right) - {{U}_{0}}\left( {\frac{r}{{{{r}_{2}}}}} \right)} \right]d\frac{r}{{{{r}_{2}}}} = \\ = \,\,r_{2}^{2}\mathop \sum \limits_{n = 1}^\infty {{c}_{n}}\left( t \right){{s}_{{2n}}}\left( {\frac{r}{{{{r}_{2}}}}} \right) - \frac{1}{2}\left( {{{r}^{2}} - r_{1}^{2}} \right){{U}_{0}}. \\ \end{gathered} $Подставляя (46) и (47) в (38), получаем напряжения для внутреннего и внешнего цилиндров:
(49)
$\begin{gathered} {{\sigma }_{{1r}}} = \frac{{ - {{\alpha }_{1}}{{E}_{1}}}}{{{{r}^{2}}}}{{I}_{1}}\left( {r,t} \right) + \frac{{{{C}_{1}}{{E}_{1}}}}{{2\left( {1 - {{\nu }_{1}}} \right)}}, \\ {{\sigma }_{{1\phi }}} = - {{\alpha }_{1}}{{E}_{1}}U + \frac{{{{\alpha }_{1}}{{E}_{1}}}}{{{{r}^{2}}}}{{I}_{1}}\left( {r,t} \right) + \frac{{{{C}_{1}}{{E}_{1}}}}{{2\left( {1 - \nu } \right)}}; \\ \end{gathered} $(50)
$\begin{gathered} {{\sigma }_{{2r}}} = \frac{{ - {{\alpha }_{2}}{{E}_{2}}}}{{{{r}^{2}}}}{{I}_{2}}\left( {r,t} \right) + {{E}_{2}}\left[ {\frac{{{{C}_{2}}}}{{2\left( {1 - {{\nu }_{2}}} \right)}} + \frac{{{{C}_{3}}}}{{\left( {1 + {{\nu }_{2}}} \right){{r}^{2}}}}} \right], \\ {{\sigma }_{{2\phi }}} = - {{\alpha }_{2}}{{E}_{2}}{{U}_{\sigma }} + \frac{{{{\alpha }_{2}}{{E}_{2}}}}{{{{r}^{2}}}}{{I}_{2}}\left( {r,t} \right) + \\ + \,\,{{E}_{2}}\left[ {\frac{{{{C}_{2}}}}{{2\left( {1 - {{\nu }_{2}}} \right)}} + \frac{{{{C}_{3}}}}{{\left( {1 + {{\nu }_{2}}} \right){{r}^{2}}}}} \right]. \\ \end{gathered} $Константы ${{C}_{1}}{\kern 1pt} --{\kern 1pt} {{C}_{3}}$ в формулах (49), (50) определяются из трех уравнений, выражающих равенство перемещений и напряжений на границе слоев, а также заданное давление среды ${{P}_{L}}$ на внешней границе второго цилиндра, равное радиальному напряжению. Эти три условия образуют систему
(51)
$\begin{gathered} {{\left( {{{u}_{r}}} \right)}_{{r = {{r}_{1}} + 0}}} = {{\left( {{{u}_{r}}} \right)}_{{r = {{r}_{1}} - 0}}},\,\,\,\,{{\left( {{{\sigma }_{r}}} \right)}_{{r = {{r}_{1}} + 0}}} = {{\left( {{{\sigma }_{r}}} \right)}_{{r = {{r}_{1}} - 0}}}, \\ {{\left( {{{\sigma }_{r}}} \right)}_{{r = {{r}_{2}}}}} = {{P}_{L}}. \\ \end{gathered} $Подставляя в (51) соотношения (46), (47), (49), (50) и учитывая, что ${{I}_{2}}\left( {{{r}_{1}},{\text{Fo}}} \right) = 0$, получаем систему уравнений для определения констант ${{C}_{1}}{\text{--}}{{C}_{3}}$:
(52)
$\begin{gathered} \frac{{\left( {1 + {{\nu }_{1}}} \right){{\alpha }_{1}}}}{{{{r}_{1}}}}{{I}_{1}}\left( {{{r}_{1}}} \right) + \frac{{{{C}_{1}}}}{2}{{r}_{1}} = \frac{{{{C}_{2}}}}{2}{{r}_{1}} + \frac{{{{C}_{3}}}}{{{{r}_{1}}}}, \\ \frac{{ - {{\alpha }_{1}}{{E}_{1}}}}{{{{r}_{1}}^{2}}}{{I}_{1}}\left( {{{r}_{1}}} \right) + \frac{{{{C}_{1}}{{E}_{1}}}}{{2\left( {1 - {{\nu }_{1}}} \right)}} = \frac{{{{E}_{2}}}}{{1 - {{\nu }_{2}}^{2}}} \times \\ \times \,\,\left( {\frac{{{{C}_{2}}\left( {1 + {{\nu }_{2}}} \right)}}{2} + \frac{{{{C}_{3}}}}{{{{r}_{1}}^{2}}}\left( {{{\nu }_{2}} - 1} \right)} \right)\frac{{ - {{\alpha }_{2}}{{E}_{2}}}}{{{{r}_{2}}^{2}}}{{I}_{2}}\left( {{{r}_{2}}} \right) + \\ + \,\,\frac{{{{E}_{2}}}}{{1 - {{\nu }_{2}}^{2}}}\left( {\frac{{{{C}_{2}}\left( {1 + {{\nu }_{2}}} \right)}}{2} + \frac{{{{C}_{3}}}}{{{{r}_{2}}^{2}}}\left( {{{\nu }_{2}} - 1} \right)} \right) = {{P}_{L}}. \\ \end{gathered} $Решая систему (52), находим значения ${{C}_{1}}{\text{--}}{{C}_{3}}$:
(53)
$\begin{gathered} {{C}_{1}} = 2\left( {1 - {{\nu }_{1}}} \right){{d}_{1}},\,\,\,\,{{d}_{1}} = \frac{1}{{r_{1}^{2}{{z}_{0}}}} \times \\ \times \,\,\left( {{{\alpha }_{1}}{{I}_{{r1}}}\left[ {{{E}_{1}}\left( {1 + {{\nu }_{2}}} \right)r_{2}^{2} + {{E}_{1}}\left( {1 - {{\nu }_{2}}} \right)r_{1}^{2} - } \right.} \right. \\ - \left. {\left. {{{E}_{2}}\left( {1 + {{\nu }_{1}}} \right)\left( {r_{2}^{2} - r_{1}^{2}} \right)} \right]} \right) + \frac{1}{{{{z}_{0}}}}\left( {2{{\alpha }_{2}}{{E}_{2}}{{I}_{{r2}}} + 2{{P}_{L}}r_{2}^{2}} \right), \\ {{z}_{0}} = {{E}_{1}}\left[ {\left( {1 + {{\nu }_{2}}} \right)r_{2}^{2} + \left( {1 - {{\nu }_{2}}} \right)r_{1}^{2}} \right] + \\ + \,\,{{E}_{2}}\left( {1 - {{\nu }_{1}}} \right)\left( {r_{2}^{2} - r_{1}^{2}} \right),\,\,\,\,{{C}_{2}} = - 2\left( {1 - {{\nu }_{2}}} \right){{d}_{2}}, \\ {{C}_{3}} = - \left( {1 + {{\nu }_{2}}} \right){{d}_{2}},\,\,\,\,{{d}_{2}} = \frac{1}{{{{z}_{0}}}} \times \\ \times \,\,\left( {\left[ {{{E}_{2}}\left( {1 - {{\nu }_{1}}} \right) - {{E}_{1}}\left( {1 - {{\nu }_{2}}} \right)} \right]{{\alpha }_{2}}r_{1}^{2}{{I}_{{r2}}} + 2{{E}_{1}}{{\alpha }_{1}}r_{2}^{2}{{I}_{{r1}}}} \right) + \\ + \,\,\frac{1}{{{{E}_{2}}{{z}_{0}}}}\left( {\left[ {{{E}_{2}}\left( {1 - {{\nu }_{1}}} \right) - {{E}_{1}}\left( {1 - {{\nu }_{2}}} \right)} \right]{{P}_{0}}r_{1}^{2}r_{2}^{2}} \right). \\ \end{gathered} $В (53) для краткости обозначено ${{I}_{{r1}}} = {{I}_{1}}\left( {{{r}_{1}},t} \right),$ ${{I}_{{r2}}} = {{I}_{2}}\left( {{{r}_{2}},t} \right)$. Теперь можно определить перемещения по формулам (46), (47) и напряжения по формулам (49), (50) для случая тонкого диска.
Далее аналогичные решения можно получить для бесконечно длинного (или зафиксированного по торцам) цилиндра. Для компонент тензора деформаций при изменяющейся только по радиусу температуре зафиксированного по торцам цилиндра ненулевыми являются только две компоненты тензора деформаций и три компоненты тензора напряжений (см. [11]):
(54)
$\begin{gathered} {{\varepsilon }_{r}} - \alpha {{U}_{\sigma }} = \frac{1}{E}\left( {{{\sigma }_{r}} - \nu \left( {{{\sigma }_{\phi }} + {{\sigma }_{z}}} \right)} \right), \\ {{\varepsilon }_{\phi }} - {{\alpha }_{T}}{{U}_{\sigma }} = \frac{1}{E}\left( {{{\sigma }_{\phi }} - \nu \left( {{{\sigma }_{r}} + {{\sigma }_{z}}} \right)} \right), \\ {{\varepsilon }_{z}} - {{\alpha }_{T}}{{U}_{\sigma }} = \frac{1}{E}\left( {{{\sigma }_{z}} - \nu \left( {{{\sigma }_{\phi }} + {{\sigma }_{r}}} \right)} \right), \\ \end{gathered} $(55)
${{\sigma }_{z}} = \nu \left( {{{\sigma }_{\phi }} + {{\sigma }_{r}}} \right) - {{\alpha }_{T}}E{{U}_{\sigma }}.$Исключение напряжений из системы (54) и подстановка в уравнение равновесия (39) после преобразований приводит к уравнению для перемещений
(56)
$\frac{\partial }{{\partial r}}\left[ {\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r{{u}_{r}}} \right)} \right] = \frac{{1 + \nu }}{{1 - \nu }}{{\alpha }_{T}}\frac{{\partial {{U}_{\sigma }}}}{{\partial r}}.$Интегрируя дважды (56) на отрезках $\left( {0,{{r}_{1}}} \right)$ и $\left( {{{r}_{1}},{{r}_{2}}} \right)$ и используя те же соображения, что и в случае тонкого диска, получаем выражения для перемещений
(57)
$\begin{gathered} {{u}_{r}} = \frac{{\left( {1 + {{\nu }_{1}}} \right){{\alpha }_{1}}}}{{\left( {1 - {{\nu }_{1}}} \right)r}}{{I}_{1}}\left( r \right) + {{C}_{1}}r,\,\,\,\,r \leqslant {{r}_{1}}; \\ {{u}_{r}} = \frac{{\left( {1 + {{\nu }_{2}}} \right){{\alpha }_{2}}}}{{\left( {1 - {{\nu }_{2}}} \right)r}}{{I}_{2}}\left( r \right) + {{C}_{2}}r + \frac{{{{C}_{3}}}}{r},\,\,\,\,{{r}_{1}} \leqslant r \leqslant {{r}_{2}}. \\ \end{gathered} $Выражения (57) позволяют получить выражения для радиальных и окружных напряжений в обоих цилиндрах:
(58)
$\begin{gathered} {{\sigma }_{{1r}}} = \frac{{ - {{\alpha }_{1}}{{E}_{1}}}}{{\left( {1 - {{\nu }_{1}}} \right){{r}^{2}}}}{{I}_{1}}\left( {r,t} \right) + \frac{{{{E}_{1}}}}{{1 + {{\nu }_{1}}}}\left( {\frac{{{{C}_{1}}}}{{1 - 2{{\nu }_{1}}}}} \right),\,\,\,\,r \leqslant {{r}_{1}}; \\ {{\sigma }_{{2r}}} = \frac{{ - {{\alpha }_{2}}{{E}_{2}}}}{{\left( {1 - {{\nu }_{2}}} \right){{r}^{2}}}}{{I}_{2}}\left( {r,t} \right) + \frac{{{{E}_{2}}}}{{1 + {{\nu }_{2}}}} \times \\ \times \,\,\left( {\frac{{{{C}_{2}}}}{{1 - 2{{\nu }_{2}}}} - \frac{{{{C}_{3}}}}{{{{r}^{2}}}}} \right),\,\,\,\,{{r}_{1}} \leqslant r \leqslant {{r}_{2}}; \\ {{\sigma }_{{1\phi }}} = \frac{{{{\alpha }_{1}}{{E}_{1}}}}{{\left( {1 - {{\nu }_{1}}} \right){{r}^{2}}}}{{I}_{1}}\left( {r,t} \right) - \frac{{{{\alpha }_{1}}{{E}_{1}}{{U}_{\sigma }}}}{{1 - {{\nu }_{1}}}} + \\ + \,\,\frac{{{{E}_{1}}}}{{1 + {{\nu }_{1}}}}\frac{{{{C}_{1}}}}{{1 - 2{{\nu }_{1}}}},\,\,\,\,r \leqslant {{r}_{1}}; \\ {{\sigma }_{{2\phi }}} = \frac{{{{\alpha }_{2}}{{E}_{2}}}}{{\left( {1 - {{\nu }_{2}}} \right){{r}^{2}}}}{{I}_{2}}\left( {r,t} \right) - \frac{{{{\alpha }_{2}}{{E}_{2}}{{U}_{\sigma }}}}{{1 - {{\nu }_{2}}}} + \\ + \,\,\frac{{{{E}_{2}}}}{{1 + {{\nu }_{2}}}}\left( {\frac{{{{C}_{2}}}}{{1 - 2{{\nu }_{2}}}} + \frac{{{{C}_{3}}}}{{{{r}^{2}}}}} \right),\,\,\,\,{{r}_{1}} \leqslant r \leqslant {{r}_{2}}. \\ \end{gathered} $Подставляя выражения (57) и (58) в условия (51), получаем систему для определения коэффициентов ${{C}_{1}}{\text{--}}{{C}_{3}}$:
(59)
$\begin{gathered} \frac{{\left( {1 + {{\nu }_{1}}} \right){{\alpha }_{1}}}}{{\left( {1 - {{\nu }_{1}}} \right){{r}_{1}}}}{{I}_{1}}\left( {{{r}_{1}},F} \right) + {{C}_{1}}{{r}_{1}} = {{C}_{2}}{{r}_{1}} + \frac{{{{C}_{3}}}}{{{{r}_{1}}}}, \\ \frac{{ - {{\alpha }_{1}}{{E}_{1}}}}{{\left( {1 - {{\nu }_{1}}} \right){{r}_{1}}^{2}}}{{I}_{1}}\left( {{{r}_{1}},F} \right) + \frac{{{{E}_{1}}}}{{1 + {{\nu }_{1}}}}\left( {\frac{{{{C}_{1}}}}{{1 - 2{{\nu }_{1}}}}} \right) = \\ = \,\,\frac{{{{E}_{2}}}}{{1 + {{\nu }_{2}}}}\left( {\frac{{{{C}_{2}}}}{{1 - 2{{\nu }_{2}}}} - \frac{{{{C}_{3}}}}{{{{r}_{1}}^{2}}}} \right), \\ \frac{{ - {{\alpha }_{2}}{{E}_{2}}}}{{\left( {1 - {{\nu }_{2}}} \right){{r}_{2}}^{2}}}{{I}_{2}}\left( {{{r}_{2}},F} \right) + \frac{{{{E}_{2}}}}{{1 + {{\nu }_{2}}}}\left( {\frac{{{{C}_{2}}}}{{1 - 2{{\nu }_{2}}}} - \frac{{{{C}_{3}}}}{{{{r}_{2}}^{2}}}} \right) = {{P}_{L}}. \\ \end{gathered} $Решение системы (59) для подстановки в (58) удобнее представить в виде
(60)
$\begin{gathered} {{C}_{1}} = \left( {2\nu _{1}^{2} + {{\nu }_{1}} - 1} \right){{d}_{1}},\,\,\,\,{{z}_{0}} = {{E}_{2}}\left( {2\nu _{1}^{2} + {{\nu }_{1}} - 1} \right) \times \\ \times \,\,\left( {r_{2}^{2} - r_{1}^{2}} \right) + {{E}_{1}}\left( {1 + {{\nu }_{2}}} \right)\left[ {\left( {1 - 2{{\nu }_{2}}} \right)r_{1}^{2} + r_{2}^{2}} \right], \\ {{d}_{1}} = \frac{{ - {{\alpha }_{1}}{{I}_{1}}}}{{\left( {1 - {{\nu }_{1}}} \right)r_{1}^{2}{{z}_{0}}}}\left( {{{E}_{1}}\left( {1 + {{\nu }_{2}}} \right)\left[ {\left( {1 - 2{{\nu }_{2}}} \right)r_{1}^{2} + r_{2}^{2}} \right] + } \right. \\ \left. { + \,\,{{E}_{2}}\left( {1 + {{\nu }_{1}}} \right)\left( {r_{1}^{2} - r_{2}^{2}} \right)} \right) - \frac{{2\left( {1 + {{\nu }_{2}}} \right)}}{{{{z}_{0}}}} \times \\ \times \,\,\left( {\left( {1 - {{\nu }_{2}}} \right){{P}_{0}}r_{2}^{2} + {{E}_{2}}{{\alpha }_{2}}{{I}_{2}}} \right),\,\,\,\,{{C}_{2}} = \left( {1 - 2{{\nu }_{2}}} \right) \times \\ \times \,\,\left( {1 + {{\nu }_{2}}} \right){{d}_{2}},\,\,\,\,{{d}_{2}} = \frac{{{{\alpha }_{2}}{{I}_{2}}}}{{\left( {1 - {{\nu }_{2}}} \right){{z}_{0}}}}\left( {\left( {2\nu _{1}^{2} + {{\nu }_{1}} - 1} \right){{E}_{2}} - } \right. \\ \left. { - \,\,\left( {1 + {{\nu }_{2}}} \right){{E}_{1}}} \right) + 2{{\alpha }_{1}}{{I}_{1}}\frac{{{{E}_{1}}\left( {1 + {{\nu }_{1}}} \right)}}{{{{z}_{0}}}} + \\ + \,\,{{P}_{0}}\left[ {\frac{{{{E}_{1}}{{{\left( {1 + {{\nu }_{2}}} \right)}}_{0}}r_{2}^{2}}}{{{{E}_{2}}{{z}_{0}}}} - r_{2}^{2}\left( {2\nu _{1}^{2} + {{\nu }_{1}} - 1} \right)} \right], \\ {{C}_{3}} = \left( {1 + {{\nu }_{2}}} \right){{d}_{3}},\,\,{{d}_{3}} = \frac{{ - {{I}_{2}}{{\alpha }_{2}}r_{1}^{2}}}{{\left( {1 - {{\nu }_{2}}} \right){{z}_{0}}}} \times \\ \times \,\,\left( {\left( {2\nu _{1}^{2} + {{\nu }_{1}} - 1} \right){{E}_{2}} - \left( {2\nu _{2}^{2} + {{\nu }_{2}} - 1} \right){{E}_{1}}} \right) + \\ + \frac{1}{{{{z}_{0}}}}\left( {2{{\alpha }_{1}}{{E}_{1}}\left( {1 + {{\nu }_{1}}} \right)r_{2}^{2}{{I}_{1}}} \right)\frac{{ - {{P}_{0}}r_{1}^{2}r_{2}^{2}}}{{{{z}_{0}}}} \times \\ \times \,\,\left[ {2\nu _{1}^{2} + {{\nu }_{1}} - 1 - \frac{{{{E}_{1}}\left( {2\nu _{2}^{2} + {{\nu }_{2}} - 1} \right)}}{{{{E}_{2}}}}} \right]. \\ \end{gathered} $Подстановка (60) в формулы (58) дает выражения для напряжений:
(61)
$\begin{gathered} {{\sigma }_{{1r}}} = \frac{{ - {{\alpha }_{1}}{{E}_{1}}}}{{\left( {1 - {{\nu }_{1}}} \right){{r}^{2}}}}{{I}_{1}}\left( {r,F} \right) - {{E}_{1}}{{d}_{1}},\,\,\,\,r \leqslant {{r}_{1}}; \\ {{\sigma }_{{2r}}} = \frac{{ - {{\alpha }_{2}}{{E}_{2}}}}{{\left( {1 - {{\nu }_{2}}} \right){{r}^{2}}}}{{I}_{2}}\left( {r,F} \right) - {{E}_{2}}\left( {{{d}_{2}} - \frac{{{{d}_{3}}}}{r}} \right), \\ {{r}_{1}} \leqslant r \leqslant {{r}_{2}}; \\ {{\sigma }_{{1\phi }}} = \frac{{{{\alpha }_{1}}{{E}_{1}}}}{{\left( {1 - {{\nu }_{1}}} \right){{r}^{2}}}}{{I}_{1}}\left( {r,F} \right) - \frac{{{{\alpha }_{1}}{{E}_{1}}{{U}_{\sigma }}}}{{1 - {{\nu }_{1}}}} - {{E}_{1}}{{d}_{1}}, \\ r \leqslant {{r}_{1}}; \\ {{\sigma }_{{2\phi }}} = \frac{{{{\alpha }_{2}}{{E}_{2}}}}{{\left( {1 - {{\nu }_{2}}} \right){{r}^{2}}}}{{I}_{2}}\left( {r,F} \right) - \frac{{{{\alpha }_{2}}{{E}_{2}}{{U}_{\sigma }}}}{{1 - {{\nu }_{2}}}} + \\ + \,\,{{E}_{2}}\left( {{{d}_{2}} - \frac{{{{d}_{3}}}}{r}} \right),\,\,\,\,{{r}_{1}} \leqslant r \leqslant {{r}_{2}}. \\ \end{gathered} $Уравнения для напряжений (61) отличаются от их первоначальной записи в виде (58) отсутствием неопределенности вида $\frac{0}{0}$ для случаев ${{\nu }_{1}} = \frac{1}{2}$ и/или ${{\nu }_{2}} = \frac{1}{2}$, которые соответствуют случаю “несжимаемости”.
ИСПОЛЬЗОВАНИЕ РЕШЕНИЙ ДЛЯ ВЕРИФИКАЦИИ РАСЧЕТНЫХ СРЕДСТВ
Как отмечалось выше, аналитические модели, как правило, менее точно отражают особенности конструкции моделируемых объектов, нежели численные. В то же время в упрощенной постановке решения будут очень близки или будут совпадать. В рамках верификации трехмерного конечноэлементного кода FENIA [12, 13], разрабатываемого как средство для сопряженного моделирования тепловых и механических процессов на различных инженерных объектах, можно сравнить аналитические решения для полей температур и механических напряжений с численными. Использование полученного аналитического решения позволяет проверить точность моделирования численным методом разрывов напряжений на границе сред, вызванных скачкообразным изменением свойств материалов. Аналитическое решение обеспечивает точное представление рассчитываемых величин вблизи разрыва, в то время как численное моделирование с высокой точностью сопряжено со значительными сложностями и может приводить к заметным ошибкам в этой области [14]. Важно отметить, что полученное решение служит для проверки лишь одного из аспектов качества расчетного кода и является малой частью матрицы верификации. Действительно, для полной верификации кода требуются трехмерные задачи. Тем не менее, данное решение хорошо соответствует области применимости кода и позволяет проверить его работу на задачах, близких к реальным.
В приближенной постановке моделируется охлаждениe бидона с остеклованными радиоактивными отходами (РАО) в специализированном хранилище ФГУП “ПО “Маяк”. Рассматривается упрощенная конфигурация хранилища, позволяющая использовать аналитическое описание. Кроме того, имеются значительные неопределенности физических свойств стекол [15, 16], радионуклидного состава и, как следствие, тепловыделения [17], а также режима охлаждения. Таким образом, приводимые результаты носят прежде всего верификационный и оценочный характер.
Специализированное хранилище состоит из бетонных отсеков (стояков), внутри которых расположены пеналы из нержавеющей стали для размещения бидонов с РАО [18]. Каждый пенал вмещает три бидона. Остеклованные отходы обладают достаточно высоким тепловыделением – до 5 кВт/м3 [19]. Охлаждение бидонов осуществляется подачей воздуха в зазор между пеналом и стояком. Схема расположения бидонов показана на рис. 3. Габаритная высота бидона составляет 1000 мм. В каждом бидоне помещается около 0.2 м3 остеклованных отходов массой ~500 кг. Внешний диаметр пенала составляет 630 мм, внутренний диаметр стояка 650 мм, т.е. ширина воздушного зазора составляет 10 мм. Скорость потока воздуха в зазоре составляет 2 м/с. Толщина стенок пенала – 8 мм, высота пенала – 3400 мм [18].
В расчете принималось, что для иммобилизации РАО применялось стекло с содержанием РАО 20 мас. %, пенал выполнен из стали 08Х18Н10. Зазор между бидоном и пеналом не моделировался, считалось, что стекло непосредственно контактирует с пеналом. Параметры расчета, приведенные в таблице, задавались на основании данных о свойствах материалов из [15, 16, 20, 21]. Температурные зависимости для них не учитывались. Тепловыделение принималось равномерным по объему бидона и равным 5 кВт/м3 (проектный предел для помещаемых в хранилище отходов [19]) в начале расчета и снижалось по экспоненциальному закону
где q0 – начальное тепловыделение, а значение параметра p таково, что тепловыделение снижается в e раз за 40 лет. Это позволяет считать, что для рассматриваемой ниже задачи, в которой механические напряжения меняются в течение одних суток, вместо (62) можно положить $q\left( t \right) = {{q}_{0}} = {\text{const}}$ и использовать формулу (34). В качестве начальной принималась температура РАО 500°С. Это связано с тем, что при более высоких температурах упругая модель неприменима вследствие размягчения стекла [17, 22]. Температура воздуха принималась равной 20°С, что соответствует температуре на входе в стояк. Коэффициент теплоотдачи в граничных условиях (3), рассчитанный по методике, приведенной в [18], взят равным 150 Вт/(м2 К). При моделировании напряженно-деформированного состояния на верхнем и нижнем торцах бидона задавалось условие нулевого перемещения вдоль оси бидона, т.е. предполагалось, что давление вышележащих бидонов препятствует вертикальным перемещениям. Начальные напряжения считались нулевыми.Свойства материалов, принятые в расчете
| Остеклованные РАО | Сталь 08Х18Н10 | |
|---|---|---|
| Коэффициент теплопроводности, Вт/(м К) | 1.5 | 18 |
| Плотность, кг/м3 | 2760 | 7800 |
| Удельная теплоемкость, Дж/(кг К) | 800 | 500 |
| Коэффициент Пуассона | 0.25 | 0.29 |
| Модуль Юнга, ГПа | 83 | 200 |
| Коэффициент линейного термического расширения, 1/К | 11.8 × 10–6 | 18 × 10–6 |
Полученные аналитически и с помощью расчетного кода FENIA температуры приведены на рис. 4. Кратковременное повышение температуры в центре бидона связано с превышением остаточного тепловыделения РАО над отводом тепла. Видно, что за время порядка двух суток тепловое состояние стабилизируется. Описанные результаты получены для нижнего бидона в стояке, для которого температура охлаждающего воздуха минимальна. Проходящий мимо расположенных выше бидонов воздух будет иметь более высокую температуру, вследствие чего для них выход на квазистационарное состояние будет более длительным, а температура более высокой.
Рис. 4.
Расчетные зависимости температуры: (а) – вдоль радиуса бидона для двух моментов времени: 1, 3 – 12 ч; 2, 4 – 45 ч; (б) – от времени в двух точках: 1, 3 – в центре; 2, 4 – на поверхности между цилиндрами; 1, 2 – аналитический расчет; 3, 4 – FENIA.
На рис. 5 приведены результаты расчета напряженно-деформированного состояния бидона. Очень высокие окружные и осевые напряжения на границе стекла и стали связаны с тем, что не моделировалась предыстория – заливка стекла в бидон с его разогревом и учетом изменения механических параметров стекла при затвердевании, а также последующее остывание бетона. При учете этих процессов расчетные напряжения будут более реалистичными. Также при моделировании не в полной мере учитывались эффекты, связанные с нагрузкой со стороны вышележащих бидонов, равно как и влияние днища и крышки бидона.
Рис. 5.
Расчетные зависимости напряженно-деформированного состояния бидона для двух моментов времени: (а) – радиальные напряжения, (б) – окружные напряжения, (в) – осевые напряжения, (г) – перемещения: 1, 3 – 12 ч; 2, 4 – 45 ч; 1, 2 – аналитический расчет; 3, 4 – FENIA.
Проведенные расчеты показали полное совпадение численного и аналитического решений, т.е. расчетный код FENIA корректно решает задачи такого рода. Сравнение показывает, что численное решение позволяет хорошо отследить “разрывы” напряжений, возникающие вследствие “разрыва” свойств материалов.
ЗАКЛЮЧЕНИЕ
Полученные аналитические решения могут использоваться как для верификации расчетных средств, так и для получения оценок температур и напряжений в различных объектах, которые могут быть приближенно представлены в виде двухслойного цилиндра. В частности, решение может использоваться для проверки точности численного моделирования разрывов напряжений на границе сред вследствие скачкообразного изменения свойств материалов.
Проведенный расчет для остывания бидона с остеклованными РАО показал, что с помощью аналитического решения можно получить достаточно надежную оценку изменения поля температур. Для моделирования же напряженно-деформированного состояния полученные результаты дают сильно завышенную оценку напряжений. Это связано не с недостатками процедуры численного или аналитического решения, а с необходимостью учета начального пространственного распределения напряжений. Его требуется рассчитывать численно, так как при остывании стекла заметно меняются его теплофизические свойства, причем при переходе через температуру стеклования его поведение меняется с вязкопластического на упругое. Кроме того, требуется учитывать неопределенности исходных данных для моделирования, которые связаны как с влиянием РАО на свойства стекла, так и с особенностями размещения бидона в стояке, от которых зависит и его охлаждение, и напряжения и деформации, вызванные давлением вышележащих бидонов.
Список литературы
Самарский А.А., Вабищевич П.Н. Вычислительная теплопередача. М.: Едиториал УРСС, 2003. 784 с.
арзов Г.П., Марголин Б.З., Швецова В.А. Физико-механическое моделирование процессов разрушения. СПб.: Политехника, 1993. 391 с
Токарев Ю.Н., Моисеенко Е.В., Дробышевский Н.И., Бутов Р.А. Точные аналитические решения нестационарных задач теплопроводности и их роль в верификации моделей тепло-гидро-механических процессов в пунктах глубинного захоронения РАО // Радиоактивные отходы. 2018. № 4(5). С. 90.
Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer. ASME, 2009. 100 p.
Roache P.J. Building PDE Codes to be Verifiable and Validatable // Comput. Sci. Eng. 2004. V. 6. Iss. 5. P. 30.
Oberkampf W.L., Roy C.J. Verification and Validation in Scientific Computing. Cambridge Univ. Press, 2010. 790 p.
Владимиров В.С. Уравнения математической физики. Изд. 4-е. М.: Наука; Глав. ред. физ.-мат. лит., 1981. 512 с.
Цветков Ф.Ф., Григорьев Б.А. Тепломассообмен. М.: МЭИ, 2001. 548 с.
Жорник В.А., Рыбинская А.А., Савочка П.А. Об одной нестационарной задаче термоупругости для двухслойного цилиндра // Изв. вузов. Северо-Кавказский регион. Естеств. науки. 2006. № S12. С. 34.
Видин Ю.В., Злобин В.С. Аналитический расчет нестационарного температурного поля плоского тела при переменном коэффициенте теплопроводности // ТВТ. 2019. Т. 57. Вып. 5. С. 790.
Тимошенко С.П., Гудьер Дж. Теория упругости. М.: Наука, 1975. 576 с.
Бутов Р.А., Дробышевский Н.И., Моисеенко Е.В., Токарев Ю.Н. Численное моделирование напряженно-деформированного состояния массива горных пород при размещении в нем пункта глубинного захоронения радиоактивных отходов и визуализация результатов // Горный информационно-аналитический бюллетень. 2019. № S37. С. 343.
Butov R.A., Drobyshevsky N.I., Moiseenko E.V., Tokarev U.N. Finite Element Code FENIA Verification and Application for 3D Modelling of Thermal State of Radioactive Waste Deep Geological Repository // J. Phys.: Conf. Ser. 2017. V. 891. 012174.
Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. 2-е изд. М.–СПб.: Физматлит; Лаб. базовых знаний, 2001. 630 с.
Алой А.С., Трофименко А.В., Кольцова Т.И., Никандрова М.В. Физико-химические характеристики остеклованных модельных ВАО ОДЦ ГХК // Радиоактивные отходы. 2018. № 4(5). С. 67.
Ремизов М.Б., Козлов П.В., Власова Н.В., Беланова Е.А., Руденко А.В., Катаев А.А., Редькин А.А., Ткачева О.Ю., Докутович В.Н., Филатов Е.С., Зайков Ю.П. Тепло- и электропроводность расплавов алюмофосфатных и боросиликатных стекол, содержащих имитаторы высокоактивных отходов от переработки ОЯТ // Физика и химия стекла. 2019. Т. 45. № 2. С. 120.
Ремизов М.Б., Козлов П.В., Борисенко В.П., Дементьева И.И., Блохин П.А., Самойлов А.А. Разработка алгоритма оценки радионуклидного состава остеклованных ВАО ФГУП “ПО “Маяк” для цели их безопасного захоронения // Радиоактивные отходы. 2018. № 3(4). С. 102.
Бельтюков В.А., Зайков Ю.В. и др. Хранение отработавших свой ресурс радиоактивных источников тепла и энергии на ПО “Маяк” // Вопросы радиационной безопасности. 2000. № 2. С. 19.
Глаголенко Ю.В., Дзекун Е.Г., Медведев Г.М., Ровный С.И., Уфимцев В.П., Захаркин Б.С. Переработка отработавшего ядерного топлива АЭС и жидких радиоактивных отходов на ПО “Маяк” // Атомная энергия. 1997. Т. 83. Вып. 6. С. 446.
ГОСТ EN 1748-1-1-2016 Стекло боросиликатное. Технические требования. Межгос. стандарт. 2018-03-01.
Кириллов П.Л., Терентьева М.И., Денискина М.Б. Теплофизические свойства материалов ядерной техники. Учеб. спр. пособ. для студ. / Под общ. ред. проф. Кириллова П.Л.; 2-е изд., испр. и доп. М.: ИздАТ, 2007. 200 с.
Nuclear Waste Conditioning / Ed. Pradel Ph. Paris: Nuclear Energy Division, 2009. 145 p.
Дополнительные материалы отсутствуют.
Инструменты
Теплофизика высоких температур
Пн.-пт.: 10.00-19.00