Использование препроцессора Femap совместно с комбинированным методом решения нелинейных конечноэлементных задач Текст научной статьи по специальности «Физика»

УДК 539.3; 539.4

DOI: 10.30987/article_5cb58f52b9a0d7.77781783

А.П. Бабин, М.В. Зернин, А.А. Просолович

ИСПОЛЬЗОВАНИЕ ПРЕПРОЦЕССОРА FEMAP СОВМЕСТНО С КОМБИНИРОВАННЫМ МЕТОДОМ РЕШЕНИЯ НЕЛИНЕЙНЫХ КОНЕЧНОЭЛЕМЕНТНЫХ ЗАДАЧ

Ранее в БГТУ разработан универсальный, но сравнительно сложный комбинированный алгоритм решения конечноэлементных задач, в которых разрешающие уравнения могут учитывать нелинейности различного типа. Предпринята попытка использовать современный интерфейс программного пакета FEMAP совместно с реализующей этот алгоритм программой.

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

Ключевые слова: метод конечных элементов, нелинейные разрешающие уравнения, итерационный алгоритм, интерфейс подготовки данных, конвертирование данных.

A.P. Babin, M.V. Zernin, A.A. Prosolovich

FEMAP PREPROCESSOR USE JOINTLY WITH COMBINED METHOD FOR SOLUTION OF NON-LINEAR FINITE ELEMENT PROBLEMS

Earlier in BSTU there was developed a universal but comparatively complex combined algorithm for the solution of finite element problems in which resolving equations can take into account different nonlineari-ty. An attempt is made to use a modern interface of the FEMAP program package jointly with the program realizing this algorithm. In fact a converter of data prepared in FEMAP formats required for the original

Постановка задачи

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

program is developed. A new program formed is tested by some examples. Now the program realizing a combined algorithm has a modern interface and it is convenient for use.

Key words: finite element method, nonlinear resolving equations, iteration algorithm, interface of data preparation, data translation.

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

Для обеспечения сходимости в таких случаях авторами разработан комбинированный алгоритм [1-5], который можно рассматривать как объединение существующих алгоритмов. Этот алгоритм реализует возможность применения в каждый конкретный момент вычисления той итерационной процедуры, которая гарантиру-

ет сходимость. Практически все конечно-элементные программные пакеты содержат следующие составляющие: препроцессор (ввод и подготовка данных); программа -решатель задачи; постпроцессор (анализ и представление в удобном виде результатов расчетов). Ранее разработанная программа [1-5] реализует комбинированный алгоритм, но стадия подготовки данных для расчетов недостаточно удобна. В настоящее время выполнена разработка алгоритма, который несет функцию подготовки данных об исследуемой конечноэлемент-ной модели, созданной в программном пакете FEMAP with NX Nastran [6]. Факти-

чески решены три задачи. На первом этапе в среде FEMAP with NX Nastran производится моделирование геометрии и конеч-ноэлементной сетки, приложение нагрузок и граничных условий, а также вывод информации о модели в отдельный документ. На втором этапе происходит поиск рационального метода обработки полученных данных и его реализация на языке С++. На третьем этапе изменена часть кода разработанной ранее программы-решателя с целью обеспечения возможности ее взаимодействия с информацией, получаемой из программного комплекса FEMAP with NX Nastran.

Краткое описание алгоритма комбиниро

В задачах механики деформируемого твердого тела может проявляться физическая или геометрическая нелинейность. В зависимости от уровней воздействия на материалы для КЭ различных участков тел могут применяться модели теории упругости (в том числе нелинейные), теории пластичности или теории ползучести. Методы решения физически и геометрически нелинейных задач достаточно хорошо разработаны в соответствующих разделах механики деформируемого твердого тела. Авторами реализованы [2; 3] различные известные итерационные методы решения нелинейных систем уравнений как сходящейся последовательности линейных задач. В качестве основных алгоритмов решения нелинейных задач рассматривались применяемые в теории пластичности метод дополнительных напряжений (МДН) и метод дополнительных деформаций (МДД). Практика расчетов показала, что ни МДН, ни МДД в отдельности не позволяют решать задачи со всеми встречающимися типами нелинейностей.

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

1нного метода

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

