
\documentclass[10pt]{amsart}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{graphicx}

\def\udcs{517.95} %Здесь автор определяет УДК своей работы
\def\mscs{35R30} %Здесь автор определяет классификаторы AMS своей работы
\setcounter{page}{144}

\newtheorem{lemma}{Лемма}
\newtheorem{theorem}{Теорема}
\newtheorem{definition}{Определение}
\newtheorem{corollary}{Следствие}
\newtheorem{proposition}{Предложение}
\newtheorem{remark}{Замечание}

\def\logo{{\bf\huge S\raisebox{0.2ex}{\hspace{0.55ex}\raisebox{0.05ex}e\hspace{-1.65ex}$\bigcirc$}MR}}

\def\semrtop
     {
  \vbox{
     \noindent\logo
     \hspace{80mm}\raisebox{1ex}{ISSN 1813-3304 }

     \vspace{5mm}

     \begin{center}
     {\huge СИБИРСКИЕ \ ЭЛЕКТРОННЫЕ} \\[2mm]
     {\huge МАТЕМАТИЧЕСКИЕ ИЗВЕСТИЯ} \\[2mm]
     {\large Siberian Electronic Mathematical Reports} \\[1mm]
     {\LARGE\tt{http://semr.math.nsc.ru}}\\[0.5mm]
%     {\small Институт математики им. С.Л.Соболева СО РАН}\\[-1mm]
%     {\small Sobolev Institute of Mathematics SB RAS}
     \end{center}
     \vspace{-3mm}
     \noindent
     \begin{tabular}{c}
     \hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa} \\
     \hline\hline
     \end{tabular}

     \vspace{1mm}
     {\flushleft\it Том 16, стр. 144--144 (2019) \hspace{65mm}{\rm\small УДК \udcs}}
     \newline
     	{\rm\small DOI~10.33048/semi.2019.16.xxx}\hphantom{aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa}{\rm\small MSC\ \ \mscs }
  }%\vbox
}

\begin{document}
\renewcommand{\refname}{References}


\thispagestyle{empty}

\title[Корректность модели сложного теплообмена]{Корректность модели сложного теплообмена на основе упрощенного метода
сферических гармоник третьего порядка}
\author{{Г.В. ГРЕНКИН}}%
\address{Gleb Vladimirovich Grenkin  % Контактные данные всех авторов и место работы указываются только на английском языке
\newline\hphantom{iii} Vladivostok State University,
\newline\hphantom{iii} ul. Gogolya, 41,
\newline\hphantom{iii} 690014, Vladivostok, Russia}%
\email{glebgrenkin@gmail.com}%

\thanks{\sc Grenkin, G.V.,
Correctness of a complex heat transfer model based on the simplified spherical
harmonics method of the third order}
\thanks{\copyright \ 2023 Гренкин Г.В}
\thanks{\it Поступила 1 января 2015 г., опубликована 31 декабря 2015 г.}%


\semrtop \vspace{1cm}
\maketitle {\small
\begin{quote}
\noindent{\sc Abstract. }
The paper deals with a system of equations that describes steady-state process of
radiative-conductive heat transfer in a bounded domain with
boundary conditions of specular and diffuse reflection of radiation
and boundary conditions of the third kind for temperature.
For the description of radiative energy field, the $SP_3$ approximation
of the simplified spherical harmonics method is used.
We establish properties of existence and uniqueness of the solution of the
boundary value problem under constraints on coefficients in boundary conditions
which are fulfilled over entire range of feasible physical data.
\medskip

\noindent{\bf Keywords:} radiative heat transfer, $SP_3$ approximation, $SP_N$ approximation, cross diffusion, unique solvability.
 \end{quote}
}


\section{Введение}


Повышение точности приближенного решения уравнения переноса излучения, совмещенного с уравнением теплопроводности, представляет практический интерес в связи с техническими и медицинскими приложениями \cite{Cheb2, LaserAbl, LaserAbl2,  Schneider, Larsen2002}. Как показало сравнение \cite{Schneider, Larsen2002, ModestLei, KovtMC}, использование диффузионного приближения позволяет неплохо описать общую картину сложного теплообмена, в то время как для достижения точности, близкой к той, которую дают более трудоемкие методы Монте-Карло, потребуются приближения метода сферических гармоник хотя бы третьего порядка.
Однако, с математической точки зрения, эти уравнения требуют дополнительного анализа корректности для более надежного применения их на практике.

В данной работе исследуется стационарная модель сложного теплообмена на основе $SP_3$-приближения уравнения переноса излучения. Ранее была доказана корректность $SP_N$-приближений без учета зависимости источников излучения от температуры \cite{ZhengHan}, а также корректность нестационарной модели сложного теплообмена применительно к анализу обратных задач \cite{TsePinnau}. Стационарные модели сложного теплообмена изучались в рамках диффузионного ($P_1$) приближения \cite{KovtCheb, SolvabKovt, Inhomog}. Отметим также работы \cite{Amosov1, Amosov2}, в которых установлена корректность моделей теплообмена с полным уравнением переноса излучения, и работы \cite{SPN2, SPN1, SPN4, KloseLarsen, SPN3, SPN5, SPN6}, посвященные теоретическому и численному анализу $SP_N$-приближений.

Важно установить корректность рассматриваемой модели для обоснования возможности использования стационарных уравнений для описания установившегося состояния полей температуры и энергии излучения. Существование решения краевой задачи для $SP_3$-модели сложного теплообмена доказывается таким же методом, как и при анализе $P_1$-модели \cite{KovtCheb}. Единственность решения получена при дополнительном условии на коэффициенты в граничных условиях, которые выполняются для всего диапазона допустимых физических данных.



\section{Вывод варианта модели. Постановка задачи}

Ранее в работе \cite{TsePinnau} были исследованы нестационарные уравнения $SP_3$-приближения вместе с уравнением теплопроводности и уравнением свободной конвекции, при этом граничные условия были записаны в общем виде, и относительно коэффициентов предполагалось условие строгого диагонального преобладания. 
В этом разделе мы выведем выражения для граничных коэффициентов и сформулируем условие, которому они должны удовлетворять, чтобы краевая задача гарантированно имела единственное решение. Здесь же мы установим, что это условие выполняется для всего диапазона допустимых физических данных в рамках варианта $SP_3$-модели с краевыми условиями зеркального и диффузного отражения с коэффициентами отражения, не зависящими от направления падающих лучей.

Установившийся процесс теплового излучения в среде, пропускающей, излучающей, поглощающей и рассеивающей тепловое излучение, заполняющей область $\Omega \subset \mathbb R^3$ с границей $\Gamma$, которая излучает, поглощает, а также зеркально и диффузно отражает тепловое излучение, моделируется уравнением \cite{Modest}
$$
\omega \cdot \nabla_x I(x, \omega) + \kappa I(x, \omega)
= \frac{\kappa_s}{4\pi}\int_S P(\omega, \omega') I(x, \omega')\,d\omega' + \kappa_aI_b(T(x)) + \frac{q(x)}{4\pi}
$$
с граничным условием
$$
I(x, \omega) = \varepsilon(x)I_b(T_b(x))+ \rho_s(x)I(x, \omega_R) + \frac{\rho_d(x)}{\pi}\int\limits_{\omega'\cdot\mathbf n > 0}I(x, \omega')\omega' \cdot \mathbf n\,d\omega', \;\; \omega \cdot \mathbf n < 0, \; x \in \Gamma.
$$
Здесь $I(x,\omega)$ -- интенсивность излучения, $T(x)$ -- абсолютная температура. Направления излучения отождествляются с точками единичной сферы $S = \{\omega \in \mathbb R^3 \colon |\omega| = 1\}$. Через $I_b(T)$ обозначена интенсивность излучения абсолютно черного тела при температуре $T$: $I_b(T) = \sigma n^2 T^4/\pi$, где $n$ -- показатель преломления, $\sigma$ -- постоянная Стефана-Больцмана. Через $T_b(x)$ обозначена температура границы области, $q(x)$ -- объемная плотность изотропных источников излучения, $\kappa_a$ -- коэффициент поглощения, $\kappa_s$ -- коэффициент рассеяния, $\kappa = \kappa_a + \kappa_s$, $P(\omega,\omega')$ -- фазовая функция рассеяния, $\varepsilon$ -- коэффициент излучения поверхности, $\rho_s$ и $\rho_d$ --- коэффициенты зеркального и диффузного отражения, $\varepsilon + \rho_s + \rho_d = 1$, $\omega_R = \omega - 2(\omega \cdot \mathbf n)\mathbf n$ -- направление отражения, $\mathbf n$ -- вектор внешней нормали к $\Gamma$.

Поле температуры описывается уравнением теплопроводности:
$$
-k\Delta T(x) = -{\mathrm {div}}\, \mathbf q_r(x) + F(x),
$$
где $k$ -- коэффициент теплопроводности, $F(x)$ -- объемная плотность источников тепла, $\mathbf q_r$ -- вектор плотности потока излучения, определяемый формулой $\mathbf q_r(x)= \int_S I(x,\omega)\omega\,d\omega$.
На границе области поставим граничное условие III рода, выражающее закон Ньютона\,--\,Рихмана:
$$
k\frac{\partial T}{\partial n} + h(x)(T(x) - T_b(x)) = 0, \;\; x \in \Gamma,
$$
где $h(x)$ -- коэффициент теплоотдачи.

Далее для вывода приближения упрощенного метода сферических гармоник рассмотрим аналогичное уравнение, описывающее перенос излучения в плоскопараллельном слое \cite{SolvabKovt, Modest, CaseZweifel}:
\begin{equation}
\label{I1-eq}
\mu I_\tau(\tau, \mu) + I(\tau, \mu) = \frac{c}{2}\int_{-1}^1 p(\mu, \mu')I(\tau, \mu')\,d\mu' + (1 - c)I_b(T(\tau)),
\end{equation}
с граничными условиями
\begin{gather}
\label{I1-bc1}
I(0, \mu) = \varepsilon_1 I_b(T_{b1}) + \rho_1^s I(0, -\mu) + 2\rho_1^d\int_0^1 I(0, -\mu')\mu'\,d\mu', \quad \mu > 0,
\\
\label{I1-bc2}
I(\tau_0, \mu) = \varepsilon_2 I_b(T_{b2}) + \rho_2^s I(\tau_0, -\mu) + 2\rho_2^d\int_0^1 I(0, \mu')\mu'\,d\mu', \quad \mu < 0,
\end{gather}
где $I(\tau, \mu)$ -- интенсивность излучения в точке с оптической координатой $\tau \in [0,\tau_0]$ в направлении, составляющим с положительным направлением оси $\tau$ угол с косинусом $\mu \in [-1,1]$, $c= \kappa_s/\kappa$ -- альбедо однократного рассеяния, $p(\mu, \mu')$ -- фазовая функция, $T(\tau)$ -- абсолютная температура.

Будем искать приближенное решение уравнения (\ref{I1-eq}) в виде отрезка ряда Фурье по многочленам Лежандра $P_k(\mu)$, а также разложим по этим же функциям фазовую функцию рассеяния:
\begin{gather*}
I(\tau,\mu) = \sum_{k=0}^3 \phi_k(\tau)P_k(\mu) = \sum_{k=0}^3 \frac{2k+1}{4\pi} \widetilde\phi_k(\tau)P_k(\mu),
\quad
p(\mu, \mu') = \sum_{k=0}^3 A_k P_k(\mu')P_k(\mu).% = \sum_{k=0}^3 \frac{2k+1}{4\pi}\widetilde A_k P_k(\mu')P_k(\mu).
\end{gather*}
Подставив разложения функций $I$ и $p$ в уравнение (\ref{I1-eq}) и в граничные условия (\ref{I1-bc1}), (\ref{I1-bc2}) и умножая полученные соотношения на $P_l(\mu)$, $l = 1, 3$, а затем интегрируя по $\mu$, получим, следуя \cite{Modest, CaseZweifel}, следующие уравнения относительно функций $\phi_k(\tau)$:
\begin{gather}
\label{eq-a0}
\frac{1}{3}\phi_1' + (1 - cA_0)\phi_0 = (1-c)I_b(T(\tau)),
\\
\label{eq-a1}
\phi_0' + \frac{2}{5}\phi_2' + (1 - cA_1/3)\phi_1 = 0,
\\
\label{eq-a2}
\frac{2}{3}\phi_1' + \frac{3}{7}\phi_3' + (1-cA_2/5)\phi_2 = 0,
\\
\label{eq-a3}
\frac{3}{5}\phi_2' + (1 - cA_3/7)\phi_3 = 0
\end{gather}
вместе с граничными условиями Маршака:
\begin{gather*}
\sum_{k=0}^3 \left[ p_{kl} - (-1)^k\rho_1^s p_{kl} - 2(-1)^k\rho_1^dp_{0l}p_{1k}\right]\phi_k(0) = \varepsilon_1I_b(T_{b1})p_{0l}, \;\; l = 1, 3,
\\
\sum_{k=0}^3 \left[ p_{kl} - (-1)^k\rho_2^s p_{kl} - 2(-1)^k\rho_2^d p_{0l}p_{1k}\right](-1)^k\phi_k(\tau_0) = \varepsilon_2 I_b(T_{b2})p_{0l}, \;\; l = 1, 3,
\end{gather*}
где $p_{kl} = \int_0^1 P_k(\mu) P_l(\mu) \, d\mu$.

Переходя к функциям $\widetilde\phi_k$, получим из уравнений (\ref{eq-a1}), (\ref{eq-a3}) соотношения
\begin{gather*}
\widetilde \phi_1 = -\frac{1}{3}(1 - cA_1/3)^{-1}\widetilde\phi_0' - \frac{2}{3}(1 - cA_1/3)^{-1}\widetilde\phi_2',
\\
\widetilde\phi_3 = -\frac{3}{7}(1 - cA_3/7)^{-1}\widetilde\phi_2',
\end{gather*}
отсюда после подстановки в (\ref{eq-a0}), (\ref{eq-a2}) получаем уравнения
\begin{gather*}
-\frac{1}{3}(1 - cA_1/3)^{-1}\widetilde\phi_0'' - \frac{2}{3}(1 - cA_1/3)^{-1}\widetilde\phi_2'' + (1 - cA_0)\widetilde\phi_0 = 4\pi(1-c)I_b(T(\tau)),
\\
-\frac{2}{3}(1 - cA_1/3)^{-1}\widetilde\phi_0'' - \frac{4}{3}(1 - cA_1/3)^{-1}\widetilde\phi_2'' - \frac{9}{7}(1 - cA_3/7)^{-1}\widetilde\phi_2'' + 5(1 - cA_2/5)\widetilde\phi_2 = 0
\end{gather*}
и граничные условия
\begin{gather*}
-\frac{1}{3}(1 - cA_1/3)^{-1} \widetilde\phi_0'(0) - \frac{2}{3}(1 - cA_1/3)^{-1}\widetilde \phi_2'(0) + \frac{\varepsilon_1}{2(2 - \varepsilon_1)}(\widetilde\phi_0(0) - 4\pi I_b(T_{b1})) +
\\
+ \frac{5\varepsilon_1}{8(2 - \varepsilon_1)}\widetilde\phi_2(0) = 0,
\\
\frac{\rho_1^d}{4}\left(-\frac{1}{3}(1 - cA_1/3)^{-1}\widetilde\phi_0'(0) -
\frac{2}{3}(1 - cA_1/3)^{-1}\widetilde\phi_2'(0) \right)
+ \frac{3}{7}(1 + \rho_1^s)(1 - cA_3/7)^{-1}\widetilde\phi_2'(0) +
\\
+ \frac{\varepsilon_1}{8}(\widetilde\phi_0(0) - 4\pi I_b(T_{b1})) - \frac{5}{8}\left(1 - \rho_1^s + \frac{1}{4}\rho_1^d \right)\widetilde\phi_2(0) = 0,
%-\frac{1}{3}(1 - cA_1/3)^{-1}\widetilde\phi_0'(0) - \frac{2}{3}(1 - cA_1/3)^{-1}\widetilde\phi_2'(0)  + \frac{12(1 + \rho_1^s)}{7\rho_1^d(1 -cA_3/7)}\widetilde\phi_2'(0) +
%\\
%+ \frac{\varepsilon_1}{2\rho_1^d}(\widetilde\phi_0(0) - 4\pi I_b(T_{b1}))
%- \frac{5}{2\rho_1^d}(1 - \rho_1^s + \rho_1^d/4)\widetilde\phi_2(0) = 0,
\\
\frac{1}{3}(1 - cA_1/3)^{-1} \widetilde\phi_0'(\tau_0) + \frac{2}{3}(1 - cA_1/3)^{-1}\widetilde \phi_2'(\tau_0) + \frac{\varepsilon_2}{2(2 - \varepsilon_2)}(\widetilde\phi_0(\tau_0) - 4\pi I_b(T_{b2})) +
\\
+ \frac{5\varepsilon_2}{8(2 - \varepsilon_2)}\widetilde\phi_2(\tau_0) = 0,
\\
%\frac{1}{3}(1 - cA_1/3)^{-1}\widetilde\phi_0'(\tau_0) + \frac{2}{3}(1 - cA_1/3)^{-1}\widetilde\phi_2'(\tau_0)  - \frac{12(1 + \rho_2^s)}{7\rho_2^d(1 -cA_3/7)}\widetilde\phi_2'(\tau_0) +\\
%+ \frac{\varepsilon_2}{2\rho_2^d}(\widetilde\phi_0(\tau_0) - 4\pi I_b(T_{b2}))
%- \frac{5}{2\rho_2^d}(1 - \rho_2^s + \rho_2^d/4)\widetilde\phi_2(\tau_0) = 0.
\frac{\rho_2^d}{4}\left(\frac{1}{3}(1 - cA_1/3)^{-1}\widetilde\phi_0'(\tau_0) +
\frac{2}{3}(1 - cA_1/3)^{-1}\widetilde\phi_2'(\tau_0) \right)
- \frac{3}{7}(1 + \rho_2^s)(1 - cA_3/7)^{-1}\widetilde\phi_2'(\tau_0) +
\\
+ \frac{\varepsilon_2}{8}(\widetilde\phi_0(\tau_0) - 4\pi I_b(T_{b2})) - \frac{5}{8}\left(1 - \rho_2^s + \frac{1}{4}\rho_2^d \right)\widetilde\phi_2(\tau_0) = 0.
\end{gather*}
из которых следуют условия
\begin{gather*}
-\frac{2}{3}(1 - cA_1/3)^{-1}\widetilde\phi_0'(0)
- \frac{4}{3}(1 - cA_1/3)^{-1}\widetilde\phi_2'(0)
- \frac{9}{7}(1 - cA_3/7)^{-1}\widetilde\phi_2'(0) +
\\
+ \frac{5\varepsilon_1}{8(2 - \varepsilon_1)}(\widetilde\phi_0(0) - 4\pi I_b(T_{b1}))
+ \left(\frac{25\varepsilon_1}{8(2-\varepsilon_1)} + \frac{75\rho_1^d}{16(1+\rho_1^s)(2-\varepsilon_1)} \right)\widetilde\phi_2(0) = 0,
\\
\frac{2}{3}(1 - cA_1/3)^{-1}\widetilde\phi_0'(\tau_0)
+ \frac{4}{3}(1 - cA_1/3)^{-1}\widetilde\phi_2'(\tau_0)
+ \frac{9}{7}(1 - cA_3/7)^{-1}\widetilde\phi_2'(\tau_0) +
\\
+ \frac{5\varepsilon_2}{8(2 - \varepsilon_2)}(\widetilde\phi_0(\tau_0) - 4\pi I_b(T_{b2}))
+ \left(\frac{25\varepsilon_2}{8(2-\varepsilon_2)} + \frac{75\rho_2^d}{16(1+\rho_2^s)(2-\varepsilon_2)} \right)\widetilde\phi_2(\tau_0) = 0.
\end{gather*}

Функция $\widetilde\phi_0 = 4\pi\phi_0$ аппроксимирует плотность энергии излучения $\int_{-1}^1 I(\tau, \mu)\,d\mu$, функция $\widetilde\phi_1 = 4\pi\phi_1/3$ аппроксимирует величину плотности потока излучения $\int_{-1}^1 I(\tau,\mu)\mu\,d\mu$, функции $\phi_2$ и $\phi_3$ не имеют простого физического смысла.

Чтобы теперь получить уравнения $SP_3$-приближения, заменим в уравнениях плоскопараллельного слоя одномерные дифференциальные операторы на многомерные и перейдем обратно от переменной $\tau$ к переменной $x$. Будем иметь
\begin{gather*}
-\alpha_{11}\Delta\widetilde\phi_0 - \alpha_{12}\Delta\widetilde\phi_2 + \kappa_a\widetilde\phi_0 = 4\pi\kappa_aI_b(T(x)),
\\
-\alpha_{21}\Delta\widetilde\phi_0 - \alpha_{22}\Delta\widetilde\phi_2 + \widetilde\kappa\widetilde\phi_2 = 0,
\\
\alpha_{11}\frac{\partial\widetilde\phi_0}{\partial n}
+ \alpha_{12}\frac{\partial\widetilde\phi_2}{\partial n}
+ \frac{\varepsilon}{2(2-\varepsilon)}(\widetilde\phi_0 - 4\pi I_b(T_b))
+ \frac{5\varepsilon}{8(2 - \varepsilon)}\widetilde\phi_2 = 0,
\\
\alpha_{21}\frac{\partial\widetilde\phi_0}{\partial n} + \alpha_{22}\frac{\partial\widetilde\phi_2}{\partial n}
+ \frac{\varepsilon}{8(2 - \varepsilon)}(\widetilde\phi_0 - 4\pi I_b(T_b))
+ \left(\frac{5\varepsilon}{8(2-\varepsilon)} + \frac{15\rho^d}{16(1+\rho^s)(2-\varepsilon)} \right)\widetilde\phi_2 = 0.
\end{gather*}
Здесь $\alpha_{11} = (3\kappa(1 - cA_1/3))^{-1}$, $\alpha_{12} = 2\alpha_{11}$, $\alpha_{21} = 2\alpha_{11}/5$, $\alpha_{22} = 4\alpha_{11}/5 + \dfrac{9}{35\kappa(1-cA_3/7)}$, $\widetilde\kappa = \kappa(1 - cA_2/5)$.

Дополним $SP_3$-приближение уравнения переноса излучения уравнением теплопроводности, введя нормированные величины
$\theta = T/T_{\max}$, $\varphi = \widetilde\phi_0 / 4\sigma n^2 T_{\max}^4$, $\xi = \widetilde \phi_2 / 4\sigma n^2 T_{\max}^4$.
Приходим к краевой задаче
\begin{gather}
\label{a1}
-a\Delta \theta + b\kappa_a(\theta^4 - \varphi) = 0,
\\
-\alpha_{11}\Delta \varphi - \alpha_{12}\Delta \xi + \kappa_a(\varphi - \theta^4) = 0,
\\
-\alpha_{21}\Delta\varphi - \alpha_{22}\Delta\xi + \widetilde\kappa \xi = 0,
\\
\label{a2}
a\frac{\partial\theta}{\partial n} + \beta(\theta - \theta_b) = 0,
\\
\alpha_{11}\frac{\partial\varphi}{\partial n} + \alpha_{12}\frac{\partial\xi}{\partial n} + \gamma_{11}(\varphi - \theta_b^4) + \gamma_{12}\xi = 0,
\\
\alpha_{21}\frac{\partial\varphi}{\partial n} + \alpha_{22}\frac{\partial\xi}{\partial n} + \gamma_{21}(\varphi - \theta_b^4) + \gamma_{22}\xi = 0.
\end{gather}
Здесь $a = k/\rho c_p$, $b = 4\sigma n^2 T_{\max}^3/\rho c_p$, $\beta = h/\rho c_p$, $\theta_b = T_b/T_{\max}$,
$$
\gamma_{11} = \frac{\varepsilon}{2(2-\varepsilon)},\;\;
\gamma_{12} = \frac{5\varepsilon}{8(2 - \varepsilon)},\;\;
\gamma_{21} = \frac{\varepsilon}{8(2-\varepsilon)},\;\;
\gamma_{22} = \frac{5\varepsilon}{8(2-\varepsilon)} + \frac{15\rho^d}{16(1+\rho^s)(2-\varepsilon)}.
$$


Чтобы избавиться от перекрестной диффузии в уравнениях, произведем преобразование координат $\begin{pmatrix} \psi_1 \\ \psi_2 \end{pmatrix} = T\begin{pmatrix} \varphi \\ \xi \end{pmatrix}$, где $T = \begin{pmatrix} 1 & \lambda_1 \\ 1 & \lambda_2 \end{pmatrix}$, выбрав, как и в \cite{Larsen2002, TsePinnau}, числа $\lambda_1$, $\lambda_2$ из условия равенства нулю внедиагональных коэффициентов в уравнениях. При этом матрицы коэффициентов уравнений $\widehat\alpha = \begin{pmatrix} \alpha_{11} & \alpha_{12} \\ \alpha_{21} & \alpha_{22} \end{pmatrix}$ и граничных условий $\widehat \gamma = \begin{pmatrix}\gamma_{11} & \gamma_{12} \\ \gamma_{21} & \gamma_{22} \end{pmatrix}$ преобразуются, соответственно, к матрицам $T\widehat\alpha T^{-1}$ и $T \widehat\gamma T^{-1}$.
Приходим к соотношениям
\begin{gather}
\label{a3}
-\alpha_1\Delta\psi_1 + \widetilde\kappa\psi_1 = \widetilde\kappa_s\varphi + \kappa_a\theta^4,
\\
\label{a4}
-\alpha_2\Delta\psi_2 + \widetilde\kappa\psi_2 = \widetilde\kappa_s\varphi + \kappa_a\theta^4,
\\
\label{a5}
\alpha_1\frac{\partial\psi_1}{\partial n} + \beta_{11}\psi_1 + \beta_{12}\psi_2 = \eta_1\theta_b^4,
\\
\label{a6}
\alpha_2\frac{\partial\psi_2}{\partial n} + \beta_{21}\psi_1 + \beta_{22}\psi_2 = \eta_2\theta_b^4,
\end{gather}
где
$$
\lambda_i = \frac{\alpha_{22} - \alpha_{11}}{2\alpha_{21}} + \frac{(-1)^i}{2\alpha_{21}}\sqrt{(\alpha_{11} - \alpha_{22})^2 + 4\alpha_{12}\alpha_{21}},
\;\;
\omega_1 = \frac{\lambda_2}{\lambda_2 - \lambda_1}, \;\;
\omega_2 = \frac{\lambda_1}{\lambda_1 - \lambda_2},
$$
$$
\alpha_i = \frac{\alpha_{11} + \alpha_{22}}{2} + (-1)^i\sqrt{(\alpha_{11} - \alpha_{22})^2 + 4\alpha_{12}\alpha_{21}} > 0,
$$
$\widetilde\kappa_s = \widetilde\kappa - \kappa_a \geq 0$,
$\varphi = \omega_1\psi_1 + \omega_2\psi_2$,
$\omega_1 + \omega_2 = 1$, $\omega_{1,2} \geq 0$,
$\begin{pmatrix}\eta_1 \\ \eta_2\end{pmatrix} = T\begin{pmatrix}\gamma_{11} \\ \gamma_{21} \end{pmatrix}$,
и поскольку
\begin{multline*}
\widehat\beta = \begin{pmatrix}\beta_{11} & \beta_{12} \\ \beta_{21} & \beta_{22}\end{pmatrix} = T\widehat\gamma T^{-1} = 
\\ = \frac{1}{\lambda_2 - \lambda_1}
\begin{pmatrix}
(\gamma_{11} + \lambda_1\gamma_{21})\lambda_2 - (\gamma_{12} + \lambda_1\gamma_{22}) & -(\gamma_{11} + \lambda_1\gamma_{21})\lambda_1 + (\gamma_{12} + \lambda_1\gamma_{22}) \\
(\gamma_{11} + \lambda_2\gamma_{21})\lambda_2 - (\gamma_{12} + \lambda_2\gamma_{22}) & -(\gamma_{11} + \lambda_2\gamma_{21})\lambda_1 + (\gamma_{12} + \lambda_2\gamma_{22})
\end{pmatrix},
\end{multline*}
то $\eta_i = \beta_{i1} + \beta_{i2}$.

Справедливы следующие свойства коэффициентов в граничных условиях при $\varepsilon > 0$.
\begin{enumerate}
\item
Диагональные коэффициенты строго положительны: $\beta_{ii} > 0$.
Коэффициент $\beta_{12}$ отрицателен при
$$
\lambda_1\frac{\rho^d}{\varepsilon(1 + \rho^s)} < \frac{2}{15}(\lambda_1^2 - \lambda_1 - 5).
$$
Коэффициент $\beta_{21}$ отрицателен при
$$
\lambda_2\frac{\rho^d}{\varepsilon(1 + \rho^s)} > \frac{2}{15}(\lambda_2^2 - \lambda_2 - 5).
$$
\item
Матрица $\widehat\beta$ имеет строгое диагональное преобладание.
\item
Выполняются неравенства $\omega_1\beta_{11} + \omega_2\beta_{21} > 0$, $\omega_1\beta_{12} + \omega_2\beta_{22} > 0$.
\end{enumerate}

Для формализации краевой задачи (\ref{a1}), (\ref{a2}), (\ref{a3})--(\ref{a6}) будем использовать пространство Соболева $V = H^1(\Omega)$. Через $(f, v)$ обозначаем значение функционала $f \in V'$ на элементе $v \in V$ и скалярное произведение в $L^2(\Omega)$, если $f, v \in L^2(\Omega)$. Будем предполагать, что исходные данные удовлетворяют условиям:

(i) \;$\beta_{ij} = \mathrm{const}$, $\beta_{ii} > 0$, матрица $\widehat \beta$ имеет строгое диагональное преобладание; $\beta \geq \beta_0 = \mathrm{const} > 0$, $\beta \in L^\infty(\Gamma)$; $0 \leq \theta_b \in L^\infty(\Gamma)$;

(ii)\;$\omega_1\beta_{11} + \omega_2\beta_{21} > 0$, $\omega_1\beta_{12} + \omega_2\beta_{22} > 0$.

Определим операторы $A \colon V \to V'$, $B_{1,2} \colon V \times V \to V'$ и функционалы $h, g_{1,2} \in V'$ по следующим формулам, справедливым для любого $v \in V$:
\begin{gather*}
(A\theta, v) = a(\nabla\theta, \nabla v) + \int_\Gamma \beta\theta v\, d\Gamma,
\\
(B_1\{\psi_1, \psi_2\}, v) = \alpha_1 (\nabla\psi_1, \nabla v) + \int_\Gamma (\beta_{11}\psi_1 + \beta_{12}\psi_2)v\,d\Gamma,
\\
(B_2\{\psi_1, \psi_2\}, v) = \alpha_2 (\nabla\psi_2, \nabla v) + \int_\Gamma (\beta_{21}\psi_1 + \beta_{22}\psi_2)v\,d\Gamma,
\\
(h, v) = \int_\Gamma \beta\theta_b v\,d\Gamma, \quad
(g_1, v) = \int_\Gamma \eta_1\theta_b^4v\,d\Gamma, \quad
(g_2, v) = \int_\Gamma \eta_2\theta_b^4v\,d\Gamma.
\end{gather*}

\begin{definition}
Тройка $\{\theta, \psi_1, \psi_2\} \in V \times V \times V$ называется слабым решением  краевой задачи (\ref{a1}), (\ref{a2}), (\ref{a3})--(\ref{a6}), если
\begin{gather}
\label{weak1}
A_1\theta + b\kappa_a(|\theta|^3\theta - \varphi) = h, \\
\label{weak2}
B_1\{\psi_1, \psi_2\} + \widetilde\kappa\psi_1 = \widetilde\kappa_s\varphi + \kappa_a|\theta|^3\theta + g_1, \quad
B_2\{\psi_1, \psi_2\} + \widetilde\kappa\psi_2 = \widetilde\kappa_s\varphi + \kappa_a|\theta|^3\theta + g_2,
\end{gather}
где $\varphi = \omega_1\psi_1 + \omega_2\psi_2$.
\end{definition}


\section{Существование решения}

В \cite{TsePinnau} установлено следующее свойство решения краевой задачи (\ref{a3})--(\ref{a6}) при заданном поле температуры $\theta$.

\begin{theorem}
\label{thm1}
Пусть выполняются условия (i), $0 \leq \theta \leq M$, $M = \|\theta_b\|_{L^\infty(\Gamma)}$.
Тогда существует единственное решение $\{\psi_1, \psi_2\}$ задачи (\ref{weak2}), причем $0 \leq \psi_{1,2} \leq M^4$.
\end{theorem}

Доказательство существования слабого решения задачи (\ref{a1}), (\ref{a2}), (\ref{a3})--(\ref{a6}) проведем методом, предложенным в \cite{KovtCheb}, а именно, построим итерационный процесс, порождающий монотонные и ограниченные последовательности приближений метода простой итерации, которые, как мы докажем, дают в пределе искомое слабое решение. Кроме того, мы покажем, что для предела $\{\theta_*, \psi_{1*}, \psi_{2*}\}$ и любого ограниченного слабого решения $\{\widehat\theta, \widehat \psi_1, \widehat\psi_2\}$ выполняются неравенства $\theta_* \leq \widehat\theta$, $\psi_{1*} \leq \widehat\psi_1$, $\psi_{2*} \leq \widehat\psi_2$. Этот факт будет использован в доказательстве единственности ограниченного решения.

Определим нелинейный оператор $F \colon L^\infty(\Omega) \to L^\infty(\Omega) \cap V$, который по заданной функции $\varphi$ дает решение $\theta = F(\varphi)$ задачи
(\ref{weak1}), а также нелинейные операторы $G_{1,2} \colon L^\infty(\Omega) \to L^\infty(\Omega) \cap V$, которые по заданной функции $\theta$ дают решения $\psi_i = G_i(\theta)$ задачи (\ref{weak2}).

Справедливы следующие свойства операторов $F$ и $G$, которые доказываются по аналогии с \cite{KovtCheb} с использованием теоремы \ref{thm1}.

\begin{lemma}
Если $\theta_1 \leq \theta_2$ и $\varphi_1 \leq \varphi_2$ п.в. в $\Omega$, то $F(\varphi_1) \leq F(\varphi_2)$ и $G_i(\theta_1) \leq G_i(\theta_2)$ п.в. в $\Omega$.
\end{lemma}

\begin{lemma}
Пусть $\{\widehat\theta, \widehat\psi_1, \widehat\psi_2\}$ --- любое неотрицательное, ограниченное решение задачи (\ref{weak1}), (\ref{weak2}). Тогда, если $\{\theta_k, \psi_{1k}, \psi_{2k}\}$ --- последовательность приближений метода простой итерации $\theta_0 = 0$, $\psi_{1k} = G_1(\theta_k)$, $\psi_{2k} = G_2(\theta_k)$, $\theta_{k+1} = F(\varphi_k)$, где $\varphi_k = \omega_1\psi_{1k} + \omega_2\psi_{2k}$, то
$0 \leq \theta_k \leq \theta_{k+1} \leq \widehat \theta$,
$0 \leq \psi_{1k} \leq \psi_{1,k+1} \leq \widehat \psi_1$,
$0 \leq \psi_{2k} \leq \psi_{2,k+1} \leq \widehat \psi_2$,
$k = 1, 2, \ldots$.
\end{lemma}

\begin{remark}
В формулировке леммы 2 можно опустить слова <<неотрицательное решение>>. Действительно, нулевую нижнюю границу можно заменить на $C_0 = \mathrm{ess\,\inf}\, \widehat\theta$ и положить $\theta_0 = C_0$. Доказательство практически не изменится.
\end{remark}

Заметим, что последовательность приближений метода простой итерации обладает свойством ограниченности: $0 \leq \theta \leq M$, $0 \leq \psi_{1,2} \leq M^4$, где $M = \|\theta_b\|_{L^\infty(\Gamma)}$. Таким образом, в пределе получаем решение, обладающее этим же свойством. 
Из лемм 1, 2 вытекает следующий результат.
\begin{theorem}
Пусть выполняются условия (i). Тогда существует решение $\{\theta, \psi_1, \psi_2\}$ задачи (\ref{weak1}), (\ref{weak2}), обладающее свойством $0 \leq \theta \leq M$, $0 \leq \psi_{1,2} \leq M^4$, где $M = \|\theta_b\|_{L^\infty(\Gamma)}$.
\end{theorem}


\section{Единственность решения}

\begin{theorem}
Пусть выполняются условия (i), (ii). Тогда существует единственное ограниченное решение задачи (\ref{weak1})--(\ref{weak2}).
\end{theorem}

\begin{proof}

Пусть $\{\theta_*, \psi_{1*}, \psi_{2*}\}$ --- предел последовательности приближений метода простой итерации, и
$\{\widehat\theta, \widehat\psi_1, \widehat\psi_2\}$ --- произвольное решение задачи (\ref{weak1})--(\ref{weak2}) из класса $[L^\infty(\Omega)]^3$.
Рассмотрим разность этих двух решений $\theta = \widehat\theta - \theta_*$, $\psi_1 = \widehat\psi_1 - \psi_{1*}$, $\psi_2 = \widehat\psi_2 - \psi_{2*}$. По лемме 2 имеем $\theta \geq 0$, $\psi_{1,2} \geq 0$.

Функции $\theta, \psi_1, \psi_2$ подчиняются уравнениям
\begin{gather}
\label{g1}
A_1\theta + b\kappa_a((|\theta_1|^3\theta_1 - |\theta_2|^3 \theta_2) - \varphi) = 0, \\
\label{g2}
B_1\{\psi_1, \psi_2\} + \widetilde\kappa\psi_1 = \widetilde\kappa_s\varphi + \kappa_a(|\theta_1|^3\theta_1 - |\theta_2|^3 \theta_2), \\
\label{g3}
B_2\{\psi_1, \psi_2\} + \widetilde\kappa\psi_2 = \widetilde\kappa_s\varphi + \kappa_a(|\theta_1|^3\theta_1 - |\theta_2|^3 \theta_2),
\end{gather}
где $\varphi = \omega_1\psi_1 + \omega_2\psi_2$.

Умножим (\ref{g2}) скалярно на $v = \omega_1$, (\ref{g3}) на $v = \omega_2$, (\ref{g1}) на $v = 1/b$ и сложим результаты. Получим
$$
\frac{1}{b}\int_\Gamma \beta\theta\,d\Gamma + \int_\Gamma (\omega_1 \beta_{11} + \omega_2\beta_{21})\psi_1\,d\Gamma
+ \int_\Gamma (\omega_1\beta_{12} + \omega_2\beta_{22})\psi_2\,d\Gamma = 0,
$$
отсюда $\theta|_\Gamma = \psi_1|_\Gamma = \psi_2|_\Gamma = 0$.

Умножив уравнение (\ref{g2}) на $b\omega_1$, уравнение (\ref{g3}) на $b\omega_2$ и сложив их сумму с уравнением (\ref{g1}), приходим к тождеству
$$
a(\nabla\theta, \nabla v) + b\omega_1\alpha_1(\nabla\psi_1, \nabla v) + b\omega_2\alpha_2(\nabla\psi_2, \nabla v) = 0,
$$
которое справедливо для любого $v \in V$.
Полагая $v = a\theta + b\alpha_1\omega_1\psi_1 + b\alpha_2\omega_2\psi_2$, получаем $a\theta + b\alpha_1\omega_1\psi_1 + b\alpha_2\omega_2\psi_2 = 0$ п.в. в $\Omega$, отсюда с учетом неотрицательности $\theta$, $\psi_1$, $\psi_2$ вытекает, что $\theta = \psi_1 = \psi_2 = 0$ п.в. в $\Omega$.
\end{proof}


\section{Примеры}

Приведем результаты вычислительного эксперимента при физических параметрах, соответствующих воздуху: $\Omega$ -- квадрат со стороной $L = 1$, $k = 0.0515$, $\kappa = 10$, $\kappa_a = 1$, $T_{\max} = 773$, $n = 1$, $h = 10$, $\varepsilon = 0.5$, $\rho^d = 0.2$, рассеяние в среде считаем изотропным. Граничная температура $\theta_b$ имеет вид $\theta_b(x, 0) = 0.5$, $\theta_b(x, L) = 1$, $\theta_b(0, y) = \theta_b(L, y) = 0.5 + y/(2L)$. Среднеквадратичное отклонение между полями температуры в $P_1$- и $SP_3$-моделях составило 0.00057 для $\theta$ и 0.0017 для $\varphi$. При этом наибольшее отклонение между полями достигается ближе к нижней границе квадрата.
%Поле температуры %и усредненной по направлениям интенсивности излучения
%для $SP_3$-модели представлено на рис. 1.% на рис. 1, 2.
Если теперь положить $\varepsilon = 0.75$, $\rho_d = 0.2$, то аналогичные характеристики составят, соответственно, 0.00051 и 0.0015.
При $\rho_d = 0$ характеристики отклонения между моделируемыми полями 0.0006 и 0.0018.

%\begin{center}
%\includegraphics[width=10cm]{1.eps}

%Рис. 1. Поле $\theta$ для $SP_3$-модели сложного теплообмена.
%\end{center}

%\begin{center}
%\includegraphics[width=10cm]{2.eps}

%Рис. 2. Поле $\varphi$ для $SP_3$-модели сложного теплообмена.
%\end{center}

Наконец, приведем результаты моделирования процесса сложного теплообмена при физических параметрах, соответствующих стеклу. Положим $L = 0.25$, $k = 1.2$, $\kappa = 100$, $\kappa_a = 100$, $T_{\max} = 773$, $n = 1.47$, $h = 10$, $\varepsilon = 0.75$, $\rho^d = 0.2$. Граничная температура вычисляется по той же формуле, что и выше. Отклонения между полями $\theta$ и $\varphi$ в $P_1$- и $SP_3$-моделях составляют, соответственно, 0.00026 и 0.00048. Если теперь положить $\varepsilon = 0.75$, $\rho_d = 0$, то эти характеристики будут равны 0.00026 и 0.00052.

Из проведенных вычислений можно сделать вывод, что наблюдается тенденция увеличения отклонения между моделируемыми полями в $P_1$- и $SP_3$ моделях при преобладании зеркального отражения.

Для вычислений применялся пакет FreeFem++ \cite{FreeFem}, использован упрощенный метод Ньютона квазилинеаризации уравнений сложного теплообмена \cite{NewtonMethod} наряду с полным методом Ньютона. Скорость сходимости упрощенного метода Ньютона в рассмотренных примерах аналогична сходимости для $P_1$-модели. Программная реализация доступна по ссылке https://github.com/grenkin/test-newton.


\bigskip
\begin{thebibliography}{1} % Ссылки на литературу приводятся только на английском языке (даже если у источника нет англоязычной версии).

\bibitem{Cheb2}
A.Yu. Chebotarev, A.E. Kovtanyuk, N.D. Botkin,
{\it Problem of radiation heat exchange with boundary conditions of the Cauchy type,}
Commun. Nonlinear Sci. Numer. Simul. {\bf 75} (2019), 262--269.

\bibitem{LaserAbl}
A.E. Kovtanyuk, A.Yu. Chebotarev, A.A. Astrakhantseva, A.A. Sushchenko, 
{\it Optimal control of endovenous laser ablation,}
Optics and Spectroscopy, {\bf 128}:9 (2020), 1508--1516.

\bibitem{LaserAbl2}
A. Kovtanyuk, A. Chebotarev, A. Astrakhantseva,
{\it Inverse extremum problem for a model of endovenous laser ablation,}
J. Inverse Ill-Posed Probl., {\bf 29}:3 (2021), 467--476.

\bibitem{Schneider}
E. Schneider, M. Seaid, J. Janicka, A. Klar, {\it Validation of simplified $P_N$ models for radiative transfer in combustion systems,} Commun. Numer. Meth. Engng., {\bf 24} (2008), 85--96.

\bibitem{Larsen2002}
E.W. Larsen, G. Th\"ommes, A. Klar, M. Seaid, T. G\"otz,
{\it Simplified $P_N$ approximations to the equations of radiative heat transfer and applications}, J. Comp. Phys., {\bf 183}:2 (2002), 652--675.

\bibitem{ModestLei}
M.F. Modest, S. Lei, {\it The simplified spherical harmonics method for radiative heat transfer,} Journal of Physics: Conference Series. IOP Publishing, {\bf 369}:1 (2012), 012019.

\bibitem{KovtMC}
A.E. Kovtanyuk, N.D. Botkin, K.-H. Hoffmann, {\it Numerical simulations of a coupled radiative-conductive heat transfer model using a modified Monte Carlo method}, Int. J. Heat Mass Transf., {\bf 55}:4 (2012), 649--654.

\bibitem{ZhengHan}
H. Zheng, W. Han, {\it On simplified spherical harmonics equations for the radiative transfer equation}, J. Math. Chem., {\bf 49}:8 (2011), 1785--1797.

\bibitem{TsePinnau}
R. Pinnau, O. Tse, {\it Optimal control of a simplified natural convection-radiation model}, Comm. Math. Sci., {\bf 11}:3 (2013), 679--707.

\bibitem{KovtCheb}
A.E. Kovtanyuk, A.Yu. Chebotarev, N.D. Botkin, K.-H. Hoffmann,
{\it Unique solvability of a steady-state complex heat transfer model,}
Commun. Nonlinear Sci. Numer. Simul., {\bf 20}:3 (2015), 776--784.


\bibitem{SolvabKovt}
A.E. Kovtanyuk, A.Yu. Chebotarev, N.D. Botkin, K.-H. Hoffmann,
{\it Solvability of $P_1$ approximation of a conductive-radiative heat transfer problem,} Appl. Math. Comput., {\bf 249} (2014), 247--252.

\bibitem{Inhomog}
A.Yu. Chebotarev, G.V. Grenkin, A.E. Kovtanyuk, {\it Inhomogeneous steady-state problem of complex heat transfer,} ESAIM Math. Model. Numer. Anal. (2017), 2511--2519.

\bibitem{Amosov1}
A.A. Amosov, {\it Nonstationary problem of complex heat transfer in a system of semitransparent bodies with radiation diffuse reflection and refraction boundary-value conditions,} Proceedings of the Seventh International Conference on Differential and Functional-Differential Equations (Moscow, August 22--29, 2014). Part 2, CMFD, 59, PFUR, M., (2016), 5--34.

\bibitem{Amosov2}
A.A. Amosov, {\it Stationary problem of complex heat transfer in a system of semitransparent bodies with boundary conditions of diffuse reflection and refraction of radiation}, Comput. Math. Math. Phys., {\bf 57}:3 (2017), 515--540.

\bibitem{SPN2}
P.S. Brantley, E.W. Larsen, {\it The simplified $P_3$ approximation,} Nuclear science and engineering, {\bf 134}:1 (2000), 1--21.

\bibitem{SPN1}
R. Backofen, T. Bilz, A. Ribalta, A. Voigt, {\it $SP_N$-approximations of internal radiation in crystal growth of optical materials,} Journal of Crystal Growth, {\bf 266}:1--3 (2004), 264--270.

\bibitem{SPN4}
A. Klar, J. Lang, M. Seaid, {\it Adaptive solutions of $SP_N$-approximations to radiative heat transfer in glass,} Int. J. Therm. Sci., {\bf 44} (2005), 1013--1023.

\bibitem{KloseLarsen}
A.D. Klose, E.W. Larsen, {\it Light transport in biological tissue based on the simplified spherical harmonics equations,} J. Comput. Phys., {\bf 220} (2006), 441--470.

\bibitem{SPN3}
M. Chu, K. Vishwanath, A.D. Klose, H. Dehghani, {\it Light transport in biological tissue using three-dimensional frequency-domain simplified spherical harmonics equations,} 
Phys. Med. Biol., {\bf 54} (2009), 2493--2509.

\bibitem{SPN5}
R.G. McClarren, {\it Theoretical aspects of the simplified $SP_N$ equations,}
Transport Theory Statist. Phys. {\bf 39} (2011), 73--109.

\bibitem{SPN6}
M.F. Modest, J. Cai, W. Ge, E. Lee, {\it Elliptic formulation of the simplified spherical harmonics method in radiative heat transfer,} Int. J. Heat Mass Transf., {\bf 76} (2014), 459--466.

\bibitem{Modest}
M.F. Modest, {\it Radiative heat transfer}, Academic Press, 2013.

\bibitem{CaseZweifel}
K.M. Case, P.F. Zweifel, {\it Linear transport theory,} Addison-Wesley, 1972.

\bibitem{FreeFem}
F. Hecht, {\it New development in FreeFem++,} J. Numer. Math., {\bf 20}:3--4 (2012), 251--265.

\bibitem{NewtonMethod}
G.V. Grenkin, 
{\it Convergence of Newton's method for equations of complex heat transfer},
Dal'nevost. Mat. Zh., {\bf 17}:1 (2017), 3--10.


\end{thebibliography}
\end{document}
