А.А. Самарский - Теория разностных схем (1989) (1160092), страница 69
Текст из файла (страница 69)
В неоднородной среде й, с, ~ могут быть разрывными функцииями х и Г (для разных веществ зависимость й, с, ~ от температуры и может быть различной). Типичным является случай, когда функции й = й(и), с = с(и), ~(и) зависят только от температуры и: с(и) д д (х(и)д ) +~(и (14) При этом итерации сходятся со скоростью геометрической прогрессии со знаменателем д,=й, так что1 о ((~рз о1(Чв )ог Так как 0(1, то итерации сходятся прн любом выборе начального приближения у. 2. Квазилииейное уравнение теплопроводности. В предыдущих главах мы рассматривали лишь линейное уравнение теплопроводности.
Между тем, для высокотемпературных процессов, протекающих, например, в плазме, коэффициент теплопроводности является нелинейной функцией температуры (и плотности), а в ряде задач, кроме того, функцией градиента температуры. Далее, источники тепла (правые части в уравнении теплопроводности) могут зависеть от температуры, если, например, тепло выделяется в результате химической реакции. От температуры может зависеть и теплоемкость среды. Таким обрааом, мы приходим к нелинейному уравнению теплопроводности (12) е ь ивлзилинввнои тглвнвнив тэлпниэгонодности 417 Вводя новую переменную и = ') й (й) оф, приведем уравнение е И4) к виду где <р(и) = ) с($) а$.
е и Если положить и= ) с($)о(3, то вместо И4) получим урево ионне 3= ( ()Р)+7(), так что ~х(и)йи = ~й(и)г)и. о е Весьма часто с(а) и й(а) являются степенными функциями температуры: с(п) сопи, й йои'. иа+~ Вводя в этом случае и= ) с($)Щ=сеа-~- и учитывая, что Ф-а э-а ди а ди " С-ади " /а+1'1а+г а+о ди ди ад о ди о ~ и ) ди" е преобразуем уравнение И4) к виду ди д I оди1 ' б — а ао/а+1'1а+г — = — (хои -)+~(и), с= —, хо=-~ — ~ до ди1 о д*3 а+2 о (, о ) 3. Некоторые аналитические решения квазилииейного уравнения теплопроводностн.
Нелинейность коэффициента теплопроводности приводит к новым физическим эффектам, главный из которых — конечная скорость распространения тепла. Убедимся в атом, найдя простейшие частные решения уравнения —, = д ~хои — ~)> хо>0, с>0, х>0. (15) Пусть при х 0 задана температура и и,'$" Иб) и требуется найти решение уравнения И5) в области л> О, г> О, 27 л. л. соиоиоииз 418 гл. шп. мвтоды гвпшния нвлннвиных ттавнвнии предполагая, что начальная температура равна нулю: и(х, 0) ~0.
ИП Будем искать решение задачи в виде бегущей волны и(х, 1) =У(Рр — х), Р = сопзЦ где Щ) — неиавестная функция.. Подставляя это вырвжение в уравнение И5) и учитывая, что ди дУ де д(7 — =Р—, — = —— де д$' дх = дВ получаем для Щ) обыкновенное дифференциальное уравнение Ри'- (яра Г)'. Отсюда находим ррр(7рУ Р0+ сопзФ. Полагая совз1 ° О, будем иметь х (Пе)' кеР'еУ =Ж7 или е =1. Во Интегрируя еще раз, получаем — У вЂ” 5+с .
Так как 0=0 при Ф х=О (4=0), то е,=О, Юо Пе / дро ~пе (е и= ц(з) = ( — В) = ~ — ) с''~1 — -Р,-) Отсюда следует, что при х = 0 и(0, р) = Сравнивая с Иб), находим Р~о и= — — =и . Ф о рре е ° Таким образом, задача И5) — И7) имеет решение в виде бегущей волны и(х, р) = иер' (1 — — ) ррр 3 е (РР. )ъ!е О, » (РР (18) и (х, г) = 0 при х)~ Рр, если выполнено условие и = 1/о; при этом скорость волны 6 ь квазилннвинов ввмвнвнив тнплотюонодности 419 определяется параметрами к„ о и и;. Р = у' к,из/а, Отсюда видно, что на фронте температурной волны х =Р1 температура и тепловой поток равны нулю при о)0, а производная ди дх <ц)1/О (щ и)1-1/и обращается в бесконечность при о~1, конечна при о 1 и равна нулю прн 0(о < 1.
Поэтому при а) 1 можно говорить лишь об обобщенном решении уравнения теплопроводности (15). Причина конечной скорости Р фронта — нелинейная зависимость коэффициента теплопроводности от температуры. Из формулы для Р видно, что для линейной теплопроводности, когда о = О, формально получаем Р =, т. е. скорость распространения тепла бесконечна. Возможен случай, когда температурный фронт неподвижен, т. е. Р = О. Такое решение существует при специальном граничном режиме вида и(0,0=-н, т>0, (19) ' р,— в)'" ' где 1,— произвольная постоянная, если предположить, что начальное услевие задано при г =— п(х, — ) =О. (20) Будем искать решение уравнения (15) методом разделения переменных, полагая н(х, Ф) о(х)Т(г).
После подстановки этого выражения в (15) и разделения переменных получим д иди 1 дг где А — параметр разделения. Отсюда находим —, = ).Т~+'. (21) (22) и является конечной. Решение вида (18) нааывается температурной иолной. Найдем тепловой поток ии+г — — (Рз — х) ' = мене /(оР) и (х, Ф). аФ+пи 420 ..
д шн и Будем искать решение уравнения (21) в виде о'= д,(х, — х)в где а и р — неизвестные пока числа, х,— произвольное число. Подставим о' в уравнение (21): в+ — () г В 1 к а ™ — [ -+() — 1)(х,— х) — Ха (х,— х) = О. з/а+с-2 1/О Е/е о ( о Отсюда находим )р взв (З+ о) Проинтегрируем теперь уравнение для Т: Т(в) (ок(с, — в))-'", где с,— постоянная интегрирования. В результате мы приходим к функции в/ ( )в1в/а .с;в-.в)те-(")'"[',' ) вв Сравнивая с граннчяым режимом (19), находим ю = 1/и, св = гю ив = и, следовательно, х,' = 2кв(о+ 2) ив/о.
Таким образом, уравнение (15) с граничным режимом (19) имеет решение / 1 — х/е ~в/е и(х, в) = ив[ ') прн О<х(хд, (23) '~« -) в и (х, в) = 0 при х) х„ где х, — ширина области прогрева. Фронт температурной волны неподвижен, таи как х, = сопев не зависит от в и зависит только от параметров задачи к„ о, и,. На фронте тепловой поток и температура обращаются в нуль при любом о) О, а производная ди/дх= прн о) 2 (на фронте бегущей волны ди/дх при о ~ 1). Решение (23) типа естоячей волныэ существует при с(в„ что связано с типом граничного режима (19) (который называется режимом с обострением), Для уравнения с тепловым источником, зависящим от температуры по степенному закону В 1 квазнлинвинов звавнвнив ткппопговодности 421 существуют решения как типа бегущей волны (при ~( о+ 1), так и типа стоячей волны (при () о+ 1).
Численное решение таких задач по схемам скеозного счета представляет значительные трудности из-за нелинейности н обращения в бесконечность производной на фронте температурной волны. Найденные в этом пункте точные решения являются хорошими тестами для проверки точности рааностных схем.
4. Разностная схема. Метод Ньютона. Перейдем к написанию разностных схем для квазилинейвых уравнений теплопроводности. Использование явных схем нецелесообразно, если Ь(и), с(и), 7(в) являются быстроменяющимися (напрнмер, степенными) функциями температуры. Условие устойчивости, явной схемы т 1 ш(ос(и) 2 шаха(и) требует мелкого шага по времеви, определяемого часто значениями функций Ь, с в небольшом числе узлов.
Поэтому применяются безусловно устойчивые неявные схемы. Рассмотрим сначала уравнение зс (25) Используем для его решения итерационный метод Ньютона ь В ай+1 з~ %+1 Ч(У)+Ч'(У)(У вЂ” У)-ту- =ЧЬ') В+1 Для определения отсюда у при граничных условиях а+1 з+г Ус = р1(гс+дс Ув = рз(гс+г) (27) (28) е краевыми условиями и(х,'6) в,(х), п(0, П = р,((), в(1, с) д,(с). Для его решения используем нелинейную относительно у'+' разностную схему ' =У-'~~с х=х = (Ь, 0<(<К, ЬХ=1. (26) г3+и (л Предположим, что ф'(у) > с,) О, !<р" (у)! < е,. Тогда нетрудно убедиться в устойчивости этой схемы и сходимости в С со скоростью 0(т + Ь'). Однако доказательство этих фактов весьма ' громоздко, и мы не будем его приводить.
Для определения решения у'+' на новом слое мы имеем нелинейное уравнение Р(У"') — У'-" = Р(У'). 422 а'л тпг мвтоды гвйпвния нллннвиньхх твйвнглин можно пользоваться прогонкой, которая устойчива при Ф (У»О. Это видно из самого уравнения, если переписать его в виде л +1 /, 2тй + т +1 к' "' ' ~Ф (У1~+ ьв/ "' + ьч "'+' 1 = 1, 21..., /Ч вЂ” 1, й й й й где ~' = Ф(у) — Ф(у) — Ф'(у) у, у = у'. Оценим скорость сходимости итераций. Для этого введем разность А+1 А+1 = У1 — У» /+ й А й+1 й+1 где у1 = у1~ . Подставляя в (27) у1 = У1+ сь уй =- у1+ с;, по+1 лучим й АА1 й+1 А й й Ф'(Ч) с — т с;„== 1р(у) — 1р(у) + 1р'(у) и+ ту-„„= й А А = Ф (У) Ф (У) + 1Р (Ч) (Ч Ч) Учитывая затем, что й й А А Ф (у) = Ф(у)+ Ф'(у) (» — у) + 0.5Ф'(ч) (у — у)', и, следовательно, А й А Ф (У) — 1Р (ч) + Ф (у) (ч — у) = Ф" (у) сй/2, й й где У=-У+Ос, 0< О <1, приходим к следующ4му уравнению: й й+1 й+1 й й 1р'(у) г — т с- =1р" (у)ий/2=с", т= 1/А, 0(1()Ч, (29) с однородными краевыми условиями А+1 й+1 с =О, сл= — О.
(30) Пользуясь принципом максимума, получаем оценку А+1 А . й й Нсбс<~0.5|!Ф" (у)/Ф'(у) )с1РНс~(Ч64с, (31) где Ч =О 51Ф" (у)/Ф'(у)1с < О 511р" (у)1с/с1 <О 5с /с, = д„так как 1р (у) ~ с, >О, 1Ф" (УП <с,. Отсюда следует, что для сходимости итераций по квадратичному закону достаточно, чтобы начальное приближение удовлет- 424 где у»= у» у» = у»», а»(и) = а(и» и и»), например, а» (»)) = 0,5 [й (»)» 1) + й (()»)), в (»<)=й( ~ ), 2й (и» ) )< (»<») "("< =1<-')+<»<' (36) (37) (38) От способа вычисления а<(»<) силъно зависит точность расчета температурной волны (см.