Вычислительный эксперимент и верификация реализаций апостериорной оценки для пластин Рейсснера – Миндлина Текст научной статьи по специальности «Математика»

Математика

DOI: 10.18721/JPM.12111 УДК 519.63

ВЫЧИСЛИТЕЛЬНЫЙ ЭКСПЕРИМЕНТ И ВЕРИФИКАЦИЯ РЕАЛИЗАЦИЙ АПОСТЕРИОРНОЙ ОЦЕНКИ ДЛЯ ПЛАСТИН РЕЙССНЕРА - МИНДЛИНА К.В. Киселев, М.Е. Фролов, О.И. Чистякова

Санкт-Петербургский политехнический университет Петра Великого, Санкт-Петербург, Российская Федерация

В работе рассматривается апостериорная оценка точности приближенных решений задачи об изгибе пластин Рейсснера — Миндлина. Оценка построена при помощи функционального подхода, основанного на строгих математических методах, в частности методах функционального анализа. Она справедлива для любых конформных аппроксимаций точных решений, что делает ее надежной. Она также является гарантированной, и неравенство не нарушается при практической реализации. Эти свойства делают данный метод контроля точности решений привлекательным для использования в инженерных расчетах, где некоторые вычислительные детали могут быть скрыты. В статье исследованы две независимые реализации оценки. Использование вычислительного эксперимента показало корректность работы алгоритмов и близость полученных результатов. Установлено, что для широкого диапазона значений толщины степень переоценки истинной величины погрешности остается приемлемой.

Ключевые слова: апостериорная оценка погрешности, метод конечных элементов, пластина Рейсснера — Миндлина

Ссылка при цитировании: Киселев К.В., Фролов М.Е., Чистякова О.И. Вычислительный эксперимент и верификация реализаций апостериорной оценки для пластин Рейсснера — Миндлина // Научно-технические ведомости СПбГПУ. Физико-математические науки. 2019. Т. 12. № 1. С. 128-141. DOI: 10.18721/JPM.12111

A POSTERIORI ERROR ESTIMATE FOR REISSNER - MINDLIN PLATES: VERIFICATION OF IMPLEMENTATIONS AND NUMERICAL TESTING

K.V. Kiselev, M.E. Frolov, O.I. Chistiakova

Peter the Great St. Petersburg Polytechnic University, St. Petersburg, Russian Federation

A posteriori error estimate for accuracy control of approximate solutions for the problem of Reissner — Mindlin plates bending has been analyzed in the paper. The estimate was constructed using the functional approach based on rigorous mathematical grounds, in particular, on methods of functional analysis. It is valid for all conforming approximations of exact solutions, and therefore, it is robust. The estimate is guaranteed in practical implementations due to reliability of the respective inequality. The above-mentioned properties of the method of error control are very desirable for engineering analysis, where some details of computations might be hidden. Our paper investigated two independent implementations of the estimate. Using specially constructed numerical tests, correctness of both implementation algorithms and similarity of the obtained results for all examples were shown. An overestimation of the true error was established to remain acceptable for a wide range of plate thickness values.

Keywords: a posteriori error estimate, finite element method, Reissner — Mindlin plate

Citation: K.V. Kiselev, M.E. Frolov, O.I. Chistiakova, A posteriori error estimate for Reissner — Mindlin plates: verification of implementations and numerical testing, St. Petersburg State Polytechnical University Journal. Physics and Mathematics. 12 (1) (2019) 128—141. DOI: 10.18721/JPM.12111

Введение

Работа посвящена апостериорному контролю точности решений в задачах теории пластин, которые нередко возникают в различных инженерных расчетах. Чем адекватнее выбранная модель описывает реальные физические процессы, происходящие при деформировании рассматриваемого объекта, тем ближе к эксперименту будут результаты проведенного расчета.

В настоящей работе исследуется проблема контроля точности метода конечных элементов (МКЭ) для модели линейно-упругих деформируемых пластин Реисснера — Миндлина [1]. Такая модель обычно применяется для пластин малой и средней толщины. В отличие от классической модели Кирхгоффа — Лява, описывающей пластины с малым отношением толщины к характерному размеру, модель Рейсснера — Миндлина позволяет рассматривать пластины средней толщины благодаря уходу от предположения о том, что волокна, перпендикулярные срединной плоскости пластин, сохраняют это свойство после приложения нагрузки.

Система уравнений равновесия для модели Рейсснера — Миндлина имеет следующий вид:

- Div(Cs(e)) = у, ^гуу = g,

в области О, (1)

у:= 2 (Уи -е),

где и — скалярное поле прогиба срединной плоскости; 0 — векторное поле поворотов вектора нормали к срединной плоскости при деформировании; у — векторное поле, соответствующее паре (и, 0); 8 — тензор деформаций; ? — толщина пластины; gfi — распределенная поперечная нагрузка; С — симметричный тензор четвертого ранга; О — плоская область, занимаемая пластиной;

