Одномерное уравнение теплопроводности. Численное решение.

[brute11k]

Доброго дня суток!

Хотелось бы подробно описать интересующую меня проблему и найти такой же подробный ответ на него. Некоторые вещи мне удалось понять самостоятельно благодаря большому количеству имеющейся в моем запасе литературы:

Цитата:

01. Barth T.Numerical Methods for Conservation Law on Structured and Unstructured Meshes 2003
02. Chung T.J. Computational fluid dynamics (CUP, 2002)(T)(1022s)
03. Ferziger J. H., Peric M. Computional Mathods for Fluid Dynamics 2002
04. John D. Anderson, Jr. Computional Fluid Dynamics 1995
05. Андерсон Д. Таннехел Д. Плетчер Р. Вычислительная гидромеханика и теплообмен-Т2-1990
06. Андерсон Д. Таннехил Д. Плетчер Р. Вычислительная гидромеханика и теплообмен-Т1-1990
07. Белоцерковский О.М. - Численное моделирование в механике сплошных сред - 1994
08. Куликовский А.Г., Погорелов Н.В., Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений Физматлит, 2001
09. Патанкар С В Численное Решение Задач Теплопроводности И Конвективного Теплообмена При Течении В К
10. Патанкар С В Численные Методы Решения Задач Теплообмена И Динамики Жидкости 1984
11. Полежаев В И Бунэ А В Верезуб Н А И Др Математическое Моделирование Конвективного Тепломассообм
12. Роуч П. Вычислительная Гидродинамика 1980
13. Самарский А.А., Попов Ю.П. - Разностные методы решения задач газовой динамики - 1992
14. Самарский А А Попов Ю П Разностные Методы Решения Задач Газовой Динамики 1992
15. Самарский А.А. - Теория разностных схем - 1989
16. Р.Рихтмайер, К.Мортон Разностные методы решения краевых задач
17. Л.И.Турчак Основы численных методов - 1987

Теперь, непосредственно само задание(ссылка на скан листа, если что-то не понятно), которое я перепишу в формате $\TeX$

http://xmages.net/storage/10/1/0/e/0/up ... ce6cb.jpeg

Здесь можно увидеть в самом верху, в рамке находится уравнение теплопроводности: $\frac{du}{dt}=\frac{d}{dx}(\kappa\frac{du}{dx})$
далее приводится рисунок, отображающий начальное условие системы при t=0:



Заранее оговаривается метод, с помощью которого мы будем численно решать это уравнение:
Сеточные методы на Эйлеровых схемах (явная схема).

Конечно, нужно будет определить устойчивость схемы и выяснить, при каких критических значениях $\Delta{t}$ эта схема будет сходиться(ответ уже найден и равен $\tau\le\frac{h^2}{2\kappa}$, где $\tau$ - шаг по времени($\Delta{t}$), $h$ - шаг по пространственной координате)

Далее приводится 3 уравнения, назначение которых мне не совсем ясно, но имеются догадки:

$1)\kappa = \kappa_1+bx; 2)\kappa = \kappa_0{*}e^{-\frac{(x-x_0)^2}{r_0^2}}; 3)\kappa = \kappa_0+du;$

Видимо, нужно еще варьировать коэффициент $\kappa$ - коэффициент теплопроводности в зависимости от этой системы уравнений, описанной выше.

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

Мои соображения:

Во-первых, нужно задать граничные условия для нашего уравнения теплопроводности(исходя из рисунка):
Допустим $U(x,0)=U_0,$ где $a{}\le{}x{}\le{}b$ (здесь подразумевается, что a и b - концы стержня)

Во-вторых, нужно ввести пространственно-временную сетку, которая будет представлять из себя матрицу, причем каждая строка её[матрицы] характеризует значения уравнения при каком-либо времени $t$. (здесь поправьте меня, если я не правильно понял смысл всех этих сеток)

В-третьих, мы переписываем уравнение теплопроводности в виде разностной схемы: $\frac{U_i^{n+1}-U_i^{n}}{\tau}=\kappa\frac{U_{i+1}^n-2U_i^n+U_{i-1}^n}{h^2},$ откуда мы можем найти $U_i^{n+1},$ при условии, что нам известно начальное распределение температуры на предыдущем временном слое(например на $n=0$).

