CC BY
199
УДК 531
Е. Е. Красновский, А. С. Шадский
ТОЧНОСТЬ РЕШЕНИЯ ЗАДАЧ МЕХАНИКИ В ANSYS ПРИ НАЛИЧИИ ИЗГИБНЫХ НАПРЯЖЕНИЙ
Исследовано влияние порядка аппроксимации, числа точек интегрирования, наличия дополнительных форм перемещений конечных элементов, применяемых в программном комплексе ANSYS, а также густоты сетки на точность решения задач механики, в которых изгибные напряжения имеют существенное значение. На примере продемонстрированы источники погрешностей shear locking и hourglassing.
E-mail: [email protected]; [email protected]
Ключевые слова: метод конечных элементов, ANSYS, погрешность, точность решения, изгиб, shear locking, hourglassing.
Современные конечно-элементные программные комплексы, особенно такие известные, как ANSYS или ABAQUS, предлагают широкие возможности по решению задач механики, в том числе большое число типов конечных элементов и их опций. Например, только для решения задач механики деформируемого твердого тела пакет ANSYS использует четыре двумерных конечных элемента (PLANE42, PLANE82, PLANE182 и PLANE183) и семь трехмерных (SOLID45, SOLID185, SOLID95, SOLID186, SOLID285, SOLID92 и SOLID187). И это без учета элементов, использующих иные (в дополнение к механическим) степени свободы, специальные свойства материалов и т.п. (например, coupled-field-элементы, такие как PLANE13, PLANE223, SOLID5, SOLID98, SOLID226, SOLID227 и др.).
Соответственно возникает вопрос — насколько густой должна быть сетка и как выбрать тип элемента, для того чтобы свести к минимуму погрешность полученного результата и при этом решить задачу за разумное время.
В настоящей работе на примере программного комплекса ANSYS исследован вопрос влияния порядка аппроксимации, числа точек интегрирования, учета дополнительных форм перемещений, а также густоты сетки на точность решения тех задач механики, в которых изгибные напряжения имеют существенное значение. Исследование проведено для 11 указанных элементов. Кроме того, продемонстрированы вызывающие погрешность эффекты shear locking и hourglassing.
Подобные исследования для двумерных и трехмерных конечных элементов в виде прямоугольников и прямоугольных параллелепипедов проводились в работах [1-3] на примере задачи об изгибе балки в двумерной и трехмерной постановках. Однако при решении этой
задачи не были рассмотрены часто применяемые элементы в виде треугольников и тетраэдров. В настоящей работе сделана попытка восполнить этот пробел.
В качестве примера рассмотрим задачу по определению напряженно-деформированного состояния стержня c заделкой (рис. 1), нагруженного поперечной силой (т.е. рассмотрим пример, подобный примеру VM16 из ANSYS Verification Manual, в котором применен элемент PLANE42).
Этот пример может быть полезен при решении задач изгиба как стержневых, так и тонкостенных конструкций. Тем более, что в современной практике часто вместо специальных балочных элементов используют трехмерные элементы, неправильное использование которых может привести к ошибкам.
Размеры стержня: длина l = 150 мм; высота сечения d = 5 мм; ширина b = 2,5 мм. Модуль упругости материала E = 70 ГПа, коэффициент Пуассона v = 0. Стержень нагружен поперечной силой P = 5 Н. Согласно формулам сопротивления материалов (используется гипотеза плоских сечений Эйлера-Бернулли) перемещение неза-
Pl3 r bd3
крепленного конца стержня равно usm = тг^т, где I = tjtt, при этом
3EI 12
имеем usm = 3,085714 « 3,086 мм.
В теории упругости [4] рассмотрено два вида закрепления стержня, удовлетворяющих условию горизонтальности начальной касательной и условию вертикальности опорного сечения. В первом случае уравнение изогнутой оси стержня и, соответственно, значение перемещения на незакрепленном конце стержня совпадает с результатом, полученным в курсе сопротивления материалов, т.е. usm. Во втором случае
Pl3 3Pl
перемещение на свободном конце равно uel = —— + —:—, где второе
3EI 2G d
слагаемое учитывает влияние сдвигов на прогиб (G — модуль сдвига). Для модельной задачи uel = usm + 0,0000064 = 3,085721 мм« usm.
В работе использованы две расчетные схемы — плоское напряженное состояние с явным заданием толщины (plane stress with thickness) (см. рис. 1) и трехмерная модель (рис. 2).
Рис. 1. Расчетная схема — плоское напряженное состояние с явным заданием толщины (plane stress with thickness)
Рис. 2. Расчетная схема — трехмерная модель
Проведено сравнение полученного с помощью А^УБ перемещения на незакрепленном конце стержня п* с решением по формулам сопротивления материалов и подсчитано отношение п* /п^м • Вычисления проводились на четырех сетках конечных элементов: 1 х 6, 2 х 12, 4 х 12 и 8 х 24, где первая цифра — число элементов по толщине, а вторая — число элементов по длине балки.
Как правило, упомянутые элементы А^УБ в виде прямоугольника или параллелепипеда являются стандартными сирендиповыми 1-го или 2-го порядка элементами, а в виде треугольника или тетраэдра используют стандартные Ь-координаты [5]. Однако в некоторых случаях функции формы могут иметь дополнительные слагаемые и тогда ссылка на литературный источник будет приведена в явном виде. В любом случае функции формы всех элементов как стандартных, так и улучшенных, приведены в работе [6]. Там же приведены схемы расположения точек интегрирования.
Перечень рассмотренных элементов и их значимых опций приведен ниже. Отметим, что часть элементов треугольной/тетраэдрической формы являют собой вырожденные варианты прямоугольников/параллелепипедов со слившимися узлами. Так, плоский 4-узловой элемент вырождается в 3-узловой треугольный, а плоский 8-узловой — в 6-узловой треугольный; 8-узловой параллелепипед вырождается в 4-узловой тетраэдр, а 20-узловой — в 10-узловой тетраэдр.
• PLANE42: плоский 4-узловой линейный элемент, четыре точки интегрирования по схеме 2 х 2.
о Опция K2 = 0: Include extra displacement shapes (используются дополнительные функции формы согласно [7] для учета линейного распределения деформаций в элементе).
о Опция K2 = 1: Exclude extra displacement shapes (сирендиповый элемент 1-го порядка).
о Вырожденный элемент (3-узловой треугольный), одна точка интегрирования (три точки для осесимметричного применения элемента).
• PLANE182: плоский 4-узловой сирендиповый элемент 1-го порядка.
о Опция K1=0: Full integration with B-bar method, четыре точки интегрирования 2 x 2.
о Опция K1 = 1: Uniform reduced integration with hourglass control, одна точка интегрирования.
о Опция K1 = 2: Enhanced strain formulation, четыре точки интегрирования 2 x 2 (вводит пять внутренних степеней свободы для предотвращения эффектов shear locking и volumetric locking).
о Опция K1 = 3: Simplified enhanced strain formulation, четыре точки интегрирования 2 x 2 (вводит четыре внутренних степени свободы для предотвращения эффекта shear locking).
о Вырожденный элемент (3-узловой треугольный), одна точка интегрирования.
• PLANE82: плоский 8-узловой сирендиповый элемент 2-го порядка, четыре точки интегрирования — редуцированная схема 2 x 2.
о Вырожденный элемент (6-узловой треугольный), три точки интегрирования.
• PLANE183: плоский 8-узловой сирендиповый элемент 2-го порядка, четыре точки интегрирования — редуцированная схема 2 x 2.
о Вырожденный элемент (6-узловой треугольный), три точки интегрирования.
• SOLID45: 3D- 8-узловой линейный элемент.
о Опция K1 =0: Include extra displacement shapes shapes (используются дополнительные функции формы согласно [7] для учета линейного распределения деформаций в элементе).
о Опция K1 = 1: Exclude extra displacement shapes. Сирендиповый элемент 1-го порядка.
о Опция K2 = 0: Full integration, восемь точек интегрирования 2 x x 2 x 2.
о Опция K2 = 1: Uniform reduced integration with hourglass control, одна точка интегрирования (K2,= 1 одновременно подавляет extra displacement shapes, т.е. автоматически задает K1 = 1).
о Вырожденный элемент (4-узловой тетраэдр, K1 = 1 устанавливается автоматически).
• SOLID185: 3D- 8-узловой сирендиповый элемент 1-го порядка. о Опция K2 = 0: Full integration with B-bar method, восемь точек
интегрирования 2 x 2 x 2.
о Опция K2 = 1: Uniform reduced integration with hourglass control, одна точка интегрирования.
о Опция K2 = 2: Enhanced strain formulation, восемь точек интегрирования 2 x 2 x 2 (вводит 13 внутренних степеней свободы для предотвращения эффектов shear locking и volumetric locking).
о Опция K2 = 3: Simplified enhanced strain formulation, восемь точек интегрирования 2 x 2 x 2 (вводит девять внутренних степеней свободы для предотвращения эффекта shear locking).
о Вырожденный элемент (4-узловой тетраэдр).
• SOLID95: 3D- 20-узловой сирендиповый элемент 2-го порядка. о Опция K11 = 0: 14 Point Rule, 14 точек интегрирования.
о Опция K11 = 1: Reduced integration option, восемь точек интегрирования 2 x 2 x 2.
о Вырожденный элемент (10-узловой тетраэдр), четыре точки интегрирования.
• SOLID186: 3D- 20-узловой сирендиповый элемент 2-го порядка. о Опция K2 = 0: Reduced integration option, восемь точек интегрирования 2 x 2 x 2.
о Опция K2 = 1: Full integration, 14 точек интегрирования. Вырожденный элемент (10-узловой тетраэдр), четыре точки интегрирования.
• SOLID285: 3D- 4-узловой линейный тетраэдр, четыре точки интегрирования.
• SOLID92: 3D- 10-узловой тетраэдр 2-го порядка, четыре точки интегрирования.
• SOLID187: 3D- 10-узловой тетраэдр 2-го порядка, четыре точки интегрирования.
В результате проведенного анализа выявлены источники погрешностей, что позволяет сформулировать следующие выводы.
1. Стандартные линейные элементы в виде прямоугольников (PLANE42 c K2-option = 1 и PLANE182 с K1-option = 0) и прямоугольных параллелепипедов (SOLID45 c K1-option= 1 и K2-option = 0; и SOLID185 с K2-option = 0) при максимальном числе точек интегрирования (full integration) — четыре точки для прямоугольников (схема 2 x 2) и восемь для параллелепипедов (схема 2 x 2 x 2) дают большую погрешность: полученное значение перемещения меньше точного на
44% даже на сравнительно густой сетке с восемью элементами по толщине. Причина состоит в эффекте сдвигового запирания (shear locking), заключающегося в том, что такой элемент 1-го порядка аппроксимации имеет повышенную изгибную жесткость - в точках интегрирования возникают паразитные сдвиговые напряжения, вызванные тем, что стороны элемента не могут изгибаться. Следовательно, если в конструкции имеются значительные напряжения изгиба, то применение таких элементов может привести к ошибкам даже на густой сетке. Для линейных треугольников (вырожденные элементы PLANE42 и PLANE182) и тетраэдров (SOLID285 и вырожденные элементы SOLID45 и SOLID185) наблюдается аналогичный результат: на сетке с восемью элементами по толщине полученное значение перемещения меньше точного для треугольников на 73 %, а для тетраэдров на 68 %.
2. Для исключения эффекта shear locking уменьшают число точек интегрирования, т.е. применяют схему reduced integration. В этом случае линейные прямоугольники (PLANE182 с K1-option=1) и прямоугольные параллелепипеды (SOLID45 c K1-option=1 и K2-option=1 и SOLID185 с K2-option= 1) имеют одну точку интегрирования и при малом (до четырех) числе элементов по толщине страдают от эффекта hourglassing (моды с нулевой энергией). Он выражается в том, что в элементе возникают паразитные деформации, при которых компоненты напряжений в единственной точке интегрирования равны нулю. Соответственно, такие деформации не влияют на энергию системы и поэтому могут привести к ошибочным решениям. В случае, когда по толщине стержня использован один элемент, решение физически неправдоподобно — отношение u*/usm порядка нескольких тысяч, так как точки интегрирования находятся на нейтральной линии стержня и поэтому модель не отражает сопротивление изгибу. В целом использование таких элементов дает завышенное значение перемещения, так как при этом уменьшается изгибная жесткость модели. Для борьбы с эффектом hourglassing в ANSYS используется искусственная жесткость, которая "может неблагоприятно повлиять на точность решения" [6].
3. Прямоугольные параллелепипеды второго порядка (SOLID95 с K11-option= 1 и SOLID186 с K2-option = 0) при одном элементе по толщине в случае уменьшенного числа точек интегрирования (восемь точек по схеме 2 x 2 x 2 тоже подвержены влиянию эффекта hourglassing, сопровождаемого появлением отрицательного значения диагонального элемента матрицы жесткости и бесконечных перемещений. В случае одного элемента по толщине такой эффект может наблюдаться и для прямоугольного элемента PLANE82 [6].
В результате анализа причин погрешностей расчета для практического применения при решении стационарных задач механики, в которых напряжения изгиба имеют существенное значение, можно дать следующие рекомендации.
1. При применении линейных элементов в виде прямоугольников и прямоугольных параллелепипедов с максимальным числом точек интегрирования (full integration) необходимо применять специальные опции include extra displacement shapes и enhanced strain formulation. Эти опции подавляют эффект shear locking и на модельной задаче позволяют получить погрешность в 0,5 % уже при одном элементе по толщине. В противном случае, так же как и при применении линейных треугольников и тетраэдров, можно ожидать значительного занижения полученных перемещений по сравнению с точным решением.
2. При использовании линейных элементов в виде прямоугольников и прямоугольных параллелепипедов с одной точкой интегрирования (reduced integration) необходимо использовать не менее четырех элементов по толщине конструкции; только в этом случае погрешность составляет не более 1,5 %.
3. В случае применения квадратичных элементов с уменьшенным числом точек интегрирования (reduced integration) в виде прямоугольных параллелепипедов необходимо использовать не менее двух элементов по толщине конструкции; в этом случае погрешность решения модельной задачи не превышает 0,1%. В противном случае можно получить физически неправдоподобные решения, а именно, бесконечные перемещения. Хотя на исследованной модельной задаче такие квадратичные элементы в виде плоских прямоугольников обеспечили погрешность в 0,03 %, рекомендуется использовать не менее двух элементов по толщине конструкции, поскольку плоские квадратичные элементы имеют редуцированную схему интегрирования 2 х 2 и, следовательно, подвержены эффекту hourglassing. Дополнительное преимущество элементов с уменьшенным числом точек интегрирования — это предотвращение эффекта объемного запирания (volumetric locking), возникающего при расчете конструкций из несжимаемых или практически несжимаемых материалов.
4. Если используются элементы второго порядка с максимальным числом точек интегрирования (full integration) в форме треугольников, прямоугольных параллелепипедов и тетраэдров, то они позволяют получать точное решение даже на грубой сетке (погрешность с одним элементом по толщине 0,6%, с двумя и более слоями элементов по толщине — 0,03 %), так как позволяют моделировать линейное распределение деформаций по толщине элемента.
СПИСОК ЛИТЕРАТУРЫ
1. ABAQUS User Manual. Getting Started with ABAQUS. Version 6.6.
2. Eric Qiuli Sun. Shear Locking and Hourglassing in MSC Nastran, ABAQUS, and ANSYS. http://www.mscsoftware.com/events/vpd2006/na/presenta-tions/tech_papers/27.pdf
3. КрасновскийЕ. Е. Точность решения плоских задач механики в ANSYS при наличии больших изгибных напряжений // Сб. тр. Седьмой конф. пользователей програм. обеспеч. CAD-FEM GmbH (Москва, 25-26 апреля 2007 г.) / Под ред. А.С. Шадского. - М.: Полигон-пресс, 2007. - С. 230-232.
4. Зенкевич О. Метод конечных элементов в технике. - М.: Мир, 1975. -541 с.
5. T a y l o r R. L., B e r e s f o r d P. J., and Wilson E. L. A Nonconforming element for stress analysis, international // J. for Numerical Methods in Engineering.
- 1976. - Vol. 10. - P. 1211-1219.
6. A N S Y S Element library.
7. Филоненко-БородичМ. М. Теория упругости. - М.: Физматлит, 1959.
- 365 с.
Статья поступила в редакцию 25.10.2011
