Вычислительная теплопередача Самарский А.А. Вабищевич П.Н. (1013606), страница 109
Текст из файла (страница 109)
1!римеры численного моделироеонил Нижняя и верхняя границы: краевое условие третьего рода РО 401= 2, Х1-1 Х! = Х1Ь + (1-1)'Н! А1(1 — 1,1) = 0.500'Н21'К1(Х1 — 0.500'Н1,Х2Ь) А1(1,1) = 0.500'Н21'К1(Х1+0.5ОО'Н1,Х2Ь) АО(1,1) = А1(1 — 1,!) + А1(1,1) + А2(1,1) + Н1'В1 В(1,1) = О.ОО А1(1 — 1,Х2) = 0.5ОО'Н21'К1(Х! — 0.5ОО'Н1,Х2К) А1(1,Х2) = 0.5РО*Н21'К1(Х1+0.5РО*Н1,Х2К) А2(1,Х2) = О.ОО АО(1,Х2) = А1(! — 1„Х2) + А!(1,Х2) + А2(1,Х2-1) + Н1'В! г(1,Х2) = О.ОО 40 СО1ЧТ1Х!1Е С С Правая граница: краевое условие третьего рода С РО 50 У = 2, Х2-1 Х2 = Х2Ь + (У вЂ” 1)'Н2 А1(Х1,2) = О.ОО А2(Х1,2 — 1) = 0.5ОО'Н12'К2(Х1ЕгК2 — 0.5РОеН2) А2(Х1,3) = 0.5РО'Н12'К2(Х1К,Х2+0.500'Н2) АО(Х1,Х) = А1(Х1 — 1,2) + А2(Х1,1 — 1) + А2(Х1,1) Э + Н2'В! Г(Х1,2) = О.ОО 50 СОХТРШЕ С С Левый нижний угол С АО(1,1) = А1(1,1) + А2(1,1) + 0.5130'Н!'В! Р(1,1) = О.ОО Левый верхний угол А2(1,Х2) = О.ОО АО(1,Х2) = А1(1,Х2) + А2(1,Х2 — 1) + 0.5ОО'Н1'В1 г(1,Х2) = О.РО 13.6, Расчел2 л2ермоетал2а С Правый нижний угол С А1(Х1,1) = О.ОО АО(Х1,1) = А1(Х1 — 1,1) + А2(Х1,1) + 0.5ОО»(Н1+Н2)»В1 Р(Х1,1) = О.ОО С С Правый верхний угол С А1(Х1,Х2) = 0.00 А2(Х1,Х2) = О.ОО АО(Х!,Х2) = А1(Х1 — 1,Х2) + А2(Х1,Х2-1) + 0.5130»(Н1+Н2)»В1 Р(Х1,Х2) = О.ОО С КЕТ1ЖХ ЕХВ Подстановка (11) в (10) даст ° М ° 2 ! 2 аа(о),~~ Х~~ сара(х) 1 121122+ 2 Х~~ оа аа~ »=1 а=! (13) ° М 1 гсдрд(х) — 1)р,(х)Ь1Ь2+ — гс» = О, а = 1,2,...,М.
Для решения системы линейных уравнений (14) используется обычный метод !пусса. Я1ВКО11Т!ХЕ ОА1133 ( М, Р, % ) С С Решение системы уравнений с плотной матрицей Р С методом Гаусса беа выбора ведущего елемента. С решение помещается в массив правых частей %. С 1МР13С1Т КЕА! »3 ( А-Н, О-Е ) Р!МЕХЯ!ОХ Р(М,М)„%(М) С ВОЗОК=1,М-1 ЭО 20! = К+1, М К = Р(1,К) / Р(К,К) Условие (6) для функционала (13) дает следующую линейную систему уравнений для определения оптимального управления нп 748 Вава 13. Примеры числеяиоеомоделирования РО 1О Х = К+1, М Р(1,1) = Р(1,1) — К'Р(К,Х) 10 СОРГПР(ЮЕ %(1) = %(1) — К'%(К) 20 СО)ШТРИХЕ 30 СОРГО(ХЕ РО501=М,1, -1 Я =ОВО 1Р ( 1 ХТ.
М ) ТНЕг( РО 40 Х = 1+1, М Я = 8 + Р(1,1) %(Х) 40 СОМПг((ХЕ Ег(Р 1Р %(1) = (%(1)-8) Х Р(1,1) 50 СОРГП1С((ХЕ С КЕТУ КХ ЕМЭ Таким образом, реализация алгоритма для рассматриваемой задачи основана на решении М разностных задач (12) и определении искомого вектора управления нз линейной системы уравнений (14). Алгоритм реализован в виде следующей основной программы.
РКООКАМ ТЕБТ06 1МР(.1С1Т КЕАЬэ8 ( А — Н, О-Х ) РАКАМЕТЕК ( 1Р1М = 40000, М = 7 ) Р1МЕ)г(31Ог( А(1РХМ), Р(М,М), %(М), Х23РАТ(М) Длина массива А должна быть достаточной для раэмещения коэффициентов симметричной матрицы раэиостной задачи, эаданной главной АО, соседней верхней А1 и удаленной верхней А2 диагоналями, векторов решения У(1) для 1 1,2, „М и правой части Г, а также векторов, участвующих в итерационном алгоритме решения ревностного уравненяя (см. алгоритм 801ЛЕ1).
С Если )Ч вЂ” число компонент вектора неиэвестных раэностной С ввдачн, то раэмещение фрагментов описывается списком С эквивалентностей: 13.6. Расчет термостата 749 С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С С ЩШУАЬЕМСЕ ( А(1), АО ), я ( А(М+2), А1 ), * ( А(2еМ+М1+1), А2 ), * ( А(7еМ+1), Х ), * ( А((7+М)еМ+1)„Р ) СОММОМ / Т06 / Х1Ь, Х12К, Х13К, Х1К, Х21., Х221., Х22К, Х23К, Х2К, С2, СЗ, В1, е М1, М2, Н1, Н2 СОММОМ / ЗОЫКСЕ / Х13, Х2$ Ввод данных аадачи: Х11., Х21. — коордяваты левого нижнего угла корпуса термостата: Х1К, Х2К вЂ” коордиваты правого верхнего угла корпуса термостата; Х11., Х221 — координаты левого нижнего угла объекта термостатяровавия; Х12К, Х22 — координаты правого верхнего угла объекта термостатировавия; Х1ЗВ, ХЗЗ — коордяиаты правого верхнего угла камеры термостата; С2, СЗ, В1 — звачеиия беаразмерных теплофизических параметров в ураввеиии и граничном условии; М вЂ” количество источников тепла в корпусе термостата; Х18 — гориеонтальиая координата ливии, иа которой расположены источники тепла; ХЗВВАТ вЂ” массив вертикальных координат источников тепла; М1, М2 — число уелов сетки ло соответствующим иаправлевиям; ЕР — требуемая относительная точность итерациоивого приближения к решению; ЕРЕМ вЂ” параметр, определяющий множитель штрафа в функционале качества.
Х1Ь = 0.1)0 Х12К = 0.300 Х13К = О.ООО Х1К = 1.ОО Х2Ь = 0.1)0 Х221, = 0.21)0 750 ядава 13. Примеры числеияогомоделирования С С С 1 С С С С С С С С С С С С С С С Х22К = 0.71)0 Х23К = 0.81)0 ХгК = 1.1)О Х18 = ОЛ)0 Х2ЯЭАТ(1) = 0.2ВО Х2Я)АТ(2) = 0.3130 Х281)АТ(З) = 0.4130 Х2$13АТ(4) = 0.51)О Х2Я)АТ(5) = 0.61)0 Х2ЯЭАТ(6) = 0.71)0 Х2Я)АТ(7) = 0.81)0 С2 = 10.130 СЗ = 0.1130 В1 = 10.130 ЕРЕХ = 1.)3+3 Х1 = 51 Х2=51 ЕРБ = 1.1) — б РЕХАЬТ = 1.1)0(ЕРЕХ'ЕРЕХ) ОРЕХ ( 06, ЛЬЕ - 'ЕОО.РТС' ) 1)0 1 1 = 1, Ю1М А(1) = 0.1)0 СОХТ)ХЦЕ Х = Х1'Х2 В следующем цякле определяются фуякции влияния отдельных источников 1решения краевых задач для уравнения с правой частью з виде соответствующей дельта-функции). 1)О10Ь=1,М Х28 = Х2ЮАТ(Ь) В подпрограмме И%06 определяются центральный Ао, правый А1 и верхний А2 коэффициенты раэностяой схемы на патиточечном шаблоне и правая часть г.
САЬЬ НЭБОб ( А(1), А(Х+2), А(2'Х+Х1+1), А((7+М)"Х+1) ) В подпрограмме ВОЬУЕ1 находится решение ревностной задачи итерационным поцеремеино-треугольным методом 13.6. Расчет термсстата 751 приближенной факторнэацин — сопряженных градиентов. Начальное приближение передается в массиве А, начиная с ((б+Ц"Х+1)-й компоненты; на атом же месте на выходе содержится итерационкое приближение относительной точности ЕРВ. САЬЬ БОЬчЕ1 ( Х, Х1, А(1), А((6+Ь)'Х+1), А((7+М)*Х+1), ЕРБ ) С 10 СОМПХЦЕ С С С С С С В следующем цикле определяются коэффициенты и правая часть системы уравнений с симметрячиой положительно определенной матрицей относительно %(Ь), 1 1, ...,М, выражающей необходимо условие минимума функцяонвла качества: ПО 30 1.
= 1, М Р(1.,1.) = К1Р(ЗУ ( А((6+1,)'Х+1), А((6+1.)'Х+1) ) + РЕХА1Т %(1.) = К1РЫ1 ( А((6+Ь)'Х+1) ) 1Р ( Ь .ОТ. 1 ) ТНЕХ ОО 20 ЬЬ = 1, 1.-1 Р(1„ЬЬ) = К)Р()ч' ( А((6+Ь)*Х+1), А((6+ЬЬ)*Х+1) ) Р(ЬЬ,)) = Р(Ь,ЬЬ) 20 СОХТ(Х()Е ЕХ)Э 1Р 30 СОХПХ(ЗЕ С С С С В подпрограмме САСЯЗ методом Гаусса решается система уравнений отиосятельно вектора оптимального управления % 1)О 50 1 = 1, Х А(1) = 0.1)0 1)О 40 Ь = 1, М А(1) = А(1) + %(Ь)"А((6+Ь)*Х+1) 40 СОЫПХЦЕ 50 СОХТ1Х()Е С %БОМ = 0.1)0 СА1 Ь ОА(3ББ ( М, Р, % ) С С Находится распределение температуры в случае испольэования С вектора оптимального управления мощностью источников: С 753 13.6. Расчет термостата СОММОЙ / Т06 / Х11., Х12К, Х13К, Х1К, Ф Х21., Х22Ь, Х22К, Х23К, Х2К, Ф С2, СЗ, В1, Ф 'с11, Х2, Н1, Н2 К1 = 0.5РО'( К(Х1,Х2-0.5РО'Н2) + К(Х1,Х2+0.5РО*Н2) ) КЕТРКй1 ЕХР РОУВЬЕ РКЕС1810й( ИЛЧСПОй1 К2 ( Х1, Х2 ) 1МРЫС1Т КЕА1 ей ( А-Н, 0 — Е ) КЕАЬе8 К С СОММ01ч' / Т06 / ХН., Х12К, Х13К, Х1К, е Х21 Х22Ь, Х22К, Х23К, Х2К, С2, СЗ, В1, е Х1, 1ч2, Н1, Н2 С К2 = 0.500*( К(Х1 — 0.5РО'Н1,Х2) + К(Х1+0.5РО'Н1,Х2) ) С КЕТЛЕ~ ЕНР РОУВЬЕ РЕЕС!ЯОР1 РР1чСПОХ К1Р1ЗЧ ( Ц Ч ) С С Вычисление скалярного проивведения двух функций на С множестве сеточных функций, ваданкых в области С объекта термостатировения С 1МР1.1С1Т КЕАЬ'8 ( А-Н, О-Е ) СОММОЙ / Т06 / Х1Ь, Х12К, Х13К, Х1К, Х21., Х22Ь, Х22К, Х23К, Х2К, Ф С2, СЗ, В1, Э Х1, Х2, Н1, Н2 Р1МЕ1ЧВ10Н ЩХ1,Х2), Ч(Х1,0(2) 1В = 2 1Е = 1ХТ ( (Х12К вЂ” Х1Ь)Н1 ) ЗВ = 2 + 1о1Т ( (Х22Ь вЂ” Х2Ь)Н2 ) 754 УЕ = 1Р1Т ( (Х22К-Х2Ь)Н2 ) Б = О.РО РО 201=!В,ЗЕ РО 1О ! = 1В, 1Е 5 = Б + 1З(11)'Ч(ЬЗ) 10 СОР1П1ЧУЕ 20 СОР1ПНЮЕ С К1РЫЧ = 5'Н!"Н2 С КЕТ11Кс! ЕНР РОУВЬЕ РКЕС1ЯОР! ИПчСПОР! К!Р!З1 ( Ы ) С С С С С Вычисление скалярного произведения заданной функции и единицы иа множестве сегочных функций, заданных в области объекта термостатирования.
1МРЫСП КЕА1е8 ( А — Н, О-л. ) С СОММОН / Т06 / Х11., Х12К, Х13К, Х1К, Ф Х21., Х221., Х22К, Х2ЗК, Х2К, Ф С2, СЗ, В1, Ф Н1, Х2, Н1, Н2 Р1МЕНЯОН 1З(М1,Р12) 1В = 2 !Е = 1ХТ ( (Х12К-Х1Ь)Н1 ) ЗВ = 2 + 1Р1Т ( (Х22Ь вЂ” Х2Ь)Н2 ) 1Е = !е1Т ( (Х22К-Х2Ь)Н2 ) К!РЩ 5еН1еН2 С Б = ОРО РО 20 1 = ЗВ, 1Е РО 10 1 = 1В, 1Е $ = Б + Щ1,1) 10 СОХПе1ЦЕ 20 СОРУП1~П1Е С 1Ъава 13.
Примеры численного моделирования 13.6. Расчет термостата 755 иЕТУЕХ Е1ч13 Здесь подпрограммы функции К, К1 и К2 используются для вычисления коэффициентов разностного уравнения по разрывному коэффициенту теплопроводности. Назначение других подпрограмм не требует каких-либо дополнительных пояснений. 13.6.4. Примеры расчетов Таблица 18.1 е, е2 ез ег В!=10 еж!02 4,239 4,230 3,849 3,428 2,864 4,10! 2,257 В! ~10 е = 10~ 5,64! 5,163 3,810 3,071 1,681 4,512 В!= 10 е = 104 -0,1872 -0,6631 13,25 7,547 2,527 0,3146 1,703 В!= !О! е = 1О' 8,022 6,694 5,042 3,356 8,663 9,099 8,839 В1=! е= 10~ 0,7207 0,6233 0,4566 0,3965 0,3496 0,53!9 0,3133 При расчете режима термостата необходимо контролировать не только достигаемую степень близости температуры в объекте термостатирования к номинальной, но и за тем, какими величинами источников тепла поставленные цели достигаются.
В частности, практическая реализация оптимального режима термостатирования должна достигаться только положительными источниками тепла. С этой целью проводится варьирование параметра штрафа е. Рассматривается вариант расчета, когда управление осуществляется семью источниками, расположенными на линии х! = 0,7, через равные интервалы от хз —— 0,2 до хз = 0,8. Толщина корпуса термостата снизу и сверху 0,2, а сбоку 0,4. Горизонтальный размер термостата равен 2, вертикальный равен 1, размеры термостатируемого объекта 0,6 и 0,5 соответственно. Принято„что в безразмерных переменных коэффициент теплопроводности в термостатируемом теле равен 10 (лз), в камере термостата 0,1 (лз).
Расчеты выполнены на равномерной сетке 51 х 51. Влияние параметра штрафа е иллюстрируется данными в следующей таблице. 756 1)в(вя 13, Примеры числеииого моделироеоиия в(к) е(к) 1.0 1.0 к кз 0.5 0,5 1.0 ' 0.0 0.0 О. .0 0.0 к, Рве. 13.23 Рие. 13.26 е(к) 1.0 0.5 0.5 0.0 Рие. 13,27 При увеличении е (уменьшении штрафа) (см. данные по задаче с В! = 10) наступает момент (е = 10"), хогда среди источников появляются отрицательные. При этом, конечно, температура в образце все более точно приближается к номинальной.
В табл. 13.1 приведены также рассчитанные величины источников для чисел Био ! и 100. Этими расчетами иллюстрируется, в частности, влияние температуры окружающей среды. На рис. 13.25 представлено рассчитанное поле температур для варианта с В! = 10 н параметра штрафа е = 1О'. Непрерывные кривые соответствуют изотермам и(х) = )е, )е = 1, 2,.... Пунктирные линии представляют собой изотермы между и = О и н = 2 с шагом 0,2. Видно, что с хорошей точностью температура в термостатируемом образце близка к необходимой температуре, равной 1. Аналогичные выводы можно 757 13.7. Восстановление внешних тепловых нагрузок сделать и относительно вариантов с В1 = 10 (рис.
13.26) и В1 = 1 (рис. 13.27). Здесь также параметр штрафа выбран равным 10з. Величины токов для этих вариантов приведены в табл. 13.1. 13.0.5. Задачи Задача 1. Проведите исследование оптимальных температурныхполей в термостате при изменении следующих вычислительных параметрот (а) числа узлов по пространственным переменным; (й) параметра штрафа е.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
