Вычислительная теплопередача Самарский А.А. Вабищевич П.Н. (1013606), страница 103
Текст из файла (страница 103)
Примеры численного моделирования,„ Система уравнений (12)-(14) дополняется граничными и начальными условиями, вытекающими из (4) — (8). Условия прилипания и непротекания (4) на границе единичного квадрата й записываются в виде !р(х,1)=0, хбдй, 0<1<Т, (15) — (х,1) = О, х б дй, О < 1 ( Т. д!Р (16) дп Начальное распределение скорости (условие (7)) дает у!(х, 0) = О, х Е й. (17) В безразмерных переменных и(О~хм!) = 1~ и(1~ х2~1) = 0~ (18) другие условия для температуры (см. (5) и (8)) остаются прежними. 13.4,3.
Рааиоетиая схема Для численного решения задачи конвекции (5), (8), (12)-(18) будем использовать безусловно устойчивую линеаризованную разностную схему из п.9.5. Для этого используется несколько необычная запись конвективных слагаемых. Зададим вектор н = (и!!, и!!) компонентами ди! д!о ги! = — и!! =-— (19) дх! ' дх! и запишем (12) в виде ди! г ! — -У( )р-~. —,-Ог — =О, д1 х-~ дх! дх! «=! (20) хай, 0<1(Т. У(й) = ~:У.(Ь), а=! ! Уа(да)х ((яи(х)я)Е + яа(х)хч ) (21) Введем в й прямоугольную равномерную сетку с шагами Ь! и Ь! по переменным х! и х! соответственно. На множестве сеточных функций Н, обращающихся в нуль на границе сетки, определим сеточный оператор Лапласа ! Лу=~~ Л у, Л у= — у;,, а=1,2.
а=! ' Для конвективных слагаемых используем аппроксимации второго порядка на основе (11): 697 13.4, 1ьонеекннн е ноеости квадратного сечение с определением д„а = 1,2 на основе центральных аппроксимаций соотношений (9) и (19). Определим сеточную функцию р(х) (граничное условие Тома для вихря скорости) соотношением 2 р(х) = ~~! р,(х ), а=1 О, Ра(ха) = 2 Р' а Ь,<х <1,— Ь„ !х = 1,2. ха ла> 1а Ла~ "— 1г( .И +Лго.
! +р(х)тр =От(в+), (23) Второй этап состоит в коррекции по конвективному переносу и граничному условию для вихря скорости: + (р(х) — У(нл)Ныл+! — галл) = 0 (24) т Приходим к несамосопряженной сеточной эллиптической задаче для определения функции тока (24). Для определения температуры во внутренних узлах сетки используется линеаризованное разностное уравнение ял,! - ял 1 + гг(тл)аль! ((вл-ь!)х,л, + (вл+!)е,л,) = О (25) т г С учетом граничных условий (5)„(18) аппроксимируется уравнение (!4) в граничных узлах сетки. При построении разностной схемы аппроксимация конвективных слагаемых в соответствии с (21) реализуется в подпрограмме т. При расчете конвективного переноса температуры используется функция тока гв' (см. (25)), а для конвективного переноса массы (см. (23), (24)) '!свовьтуетсл вихрь скорости в соответствии (19) Для решения системы уравнений (13), (14), (20) с соответствующими граничными и начальными условиями используется следующая линеаризованная схема.
Сначала рассчитывается вихрь скорости !ила!П вЂ” — ЛтР„+!П из уравнения " — У(нл)грл+Л грль!!!+ р(х)фл = Ог(яль!)г, (22) Тем самым, конвективный перенос, граничное условие для вихря берутся с предыдущего временного слоя. С учетом введенных обозначений (22) записывается в виде 698 И~ава 13. Примеры численного моделирования, Я)ВКО1)ПХЕ ч' ( О, А1Ь, А1К, А2Ь, А2К ) С С С С С С С С В подпрограмме ч вычисляются четыре диагоналя: А1Ь, А1В, АЗЬ.
АДВ матрицы ревностного оператора коивективного переноса в уравнении для температуры (по ааданным нянчениям функции тока (1 РБ1) или в уравнения для функции тока (по ааданным значениям вихря скорости В - ОМС). 1МРЫС1Т КЕАЬч8 ( А — Н, О-Е ) Р1МЕХБ1ОХ 0(Х1,Х2), А1Ь(Х1,Х2), А1К(Х1,Х2), А2Ь,(Х!,Х2), А2К(Х1,Х2) СОММОХ / Т04 / Н181, Н281, Н1412, Н2412, НН81, ТА(31, РК1, ОК2Н1, Х1, Х2 РО 20 1 = 2, Х2-1 РО 10 1 = 2, Х1 — 1 А11.(1,1) = — НН81'( 0(1-1.!+Ц 0(1 1 ) Ц + 0(1,1+Ц вЂ” 0(1,1-Ц ) А1К(1,1) = НН8!'( 0(1+1,1+Ц вЂ” 0(1+1,1- Ц + 0(1,2+Ц вЂ” 0(1,1-Ц ) А2Ь(1,1) = — НН81'( 0(1 1,.1 — Ц вЂ” 0(1+1,1 — Ц + 0(1-1,1) — 0(1+1,Ц ) А2К(1„1) = НН8!'( 0(1-1,1+Ц 0(1+1,.1+Ц + 0(1-1,1) — 0(1+1,1) ) 10 СОЯПХЦЕ 20 СОХПХЮЕ Формирование сеточной несамосопряженной задачи для определения температуры (разностное уравнение (25) с соответствующими граничными условиями) на новом временном слое проводится в подпрограмме РРЯЬ Разностная задача записывается на обычном пятиточечном шаблоне в виде -А2!гу,у(4, — АИ! ~уу; ~4+ А01; у!.— — А)ге уг, ~4 — А2гбу;4+~ — — Л~, ! =1,2,...,)Ун у = 1,2,...,дгп (26) В подпрограмме проводится вычисление коэффициентов разностного .у равнения.
699 13.4. Копввкдил в полости квадратпгмо сечения Я1ВКОЦТ!ХЕ Р!лЯ1 ( АОО, А11., А1К, А21., А2К, 11, Р ) В подпрограмме ИВС вычнеляютея пять диагоналей: АОО, А1Ь, А1В, А2Ь, А2В, несимметричной матрицы и правая часть У разноетной схемы уравнения для температуры. В массиве У передаатея еначенне температуры на предыдущем слое по времени. на входе в массивах А1Ь, А1В, А21., А2В передаютея коэффициенты ревностного оператора конвективного переноса при эаданных значениях функции тока. 1МРЫСГТ ВЕАЬе8 ( А-Н, 0 — Е ) Р1МЕХЕ1ОХ АОО(Х1,Х2), А11.(Х1,Х2), А1К(Х1,Х2), Ф А2ЦХ1,Х2), А2К(Х1,Х2), Ф ЩХ1,Х2), Р(Х1,Х2) СОММОХ / Т04 / Н181, Н2Я, Н1412, Н2412, НН81, Ф ТАШ, РК1, ОК2Н!, Ф Х1, Х2 С С Внутренние уалы С !30 20 1 = 2, Х2-1 !30 10 1 = 2, Х1-1 АОО(1,1) = ТАЫ! + РК!е( Н1Я + Н1$1 + Н2Я + Н281 ) А11.(1,!) = РКРН1Я вЂ” А1Ц1,!) А1К(1,!) = РК!"Н1$1 — А1Щ1„1) А2Ь(1,!) = РК1еН2Я вЂ” А2Ц1„!) А2Щ1,1) = РК1еН2$1 — А2К(1,!) Р(1,1) = ТАШ'!1(1,!) 10 СОХТП~ПЗЕ 20 СОХТ1ХЫЕ С С Нижняя граница, краевое условие второго рода С 130 30 1 = 2, Х1-1 АОО(1,1) = 0.500'ТАШ + РК1е( Н1$1 + Н2$1 ) А11.(1,1) = 0,5!10еРК!еН1Я А1Щ1,1) = 0.500еРК!"Н1Я А2К(1,1) = РК1'Н2Я Р(1,1) = 0.5110'ТА!1!ЧЗ(1,1) 30 СОХТ!Х!1Е о Игала 13.
Примеры численного моделирования С Верхняя граница, краевое условие второго рода С РО 40! = 2, Х1-1 АОО(1,Х2) = 0.5РО'ТАЫ1 + РК1ч( Н1Я + Н2Я ) А1Ь(1,Х2) = 0.5РОчРК!'Н1Я А1К(1,Х2) = 0.5РО'РК1чН1$1 А2Ь(1,Х2) = РК1'Н2$1 г(1,Х2) = 0.5РО'ТАЫ1ч11(1,Х2) 40 СОХП1~П1Е С С Левая граккца, краевое условие первого рода С РО501=1, Х2 Ы(1,1) = 1.РО АОО(1,3) = 1.РО А1К(1,1) = О.РО г(2,1) = г(2,1) + А1Ь(2,1) А1Ц2,1) = О.РО А2Ц1,1) = О.РО А2К(1,1) = О.РО г(1,1) = 1.РО 50 СОХПХ11Е С С Правая граница, краевое условие первого рода С РО601=1,Х2 Ы(Х1,1) = О.РО АОО(Х1,1) = 1.РО А11.(Х1,1) = О.РО А1К(Х1 — 1,1) = О.РО А2ЦХ1,1) = О.РО А2К(Х1,1) = О.РО г(Х1,1) = О.РО бО СОХПХ15Е С С Формирование разностной задачи для определения вихря скорости на промежуточном слое по времени проводится в соответствии с (23) в подпрограмме рРЯОМН чала~отсе коэффициента ггто л1ог ', .ллг 701 13.4. Конввкяин в полости квадратного сечения БПВКОЮТ1ХЕ РПБОМСг ( АО, А1, А2, ОМб, Р, Ф АОО, А11., А1К, А2Ь, А2К, РЯ, 11 ) С С С С С С С С С С С С С С В подпрограмме УЭЗОМО вычисляются тря диагонали: АО, А1, А2, симметричной матрицы и правая часть Р рааностной схемы уравнения для вихря скорости.
Массиву АОО прнсваиваютея эначения сеточной функции т(ОМО)РВ1, выступающей в качестве одного нз слагаемых в правых частях уравнений для вихри скорости и функции тока. В массивах ОМО, Р31 н 1Г передаются аиаченна соответственно вихри скорости и функции тока с предыдущего слоя по времеви и температуры с текущего слоя. В массивах А1Ь, А1В, А2Ь. А2В передаются коэффициенты равностного оператора конвектнвного переноса прн еадаиных аначеннях вихря. 1МРЫС1Т КЕАЬе8 ( А — Н, Π— Е ) 0!МЕХИ!ОХ АО(Х1,Х2), А1(Х1,Х2), А2(Х1,Х2), Ф ОМО(Х1,Х2), Р(Х1,Х2), АОО(Х1,Х2)„А1Ь(Х1,Х2), А1К(Х1„Х2), Ф А21.(Х1,Х2), А2К(Х1,Х2), е РБ1(Х1,Х2), Ь1(Х1„Х2) С СОММОХ / Т04 / Н1Я, Н281, Н1412„Н2412, НН81, » ТАШ, РК1, бК2Н1, Ф Х1, Х2 С С Внутренние уааы С 90 20 1 = 2, Х2-1 ВО 10 1 = 2, Х1-1 АО(1,1) = ТАШ + Н181 + Н1Я + Н281 + Н2Б1 А1(1,!) = Н1Я А2(1,!) = Н2Б! АОО(1,!) = А1Ь(1,!)"РБ1(1 — 1,!) + А1К(1,!)еРБ1(1+1,Х) + А2Ь(1,!)'РБ1(1,1-1) + А2К(1,1)"РБ1(1,!+1) Р(1,1) = АОО(1,!) + ТАШ'ОМьт(1,!) Ф + аК2Н!Ф( Ц(!+1,!) — 11(1 — 1,1) ) с самосопряженным оператором на пятнточечном шаблоне (см.
п. 13.1). Предварительно рассчитывается конвективный перенос с помощью описанной ранее подпрограммы У. Текст подпрограммы снабжен достаточно подробными комментариями и в дополнительных пояснениях не нуждается. УОЗ 13.4. Конвенция в тыости квадратного сечения Левый ияжний угол АО(1,1) = 1.ОО А1(1,1) = 0 ОО А2(1„1) = 0.1Ю Р(1,1) = ОМО(1,1) С С С Правый нижний угол АО(Х!,1) = 1.ОО А2(Х1,1) = О.ЕЮ Р(Х1,1) = ОМО(Х1,1) С С С Левый верхний угол АО(1,Х2) = 1.ОО А1(1,Х2) = 0.1Ю Р(1,Х2) = ОМО(1,Х2) С С С Правый верхняй угол АО(Х1,Х2) = 1.ОО Р(Х1,Х2) = ОМО(Х1,Х2) Подпрограмма ИЭБРЯ1 предназначена для формирования пятндиагональной сеточной задачи для определения функции тока на новом временном слое и реализует полушаг коррекции (24). Сеточная задача снова записывается в виде (2б).
С С С С С С С С С С Я1ВК01УТНЧЕ НЭБР31 ( АОО, А1Ь, А1К, А21., А2К, е РБ1, Р, ОМО ) В подпрограмме Н13Р31 вычисляются пать диагоналей: АОО, А1Ь, А1К, А2Ь, А2К, несимметричной матрацы н правая часть Р разностной схемы уравнения для функции тока. В массивах Р31 н ОМО передаются значения соответственно фующии тока и вихря скорости на предыдущем слое по времени, в массиве АОО передается величина Ч(ОМО)РВЬ вычисленная ранее в подпрограмме И13ОМО.
33!ава 13. Примеры численного моделирования 1МР3.!С1Т КЕА3 е8 ( А-Н, О-г. ) 33!МЕХЕ!ОХ АОО(Х1,Х2), А1ЦХ1,Х2), А1К(Х1,Х2), Ф А2ЦХ1,Х2), А2К(Х1,Х2), Р81(Х1,Х2), Р(Х1,Х2), ОМчл(Х1,Х2) СОММОХ / Т04 / Н181, Н281, Н1412, Н24!2, НН81, Ф ТАШ, РК1, чзК2Н!, Х1, Х2 С С Внутренние узлы С ЭО 20 3 = 2, Х2-1 130 10 ! = 2, Х! — 1 Р(1,3) = ТА!3!'ОМб(1,3) — АОО(1,3) АОО(1,3) = ТА!3!е( Н1Б1 + Н181 + Н2Я + Н2Я ) А1Ц1,3) = ТАШ'Н181 + А1Ц1,3) А1К(1,3) = ТАШеН181 + А1К(1,3) А21.(1,3) = ТА!3!'Н28! + А2Ц1,3) А2К(1,3) = ТАШеН281 + А2К(1,3) 10 СОЮПХЦЕ 20 СОХТЗХ13Е С С Левая граница, краевое условие первого рода С 330 30 3 = 2, Х2-1 РЯ(1,3) = 0.130 АОО(1,3) = 1.330 А1К(1,3) = 0.330 А23.(1,3) = О.!30 А2К(1,3) = 0.330 Р(1,3) = 0.130 АОО(2,3) = АОО(2,3) + Н!412 А13.(2„3) = 0.130 Р(2,3) = Р(2,3) + Н1412'РБ!(2,3) 30 СОМПХ13Е С С Правая граница, краевое условие первого рода С 33О 40 3 = 2, Х2-1 Р81(Х1,3) = 0.130 АОО(Х1,3) = 1.ОО А11.(Х1,3) = 0.130 А2ЦХ1,3) = 0.330 13.4, яонеекцил е лолосжи квадро голого сечения 705 А2К(Х1,1) = О.ОО Р(Х1,1) = 0.330 АОО(Х1 — 1,1) = АОО(Х1 — 1,1) + Н1412 А1К(Х1 — 1,1) = 0.00 Р(Х1 — 1,1) = Р(Х1 — 1,1) + Н1412"РБ!(Х1 — 1,1) 40 СОХПХ!1Е С С Нижняя граняца, краевое условие первого рода С 00501=2,Х1 — 1 Р81(1,1) = О.ЭО АОО(!,1) = ЬРО А1Ц!,1) = О.ОО А1К(1,1) = ОЛ)0 А2К(1,1) = 0.00 Р(1,1) = О.ОО АОО(1,2) = АОО(1,2) + Н2412 А2Ь(1,2) = О.ОО Р(1,2) = Р(1,2) + Н2412*РБ!(1,2) 50 СОХПХ!1Е С С Верхняя граница, краевое условие первого рода С РО 60 1 = 2, Х1-1 РЯ1(1,Х2) = О.ОО АОО(1,Х2) = 1.ОО А1Ь(1,Х2) = О.ОО А1К(1,Х2) = 0.00 А2Ц1,Х2) = О.ОО Р(1,Х2) = О.РО АОО(1,Х2 — 1) = АОО(1,Х2-1) + Н2412 А2К(!,Х2 — 1) = О.ОО Р(!,Х2-1) = Р(1,Х2 — 1) + Н2412'РЯ1(1,Х2-1) 60 СО1ЧТ1ХЦЕ С С Левый нижний угол С РЯ(1,1) = О.ОО АОО(1,1) = 1.ОО А1К(1,1) = О.ОО А2К(1,1) = 0.110 Р(1,1) = О.ОО 706 Пгава 13.
Примеры численного моделирования Правый нижний угол РЯ(Х1,1) = О.ОО АОО(Х1,1) = 1.ОО А11.(Х1,1) = О.ОО А2К(Х1,1) = О.ОО Р(Х1,1) = О.ОО Левый верхний угол Р31(1,Х2) = О.ОО АОО(1,Х2) = 1.ОО А1К(1,Х2) = О.ОО А21.(1,Х2) = О.ОО Р(1,Х2) = О.ОО Правый верхний угол Р$1(Х1,Х2) = О.ОО АОО(Х1,Х2) = 1Ш АН (Х1,Х2) = О.ОО А21,(Х1,Х2) = О,ОО Р(Х1,Х2) = О.ОО КЕПЛ~Х ЕХР По найденной функции тока численным дифференцированием находится вихрь скорости. Для этого используется подпрограмма ОХОРБ1. $11ВКОЦТ1ХЕ ОМОР31 ( ОМО, Р31 ) В подпрограмме ОМОР31 расчнтывается значение вихря ОМО на текущем слое по значениям функция тока РЯ.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
