XIII Власова Б.А., Зарубин B.C., Кувыркин Г.Н. Приближенные методы математической физики (1081417), страница 64
Текст из файла (страница 64)
Разностная схема (8.71) и ее шаблон являются шеститочечными. Для одномерной равномерной сетки он изображен на рис. 8.4. 464 8. ОДНОМЕРНЫЕ КРАЕВЫЕ ЗАДАЧИ Процесс вычислений по (8.72) будет устойчивым, т.е. погрешность вычисления температур Т„", п = О, А', по абсолютной величине не будет превышать наибольшей из погрешностей, возникших при вычислениях на предыдущем интервале времени, если сумма модулей коэффициентов при значениях температуры в правой части (8.72) не превышает единицы: 1 — — Ь1~)+ Ь1ь < 1, п = О, Ж. 1- Ь„~ (и„(+ )с„( л» л» Отсюда в силу частичного диагонального преобладания матрицы А в (8.71) и положительности ее диагональных элементов имеем )а„) + Ь„+ )с„! 2Ь„ Лавашах " " "< шах — "<2, »ыо, Н в» »=О,Н Л» что не противоречит условию (8.62), но приводит в общем случае к несколько менее жесткому ограничению л» Ыь < пп'и— »»0,1Ч Ь» (8.73) выбора интервала времени Ь1ь.
Если в (8.71) принять о = 1, то получим неявную двухслойную разностпную схеляу в виде СЛАУ (Б+ АЬЬ,)ть =БТь ~ + ~ьД1ь (8.74) с трехдиагональной симметрической матрицей 5, = Я+ АЫь имеющей положительные диагональные элементы и обладающей частичным диагональным преобладанием. Четырехточечный шаблон этой схемы для одномерной равномерной сетки изображен на рис. 8.5. Иэ (8.74) не удается получить явного выражения вида (8.72) для вычисления узловых значений Т», но такую СЛАУ можно решить обычным неоюдо и нроеонки, причем разностинал схеми (8.74) в сочетании с алгоритмом метода прогонки будет корректной, Отметим, что при Ь1ь = Ь1= = сопяс коэффициенты в„в (8.86) не будут зависеть от номера 8.3.
Неетационарнал задача тенлонроаодноети 4б5 й интервала времени. Это существенно уменьшает число арифметических операций при решении нестационарной краевой задачи. и+1 Рис. в.в Выбор в (8.71) 11 = 1/2 приводит к так называемой двухсдойиой симметричной разностной схеме. Она также является неявной и соответствует СЛАУ (.5'+ А — )Тк = (5 — А — )Ть 1+ (~ь 1+ ~я)— относительно вектора Те узловых значений температуры в конце интервала времени 1а11ь. Трехднагональная симметрическая 1 матрица л.' = о + -А1а1ь этой СЛАУ также имеет положительные диагональные элементы и частичное диагональное преобладание, что обеспечивает корректность разностной схемы при 11 = 1/2 в сочетании с методом прогонки.
Поскольку эта схема эквивалентна (8.65), то она устойчива при любых значениях Ь1ы а погрешность аппроксимации этой схемой производных по времени пропорциональна (Ь1ь)~. Шеститочечный шаблон этой схемы аналогичен изображенному на рис. 8А. Замечание 8.4. Отметим, что формально устойчивая разностная схема может подавлять возникающие вычислительные погрешности двумя путями: затухающими колебаниями и монотонно. Первый путь не отвечает физическому содержанию задачи теплопроводности, так как при отсутствии внешнего теплового воздействия температура любой точки тела стремится к установившемуся значению монотонно, без колебаний 466 8. ОДНОМЕРНЫЕ КРАЕВЫЕ ЗАДА ЧИ около этого значения.
Поэтому для получения при численном решении задачи нестационарной теплопроводности с применением (8.71) физически оправданных результатов следует выбирать 1з1» иэ условия монотонного затухания погрешности. Можно показать', что схема (8.71) в случае равномерной пространственной сетки с шагом Л обладает таким свойством при Л Производную по времени в (8.56) можно аппроксимировать не только при помощи разности узловых значений Т„в начале и конце Й-го интервала времени Ь1», но использовать для этого известные значения Т„в предшествующие моменты времени.
Такал аппроксимация позволяет повысить порядок погрешности по Ь1», но требует в процессе решения задачи хранить в памяти ЭВМ результаты вычислений на нескольких предыдущих интервалах. Если использовать центральную конечную разность и скорость изменения температуры в пределах удвоенного интервала времени 2Ь1=1» — 1» э вычислить в момент времени 1» 1, то из (8.56) с учетом (8.53) получим лвную трехслойную разностную схему Т„" =Т„~+2(а Т» ' — Л„Т~ '+с„Т~+1'+~» ')»ь1. (8.75) Эта схема неустойчива при любых значениях Ь1, так как сумма модулей коэффициентов при значениях температур в правой части (8.75) превышает единицу.
Замена Т" ' -(Т»+Т„" а) в (8.75) приводит к явной трех- 1 2 слойной схеме (1 — 2Б„М) Т» э + 2(а„Т„1 + с„Т„"+1 + )„) Д» 1 + 25„1.'11 устойчивой в силу (8.54) и (8.55) при любых значениях Д»1, поскольку теперь сумма модулей коэффициентов при значениях 'См., например: Зарубин В.С., Селиванов В.В. 467 а4. Некоторые динамические задачн температур в правой части не превышает еднннцы. Однако эта схема прн больших значениях Ь| становится неточной, узловые значения температуры от интервала к интервалу нзменяются немонотонно даже прн отсутствии внешнего теплового воздействия, н результаты расчета теряют физический смысл (см. замечанне 8.4).
С учетом физических соображений выбор больших значеннй Л1 допускает нелвнал трехслойная ризностнол схема, в матрнчной запнсн имеющая внд 2Ьс 3 = — А(Те+Та 1+Та г)+1д 1. (8,77) 8.4. Некоторые динамические задачи Трехслойные разностные схемы прнменяют прв решении методом конечных разностей (МКР) динамических краевых задач, описываемых дифференциальными уравнениямн со второй производной искомой функции по времени 1. Напрнмер, прн распространении в идеальной (невязкой) среде волн малой амплитуды функция и(1,х) перемещения частиц среды удовлетворяет волновому уравнению дзи(1,х) Дзи(1,х) — а, 1>0, хб(01), д1з дхз (8.78) Она удобна, в частности, прн решении нелннейных задач, когда с, д н! е в (8.49) произвольным образом эавнсят от температуры.
В этом случае элементы матрицы А в координаты вектора ~ь 1 вычисляют по известным узловым значениям температуры в момент времени 1ь 1 н, используя метод прогонки, находят узловые значения температуры в момент времени 1е. Однако прв расчете на первом интервале по временн (Й = 1) вектор Тл з нензвестен. Поэтому применение трехслойной разностной схемы прн решении задачи теплопроводностн возможно лишь со второго интервала, а расчет на первом интервале следует провести по какой-либо двухслойной разностной схеме. 468 8. ОДНОМЕРНЫЕ КРАЕВЫЕ ЗАДАЧИ где а — скорость звука в этой среде. Примем для определенности на концах отрезка [0,1) граничные условия в виде и(С,О) = Хо(С) и и(2,!) = Л(2), а начальные условия и(0, х) = и'(х),, = о'(х) ди(1,х) ) ~=о (8.79) согласуем с граничными: и (0) = 7о(0), и (1) = Л(0), о (О) = А( Нс=о' о ( ) = И' Ис=о' (8.80) На равномерной одномерной сетке с шагом 6 и узлами х„= и12, и= О, Х, при постоянном интервале времени Ь1= т = = сопеС можно аппроксимировать (8.78) нзрехслойиой симметричной разностпной схемой + а2г2 62 Л-2 Ь-2 Л-2 Ь-1 Л 2 Ь-1 +ц " +(1 — 20) " " ", (88!) и„, — 2и„+ и„+, и„, — 2и„' + и„ ди(1, х„) и1 — и„' д1, о 2г н=О,Х, где верхним индексом — 1 отмечены узловые значения и„функции и(1,х) на фиктивном слое, соответствующем моменту вре- и = 1, )1' — 1, с весовыми коэффициентами ц и 1 — 221.
Кроме того, из граничных условий имеем ио~ — — ~о(еь) и ин~ — — ~~(4), где Е~, = = кг, а б г(, а из первого условия (8.79) — ио = и'(х„), н = О, Ж. Так как в (8.81) аппроксимация вторых частных производных искомой функции и(1,х), входящих в (8.78), имеет второй порядок погрешности (см 1.2), то целесообразно, чтобы аппроксимация производной в (8.79) имела бы такой же порядок погреШности.
Для этого используем центральную конечную разность и запишем 469 ВЛЬ Некоторые динамические задачи мени ! = — г. Отсюда находим и„' = и„' — 2йт и после подста- новки в (8.81) при и = 1 получаем СЛАУ (2~. ~+ 1) „2 1 = и„+ и„г+ у (и„, — 2и„+ и„+,)— оо1292ооо 2 2( о 2 а+ о ) (882) относительно неизвестных значений и„', и = 1, 1ч' — 1. СЛАУ (8.82) имеет трехдиагональную симметрическую матрону с одинаковыми диагональными элементами 6к =к 2!!т2+ 1, ат где у = — „, а все ненулевые внедиагональные элементы также одинаковы и равны -21у2, т.е. при у > 0 и любых значениях у она является матриией с диагональным преобладанием. Следовательно, СЛАУ (8.82) имеет единственное решение, которое нетрудно получить с учетом известных из граничных условий значений ис!, = Ят) и и!!о —— /!(т) обычным методом прогон!си.
Затем тем же методом можно при и > 2 решить СЛАУ 2 Л 1 (2 2» 1) !с 2 й +(! 2 ) 2( 'с-! 2 ~-!+ ~ — !)+2 ~ — ! ~ — 2 которая следует из (8.81). СЛАУ (8.83) также имеет матрицу с диагональным преобладанием, что обеспечивает вычислительную устойчивость алгоритма метода прогонки. Можно показать*, что разностная схема (8.81) в сочетании 2 ат 1 с методом прогонки при условии т = — < обладает и Лз 1 — 4с! устпойчивостью по входным данным. При 21= 1/4 это условие выполняется для любых значений интервала времени г.
В частном случае у = 0 (8.82) и (8.83) явно разрешимы относительно 'Смэ Самарский А.А. 470 а ОднОмеРные кРАеВые 3АЯАчи искомых узловых значений и~, и = 1, Ф-1, й Е М, но при этом Л возникает ограничение г < — на выбор интервала времени. а При невыполнении равенств (8.80), выражающих согласование начальных и граничных условий, волновое уравнение (8.78) может иметь негладкие решения н(1,х) с разрывами частных производных на линиях х ~ а1 = сопв1, называемых характеристиками волнового уравнения (ХН].
Негладкие решения, имеющие разрывы не только производных, но и самих искомых функций, характерны для задач газовой динамики. Для построения таких решений используют специальные разностные схемы*. Рассмотрим особенности построения разностной схемы для динамической задачи, описываемой системой уравнений с несколькими неизвестными функциями.
В качестве примера обратимся к одномерной математической модели магнитной гидродинамики (см. пример 3.3), содержащей систему нелинейных уравнений др дри ~(п д г ро гл дН диН вЂ” + — =О, р — + — ~'р(р)+ — Н'~ =0, — + — =0 д! дх1 ' с(1 дх1( 2 ) ' дг дх~ с тремя неизвестными функциями: плотностью р, скоростью среды о и напряженностью магнитного поля Н, зависящими от времени 1 и координаты х~ (ро -- магнитная постоянная, зависимость Р(р) давления среды от плотности считаем заданной). Используя первое уравнение, преобразуем первое слагаемое во втором уравнении: с(и до до дро дрог р — =р — +р — = — + —. й д! дх~ д1 дх1 ' т Для сокращения записи введем матрицу-столбец Н = (р ро Н) размера 3 х 1 и представим систему уравнений в матричном 'См., например: Самарский А.А., Попов Ю.П., а также: Пирумов У.Г., Росляков Г.С.
или Риятмайер Р., Мортон К. 471 8.4. Некоторые динамические задачи виде де/ дети дИ' + (8.84) д1 дх! дх! Н' где И' = (О р(р) + /то — 0) — матрица-столбец. 2 Для решения системы (8.84) при помощи МКР можно использовать как явные, так и неявные разностные схемы. Однако в случае неявной схемы неизвестные величины в конце каждого интервала по времени в силу нелинейности системы приходит~я находить последовательными приближениями.