Bryansk, Bryansk, Russian Federation
Bryansk, Bryansk, Russian Federation
Bryansk, Bryansk, Russian Federation
GRNTI 55.01 Общие вопросы машиностроения
GRNTI 55.13 Технология машиностроения
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 nonlinearity. An attempt is made to use a modern interface of the FEMAP program package jointly with the program realizing this algorithm. In fact a convertor 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.
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 |
|
1'12' ’ |
|
2 |
|
2' |
|
3 |
|
3' |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
где 
- матрица жесткости; 
- вектор узловых перемещений; 
- вектор узловых сил от внешних воздействий; 
- вектор узловых сил от дополнительных деформаций; 
- вектор узловых сил от дополнительных напряжений.
Опишем итерационную схему комбинированного алгоритма по рис. 1. На первой итерации вычисляются упругие напряжения и деформации от приложенной внешней нагрузки. При этом начальные деформации и начальные напряжения принимаются равными нулю. Допустим, что полученная точка находится выше кривой (точка 1 на рис. 2). Тогда определяем напряжения , которые соответствуют рассчитанным деформациям и нелинейной характеристике (точка 1'), и вычисляем приращения дополнительных напряжений: 
=-. Начальные напряжения следующей итерации изменяются на величину приращений дополнительных напряжений . На следующей итерации находим деформации и напряжения (точка 2). Таким образом на первой итерации реализуется схема МДН.
Точка 2 находится ниже кривой . Следовательно, 
- те деформации, которые соответствуют полученным напряжениям и нелинейной характеристике . Далее определяются приращения дополнительных деформаций 
и начальные деформации для следующей итерации. Таким образом реализуется схема МДД. Итерационный процесс заканчивается, когда приращения дополнительных деформаций и напряжений становятся меньше заранее заданной допустимой погрешности этих параметров. В итоге процесс решения в каждом нелинейном КЭ может сходиться либо по схеме МДД, либо по схеме МДН, либо по их комбинации (каждый из этих методов используется на некоторых этапах решения), независимо от схемы сходимости в других КЭ.
В некоторых случаях на диаграммах свойств КЭ проявляется скачкообразное изменение этих свойств. Рассмотрим в качестве примера нелинейые свойства слоев в зоне контактирования поверхностей (так называемого «третьего тела», или «контактной псевдосреды» [1-5]). В общем случае нормальные деформации 
контактной псевдосреды нелинейно зависят от нормальных напряжений в контакте (контактных давлений). При сухом контактировании эта взаимосвязь качественно отображается кривой линией, показанной на рис. 2а. Если после нагружения определенными давлениями приложить к контактной псевдосреде касательные усилия и измерить ее касательные деформации , то получим нелинейные зависимости, изображенные на рис. 2б. Сначала контактный слой нелинейно деформируется в пределах предварительного смещения. После достижения некоторого предельного значения касательного напряжения начинается относительное движение поверхностей при касательном напряжении (сопротивление трению скольжения меньше сопротивления трению покоя). Поэтому на рис. 2б имеется зона скачкообразного изменения свойств (вертикальные участки линий 
), соответствующая переходу от состояния покоя в условиях предварительного смещения к процессу относительного скольжения поверхностей, которому соответствуют горизонтальные участки линий .
Если моделируются односторонние связи контактирующих поверхностей, то координаты можно отсчитывать только в сторону сжимающих нормальных напряжений и деформаций (рис. 2а). Силами сцепления поверхностей (например адгезионными) обычно пренебрегают из-за их малости по сравнению со сжимающими контактными давлениями. В более общем случае учитывают для контактной псевдосреды различные модели взаимодействия поверхностей при их прижатии друг к другу и при приложении растягивающих нагрузок. В этом случае координаты отсчитываются в двух направлениях, причем график положительных (растягивающих) нормальных напряжений обрывается после достижения некоторого их предельного значения 
. Такие нелинейные свойства контактной псевдосреды можно задавать в МКЭ как нелинейные свойства контактных КЭ [1-5].
Применяемая в сочетании с МКЭ дискретизация временной оси позволяет рассматривать практически любой характер изменения внешних воздействий на объект во времени. Высокопроизводительные ЭВМ позволяют прикладывать внешние нагрузки по малым шагам и отыскивать решение на каждом шаге. Таким образом, фактически моделируется реальный процесс нагружения с учетом истории изменения внешнего воздействия. Можно моделировать процессы разгрузки, появления остаточных напряжений и деформаций, последующие этапы повторного нагружения [1-5].
|
|
|
|
|
|
|
h |
|
H |
Разрывы на характеристиках КЭ могут существовать и для их исходного состояния. В частности, зазоры между контактирующими поверхностями можно задавать как разрывные свойства контактных КЭ. Стадия определения границ площадки контакта может быть реализована в рамках механики контактной псевдосреды [1-5] как решение задачи со специфическими свойствами контактных КЭ, учитывающих не только нелинейное деформирование слоя, но и наличие исходного зазора между поверхностями. Можно не использовать прием выключения-включения контактных КЭ в зону контактирования, требующий применения внешнего итерационного цикла для поиска площадки контакта. В рамках механики контактной псевдосреды используются контактные КЭ, которые всегда включены в схему расчета. Но они наделены индивидуальными свойствами в зависимости от их геометрического положения: начальный зазор включен в свойства контактных КЭ. Подобные элементы в пределах сближения не оказывают сопротивления внешней нагрузке. Такой подход сводит задачу о контактировании тел с различной геометрией (внешней нелинейностью) к задаче с внутренней нелинейностью, определяемой свойствами контактных КЭ.
В итоге построен итерационный алгоритм [1-5] для решения конечноэлементных задач, в которых разрешающие уравнения могут иметь нелинейности различного типа. Этот алгоритм можно рассматривать как объединение существующих алгоритмов: он реализует возможность применения при расчете того из них, который гарантирует сходимость итерационной процедуры. В одной и той же расчетной схеме, но в разных КЭ с различными нелинейностями одновременно могут использоваться различные варианты итерационных процедур. Кроме того, даже в одном и том же КЭ могут использоваться различные варианты итерационных процедур, если нелинейные характеристики для разных компонент НДС соответствуют нелинейностям различного типа.
Алгоритм позволяет решить многие задачи с нелинейностями различных типов, в том числе и с несколькими типами нелинейностей. Но встроить этот сложный алгоритм в известные промышленные программные пакеты, такие как Nastran и Ansys, очень сложно. Нами предпринята попытка использовать возможности препроцессора FEMAP для подготовки исходных данных для оригинальной программы-решателя [1-5]. В принципе пакет FEMAP является независимым программным продуктом и взаимодействует с пакетами МКЭ Nastran, Ansys и другими, передавая им подготовленные файлы данных и получая файлы результатов для дальнейшего анализа. Нами преобразованы форматы подготовленных данных и оригинальная программа таким образом, чтобы можно было решать нелинейные задачи с различными типами нелинейностей.
Разработан алгоритм обработки информации о созданной в программном комплексе FEMAP with NX Nastran модели. В алгоритме учтены все особенности представления данных, обязательные для корректной работы программы расчета контактных задач (перенумерация узлов и элементов, ориентирование контактных элементов относительно контактирующих тел и др.). Разработанный алгоритм реализован на ЭВМ в виде программы. Изменена часть кода программы решения контактных задач, отвечающая за создание конечноэлементной модели исследуемого объекта. Изменения несут в себе функцию взаимодействия этой программы с обработанными данными о модели, созданной в программном комплексе FEMAP with NX Nastran.
Результаты решения некоторых нелинейных задач
Решено несколько простейших линейных задач, таких как изгиб бруса и другие. Эти задачи по новой версии программы решены без погрешностей и здесь не рассматриваются.
В качестве простого примера нелинейной задачи была решена задача о многократном нагружении шероховатого слоя и поиске площадки контакта с предположением, что центральная часть контактного слоя предварительно была деформирована жестким штампом [7] и в ней возникли остаточные деформации (рис. 4). Зазор, образованный первоначальным воздействием штампа, в модели задается как свойство контактного конечного элемента (по рис. 3), поэтому контактная среда смоделирована сплошным слоем элементов с соответственно различными свойствами. Информация о свойствах материалов и КЭ задается непосредственно в программе-решателе.
Эта же задача была решена по ранее разработанной программе [1-5]. Результаты решения в зависимости от величины нагрузки F приведены на рис. 5. Здесь видно, что при приложении малой нагрузки нормальные напряжения в центре равны нулю (нагрузку эти элементы не воспринимают) (рис. 5а). При увеличении нагрузки в контакт входят не только крайние, не деформированные ранее элементы, но и элементы, находящиеся в центре (рис. 5б). При еще большем увеличении нагрузки элементы на краях и в центре деформируются упругопластически (рис. 5в). При увеличении нагрузки упругопластическая зона растет, а зона упругого деформирования сокращается, затем вся зона в центре начинает деформироваться упругопластически. В итоге распределение напряжений в контакте становится плавным, без скачков и разрывов (рис. 5г).
Опишем подробнее, как организованы взаимосвязи препроцессора и решателя. После того как модель построена, заданы граничные условия и условия нагружения, нужно выгрузить эту информацию в отдельный документ. Во вкладке List нужно выбрать пункт Destination, здесь указать путь к текстовому файлу, в который будет выгружена информация. Далее необходимо задать, какая конкретно нужна информация. В той же вкладке List выбирается пункт Model, поочередно из списка нужно выбрать пункты Node, Element, Constrain – Individual, Load – Individual. Таким образом, выводится информация поочередно об узлах, элементах, условиях закрепления и нагрузках. Соблюдение именно такого порядка строго обязательно, так как именно в такой последовательности программа-решатель обрабатывает входные данные.
|
|
||
|
|
|
Была решена задача о поиске площадки контакта пластины и кольца с расположенным между ними тонким слоем с нелинейными свойствами. Характер свойств контактной среды продемонстрирован на рис. 2. Геометрия модели представлена на рис. 6. По этим данным в среде FEMAP with NX Nastran была построена конечноэлементная модель. Задача решена для различных внешних нагрузок, а именно: 12, 24, 48, 96, 192 Н.
В результате работы программы были получены данные о перемещении всех узлов конечных элементов исследуемой модели на каждой итерации. По этим данным в программном комплексе Mathcad 15.0 были построены диаграммы деформирования контактного слоя для всех выбранных уровней нагрузки.
На рис. 7 приведена диаграмма состояния контактного слоя в направлении нормали (направление Y) без приложения внешней нагрузки. Здесь же приведена линия определения свойств контактной среды в пределах сближения. До этой линии контактный КЭ не оказывает сопротивления внешним воздействиям, так как выбирается зазор между контактирующими поверхностями при нулевых давлениях в контактных КЭ. Если деформация контактного конечного элемента проходит в зоне, отсеченной этой линией и поверхностью пластины, контактная среда приобретает свойства, соответствующие рис. 2а.
|
|
|
I |
Рис. 7. Диаграмма состояния контактного слоя до приложения нагрузки: I - линия,
очерчивающая контактирующую поверхность кольца; II - линия поверхности пластины;
III - линия, соответствующая верхней поверхности слоя
Далее приведена серия диаграмм деформирования контактных конечных элементов в зависимости от уровня приложенных внешних нагрузок: 12 Н (рис. 8а), 24 Н (рис. 8б), 48 Н (рис. 8в), 96 Н (рис. 8г), Здесь показаны только половины симметричных рисунков. Картина деформирования при максимальной нагруззке (192 Н) показана на рис. 9. На рис. 8 и 9 видно, как при увеличении нагрузки ширина площадки контакта увеличивается.
По данным расчета построена диаграмма, отражающая зависимость ширины площадки контакта в проекции на ось Y от уровня внешних нагрузок (рис. 10).
а) б)
в) г)
Рис. 8. Деформация контактных конечных элементов при нагрузках: а - 12 Н; б - 24 Н; в - 48 Н; г - 96 Н
Рис. 9. Деформация контактных конечных элементов при Рис. 10. Зависимость ширины площад-
нагрузке 192 Н ки контакта от величины нагрузки
Заключение
Разработан алгоритм обработки информации о созданной в программном комплексе FEMAP with NX Nastran модели. Алгоритм реализован на ЭВМ в виде программы. Программа была протестирована, результаты работы согласуются с ожидаемыми. Изменена часть кода программы решения нелинейных задач, отвечающая за создание конечноэлементной модели исследуемого объекта. Программа была протестирована. Ход расчета сопоставлялся с ходом расчета программы без внесенных изменений. В итоге была достигнута идентичность результатов, получаемых прежней версией программы и новой.
1. Babin, A.P. Konechnoelementnyy algoritm resheniya kontaktnyh zadach s uchetom nelineynyh effektov / A.P. Babin // Dinamika, prochnost' i nadezhnost' transportnyh mashin. - Bryansk: BGTU, 2002. - S. 138-148.
2. Zernin, M.V. K issledovaniyu kontaktnoy zhestkosti s ispol'zovaniem modeli mehaniki kontaktnoy psevdosredy / M.V. Zernin, A.P. Babin // Zavodskaya laboratoriya. - 2001. - T. 67. - № 6. - S. 51-54.
3. Babin, A.P. Uchet vliyaniya nelineynyh svoystv poverhnostnyh sloev pri konechnoelementnom reshenii zadach o kontaktnom vzaimodeystvii deformiruemyh tel / A.P. Babin, M.V. Zernin // Trenie i smazka v mashinah i mehanizmah. – 2008. - № 3. - S. 3-16.
4. Babin, A.P. Konechnoelementnoe modelirovanie kontaktnogo vzaimodeystviya s ispol'zovaniem polozheniy mehaniki kontaktnoy psevdosredy / A.P. Babin, M.V. Zernin // Izv. RAN. Mehanika tverdogo tela. – 2009. - № 4. - S. 84-107.
5. Zernin, M.V. Algoritm kombinirovannogo metoda resheniya konechnoelementnyh zadach s nelineynostyami razlichnogo tipa / M.V. Zernin, A.P. Babin // Vestnik Bryanskogo gosudarstvennogo tehnicheskogo universiteta. - 2009. - № 4. - S. 57-64.
6. Rudakov, K.N. UGS Femap 9.3. Geometricheskoe i konechnoelementnoe modelirovanie konstrukciy / K.N. Rudakov. - M.: DMK Press, 2009. – 296 s.
7. Ponomarev, V.S. Vnedrenie zhestkogo shtampa v elastichnyy material. Mehanika deformirovannogo tverdogo tela / V.S. Ponomarev. – Tomsk, 1987. – S. 147-152.