В-четвертых, следует учитывать еще эти три уравнения, связанные с самим коэффициентом теплопроводности. Таким образом задача усложняется - требуется понять, каким именно образом следует организовать эту связь и как мы можем вообще их применить либо хотя бы выразить. Я предполагаю, что начальный коэффициент теплопроводности (та самая $\kappa_0$) должен задаваться в начале, а затем уже выводятся и находятся остальные значения уравнений.

Возникает некоторое количество логичных вопросов:
1. На графике нарисован "стержень" - с равномерной начальной температурой по всей его длине. Что будет означать профиль графика в последующие моменты времени - нагрев окружающей среды и остужение самого стержня?
2. Достаточно ли я описал задачу(для описания её непосредственно в дискретной форме, в ЭВМ) и нет ли где "подводных" камней, либо неопределенности касаемо условий задачи?
3. Каким образом проще и понятнее ввести дополнительные условия для коэффициента теплопроводности(имеются в виду 3 вышеприведенных уравнения)?
4. Мелкие, но необходимые для целостной картины вопросы, описанные выше, в тексте.
5. Правильность моих рассуждений.

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

[brute11k]

Еще вот возник вопрос по поводу граничных условий: если мы задаем на концах стержня(a,b) температуру U=0, то значит ли это, что она там постоянна, или это где-то в другом месте меняется?

[Munin]

brute11k в сообщении #383865 писал(а):

В-третьих, мы переписываем уравнение теплопроводности в виде разностной схемы

С этим шагом вы поторопились. Он обычно один из центральных. Нужно узнать, какая именно сетка требуется (существует много разных вариантов), и аккуратно преобразовать дифур в нужный вид. У вас, кроме того, в нескольких вариантах $\kappa$ зависит от $x,$ так что сначала надо вынести его из-под производной, и только потом преобразовывать в разностную схему. И в-третьих, вам повезло, что схема у вас явная.

brute11k в сообщении #383865 писал(а):

1. На графике нарисован "стержень" - с равномерной начальной температурой по всей его длине. Что будет означать профиль графика в последующие моменты времени - нагрев окружающей среды и остужение самого стержня?

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

brute11k в сообщении #383865 писал(а):

2. Достаточно ли я описал задачу

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

Самым коварным (или кошмарным) мне видится вариант 3 (кстати, вы при переписывании упустили важную точку, $du$ - это совсем не то же самое, что $d\cdot u$), поскольку если у вас $\kappa$ зависит от $u,$ то само уравнение становится нелинейным, полностью меняется вид уравнения, и способ, как его решать. Даже условие устойчивости может измениться, причём перестанет быть условием, которое можно рассчитать заранее, а потребует пересчёта в ходе решения.

[brute11k]

Неужели всё настолько сложно в моём случае? Мне всего лишь нужно показать, какова устойчивость явной схемы эйлера и, чуть усложняя задачу, зависимость $\kappa{}(x)$ и $\Delta{}t$ (на самом деле, мне придется еще раз уточнить задание, ибо я не понял, причем тут эта зависимость[$\kappa{}(x)$ и $\Delta{}t$], которая всё только усложняет).

(Оффтоп)

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

Munin в сообщении #384998 писал(а):

С этим шагом вы поторопились. Он обычно один из центральных. Нужно узнать, какая именно сетка требуется (существует много разных вариантов), и аккуратно преобразовать дифур в нужный вид. У вас, кроме того, в нескольких вариантах $\kappa$ зависит от $x,$ так что сначала надо вынести его из-под производной, и только потом преобразовывать в разностную схему. И в-третьих, вам повезло, что схема у вас явная.

По поводу выбора сетки - я думаю можно для начала выбрать самую простую, статическую сетку(есть еще лагранжевая - она нам не нужна здесь), а шаблон(если я правильно понял, это и есть разностная схема) взять самый простой:


По поводу "$\kappa$ зависит от $x$" - выносить из под производной нужно в начальном уравнении ($\frac{du}{dt}=\frac{d}{dx}(\kappa\frac{du}{dx})$)? Если так, тогда что мешает нам переписать его в виде:
$\frac{du}{dt}=\kappa^2\frac{d^2u}{dx^2}$ - получается, что коэффициент $\kappa$ вынесен за производную. Ведь это нам и нужно? Разностная схема по-идее остается та же самая: $\frac{U_i^{n+1}-U_i^{n}}{\tau}=\kappa\frac{U_{i+1}^n-2U_i^n+U_{i-1}^n}{h^2}$ (или я не прав?)

