Fletcher-1-rus (1185917), страница 59
Текст из файла (страница 59)
Что касается дискретизации многомерных уравнений, проводимой на основе конечно-элементного и конечно-разностного подходов, то метод конечных элементов моделирует пространственные производные по данным, распределенным в пределах некоторой площадки, тогда как метод конечных разностей моделирует те же производные по данным вдоль определенного пространственного направления. Вообще говоря, конечно- элементная формулировка более точна, но вместе с тем и ме- направлениях х и у, то после деления всех членов на ЛхЛу получим ГдТ1 М»ЭМи~д ~ =ахМ„ЭЬххТП»+аиМ»Э1.„„Т, », (8.29) ). дг )г» $8.3.
Схемы расщепления 338 нее экономична, так как связана с расчетом большего числа членов. Если вместо линейных лагранжевых элементов используются квадратичные лагранжевы элементы, то каждый из массовых и разностиых операторов будет содержать уже по пять членов, таи что тензорное произведение типа М,З 1.екТ), з бу- п41 1-1 1 1+1 Прогонка в направлении х (постоянное)г) (а) Конечно-злементная формулировка и+1 — — Рй+1 г г г)г формулировка Рис. 8.3. Активные узлы для схем расщепления по уравнениям (8.93), (824) или (8.39), (8.40), дет содержать 25 членов. Однако такое выражение согласуется, как правило, с представлением производной дзТ/дуз с точностью до членов порядка 0(буз).
8.3,1. Приемо( расщепления конечно-элементных уравнений Здесь будет изложено два способа расщепления конечно- элементных уравнений. Первый из этих способов приводит н алгоритму, аналогичному схеме НПН (8.14), (8.15). Второй алгоритм сводится к выражению решения через поправки ЬТА м как по схеме (8.23), (8.24). Так же как и при исследовании одномерного уравнения диффузии, в соотношение (8.29) будет вводиться конечно-разност- и )г ! 1-1 1 1+1 Прогонке в направлении х (постоянное а) [Ы Конечно-Пазностная — — — -Рй+! г х ° г-— — — уй г l l )г-1 1-1 1 1+1 Прогонка в направлении у (постоянное1) --- — Рй+! г г — е(( х г (г-1 1 1 1 1+1 Прогонка в направлении У (постоянноеу) Гл. 8. Многомерное урзвненне диффузии ное представление производной дТ(д1, а пространственные члены будут моделироваться средневзвешенными значениями для и-го и (и+ 1)-го слоев по времени.
Таким образом, соотношение (8.29) заменяется выражением дтп+' М„Э̄— '=(а„М„Эй„„+ а„М„ЭЬнн) [(1 — 8) Т;, в+ 6Тг,з]. (8.33) Уравнение (8.33) представляет собой конечно-элементный эквивалент уравнения (8.17). Расщепление типа НПН с конечными элементами осуществляется путем введения в левую часть уравнения (8.33) дополнительного члена Ма„а 1.„„Э 1.нн Ь6~Т"+ — (1 — 8) Т"1, так что это уравнение принимает вид (Мл — ахб а11.„) Э (Мн — аФ Мьнн) Т7Х = = [М, + а„(1 — [1) ЫЕ„] Э [Мн + ан (1 — р) ог Ьнн] Т~, ы (8.34) Решение уравнения (8.34) может быть весьма эффективно проведено в два этапа.
На первом из них соотношение (М, — ал6 Ш.л) Т~, л = [М„+ ан (! — [)) Магно] Ть е (8.35) позволяет получить трехдиагональную систему на каждой из сеточных линий, параллельных оси х. На втором этапе соотношение (Мн — анр Жанн) Т7,"»' = [М, + а, (1 — [)) Л(з.„] Тг, е (8.36) приводит к трехдиагональной системе на каждой из сеточных линий, параллельных оси у. При р = 0.5 единственным различием между конечно-элементной конструкцией типа НПН в форме (8.35), (8.36) и конечно-разностной конструкцией НПН вида (8.14), (8.15) будет появление в первой из них массовых операторов. Конечно-элементный алгоритм по схеме НПН является экономичным, имеет второй порядок точности по времени и пространству и обладает безусловной устойчивостью прп р ) 0.5+(6 — 0.25)/з, где 6 определяется по формуле (8.44), а вместо з можно использовать з„или з .
Расщепление может осуществляться и другим способом. При этом член Тг~е, находящийся в правой части уравнения (8.33), разлагается в ряд Тейлора в окрестности и-го временнбго слоя с отбрасыванием членов 0(ЛР). После перегруппи- 4 8.3. Схемы расщепления ровки результат может быть записан в виде 1Мл Э Мц — ]3 Лг (а,Мц Э Елл + ацМ, Э 7.цц)] ЛТт "а = = Л( (алМц Э ~~л + ацМл Э1 цц) Ту. а. (8.37) Если в левую часть уравнения (8.37) ввести дополнительный член р~ Лт алацЕ„Э Ецц ЛТт, а, то получится расщепляющееся уравнение тМл — р Л!алТ.лл) Э (Мц — 8 Мац1.цц) ЛТ~,+~ = = Ж(алМц Э).ли+ ацМ,Э Лиц) Т7, ы (8.38) причем реализация расщепления включает соотношение (̄— 8 Л!ал 1.„) ЛТь а = Ж (а,Мц Э Ьлл + ацМ, Э 7.цц) Тт", а (8.39) для всех сеточных линий, идущих в направлении х, а также соотношение (Мц — р Лгаце.цц) ЛТ!~а = ЛТ~, а (8.40) для всех сеточных линий, идущих в направлении у. Уравнения (8.39), (8.40) можно сравнить с уравнениями (8.23), (8.24).
Правая часть уравнения (8.39) требует большего времени для своего расчета, чем правая часть (8.23). Расчет правых частей осуществляется на временнбм слое и. Как можно видеть из рис. 8.3, схема метода конечных элементов охватывает девять активных узлов на п-м слое по времени по сравнению с пятью активными узлами для схемы конечно-разностного метода. Однако процесс расчета членов уравнения (8.39), а также реализация решения трехдиагональных систем уравнений для обеих формулировок различаются несущественно. Как свидетельствуют вычислительные эксперименты, метод конечных элементов работает, вообще говоря, медленнее метода конечных разностей, но обычно обеспечивает более точное решение. Вычислительный алгоритм (8.39), (8.40) согласуется с уравнением (8.1) при ошибке аппроксимации порядка 0(Л(т, Лх', Лу'), если р = 0.5.
Такие же уравнения по схеме с расщеплением будут получены, если вместо линейных лагранжевых элементов использовать квадратичные лагранжевы элементы (п. 5.3.2); при этом, однако, операторы М„Т.цц и т. д. будут содержать максимальное число членов, т. е. пять вместо трех [Р!е1сйег, 1984]. 22 К. Флетчер, т. ! ЗЗ8 Гл. 8. Многомерное уравнение днффуанн 8.8.2, Т%111Р: обобщенная схема для реализации расчетов с конечногми разностями или конечными элементами Различные схемы расщепления, рассмотренные в п.
8.2.2, 8.2.3 и 8.3.1, могут быть объединены в одну обобщенную трехслойную схему, охватывающую как конечно-разностную, так и конечно-элементную формулировки. Мы покажем это здесь для области, изображенной на рис. 8.!. Распределение температуры в области 0 < х < 1, 0 < у < < ! определяется уравнением (8.1) со специальными граничными условиями Т(0, д, С) = 20+ 80у, Т(1, у, 1) = 20+ 80 [у — е-'""" яп (0.5лу)), Т (х, О, 1) = 20, Т (х, 1, 1) = 20+ 80 [у — е олон*г яп(0 5лх)] (8.4!)~ Для данного примера существует следующее аналитическое ре- шение: Т=20+ 80 [у — е-ее'"' яп (0.5лх) яп (О.блу)).
(8.42т (!+У)ат",+,' — тат," а Мх Э Ма = а (М~ Э е + М, Э Е„,) [(1 — Р) Ть е + 8Т", е 1 (8.43) т. е. соотношение, эквивалентное (8.28). Однако здесь массовые операторы (8.30) по координатным направлениям будут представлены в обобщенном виде М,=Мед= [6, 1 — 26, 6).
(8.44) Ясно, что конечно-разностная схема (п. 8.2.2 и 8.2.3) соответствует значению 6 = О, тогда как конечно-элементная схема (п. 8.3.1) соответствует значению 6 = !!6. После приближенной факторизации уравнения (8.43) получается следующий двухэтапный алгоритм, заменяющий (8.39), (8.40): Обобщенная схема может быть введена за счет трехслойной конечно-разностной дискретизации уравнения (8.29) прн а=а,=а„, чтодает $8.3. Схемы расщеплении Первый этап ~М вЂ” ( ~н ) Т,„~АТ! а=КН5" (8.45) тде ахН5" =(аАГ(МнЭЬаа+ Ма Э 1.„н) ТК а+ гМаЭМдАТ!', а1/(! + т). (8.46) Второй этап !(̄— ( ~и ) 7.„„1 АТ~,+' = АТ;, и. (8.47) Для уравнений, формируемых в точках вблизи границ, нужмо задать граничные значения величин АТ' в (8.45) и АТ"е' в (8.47). Надлежащий вид этих граничных значений проще всего определить путем свертывания (8.45) и (8.47) в одну комбинированную схему, получаемую за счет исключения АТ'. Тогда становится ясно, что требуемые для решения (8.47) граничные значения АТ" могут быть получены непосредственно из формул (8.41), тогда как граничные значения для АТ" должны вычисляться также из (8.41) через посредство обращенного соотношения (8.47); например, в точке ! =ИХ (х= !) имеем АТмх,а=[Ма — (! „!)Т.не)АТмх~е.
(848) Если граничные условия (8.4!) не зависят от времени, то для АТ* и АТ"+' получаются нулевые граничные значения. Известные граничные значения зависимых переменных АТ" и АТ"+' .переносятся в правые части уравнений (8.45) и (8.47) соответственно до решения уравнений. Реализация вышеописанного алгоритма проводится с помощью программы Т%01Р (рис. 8.4).