Х =1 Ek/(1 + у) 2

(Е — модуль Юнга, V — коэффициент Пуассона, к — корректировочный коэффициент; во многих расчетах берется к = 5/6).

В настоящее время существует ряд известных коммерческих пакетов для проведения расчетов в рамках данной модели при помощи МКЭ. Поскольку программная реализация подобных продуктов, таких как ЛК8У8 (популярный среди инженеров

вычислительный пакет), является закрытой, оказывается, как правило, невозможным узнать во всех деталях, каким образом производился необходимый расчет. В итоге трудно оценить величину ошибки, если не разработан соответствующий строгий математический аппарат или не реализована оценка погрешности, полученная с его помощью. Существенные отклонения от экспериментальных данных могут быть обусловлены и неправильно выбранной моделью, и ошибками при ее применении, и низким качеством расчета при вполне адекватной модели.

Таким образом, возникает необходимость получения и численного исследования апостериорной оценки точности приближенного решения, которая являлась бы гарантированной и надежной, то есть не допускала бы недооценки погрешности при практической реализации и не зависела бы от скрытых деталей процедуры вычисления.

Первый вариант такой оценки был получен в работе [2] на основе функционального подхода и модифицирован в работе [3]. Обзор последних результатов, касающихся ее численной реализации, содержится в работе [4]. Историю развития функционального подхода с 1996 года по настоящее время можно проследить по источникам, указанным в монографиях [5 — 7] и недавно опубликованных статьях [8, 9]. Обзоры других подходов в применении к рассматриваемой задаче можно найти, например, в статьях [10, 11] и цитируемой там литературе.

Общее представление подобных апостериорных оценок имеет следующий вид:

II и - и || <

< М (и, D, у, у 2,..., с, С2,...),

(2)

где в левой части рассматривается норма отклонения неизвестного точного решения и от полученного приближенного решения и ; она выбирается исходя из постановки задачи; в правой находится функционал М — мажоранта отклонения. Аргументами функционала являются приближенное решение и, параметры задачи Д набор постоянных (СрС^,...). Важно, что значение этих констант определяется свойствами задачи, но не свойствами разбиения. Наконец, (у, у2,...) представляет собой набор свободных элементов.

Функционал М должен удовлетворять естественному условию — давать нулевое значение тогда и только тогда, когда приближенное решение совпадает с точным решением.

В то же время знание не только величины глобальной ошибки

приближенного решения, но также локального распределения ошибки на элементах разбиения представляется не менее полезным и важным. Оно позволяет выделить области с наибольшим отклонением приближенного решения от точного. Далее, вместо дробления всей сетки, производится необходимое улучшение лишь на тех элементах, где ошибка велика. Подобный подход называется адаптивным; он способствует сокращению вычислительных ресурсов, требуемых для получения решения необходимого качества.

Таким образом, полный цикл проведения расчета с использованием адаптивного подхода принимает следующий вид:

начальное разбиение области ^

получение решения ^

оценка глобальной ошибки

(выход из цикла расчета

при достижении требуемой

точности) ^ (3)

определение областей

с наибольшими локальными

ошибками ^

дробление сетки

в выделенных областях ^

расчет на новой сетке ^

получение более точного решения.

Для внедрения данного цикла в процесс решения реальных инженерных задач необходима уверенность в свойствах исследуемой оценки, а также определение набора параметров, при котором достигается надежный результат.

Ранее авторами данной статьи М.Е. Фроловым, а затем О.И. Чистяковой были выполнены две независимые реализации апостериорной оценки для пластин Рейсснера — Миндлина (см. обзор в работе [4]): на языке FORTRAN и в пакете MATLAB. Было проведено начальное исследование надежности результатов работы каждого из алгоритмов, дающих численную величину оценки. На данный момент необходимо расширить исследование поведения апостериорной

оценки путем рассмотрения большего числа примеров, чтобы оценить корректность использования данного инструмента, достоверность результатов и создать основу для развития метода.

Целью настоящей работы является сравнительный анализ двух реализаций вычисления апостериорной оценки точности приближенного решения задачи об изгибе пластины Рейсснера — Миндлина. Анализ включает построение примеров, имеющих различные особенности, проведение расчетов для разных значений толщины и сравнение как величин оценки, так и сопутствующих вспомогательных параметров.

Математическая постановка задачи

Обратимся к системе уравнений, описывающей модель Рейсснера — Миндлина для линейно-упругих пластин толщины t. Постановка записывается в терминах скалярной функции u=u(x) — прогиба срединной плоскости в точке x и векторной функции 0 = 0(х) — углов поворота вектора нормали к срединной плоскости. Срединная плоскость изначально занимает ограниченную связную область О ^ К2 с границей Г, непрерывной по Липшицу.

Постановка задачи: Найти пару (и,0) и соответствующее этой паре векторное поле у = у(х), удовлетворяющие системе уравнений (1).