Munin в сообщении #384998 писал(а):

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

Ясно. [a;b] - длина стержня, а сам стержень представлен как бесконечно тонкий, одномерный кусок "провода", а величина $l$ - показывает какая часть(длина) стержня в данный момент имеет температуру $U_0$. Профили в последующие моменты времени покажут распределение температуры по нему[стержню] - причем ясно, что со временем она будет убывать(ведь в нашей задаче нет источника тепла - только начальная температура). Для наглядности, что я имею ввиду, покажу картинку:

Munin в сообщении #384998 писал(а):

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

Мне кажется, граничные условия - это как раз таки и есть длина стержня [a;b] (то есть мы ищем решение не на всей числовой прямой $-\infty\le{}x\le\infty$, а лишь в пределах $a\le{}x\le{}b$). А начальные условия - указаны на графике(границы стержня [a;b], область, где температура отлична от нуля). Даже если это не так, и на самом деле они задаются по-другому, то что мешает нам взять именно такие начальные данные(например, в учебных и самых простых примерах)?
Таким образом, наша начальная сетка(в программе, соответственно двумерный массив) будет заполнена начальными значениями профиля. Последующие вычисления производятся уже с помощью этой самой разностной схемы, описанной выше.

Munin в сообщении #384998 писал(а):

(кстати, вы при переписывании упустили важную точку, $du$ - это совсем не то же самое, что $d\cdot u$)

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

Munin в сообщении #384998 писал(а):

поскольку если у вас $\kappa$ зависит от $u,$ то само уравнение становится нелинейным, полностью меняется вид уравнения, и способ, как его решать. Даже условие устойчивости может измениться, причём перестанет быть условием, которое можно рассчитать заранее, а потребует пересчёта в ходе решения.

Тогда, может, стоит сначала рассмотреть задачу без учёта последнего, третьего уравнения? Для наглядности, простоты и просто тренировки. А затем, когда всё будет понятно, добавить это условие и рассматривать задачу, уже отталкиваясь от накопленного опыта.

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

[Munin]

brute11k в сообщении #385121 писал(а):

Неужели всё настолько сложно в моём случае?

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

brute11k в сообщении #385121 писал(а):

По поводу "$\kappa$ зависит от $x$" - выносить из под производной нужно в начальном уравнении ($\frac{du}{dt}=\frac{d}{dx}(\kappa\frac{du}{dx})$)? Если так, тогда что мешает нам переписать его в виде:
$\frac{du}{dt}=\kappa^2\frac{d^2u}{dx^2}$ - получается, что коэффициент $\kappa$ вынесен за производную. Ведь это нам и нужно?

Если $\kappa=\mathrm{const},$ ничего не мешает (только квадрата там не будет, разумеется). Но как раз в ваших вариантах 1, 2, 3 это не так. Надо честно дифференцировать.

brute11k в сообщении #385121 писал(а):

Профили в последующие моменты времени покажут распределение температуры по нему[стержню] - причем ясно, что со временем она будет убывать(ведь в нашей задаче нет источника тепла - только начальная температура).

Не убывать, а разглаживаться (там где был нуль, возрастает). Остальное всё верно.

brute11k в сообщении #385121 писал(а):

Мне кажется, граничные условия - это как раз таки и есть длина стержня [a;b] (то есть мы ищем решение не на всей числовой прямой $-\infty\le{}x\le\infty$, а лишь в пределах $a\le{}x\le{}b$).

Нет. Граничные условия - это не просто длина стержня, а те дополнительные уравнения, которые задаются на концах стержня, то есть в точках $a$ и $b.$ Эти уравнения физически могут отвечать либо теплоизолированному концу стержня, либо тепловому контакту с некоторым резервуаром с фиксированной температурой, либо ещё какой-то более сложной ситуации. Задавать их обязательно, поскольку без них задача не будет иметь решения (точнее, будет иметь бесконечно много решений, а из них надо выбрать одно и только одно). Существуют правила, какие комбинации начальных и граничных условий соответствуют существованию и единственности решения - они называются корректными постановками краевых задач. Короче, имейте в виду, что без граничных условий ваша задача не имеет решения.

