Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 64
Текст из файла (страница 64)
В точке х = 0: а) функция и(г, х) непрерывна, т.е. и(г, — 0) = и(г, + 0); б) непрерывен тепловой поток, т.е. к,их(г, — 0) = кзих(г, + 0). Введем (временно) в точке х = 0 температуру и и запишем разностную аппроксимацию условия непрерывности теплового потока (справа и слева от разрыва функция и гладкая, только точка х = 0 является точкой нарушения гладкости) в виде »» х-ш "пг хо к~ — Кд — = хг — д7з- (Временной индекс не пишем, он может быть и л, и и+1.) Отсюда Опуская простые выкладки, вычисляем П»: И», ьз~ П =2(и — и )! — + — .
о пз -Пз/~х х~' 1 2 (б) Такую аппроксимацию теплового потока в точке контактного разрыва иногда называют «наилучшей». (Ниже будет показано, как опасно придавать этому термину слишком универсальное значение.) Эта же формула применяется не только в случае разрыва коэффициента х, но и при переменном коэффициенте теплопроводности. В частности, в рассмотренной выше задаче «газодинамика + теплопроводность», коэффициент х зависит от термодинамическнх величин и, р (р — плотность вещества), а эти величины определены в полуцелых точках. Таким образом, в каждой «целой» точке т теп- ловой поток аппроксимнруется формулой П =2(и +, — и», )г ~ +и'+ х,„+пз х (7) Аппроксимация при расчете тепловой волны. Выведенная выше «нанлучшая» формула (7), однако, не пригодна для расчета тепловой волны.
И вот почему. Пусть в начальных данных фронт тепловой волны находится в точке хп т.е. и»из = 0 при ш» 1. В этом случае П = 0 при гл > 1 (так как й +и lи«»из= 1/О). Таким обра- НРНБли!Еенные метОды ВычислительнОЙ Физики 1ч. и 3!4 зом, тепловой поток через точку х, равен нулю и температура при х и х, останется нулевой и в дальнейшем, т,е. фронт тепловой волны не продвигается, а застревает в начальном положении. Для правильного расчета тепловой волны используегсм аппроксимация, учитывающая структуру функции и(1, х) в окрестности фронта тепловой волны (3).
Обратим внимание на то, что и«(1, х) — линейная функция х в окрестности фронта. Учнтывам этот факт, в расчетах применяют аппроксимацию типа ь Л " -ш+ " ~!а ~ +!и " -!а и„— (8) Длм линейной функции и«линейная интерполяцим является точной. В свое время автору приходилось решать задачи гидродннамнки с теплопроводностью в ситуацни, когда коэффициент теплопроводности и имел внд х(и, р, х) = /(х, р)и», причем функция ~ была разрывной по х; искомая функцнм р(г, х) тоже была разрывной (контактные разрывы).
Прн этом были плохи обе формулы (7) и (8): первая препятствовала правильному расчету тепловой волны (а это явление играло важную роль в проводившихсм расчетах), вторая приводила к погрешностям на контактных разрывах. Решение было найдено в виде компромиссной формулы П "!а и" +и" « -и (р) 2 05ДЛ У)„, !и+(А У)„,«!Н1' в которой разрывная часть коэффициента теплопроводностн учитывалась так, как это рекомендуется теорией для разрывного коэффициента, а множитель и", ответственный за фронт тепловой волн»!, усреднмлся с учетом типичного )рафика и(г, х) в окрестности фронта.
Уравнение теплопроводностн с нелинейным источником. Рассмотрим уравнение (1) в случае, если х = н(и), О= Д(и), Допустим, мы используем неявную схему. Возникает вопрос: что делать с нелинейносгями в н и Д? Есть два варианта. Можно оставить их «на ин:кием слое» и получить простую схему «+! а ««! +1 ««! ~+!1 где и" +пз — — н(и«+ и" +!)/2. Уравнение «на верхнем слое» (для и"+') линейное; оно решается прогонкой. Второй вариант (нелинейность «на верхнем слое») отличается от схемы (10) только тем, что в нем используются значения Д(и" +').
В этом случае уравнения «на верхнем слое» нелинейные. Их прихо- 3гз нзлинкйног.лАвнвнив тзплопговодности 12П дится решать итерациями с линеаризацией по методу Нъютона и прогонкой для линеаризованных уравнений. Это, конечно, намного сложнее, чем при аппроксимации (10), Из общих соображений трудно понять, зачем нужна такая трудно реализуемая схема. Однако в литературе часто встречаются указания на предпочтительность именно более сложной схемы. В чем дело? Попробуем немного прояснить этот вопрос.
Все дело в характере нелинейности и в шаге т по времени. Грубо говоря, дело обстоит так, Если в рассчитываемом процессе шаг т таков, что [тД„(и) [ ~1, то обе схемы.более или менее равносильны, н следует отдать предпочтение более простой схеме (10).
Поясним это положение следующими оценками. За один шаг т температура изменится на тД (предполагаем, что хи,„— величина того же порядка), т.е. и«»' ж и«+ тД. Вычислим Д(п«+1) ж Д(л«+ тД) Д(в«) ч тД Д (1 1 тД )Д(и«) Если [тЦ„[ ~1, то Д(и«) и Я(и«+') почти совпадают и схема, в которой Д вычисляется на верхнем слое, мало чем отличается от схемы с вычислением Я на нижнем слое. Но бывают задачи, например связанные с расчетом тепловых явлений в звездах и близких к ннм обьектах, когда вычисления с шагом т, таким, что т[Д„[«м1, немыслимы.
Это слишком малый шаг. С таким шагом т за приемлемое время работы ЭВМ не удается провести расчет на заданном интервале времени [О, Т[ (число шагов Т/т слишком'велико) и нужно считать с шагом т:Ф.1/[ Д„[. При определенных условиях выходом нз положения является аппроксимация источника на верхнем слое. Нужно сказать, что эта ситуация внешне и по существу очень близка к тому кругу вопросов, которые мы обсуждали при описании жестких систем и методов их интегрирования (там тоже решающую роль играли неявные схемы). Естественным условием применимости схемы с болъшим шагом т является «усгойчизость»; Д„< О. Описанная выше ситуация часто встречается в задачах астрофизики, когда Я= (А — Д, где Ц, н Д определяют выделение (за счет ядерных реакций,.например) и поглощение энергии соответсгвенно.
Оба этих процесса очень интенсивны и «почти сбалансированы», т.е. [Д, — Д [ к [ф[+ [Щ. Другими словами, выделившаяся в какой-ю точке энергия поглощается почти в этом же месте. Разумеется, термин «почти в этом же месте» означает, что энергия поглощается на расстоянии от места выделения, менъшем шага счетной сетки. В задачах, связанных с расчетом процессов в звездах, когда по радиусу звезды вводится 1Оз «- 105 точек, шаг Ь достигает тысяч километров. ПРИБЛИЖЕННЫЕ МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ ФИЗИКИ 1ч, и 3!6 Если в описываемой ситуации Д„> О, нужно считать с шагом Ежа„! или придумывать что-то другое, более сложное, чем переход к неявной схеме. Разложим решение в ряд Фурье: и" = ~х, с~ з1п (Лингй). е Тогда для коэффициентов Фурье се«без труда можно получить разностиое уравнение: +! 4 г ЛЛЛ «+! = — — ып — е" Лг 2 В ! -« т.е.
с„= с„! 1 + 4 — з1п — ) В/ т . глнл~ Лг 2 ) Точное решение дифференциальной задачи имеет вид и(г, х) = ~ с'!е л«! Егп (Лхх), Е О точности разностиого решения можно судить, сравнивая величины е «"' и !1+4 — з!и — ) . При )Лхй~ !С1 можно положить -Л / Т 2~"! Л 2) 4 т згп2 х 1 4 4 т ЛгигЛЗФ« 1 4- ТЛгиг Лг Зачем нужны неявные схемы? Ответ на этот вопрос кажется совершенно очевидным.
Неявные схемы нужны для того, чтобы считать задачи с шагом по времени т, существенна ббльшим, нем это позволяет известное условие Куранта. Для уравнения и, = и„„, Оих ч 1, /~0, и(2,0) == и(г, 1) =О, на примере которого мы ниже рассмотрим некоторые вопросы, условие Куранта, как известно, имеет вид т ~ 0.5Лг. Итак, неявные схемы позволюот проводить устойчивый счет при т м Лг/2, например при т = Л. Но всегда ли следует пользоваться этим преимуществом? Обсудим это.
Результат будет примерно такой. Для уравнения теплопроводности соотношение т ж Лг, в известной мере, естественное. Его не следует нарушать очень уж сильно, счет с т м» Лг требует большой осторожности. Это связано с другим важным понятием — с фактической погрешностью аппроксимации. Перейдем к конкретному анализу. Допустим, мы проводим расчет с шагом Л = 1/100. Имеет ли смысл расчет с шагом т = Л = 1/100? Вообще говоря, иет. Используем неявную схему: «Ф! ««+! «+!+ «+! — лг=1,2,... М вЂ” 1. Т Лг нклннкйнок !тАвнкннк ткплопговодностн в 211 Если еще принягь, что г/!гхг~ 1, то 1 + гхгяг м е'к и и 1+4 — зш — же < г г йий'! Зг 2! Этой простой оценкой мы выяснили следующие обстоятельства.
Достаточно правильно вычисляются те Фурье-компоненты решения (т.е. те члены ряда х~; с" з!и (кпн!/!)), у которых Ы/1 !с1. Так как /! = 0.01, то это, грубо говоря, первые 20 гармоник (20п /г 0.6). В общем случае, при Ь 1/К, это примерно 1/5 всех присутствующих в разностном решении гармоник, так как (11//5)хл ~ 0.6 . Считаем, что О.б«м1, тзк как з1пга/аг 0.89 при а = 0.6.