XIII Власова Б.А., Зарубин B.C., Кувыркин Г.Н. Приближенные методы математической физики (1081417), страница 68
Текст из файла (страница 68)
6), а для численного решения применимы метод конечных резкостей (МКР), метод конечных элементов и метод граничных элементов. Можно также в области решения многомерной нестациоиарной задачи ввести пространственную сетку и на этой сетке аппроксимировать производные искомых функций по пространственным координатам. В результате получим систему обыкновенных дифференциальных уравнений (ОДУ) относительно изменяющихся во времени узловых значений этих функций. Такая система в сочетании с заданными начальными условиями составит математическую формулировку задачи Коши, для численного решения которой можно использовать один из вариантов метода Рунге — Кутты.
Как и в случае одномерных иестационариых задач, третий подход к решению многомерных задач объединяет первые два и связан с переходом к дискретной математической модели Рассматриваемого физического процесса как во времени, так и в пространстве. Эта модель на каждом к-м интервале приведет 492 9, МНОГОМЕРНЫЕ ЗАДА ЧИ к системе конечных уравнений относительно узловых значений искомых функций в момент времени ге в конце этого интервала. Если исходная краевая задача является линейной, то конечные уравнения будут также линейными. 9.2.
Двумерная и трехмерная задачи теплопроводности т(г,Р) =~(1,Р), Р~Г; Т(0,М) =Т'(М), М~Р, (9.2) где à — граница области Р (рис. 9.1). Введем в Р двумерную пространственную сегвку Рк с шагами Й;, г' = 1, 2, постоянными вдоль каждой из координатных осей (на рис. 9.2 принято 61 — — вз — — Ь = сопв$). При этом в узды, вообще говоря, не принадлежащие Г, но расположенные наиболее близко к границе и составляющие множество Гь, перенесем из ближайшей к каждому узлу точки Р Е Г заданные граничные значения (9.2) температуры. Отметим, что возникающая при этом погрешность является основной причиной, ограничивающей применение МКР к решению задач в многомерных областях произвольной конфигурации.
Равенством Т„, — 2Т„+ Т„+1 ЛИТ„= а 6, 1 1=1,2 (9.3) Третий подход к решению многомерных задач (см, 9,1) рассмотрим сначала на достаточно простом примере двумерной нестационарной задачи теплопроводности в твердом теле, описываемой уравнением От(г, М) !О'т(г, М) О'т(г, М) 1 =а~ О э + О э, Г>0, Мег', (91) хэ ( где а — коэффициент температуропроводности, а температура Т(г,м) зависит от времени г и декартовых координат яы хз точки М Е Р двумерной области Р решения задачи. Граничные и начзльные условия примем в виде 9.2.
Двумерная н трехмерная задачи тепяопроводноети 493 Рис. 9.1 Рис. 9.2 где и — номер внутреннего узла двумерной сетки Рь, отсчитываемый вдоль координатной оси Ох;, введем размосттамьяб отзеРоттаоР Л„. ТогДа ДлЯ любого внУтРеннего Узла Мо Е Ра зтой сетки, имеющего температуру 7„(1), с учетом (7.3) запишем = Лит„(2)+О(ь;). дх~ м м Опуская нижний индекс в обозначении узлового значения температуры и аппроксимируя в (9.1) производную по времени на й-м интервале, для каждого внутреннего узла сетки К, получаем Т" — Т" = ~) Лп(11;Т +(1 — т1;)Т" '), т1, б [0,1). (9.4) 1м1,2 В более подробной записи (9.4) имеет вид Т „— Т„,„Т -1 о — та+7 +1, — 27е Ь2 =91 ' 62 ь 1 + 1 7~ „1 — 27~„+7~ „+1 Т~~ „' 1 — 27~„'+7~ „'+1 + т12 + (1 тт2) ~2 ~2 а мнОГОЧЕРнЫЕ ЗАДА НИ 494 Шаблон этой разностной схемы изображен на рис. 9.3.
Она является двумерным аналогом двухслойной схемы (8.71) с весами, имеет при т1! — — т1э = 1/2 погрешность 0((Ь4)~, Ь~!+ Ьз~) и приводит к системе линейных алгебраических уравнений (СЛАУ) относительно значений Т" во внутренних узлах сетки Рю Для решения этой СЛАУ необходимо использовать известные значения Т в узлах, принадлежащих множеству Гь, с и значения Т~ ' во всех узлах сетки Еь, вычисленных на предшествующем интервале времени. Ясно, что при й = 1 значения Т~ ' определены начальными условиями (9.2). Т" ил+1 --Й а Т „! т, Т" ' ~в~-3л а+1,и Рис.
9.3 Из (9.4) следует, что квадратная матрица СЛАУ в каждой строке содержит пять ненулевых элементов, причем один из них расположен на диагонали, а расположение остальных зависит от соответствия между нумерацией узлов и уравнений системы. Однако при любом варианте нумерации эту матрицу не удается привести к плтидиагональной, что не позволяет использовать для решения СЛАУ соответствующий вариант метода прогонки (см.
Д.8.1). Решение этой СЛАУ можно получить при помощи матричной прогонки, обеспечивающей заметное снижение общего числа арифметических операций, если область г' является прямоугольной. В случае произвольной области г' одним из экономичных по числу арифметических операций методов решения рассматриваемой СЛАУ является !зродольно-!зо!зереннал прогонки, цг. двумерная и трехмерияа задачи топлопроаодиоохи 495 Она основана на переходе от разностной схемы (9.4) к разност- ной схеме Т" '1г — Тя ' Ыь — гЛЛиТ + (1 — Ъ)Лггт" ', Т" — Т" = (1 — г1г)Л Ть '1г+ г1гЛг Т", Ыь (9.5) В; = 1 — г1;ЛнЬгь, С; = 1+ (1 — г1;)ЛнЬгь, г = 1, 2, (9.6) где 1 — тохсдественный оператор. Тогда (9.5) примет вид В,Т~ ~1~ =СгТ~ В,Т" = С, Т'-'1', откуда получаем ВгВгТ" =С~СгТ" ', т.е.
с учетом (9.6) при- ходим к разностной схеме г Т~ Т (1+ ~,9гЛ„Лгг(Ь1,)') Ьсь = (1 — Ъ вЂ” Чг)Ль ьЛггТ Ь1ь+ + ~~~ Лн(т1;Т +(1 — гу,)Т~ ~). (9.7) Каждое из двух соотношений этой схемы представляет совокупность одномерных неявных двухслойных разностных схем по осям х, и хг соответственно, причем первое из них содержит неизвестные узловые значения Т" '1г температуры в момент времени 1я,1г в середине Й-го интервала, а второе — значения Т" в тех же узлах в момент гя в конце этого интервала. Отдельно взятая одномерная схема, следующая из первого соотношения (9.5), включает значения Т" '1г в цепочке узлов при хг = совая, а следующая из второго соотношения — искомые значения Т" в цепочке узлов прн хг = сопвг. Все эти узловые значения могут быть найдены при помощи обычного ляетода прогонки.
Выясним условия эквивалентности (9.4) и (9.5). Для этого исключим Т" '1г из (9.5), введя разностные операторы 9. МНОГОМЕРНЫЕ ЗАДА ЧИ ЬЯ аЬФь < г э' 1+ 2 Рассмотренные свойства разностной схемы (9.5) можно распространить на случаи переменных (в том числе разрывных) коэффициентов в уравнении теплопроводности, наличия в теле источников энерговыделения и двумерной сетки с переменным шагом в направлении каждой иэ координатных осей, а также — на случаи двумерных задач в цилиндрической системе координат и осесимметричных задач в сферических координатах. Однако идею продольно-поперечной прогонки не удается в общем виде перенести на трехмерные задачи нестационарной теплопроводности.
Но для задачи в произвольной области с размерностью а', описываемой дифференциальным уравнением дТ ~ д'Т д = .,а.,' (9,8) Сравнивая уравнение (9.7) с (9.4), заключаем, что при п1 = = пз = 1/2 рвзностнал схема (9.7), а значит, и (9.5) аппроксимируют дифференциальное уравнение (9.1) с погрешностью 0((Ь1ь)з,6~+ йэ2). Отметим, что несмотРЯ на совпаДение порядков погрешности разностных схем (9.4) и (9.5) результаты расчетов по ним будут различны, что объясняется наличием в (9.7) дополнительных слагаемых. В частном случае и, = п2 — — 1 каждое иэ соотношений (9.5) приводит к неявной двухслойной разноетной схеме, обеспечивающей вычислительную устойчивость алгоритма метода прогонки при любых значениях Ь|ь > О. В этом случае (9.5) называют локально-одномерной разностпной схемой. Ясно, что переход в (9.5) при выборе п1 = пз = 0 к явной схеме не имеет смысла, поскольку такой выбор приводит к явной двухслойной разностной схеме непосредственно из (9.4).
Можно показать, что при этом алгоритм вычислений устойчив при выполнении условия Нг. Двумерная и трехмерная эадачи тендопроводности 497 можно построить локально.одномерную разностную схему ТЯ 1+'1Я вЂ” Тя '+11 %Я = Лп(цТ ''~ч + (1 — т~я)Т +1' у"), (9.9) которая при каждом значении 1'= 1, Ы позволяет использовать обычный метод прогонки в отличие от разностной схемы Т~ — Т" ' = ~~1 Л„(0;Т + (1 — т1;)Т" '), (9.10) 1=1 обобщающей (9.4) и аппроксимирующей (9.8) непосредственно.
При и; = 1/2 для всех значений 1= 1, Н разностпая схема (9.10) аппроксимирует (9.8) со вторым порядком погрешности как по Ь1ь, так и по шагам 6; пространственной сетки, а в остальных случаях порядок погрешности по Ыя уменьшается до единицы. Выясним порядок погрешности схемы (9.9). Для зтого исключим из (9.9) промежуточные значения Т~ 1+'~Я, 1 = = 1, И вЂ” 1, и с учетом (9.6) запишем Отсюда, например, для случая Н = 3 получим ТЯ Ть 1 3 = ~ Лп(~,Т" +(1-Ч;)Т"-')+ ям 1 +Л11ЛггЛзз(11727зТ +Х1Х2ХзТ" ')(ЬА) — (171112Л11Лгг+ Ъ0зЛггЛзз+ ЧзЧ1ЛззЛ11)Т (сясь) + + (Х1ХгЛ11Лгг+ ХгХзЛггЛзз+ ХзХ1 ЛззЛы) Т 1т1ы 498 9.