Fletcher-1-rus (1185917), страница 58
Текст из файла (страница 58)
Однако для достижения суммарной ошибки аппроксимации порядка 0(Л/т) необходимо для промежуточного решения Т вводить граничные значения так, чтобы они были совместимы с алгоритмами для внутренних точек (8.14) и (8.15). Например, если ставятся граничные условия Дирихле, то задание Тих, а =Ьа+'" при х =! приводит к алгоритму с ошибкой аппроксимации порядка 0(Л1). Чтобы ошибка аппроксимации .получила порядок 0(Лгз), необходимо задать Т'„а с помощью 8ЗО Гз.
8. Многомерное уравнение диффузии соотношений (8.12) и (8.13), т. е. по формуле Тнх, е = 0.5 (Ье + Ьеи ) — 0 25 А!Еде (Й вЂ” бе). (8. 17) Аналогичная проблема возникает и при использовании схем с приближенной факторизацией; надлежащие формы задания граничных условий обсуждаются в п. 8.3.2 и в 8.4. Таким образом, можно прийти к выводу, что схема НПН для двумерного случая обладает требуемыми свойствами: безусловно устойчива, имеет второй порядок точности и допускает экономичный метод решения. Схема НПН распространяется и на трехмерное уравнение диффузии, при решении которого вместо алгоритмов (8.14) и (8.!5) необходимо вести счет на трех дробных шагах, каждый из которых занимает время Аг/3.
В случае трех измерений схема НПН оказывается экономичной, имеет второй порядок точности по пространству, но является лишь условно устойчивой. Для ее устойчивости необходимо, чтобы з„, з„, ае « ' 1.5, где ее = а,йг(аа'. 8.2.2. Обобщенная двухслойная схема Попытаемся теперь дать обобщение идеи расщепления. Неявная двухслойная конечно-разностная схема общего вида применительно к уравнению (8.1) может быть записана в виде е-~-1 ат е и — — (1 — й)ь',а,.Е, Т! е+ а„ЕнвТ! в)— — 8 (а„1.„Т!г;ве' + а1,„вТ!", е) = О, (8.18) где Иг,+,'=Т",,", — Т,", Здесь величину ЬТ!, е можно рассматривать как поправку к решению на временнбм слое л, требуемую для продвижения на временнбй слой (н + 1).
Для минимизации накапливающейся ошибки округления полезно, чтобы величина АТ!,+е' фигурировала в явном виде в тексте компьютерной программы. Роль коэффициента р в уравнении (8.18) состоит в задании определенных весовых множителей при членах, соответствующих временным слоям п и (я+ 1); та же идея была использована н при выводе уравнения (7.24). С помощью перестановки членов уравнения (8.18) можно получить неявный алгоритм для определения ЬТ!,е. Вначале член Т!",+» разлагается в ряд Тейлора в окрестности а-го временнбго слоя в соответствии с формулой (3.16), а именно Тг, е =Т";, е+ Лг~ в! ~ + 0.5К! ~вге~ + !, е !. е $3.2.
Методы расщепления для многомерных задач 331 что можно аппроксимировать формулой Т ат" (8.19) Подстановка формулы (8.19) в уравнение (8.18) дает я+1 ат ч и т г и+! и+!т — (а 1„„Т! я+ ацЕццТ! е) — (1(а Е„хЬТ! е+ацЕццЬТ; е )=О, (8.20) или после некоторой перегруппировки !1 — 8 Ь1 (аяЕлл + ац1 цц)) ЬТг,+а' = Ж (ахйхл + азйцц) Т7, е. (8 21) В левой части соотношения (8.21) фигурируют алгебраические операторы, действующие как в одном, так и в другом направлениях.
Чтобы можно было воспользоваться алгоритмом Томаса, соотношение (8.21) заменяется другим соотношением, получаемым в результате приближенной факторизации, а именно Как показывает сравнение соотношений (8.21) и (8.22), в ле- вой части последнего содержится дополнительный член 3 2 ч+! 8 М ахацЕхцЕццЛТ!,ц. Отсюда следует, что соотношение (8.22) аппраксимирует (8.21) с точностью 0(Л(з). На каждом шаге по времени выполнение соотношения (8.22) может быть реализовано в форме двухэтапного алгоритма. Первый этап сводится к решению следующей системы уравнений, справедливых на каждой из сеточных линий, направленных параллельно оси х (лииии постоянного й на рис. 8.2): (1 — 8 Ца,Е ) ЬТгд ц = Ж (а„Езх + азЕцц) Тг, ь (8.23) Это уравнение служит для определения ЬТгп д что можно рассматривать как промежуточное приближение для Ы'!,ц~'.
Если Е,„ представляет собой трехточечный оператор с центральной разностью, то уравнения (8.23) будут иметь форму трехдиагональной системы, которая может быть эффективно решена с помощью алгоритма Томаса (п. 6.2.2). На втором этапе решается следующая система уравнений: (1 — б МацЕцц) ЬТг,"к' = ЛТгд ь (8.24) (1 — бала,Е„)(1 — 8ИацЕцц)ЛТ!",е' =М(а,Е„, +ацЕц ) Тг, ь (8.22) Гл. 8. Многомерное уравнение диффузии справедливая на каждой из сеточных линий, параллельных оси у (линии постоянного 1' на рис. 8.2).
Структура уравнений (8.23) и (8.24) аналогична структуре уравнений (8.!4) и (8.15), характеризующих схему НПН. В случае уравнений более сложной формы основной вклад в суммарное время исполнения алгоритма вносит расчет всех пространственных членов в правых частях уравнений (8.23). Реализация, рассматриваемая здесь, требует лишь однократного расчета пространственных членов на одном шаге по времени. В отличие от этого схема НПН требует проведения двух подобных расчетов.
Алгоритм (8.23), (8.24) предложен в работе [Ровд!аз, бнпп, 1964]. Двухэтапный алгоритм (8.23), (8.24) является безусловно устойчивым при р ) 0.5 и имеет ошибку аппроксимации порядка О(Лге, Лх', Лу'), если р=0.5. Построение, связанное с приближенной факторизацией, обобщается и на случай трех измерений, причем в отличие от схемы НПН соответствующая схема безусловно устойчива при р ) 0.5. Если изложенная здесь методика используется для построения решений стационарных задач в форме предельного стационарного состояния решений для соответствующих нестационарных задач ($6.4), то полезно в уравнение (8.23) ввести определение ЙНЯ =(а„1лл+ а„7еи)Ту е, (8.25г.
По мере приближения к стационарному решению величина ВНБ стремится к нулю; следовательно, контроль за изменением КН5 даст возможность установить близость к стационарному решению. 8.2.3. Обобщенная трехслойная схема Обобщенная трехслойная неявная схема для решения одйомерного уравнения диффузии задавалась в форме (7.25). Соответствующая трехслойная схема для случая двух измерений записывается в виде (1+ т) лт",+„' т дт,", — = (1 — 6) (а,7.
+ аи7-ен)Тг", е + + 6 (ал7.„л + ад1 „и) Т1, е (8.26) где Если применить то же самое построение, которое использовалось для перехода от соотношения (8.18) к уравнениям (8.23) и (8.24), то из (8.26) получится следующий двухэтапный ал- $8.3. Схемы расщепления горитм. На первом этапе соотношение (1 —, + Ь!ая7.„„) ЬТь я = М е т е = ( + ! (а,Е „+а„7.„я)Т! я+ + ЬТ~ е (8.27) дает возможность получить трехдиагональную систему уравнений, относящихся к каждой сеточной линии, параллельной оси х. На втором этапе каждого шага по времени используется следующее уравнение: 1 — + ЬГая7.„„~~ ЬТ, = ЬТ ( !!+ ! я+1 !! + т) (8.28) 9 8.3. Схемы расщепления и метод конечных элементов Здесь мы будем применять метод Галеркина с конечными элементами (9 5.3) к решению двумерного уравнения диффузии (8.1) с граничными и начальными условиями, заданными в форме (8.2) и (8.3), и попытаемся определить, нуждаются ли в модификации схемы расщепления, разработанные в $2.8, чтобы охватить конечно-элементную разновидность дискретизированных уравнений.
Воспользуемся элементами прямоугольной формы с билинейными интерполяционными функциями типа (5.59) внутри каждого элемента. Если метод Галеркина с конечными элементами применяется на сетке, однородной в При выборе специальных значений р = 1, у = 0.5 двухэтапный алгоритм, заданный с помощью уравнений (8.27) и (8.28), согласуется с уравнением (8.1) с ошибкой аппроксимации 0(Ь!е, Лх', Ьу') и является безусловно устойчивым. При реализации первого шага по времени (и = О) необходимо воспользоваться двухслойной схемой, подобной (8.23), (8.24). Рассмотренные в этом параграфе схемы расщепления естественным образом обобшаются и на случай трех измерений [М!!сне!1, Пг!!!!йз, !980]. Современный подход к процессам расщепления, при котором неявное уравнение дополняется некоторым членом, как правило, порядка 0(ЬР) с целью осушествлеиия факторизации, подробно обсуждается в работе [Понг!ау, 1977).
При решении уравнения диффузии в двух или трех измерениях имеется возможность применения схем расщепления повышенного порядка. Некоторые из таких схем рассматриваются в книге [М!!сне!1, Пг!!!!!йз, 1980). 334 Гн. 8. Многомерное уравнение диффузии где знак Э обозначает тензорное (или внешнее) произведение (Мазе, 19711; М, и Ми — массовые операторы соответствующих направлений, а Е„„и Е„„— разностные операторы тех же направлений (см. приложение А.2). Массовые операторы по направлениям имеют вид М (1/6 2/3 1/6) Ми (1/6 2/3 1/6)г (8 30) а разностные операторы по направлениям — вид (8.31) Таким образом, член Ма Э /.„,Тп» дает девятиточечное представление производной д'Т/дх'.
Обращаясь к рис. 8.3, получим М Э/. т = — 1( '-"+' ',"'+ "'"' ~+ И ахб» 6( Ьх» 2 ( Т...-2Т,, +Т,+,, ) +61 а* 1+ г — + г Т 2Т .»-Т 6 Ьх' (8,32) Введенное ранее конечно-разностное представление, например в форме (8.4) и следующих за этим соотношений, может рассматриваться в рамках соотношения (8.29), если определить конечно-разностные массовые операторы по направлениям с помощью выражений М, = (М» ) =(О, 1, 0), так что конечно-разностные выражения будут содержать только три члена.