Здесь {а} - вектор напряжений; [о] -матрица упругих коэффициентов; {£■} -вектор деформаций; {а*} и {£■*} - соответственно векторы дополнительных напряжений и деформаций, в которые сведены все нелинейные составляющие физического уравнения. Переопределяя такие дополнительные (псевдоначальные) параметры на каждой итерации, можно решить нелинейную задачу как рекуррентную последовательность линейных задач. Уравнения (1) аналогичны уравнениям упругой задачи при наличии одновременно начальных напряжений и начальных деформаций. Для удовлетворения физических зависимостей на каждой итерации определяются новые значения или дополнительных напряжений, или дополнительных деформаций.

Рис. 1 демонстрирует, как осуществляется выбор варианта итерационного процесса в соответствии с комбинированным методом (МДДН). Основную идею этого алгоритма можно сформулировать следующим образом. Если точка упругого решения на текущей итерации находится

над кривой с(б) (точка 1), то следующее решение отыскивается по схеме МДН, а если ниже (точка 2), то по схеме МДД.

Конечноэлементная система разрешающих уравнений имеет вид

[К]{и} = / [в] № *}& -1 [в] {с *}&

где [К] = Дв]7 [о]{В>Ж - матрица жестко-

V

сти; {и} - вектор узловых перемещений; {Р> - вектор узловых сил от внешних

воздействий; - вектор

V

узловых сил от дополнительных деформаций;

- вектор узловых сил от

V

дополнительных напряжений.

с

- б

Рис. 1. К выбору варианта реализации итерационного процесса в методе дополнительных деформаций и напряжений

Опишем итерационную схему комбинированного алгоритма по рис. 1 . На первой итерации вычисляются упругие напряжения и деформации от приложенной внешней нагрузки. При этом начальные деформации {б *} и начальные напряжения {с *} принимаются равными нулю. Допустим, что полученная точка находится выше кривой с(б) (точка 1 на рис. 2).

Тогда определяем напряжения с *(1), которые соответствуют рассчитанным деформациям и нелинейной характеристике с(б) (точка 1'), и вычисляем приращения дополнительных напряжений:

ёа*(1) = а*(1)-а(1). Начальные напряжения следующей итерации {с *} изменяются на величину приращений дополнительных напряжений ёс *(1). На следующей итерации находим деформации и напряжения (точка 2). Таким образом на первой итера-

ции реализуется схема МДН.

Точка 2 находится ниже кривой с(б). Следовательно, б *(2) - те деформации, которые соответствуют полученным напряжениям и нелинейной характеристике с(б). Далее определяются приращения

дополнительных деформаций ёб *(2) и начальные деформации {б *} для следующей итерации. Таким образом реализуется схема МДД. Итерационный процесс заканчивается, когда приращения дополнительных деформаций и напряжений становятся меньше заранее заданной допустимой погрешности этих параметров. В итоге процесс решения в каждом нелинейном КЭ может сходиться либо по схеме МДД, либо по схеме МДН, либо по их комбинации (каждый из этих методов используется на некоторых этапах решения), независимо от схемы сходимости в других КЭ.

В некоторых случаях на диаграммах

свойств КЭ проявляется скачкообразное изменение этих свойств. Рассмотрим в качестве примера нелинейые свойства слоев в зоне контактирования

поверхностей (так называемого «третьего тела», или «контактной псевдосреды» [15]). В общем случае нормальные деформации £п контактной псевдосреды нелинейно зависят от нормальных напряжений оп в контакте (контактных давлений). При сухом контактировании эта взаимосвязь качественно отображается кривой линией, показанной на рис. 2а. Если после нагружения определенными давлениями приложить к контактной псевдосреде касательные усилия т и измерить ее касательные деформации у, то получим нелинейные зависимости, изображенные на рис. 2б. Сначала контактный слой нелинейно

деформируется в пределах

предварительного смещения. После достижения некоторого предельного значения касательного напряжения т™ах

начинается относительное движение поверхностей при касательном

напряжении т^ < т™* (сопротивление

трению скольжения меньше

сопротивления трению покоя). Поэтому на рис. 2б имеется зона скачкообразного

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

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

(растягивающих) нормальных напряжений обрывается после достижения некоторого их предельного значения ст^3*. Такие

нелинейные свойства контактной псевдосреды можно задавать в МКЭ как нелинейные свойства контактных КЭ [1-5].

а)

У

б)

Рис . 2. Схематическое изображение графиков взаимосвязи напряжений и соответствующих перемещений контактного слоя в нормальном (а) и касательном (б) направлениях

Применяемая в сочетании с МКЭ дискретизация временной оси позволяет рассматривать практически любой харак-

тер изменения внешних воздействий на объект во времени. Высокопроизводительные ЭВМ позволяют прикладывать внеш-

ние нагрузки по малым шагам и отыскивать решение на каждом шаге. Таким образом, фактически моделируется реальный процесс нагружения с учетом истории изменения внешнего воздействия. Можно моделировать процессы разгрузки, появления остаточных напряжений и деформаций, последующие этапы повторного на-гружения [1-5].

H

Рис.3. Свойства контактного конечного элемента в месте зазора поверхностей

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

ностью) к задаче с внутренней нелинейностью, определяемой свойствами контактных КЭ.

