Численное решение двумерного нестационарного уравнения тепло

[Munin]

svv в сообщении #1206085 писал(а):

Выбор $k_1\neq k_2$ означает, что коэффициент теплопроводности анизотропен — различен в разных направлениях. Но если так, то, скорее всего, всё ещё хуже

Такое уравнение может получиться из изотропного преобразованием координат. А координаты в технических вычислениях стараются подбирать пусть и криволинейные, но локально прямоугольные (например, сферические, цилиндрические, эллиптические). Так что часто анизотропия будет, но без перекрёстных коэффициентов.

Впрочем, это я фантазирую.

[svv]

Допустим даже, коэффициент тензорный, но от координат (декартовых) не зависит. Можно было бы поворотом системы координат привести его к диагональному виду. Но если прямоугольная область задачи имела границы, параллельные осям, а после поворота этого не будет, стоит ли это делать?

[Munin]

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

[SteelRend]

А ещё возник такой вопрос (жаль, что я не просчитал эту проблему заранее). В той схеме, что я написал в заглавном сообщении, имеются значения $a^{(1)}_{i+1,j}, a^{(1)}_{i,j}$ ($a^{(2)}_{i,j+1}, a^{(2)}_{i,j}$). При аппроксимации граничных условий через фиктивные узлы будем иметь дело с $a^{(1)}_{0,j}, a^{(1)}_{N_1,j}$ ($a^{(2)}_{i,0}, a^{(1)}_{i,N_2}$), и если $a^{(1)}_{i,j} = k_1((i-\frac{1}{2})h_1,j h_2)$, то нам понадобятся значения $k_1(-0,5 h_1, j h_2), k_1((N_1+0,5) h_1, j h_2)$, которые не были заданы, т.к. находятся за границами рассматриваемой области. Как в этом случае следует поступить?

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

[GAA]

Добавил краткую аннотацию.

SteelRend в сообщении #1206130 писал(а):

Верно?

У Вас переменные коэффициенты. Если нужно, чтобы сетка не выходила за границы области, то в случае постоянных коэффициентов $k_1$ и $k_2$ легко получить и так (и это изложено в учебниках).

SteelRend в сообщении #1206375 писал(а):

Получается ... ?

В случае переменных коэффициентов можно начать с попытки аппроксимации производной с использованием самого уравнения (хоть Тейлор тут и не к месту, но перед чтением теории пользительно). Пусть в одномерном случае (тут всё просто)
$ u_t = (ku_x)_x$ (1)
$\frac{u_1-u_0} h = u_x(0,t) +\frac h 2 u_{xx}(0,t)+O(h^2)$ (2)
Теперь, используя (1) в (2), пытаемся получить аппроксимацию $u_x(0,t)$ со вторым порядком по $h$ и первым порядком по $\tau$ ($\tau$ — шаг по времени). Здесь $u_t$ — производная по $t$, $u_x$ — производная по $x$.

Более простого способа я не знаю. Может другие участники подскажут.

-- Вт 04.04.2017 09:02:38 --

(Сразу забыл написать.) Шаг постоянный, узел 0 лежит на левой границе (как в начальном сообщении ветки).

-- Вт 04.04.2017 09:25:11 --

И ещё. Если коэффициент $k$ не постоянный, то условие второго рода, пусть на левой границе, обычно записывают в виде $k(0)u_x(0,t)= g(t)$. Тогда левая часть (2) изменяется соответствующим образом и может будет проще.

[GAA]

Вечером не будет времени. Поэтому напишу подробней (хоть и банальные вещи).
(Если условие не на поток, а на производную, то умножим обе части граничного условия на $k(0)$ и получим в левой части граничного условия поток. Конечно, предполагается, что $k(0) > 0$.)

1.

