Вычислительная теплопередача Самарский А.А. Вабищевич П.Н. (1013606), страница 99
Текст из файла (страница 99)
(5) дхз о~11 В1~ = —, й В безразмерных переменных уравнение (1) принимает вид (1+Ягео(и)) — — ~~~, — = О, х Е Й, 1) О. (7) ди ди 01 дхз 1)>аничное и начальное условия (2), (6) остаются, а (3) — (5) дают ди — +Вй(и+1)=О, х =1, 0<хз<1н (8) дх Будем считать, что расплав заливается при некоторой температуре ио (ио = сопзг ) и'), т. е. начальное условие имеет вид и(х, О) = и, х б Й.
(6) В принятом приближении Стефана граница фазового перехода Я = Я(1) в каждый момент времени определяется следующим образом: Я = (х1х 0 Й, и(х,1) = и'). Уравнение (1) с условиями (2)-(6) полностью описывает температурное поле при затвердевании расплава. Сформулируем поставленную задачу в безразмерных переменных. Клк обычно, для безразмерных величин используются те же обозначения что и для размерных. Обезразмеривание проведем на 1ы безразмерную температуру опрелелим отношением (и-и*)/(и' — и,). Задача характеризуется безразмерными числами Стефана и Био: Л оз11 Вге = В1з = —. с(и' — и,) й 664 Рава 13. Примеры численногомоделироеания, ди — — + В!1(в + 1) = О, дх2 ди — + В!!(и + 1) = О, дхз хт = О, 0 < х! < 1, (9) хз = !и 0 < х! < 1.
(10) Поставленная задача характеризуется параметрами !з, й, 8!е, В!м В!з. 13.2.2. Рааиоетиая схема Для численного решения задачи (2), (6), (7)-(10) используем разностную схему сквозного счета со сглаживанием коэффициентов, рассмотренную в п. 7.2. Вместо (7) используется уравнение ди д~и (! + $1е б(и, Ь)) — — ~~ — = О, х б й, ! > О, (11) дг, дх~ где б(и, Ь) аппраксимирует б-функцию. В приведенных ниже расчетах использовалась 1 — !и! <Ь, б(и, Ь) = (12) О, !в!>Ь, Для решения задачи (6), (7) — (12) на обычной равномерной разност- ной сетке в прямоугольнике й применялась линеаризованная разностная схема с весами: Ь(у„) "+' " + й(аун и + (! — а)у„) = рн, T хбй, н=0,1,.... (13) В (13) Ь(у) = ! + Все б(у, й ), (14) а разностный оператор й во внутренних узлах сетки задается следующим образом: — Е+ ай)унм = 1 — "Š— (1 — а)й у~+'рн ь(у„) ~ /ь(у.) (15) в=0,1, хбй, 2 Ау = — ~~~ уе,л., х б ы.
а=! Для граничных узлов используются аппроксимации, аналогичные приведенным в п. 13.1. Правая часть (13) отлична от нуля только в граничных узлах и обусловлена неоднородными граничными условиями (8)-(10). Для определения решения на верхнем слое из (13) получим сеточное эллиптическое уравнение 13.2. Затвврдввание расплава в полости прямоугольной формы бб5 БХХВКОХГПХЕ ГХХБ02 ( АО, А!, А2, У, Р ) С С С С С С С ИЗВ02 — Подпрограмма формирования системы линейных уравнений для раэностной задачи на каждом временном слое.
В!У) — Подпрограмма вычисления значений эффективной теплоемкости. 1МРХЛС1Т ЙЕАХ,ьЗ ( А — Н, Π— Е ) ХХХМЕХБ1ОХ АО(Х1,Х2), А1(Х1,Х2), А2(Х1,Х2), У(Х1,Х2), Г(Х1,1Ч2) СОММОХ / Т02НЭБ / Х1, Х2, Н1, Н2, В11, В12, Б1СМА, Б21, Б12, Ф БС, БС21, БС!2 Внутренние узлы РО 20 Х = 2, Х2 — 1 ОО 1О 1 = 2, !Ч1 — 1 ВН = В(У(1,1)) А1(1-1,1) = Б21 А!(1,1) = Б21 А2(1,Х вЂ” 1) = Б12 А2(1,1) = Б12 АО(1,Х) = А1(!-1 Х) + А1(1,1) + А2(1,Х вЂ” 1) + А2(1,Х) + В11 Г(1,Х) = БС21'( У(1+1,1)+У(1-1,1) — 2,130'У(11) ) + БС!2'( У(1,Х+1)+У(1,1-1)-2,00'У(1,Х) ) * + В11эУ(1,1) 10 СОМП1Ч11Е 20 СОМПХПЕ С С Левая граница: краевое условие второго рода С !ХО 30 Х = 2, !Ч2 — ! В1Х = 0.5ОО'В(У(1,Х)) А2(1,1-1) = 0.5130'Б12 Зта задача записывается на пятиточечном шаблоне в симметричном виде (см.
(!6) из и. !3.1). Формирование коэффициентов этой сеточной задачи проводится в подпрограмме РОБ02, по тексту которой легко прослеживаются используемые аппроксимации, в частности, на различных сторонах границы и в угловых точках расчетной области. Сеточная задача (15) на каждом шаге по времени решается итерационно. Ниже приведен текст этой подпрограммы. 666 Р~ава 13.
Примеры висленлого моделирования, А2(1,2) = 0.5РО*Б12 АО(1,2) = А1(1,2) + А2(1,! — 1) + А2(1,!) + В1! г(1,2) = БС21'( У(2,!)-У(1,!) ) + 0.5РОьБС12"( У(1,!+1)+У(1,! — 1)-2.РО'У(1,.!) ) + В1ХьУ(1,2) 30 СОХПХ!!Е С С Нижняя граница: краевое условие третьего рода С РО 40 1 = 2, Х! — 1 В11 = 0.5РО'В(У(1,1)) А1(1 — 1,!) = 0.5РОьБ21 А1(1,1) = 0.5РО Б21 АО(1,1) = А1(1- 1,1) + А1(1,1) + А2(1,1) + В11 + Н1ьВ11 Б10МА В(1,1) = 0.5РО'БС21'( У(1+1,1)+У(1 — 1,1) — 2.РО'У(1,1) ) + БС12'( У(1,2)-У(1„1) ) Ф + В11 "У(1,1) — Н1'ВИ'( БС"У(1,1)+1.РО ) 40 С01ЧТ!ХЦЕ С С Верхняя граница: краевое условие третьего рода С РО 50 1 = 2, Х1 — 1 В1Х2 = 0.5РО'В(У(1,Х2)) А!(1-1,Х2) = 0.5РОьБ21 А1(1,Х2) = 0.5РОьБ21 А2(1,Х2) = О.РО АО(1,Х2) = А1(1-1,Х2) + А!(1,Х2) + А2(1,Х2-1) + В!Х2 + Н1*В12ьБ!ОМА г(1,Х2) =0.5РО'БС21'( У(1+1,Х2)+У(1 — !,Х2)-2.РО'У(1,Х2) ) + БС12*( У(1,Х2-1) — У(1,Х2) ) + В1Х2'У(1,Х2) — Н1'В12'( БС'У(1,Х2)+!.РО ) 50 СОХПХбЕ С С Правая граница: краевое условие третьего рода С Во 60 2 = 2, Х2-1 ВХ12 = 0.5РО'В(У(Х1,!)) А1(Х1,2) = 0 РО А2(Х1,Х вЂ” 1) = 0.5РО'Б12 А2(Х1,!) = 0.5РО'Б12 АО(Х1,!) = А1(Х1 — 1,!) + А2(Х1,,! — !) + А2(Х1,!) + ВХ1! + Н2'В!1'Б!ОМА 13.2.
ЗатверР.'ванне расплава в полости прямоугольной формы 667 Р(Х1,3) = БС21'( У(Х1 — 1!)-У(Х1,2) ) + 0.5РО'БС12'( У(Х!,2+1)+У(Х1,1-1)-2.РО'У(Х1,!) ) + ВХ!геУ(Х1,2) — Н2'В!1'( БС*У(Х1Д+1.РО ) 60 СОХТ1ХРЕ С С Левый нижний угол С В11 = 0.25РО'В(У(1,1)) АО(1,1) = А1(1,1) + А2(1,1) + В11 + 0.5РО'Н!еВ!1*Б1ОМА В(1,1) = 0.5РО'БС21'( У(2,1) — У(1,!) ) + 0.5РО*БС12'( У(1,2)-У(1,1) ) + В11'У(1,1) — 0.5РО*Н!'В!1*( БС*У(1,1)+1.РО ) С С Левый верхний угол С В1Х2 = 0.25РО*В(У(1,Х2)) АО(1,Х2) = А1(1,Х2) + А2(1,Х2 — 1) + В1Х2 + 0.5РО'Н1*В12'Б!ОМА Г(1,Х2) = 0.5РО'БС21'( У(2„Х2) — У(1,Х2) ) + 0.5РО*БС12*( У(1,Х2 — 1) — У(1,Х2) ) + В1Х2'У(1,Х2) — 0.5РО'Н1*В12'( БС*У(1,Х2)+1.РО ) С С Правый нижний угол С ВХ11 = 0.25РО'В(У(Х1,1)) АО(Х1,1) = А1(Х1-1,1) + А2(Х1,1) + ВХ11 + 0.5РО*(Н1+Н2)'В11'Б1ОМА Г(Х1,1) = 0.5РО'БС21'( У(Х1-1,!) — У(Х1,1) ) + 0.5РО'БС12'( У(Х!,2)-У(Х1,!) + ВХ11'У(Х1,1) ) — 0,5РО'(Н!+Н2) В!!е( БС'У(Х1,1)+1.РО ) С С Правый верхний угол С ВХ1Х2 = 0.25РО*В(У(Х1,Х2)) АО(Х1,Х2) = А1(Х1-!,Х2) + А2(Х1,Х2-1) * + ВХ1Х2 + 0.5РО'( Н!'В12+Н2'В11 )'Б1ОМА г(Х1,Х2) = 0.5РО'БС21*( У(Х1 — 1,Х2)-У(Х1,Х2) ) + 0.5РО'БС12'( У(Х1,Х2 — 1)-У(Х1,Х2) ) Ф + ВХ1Х2*У(Х1,Х2) — 0.5РО'( Н1'В12+Н2'В!1 ) е '( БС'У(Х1,Х2)+1.РО ) 668 (дава 13.
Пршиеры численного моделирования, РКООКАМ ТЕЯТ02 С С ТЕЯТО2 — ззтвердеввние расплава в полости С прямоугольной формы. С 1МРЕ1С1Т КЕА1ь8 ( А — Н, 0 — Е ) С РАКАМЕТЕК ( Н31М = 50000 ) С Т)1МЕХЯ10Х А(101М) С С С С С С С С С Способ размещения указанных фрагментов в массиве А, С описываемый ниже„позволяет избежать в алюритме С ЯОЬЧЕ1 специального рассмотрения компонент векторов С и строк матриц, соответствующих граничным узлам С сеточного прямоугольника. С С Еелн Х вЂ” число компонент вектора неизвестных разностной С задачи, то размещение фрагментов описывается списком С эквивалентностей: С С С С С С 'С СОММОХ / Т02Н)Я / Х1, 1Ч2, Н1, Н2, В11, В12, Длина массива А должна быть достаточной для размещения коэффициентов симметричной матрицы ревностной задачи, заданной главной АО, соседней верхней А1 н удаленной верхней А2 днагоналямн, векторов решения У и правой части Р, в также векторов, участвующих в итерационном алгоритме решения ревностного уравнения (ем.
алгоритм ЯО(ЛЕ1). Е(2ШЧА(ЕХСЕ ( А(1), АО), * ( А(Х+2). А1 ), * ( А(2*Х+Х1+1), А2 ), * ( А(увХ+1), У ), * ( А(йеХ+1), Р ) 18.2.8. Программа Реализация лннсаризованной схемы (13) для приближенного решения задачи Стефана (2), (6)-(10) осушествлястся программой ТЕЯТ02. На каждом временном слое решается сеточная эллиптическая задача (14). Для этого используется подпрограмма Я01.ЧЕ1, реализующая попеременно-треугольный итерационный метод приближенной факторизации— сопряженных градиентов. Подробное описание алгоритма и подпрограммы Я01ЧЕ1 дано в п. 13.1. Б1ОМА, Б21, $12, Ф БС, БС21, БС!2 СОММОХ / ТО2В / ОЕЬ, ННТ, БННТЕ С С Ввод данных задачи: С С С С С С С С С С С С С С Х11ь Х2Ь вЂ” координаты левого нижнего угла расчетной прямоугольной области; Х1В, Х2К вЂ” коордииаты правого верхнего угла; В11, В12 — Числа БИО; БТŠ— Число Стефана; 11ЕЬ вЂ” полуширина сглаживаиия дельта-функции; Б10МА — весовой множитель в разностиой схеме; Х1, М2 — число узлов сетки по соответствующим направлениям; ЕРБ — требуемая относительная точность итерациовиого приближения к решению; ТАЮ вЂ” шаг по времени.
Х11. = О.ОО Х1К = 1.ОО Х21. = 0.130 Х2К = 1.ОО В11 = 1ОО.ОО В12 = 0.5ОО БТЕ = 0.25ОО ОЕЬ = 0.0400 Х1 = 51 Х2 = 51 ЕРБ = !.Π— 2 Б1ОМА = 1.ОО ЫО = 0.25ОО ТА11 = 0.00100 М = 25 ОРЕХ ( 01, Р1ЬЕ = 'В02.0АТ' ) ОРЕХ ( Об, Р1ЬЕ = 'К02.РТС' ) ВО!01 = 1, 1Р1М А(1) = О.ОО ' 10 СОХТ1ХЮЕ С Х = Х1*Х2 13.2. Затвердевапие расплава в полости прямоугольной формы 669 б70 1)О201= 1,Х А(7'Х+1) = [30 20 СОХТ)Х1)Е С ХТАЫ = 0 МТАЫ = 0 С ЗО СОХТ1ХСЕ С С С С С С С С С С С С С С С С С С Бгава 13. Примеры численного моделирования Н1 = (Х1К-Х1Ь) / (Х1 — 1) Н2 = (Х2К-Х21.) / (Х2 — 1) Б21 = Б1ОМА*Н2/Н1 812 = Б1ОМА'Н1/Н2 БС = 1.ЫΠ— Б1ОМА БС21 = БС*Н2/Н1 БС12 = БС'Н1/Н2 ННТ = Н1'Н2 / ТАЬ) БННТЕ = ННТ'( !.1)0+0.51)0'БТЕ/1)ЕЬ ) В подпрограмме И)302 определяются центральный АО, правый А1 и верхний А2 коэффициенты ревностной схемы на пятиточечном шаблоне и правая часть Р.