В итоге построен итерационный алгоритм [1-5] для решения конечноэлементных задач, в которых разрешающие уравнения могут иметь нелинейности различного типа. Этот алгоритм можно рассматривать как объединение существующих алгоритмов: он реализует возможность применения при расчете того из них, который гарантирует сходимость итерационной процедуры. В одной и той же расчетной схеме, но в разных КЭ с различными нелинейностями одновременно могут использоваться различные варианты итерационных процедур. Кроме того, даже в одном и том же КЭ могут использоваться различные варианты итерационных процедур, если нелинейные характеристики для разных компонент НДС соответствуют нелинейностям различного типа.

Алгоритм позволяет решить многие задачи с нелинейностями различных типов, в том числе и с несколькими типами нелинейностей. Но встроить этот сложный алгоритм в известные промышленные программные пакеты, такие как Nastran и Ansys, очень сложно. Нами предпринята попытка использовать возможности препроцессора FEMAP для подготовки исходных данных для оригинальной программы-решателя [1-5]. В принципе пакет FEMAP является независимым программным продуктом и взаимодействует с пакетами МКЭ Nastran, Ansys и другими, передавая им подготовленные файлы данных и получая файлы результатов для дальнейшего анализа. Нами преобразованы форматы подготовленных данных и оригинальная программа таким образом, чтобы можно было решать нелинейные задачи с различными типами нелинейностей.

Разработан алгоритм обработки информации о созданной в программном комплексе FEMAP with NX Nastran модели. В алгоритме учтены все особенности представления данных,

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

задач, отвечающая за создание конечноэлементной модели исследуемого объекта. Изменения несут в себе функцию взаимодействия этой программы с обработанными данными о модели, созданной в программном комплексе FEMAP with NX Nastran.

Результаты решения некоторых нелинейных

Рис. 4. Модель контактирования двух деталей через нелинейный шероховатый слой

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

В качестве простого примера нелинейной задачи была решена задача о многократном нагружении шероховатого слоя и поиске площадки контакта с предположением, что центральная часть контактного слоя предварительно была деформирована жестким штампом [7] и в ней возникли остаточные деформации (рис. 4). Зазор, образованный первоначальным воздействием штампа, в модели задается как свойство контактного конечного элемента (по рис. 3), поэтому контактная среда смоделирована сплошным слоем элементов с соответственно различными свойствами. Информация о свойствах материалов и КЭ задается непосредственно в программе-решателе.

задач

Эта же задача была решена по ранее разработанной программе [1-5]. Результаты решения в зависимости от величины нагрузки F приведены на рис. 5. Здесь видно, что при приложении малой нагрузки нормальные напряжения в центре равны нулю (нагрузку эти элементы не воспринимают) (рис. 5а). При увеличении нагрузки в контакт входят не только крайние, не деформированные ранее элементы, но и элементы, находящиеся в центре (рис. 5б). При еще большем увеличении нагрузки элементы на краях и в центре деформируются

упругопластически (рис. 5в). При увеличении нагрузки упругопластическая зона растет, а зона упругого деформирования сокращается, затем вся зона в центре начинает деформироваться упругопластически. В итоге распределение напряжений в контакте становится плавным, без скачков и разрывов (рис. 5г).

