Вычислительная теплопередача Самарский А.А. Вабищевич П.Н. (1013606), страница 100
Текст из файла (страница 100)
САЬЬ ГТ)802 ( А(1), А(Х+2), А(2'Х+Х!+1), А(7'Х+1), А(8эХ+!) ) В подпрограмме ЯОЬУЕ1 находится решение разностной задачи итерационным попеременно-треугольным методом приближенной факторизации — сопряженных градиентов. начальное приближение передается в массиве А, начиная с 17"Х+1)-й компоненты; на атом же месте на выходе содержится итерационное приближение с относительной точностью ВРИ.
САЬЬ БОЬУЕ1 ( Х, Х1, А(1), А(7'Х+!), А(8'Х+1), ЕРБ ) ХТА1) = ХТАЬ) + 1 МТАЬ) = МТА1) + 1 Запись решения в момент времени Т 1Е ( МТАЬ) .ЕО. М ) ТНЕХ Т = ХТА1)*ТАЫ 13.2, Затвердевание расплава в яовооии прямоугольной формы 671 %К!ТЕ ( 01,' ) Т %К!ТЕ ( 01,' ) (А(7*Х+1), 1=1,Х) МТАБ = 0 ЕХП 1Г 00401=!Х 1Р ( А(7ьХ+!),ОЕ. 0.00 ) ОО ТО 30 40 СОМПХ1!Е С С Запись решения на момент полного затвердевания С Т = ХТА1! ТАП %КПЕ ( 01,* ) Т ттК1ТЕ ( 01,' ) (А(7*Х+1), 1=1,Х) С1.ОБЕ ( 01 ) С СЕОБЕ ( 06 ) С ВТОР ЕХП 001!В1Е РКЕС1210Х И!ХСПОХ В ( 17 ) КЕАЬ'8 Ц ПЕ1., ННТ, 8ННТЕ СОММОХ Т02В ОЕ1., ННТ, БННТЕ В =ННТ 1Р ( 0АВ8(!!) Л.Е.
0Е1, ) В = 8ННТЕ Подпрограмма-функция В используется для вычисления эффективного коэффициента теплоемкости в соответствии с (! 4). 13.2.4. Примеры расчетов Приведем некоторые данные характерных расчетов, лля выполнения которых использовалась описанная программа. В рассматриваемой задаче Стефана интерес, прежде всего, представляет динамика границы фазового перехода и, в частности, время полного затвердевания.
Бгава 13. Примеры численного моделирования 672 и(х,г) и(х,г) !.О хз кз 0.5 0.5 0.0 0.5 1,0 0.0 0.0 х Рис. 13Л ОД 1.О х, Рие. 13.3 иЬл) и(хд) !.0 хз хз 0.5 0.5 1.0 0.0 0.0 0.5 !.0 к, х~ Рие. 13.10 О.О 0.0 0.5 Рис. 13.0 На рис. 13.7 представлены изотермы в(х,() = сопя! через би = 0,1 на момент времени ! = 0,05 лля задачи (2), (6)-(10) при В1~ —— 100, В11 = 0,5, 8!е = 0,25. Начальная температура расплава й = 0,25. Расчеты выполнены на сетке 51 х 5! для шага по времени т = 0,001, полуширине сглаживания б-функции Ь = 0,04. Использовалась чисто неявная линеаризованная схема, т.
е. в (13) о = 1. Для этого же варианта на рис. 13.8 представлена граница фазового перехода на различные моменты времени, через б( = 0,025. Время полного затвердевания расплава а 0,18, Влияние энтальпии фазового перехода (числа Стефана) иллюстрируется рис. !3.9. Здесь, в отличие от рис. 13.8, 8(е = 0,5. Изменение энтальпии фазового перехода существенно влияет на процесс затвердевания, при увеличении 8!е возрастает время затвердевания. 673 13.3.
Переизлучение в твердом теле с невыпуклым сечением 13.2.6. Задачи Задача 1. Проведшие исследование времени затвердевания при изме- нении следующих вычислительных параметров: (а) числа узлов по пространственным переменным„' (б) шага сетки по времени; (в) веса о' разностной схемы; (г) полуширины сглаживания б-функции с3; (д) относительной точности итерационного процесса е. Задача 2.
Исследуйте влияние на процесс затвердевания: (а) высоты полости 1з! (б) начальной температуры расплава й; (в) знтальпии фазового перехода (Бге); (г) условий охлазкдения (В!с, В1з). 13.3. Переизлучение в твердом теле с невыпуклым сечением 13.3.1. Постановка задачи Рассмотрим модельную (см. и. 8.3, 8.4) задачу теплообмена теплопроводиостью и излучением для твердого тела с невыпуклым Ь-образным сечением й. Расчетная область й состоит (рис. 13.11) из двух прямоугольников: й = й' ! ! й", й' = (х ~ х = (хп хз), йн = (х ~ х = (хн хз), Невыпуклую часть границы обозначим с ч 7=7 г!7 0<х,<1', аш1,2), 0<х <1", а=1,2).
где у =(х!хбдй, х!=1,), у =(х)хбдй, хг=1з) и Г= бй'!у. Представляет интерес исследование влияния условий охлаждения. На рис. 13.10 представлены границы затвердевания для задачи, когда в отличие от рис.! 3.8 В!з = 0,25. Уменьшение интенсивности теплообмена с верхней границы замедляет образование затвердевшей корки на этом участке границы. Приведенная программа позволяет провести численное исследование процесса криствтлизации или плавления в различных технологических условиях. 674 Глава 13. Примеры численного моделиргманил Будем исследовать установившуюся температуру твердого тела с сечением й в условиях конвективного теплообмена с окружающей средой постоянной температуры и„когда нагрев тела обусловлен объемными ис- 1' точниками.
Основной вопрос связан с учетом переизлучения на невыпуклой части границы у. Стационарное поле в изотропном 1г твердом теле описывается уравнением хай, к Рис. 13.11 где 1(х) — плотность объемных источников тепла. Потери тепла на границе обусловлены конаекгивным теплообменом с окружающей средой и излучением. Пренебрегая излучением окружающей среды, получим ди 1 — + а1(и — и,) + д(х) = О, х б дй, (2) где о, — коэффициент конвективного теплообмена, а д(х) определяет радиационный поток. На выпуклой части границы имеем д(х) = до(х, и), (3) хЕГ, где чо(х, и) = ази (х), (4) а о з — излучательная способность тела.
На невыпуклой части границы 7 имеем д(х) = о (х) — д+(х), где о (х) — поток испускаемого излучения, а д~(х) — поток падающего излучения (д(х) — результирующий поток). При заданных температурах д (х) определяется из интегрального уравнения д (х) — у С(х,~)д (С)НС =до(х,и), х Е у, (5) где т — коэффициент отражения. Ядро этого интегрального уравнения определяется (см.
п. 2.5, 8.3) выражением 1 б(х, Д = соз (и(х), г) соз (и(~), г), тг(х, е) (6) б75 13.3. Переизлучеиие в твердом теле с иевыиуилым еечеиием где г(х, 4) — расстояние между точками х и 4, а соз (я(х), г) — косинус угла между нормалью к 7 в точке х и отрезком, соединяющим * и С.
Для д+(х) имеем выражение д+(х) = б(х,1)д (с)дс, х Е у. 7 Сформулируем задачу (1)-(6) в безразмерных переменных. Будем считать, что твердое тело однородно (1г(х) = й = сопя!), а тепловые источники распределены равномерно (у(х) = 7 = сопя!). В качестве характерного линейного размера возьмем 1",, безразмерную температуру опрелелим как отношение и/и,. Тогда в безразмерных переменных уравнение теплопроводности примет вид дзи — — 2=02, хай, , дх2 (7) где Оз — число Остроградского: Ох=в У (1н) 2 йи, Интенсивность контактного теплообмена характеризуется числом Био: 1н В! = —. й Мерой отношения переноса тепла излучением к потоку тепла за счет теплопроводности служит число Кирпичева и21! ин и 3 К! = —. й Граничное условие (2) в безразмерных переменных примет вид ди — + В!(и — 1) + К1 д(х) = О, х Е дй, д где радиационный поток определяется согласно (3) на Г и из интеграль- ного уравнения (5) на у, причем теперь йв(х, и) = и (х).
(9) Поставленная задача (3), (5), (7)-(9) является нелинейной и характе- ризуется следующими основными безразмерными параметрами: Оа, В1, Кз, л. Ее основная особенность заключается в нелокальности граничнык условий на у. 676 Птава 13. Примеры чиелеииоео моделирооаиия 13.3.2. Разностная задача Для численного решения задачи введем равномерную сетку с шагами йт, 'аг по пространственным переменным хт и хг соответственно.
Будем считать, что сетка согласована с границей расчетной области, т. е. узлы лежат на границе дй. Простейший пример такой сетки приведен на рис. 13.12. Приближенное решение для температуры обозначим у(х), (" г а радиационный поток — е(х). Пусть ы — множеспю внутренних узлов сетки. Уравнение (7) аппроксимируется на обычном пятиточечном шаблоне с помощью разиостного уравнения Ряс. 18.12 увил — уе,*, = Оз, х Е ы. (10) Сеточная задача в граничных узлах записывается с учетом граничных условий третьего рода (8). пусть ть', 7'„' — множество узлов, лежащих на У' и Уо соответственно (Уь — — ть О тьи) без Узлов (1'„1',), (1",,1,"). На нижней части границы области й, аппроксимнруя уравнение (7) и условие (8), получим 2 2, 2 — (-аг(х)у„+ В1 у) — (от у;,),, + — Кл е(х) = — В1+ Оз, йг йг йг (11) хг — — О, 0<хт <1т.
Аналогичные аппроксимации используются на других участках границы и, в частности, в угловых точках. Разностное уравнение лля температуры запишем в виде Лу+ Ро(х)е = Рт(х), х Е то, (12) где й — линейный сеточный оператор, определяемый в соответствии с (10), (11) и аналогичными соотношениями в других граничных узлах а Е ды = й ~ ы = Гь О уь. Сеточная функция Ро(х) отлична от нуля только в граничных узлах. Сформулируем теперь разностную задачу для радиационных потоков.
Из (3), (9) имеем е(х) = Рг(у), х Е Гю (13) где Рг(у) = у~(х). (14) Аппроксимацию интегрального уравнения (5) проведем на основе интегро-интерполяционного метода (более подробно см. п. 8.3). Упорядочим гРаничные Узлы хт Е тю г = 1,2,...,тл и свЯжем с каждым 13,3. Лереизлучение в твердом теле с невмлуквмм сечением 677 таким Узлом окРестность гРаницы уы.
К х; отнесем половины отРезков, соединяющих соседние граничные узлы. При интегрировании интегрального уравнения (5) по отмеченным окрестностям получим з (х!) — х~ ~!рбз (ху) = Рг(у(х!)), ! = 1,2,...,пг. (15) г=! Для угловых коэффициентов в (15) используются расчетные формулы следующего вида: шеез 7м аэ Их. (1б) — / 6(х,б)Щдх, 1 щезз7л ./ з(х;) =з (х;) — ~ р; з (х.), ! = 1, 2,..., пг. Задачу (13)-(1б) по вычислению радиационного потока запишем в сле- дующем операторном виде (17) Здесь оператор Я отличен от нуля только на множестве граничных узлов 7л. Для приближенного решения нелинейной разностной задачи будем использовать простейший итерационный метод последовательного уточнения потерь на излучение.