Вычислительная теплопередача Самарский А.А. Вабищевич П.Н. (1013606), страница 97
Текст из файла (страница 97)
Аппроксимация граничного условия (3) на левой границе дает у(0, хз) = О, 0 < хз < 1,. (13) На нижней и верхней границах (см. (3)) имеем 1 -(а~ух,)„= О, 0<х~ <1, хз=О, (14) 1 -(а,уе,)„= О, 0 < х~ < 1, хз = 1т. При аппроксимации граничного условия (9) получим 1 В1 1 В1 — а,(х)у;, + — у - - (азуе,)„ = †, х~ = 1, О < хз < 1ь (15) Аналогично аппроксимируется уравнение в угловых точках расчетной области й = й о дй. После умножения на произведение шагов сетки Ь,Ьз разностная задача (10) — (15) записывается на пятиточечном шаблоне в следующем симметричном виде -А2ьу ~уд ! — А1г,ду;,д+АООуб — А10уьмд — А20уд+~ — — ф, (16) 1 = 1, 2,..., Ьгп у = 1, 2,..., Ьгн В такой записи предполагается, что любая сеточная функция вне сеточной области обращается в ноль.
Здесь предполагается, что нумерация узлов по каждому направлению начинается с единицы. 1 — з(*)у*,— Ьз 1 — аз(х)у;,— 2 1лава 13. Примеры численного моделирования 13.1. Стационарное телаолроводность в кусочно-однородном теле 647 Я1ВК06Т1ХЕ Р11$0! ( АО, А1, А2, Г ) С С С С С ИЖ01 — Подпрограмма формирования системы линейных уравнений для разиоствой задачи на патиточечном шаблоне. 1МРЕ!С1Т КЕА1 '8 ( А — Н, 0 — Е ) КЕА1'8 К1, К2, КАРРА 0!МЕХЕ!ОХ АО(Х1,Х2), А1(Х1,Х2), А2(Х1,Х2), Г(Х1,Х2) СОММОХ / Т01 / ХНо Х111, Х1К, е Х21., Х2Тг, Х2К, е КАРРА, В1, е Х1, Х2, Н!, Н2 С С Шаги сетки С Н! = (Х1К-ХИ.) / (Х1-!) Н2 = (Х2К вЂ” Х2Е) / (Х2 — 1) Н12 = Н1 / Н2 Н2! = Н2 / Н1 С С Внутренние узлы С 110 20 1 = 2, Х2-1 Х2 = Х21.
+ (1-1)'Н2 ОО 10 1 = 2, Х1-1 Х1 = Х1Е + (1 — !)'Н1 А1(1-1,1) = Н21'К1(Х! — 0.500еН1,Х2) А1(1,1) = Н21'К1(Х1+0.5110'Н1,Х2) А2(1,1-1) = Н12'К2(Х1Х2 — 0.500'Н2) А2(1,1) = Н12'К2(Х1,Х2+0.5110'Н2) АО(1,1) = А1(1-1,1) + А1(1,1) + А2(1,1 — 1) + А2(1,1) Г(1,1) = 0.110 10 СОХТ1Х%3Е 20 СОХПХЦЕ Формирование разностной задачи (массивов в (1б)) проводится в подпрограмме И)301, Назначение и структура втой подпрограммы с учетом приведенного ниже текста в дополнительных комментариях не нуждается. Р!ава 13.
Примеры численного моделирования Левая граница: краевое условие первого рода РО 30 Х = 2, Х2 — 1 А1(13) = О.РО Аг(1,1) = О.РО АО(13) = 1.РО Р(1,У) = О.РО 30 СОМПХЮЕ С С Нижняя к верхняя границы: краевое условие второго рода С РО 40 1 = 2, Х1-! х! = хьь + (1-1) и! А1(1 — 1,1) = 0.5РОеН21чК!(Х! — 0.5РОеН1,Хгь) А1(!,1) = 0.5РО'Н21*К1(Х1+0.5РО'Н1,Х2Ь) АО(!,1) = А1(1-1,1) + А1(1,1) + А2(1,!) В(1,1) = О.РО А1(1-1,Х2) = 0.5РО'Н21'К1(Х1-0.5ВО'Н!,хгВ) А1(1,Х2) = 0.5ВО*Н21'К1(Х1+0.5РО'Н1,хгК) А2(1,Х2) = 0 РО АО(1,Х2) = А1(1-1,Х2) + А!(1,Х2) + А2(1,Х2 — !) В(1,Х2) = О.ВО 40 СОХ'ПСЕ С С Правая граница: краевое условие третьего рода С РО 50 ! = 2, Х2 — 1 хг = хгь + (1-!)*Нг А1(Х1,1) = О.РО А2(Х1,1 — 1) = 0.5РО*Н12'Кг(Х1В,Х2 — 0.5РО'Н2) А2(Х1,1) = 0.5РО*Н!2'Кг(Х1Е,Х2+0.5РО*Н2) АО(Х1,1) = А1(Х1-1,1) + Аг(Х1,1-1) + Аг(Х1,1) Ф + Н2*В! Р(Х1,1) = Н2*В1 50 СОМПХЬЕ С С Левый нижний угол С А1(1,1) = О.РО А2(1,1) = О.РО АО(1,1) = 1.РО В(1,1) = О.РО 13.1.
Стационарная телеопрооодность в кусочно-однородном теле 649 С С Левый верхний угол С А1(!,М2) = О.!70 А2(1,Х2) = О.!70 АО(1,Х2) = 1.00 Р(!,Х2) = О.РО С С Правый нижний угол С А1(Х1,1) = О.РО АО(М1,1) = А1(Х1 — 1,1) + А2(М1,1) + 0.5ВО'Н2'В! г(Х1,1) = 0,5ВО'Н2'В! С С Правый верхний угол С А1(М1,М2) = 0.00 А2(Х1,Х2) = О.!70 АО(Х1,Х2) = А1(Х! — 1,М2) + А2(М1,Х2 — !) + 0.500'Н2'В1 Г(Х1,Х2) = 0.5РО*Н2'В1 13.1.4.
Итерационное решение сеточной задачи Для решения разностной задачи (10) — (15) (индексная запись (1б)) используется попеременно-треугольный метод (см. п.4.7). Запишем разностную задачу в виде (17) Ау = 7, причем (18) где 27 — диагональная, а Ь вЂ” нижне-треугольная матрицы. В попеременно-треугольном итерационном методе оператор перехода с одного итерационного шага на другой выбирается в виде: В = (а+ ь)а-! (а+ ь') (19) где гл — диагональная матрица с положительными элементами. Некото- рые способы выбора б отмечались выше в п. 4.7.
650 Б»ава 13. Примеры численного моделирования При реализации попеременно-треугольного метода удобно преобразовать исходную задачу (17), (18) к задаче, для которой оператор перехода имеет наиболее простой вид. Домножим слева исходное уравнение (17) на С '/» и обозначим У = а'/'у, К= а '"1. Тогда уравнение (17),(18) примет вид Ау =,г, (2!) где теперь (22) (24) Чо = »ео, 9» = ге»+/3»9»-ц /о = 1,2,..., где теперь т» = ~ Ау», ге» = В т». (25) Итерационные параметры а» и /3» определяются соотношениями /5» =, /» =1,2,..., (т», »е») (т»-пи» ~) а»= ', 9=0,1,....
(», ») (Оы Ач») При этом для вычисления т» используется рекуррентная формула т» — — т» ~ — а» ~А9» н й= 1 2,.... (27) Модификация попеременно-треугольного метода — сопряженных градиентов позволяет исключить при реализации умножение матрицы А на вектор и не вычислять»е» в (25)-(27). Принимая во внимание (23), из (25) получим (2б) ге» =У '(1/*) 'т». (28) Для преобразованной задачи (20) — (22) оператор перехода попеременно- треугольного метода В = б '/»ВС '/' записывается следующим образом: В=У*У, У=Е+Х. (23) Вычислительные затраты на обращение этой матрицы В меньше, чем для матрицы исходной задачи (17), (18) (см. (19)).
Для приближенного решения уравнения (21), (22) использовался итерационный метод сопряженных градиентов. При реализации используются (ср. с п. 4,б) расчетные формулы: у»~~ =у»+а»9», /» = О, 1,..., 13.1. Стационарная тенлонроеодноеть е нусочно-однородном теле 651 (31) й=1,2..., 3 Г ) ( 4) й=0,1,... (д», 1») Приведенные расчетные формулы (31)-(34) показывают, что переход с одной итерации на другую осуществляется фактически на основе решения двух систем уравнений с матрицами У и У* (при вычислении Т» согласно (33)) н вычисления двух скалярных произведений при расчете итерационных параметров согласно (34).
В подпрограмме Я01УЕ! реализован вариант попеременно-треугольного итерационного метода — метод приближенной факторизации (см. п.4,7). Алгоритм реализации попеременно-треугольного метода — сопряженных градиентов в соответствии с (30) — (34) состоит в следующем. 1. Вычисление элементов диагональной матрицы С У . 2. Переход от задачи (17), (18) и задаче (21), (22).
Р:=С ~РС ьн — 2Е, С-! /2 ~ ~ С-1/2 С учетом (28) введем новые вектора у» = Уу», г» = (У") 'г» = У»о», д» = Уд», Т» = Я") ~Ад», й = О, 1,. Принимая во внимание (22) н (23), для 1» имеем 1» = (17')-'(х,+Р+Т,*)17 'д, = =(17') '(17'+(Р— 2Е)+17)17 'д» = =У 'д»+Ш') '(д»+(Р— 2Е)(7'д»). (30) С учетом введенных обозначений из (24) умножением слева на У получим у»ч = у» + а»дм й = О, 1,..., до = гь, д» = г» +)3»д»-и й = 1, 2,....
Аналогично, домножая (27) слева на (У') ~ „получим гь=у» ~ — а» А и й=1,2,.... (32) Для вычисления вектора Т» из (30) имеем 1»=Р 'д»+(П')-'(д»+(й-2Е)1Г'д»), й=0,1,.... (33) Для итерационных параметров Ц и а» из (26) получим (р», у») (у» и г» ~) 652 17гвва 13. Примеры численного моделирования 3. Начальное приближение. 1/2 Уо: = а Ур У: = Уо = Пуо. 4, Начальная невязка. Ь=(П*) '(О '"~-Пу -У)-1, (то, то) = (У, У)> ттО, О, О.
ттО: = тт: = ттй: = пЫ: = 5. Начало очередной итерации. Д .'= тт о тто, тт1:= 1|тт. 6. Решение системы уравнений с верхне-треугольной матрицей. О:= Оь =,7+ бь ь О, 1:=П д. 7. Решение системы уравнений с нижне-треугольной матрицей. 1:= оь = (П') '(9+271)+ 1, гч:=(гь,дь) =(г ч) аь:= тт119. 8. Пересчет итерационного приближения к решению модифиццрованной системы и невязки. где вектор е = (1, 1,..., 1).
У:=У+сьоЯ 7;= то = у — аь ь1, тт:=(ть,ть) =(г,г), нИ:= в11+1, при необходимости продолжить итерации бо го 5. 9. Обратное масштабирование (искомое приближенное решение). у:= С ч~П 'у. Ниже приведен текст программы БОТУЕ1, реализующей описан- ный алгоритм итерационного попеременно-треугольного метода — со- пряженных градиентов в варианте приближенной факторизации. Выбор матрицы С в (19) подчинен (см. п.4.7) равенству строчных сумм Ве = Ае, 13.1. Стационарнав тенлолроеодность о кусочно-однородном теле 653 С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С Б()ВК01ХПХЕ 801МЕ1 ( Х, 1., А, Х, Р, ЕР8 ) 1МРЬ1С1Т КЕА1е8 ( А — Н, Π— Е ) 1)1МЕХВ10Х А(1), У(Х), Р(Х) Параметры Х— Ь— А— У— У— ЕРЗ Описание фрагментов массива А: Элемевты главной диагонали располагаются в компонентах с 1"й по 1ч-ю.
11К=Х+1 Элементы соседней верхней диагонали располагаются порядок системы (число узлов сегочвого прямоугольника); полуширина ленты пятидиаговальной матрицы системы (число узлов сеточного прямоугольника по первому направлевию); массив, в котором задаются элементы симметричной пятидиагональной матрицы тремя диагоналями: главной, верхней соседней и верхней удаленной, что соответствует коэффициентам ревностных уравнений, действующих на компоненты сеточных функций в центральном, соседнем справа и соседнем сверху узлах пятиточечного шаблона; вслед за коэффициентами в массиве а будут расположены векторы, используемые ва итерациях; массив, содержащий на входе начальное приближение и на выходе — итоговое итерационное приближение к решению; массив, содержащий на входе правую часть системы; — параметр, задающий относительную точность итогового итерацяонного приближения к решению. ЗО1.ЧЕ1 — Подпрограмма решения системы линейных уравнений пятиточечиой ревностной схемы для самосопряжениого двумервого эллиптического уравнения второго порядка с переменными коэффициентами в прямоугольнике.
Реализован поперемевно-треугольвый метод приближенной факторнзации — сопряженных грздиевтов. 654 Т)шва 13. Примеры численного моделирования С в компоневтах с (11К+1)-й по (11К+г(-1)-ю: С Ь2К = 2Х+ Ь С С Элементы удаленной верхней диагонали располагаются С в компонентах с (12В+1)-й по (12К+К-Ь)-ю: С 10 = 3'Х + Ь 1Т =10+ Х 1ТЬ= 1Т+ Ь 1() = 1Т + Х + 1- С С С С С С Реализация алгоритма ЯОЬУЕ1. С С 1.