Опишем подробнее, как

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

необходимо задать, какая конкретно нужна информация. В той же вкладке List выбирается пункт Model, поочередно из списка нужно выбрать пункты Node, Element, Constrain - Individual, Load -Individual. Таким образом, выводится информация поочередно об узлах, элементах, условиях закрепления и

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

'/.'// //// //////////////////////////

б

Рис. 6. Модель контактирования кольца и основания через нелинейный шероховатый слой

программа-решатель обрабатывает

входные данные.

4 6

Длина, мм

6)

Длина, мм

г)

Была решена задача о поиске площадки контакта пластины и кольца с расположенным между ними тонким слоем с нелинейными свойствами. Характер свойств контактной среды

продемонстрирован на рис. 2. Геометрия модели представлена на рис. 6. По этим данным в среде FEMAP with NX Nastran была построена конечноэлементная модель. Задача решена для различных внешних нагрузок, а именно: 12, 24, 48, 96, 192 Н.

В результате работы программы были получены данные о перемещении всех узлов конечных элементов исследуемой модели на каждой итерации. По этим данным в программном комплексе Mathcad 15.0 были построены диаграммы деформирования контактного слоя для всех выбранных уровней нагрузки.

На рис. 7 приведена диаграмма состояния контактного слоя в направлении нормали (направление Y) без приложения

Рис. 5. Графики давлений в шероховатом слое при различных уровнях его нагружения

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

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

№ соответственных узлов контакта Рис. 7. Диаграмма состояния контактного слоя до приложения нагрузки: I - линия, очерчивающая контактирующую поверхность кольца; II - линия поверхности пластины; III - линия, соответствующая верхней поверхности слоя

Далее приведена серия диаграмм деформирования контактных конечных элементов в зависимости от уровня приложенных внешних нагрузок: 12 Н (рис. 8а), 24 Н (рис. 8б), 48 Н (рис. 8в), 96 Н (рис. 8г), Здесь показаны только половины симметричных рисунков. Картина деформирования при

максимальной нагруззке (192 Н) показана

на рис. 9. На рис. 8 и 9 видно, как при увеличении нагрузки ширина площадки контакта увеличивается.

По данным расчета построена диаграмма, отражающая зависимость ширины площадки контакта в проекции на ось Y от уровня внешних нагрузок (рис. 10).

t» к

я

£

=

S:

с с.

ыо

1x10

: 4 б

№ соответственных узлов контакта

а)

- :xio

0 с

= 1*10"

1 -

с

м

2 4 6

№ соответственных узлов контакта б)

2*10

L> О О

- 1x10"

5 !r К

ь

О

м

3 2^10 *

и

о О

5 1x10"

й

ё X К

&

О С

м

: 4 б

Л'1 соответственных узлов контакта

в) '

: 4 б

№ соответственных узлов контакта

г)

Рис. 8. Деформация контактных конечных элементов при нагрузках: а -12 Н; б - 24 Н; в - 48 Н; г - 96 Н

- >10' к

о С

о

- 1x10"

й &

=

к с.

О

с М

-4

2 4 6

№ соответственных узлов контакта

Рис. 9. Деформация контактных конечных элементов при нагрузке 192 Н

Рис. 10. Зависимость ширины площадки контакта от величины нагрузки

Заключение

Разработан алгоритм обработки информации о созданной в программном комплексе FEMAP with NX Nastran модели. Алгоритм реализован на ЭВМ в виде программы. Программа была протестирована, результаты работы согласуются с ожидаемыми. Изменена часть кода программы решения нелинейных задач, отвечающая за

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

1. Бабин, А.П. Конечноэлементный алгоритм решения контактных задач с учетом нелинейных эффектов / А.П. Бабин // Динамика, прочность и надежность транспортных машин. - Брянск: БГТУ, 2002. - С. 138-148.

2. Зернин, М.В. К исследованию контактной жесткости с использованием модели механики контактной псевдосреды / М.В. Зернин, А.П. Бабин // Заводская лаборатория. - 2001. - Т. 67. -№ 6. - С. 51-54.

