CC BY
42
_____________УЧЕНЫЕ ЗАПИСКИ КАЗАНСКОГО УНИВЕРСИТЕТА
Том 153, кн. 4 Физико-математические пауки
2011
УДК 519.6
ПРОСТЫЕ АЛГОРИТМЫ ВЫЧИСЛЕНИЯ КЛАССИЧЕСКИХ АПОСТЕРИОРНЫХ ОЦЕНОК ПОГРЕШНОСТИ ЧИСЛЕННЫХ РЕШЕНИЙ ЭЛЛИПТИЧЕСКИХ УРАВНЕНИЙ
В. Г. Корнеев
Аннотация
Как известно, «классический» подход к апостериорной оценке погрешности приближенных решений задач механики твердого тела вытекает из встречных вариационных принципов Лагранжа и Кастильяпо. Если, например, задача липейпая и ее приближенное решение удовлетворяет геометрическим условиям, то потенциальная энергия ошибки оценивается потенциальной энергией разности тензора напряжений приближенного решения и любого тензора, удовлетворяющего уравнениям равновесия. Мы показываем, что для многих задач вычисление уравновешенных тензоров требует асимптотически оптимального числа арифметических действий. Улучшены также известные апостериорные оценки посредством произвольных неуравновешенных тензоров напряжений. Численные эксперименты показывают, что паши оценщики погрешности обеспечивают весьма хорошие индексы эффективности, как правило, сходящиеся к единице, имеют лилейную сложность и робастны.
Ключевые слова: апостериорные оценки, погрешность приближенных решений, метод конечных элементов.
1. Введение
Гарантированная апостериорная оценка оценивает погрешность сверху, используя только данные задачи и полученное приближенное решение, то есть дает надежный ответ, удовлетворяет ли приближенное решение требованиям к точности. Процедуры апостериорной оценки погрешности, или. как их еще называют. апостериорные оценщики погрешности, становятся важным компонентом программ численного решения краевых задач, встречаемых в научных исследованиях и инженерных приложениях. Индикаторы, погрешности, дающие представление по крайней мере о порядке погрешности, имеются в ряде коммерческих пакетов программ, таких, например, как ANSYS, ABACUS. FLUENT и др. Отчасти это связано с тем. что современный заказчик численного исследования, как правило, ие является специалистом в области вычислений. Таковым нередко не является и непосредственный пользователь пакетов программ. Потому создание быстрых алгоритмов апостериорного оценивания погрешности одно из наиболее актуальных направлений современной вычислительной математики. Хотя основы построения апостериорных оценок погрешности известны давно и восходят к встречным, экстремальным принципам минимума потенциальной энергии Лагранжа и максимума дополнительной работы Кастильяпо. интенсивное развитие этого направления началось в конце прошлого века. Это произошло в связи с резко возросшей
ролыо численного анализа и математического моделирования в науке и инженерной практике и увеличением мощности компьютеров. Имеется несколько сложившихся техник вывода, позволивших получить весьма эффективные апостериорные оценки. Центральная среди них основывается на применении уравновешенных (сбалансированных), то есть удовлетворяющих уравнениям равновесия (баланса), полей напряжений (потоков). Ее примером является метод равновесных невязок, который успешно применяли Эйнсворт. Демкович и Ким [1]. Люс и Вольмут [2]. Вейчодски [3]. Брасс и Шёберль [4], в работах которых имеется обширная библиография. Однако, как и в других техниках, равновесность понималась в дискретном смысле, а применение полей, удовлетворяющих уравнениям равновесия в точности, избегалось под предлогом, например, вычислительных и других трудностей. В настоящей статье мы показываем, что прямое получение уравновешенных, то есть в нашем понимании удовлетворяющих уравнениям равновесия в точности, тензоров напряжений может быть дешевым по вычислительной стоимости и позволяет получить простые гарантированные апостериорные оценки погрешности.
Рассматриваются две техники получения таких оценок. В одной на основе значений решения методом конечных элементов (МКЭ) в точках суперсходимости определяется достаточно хорошее из возможных дифференцируемое приближение г тензора напряжений а задачи. Компоненты тензора г могут находиться, например. как кусочно-полиномиальные функции конечно-элементных пространств класса С. Тензор г корректируете затем до уравновешенного тензора т = г + + 5т, который и применяется в апостериорной оценке погрешности. При первом краевом условии (на границе тела заданы перемещения) вычисление корректирующего тензора 5т сводится к вычислению одномерных интегралов от невязки в уравнениях равновесия. Поскольку в изложенном алгоритме не требуется решать какие-либо вспомогательные задачи, в частности системы алгебраических уравнений (СЛАУ), мы называем его алгоритмом прямого вычисления апостериор-
т
иутем численного решения методом Галеркина двойственной задачи, выражающей принцип Кастильяно максимума дополнительной работы. В качестве координатных применяются специальные локализованные самоуравновешенные тензоры,, удовлетворяющие уравнениям равновесия при отсутствии нагрузки. Они позволяют использовать МКЭ для решения двойственных задач. В случае двумерных эллиптических уравнений второго порядка для решения двойственной задачи оказывается возможным применять тот же МКЭ, которым решалась основная задача. Таким образом, для получаемых вспомогательных СЛАУ применим богатый набор алгоритмов, разработанных для МКЭ.
Наш подход позволяет также получить новые общие гарантированные вычисляемые апостериорные оценки погрешности, не содержащие констант, отличных от единицы. Примером являются оценки (2.5), (2.6) теоремы 1 и оценка (2.12). Их преимущества в том, что они наиболее близки наиболее точным, но практически не вычисляемым известным оценкам н в то же время алгоритмы их вычисления просты и во многих случаях могут быть оптимальными по вычислительной работе.
Мы приводим результаты численных экспериментов, подтверждающих эффективность алгоритмов, основанных на излагаемом подходе. Хотя краевые задачи в этих экспериментах относительно просты, они служили для тестирования алгоритмов апостериорных оценок погрешности и в других работах (см. [5 9]). Дополнительные результаты численного тестирования наших алгоритмов апостериорных оценок погрешности можно найти в [6, 7]. Отметим, что ради экономии места мы везде предполагаем все функции достаточно гладкими, если требования к гладкости опущены.
2. Апостериорные оценки погрешности для приближенных решений задач линейной теории упругости
Пусть Е - модуль Юнга, V - коэффициент Пуассона, удовлетворяющие неравенствам
0 < С1 < Е(х) < С2, 0 < < v(x) < С2, и < 0.5
\Т
С1,
С2,
{£м}
м(х)
3
(и1(х),и2(х),из(х)) Т - вектор переме-симметричные тензоры напряжений и
С1 С2
щений, а = {ам}к, г=1 “ - I~*,чк,1=1
деформаций, / и 4 - векторы объемных и поверхностных нагрузок, действующих на трехмерное тело П, х = (х1? х2, х3) . Верхний индекс Т означает транспонирование.
Краевая задача линейной теории упругости для изотропного материала в области П со смешанными краевыми условиями на границе д П = Ги У Гст, Ги р| Гст = 0, как известно, описывается системой уравнений
ё1у а = /, ап
£(м) = {ем(м)} к;=1 Е
1 /диь дил
= г (з^7 +
дм;
(2.1)
(2.2)
ак,к —
(1 - ^(1 - 2v)
(1 — Ь')£и ,к+^{£к+1 ,й+1+£к+2,к+2) , Ск,пг — ;---£к,пи (2.3)
] 1 + V
Е
где к = ш, ап - напряжение в точке границы на площадке с внешней нормалью дП
П и части границы Ги и Гст везде считаем односвязными.
Пусть Н1(П) = [И ^П)] , И к (П) - пространство С.Л. Соболева функций, имеющих суммируемые с квадратом производные до порядка к включительно, Уф =
{V е Н 1(П) :
= ^ } > <3/, г — множество тепзоров а, удовлетворяющих
уравнениям равновесия (2.1) в обобщенном смысле, _ энергетические
нормы, выраженные через перемещения и напряжения, соответственно. Нормы функции V в пространствах Ь2(П), И к(П) и квазинорму порядка к обозначим через |Н|о,п, ||V^к,п и ^|к,п соответственно.
хк
границу области не более чем дважды, сначала часть границы Гк, а , а затем часть границы Гк, ь> и пусть хк = а(хк+1, хк+2) - уравнение поверхности Гк,а.
Обратимся для простоты к случаю первой краевой задачи, когда Г и = д П; используя в этом случае вместо 3/, г обозначение 3 /. Конструктивно множество тензоров т е 3/ можно определить следующим образом. Задаем произвольные достаточно гладкие касательные напряжения тк,; = т;,к, к = I, и функции $к,к(хк+1, хк+2) = ^к,к (вк(хк+1, хк+2), хк+1, хк+2 ) значений нормальных напряжений на поверхностях Гк,а. Напряжения тк,к вычисляем по формулам
тк,к ^к,к (хк+1? хк+2 ) +
/к -
дтк,
к+1
дтк,
к+2
дх
к+1
дх
к+2
(Пк,хк+1,хк+^ ^Пк. (2.4)
аь(жь+1,жь+2)
Имеется ряд других способов определения тензоров т е 3 (см. [6, 7, 10]). Можно, например, задать сначала произвольный тензор я и по нему определить т е 3 /.
Теорема 1. Пусть V - произвольный вектор) из Уо, <г(^) - тензор напряжений, получаемый по и = V согласно (2.2), (2.3), г - произвольный симметричный тензор с компонентами из Нх(0). Тогда справедливы оценки
[и - V]и < [о'Ы - т]ст, [и - V] и < |_^) - г]ст + |/т_и, (2-5)
[и - V] и < |_^) - г]ст +
Хк
^Г- \ 1 (; 0гк,к+1 огкк+2\, ,
+ ]> —(Ік-—,----------------^----------^----- )(т,хк+1,хк+2)йщ , (2.6)
УЕ ' ®Хк &хк-\-1 иХк+2 ' О,
Й-1,2>3 ак
где а& = ай(хй+1,Хй+2) ? т - тензор из е компонентами ^ для к = I,
т&,& г&,&(а&? х&+1 ,х^+2) + хк
+ / (*- (27) а к ( х к +1 ?х к + 2 }
£т = diag [ ^т1;1, £т2,2, ^тз,з ] - диагональный тензор с компонентами
Хк
л- (і 0гк,к+1 0^,й+2\,
5тк,к= 1к-----^-------^------------^----)(щ,хк+1,хк+2)ащ. (2.8)
.) ^ джЙ ожЙ+1 ожЙ+2 /
ак(хк+1,хь+2)
Доказательство. Нетрудно проверить, что определения тензора т посредством (2.4) и (2.7) эквивалентны и т € ^/. Поэтому первая оценка (2.5) есть известная, называемая в литературе классической, оценка, которую можно найти у С.Г. Михлина в [11]. В настоящее время ее доказательство требует нескольких строк. Действительно, из (2.1)-(2.3) заключаем, что для любых тензоров а € и сто € ^о,о имеем
J [а : е(ад) + f • ад] йж + J £ • ад йв = 0, J ао : е(ад) йж = 0 V ад € V о. (2.9)
п гст п
Первое тождество получаем умножением уравнений равновесия на V и интегрированием по частям с использованием остальных соотношений (2.1) (2.3). Оно является обобщенной формулировкой уравнений равновесия. Второе тождество частный случай первого и выражает известный в механике принцип: работа виртуального изменения напряженного состояния на любом виртуальном изменении деформированного состояния равна нулю. Так как (м — V) € Vо и (а(м) — т) € € ^о,о, то второе тождество (2.9) приводит к равенству
1_ан -т ] I = |_аи -а(и)] ^ + Им)—т ] ^ = 1_м—vJ ^ + Им)—т ] и
из которого вытекает первая оценка (2.5). Другие оценки требуют лишь применения неравенства треугольника, неравенства Коши и элементарных оценок коэффициентов матрицы упругих коэффициентов. □
Чтобы сравнить с известными апостериорными оценками, рассмотрим для простоты задачу Дирихле для уравнения Пуассона: Дм = f(ж), ж = (жх,ж2) € О, =0
дП
|У(г,- г0||о,„ < (1+£)||^г1- 2/||^Г2+ (і + | ) ||/ ~ V • у|| " ± (2.10)
||Щи - г.)|| I п < (1 + е) || V, - у\\ + СП (1 + I ) ||/ - V • у\\I п, (2-11)
в которых с а - постоянная го неравенства Фрндрнхса, е - произвольное положительное число, ЦадЦ -1 п - норма фупкции ю в пространстве Н-1(О). За историей вопроса и оценками (2Л0), (2.11) мы отсылаем к [12 14, 15]. Первая не применяется на практике, поскольку вычисление нормы ^•|_1 а затруднительно. Вторая оценка является вычисляемой, но она содержит постоянную, зависящую от области. Так же, как (2.6) (см. [6, 7]), мы получаем оценку
11У(м-”)|| О,а ^ (1+е) ||У-у\\0,а '
хк
I 2
+ (1 + ”)|| 53 ак(/~ у 'У)(т,хъ-к)<1щ (2.12)
к=1,2 , ч - ’
ак (х3 — к)
где - произвольные функции, удовлетворяющие равенству «1 + «2 = 1 • Она вычисляема, причем достаточно просто, поскольку требует лишь вычисления двух одномерных н двух двойных интегралов. В то же время, как и (2.10), она не содержит констант, отличных от 1. В сравнении с (2.11) оценка (2.12) представляется
са
невязку, а интегралы от нее.
3. Получение сбалансированных потоков путем решения двойственных задач
Численное решение двойственной задачи может быть не более сложным, чем решение исходной задачи. Рассмотрим уравнение Пуассона в двумерной области
—Дм = /(х), х = (ж1, Ж2) € О,
дм
^ ? "о
г в дп
(3.1)
Г N
где дО, = ГдУГдг, Гдр|Гдг = 0. Для этой задачи уравнения баланса имеют вид
ди_ сН2_
Я ' Я */ ’ 7
дх1 дх2
Г N
и следовательно, множество сбалансированных потоков £ имеет вид Q /,д = | £ : —J£ •У^х = ^(^) VV € Vо(О) |,
где
^(V) = J /V ^х + J gv ^в, Vф(О) = | V € Н 1(О) : V = ф |.
а ГN
Обозначение Q/ сохраняем для пространства, соответствующего случаю Гд = дО. Пусть также в - координата па границе дО и во - ее наименьшее значение на Г^. Поскольку уравнение (3.1) можно свести к однородному уравнению, сформулируем результат для этого случая. Для простоты считаем также множества Гд, Г N односвязными, а фу нкцию ф равной нулю в точках их раздела.
Теорема 2. Пусть данные задачи (3.1) достаточно гладкие и / = 0, так что функции
ф*(в)= ! дсЬ, веГдг, д* = ^(«),
в € Г
определены. Пусть также ю - решение задачи —Дю = 0, х = (ж1, Ж2) € О, ю
Тогда
а* ди’
= Ф , тг-
г N дп
= 0*. (3.2)
Г в
ди ди> ди ди> (3 3)
дх1 дх2 дх2 дх1
Доказательство. Решение задачи (3.1) минимизирует функционал
1 [
J(v) = — а(г>, г’) — -Р(г’) Уг’£Уф, где а(г’,ги) = Х71’-Х7гис1,х, (3-4)
а
в то время как решение я двойственной задачи минимизирует функционал J*(£):
7*(г) = шш *(£), 7*(£) = I £ • £йх + I (£ • п) фйв. (3.5)
/,а 7 .1
а Гв
Взяв любой фиксированный вектор £ /,д € Q/,д, приходим к эквивалентной формулировке относительно я = я о,о + £/,д в виде тождества
^( яо,о + / ) • £ о, о йх +J ф (£о,о • п) йв = 0 V £о,о € Qо,о. (3.6)
а Гв
Нетрудно убедиться в том, что линейное пространство Q 0,0 можно задать посредством производящих функций как
<Зо,о = | ^ : +к = (-1) Й_1 —- Vг<; £ Т-Уо.г^ }; (3-7)
= {ю € Н 1(О) : ю = х },
Г N
и переформулировать (3.6) в виде: найти функцию ю € Ж^*,^ , удовлетворяющую тождеству
J (Ую + £ /,о) •У^йх + J ф (У-0 • п) йв = 0 V ^ € Ж0, ГN. (3.8)
а Гв
Проинтегрировав второе слагаемое по частям, придем, как можно убедиться, к тождеству
J Ую •У^йх = ^*(^) ^ У 0*^йв V^ , (3.9)
а Гв
которое при / = 0 и есть слабая формулировка задачи (3.2). □
Заметим, что конструктивное определение вектора £/,д достаточно простое.
Можно, например, взять произвольный вектор £/ = (£1,/ , £2,/)^ € Q/ с компонентами
х к
= 0к^(ак,х^-к) + ^ J/(??*;, жз-й) <*?*;, (3.10)
ак
вычислить поток £„,/ = / • п в направлении внешней нормали в точках границы Г« и функции
в
Ф/ = J *П,/ Ф* = Ф* - Ф/. (3.11)
в 0
После этого можно принять
, _ , , , , _ / дф* дф*\
+ Ъ.Г,9- [дхз, дх1),
где Ф* - функция, заданная на П и совпадающая с Ф* из (3.11) на Г«.
Следствие 1. Равновесные поля для апостериорных оценок погрешности
можно получать посредством МКЭ. При этом для решения основной и двой-
ственной задач .можно, вообще говоря, применять один и тот же МКЭ, а для решения СЛАУ МКЭ основной и двойственной задачи богатый набор быстрых алгоритмов, разработанных для МКЭ. В случае регулярных тлшсопряж-енных эллиптических уравнений различия вычислительных схем обусловлены лишь различием, краевых условий основной и двойственной краевых задач. Наличие вырождения коэффициентов и других нерегулярностей м.ожет потребовать различных способов их учета в схемах МКЭ для указанных краевых задач.
Пусть - пространство МКЭ с криволинейными вблизи границы конечными элементами, ассоциированными с треугольным базисным элементом порядка р, и удовлетворяются соответствующие условия квазиоднородности. Такие МКЭ-модели предложены и исследованы в [16, 17]. Достаточные условия квазиоднородности можно переформулировать следующим образом. Отображение х = X (£) : 6о ^ 6 базисного элемента 6о на конечный элемент ансамбля 6 представимо в виде суперпозиции X(£) = У(2(£)) двух отображений, где ? = 2(£) : 60 ^ 6^ -аффинное отображение и 6^ — треугольник с прямолинейными сторонами, верши-
6
отображение х = У (?) : 6^ ^ 6 имеет ограниченные производные порядка не выше р +1. При этом ансамбль треугольников 6^ удовлетворяет обычным условиям квазиоднородности, а ||У(?)||то < с, где || • ||то - норма в пространстве С.Л. Соболева
С+1ы.
Для простоты ограничимся формулировкой результата о сходимости приближенных решений двойственной задачи отдельно для случаев задач Дирихле и Неймана.
Пусть П1 = ид, и^ - решения МКЭ из прострапства ир(П) соответствующих двойственных задач (3.9), и = (21,^, Приближенные значения компонент
потока находятся по формуле
- = (-1) 6-1 +**,/,о- (3.12)
дх^ дхз_д.
Теорема 3. Пусть решения и = ид, задач Дирихле (Гд = дП) и Неймана
(Г^ = дП) и правая часть / удовлетворяют условиям
х к д {
/> ч|мдг | 1,П < С ( ||/|| р-1,П + |Ы| р-0.5,дп ), дхз_й ,]
а к
|ид | 1,п < С ( ||/1| р-1,П + ||Ф|| р+0.5,дП ), с = COnst.
Тогда
X k
II v« - zh II o,n < hp | ^2 —--- I f{in,x3-k)di]k
1 „ -~3-fc „ ,2 a k
Мы опускаем доказательство этой теоремы, которая есть прямое следствие теоремы 1 и известных результатов о сходимости МКЭ (см. [16. 17]).
Замечание 1. Вектор удовлетворяет уравнению баланса в области. Однако краевому условию Неймана на rN он удовлетворяет только в том случае, если ф* G UP(^)|rN j то есть ф* принадлежит пространству следов на rN КЭ функций. Если это не имеет места, то фактически G Q f gh, где соответствует функции ф^, получаемой в результате аппроксимации в схеме МКЭ функции ф*. При этом можно воспользоваться апостериорной оценкой (14) из [10] с дополнительным слагаемым вида c ||g — || o,rN > c = const, в ее правой части. Дополнительное
слагаемое имеет, как следует из оценки (14) работы [10], более высокий порядок малости по сравнению с основным членом.
4. Нелинейные задачи
Уравнения равновесия не зависят от типа закона Гука и сохраняют свой вид, например, для ортотропных, трансверсально изотропных и других типов упругих тел. Они также не изменяются для широкого круга физически и геометрически нелинейных тел. Следовательно, способы получения уравновешенных и самоурав-новешенных тензоров напряжений, рассмотренные в настоящей статье, применимы для широкого круга задач механики сплошной среды. При этом прямые способы нахождения указанных тензоров одинаковы для линейных и нелинейных задач. Перечисленные выше факторы влияют только на относительно дешевые операции вычисления тензора напряжений для приближенного решения задачи и энергетической нормы нлн функционалов потенциальной энергии в правой части апостериорной оценки.
Обратимся к общему случаю нелинейных задач механики твердого тела, для которых справедливы принцип Лагранжа минимума потенциальной энергии и принцип Кастпльяно максимума дополнительной работы, что позволяет получить оценки сверху и снизу действительной потенциальной энергии тела. Как хорошо известно, в этом случае апостериорная оценка погрешности имеет вид
L(v) — L(u) < L(v) — С(т), (4.1)
где L(v) - полная потенциальная тела на перемещениях v, удовлетворяющих всем геометрическим условиям, и - точное решение задачи, минимизирующее функционал L, т - любой тензор напряжений, удовлетворяющий уравнениям равновесия, C - функционал дополнительной работы на тензоре т. Пусть V - пространство перемещений с конечной энергией. Если || • ||у — норма в V, для которой выполняется неравенство
в ||v — u ||v < L(v) — L(u), в = const, (4.2)
то из (4.1) следует апостериорная оценка погрешности
в ||v — u ||у < L(v) — C(т). (4.3)
Оценка (4.1) выражает основополагающие свойства встречных вариационных
принципов Лагранжа и Кастильяно (см., например, [18 22]). Оценка (4.2) может
иметь место, например, для задач с монотонными дифференциальными операторами в частных производных. Математические аспекты оценок (4.1), (4.2) можно
найти в монографиях [23 25] и др. Они относятся к теории двойственности и теории монотонных/коэрцитивных дифференциальных операторов, позволяющей, в частности, сформулировать условия гладкости к исходным данным и условия на тип нелинейности, при которых справедлива оценка (4.2).
т
ближенных решений (например. МКЭ) нелинейных задач те же самые, что и для соответствующих линейных задач. Таким образом, при вычислении апостериорной оценки отличия связаны только с несколько более сложными по сравнению с линейным случаем вычислениями функционалов L(v) и С(т) для заданных v т
мер. апостериорный оценщик, оптимальный по вычислительной работе для линейной задачи теории упругости, остается оптимальным по порядку и для физически нелинейной задачи теории упругости.
Остановимся кратко также па вычислении т и одновременно функционала С(т) путем решения двойственной задачи, эквивалентной задаче минимизации функционала дополнительной работы. Основная задача заключается в нахождении u G U такого, что
L(u) = inf L(v), U С V, (4.4)
LV
хово пространство с нормой || • ||V, U - замкнутое, выпуклое подмножество V.
Тензор т и функционал дополнительной работы С(т) являются соответственно двойственной переменной н двойственным функционалом по отношению к основной задаче (4.4). При этом тензор т принадлежит множеству Qf,t тензоров а, удовлетворяющих уравнениям равновесия (2.1) в обобщенном смысле, и представим в виде т = тд( + т0,0, т0,0 G Qo,o > где Q0,0 - пространство самоуравновешенных тензоров, то есть удовлетворяющих уравнениям равновесия при f = 0, t = 0. Двойственная задача заключается в нахождении тензора а, для которого
С(а) = iQf С(т)= iQf С(т^ + 7). (4.5)
т eQf.t TEQo.o
Если С (тд( + y) - выпуклый, полунепрерывный снизу функционал от 7, коэрцитивный на рефлексивном банаховом пространстве Q0j0, то (см., например, [25])
С(т) ^ С(а) = L(u) < L(v) Ут G Q/>t, V v G U,
откуда следует (4.3).
Ряд авторов считает, что использование апостериорных оценок (4.3) невыгодно по двум причинам: 1) сложность прямого вычисления по v тензора G Q/,t, аппроксимирующего а, и 2) сложность, особенно в случае нелинейной задачи, решения двойственной задачи (4.5). Первая из причин уже обсуждалась выше, поэтому добавим лишь несколько слов по поводу второй. Как и в линейном случае, можно получить вычислительные алгоритмы со свойствами, подобными МКЭ. При этом вопреки 2) численное решение двойственной задачи может оказаться даже более простым по сравнению с основной задачей. Примером может служить задача с коэффициентом Пуассона, близким в некоторых подобластях к 0.5.
5. Численные эксперименты
Результаты ряда численных экспериментов, включая примеры 5.1, 5.3, приведены в работах Ануфриева, Корнеева и Костылева [6, 7], где можно найти их более полное обсуждение. Позднее Костылевым были выполнены расчеты для других задач, в частности для второго из обсуждаемых ниже примеров.
5.1. Эллиптическое уравнение с разрывным коэффициентом. Рассматривалось уравнение V • pVu = f в квадрате (0,1) х (0,1) с р = р1 = 10-2 и р = р2 = 102 в его левой и правой половинах соответственно. Правая часть и краевые условия соответствовали точному решению
Зтг (‘г’2 + 1; Х <
w = (cos(27rx) — l) cos —— у I
2 1.25 + 0.25^-(*-i)2^, X>1
р2 р2 2
При этом иа левой и иижией гранях ставилось краевое условие Неймана, а на
(0, 1) х (0, 1)
вался равномерной квадратной сеткой с четным числом шагов, которое изменялось от 8 до 1024. Пространством МКЭ служило пространство непрерывных кусочнополилинейных функций. Для рассматриваемой задачи апостериорная оценка погрешности имеет внд
||V(u — ufem)|| р < ||р Vufem — t|| р_, (5.1)
где для z = (zi, ^2)^ и ш > 0 имеем
I и 2 _
У w(z2 + z|) dx.
Для решений МКЭ, полученных на разных сетках, вычислялись
е = || V (и — и(ет )| И П = ||р — ^| р _,
где е - норма точной погрешности решения МКЭ и п _ апостериорная оценка погрешности, сравнение которых позволяет судить об эффективности апостериорного оценщика. Ниже приводим весьма простой алгоритм вычисления апостериорной оценки погрешности, более подробное описание которого можно найти в [6, 7] (алгоритм 2.4).
Алгоритмы прямого вычисления апостериорной оценки погрешности
Шаг 1. Для каждого узла ж(г) = (*1 /г, г) € П, * = (*ь *2), вычисляется значение сеточной функции = {"У2д(*)} о 1 п, которое является конечно-
разностной аппроксимацией в узле производной д(рдм/дх1)/дж1, вычисленной по конечно-элементному решению и (ж). Для внутренних узлов на горизонтальных сеточных линиях применяются формулы
«2,1(*) = г-2 [«1д(*1 + 1, *2) — «1д(*)],
«1,1(*) = р(г(*1 — 0.5,*2)) [^т(г*) — (г(*1 — 1,^))]•
Для узлов * = (0, *2) та ос и ж2 полагаем
«2,1(0, *2) = Р(0, *2)д2и;п1;/дж!(0, *2),
где ^ тг (ж 1) - интерполяционный полином Лагранжа 3-го порядка по значениям и(ет в узлах ж(г), *1 = 0,1, 2, 3.
Аналогично, для узлов * = (п, *2)
«2,1(п, *2) = р(п, *2)д2^г/дж1(п, *2),
Табл. 1
N е V 1сК
64 9.14308 1.02954 • 101 1.12604
256 4.44483 4.77467 1.07421
1024 2.20395 2.25318 1.02234
4096 1.09958 1.10628 1.00609
16384 5.49486 • 10 _1 5.50364 • 10 _1 1.00160
65536 2.74705 • 10 _1 2.74818 • 10 _1 1.00041
262144 1.37348 • 10 _1 1.37362 • 10 _1 1.00010
1048576 6.86734 • 10 ~2 6.86751 • 10 ~2 1.00003
где и'шь (х 1) - интерполяционный полином Лагранжа 3-го порядка по значениям В узлах Ж(г) , *1 = п, п — 1, п — 2, п — 3.
Найденная сеточная функция однозначно определяет кусочно-билинейную функцию /т1;(ж) € V(О) равенствами *(х(г)) = «2д(*), ж(г) € О.
Шаг 2. Находим поток £ = (^ 1?^2)Т, удовлетворяющий уравнению баланса, путем вычисления интегралов
®1 Х2
^(х) = ^(^1,Х2) йх*1, ^2(х)^У [ /(Х1,^2)+ ^(Х1,^2)] ^2-
0 0
Шаг 3. Вычисляем правую часть ||рУи£ет — ^|р-1 оценки (5.1).
Замечание 2. Имеется много других вариантов алгоритмов, основанных на построении по решению МКЭ потока, в точности удовлетворяющего уравнению баланса. Например, можно интерполировать по средним значениям в узлах непосредственно одну из компонент ^, а другую компоненту находить интегрированием / — д^/джд;. Несколько из них можно найти в [6, 7, 10].
Основным показателем эффективности апостериорного оценщика в отношении точности является индекс эффективности 1е« = п/е, где левая и правая
части в (5.1), то есть точное значение энергетической нормы погрешности решения МКЭ н ее апостериорная оценка. Результаты вычислений показаны в табл. 1, в которой N - число неизвестных. Из нее видно, что апостериорный оценщик весьма эффективен и индекс эффективности быстро сходится к 1. Он также является быстродействующим, что видно из рис. 1, где показана зависимость от N затрат времени РС в миллисекундах. Верхняя линия отвечает времени решения СЛАУ МКЭ многосеточным методом, нижняя времени вычисления апостериорной оценки, которое оказывается пропорциональным числу неизвестных. Таким образом, апостериорный оценщик оптимален по вычислительной работе (имеет линейную вычислительную сложность) н менее трудоемок, чем многосеточный метод решения СЛАУ МКЭ. Он, кроме того, является робастным, так как скачок коэффициента р, иропорциональный 104, не привел к ухудшению качества алгоритма: индекс эффективности апостериорной оценки, как и в других примерах, быстро сходится к единице, то есть апостериорная оценка весьма близко следует за точным значением погрешности.
Проводились также численные эксперименты для задачи плоской теории упругости. Рассматривалась первая краевая задача и = (и1;и2)т = плоской деформации в единичном квадрате О = (0,1) х (0,1), которая решалась методом
Рис. 1
Табл. 2
N е ?7(1) /(1) ?/2) 1—1
16 3.29126 10-1 2.11430 6.42399 2.45022 7.44462
64 1.63481 ю-1 5.54014 • ю-1 3.38885 3.54162 • 10-1 2.16638
256 8.17444 10~2 1.94481 • 10-1 2.73914 1.46349 • 10-1 1.79033
1024 4.08766 10~2 8.62833 • 10~2 2.11083 5.45946 • 10~2 1.33560
4096 2.04389 10~2 4.18226 • 10~2 2.04622 2.56000 • 10~2 1.10387
16384 1.02196 10~2 2.07433 • 10~2 2.02976 1.05054 • 10~2 1.02797
65536 5.10979 10~3 1.03501 • 10~2 2.02554 5.14644 • 10~3 1.00717
262144 2.55490 10~3 5.17229 • 10~3 2.02446 2.55953 • 10~3 1.00181
1048576 1.27745 10~3 2.58580 • 10~3 2.02419 1.27804 • 10~3 1.00046
конечных элементов на равномерной квадратной сетке с п2 ячейками. Принималось, что Е =100, V = 0.2. Правая часть и функция соответствовали точному
решению
и 1_(ж) = 8ш(пЖ1 ) 81п(2пЖ2) + Х1 + Х2,
и2(х) = 81п(2пхх) вт(пх2) + 0.25(х1 + 1)(х2 + 1).
Апостериорные оценки вычислялись по алгоритмам 3.1, 3.2, описанным в [6, 7]. В первом алгоритме вычислялась и интерполировалась первая производная дт1,2/дх2 напряжения т^, полученного методом конечных элементов. По ней посредством операций вычисления одномерных интегралов и дифференцирования вычислялся тензор напряжений, удовлетворяющий уравнениям равновесия. Во втором алгоритме с целью получения симметричного относительно осей хк алгоритма вычислялась и интерполировалась вторая производная д 2т12/дх1 дх2 напряжения Т1,2 , полученного МКЭ. Соответственно, равновесные тензоры напряжений находились путем вычисления двойных интегралов. Численные результаты приведены в табл. 2, в которой N - число узлов сетки, е - точная погрешность решений МКЭ в энергетической норме, г/(к\ к = 1, 2, - апостериорные оценки погрешности решений МКЭ в энергетической норме, /^ - соответствующие индексы эффективности.
Рис. 2
Табл. 3
N Ьп-С ггог е 1аа
25 0.00149957 0.06471507 5.75502518
100 0.01604613 0.14359042 1.5897243
400 0.00088007 0.03620792 3.59732922
1600 0.00021309 0.01927863 4.28158306
6400 6.4957 • 10-6 0.00956583 3.37412084
25600 1.7082 • 10-6 0.00467778 2.11054348
102400 4.3256 • 10 ~6 0.00232179 1.39053788
409600 1.0849 • 10 ~6 0.00115862 1.11269864
1638400 2.7144 • 10-7 0.00057902 1.02947749
На основании табл. 2 можно сделать следующие выводы. Апостериорные оценки погрешности чувствительны к способу их вычисления. Второй алгоритм оказался заметно более эффективным. Сравнивая с табл. 1. заключаем, что разрывный коэффициент в уравнении ухудшает погрешность значительнее (примерно на порядок). чем усложнение уравнений при гладком решении. Численные эксперименты подтверждают оптимальность предлагаемых апостериорных оценщиков погрешности плоской задачи теории упругости по вычислительной трудоемкости относительно числа неизвестных (см. дополнительную информацию в [6. 7]).
5.2. Уравнение Пуассона с решением, имеющим большие градиенты. В этом пункте изучалась задача Дирихле для уравнения Пуассона —Ап = f с точным решением вида
п = ж(1 — ж)у(1 — у) ехр(—1000[(ж — 0.5)2 + (у — 0.117)2]).
рассматривавшаяся также в [5. 8. 9]. На всей границе ставилось однородное краевое условие Дирихле. Точное решение показано на рис. 2 и имеет значительные градиенты. Применялся тот же МКЭ на равномерной квадратной сетке и тот же апостериорный оценщик, что и в первом примере.
Численные результаты сведены в табл. 3. В ее втором столбце приведена норма точной ошибки МКЭ в пространстве Ь2(0.), а в третьем - норма е (при р = 1).
Индекс эффективности снова сходится к единице, но значительно медленнее, чем в первом примере. В основном это связано с тем. что для решения задачи со значительными градиентами ради упрощения вычислений применялась равномерная сетка. Отметим, что быстродействие апостериорного оценщика такое же по порядку. как в предыдущем примере, так как оно слабо зависит от правой части, то есть апостериорный оценщик имеет линейную сложность. В то же время при сравнении с другого типа апостериорными оценщиками работ [5. 8. 9] результатов для одинаковых точных значений норм погрешности он обеспечивает более точные апостериорные оценки погрешности н. как следствие, лучшие индексы эффективности.
5.3. Алгоритм, основанный на решении двойственной задачи. Рассматривалась задача Дирихле для уравнения Пуассона в единичном квадрате О = = (0,1) х (0,1)
—△и = /(х), х Є О,
д П
Ч>-
(5.2)
для решения которой применялся МКЭ па равномерной квадратной сетке хк = = хкк^ = гк Н, Н = 1/п, %к = 0,1,..., п, с билинейными конечными элементами. Пусть V(О) — пространство непрерывных на О и билинейных на каждой ячейке сетки функций и {ф(г) } 1 ~ базис в нем, каждая функция ф(г) которо-
го равна 1 в узле х(г) = гН, г = (г1; г2), и 0 - в остальных узлах. При ^ = 0. что для простоты предполагается ниже, решение МКЭ ищется в подпространстве ^(О) С V(О) функций, обращающихся в 0 на д О. В соответствии со схемой, изложенной в разд. 2, для решения двойственной задачи можно использовать пространство V(О). В этом случае в качестве сеточного подпространства пространства Q0 самосбалансированных потоков имеем
Qh,0 = |Ъ : t = (дv/дх2, —дv/дxl )Т VV е V(О) |,
базисом в котором является |Ъ(г) = (дф(г)/дх2, — дф(г)/дх1 )Т | . Как следует
I ) г\ ,г2 = 1
из разд. 1 и 2, сбалансированный поток / = (Ъ1,/ , Ъ2,/) можно определить многими способами. В численных экспериментах применялись два наиболее простых способа, когда вычисления проводятся по формулам
хк
Ч/ =0, Ъ2,/ = J I(хь&) ^2, (5.3)
0
а также по инвариантным относительно наименования осей формулам
Хк
Ък,/ =0.^ I(^к,хэ-к) ^к. (5.4)
Таким образом, приближенное решение +Ъ/, е Q0, находилось путем
решения (п + 1)2 уравнений
ік,! = 0.5 J(^о + (їх = 0,
*1, *2
0,1,
(5.5)
п.
Табл. 4
N е 1аа -йг
16 1.59912 2.93605 2.05747 1.20610
64 7.60498 • 10-1 1.95927 2.16500 1.23228
256 3.70781 • 10 _1 1.26742 2.21586 1.24493
1024 1.84023 • 10 _1 1.07113 2.23084 1.24869
4096 9.18344 • 10 ~2 1.01826 2.23475 1.24967
16384 4.58949 • 10 ~2 1.00460 2.23574 1.24992
65536 2.29446 • 10 ~2 1.00115 2.23599 1.24998
262144 1.14720 • 10 ~2 1.00029 2.23605 1.24999
1048576 5.73594 • 10 ~3 1.00007 2.23606 1.25000
По теореме 2 система уравнений (5.5) есть система уравнений МКЭ для задачи Неймана
= 0. (5.6)
8П
При этом если г есть ее решение МКЭ, то так же, как и в (3.12), имеем
= (до/дх2 + tlJ, -дш/дх1 + ^2,/ )Т.
Ниже для наглядности приводим результаты вычислений для простой задачи с правой частью 5п2 вт(2пх1) вш(пх2) и точным решением и = вт(2пх1) вш(пх2). Две апостериорные оценки, соответствующие вычислению вектора I/ согласно (5.3) и (5.4), имеют два индекса эффективности, обозначаемые 1^- и IДля сравнения применялся также алгоритм прямого (то есть без решения вспомогательных задач) вычисления апостериорной оценки погрешности МКЭ, описанный в п. 4.1. Индекс эффективности этого алгоритма обозначаем 1ед. Значения нормы е = 1е|1,п погрешности МКЭ и перечисленных индексов эффективности в зависимости от числа неизвестных N = (п + 1)2 приведены в табл. 4.
Численные результаты позволяют сделать следующие выводы.
1. Сбалансированные потоки, найденные путем решения методом конечных элементов двойственной задачи, позволяют получать апостериорные оценки погрешности с весьма хорошими индексами эффективности, не превосходящими 1.3.
2. Сравнение 1е^ и !2№ показывает, что выбор вектора / влияет па точность апостериорной оценки. При этом вычисление / по инвариантным относительно наименования осей формулам заметно увеличивает точность.
3. Индексы 1^ те сходятся при Н ^ 0 к единице, в то время как 1ей- сходится. Отсутствие сходимости означает, что отношение погрешностей решения исходной задачи и двойственной при Н ^ 0 не уменьшается. Это объясняется тем, что для решения основной н двойственной задач применялись МКЭ одного порядка точности, точнее, один и тот же МКЭ. Чтобы обеспечить сходимость индекса эффективности Iф-, для решения двойственной задачи, видимо, нужно применять более точный по порядку МКЭ, чем для решения основной задачи. В то же время более простой алгоритм, не требующий решения вспомогательных задач, но основанный на использовании точек суперсходимости решений МКЭ, обеспечивает сходимость 1е® ^ 1 при Н ^ 0 и, следовательно, большую точность апостериорных оценок
Н
ложеннй. Он свидетельствует о том, что усложнение алгоритмов апостериорных
оценок погрешности приближенных решений, как правило, мотивированное стремлением повышения их точности, может привести к неоправданному усложнению алгоритмов.
Работа выполнена при финансовой поддержке РФФИ (проект X- 01-11-00667).
Summary
V.G. Korneev. Simple Algorithms for Calculation of the Classical A Posteriori Error Estimates of Numerical Solutions of Elliptic Equations.
As is well-known, for the problems in solid mechanics, “classical” approach to a posteriori error estimation stems from the Lagrange and Cast.igliano variational principles. If the problem is linear and an approximate solution satisfies geometrical restrictions, then potential energy of the error is estimated by the potential energy of the difference of the stress tensor corresponding to the approximate solution and any stress tensor satisfying the equations of equilibrium. We show that in many cases, construction of equilibrated stress fields can be done for a number of arithmetic operations, which is asymptotically optimal. This approach allows us also to improve known a posteriori estimates by means of arbitrary lionequilibrat.ed tensors. Numerical experiments show that our a posteriori error estimators provide rather good efficiency indices, which often converge to unity, have linear complexity, and are robust.
Key words: a posteriori estimates, error in approximate solutions, finite element method.
Литература
1. Ainsworth М., Denikowicz L., Kim C.-W. Analysis of the equilibrated residual method
for a posteriori estimation on meshes with hanging nodes // Comp. Met.li. Appl. Math.
Engrg. 2007. V. 196, No 37 40. P. 3493 3507.
2. Luce R., Wohlmuth B. A local a posteriori error estimator based on equilibrated fluxes //
SIAM J. Num. Anal. 2004. V. 42, No 4. P. 1394 1414.
3. Vejchodsky T. Local a posteriori error estimator based on the liypercircle method //
Neittaaiimaki P., Rossi Т., Korot.ov S., Onat.e E., Periaux J., Knorzer D. (eds.)
European Congress on Computational Methods in Applied Sciences and Engineering ECCOMAS 2004. Yavaskyla, Finland, 2004. 16 p. URL: http://www.imamod.ru/
~serge/arc/conf/ECCOMAS_2004/ECCOMAS_V2/proceedings/pdf/769.pdf, свободный.
4. Braess D., Schoberl J. Equilibrated residual error estimator for MaxsweU’s equations // Mat.h. Comp. 2008. V. 77. P. 651 672.
5. Gockenbach M.J. Understanding and implementing the finite element method. SIAM,
2006. 363 p.
6. Anufriev I.E., Korneev V.G., Kostylev V.S. Exactly equilibrated fields, can they be efficiently used for a posteriori error estimation? // Учен. зап. Казап. уп-та. Сер. Физ.-матем. науки. 2006. Т. 148, кп. 4. С. 94 143.
7. Anufriev J., Korneev V., Kostylev V. A posteriori error estimation by means of the exactly
equilibrated fields. Ricam Report .V 2007-07. Linz, Austria: Johann Radon Institute for Computational and Applied Mathematics, 2007. 54 p.
8. Tomar S.K., Repin S.I. Efficient computable error bounds for discontinuous Galerkin
approximations of elliptic problems. RICAM-report. .V® 2007-39. Linz, Austria: Johann Radon Institute for Computational and applied Mathematics, 2007. 20 p.
9. Repin S.I., Tomar S. A posteriori error estimates for nonconforming approximation of elliptic problems based on Helmholtz type decomposition of the error. RICAM-report.
2007-41. Linz. Austria: Joliann Radon Institute for Computational and applied Mathematics, 2007. 16 p.
10. Корнеев В.Г. В развитие классического подхода к апостериорным оценкам погрешности приближенных решений краевых задач // Сеточные методы для краевых задач и приложения: Тр. 6-го Всерос. семинара. Казань: Изд-во Казан, ун-та, 2007.
С. 162 167.
11. Михлии С.Г. Вариационные методы в математической физике. М.: Наука, 1964. 512 с.
12. Ainsworth М., Oden J.T. A posteriori estimation in finite element analysis. N. Y.: John
Wiley & Sons, Inc., 2000. 240 p.
13. Babuska I., Strouboulis T. Finite element method and its reliability. N. Y.: Oxford Univ.
Press, 2001. 802 p.
14. Neittaanmaki P., Repin S.I. Reliable methods for computer simulation: Error control and
a posteriori estimates. N. Y.: Elsevier, 2004. 305 p.
15. Репин С.И., Фролов М. Об апостериорных оценках точности приближенных решений краевых задач для эллиптических уравнений // Жури, вычисл. матем. и матем. физ. 2002. Т. 42, Л» 12. С. 1774 1787.
16. Корнеев В.Г. О построении вариациошю-разпостпых схем высокого порядка точности // Вестп. Лепипгр. уп-та. 1970. Т. 25, Л'! 19. С. 28 40.
17. Корнеев В,Г, Схемы метода конечных элементов высоких порядков точности. Л.:
Изд-во Лепипгр. уп-та, 1977. 255 с.
18. de Vebeke F. Displacement and equilibrium models in the finite element method // Stress Analysis / Eds. O.C. Zienkiewicz, G.S. Holister. London-N. Y.: Wiley, 1965. P. 145 197.
19. Arthurs A.M. Complementary variational principles. Oxford: Clarendon Press, 1980. 154 p.
20. Мосолов П.П., Мясников В.П. Механика жесткопластических сред. М.: Наука,
1981. 207 с.
21. Washizu К. Variational methods in elasticity and plasticity. N. Y.: Pergamon Press,
1982. 630 p.
22. Бердичевский B.JI. Вариационные принципы механики сплошной среды. М.: Наука,
1983. 448 с.
23. Glowinski R. Numerical methods for nonlinear variational problems. N. Y.-Berlin-
Heidelberg-Tokio: Springer-Verlag, 1984. 493 p.
24. Duvaut G., Lions J.-L. Inequalities in mechanics and physics. Berlin-Hedelberg-N. Y.:
Springer-Verlag, 1976. 397 p.
25. Ekeland I., Temam R. Convex analysis and variational problems. Amsterdam: N. Y.:
North-Holland Pub. Co., 1976. 402 p.
Поступила в редакцию 24.10.11
Корнеев Вадим Глебович доктор физико-математических паук, профессор, профессор Санкт-Петербургского государственного политехнического университета, Санкт-Петербургского государственного университета Е-шаП: УаdimKoтееь вуаЬоо. сот