Решение ищется в обобщенном смысле. Считаем, что g еЬ2 (О), тензор C симметричен, существует пара констант \ и такая, что выполняется двусторонняя оценка

Х|к|2 < Ск:к < Х2|к|2

Vк е М2^2 , | к |2 = к : к, где М 2т2 - пространство симметричных тензоров ранга 2 размерности 2.

В данном исследовании рассматриваются пластины из однородного изотропного материала, хотя это и не является принципиальным ограничением.

Рассмотрим границу Г области О и представим ее в виде двух составляющих: Гр, Гв ^ 0, на которой пластина жестко заделана, т. е. заданы условия на перемещение u = 0 и поворот 0 = 0, и Г5 = Г \ Гр, где пластина имеет свободные граничные условия (можно рассматривать и другие типы граничных условий).

Обобщенная постановка задачи об изгибе пластин Рейсснера — Миндлина:

Найти такую тройку (и, е, у) е Их ©х 0, где

и = {ше Ш2 (О)|ш = 0 на Г0 }, <0 = {феШ2 (О) х Ш2 (О) | ф = 0 на Г0 }, 0 = Ь2 (О) х Ь2 (О),

которая удовлетворяет соотношениям

| С е(е):е(ф^ О-

О

-|у-ф d О = 0, Уфе©;

О

|у-Уш d О =

О

= | gшdО, Уш е И;

О

|(К-1г2у-(Уи -е))-т dО = 0,

О

Уте

(4)

3 (и, е) = |Г 2 Се(е): е(е) +

+1 КГ21 Уи-е|2 -gu) d О.

(5)

еи = и - г/,

ее =е-е, еу =у-у.

Введем квадрат ошибки, выраженной через разность значений функционала энергии:

ё2 = 3(и/, е) - 3(и, е).

Показано (см. работу [2]), что У (и, 0) е S имеет место соотношение 1

ё2 = -(|||е'112 ^-1*2

где

Функционал энергии в задаче имеет следующий вид:

+К-1г 2||е?||О ),

|||ее|||2:= | Се(ее): е(ее) d О,

О

||е¥||О := /1 е? 12dО.

О

Далее введем в рассмотрение необходимые для построения апостериорной оценки свободные элементы: векторное поле у и несимметричный тензор к,_ представленный в виде его составляющих к = ^к1, к2 J. Они выбираются так, что

у, к1, к2 е Н(О^у) := : = {у е Ь2 (О) Ь2 (О) | divy е Ь2 (О)}.

После тождественных преобразований и применения известных неравенств, представленных в работе [3], итоговое неравенство принимает вид:

Соответствующая задача минимизации ставится на пространстве

8 := и х в.

Перейдем непосредственно к виду исследуемой апостериорной оценки. Пусть имеется пара конформных аппроксимаций (и, е) точного решения задачи (^ 0 ) в 8, тогда векторное поле

у = —-2 (У и -е)

аппроксимируется в пространстве Q элементом

У = —-2 (у и -е).

Введем в рассмотрение тройку отклонений каждой из компонент решения:

|||ее|||2 + К-1г2||е?||О < а2 + К-1г2Ь2, (6)

где

а = ||| СЛут(к) -ё(е)||| +

+с 11 skew(K) ||О + с2с3 (| О | || g + divy |

ч ч1/2

+ | Г5 | || у-п ) +

-с4 (| О | || у + [divK1, divK2] ||О

+1 Г5 | || [к1 -п, к2 -п]|Г)1/2,

Ь = У у-У||О +

+ сз ( | О | 11 g + divy ||О + 2 \1/2

+ |Г^ ||у-п||г,) ,

(7)

Константы в последней строке зависят только от геометрии области О, граничных условий и свойств материала, и возникают из следующих неравенств:

II Уф 110 < с2!!! ф |||2 !1ф|1О < III ф |||2:

ОI

О |

I ш !!о +;

1

< с.

ФI

I х 5

Уш 1

ш

2

О,

!г2 <

11 V

Ф112 <

< с42 ||| Ф |

Методология реализации вычислительного эксперимента

Рассмотрим детально шаги адаптивного алгоритма (3), который обычно используется для более эффективного уточнения приближенного решения.

На первом этапе задаются геометрия исследуемой области, толщина пластины, параметры материала, строится начальное разбиение. Далее производится расчет в подходящем вычислительном пакете, например А№У8.

После получения значений перемещения и поворотов в узлах сетки можно переходить к следующему блоку — применению исследуемой апостериорной оценки и определению величины глобальной ошибки. Она на всем решении определяется как корень из суммы квадратов локальных ошибок на каждом элементе, именуемых индикаторами. Если значение глобальной ошибки неудовлетворительно в рамках поставленной задачи, то рассматривается распределение индикаторов на элементах разбиения.