С начальными условиями вы разобрались правильно.

brute11k в сообщении #385121 писал(а):

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

Никакого дифференциала нет. У вас на бумажке написано $\kappa=\kappa_0+d\cdot u$ - здесь $d$ - просто число. А вот вы, переписывая, написали $du,$ что как раз читается как "дифференциал $u$", и в данном контексте нелепо.

brute11k в сообщении #385121 писал(а):

Тогда, может, стоит сначала рассмотреть задачу без учёта последнего, третьего уравнения? Для наглядности, простоты и просто тренировки.

Тогда вообще решите её для случая $\kappa=\mathrm{const},$ добейтесь, чтобы программа выводила картинку, похожую на ту, которую вы привели, с разглаживающимся профилем. А потом переделывайте её под заданные варианты.

brute11k в сообщении #385121 писал(а):

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

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

[ewert]

brute11k в сообщении #383865 писал(а):

В-третьих, мы переписываем уравнение теплопроводности в виде разностной схемы:

Это правда, но только если $\kappa$ не зависит от $x$. А если зависит, то правую часть надо переписать в виде

$\dfrac{1}{h^2}\big(k_{i+{1\over2}}(U_{i+1}-U_{i})-k_{i-{1\over2}}(U_{i}-U_{i-1})\big)$

(где, естественно, $\kappa_{i\pm{1\over2}}\equiv\kappa(x_i\pm{h\over2})$), что общих свойств схемы принципиально не меняет. Только тогда в Вашем замечательном условии устойчивости в знаменатель надо ставить, конечно, максимум $\kappa(x)$ по всем иксам.

[Munin]

ewert
А что вы скажете о варианте $\kappa=\kappa(u)$?

[ewert]

Ничего не скажу (насчёт устойчивости) -- я с нелинейными уравнениями дела не имел.

[egor20]

Там где у вас во первых указано не граничное, а начальное условие

[ewert]

egor20 в сообщении #385773 писал(а):

Там где у вас во первых указано не граничное, а начальное условие

Не понял вопроса. Но в любом случае он явно не по делу: граничные и начальное условия ставятся независимо друг от друга, и все необходимы.

[mitia87]

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

[eugrita]

Выскажу еще свой взгляд на все это. хоть я и не узкий специалист в этой области
Обращаю внимание на следующее.
1)нигде в литературе я толком не нашел ответа на следующий вопрос : оценка времени переходного процесса, т.е перехода нестационарного температурного поля в стационарное. (даже для 1-мерного)
2)интересно было бы на одной модели сравнить динамику температурных профилей при разных значениях коэф-та температуропроводности a , а также сделать выводы опять же о скорости затухания переходного процесса для постоянных граничных условий, а также "инерционности" профилей температуры для нестационарных граничных условий, например вида $Tk=A+Bsin(wt)$
где скажем w физически соответствует суточным колебаниям температуры
Тут действительно уместно бы применить прием компьютерной анимации
3)неплохо бы на модели решать достаточно практические вопросы скажем 2- или 3-слойной теплоизоляции стен. В этом случае стыкуются 2-3 участка с разными a
4)случай нелинейного уравнения конечно интересен, но я не знаю практических его применений,
пусть другие подскажут

[Munin]

eugrita в сообщении #531741 писал(а):

1)нигде в литературе я толком не нашел ответа на следующий вопрос : оценка времени переходного процесса, т.е перехода нестационарного температурного поля в стационарное. (даже для 1-мерного)

Если использовать метод Фурье, то там в явном виде появляются экспоненты $e^{-t/\tau}.$ Это и даст вам искомую оценку. Метод Фурье изложен во всех учебниках по УМФ.

eugrita в сообщении #531741 писал(а):

2)интересно было бы на одной модели сравнить динамику температурных профилей

Это всё студенческие упражнения. Большую часть можно прикинуть аналитически.

eugrita в сообщении #531741 писал(а):

3)неплохо бы на модели решать достаточно практические вопросы скажем 2- или 3-слойной теплоизоляции стен. В этом случае стыкуются 2-3 участка с разными a

Для практических вопросов важны всякие щели, а это трёхмерка, да и параметры не табличные, а как герметик на душу лёг. Оно вам надо?

[Maik2013]

brute11k

У Вас все получилось?