создание конечноэлементной модели исследуемого объекта. Программа была протестирована. Ход расчета

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

3. Бабин, А.П. Учет влияния нелинейных свойств поверхностных слоев при конечноэлементном решении задач о контактном взаимодействии деформируемых тел / А.П. Бабин, М.В. Зернин // Трение и смазка в машинах и механизмах. -2008. - № 3. - С. 3-16.

4. Бабин, А.П. Конечноэлементное моделирование контактного взаимодействия с использованием положений механики контактной псевдосреды / А.П. Бабин, М.В. Зернин // Изв. РАН. Механика твердого тела. - 2009. - № 4. - С. 84-107.

5. Зернин, М.В. Алгоритм комбинированного метода решения конечноэлементных задач с не-линейностями различного типа / М.В. Зернин, А.П. Бабин // Вестник Брянского государственного технического университета. - 2009. - № 4. -С. 57-64.

1. Babin, A.P. Finite element algorithm for contact problem solution taking into account nonlinear effects / A.P. Babin // Dynamics, Strength and Reliability of Transport Machines. - Bryansk: BSTU, 2002. - pp. 138-148.

2. Zernin, M.V. To contact stiffness investigation using models of contact pseudo-environment mechanics / M.V. Zernin, A.P. Babin // Factory Laboratory. 2001. - Vol.67. - No.6. - pp. 51-54.

3. Babin, A.P. Account of surface layer nonlinear properties impact at finite element solution of problems on contact interaction of bodies deformable / A.P. Babin, M.V. Zernin // Friction and Lubrication in Machines and Mechanisms. - 2008. - No.3. - pp. 3-16.

4. Babin, A.P. Finite element simulation of contact

6. Рудаков, К.Н. UGS Femap 9.3. Геометрическое и конечноэлементное моделирование конструкций / К.Н. Рудаков. - М.: ДМК Пресс, 2009. - 296 с.

7. Пономарев, В.С. Внедрение жесткого штампа в эластичный материал. Механика деформированного твердого тела / В.С. Пономарев. -Томск, 1987. - С. 147-152.

interaction using positions of contact pseudoenvironment mechanics / A.P. Babin, M.V. Zernin // Proceedings of RAS. Mechanics of Rigid Body. -2009. - No.4. - pp. 84-107.

5. Zernin, M.V. Algorithm of combined method for finite element problem solutions with different non-linearity / M.V. Zernin, A.P. Babin // Bulletin of Bryansk State Technical University. - 2009. - No.4. - pp. 57-64.

6. Rudakov, K.N. UGS Femap 9.3. Geometrical and Finite Element Simulation of Constructions / K.N. Rudakov. - M.: MDK Press, 2009. - pp. 296.

7. Ponomaryov, V.S. Rigid Die Introduction into Elastic Material. Deformed Solid Mechanics / V.S. Ponomaryov. - Tomsk, 1987. - pp. 147-152.

Статья поступила в редакцию 10.01.19 Рецензент: к.т.н., доцент Брянского государственного

технического университета Болдырев А.П. Статья принята к публикации 22. 03. 19.

Сведения об авторах:

Бабин Александр Павлович, к.т.н., доцент кафедры «Информатика и программное обеспечение» Брянского государственного технического университета, e-mail: apbabin@yandex. ru. Зернин Михаил Викторович, к.т.н., доцент кафедры «Механика и ДПМ» Брянского государст-

Babin Alexander Pavlovich, Can. Sc. Tech., Assistant Prof. of the Dep. "Informatics and Software", Bryansk State Technical University, e-mail:

apbabin@yandex. ru.

Zernin Michail Victorovich, Can. Sc. Tech., Assistant Prof. of the Dep. "Mechanics and DPM", Bryansk State Technical University, e-mail: zerninmv@mail. ru.

венного технического университета, e-mail: zerninmv@mail. ru.

Просолович Алексей Александрович, магистрант кафедры «Механика и ДПМ» Брянского государственного технического университета, e-mail: palpatinec@mail. ru.

Prosolovich Alexey Alexandrovich, Master degree student of the Dep. "Mechanics and DPM", Bryansk State Technical University, e-mail:

palpatinec@mail.ru.