При наличии картины распределения локальных ошибок появляется возможность выделить области с наибольшим значением погрешности. Существует несколько стратегий выделения элементов (см., например, работы [12, 13]); они сводятся к определению порога значения индикатора, при котором глобально не допускается выделения слишком малого или слишком большого числа элементов.

Финальной частью каждого цикла адаптивного алгоритма выступает дробление сетки (а в общем случае — и ее укруп-

нение) в выбранных областях. Задача не является тривиальной и второстепенной, поскольку качество сетки влияет на точность расчета. Необходимо сохранить конформность сетки; углы сетки должны лежать в определенном диапазоне — наличие слишком малых и слишком больших углов негативно влияет на точность расчета; при наличии симметрии также крайне желательно сохранить ее характер. Варианты алгоритмов, решающих данную проблему для четырехугольных элементов рассмотрены, в частности, в работе [14].

Далее производится перерасчет решения задачи на новом разбиении. Данный адаптивный итерационный алгоритм подразумевает, что результат, полученный на измененной сетке, снова может быть оценен с точки зрения точности и распределения погрешности по расчетной области

Численные результаты и их обсуждение

В данном исследовании сравнительный анализ проводился на четырех примерах, проявляющих различные вычислительные особенности: круглаяпластина, косоугольная пластина, квадратная пластина с крупным квадратным отверстием, квадратная пластина с малым круглым отверстием.

Каждый рассмотренный случай обладает своим набором особенностей, которые могут выявить несовершенства реализаций. Например, квадратная пластина с большим квадратным отверстием имеет области концентраций напряжения вблизи вершин отверстия. Для круглой пластины известно аналитическое решение. Во всех примерах рассматривается равномерно распределенная нагрузка. Косоугольная пластина также встречается в литературе по апостериорным оценкам [15, 4]. В нашем исследовании она использовалась для демонстрации поведения ошибки при наложении жестких граничных условий только на верхнюю и нижнюю грани (две оставшиеся свободны), в отличие от других примеров, где рассматривается заделка по всей границе области.

Перейдем к описанию двух примеров. Отметим, что хотя объем статьи не позволяет включить результаты всех проведенных исследований, выводы, сделанные в итоге, ими полностью подтверждаются.

Создание геометрии области, наложение граничных условий, построение сетки и непосредственное решение производились

О

в пакете ANSYS Mechanical (APDL). Для проверки также использовался решатель, реализованный в MATLAB.

Для анализа эффективности оценки вводится понятие индекса эффективности IeJ, общепринятое в теории апостериорного контроля точности и определяемое как отношение мажоранты к норме отклонения приближенного решения от точного,

Ieff = M/1| u —u ||.

Таким образом, If показывает степень переоценки ошибки. Соответственно, чем его значение ближе к единице справа, тем лучше. Если аналитическое решение задачи неизвестно, за точное решение u принимается приближенное решение, вычисленное на очень мелкой сетке. Этот прием всегда чуть завышает значение Ieff) и в действительности оценивание происходит даже лучше, чем мы можем это показать.

Реализации сопоставляются по наиболее важным параметрам. Сначала рассматриваются константы c1, c2, c3, c4. Для вычисления c1 дополнительно рассматривается константа в неравенстве Корна ck. На важную роль данной константы в апостериорных оценках для задач теории упругости указывается, в частности, в монографии [6]. Тогда

С = скУ112(1 + v) / E

(для задач, где Г s = 0, ск = V2).

Далее сравнивается ошибка e — корень квадратный из левой части оценки (6); мажоранта M — корень квадратный из правой части оценки (6) и ее составляющие D, S, R, где

D = ||| СЛут(К)— e(9)||| +

II У — Y iin, S = ci ||skew(K)||n,

R = c2c3

(| Q | || g + divy

\ 1/2 ) +

+

+ X-1/2tc

^S 1 11 y •n

Q | || g + divy

1Гs 1 11 y •n ||2

1/2

+

+С

2n ii 2 Q

4 (| Q | || y + [divK1, divK2 ]

~1 ~9 2 \ 1/2 +1 Ts I || [K1 • n, к2 • n]||r2s) ;

последняя сравниваемая величина — индекс эффективности If

Степень дробления сетки обозначается как R., i е{0, 1, 2, 3, 4, 5}, где нулевой индекс соответствует начальному разбиению.

Для вычисления мажоранты М первоначально запускается процесс приближенного расчета констант (7). Их вычисление происходит один раз для каждого примера и далее может использоваться при получении распределения ошибки на всех сетках и толщинах. Вычисление констант носит асимптотический характер относительно дробления сетки Ri, а удовлетворительные значения получаются даже на относительно грубых сетках, что важно для перспектив практического применения метода.

Вспомогательный второй блок включает в себя поиск эталонного решения на мелких сетках и вычисление на его основе значения функционала энергии. Этот расчет производится только в целях постановки вычислительного эксперимента и в реальных инженерных задачах не требуется. Он просто заменяется расчетом функциональной апостериорной оценки на текущей сетке, что в вычислительном плане несравнимо менее трудоемкая задача.

