Теплопроводность в сферических координатах (численно)

[cool.phenon]

Здравствуйте, уважаемые участники форума. Возник вопрос: есть краевая задача

$$\frac{\partial u}{\partial t}= \operatorname{div}\left(a^2 E \nabla u\left(r,\varphi,\psi \right) \right) $$

в шаре

$$B_{0}(1)=\left\{x \in {{\mathbb{R}}^{3}: \left\| x \right\| < 1 \right\} $$

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

$$\frac{\partial u}{\partial t}=a^2\left(\frac{1}{r^{2}} \frac{\partial}{\partial r}\left(r^2 \frac{\partial u}{\partial r}\right)+\frac{1}{r\sin\varphi}\frac{\partial }{\partial \varphi}\left(\frac{\sin\varphi}{r} \frac{\partial u}{\partial \varphi} \right)+\frac{1}{r\sin\varphi}\frac{\partial}{\partial \psi}\left(\frac{1}{r\sin\varphi}\frac{\partial u}{\partial \psi}} \right)\right). $$

Эту краевую задачу нужно решить численно. Для этого можно написать неявную схему

$$\frac{u^{n+1}_{i,j,k}-u^{n}_{i,j,k}}{\Delta t}=a^2 \left( \frac{2}{r_{1}+(i-1)\Delta r}\frac{u^{n+1}_{i+1,j,k}-u^{n+1}_{i-1,j,k}}{2\Delta r}+\frac{u^{n+1}_{i+1,j,k}-2u^{n+1}_{i,j,k}+u^{n+1}_{i-1,j,k}}{\Delta r^2}\right) +$$ $$+a^2\left(\frac{\ctg \left(\varphi_{1}+(j-1)\Delta \varphi \right) ) }{\left(r_{1}+(i-1)\Delta r\right)^2}\frac{u^{n+1}_{i,j+1,k}-u^{n+1}_{i,j-1,k}}{2\Delta \varphi}+\frac{1}{\left(r_{1}+(i-1)\Delta r\right)^2}\frac{u^{n+1}_{i,j+1,k}-2u^{n+1}_{i,j,k}+u^{n+1}_{i,j-1,k}}{\Delta \varphi^{2}} \right) +$$ $$+a^2\left(\frac{1}{\left(r_{1}+(i-1)\Delta r\right)^2 \sin^2\left(\varphi_{1}+(j-1)\Delta \varphi \right)}\frac{u^{n+1}_{i,j,k+1}-2u^{n+1}_{i,j,k}+u^{n+1}_{i,j,k-1}}{\Delta \psi^2} \right).$$

Проблема состоит в её реализации. Дело в том, что величины $r_{1},\varphi_{1}$ -- суть нулевые, поэтому при вычислении соответствующих значений возникает деление на нуль. Подскажите, пожалуйста, как можно переписать эту схему, чтобы избежать описанных неопределённостей?

[ewert]

В схему не вникал, но стандартный подход для этой ситуации -- использование полуцелых узлов. Т.е. $r_i=h(i+\frac12)$. Тогда все особенности в нуле автоматически исчезнут.

[cool.phenon]

То есть, если раньше мы принимали первые значения $r_{1}=0$ и $\varphi_{1}=-\frac{\pi}{2}$, то сейчас нужно изменить эти значения на полшага и принять $r_{1}=\frac{\Delta r}{2}, \varphi_{1}=\frac{\Delta \varphi}{2}-\frac{\pi}{2}$?

[sergei1961]

ewert-интересный приём, не встречал. Это где-то написано в канонических текстах?

[cool.phenon]

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

[sergei1961]

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

[ewert]

cool.phenon в сообщении #751442 писал(а):

а как быть с краевым условием на шаре, который возник при отступлении полушага?

А кто заставляет отступать на полшага именно от края?...

[cool.phenon]

ewert
Но мы отступали-то не от края, а от центра. С краевым условием на поверхности шара разобраться можно. Но если отступить от центра, то получим небольшую (на расстоянии как раз в полшага) пустоту в виде сферы, за неё я и говорил.

sergei1961
С этим методом я не знаком, но его реализация видится такой же, как и предыдущего метода. Размышлял я так: если из этого уравнения с помощью асимптотики (или каких-то других методов) вытащить решение непосредственно в одной точке -- центре сферы (естественно, что там будет не нуль) и задать значение в центре как "точечное" условие, то все проблемы уйдут. Фактически, это же условие можно было бы задать и на той пустоте, но действовать там нужно было бы осторожнее, ведь нужно было бы находить не значение в точке, а на сфере малого радиуса.

Основная трудность, с которой я столкнулся -- как именно получить значение в центре сферы? Дело в том, что с асимптотиками я очень мало работал, как именно их применять в краевых задачах, я не знаю.
sergei1961
Если Вам не сложно, подкиньте, пожалуйста, идею насчёт асимптотического подхода.

[sergei1961]

В центре сферы-а на самой сфере известны данные?

[cool.phenon]

sergei1961
Да, краевое условие известно.На поверхности сферы функция нулевая. Известно также начальное условие, но это произвольная функция.

[sergei1961]

Для уравнений типа Лапласа подошло бы тогда использование теоремы о среднем: грубо говоря, значение в центре сферы равно некоторому интегралу по сфере. Для Вашего типа уравнения может быть попробовать доказать аналог теоремы о среднем? Или поискать готовый-это очень обширная область, теоремы о среднем.

[cool.phenon]

Я понял, что Вы имели в виду. Использовать утверждение: если функция непрерывна в сфере, то

$$\int\limits_{B_{1}(0)}f(x,y,z)dxdydz=V(B)f(\xi) $$

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

[ewert]

cool.phenon в сообщении #751689 писал(а):

нужно знать значение в центре сферы со временем,

Не нужно. Всё, что Вам нужно -- это просчитать поле в целом, значение же для центра круга в разностной схеме со сдвигом на полшага нигде не используются.

[cool.phenon]

ewert
Если значение в центре известно -- тогда всё нормально. Если сдвинуться на полшага, возникает новая граница. (та, о которой я говорил выше)

[ewert]

cool.phenon в сообщении #751704 писал(а):

Если значение в центре известно -- тогда всё нормально

Тогда всё ненормально, так краевые задачи не ставятся.

cool.phenon в сообщении #751704 писал(а):

Если сдвинуться на полшага, возникает новая граница.

Никакой новой границы не возникает. Вам нужно рассчитать массив $u_{1/2},u_{3/2},u_{5/2},\ldots$. В разностном уравнении для начального узла $r_{1/2}$ формально присутствуют неизвестные $u_{-1/2},u_{1/2},u_{3/2}$. Однако фактически разность $u_{1/2}-u_{-1/2}$ в этом уравнении умножается на $r_0=0$, и потому значение $u_{-1/2}$ не требуется.