Вычислительная теплопередача Самарский А.А. Вабищевич П.Н. (1013606), страница 98
Текст из файла (страница 98)
Вычисление алементов диагонального множителя С матрицы приближенной факторизации; С Указатели 1С, )т, 1ТЬ, Щ используются для обрюцевия к фрагментам массива А, расположенным вслед за коаффициентами системы. РО 201= 1, Х А(1) = А(10+1)'А(1)'А(10+1) — 2.РО А(ПК+1) = А(10+1)'А(11К+1)'А(КО+1+1) А(12К+1) = А(10+1)'А(12К+1)еА(10+1+Ь) 20 СОХП1~П)Е С С 2. Масштабирование начального приближения С верхне-треугольной матрицей; С РО301=Х,1, — 1 А(1ТЬ+1) = У(1) / А(10+1) РО(01=1,Х 0 = А(1) — А(11К+1 — 1)'А(1Т+1-1)'(А(11К+1 — 1)) + А(12К+1 — 1)) Ф вЂ” А(12 К+1- Ь) *А(П+1 — Ь) '(А(11К+1 — Ц) + А(12К+1 — Ц) А(1Т+1) = 1.00 / О А(10+1) = РЯ()КТ(А()Т+1)) 10 СОХПХРЕ С С 2. Масштабирование матрицы системы диагональной матрицера С 13.1.
Стационарная теллолроводность в кусочно-однородном теле 655 У(Х) = А(1ТЬ+1) — А(11К+Х) "А(1Т1.+1+1) — А(12К+1) ьА(ХТЬ+Х+Ц 30 СО1чТ!г16Е С С 4. Вычисление начальной невнзки на месте правой части; С ККО = О.РО РО 40 Х = Ь11 А(1Т+Х) = Г(1)'А(1О+Х) — А(1)ьА(1ТЬ+Х) — У(Х) + А(11К+1-1)'А(1Т+Х-1) + А(12К+Х-Ь)'А(1Т+Х-Ц Г(Х) = А(1Т+Х) — А(1Т1,+Х) ККО = ККО + Р(Х)'Р(Х) 40 СОУП1чХФХЕ С г11Т = 0 КК = ККО КК1 = О.РО С С б. Выполнение очередной итерации; С 50 СО1ьХП1ч!ХЕ С ВК = ККьККХ КК1 = 1.РО / КК С С б. Решение системы с верхне-треугольной матрицей; С РО60Х=Х,1, -1 А(10+1) = Р(Х) + ВК*А(10+1) А(1ТЬ+Х) = А(Ь0+Х) + А(11К+Х)'А(1ТЬ+Х+1) + А(12К+Х)'А(1ТЬ+Х+Ь) 60 С01ьХП)ч!ХЕ С С 7.
Решение системы с нижве-треугольной матрицей; С Т0 = О.РО РО701=1, 1Ч О = А(1ТЬ+Х) А(1Т1.+Х) = А(10+1) + А(Х)'О + А(11К+Х вЂ” 1)*А(1Т1.+Х вЂ” 1) + А(12К+Х вЂ” Ь)'А(1ТЬ+Х вЂ” Ц А(1Т+Х) = О + А(1ТЬ+Х) Т0 = Т0 + А(1Т+1)ьА(10+1) 70 СОМП1~ПЗЕ 65б 1Лава 13. Лрциерм численного моделирования С АК= КК/ТО С С 8. Пересчет итерационного приближения к решению С и невязке; КК = О.ВО РО 80 1 = 1,11 У(1) = У(1) + АКеА(1О+3) Р(1) = Г(Х) — АК*А(1Т+1) КК = КК + Г(3)'Г(1) 80 СО1ЧТ11~ПЖ С С Выполнение очередной итерации завершено.
С Х1Т = Х1Т + 1 ЕРБМ1Т = 1331)КТ( КК/ККО ) С АНЕТТЕ ( 06, 100 ) Х1Т, ЕРБХ1Т С 1Р ( ЕРБЯ1Т .ОТ. ЕРБ ) ОО ТО 50 С С 9. Обратное масштабирование решения С верхне-треугольной матрицей; С РО901=М,1, — 1 А(1ТЕ+1) = У(1) + А(1!К+1)'А(1ТЕ+1+1) + А(12К+1)еА(1Т1+1+1.) У(1) = А(1О+1)'А(!Т1.+1) 90 СОХПХСЕ С 100 ГОКМАТ ( * Х1Т = ', 13, ' ЕРЯг11Т = ', 1312.3 ) С КЕП1Кг1 ЕХР 658 Глава 13. Примеры численногомоделирооаиил РКООКАМ ТЕЯТО! С С ТЕЯТ01 — Стационарная теплопроводноеть в куеочно- С однородном теле прямоугольного сечения.
С 1МР?.1СГГ КЕА1з8 ( А-Н, Π— Е ) КЕА)'8 КАРРА С РАКАМЕТЕК ( ГР1М = 50000 ) С ))1МЕ)ч)Я10)ч) А(ГР1М) С С С С С С С С С Способ размещения указанных фрагментов в массиве А, С описываемый ниже. позволяет взбежать в алгоритме С ЯОЬЧЕ1 специального рассмотрения компонент векторов С и строк матриц, соответствующих граничным узлам С сеточного прямоугольника.
С С Если г! — число компонент вектора неизвестных рээноегной С задачи, то размещение фрагментов опнсываегея списком С эквивалентностей: С С С С С С С СОММО)ч) / ТО! / Х!1., Х!Р, Х!К, ь Х21., Х2Р, Х2К, Длина массива А должна быть достаточной для размещения коэффициентов симметричной матрицы разностной задачи, заданной главной АО, соседней верхней А1 и удаленной верхней А2 диагоналями, векторов решения У н правой части Р, а также векторов, участвующих в итерационном алгоритме решения раэностного уравнения (см. алгоритм ЗОЕЧЕ1). Е()ШЧА(ЕНСЕ ( А(1), АО ), Ф ( А()э+2), А1 ), * ( А(2зн+)Ч1+1)„А2 ) * ( А(зе)я+1), У ), ( А(8"')Ч+1), Р ) 13.1.0.
Программа Для решения задачи теплопроводностн в кусочно-однородном теле использовалась программа ТЕЯТ01. Она основывается на использовании подпрограммы РРЯ01 для формирования сеточной задачи н подпрограммы ЯРЧЕ! для итерацнонного решения сеточной задачи. Поэтому ограничимся лишь приведением самого текста расчетной программы. Ф КАРРА, В1, з Х1, Х2, Н1, Н2 с С Ввод данных задачи: С С С С С С С С С С С С С С Х11, Х2Ь вЂ” координаты левого нижнего угла расчетной прямоугольной области; Х1В, Х2К вЂ” координаты правого верхнего угла; ХЮ, Х2 — координаты правого верхнего угла подобласти; КАРРА — безразмерное отношение коэффициентов теплопроводности; В1 — число Био; Х1, Х2 — число узлов сетки по соответствуюшим направлениям; ЕРБ — требуемая относительная точность итерационного приближения к решению.
Х!Ь = ОЛ?0 Х11? = 0.51?0 Х1К = Ь1?0 Х2Ь = 0.1?0 Х21? = 0.500 Х2К = Ь00 КАРРА = 100.1?0 В1 = 10.00 Х1= 51 Х2= 5! ЕР8 = 1.1?-б ОРЕХ ( Об, Н1.Е = 'К01.РТС' ) 1?О 1 1 = 1, 1Ь?1М А(1) = 0.00 1 СОХТ1ХЬ?Е С Х = Х1'Х2 С С С С С В подпрограмме И%01'определяются центральный АО, правый А1 и верхний А2 коэффициенты раэностной схемы на пятиточечном шаблоне и правая часть Р. САЬЬ НЭВО! ( А(1), А(Х+2), А(2'Х+Х1+1), А(8"Х+1) ) 13.1. Стационарная темояроводноетз в кусочно-однородном теле б59 660 ОРЕХ ( 01, Р11 Е = 'К01.РАТ' ) %К1ТЕ ( 01,' ) (А(7*Х+1), 1=1,Х) СЬОБЕ ( 01 ) С1.ОБЕ ( 06 ) БТОР ЕХП Подпрограммы-функции К„К1, К2 используются для вычисления козффициентов разностной задачи в соответствии с (11), (12).
13.1АЬ Примеры расчетов По приведенной программе выполнены методические расчеты, которые направлены на исследование как самих вычислительных алгоритмов, так и зависимости самого решения от физических параметров. На рис. 13.2 представлена зависимость погрешности приближенного решения задачи (1), (3), (6)-(9) в единичном квадрате от номера итерации. Данные приведены для варианта с к = В) и В1 = 10 на трех различных сетках. Начальное приближение в расчетах уе(я) = О, я 6 й.
Величина относительной погрешности на итерации определяется отношением На рис. 13.2 кривая 1 соответствует расчетам на сетке 26 х 26, кривая 2— на в два раза более подробной сетке (51 х 51), кривая 3 — 101 х 1О1. Данные демонстрируют теоретическую зависимость числа итераций попеременно-треугольного метода от числа узлов сетки (см. п.4.7)— при увеличении числа узлов по каждой переменной в четыре раза, число итераций увеличивается примерно в два паза. С С С С С С С С 1Лава 13. Примеры численного моделиросалия В подпрограмме 801ЛЕ1 находится решение разноствой задачи итерационным попеременно-треугольным методом приближенной факторизации — сопряженных градиентов.
Начальное приближение передается в массиве А, начиная С (7*гч+1)-й компоненты; на этом же месте на выходе содержится итерационное приближение С относительной точностью Еря. СА1Ь БО?УЕ1 ( Х, Х1, А(1), А(7'Х+1), А(8'Х+1), ЕРБ ) е = е(й) = —, гь = (11 ) 6 7 (1 — Ауь). !!~ й!! ~ — 1 -1!2 !!гс!!' 13.1. Стационарная теплопроводность в кусочно-однородном теле 661 1оде 5 1.0 -10 0.0 05 Рие. 13.3 1.0 Рис. 13.2 и 00 и(я) 1.0 яа 0.5 .5 1.0 0.0 0.0 05 0.0 0.0 1.0 я! я, Рис. 13.3 Рис. 13.4 На других рисунках представлены изотермы я(х) = сопа1 с шагом дв = 0,05 при различных значениях определяюших параметров задачи. Расчеты выполнены на квадратной сетке 51х 51.
На рис. ! 3.3 представлено поле температур в задаче с к = 10 и В! = 10. Влияние изменения безразмерного числа к (теплопроводности в подобласти й') иллюстрируется рис. 13.4, когда к = 0,1 и Вг = 10. При к — О реализуется граничное условие первого рода на внутренней границе включения й'. Этот случай соответствует использованию метода фиктивных областей (см. п. 4.8) при решении задачи Дирихле в нерегулярной области й". Аналогично при к — 0 реализуются однородные условия Неймана на внутренней границе включения. Изменение теплового режима на границе прослеживается на рис.
13.5, где к = !О и В! = 0,5. 662 Евана 13. Примеры численного моделирования 13.1.г. Задачи Задача 1. На основе приведенной программы проведите вычислительные эксперименты по исследованию скорости сходимости итерационного попеременно-треугольного метода — сопряженных градиентов от: (а) изменения безразмерного коэффициента теплопроводности к; (б) изменения числа Био; (в) изменения числа узлов по отдельным направлениям. Задача 2. Проведите численное исследование теплового состояния твердого тела с включением в широких пределах изменения: (а) геометрии тела (1з) и его включения; (б) теплофизических характеристик включения (к); (в) граничных режимов (В1).
13.2. Затвердевание расплава в полости прямоугольной формы 13.2.1. Постановка задачи Рассмотрим задачу затвердевания чистого материала, который заполняет полость прямоугольного сечения, Будем считать, что полость поддерживается при постоянной температуре окружающей среды, верхняя поверхность расплава контактирует с атмосферой (рис. 13.6). Тепловое состояние затвердевающего расплава с учетом энтальпии фазового перехода описывается уравнением теплопроводности в виде (см. п.
2.3 и главу 7) Р . 1З.Е дн д l ди 1 (с(н) + Лб(а — и')) — — ~~> — 1 к(а) — ) = О, (1) 31, Ох. ~ дх.) = где н' — температура, а Л-знтальпия фазового перехода. Для простоты ограничимся случаем, когда коэффициенты теплоемкостн и теплопроводности постоянны в твердой и жидкой фазах и не меняются при фазовом превращении. В силу симметрии достаточно рассмотреть задачу в половине сечения й = (х ( х = (хн хз), О < х„< 1„, а = 1, 2 ). 13.2. Затвердевание расплава в полости прямоугольной формы 663 Поэтому на левой границе имеем ди — =О, х, =О, О<х,<1,. (2) дх Допуская неидеальный контакт с полостью, на правой и нижней границах зададим граничное условие третьего рода".
ди й — + о,(и — и,) = О, х~ — — 1п О < хз < 1н (3) дх~ ди -й — +о,(и — и,)=0, хз=О, 0<х, <1ь (4) дхз где и, = соим < и' — температура окружающей среды. В (3), (4) о, ~ со соответствует случаю неидеального контакта затвердевающего расплава с границей полости. На верхней границе имеет место аналогичное граничное условие с коэффициентом конвекгивного теплообмена оз « о1. ди й — +о (и — и,) =О, ~=1, О<х~ <1,.