Заключительный третий блок отвечает за непосредственное вычисление мажоранты М, где на вход поступают определенные c1, c2, c3, c4, которые несколько переоцениваются для надежности.

Приведем полученные результаты сравнения на примере круглой пластины и квадратной пластины с крупным квадратным отверстием.

Пример 1. Круглая пластина. Начнем с классического примера (см. ANSYS Verification Manual), для которого известно аналитическое значение прогиба в центре пластины.

На рис. 1 представлены результаты расчета деформации круглой пластины средней толщины t в пакете ANSYS на сетке R1. Параметры модели приведены в табл. 1.

Табл. 2 содержит значения вычисленных вспомогательных констант и позволяет сравнить результат работы обеих реализаций на различных сетках. Как было указано выше, наблюдается сеточная сходимость значений констант при уменьшении размеров элементов сетки. В силу нулевых граничных условий, точное значение константы ck известно и равно V2. Для дальнейших расчетов используются константы,

полученные на сетке R2, для надежности оцененные сверху с запасом (см. графу R2*).

Табл. 3 содержит основные расчетные результаты для сравнения, а именно мажоранту М и ее составляющие. Обратимся к верхней части таблицы, т. е. к случаю, когда Rad/t = 250. Видно, что определяющим является слагаемое D, так как в эту компоненту входят те части оценки, которые соответствуют слагаемым ошибки. Компоненты S и R должны быть близкими к нулю на полях, близких к точным, что и наблюдается. В шестой и двенадцатой строках таблицы приводятся значения /_, позволяющие оценить степень переоценки для обеих реализаций. Видно, что значения, вычисленные по двум алгоритмам (в шестой строке If = 1,60), совпали, а само это значение оказывается удовлетворительным, т. е. критической переоценки (на порядок или в несколько раз) не наблюдается.

Рис. 2 иллюстрирует распределение локальных индикаторов ошибки в области круглой пластины на сетке R2. Индикаторы идентичны и указывают на то, что ошибка распределена не только в зоне максимального прогиба пластины, но и вблизи границы области. Это соответствует структуре энергетической нормы в данной задаче (эта норма определена через градиенты, которые у края рассчитываются с меньшей точностью). Отсутствие полной симметрии связано с особенностями построения сетки.

Далее приведем расчетные результаты для тонкой круглой пластины толщиной t = 0,00005 м (Rad/t = 5000), где параметры материала и тип граничных условий остались прежними (см. табл. 1). Это позволяет оценить эффективность метода в зависимости от ключевого в данном случае параметра — толщины пластины.

На рис. 3 представлено приближенное решение. Видно, что максимальное значение прогиба срединной плоскости увеличилось, по сравнению с таковым на рис. 1. Изменение толщины пластины не влияет на геометрию области, следовательно, и на результат вычисления констант (7). Из сравнения значений в шестой и двенадцатой строках табл. 3 можно также заключить, что для тонкой пластины значение индекса If снизилось с 1,60 до 1,42, т. е. приблизилось к оптимальному — единице.

Расчетные распределения локальных индикаторов ошибки по области круглой пластины малой толщины t = 0,00005 м

Рис. 1. Результаты расчета деформации круглой пластины средней толщины в пакете

ANSYS на сетке R1; uz — распределение прогиба срединной плоскости; 9x, — распределения полей поворота вектора нормали к ней. Rad/t = 250 (см. данные табл. 1)

Таблица 1

Описание модели круглой пластины

Параметр Значение Особенность

Модуль Юнга E, Н/м2 2' 1011 1. Пластина жестко закреплена по всей границе 2. Построенная на области сетка имеет структуру общего вида

Коэффициент Пуассона v 0,3

Нагрузка, Н/м2 6585,175

Радиус Rad, м 0,25

Толщина t, м 10-3 ; 5 ' 10-5

Таблица 2

Результаты вычисления констант для круглой пластины в реализациях FORTRAN и MATLAB

Сетка Число Значение константы

элементов С1 C2 c3 C4

R0 32 1,39 1,23e - 05 9,15e - 07 2,09e - 06

1,42 1,25e - 05 9,01e - 07 2,06e - 06

R1 128 1.41 1.42 1,25e - 05 9,36e - 07 9,32e - 07 0,23 2,12e - 06 2,11e - 06

R2 512 1,42 9,42e - 07 9,41e - 07 2,13e - 06

r; 512 - 1,30e - 05 9,50e - 07 0,30 2,20e - 06

Примечания. 1. Верхние числа в клетках таблицы относятся к реализации FORTRAN, нижние — к MATLAB; при совпадении данных здесь и далее указывается одно число. 2. Данные для сетки, помеченной звездочкой, относятся к оценке сверху с запасом (для надежности). 3. Rad/t = 250.