$$k(h/2)\frac {u_1^{n+1}-u_0^{n+1}} h \equiv (k(0)+k'(0)\frac h 2 + O(h^2))(u_x(0, t)+u_{xx}(0,t) \frac h 2 + O(h^2)) \equiv $$

$$\equiv k(0) u_x(0, t) + [k(0) u_{xx}(0,t)+ k'(0) u_x(0, t)] \frac h 2 + O(h^2) \equiv k(0) u_x(0, t) + (k(0) u_x(0, t))_x \frac h 2 + O(h^2).$$

Здесь слой $n+1$ соотвнетсвует моменту времени $t$ (при постоянном шаге по времени $t = (n+1)\tau$).
2. Используя уравнение, заменяем $(k(0) u_x(0, t))_x$ на $u_t(0, t)$.
3. Для $u_t(0,t)$ используем $\frac {u_0^{n+1}- u_0^n} \tau$.

Таким образом, на левой границе

$$k(h/2)\frac {u_1^{n+1}-u_0^{n+1}} h - \frac {u_0^{n+1}- u_0^n} \tau \frac h 2 = k(0) g_1(t).$$

[ewert]

Зачем всё это?... $u'_x(0)=\frac{-3u_0+4u_1-u_2}{2h}+O(h^2)$, вот и всё.

[SteelRend]

Спасибо за вашу помощь.
Мне не понятно, почему $k(0)$ мы заменяем выражением $k(\frac {h}{2})$, ведь если принять $a_i = k_{i-0,5} = k_i + O(h^{2}), \frac{a_{i+1}-a_i}{h}=\frac{k_{i+1/2}-k_{i-1/2}}{h} = k'_i+O(h^2)$, то получается, что мы аппроксимируем не $k(0)$, а $k(h)$...

ewert в сообщении #1206467 писал(а):

Зачем всё это?... $u'_x(0)=\frac{-3u_0+4u_1-u_2}{2h}+O(h^2)$, вот и всё.

мне учитель сказал, что так делать НЕЛЬЗЯ (хотя я читал, что можно), мы должны учитывать вид уравнения. Да и разве получится тогда использовать прогонку в методе переменных направлений?

[GAA]

ewert в сообщении #1206467 писал(а):

Зачем всё это?...

Из-за очевидности. И сразу имеем нужную аппроксимацию гран условия.

-- Вт 04.04.2017 14:58:51 --

SteelRend в сообщении #1206468 писал(а):

Спасибо за вашу помощь.
Мне не понятно, почему $k(0)$ мы заменяем выражением $k(\frac {h}{2})$, ведь если принять $a_i = k_{i-0,5} = k_i + O(h^{2}), \frac{a_{i+1}-a_i}{h}=\frac{k_{i+1/2}-k_{i-1/2}}{h} = k'_i+O(h^2)$, то получается, что мы аппроксимируем не $k(0)$, а $k(h)$...

Мы не аппроксимируем $k(0)$. Мы аппроксимируем поток.

[SteelRend]

Вернее, я был неправ.

$\frac{a_{i+1}+a_i}{2}=\frac{k_{i+1/2} + k_{i-1/2}}{2}=k_i+O(h^2)$, т.е. $k(0) \approx \frac{k_{1/2} + k_{-1/2}}{2}$, и тогда опять получается, что нужно знать $k(-0,5 h)$...

-- 04.04.2017, 17:05 --

GAA в сообщении #1206469 писал(а):

Мы не аппроксимируем $k(0)$. Мы аппроксимируем поток.

Но ведь поток $k(x)u_x(x,t)$ берётся в точке $x=0$, и поэтому сомножители мы аппроксимируем в этой точке, разве нет? Почему всё же там записано $k(\frac{h}{2})$? В Самарском тоже подобное видел, но не могу понять, почему аппроксимация берётся при $x=h$, а не при $x=0$.

[GAA]

Ещё раз. Мы не аппроксимируем каждый из сомножителей. Мы аппроксимируем поток. Я привел детальные наивные выкладки выше. Не хотите $k(h/2)$, возьмите $k(0)$, но не берите полусумму и посмотрите, что получится в рамках наивного подхода.

У Самарского в «Теории разностных схем» можно найти подробную теорию.
На всякий случай отмечу, что Вы ещё не указали ничего о $k(x)$, а я производными балуюсь.

[GAA]

SteelRend, может картинка облегчит понимание Самарского.
(На самом деле выписываются интегралы, но для простоты без них. Детали и строгость можно найти в книжках.)
Идея получения граничного условия. (За знаками потоков я не следил. Нужно перепроверить.)

Пусть для простоты рисунка $N=3$. На обеих границах заданы условия второго рода, т.е. заданы потоки. Тело разбиваем на 4 ячейки. Две ячейки половинной ширины — граничные.

Вложение:

BC.PNG [ 3.55 Кб | Просмотров: 734 ]

Запишем баланс для левой крайней ячейки (для правой аналогично).
Изменение количества тепла (теплоёмкость равна 1) равна объёму на изменение температуры. $S \frac h 2 (u_0^{n+1}- u_0^{n})$, где $S$ площадь "бокового" сечения ячейки. С другой стороны это изменение равно плотности потоков умножить на площадь боковой поверхности и на промежуток времени $\tau$. Учитывая, что поток на левой боковой поверхности есть $k(0)g_1(t)$, а на правой боковой $k(h/2)\frac {u_1^{n+1} - u_0^{n+1}} h$ и сокращая на $S$, получим

$$\frac h 2 (u_0^{n+1}-u_0^n) = \left[ k(h/2) \frac {u_1^{n+1} - u_0^{n+1}} h - g_1k(0)\right]\tau.$$

Как-то так в нулевом приближении.

[SteelRend]

Но если $k(x)$ коэффициент теплопроводности, тогда граничное условие второго рода нужно записать как $k(0)\frac{\partial u(0,t)}{\partial x}=g_1(t)$, тогда в скобке справа нужно записать $g_1(t)$ вместо $k(0)g_1(t)$, поскольку именно этой величие в соотв. с граничным условием равен тепловой поток на границе?

Возможно не по теме, но почему теплоёмкость равна 1 и куда в этом случае "девается" плотность $\rho$?

[GAA]

Да. Если $k(0)\frac{\partial u(0,t)}{\partial x}=g_1(t)$, тогда в скобке справа нужно записать $g_1(t)$. Но в начальном сообщении было условие $\frac{\partial u(0,t)}{\partial x}=g_1(t)$.