CC BY
15
20. Дорофеев, А.Н. Авиационные боеприпасы |Текст| / А.Н. Дорофеев, А.П. Морозов, P.C. Саркисян,— М.: Изд-во Военно-возд. инж. академии им. Н.Е. Жуковского,— 1978,— 446 с.
21. Вологжанин, О.Ю. Границы областей применения моделей поведения преграды при воздействии ударников [Текст] / О.Ю. Вологжанин, О.Ю. Вшивков, Н.А. Рыбаков // Вестник ИжГТУ,— 2010,- № 2 (46).- С. 16-18.
УДК 51 9.63: 539.3
М.Е. Фролов
РЕАЛИЗАЦИЯ ФУНКЦИОНАЛЬНОГО ПОДХОДА К АПОСТЕРИОРНОМУ КОНТРОЛЮ ТОЧНОСТИ РЕШЕНИЙ ТРЕХМЕРНЫХ ЗАДАЧ ТЕОРИИ УПРУГОСТИ
Значительный интерес к проблеме построения надежных методов, позволяющих апостери-орно контролировать погрешность вычисленных приближенных решений, возник к началу 80-х годов XX века. Он был вызван в первую очередь необходимостью проверки результатов, получаемых в процессе математического и инженерного моделирования.
Уже к концу XX века в вычислительной математике сформировалось отдельное направление, в рамках которого на основе метода конечных элементов было предложено несколько подходов к построению апостериорных оценок погрешности (см. монографии [1—6] и цитируемую там литературу). Однако подходы, ставшие за эти десятилетия классическими, основаны на дополнительных предположениях о свойствах приближенного решения, в частности, это должно быть точное решение соответствующей конечномерной задачи (галёркинская аппроксимация). Такие ограничения лишают их необходимой общности.
Идея создания более универсального подхода на основе привлечения методов математической физики, вариационного исчисления и функционального анализа, нашедшая отражение еще в работах С.Г. Михлина, на современном этапе получила свое существенное развитие в исследованиях профессора С. И. Репина и его коллег в 1996—2010 гг. Для многих задач получены мажоранты погрешности, удовлетворяющие важным требованиям — они дают надежные (гарантированные) оценки точности интегральной величины погрешности и подходят для широкого класса приближенных решений, построенных
различными методами. Теоретические основы подхода, названного функциональным, интенсивно развиваются и обобщаются уже более десяти лет (см. монографию [6] и цитируемую там литературу). Тем не менее вопрос его практического применения для контроля точности решений задач механики исследован недостаточно. Особенно это касается реализации предлагаемых методов как независимой надстройки над коммерческими пакетами для инженерных вычислений, которые широко используются в современной инженерной практике. Надежный контроль погрешности таких решений другими (нефункциональными) методами для большинства задач вообще видится крайне затруднительным. Основная причина как раз заключается в том, что стандартные подходы к построению апостериорных оценок основаны на специальных предположениях о свойствах приближенного решения. Гарантировать выполнение всех требований в случае коммерческого программного обеспечения с закрытым кодом невозможно. С вычислительной точки зрения пакеты все-таки выглядят как «черный ящик». Отдельные детали реализации скрыты или описаны не в полном объеме.
Данная работа связана с реализацией функционального подхода к оценке точности решений трехмерных задач линейной теории упругости, полученных в пакете ANSYS. Для этого на языке программирования FORTRAN создан внешний вычислительный модуль, а эффективность применения методики проверена численными экспериментами.
Математическая постановка задачи
Исходная классическая задача заключается в определении в трехмерной области ^ поля перемещений и, поля напряжений а и поля деформаций е, которые связаны системой соотношений
сИУСТ + Т'^О;
где
а =
е(и) = ^ (Vu + (Vuf),
(1)
и = м0 на Г
и
(2)
¡Щи): e(w)rffi = ¡F -wdX -
X X
+ J g-wdT,
(4)
v0 = jwe^2'(Q,Ä3)
u = u0 на Г,
и i
w = 0 на Г,
и |
Функционал энергии, соответствующий задаче, выглядит следующим образом:
1
где Г— плотность объемных сил; Ь — тензор упругих модулей, связывающий напряжения и деформации.
Будем предполагать, что среда однородна и изотропна; тогда в линейном случае этот тензор определяется двумя постоянными, связанными со свойствами материала. Это могут быть, например, коэффициент Пуассона и модуль Юнга. К системе уравнений (1), выполняющихся в области, занимаемой упругим телом, добавляются граничные условия в терминах перемещений
J(u) = — ¡Щи): e(u)dX --¡F -udO- J g-udT.
(5)
Если рассмотреть произвольное приближенное решение V е V задачи (4), то несложно показать, что энергетическая норма погрешности этого решения выражается через разность значений функционала на приближенном и точном решениях, а именно
vll|2 := - v) : е(м - v)dX =
= 2(J(v)-J(u)).
(6)
и в терминах напряжений а
где Гц, ГА, — две непересекающиеся части границы, объединение которых представляет собой всю границу области П; м0, g— перемещения, заданные на первой части границы, и поверхностные силы, заданные на второй части границы, соответственно; п — нормаль к поверхности тела.
Основой для вывода апостериорной оценки для задачи (1) — (3), как и для применения метода конечных элементов, служит ее обобщенная (слабая) постановка. При этом рассматривается решение так называемой прямой задачи, в которой основной неизвестной характеристикой является поле перемещений, а поля деформаций и напряжений определяются по нему. Эта задача выглядит следующим образом:
найти элемент и е V такой, что
Если подставить точное решение и в функционал /, то в силу слабой постановки (4), мы получим функционал в виде
J(u) = -— ¡Щи): г(u)dX =
= — ¡c(u):e(u)dQ.
(7)
Таким образом, значение функционала на точном решении дает значение энергии деформации с обратным знаком. Таким же свойством обладают галёркинские аппроксимации — точные решения конечномерных задач, соответствующих задаче(4).
В работе [7] с привлечением методов теории двойственности вариационного исчисления A.B. Музалевским и С.И. Репиным получена апостериорная оценка функционального типа для задачи (4) и ее двумерных аналогов (плоское напряженное состояние, плоская деформация, осесимметричная задача). В ней также приведены два примера численного расчета для случая плоских задач.
Если ввести обозначения для невязок соотношения (3) и первых двух соотношений (1) как
Z?(a):=diva + /' в X; 6(v,a):=a-Ie(v) в X; (8)
r(a)\=an-g на Г^,
то оценка погрешности, предложенная в работе [7], может быть преобразована к виду
где
D(v,c):= Jr1 Ю(у,о): ^(v,c); (10) п
cX — константа, не зависящая от разбиения области конечными элементами и удовлетворяющая неравенству
H|X+|rf <Cn|||w|||2, VweV0;
символами ||...||п, \\-\\Гу обозначены стандартные нормы в пространствах функций, суммируемых с квадратом в смысле интеграла Лебега, соответствующих области и части ее границы.
С точки зрения практической реализации более удобной является форма оценки погрешности, в которой стоят квадраты норм без извлечения корня. С учетом того, что поиск значения константы сп , близкого к оптимальному, для произвольной области является непростой задачей, в численных исследованиях упор сделан на
вычисление индикатора погрешности <jD(v,a). При этом уточненное поле напряжений, необходимое для его вычисления, получается из минимизации функционала следующего вида:
D(v,g) + a, ||i?(c)||X + а2 ||r(c)|^ ,
где а,, а2 - малые весовые множители, позволяющие учесть в функционале независимо уравнения равновесия и граничные условия в напряжениях.
Численные результаты
Важное место в применении любых вычислительных технологий занимает их верификация. Поскольку коммерческие пакеты (такие как ANSYS) разрабатываются большим коллективом профессионалов — инженеров, вычислителей,
программистов, и подвергаются тестированию в промышленных масштабах, то их можно брать за основу при проверке собственного программного кода. Поэтому в процессе реализации метода конечных элементов, необходимого для вычисления функционала оценки погрешности, была проведена отладка кода при помощи сравнения по ряду основных характеристик решения исходной задачи в пакете с результатами работы авторской программы. В качестве характеристик были выбраны узловые значения перемещений, которые являются базовыми неизвестными величинами, подлежащими определению, а также значение энергии деформации (как производная величина, включающая градиент перемещений). Для эффективной и быстрой отладки программы были разработаны простые тестовые задачи, включающие в себя определенные характерные особенности.
В силу ограниченности объема публикации ниже представлены только три примера, которые, тем не менее, позволяют сделать и обосновать все необходимые выводы. Примеры направлены в первую очередь на оценку эффективности методики апостериорного контроля погрешности. Для этого используется равенство (6), в котором точное решение в правой части заменено на приближенное решение, вычисленное на более мелкой сетке. Такая замена позволяет достаточно хорошо оценить истинную погрешность; при этом мы лишь незначительно занижаем ее величину. С другой стороны, такой подход дает возможность рассматривать любые примеры и даже те, в которых гладкое аналитическое решение невозможно построить или его не существует; а таких задач, как известно, большинство.
Отметим, что в пакете ANSYS реализован один из классических методов построения апостериорных оценок, который был предложен в работе [8] и нашел достаточно широкое применение в инженерной практике для случая простых задач. По этой причине кажется разумным провести сравнение реализованной функциональной методики оценки погрешности и стандартного средства пакета.
Пример 1. В данном примере рассматривается задача, в которой присутствуют только условия в перемещениях в области, изображенной на рис. 1. Известно, что наличие входящих углов в области решения задачи может влиять на качество оценки погрешности. На рис. 1 также пред-
Рис. 1. Геометрия расчетной области
и начальное разбиение конечными элементами (пример 1)
ставлено начальное разбиение конечными элементами. Последующие сетки получены из предыдущих разбиений путем их равномерного измельчения по каждой стороне в два раза (Л1 из 7?0, Я2 из и т. д.). Таким образом, число элементов каждый раз увеличивается в восемь раз.
Для оценки эффективности методики оценки погрешности вычисляется решение на мелкой сетке, обозначенной REF. Основные параметры сеток указаны в таблице. Отметим, что максимальная размерность задачи, которую пришлось решить при оценке погрешности, составила почти 0,5 млн. уравнений, а решение задачи на мелкой сетке связано уже с системой порядка 2 млн. уравнений.
Данные для рассматриваемой задачи следующие: перемещения всех граничных узлов фиксированы, ускорение az = 9,8 м/с2 (учет собственного веса конструкции), модуль Юнга Е= 10,7 ГПа, коэффициент Пуассона v = 0,15; плотность р = 1,83-103 кг/м3.
Результаты, полученные при помощи реализованного решателя, практически совпадают с теми, что дает решатель ANSYS. Отклонение в значениях поля перемещений составляет менее 0,01 %. Данные об оценке погрешности также приведены в таблице. В ней представлена оценка ошибки снизу, полученная при помощи решения, вычисленного на мелкой REF-сетке (столбец 4). Поскольку такая оценка практически совпадает
Данные о сетках и эффективности оценок
Разбиение Число Истинная погрешность, % Индикатор погрешности, %
узлов элементов Функциональный ANSYS
Пример 1
m 1820 1404 15,6 16,5 22,2
m 11200 9672 9,4 10,8 14,3
R2 81095 75168 4,2 5,5 7,9
REF 732253 (~2 млн. уравнений) 707808
Пример 2
R0 3345 2368 9,6 10,8 13,1
m 22689 18944 5,3 5,6 7,7
REF 166209 (~0,5 млн. уравнений) 151552
Пример 3
R0 369 160 69,4 50,7 75,8
RI 2025 1280 43,5 39,6 59,4
R2 13041 10240 23,6 23,4 37,3
REF 698049 (~2 млн. уравнений) 655360
с истинной погрешностью, мы пренебрегаем этой разницей и сравниваем остальные данные с ней; тем более, что такая замена может привести только к недооценке , но никак не к переоценке эффективности исследуемой методики. В последнем столбце дается значение индикатора относительной точности решения, которое выдает пакет. Вычисляется этот индикатор при помощи методики, предложенной в работе [8].
Из данных таблицы можно сделать вывод, что наличие входящих углов незначительно влияет на качество полученных оценок погрешности. ANSYS несколько завышает истинную величину погрешности. Отметим, что хотя нарушения неравенства «истинная погрешность < оценка» не происходит, используемые индикаторы не являются гарантированными верхними оценками погрешности. Они стали бы гарантированными, если применять неравенство (9), но с практической точки зрения это связано с некоторыми затруднениями.
Пример 2. Во втором примере рассматривается стальная конструкция, закрепленная с одного конца и нагруженная с другого постоянным давлением (рис. 2). Фиксировано перемещение точек торцевой грани X = 0. Давление на отмеченные участки поверхности составляет 1 МПа. Учитывается также и собственный вес конструкции а у = —9,8 м/с2. Параметры материала: Е= 207 ГПа; v = 0,30; Р = 7,83-103 кг/м3.
Параметры конечно-элементных разбиений и результаты расчетов для этого примера также представлены в таблице. Из анализа приведенных данных можно сделать вывод о том, что качество оценки погрешности оказывается высоким, поскольку используемые объемные конечные элементы (SOLID) хорошо подходят для решения данной задачи. Поэтому и точность приближенного решения быстро достигает отметки почти в 95 %, считающейся в инженерной практике вполне достаточной.
Пример 3. Заключительный пример — это по сути незначительно модифицированный VM184(ANSYS 12.0 Verification Manual), в котором рассматривается стержень малого прямоугольного сечения. Его размеры вдоль каждой из осей составляют, соответственно: 1К= 152,4 мм, lY= 5,08 мм, lz= 2,54 мм. Фиксированы перемещения узлов одного из торцов, а вблизи свободного торца приложено давление 0,23 МПа в направлении оси 0Z. Свойства материала еле-
..............
о
Рис. 2. Геометрия области и граничные условия для примера 2 (стрелка вверх указывает зону закрепления, стрелка вниз — зону и направление приложения давления)
дующие: модуль Юнга Е— 68,95 ГПа, коэффи-v
случае не требуется). При таких параметрах прогиб левого края должен быть близок к значению, приведенному в оригинале (VM184), а оно составляет 10,97 мм.
В данном примере обращает на себя внимание существенная недооценка погрешности на грубой сетке (см. таблицу — сетка R0 160 элементов). Это вполне объяснимо, поскольку для вычисления индикатора нам необходимо решить систему уравнений и определить уточненное поле напряжений. Но если сетка грубая, то это поле не может достаточно точно удовлетворить и граничным условиям в напряжениях, и уравнениям равновесия внутри области. Видно, что с измельчением сетки ситуация исправляется. Отметим, что для задачи об изгибе стержня Тимошенко в работе автора [9] получена и численно протестирована другая оценка функционального типа. Если исходную задачу можно свести к одномерной, то следует решать ее посредством балочных элементов (ВЕАМ) и применять соответствующую оценку погрешности. В работе [9] решение в ANSYS получается припомощи ВЕАМ-элементов, и также контролируется его точность посредством внешней программы, реализованной на языке FORTRAN. Результаты указывают на высокую эффективность методики и точность гарантированных верхних оценок даже на грубых сетках.
Итак, в результате исследования был реализован на языке FORTRAN и отлажен вычисли-
тельный модуль оценки погрешности решений трехмерных задач линейной теории упругости. Методика была применена к приближенным решениям, вычисленным в пакете AN SYS 12.0. При сравнении со стандартным индикатором погрешности ANSYS также следует учесть тот факт, что предлагаемые в рамках функционального подхода методы могут быть обобщены на случай нелинейных задач и задач, в которых нарушается требование соответствия точному решению конечномерной задачи. В этих ситуациях стандартный индикатор пакета недоступен.
При помощи разработанной серии численных экспериментов различной направленности модуль оценки погрешности был протестирован
СПИСОК J
ffirth, R. A review of a posteriori error estimation and adaptive mesh-refinement techniques furth.— Chichester, Stuttgart: John Wiley & Sons, B.G. Teubner, 1996,- 127 p.
2. Ainsworth, M. A posteriori error estimation in finite element analysis [Text] / M. Ainsworth, J.T. Oden.— New York: John Wiley & Sons, 2000,- 240 p.
3. Babuska, I. The finite element method and its reliability [Text] / I. Babuska, T. Strouboulis.— New York: The Clarendon Press Oxford University Press, 2001,-802 p.
maki, P. Reliable methods for computer simulation. Error control and a posteriori
maki, S.I. Repin. — Studies in Mathematics and its Applications 33.-Amsterdam: Elsevier, 2004.— 305 p.
5. Ladeveze, P. Mastering calculations in linear and nonlinear mechanics [Text] / P. Ladeveze, J.-P. Pelle.— Mechanical Engineering Series.—New York: Springer, 2005.-413 p.
на задачах, анализ которых приводит к решению систем до 0,5 млн. уравнений. При этом решение системы, связанной с исходной задачей, на мелкой сетке, необходимое нам для оценки эффективности метода, потребовало существенных вычислительных затрат, поскольку число уравнений доходило до 2 млн.
Результаты подтвердили высокую эффективность предлагаемых методик расчета апостериорных оценок функционального типа, по крайней мере, в области их обоснованного применения.
Исследование выполнено при поддержке Правительства Санкт-Петербурга в рамках конкурса грантов 2010 года для молодых ученых и молодых кандидатов наук (диплом ПСП N910699).
6. Repin, S.I. A posteriori estimates for partial differential equations [Text]: Radon Series on Computational and Applied Mathematics 4 / S.I. Repin.— Berlin: de Gruyter, 2008.— 316 p.
7. Muzalevsky, A.V. On two-sided error estimates for approximate solutions of problems in the linear theory of elasticity [Text] / A.V. Muzalevsky, S.I. Repin // Russian Journal of Numerical Analysis and Mathematical Modelling.— 2003.— Vol. 18, № 1.— P. 65-85.
8. Zienkiewicz, O. C. A simple error estimator and adaptive procedure for practical engineering analysis [Text] /O.C. Zienkiewicz, J.Z. Zhu // Internat. J. Numer. Methods Engrg.— 1987,— Vol. 24,— № 2,- P. 337-357.
9. Frolov, M.E. Functional a posteriori error estimates for certain models of plates and beams [Text] / M.E. Frolov // Russian Journal of Numerical Analysis and Mathematical Modelling.— 2010,-Vol. 25,-№2,-P. 117-129.
УДК 539.3
В.В. Тихомиров
ТРЕЩИНА ПРОДОЛЬНОГО СДВИГА, ЧАСТИЧНО ПРОНИКАЮЩАЯ В УПРУГОЕ КРУГОВОЕ ВКЛЮЧЕНИЕ С ПОКРЫТИЕМ
Проблема взаимодействия трещин с включениями различной формы в рамках линейной механики разрушения привлекала и привлекает
внимание многих исследователей. С одной стороны, интерес к ней вызывается потребностями механики неоднородных сред (геомеханики)