Таблица 3

Контрольные значения неравенства (6) для круглых пластин средней и малой толщины в двух реализациях на сетке R

Составляющая Контрольное значение

мажоранты (и мажоранта) FORTRAN MATLAB

Rad/t = 250

e 0,712е + 02

M 0,114e + 03

D 0,111e+ 03

S 0,271e + 01 0,272e + 01

R 0,723e + 00

1eff 1,60

Rad/t = 5000

e 0,114e+ 12

M 0,162e + 12

D 0,162e + 12

S 0,217e + 08

R 0,568e + 07

1eff 1,42

(Rad/t = 5000) на сетке R имеют вид, близкий представленному на рис. 2, поэтому в данной статье не приведены. На основе полученных данных можно сделать вывод о корректности работы алгоритмов для достаточно широкого диапазона значений толщины, что подтверждается и другими выполненными исследованиями.

Пример 2. Квадратная пластина с крупным квадратным отверстием. Следующий пример (табл. 4) был рассмотрен для того, чтобы проследить поведение оценки в областях концентрации напряжений. Естественно, что при данной конфигурации пластины подобные области образуются в районе вершин внутреннего отверстия.

FORTRAN indicator (а)

-0.2 -0.1 0 0.1 X, m

MATLAB indicator (b)

-0.2 -0.1 о 0.1 x, m

uz, m

■0.2 -0.1 о 0.1 x, m

-0.2 -0.1 о 0.1 x, m

■0.2 -0.1 о 0.1 x, m

Рис. 2. Расчетные распределения локального

индикатора ошибки на элементах области Рис. 3. Расчетные результаты, аналогичные

круглой пластины в реализациях FORTRAN представленным на рис. 1, но для круглой

(а) и MATLAB (b) на сетке R2; Rad/t = 250 пластины малой толщины; Rad/t = 5000 (см. (см. Пример 1) табл. 1).

Таблица 4

Описание модели квадратной пластины с крупным квадратным отверстием

Параметр Значение Особенность

Модуль Юнга Е, Н/м2 2' 109 1. Пластина жестко закреплена по внешней и внутренней границе 2. Элементы сетки имеют форму квадрата 3. Решение имеет зоны концентрации напряжения

Коэффициент Пуассона V 0,3

Нагрузка, Н/м2 6585,175

Сторона квадрата, м внешняя А внутренняя а 12' 10-3 4'10-3

Толщина 1, м 10-4

Таблица 5

Результаты вычисления констант в двух реализациях для квадратной пластины с крупным квадратным отверстием

Сетка Число элементов Реализация Значение константы

С1 С2 c3 C4

512 FORTRAN MATLAB 141 1,42 1,25e - 04 1,46e - 07 0,12 1,29e - 05 1,28e - 05

R* - 1,30e - 04 1,50e - 07 0,15 1,30e - 05

Примечания. 1. Данные для сетки, помеченной звездочкой, относятся к оценке сверху с запасом (для надежности). 2. Геометрический параметр А/? = 120.

Таблица 6

Контрольные значения неравенства (6) в двух

реализациях на сетке Я2 для квадратной пластины с крупным квадратным отверстием

Составляющая мажоранты (и мажоранта) Контрольное значение

FORTRAN MATLAB

е 0,313e + 04

М 0,529e + 04 0,534e + 04

В 0,519e + 04 0,515e + 04

£ 0,413e + 02 0,348e + 02

Я 0,567e + 02 0,161e + 03

1,69 1,70

а)

b)

Рис. 4. Результаты расчета деформации квадратной пластины с крупным квадратным отверстием; использован пакет ЛК8У8 на сетке Я.

Рис. 5. Расчетные распределения локального индикатора ошибки на элементах области в реализациях FORTRAN (a) и MATLAB (b) на сетке R2; A/t = 120 (см. Пример 2)

На рис. 4 изображены результаты расчета деформации квадратной пластины с крупным квадратным отверстием. Расчет выполнен в пакете ЛК8У8 на сетке Я1. Представлены распределения величин и, 9х и по области пластины. Как и в Примере 1, иг — прогиб срединной плоскости; 0^ , 0^ — распределения полей поворота.

Табл. 5 содержит значения вычисленных вспомогательных констант. Для дальнейших расчетов используются константы, полученные на сетке Я1, которые берутся с запасом (данные в строке Я1*).

Рассмотрим результаты сравнения контрольных значений оценки на сетке Я2. Из

Обозначения величин те же, что на рис. 1 данных табл. 6 видно, что значения 1^ для

двух реализаций практически совпали, су-

щественной переоценки не наблюдается.

Обратимся к данным, изображенным на рис. 5. Результаты реализаций совпадают, распределения индикатора практически идентичны. Зоны с минимальной ошибкой располагаются в углах внешней границы. В силу граничных условий, в этих местах узлы имеют минимальное перемещение. Зоны с максимальной ошибкой расположены в областях, где напряжения при приложении нагрузки имеют сингулярность.

Заключение

В рамках представленного исследования, для пластин Рейсснера — Миндлина было проведено сравнение результатов, полученных при помощи существующих на данный момент реализаций алгоритмов вычисления функциональной апостериорной оценки [3]. Проведен вычислительный эксперимент в применении к моделям пластин различных геометрических конфигураций области, значений толщины, параметров

материалов и структуры сетки.

На основе полученных результатов сделаны следующие выводы.

1. Для всех примеров с краевыми условиями жесткой заделки по всей границе, обе реализации дали близкие значения индекса эффективности; при этом величина переоценки ошибки удовлетворительная.

2. Распределения индикатора ошибки по расчетной области оказались практически идентичными для обеих реализаций, что также является подтверждением их корректности.

В дальнейшем представляет интерес применить подход к адаптивному решению прикладных задач с использованием ресурса суперкомпьютерного центра «Политехнический».

Работа М.Е. Фролова и О.И. Чистяковой выполнялась при финансовой поддержке гранта Президента Российской Федерации МД-1071.2017.1.

СПИСОК ЛИТЕРАТУРЫ

1. Falk R.S. Finite elements for the Reissner — Mindlin plate // Mixed finite elements, compatibility conditions, and applications. D. Boffi, L. Gastaldi (Eds). Berlin, Heidelberg: Springer, 2008. Pp. 195-232.

2. Репин С.И., Фролов М.Е. Об оценке отклонений от точного решения задачи о пластине Рейсснера — Миндлина // Краевые задачи математической физики и смежные вопросы теории функций. Записки научных семинаров ПОМИ. СПб, 2004. Т. 310. С. 145— 157.

3. Фролов М.Е. Надежный апостериорный контроль точности решений задач об изгибе пластин Рейсснера — Миндлина // Сеточные методы для краевых задач и приложения. Материалы Десятой международной конференции. Казань, 2014. С. 610—615.

4. Frolov M., Chistiakova O. Adaptive algorithm based on functional-type a posteriori error estimate for Reissner-Mindlin plates // Advanced Finite Element Methods with Applications. Selected papers from the 30th Chemnitz Finite Element Symposium, Lecture Notes in Computational Science and Engineering (LNCSE). 2019. Vol. 128. Ch. 7.

5. Neittaanmaki P., Repin S. Reliable methods for computer simulation: Error control and a posteriori estimates. Amsterdam: Elsevier, 2004. 304 p.

6. Repin S. A posteriori estimates for partial differential equations. Berlin: de Gruyter, 2008. 316 p.

7. Mali O., Neittaanmäki P., Repin S. Accuracy verification methods: Theory and algorithms. Springer Science+Business Media Dordrecht, 2014. XIII+355 p.

8. Langer U., Repin S., Samrowski T. A posteriori estimates for a coupled piezoelectric model // Russian Journal of Numerical Analysis and Mathematical Modelling. 2017. Vol. 32. No. 4. Pp. 259-266.

9. Repin S., Valdman J. Error identities for variational problems with obstacles // ZAMM — Journal of Applied Mathematics and Mechanics. 2018. Vol. 98. No. 4. Pp. 635—658.

10. Bösing P.R., Carstensen C. Discontinuous Galerkin with weakly over-penalized techniques for Reissner — Mindlin plates // Journal of Scientific Computing. 2015. Vol. 64. No. 2. Pp. 401—424.

11. Kumar B., Somireddy M., Rajagopal A.

Adaptive analysis of plates and laminates using natural neighbor Galerkin meshless method // Engineering with Computers. 2018. Vol. 35. No. 1. Pp. 201—214.

12. Dörfler W. A convergent adaptive algorithm for Poisson's equation // SIAM Journal on Numerical Analysis. 1996. Vol. 33. No. 3. Pp. 1106—1124.

13. Verfurth R. A review of a posteriori error estimation and adaptive mesh-refinement techniques. John Wiley & Sons, 1996. 127 p.

14. Караваев A.C., Копысов С.П. Перестроение неструктурированных четырехугольных и смешанных сеток // Вестник Удмуртского университета. Серия Математика. Механика.

Компьютерные науки. 2013. Вып. 4. С. 62—78. 15. Carstensen C., Xie X., Yu G., Zhou T. A

priori and a posteriori analysis for a locking-free low order quadrilateral hybrid finite element for Reissner — Mindlin plates // Computer Methods in Applied Mechanics and Engineering. 2011. Vol. 200. No. 9-12. Pp. 1161-1175.

Статья поступила в редакцию 29.12.2018, принята к публикации 02.02.2019.

СВЕДЕНИЯ ОБ АВТОРАХ

КИСЕЛЕВ Кирилл Владимирович — магистрант кафедры прикладной математики Института прикладной математики и механики Санкт-Петербургского политехнического университета Петра Великого, Санкт-Петербург, Российская Федерация.

195251, Российская Федерация, г. Санкт-Петербург, Политехническая ул., 29 kvladimirovich10@gmail.com

ФРОЛОВ Максим Евгеньевич — доктор физико-математических наук, директор Института прикладной математики и механики, заведующий кафедрой прикладной математики Санкт-Петербургского политехнического университета Петра Великого, Санкт-Петербург, Российская Федерация.

195251, Российская Федерация, г. Санкт-Петербург, Политехническая ул., 29 frolov_me@spbstu.ru

ЧИСТЯКОВА Ольга Игоревна — аспирантка кафедры прикладной математики Института прикладной математики и механики Санкт-Петербургского политехнического университета Петра Великого, Санкт-Петербург, Российская Федерация.

195251, Российская Федерация, г. Санкт-Петербург, Политехническая ул., 29 Chistiakova.Olga@gmail.com

REFERENCES

[1] R.S. Falk, Finite elements for the Reissner — Mindlin plate, Mixed finite elements, compatibility conditions, and applications, Springer, Berlin, Heidelberg (2008) 195-232.

[2] S.I. Repin, M.E. Frolov, Estimation of deviations from the exact solution for the Reissner — Mindlin plate problem, Journal of Mathematical Sciences. 132 (3) (2006) 331—338.

[3] M.E. Frolov, Reliable a posteriori error control for solutions to problems of Reissner — Mindlin plates bending, Proceedings of 10th International Conference Mesh methods for boundary-value problems and applications, Kazan University Publishing House, Kazan, Russia (2014) 610—615.

[4] M. Frolov, O. Chistiakova, Adaptive algorithm based on functional-type a posteriori error estimate for Reissner-Mindlin plates // Advanced Finite Element Methods with Applications. Selected papers from the 30th Chemnitz Finite Element Symposium, Lecture

Notes in Computational Science and Engineering (LNCSE). 128, Ch. 7 (2019).

[5] P. Neittaanmaki, S. Repin, Reliable methods for computer simulation: Error control and a posteriori estimates, Elsevier, Amsterdam, 2004.

[6] S. Repin, A posteriori estimates for partial differential equations, de Gruyter, Berlin, 2008.

[7] O. Mali, P. Neittaanmaki, S. Repin, Accuracy verification methods: Theory and algorithms, Springer Science+Business Media Dordrecht, 2014.

[8] U. Langer, S. Repin, T. Samrowski, A posteriori estimates for a coupled piezoelectric model, Russian Journal of Numerical Analysis and Mathematical Modelling. 32 (4) (2017) 259— 266.

[9] S. Repin, J. Valdman, Error identities for variational problems with obstacles, ZAMM — Journal of Applied Mathematics and Mechanics. 98 (4) (2018) 635—658.

[10] P.R. Busing, C. Carstensen, Discontinuous Galerkin with weakly over-penalized techniques for Reissner — Mindlin plates, Journal of Scientific Computing. 64 (2) (2015) 401—424.

[11] B. Kumar, M. Somireddy, A. Rajagopal, Adaptive analysis of plates and laminates using natural neighbor Galerkin meshless method, Engineering with Computers. 35 (1) (2018) 201— 214.

[12] W. Dörfler, A convergent adaptive algorithm for Poisson's equation, SIAM Journal on Numerical Analysis. 33 (3) (1996) 1106—1124.

[13] R. Verfurth, A review of a posteriori

error estimation and adaptive mesh-refinement techniques. John Wiley & Sons, 1996.

[14] A.S. Karavaev, S.P. Kopysov, A refinement of unstructured quadrilateral and mixed meshes, Vestnik Udmurtskogo Universiteta. Matematika. Mekhanika. Komp'yuternye Nauki. (4) (2013) 62-78.

[15] C. Carstensen, X. Xie, G. Yu, T. Zhou, A

priori and a posteriori analysis for a locking-free low order quadrilateral hybrid finite element for Reissner-Mindlin plates, Computer Methods in Applied Mechanics and Engineering. 200 (9-12) (2011) 1161-1175.

Received 29.12.2018, accepted 02.02.2019.

THE AUTHORS

KISELEV Kirill V.

Peter the Great St. Petersburg Polytechnic University

29 Politechnicheskaya St., St. Petersburg, 195251, Russian Federation

kvladimirovich10@gmail.com

FROLOV Maxim E.

Peter the Great St. Petersburg Polytechnic University

29 Politechnicheskaya St., St. Petersburg, 195251, Russian Federation

frolov_me@spbstu.ru

CHISTIAKOVA Olga I.

Peter the Great St. Petersburg Polytechnic University

29 Politechnicheskaya St., St. Petersburg, 195251, Russian Federation

Chistiakova.Olga@gmail.com

© Санкт-Петербургский политехнический университет Петра Великого, 2019