Содержание к диссертации
Введение
ГЛАВА 1. Состояние вопроса, обоснование научной проблемы 17
1.1. Обоснование научной проблемы применения контактной задачи при анализе динамики сборных роторов турбомашин 17
1.2. Выбор методов исследования 22
1.3. Выводы по главе 28
Глава 2. Зависимости метода конечных элементов для расчета динамики роторов турбомашин 29
2.1. Вариационно-энергетический принцип метода конечных элементов 29
2.2. Модель объемного напряженно-деформированного состояния деталей турбомашин 33
2.3. Применение несовместных функций формы при моделировании изгиба толстых пластин объемными конечными элементами 54
2.4. Вспомогательные конечные элементы 56
2.4.1. Балочный конечный элемент 56
2.4.2. Конечный элемент невесомый стержень 63
2.4.3. Конечный элемент сосредоточенная масса 64
2.4.4. Конечный элемент вязкий демпфер 64
2.4.5. Преобразование координат 65
2.5. Решение глобальной системы алгебраических уравнений равновесия 68
2.6. Выводы по главе 75
ГЛАВА 3. Математическая модель решения контактной задачи динамики сборных роторв турбомашин 77
3.1. Математическая модель контактной задачи для расчета статического
напряженно-деформированного состояния сборного ротора 77
3.1.1. Контактный конечный элемент 79
3.1.2. Алгоритм решения статической контактной задачи 88
3.1.3. Пример решения задачи о контакте двух стержней 89
3.1.4. Пример изменения условий сопряжений в сборных узлах ротора при деформировании в условиях рабочего нагружения 92
3.2. Математическая модель динамики сборной конструкции ротора 96
3.3. Алгоритм решения контактной задачи динамики сборного ротора 111
3.4. Подбор величины ускорения раскрутки ротора 113
3.5. Анализ достоверности численного решения динамической задачи 118
3.5.1. Одномассовая и двухмассовая виброударные системы 118
3.5.2. Ротор на анизотропных опорах 121
3.6. Выводы по главе 126
ГЛАВА 4. Численный эксперимент по анализу динамики реального сборного ротора авиационного гтд с учетом контактного взаимодействия деталей.. 128
4.1. Построение конечно-элементной модели сборного ротора для динамического анализа 128
4.2. Решение задачи о вынужденных колебаниях сборного ротора 137
4.3. Выводы по главе 155
Заключение 157
Библиографический список
- Выбор методов исследования
- Применение несовместных функций формы при моделировании изгиба толстых пластин объемными конечными элементами
- Алгоритм решения статической контактной задачи
- Решение задачи о вынужденных колебаниях сборного ротора
Введение к работе
Совершенствование современных турбомашин идет по пути дальнейшего увеличения удельных мощностей, при одновременном ужесточении требований по надежности и ресурсу. Решение этих проблем сопровождается применением принципов многокаскадности, модульности, развитием геометрических форм и других мероприятий, предполагающих использование в турбомашинах конструкций роторов сборного типа. Конструктивно-технологический фактор сборки ротора решает проблему компоновки турбомашины, однако, в этом случае, он имеет в значительной степени более низкий уровень надежности по сравнению со своим гипотетически монолитным аналогом.
Данное обстоятельство объясняется прежде всего тем, что свойство жесткости деформирования сборной конструкции ротора носит сложный концептуальный характер. Оно состоит здесь из двух составляющих: жесткости конструкционной (отдельных деталей) и жесткости контактной. Представленное утверждение сформулировано в ряде расчетно-экспериментальных работ [75, 99, 117 и др.], где под контактной жесткостью понимается свойство односторонней сопротивляемости контактных стыков деформированию сборной конструкции. Величины контактных сил сопротивления зависят как от внешних силовых воздействий, так и от дополнительных нагрузок, определяемых условиями сопряжений деталей: посадками с натягом (зазором), усилиями стягивания и другими факторами.
Проблема заключается в том, что контактная жесткость сборного ротора, в сравнении с конструкционной жесткостью отдельных деталей, имеет более высокую степень влияния на параметры деформирования. По экспериментальным данным [117], контактные деформации сборной конструкции могут составлять 80% и более от величин ее общих деформаций. Кроме того, контактная жесткость подвержена изменению в процессе рабочего нагружения.
Представленные обстоятельства отражаются на динамических характеристиках сборного ротора. Получение достаточно точной оценки этих в действительности, двухстороннюю связь при деформировании, что является некорректным по отношению к стыку; - упругие свойства элементов подбираются на основе отдельного расчета или результатов натурных испытаний и не меняются в зависимости от изменений эксплуатационных нагрузок, тогда как в действительности условия сопряжений деталей в процессе нагружения изменяются. Также, недостатком представленного подхода в проектировании роторов турбомашин, наряду с высоким уровнем материальных и временных затрат, является низкий уровень информативности об объекте, в особенности на начальных стадиях его создания. С чем связан риск принятия ошибочных конструктивно-технологических решений, например, для таких энергоемких и сложных механических систем, как авиационный газотурбинный двигатель (ГТД).
Изложенная проблема может быть преодолена с переходом к современной концепции инженерного анализа, которая характеризуется усилением роли математического и компьютерного моделирования с применением высокоэффективных численных решений, таких, как метод конечных элементов (МКЭ).
В проектировании сборных роторных конструкций турбомашин представленный переход обозначен применением контактной задачи механики твердого деформируемого тела при анализе роторной динамики (рис. 1, в). В этом случае в роторе, наряду со сложностью его конструктивных форм и внешнего воздействия, моделируются условия контактных сопряжений деталей и изменение этих условий: от нерабочих состояний, когда нагрузка определяется только предварительными контактными силами сопряжений, до состояний, определяемых процессами дополнительного рабочего нагружения.
Таким образом, разработка методики математического моделирования динамики сборных конструкций роторов на основе контактной задачи механики твердого деформируемого тела является актуальным научным направлением, имеющим важное народнохозяйственное значение.
Этому направлению в существующей практике отечественных авторов практически не уделено внимание, а доступная информация по тематике роторной динамики, приходящая из-за рубежа, носит рекламный характер.
Цель работы:
Разработка и реализация математической модели и методики анализа роторной динамики с учетом контактного взаимодействия деталей и его изменения в рабочих условиях нагружения сборного ротора.
Идея работы заключается в применении определенной концепции контактной задачи механики твердого деформируемого тела при решении уравнения динамики движения сборного ротора, позволяющей моделировать любые конструкции сопряжений деталей.
Для достижения поставленной цели в работе решаются задачи:
1. Разработать на основе МКЭ в полярной цилиндрической системе координат математическую модель динамического поведения сборного ротора с применением контактной задачи механики деформируемого тела.
2. Провести разработку, программную реализацию и отладку алгоритмов выше представленной математической модели.
З.На основе замкнутых аналитических решений, а также имеющихся результатов натурных испытаний, исследовать достоверность численного решения МКЭ.
4. Разработать конечно-элементную (КЭ) модель реальной конструкции сборного ротора авиационного ГТД.
5. Провести численное исследование динамики реального ротора на трех его КЭ моделях: монолитной; монолитной с подбираемыми упругими элементами на сопрягаемых поверхностях; и сборной, реализующей контактную задачу.
6. Обосновать правомерность применения контактной задачи сопоставлением результатов расчетов перечисленных КЭ моделей.
7. Сформулировать основные конструктивные рекомендации, направленные на улучшение динамических характеристик сборного ротора.
Методы исследования:
Основные результаты работы получены с применением вариационно-энергетического принципа при формировании функционала рассматриваемой физической задачи. И численного решения МКЭ, которое включает полный набор математического аппарата теории матриц, алгебраической сплайн аппроксимации и численного интегрирования. Зависимости МКЭ, используемые для моделирования напряженно-деформированного состояния (НДС) и динамики сборного ротора, построены в полярной цилиндрической системе координат. Контактная задача расчета НДС сборной конструкции решается с использованием модифицированного подхода в методе перемещений теории твердого деформируемого тела. Моделирование контактного взаимодействия деталей осуществляется посредством решения вариационного неравенства методом штрафных функций. Решение уравнения динамики движения сборного ротора построено на основе шагового метода прямого численного интегрирования Ньюмарка. Решение глобальных систем алгебраических уравнений осуществляется прямым методом . Холецкого, адаптированным для работы с разреженными матрицами.
Для разработки программного модуля, реализующего выше представленный анализ и его составляющие, использован алгоритмический язык Fortran и персональный компьютер на базе процессора Pentium IV. Подготовка КЭ моделей, куда входят: геометрическая и дискретная модели объекта, данные по внешним воздействиям и другие параметры; а также интерпретация и обработка результатов анализа, проведены с использованием программного комплекса MSC.Patran. Дополнительное тестирование разработанных математических моделей физических задач и реализованных для них алгоритмов, проведено с использованием программного комплекса MSC.Nastran. Эти комплексы предоставлены учебно-научным центром Иркутского государственного технического университета (ИрГТУ), лицензия ЕС 1916 от 19.08.1998 IrGTU.
Научная новизна работы заключается в следующем:
1. Создана и программно реализована, на основе МКЭ, математическая модель и методика анализа динамики роторов сборных конструкций с учетом контактного взаимодействия деталей конструкционного типа. Позволяющие моделировать сборные конструкции роторов с неограниченным количеством независимых пар контактных поверхностей, имеющих различные условия сопряжений (посадки с натягом или зазором, усилия стягивания), и изменения этих условий в процессе внешнего силового и температурного воздействий.
2. Разработан способ моделирования сборных болтовых соединений, позволяющий значительно уменьшать размерность динамической задачи МКЭ (редуцировать степени свободы КЭ модели) за счет применения различных типов КЭ: объемных - для моделирования деталей ротора и балочных - для стяжных элементов. При этом повышается скорость вычислительного процесса без ущерба в адекватности.
3. Разработана методика подбора ускорения раскрутки ротора при моделировании переходного процесса прохождения им критической частоты. Реальный физический процесс раскрутки ротора является слишком медлительным для моделирования посредством шагового численного интегрирования с применением МКЭ - компьютерное время решения задачи оказывается слишком большим для ее практической реализации (может исчисляться неделями и месяцами). Представленная методика позволяет подобрать ускорение, отвечающее требованиям приемлемого, с практической точки зрения, компьютерного времени расчета и максимальному приближению амплитудно-частотных характеристик (АЧХ) переходного процесса к АЧХ установившихся стационарных режимов частот вращений.
4. Проведен расчет амплитудного отклика сборного ротора в областях резонанса посредством моделирования переходного процесса прохождения критической частоты.
5. Впервые расчетным путем получены АЧХ и картина движения сборного ротора при вынужденных колебаниях с учетом условий сопряжений и их изменений в процессе нагружения. 6. На основе выполненных исследований сформулированы конструктивные рекомендации по улучшению динамических параметров сборных роторов.
Достоверность полученных результатов обеспечена применением расчетно-аналитической базы, отвечающей современному уровню развития и совершенствования турбомашин. Достоверность результатов, полученных численным решением МКЭ, доказывается высокими характеристиками сходимости с аналитически замкнутыми решениями, а так же относительно имеющихся результатов натурных испытаний.
Практическая ценность работы заключается в реализации разработанной математической модели и методов в виде вычислительной программы, позволяющей проводить анализ динамики концептуально сборных роторов турбомашин. Эта программа позволяет на стадии проектирования проводить численные исследования компоновочных решений и их влияние на общий уровень вибраций турбомашины. Что дает возможность повысить надежность работы сборных роторов, увеличить долговечность и удельную мощность, сократить временные и материальные затраты на доводку изделия. Также вычисленные динамические параметры сборного ротора позволяют давать качественное и количественное толкование результатов натурных испытаний.
Результаты, полученные в работе, использованы в процессе реального проектирования роторов ГТД и внедрены на ведущем предприятии авиадвигателестроительной отрасли НТЦ имени А. Люльки НПО "Сатурн", г. Москва.
На защиту выносятся следующие основные положения работы:
1. Разработанные математическая модель и методика • анализа динамики сборных роторов турбомашин.
2. Результаты анализа динамических характеристик реального сборного ротора авиационного ГТД в сравнении с его модельными аналогами.
3. Доказательство невозможности дальнейшего улучшения динамических характеристик жесткого ротора, относительно базовой конструкции, путем подбора упругостей опор.
4. Правомерность использования контактной задачи механики деформируемого тела при теоретическом анализе динамики сборных роторов.
Апробация работы:
Основные результаты выполненных исследований и разработок представлялись и обсуждались на: научной конференции кафедры конструирования и стандартизации в машиностроении ИрГТУ, 2003 г.; IX Всероссийской научно-практической конференции с международным участием "Безопасность-04" кафедры пром. экологии и БЖД ИрГТУ, 2004 г.; региональной научно-практической конференции "Информационные технологии" кибернетического факультета ИрГТУ, 2004 г.; научно-технических конференциях факультета транспортных систем ИрГТУ, 2003-2007 г.г.; университетских смотрах-конкурсах НИР-НИРС-НТТМ ИрГТУ, 2004,2006 г.г.; Всероссийской научно-практической конференции посвященной 75-летию ИрГТУ, 2005 г.; выставках "Инновации для экономики и социальной сферы", проходивших в выставочном центре "Сибэкспоцентр" г. Иркутск, 2005, 2007 г.г.; Всероссийских конкурсах Министерства Образования РФ на лучшую научную студенческую работу года (первые места в 2003, 2004 г.г. по разделам "Авиация, авиастроение, воздушный транспорт" и "Математическое моделирование"); Всероссийском конкурсе выпускных квалификационных работ Совета УМО ВУЗов РФ по образованию в области эксплуатации авиационной и космической техники (второе место в 2006 г.); VII-IX Всероссийских конференциях "Комплексные технологии виртуального моделирования и инженерного анализа", проводимых международной корпорацией MSC.Software, г. Москва, 2004-2006 г.г. (доклады заняли первое и вторые места в конкурсном отборе докладов ведущих промышленных предприятий и научных организаций России и стран СНГ в номинации "Аэрокосмическая промышленность"); научно-практическом семинаре отдела прочности НТЦ им. А. Люльки НПО "Сатурн", г. Москва, 2004 г.; научно техническом семинаре НПО "Энергомаш" им. В. П. Глушко, г. Химки, 2006 г.; научно-методическом семинаре кафедры строительной механики МАДИ (ГТУ), г. Москва, 2006 г.; научно-методическом семинаре кафедры конструкции и проектирования двигателей МАИ (ГТУ), г. Москва, 2006 г.; научном семинаре НТЦ МКБ "Гранит", г. Москва, 2006 г.. Личный вклад соискателя:
1. Обзор и анализ информации об исследованиях, проведенных ранее.
2. Разработка и программная реализация алгоритма нелинейного анализа динамики сборных роторных конструкций с учетом контактного взаимодействия деталей, включающие следующие блоки:
- реализация шести типов конечных элементов (объемного НДС, контактного, балочного, невесомого стержня, сосредоточенной массы, вязкого демпфера);
-реализация алгоритма решения глобальной системы алгебраических уравнений равновесия, адаптированного для работы с разреженными матрицами;
-учет граничных условий, моделирования стационарных (инерционных, температурных) и динамических (дисбалансных) нагрузок, внутренних нагрузок контактных взаимодействий деталей и условий их сопряжений (отслеживания изменений этих условий в ходе рабочего нагружения);
- моделирование эффектов демпфирования; моделирование дополнительной жесткости ротора, прецессирующего в поле центробежных сил;
- реализация алгоритма прямого численного интегрирования динамического уравнения сборного ротора на основе метода Ньюмарка.
3. Разработка методики выбора исходных численных параметров для успешного решения поставленной задачи.
4. Анализ достоверности численного решения на основе сравнения с точными решениями задач: об изгибе толстой пластины, динамики простейших виброударных систем и критических частотах вращения монолитного однодискового ротора на анизотропных опорах и данными натурных испытаний реального сборного ротора.
5. Создание КЭ модели реальной конструкции сборного ротора ГТД.
6. Численное исследование динамики изучаемого объекта.
7. Обработка и анализ полученных результатов; формулировка положений диссертации и выводов по результатам исследований.
Все приводимые в работе результаты исследований получены автором лично. Отмечающиеся в тексте результаты других исследователей, а также результаты совместных исследований с соавторами, снабжены ссылками на соответствующие источники.
Публикации:
Основные результаты работы освещены в 7 публикациях, в том числе: одна монография; одни тезисы доклада; 5 статей, из которых 2 опубликованы в центральных периодических изданиях. Личный вклад автора в опубликованные в соавторстве работы составляет не менее 50%.
Структура и объем работы:
Диссертационная работа состоит из оглавления, списка принятых сокращений, введения, четырех глав, заключения, библиографического списка и 3 приложений. Диссертация изложена на 174 страницах, включает 91 рисунок и 15 таблиц. Библиографический список охватывает 156 источников.
Работа выполнялась в 2003-2007 годах в Иркутском государственном техническом университете на кафедре самолетостроения и эксплуатации авиационной техники.
Содержание работы:
Во введении определяются объект и предмет исследования, формулируется цель работы, задачи и методы их решения; приводятся выносимые на защиту положения и краткое содержание работы по главам.
В первой главе обосновывается актуальность применения контактной задачи механики деформируемого тела при теоретическом анализе роторной динамики, позволяющей учесть различные конструктивно-силовые факторы, влияющие на динамические свойства сборной конструкции в целом.
Проводится обзор подходов к решению контактных задач. Указаны основные принципы и методы выполнения исследований.
Вторая глава посвящена математическому аппарату МКЭ, применяемому для моделирования сборного ротора и построенному на основе алгебраической сплайн аппроксимации и вариационно-энергетического принципа метода перемещений теории упругости. В главе представлены зависимости конечных элементов (КЭ), используемых в работе, основного - КЭ объемного НДС и вспомогательных КЭ: балочного, невесомого стержня, сосредоточенной массы, вязкого демпфера. Адаптация перечисленных конечных элементов для задачи моделирования ротора, заключается в применении полярной цилиндрической системы координат. Анализ точности и сходимости численного решения проведен для основного КЭ объемного НДС на упругой задаче об изгибе толстых пластин, имеющей аналитическое решение. Также в главе рассмотрен алгоритм решения глобальной разреженной системы алгебраических уравнений равновесия.
Третья глава содержит описание математического аппарата контактной задачи, используемого для анализа динамики сборной конструкции ротора. Алгоритм решения задачи построен на основе модифицированного вариационно-энергетического подхода, реализуемого относительно невязки поля перемещений сопрягаемых поверхностей. Контактное взаимодействие деталей моделируется посредством специального двух узлового контактного КЭ. Алгоритм работы контактного КЭ базируется на принципе применения штрафной жесткости для его узлов, принадлежащих одновременно и сеткам контактирующих поверхностей деталей. Величина этой жесткости адаптирована соответственно фактора непроникновения деформируемых тел друг в друга, определяемого условиями их сопряжений. В главе приведена постановка задачи роторной динамики. Представлены алгоритмы и методика выбора исходных численных параметров для успешного решения поставленной задачи (методика подбора ускорения раскрутки ротора). Для решения динамического уравнения сборного ротора используется шаговый метод прямого численного интегрирования Ньюмарка. Достоверность численного решения доказывается сравнением с аналитическими решениями для простейших виброударных систем и монолитного однодискового ротора на анизотропных опорах.
В четвертой главе последовательно описана методика численного исследования роторной динамики на примере анализа реального ротора авиационного ГТД. Подробно изложены особенности и этапы создания КЭ модели изучаемого объекта. Предложен способ моделирования сборных болтовых соединений, значительно редуцирующий общее число степеней свободы КЭ модели без существенной потери ее адекватности. На основе задачи о собственных частотах и формах колебаний доказывается невозможность дальнейшего улучшения динамических характеристик жесткого ротора путем подбора упругостей опор. Сопоставлением результатов расчетов трех КЭ моделей (монолитной, сборной с контактными элементами и "сборной" с подбираемыми упругими элементами) обосновывается правомерность использования контактной задачи при теоретической оценке динамических параметров сборного ротора. Достоверность получаемых теоретических результатов для сборного ротора подтверждается сравнением с данными натурных испытаний. В конце главы сформулированы конструктивные рекомендации, направленные на улучшение динамических характеристик концептуально сборного ротора.
Заключение содержит общую характеристику диссертационной работы и основные выводы по результатам.
В первом приложении приведена система единиц, принятая в выполненных в работе расчетах. Второе приложение содержит акт внедрения результатов работы на предприятии НТЦ имени А. Люльки НПО "Сатурн", г. Москва. В третьем приложении на компакт-диске представлена анимация результатов динамического анализа сборного ротора.
Автор выражает благодарность, глубокое почтение и признательность научному руководителю д.т.н., профессору кафедры самолетостроения и эксплуатации авиационной техники ИрГТУ А. А. Пыхалову.
Выбор методов исследования
Заметный подъем в развитии теоретических методов решения контактных задач механики деформированного тела наблюдался в конце 60-х начале 70-х годов прошлого столетия. Это связано с бурным внедрением в расчетно-инженерную практику ЭВМ повышенной производительности и, соответственно, более широким использованием численных методов решения. Несмотря на то, что численные методы, по сути, являются приближенными, достоверность решения сложных задач (при успешной реализации метода) значительно выше, нежели при использовании классических аналитических методов. Достижение высокой достоверности результатов обеспечивается отсутствием упрощений по геометрии объекта исследования, точным соблюдением граничных условий и закономерностей приложения нагрузок, учетом изменения физико-механических свойств исследуемой области.
В особенности это относится к одному из наиболее эффективных численных методов - методу конечных элементов (МКЭ). Этот метод оказал существенное влияние на развитие всей расчетно-аналитической базы механики твердого деформированного тела. Основополагающие работы в этом научном направлении принадлежат как зарубежным, так и отечественным авторам: Зенкевичу О., Бате К., Вильсону Е., Галлагеру Р., Гордону М. А., Деклу Жану, Ершову Н. Ф., Квитка А. Л., Дж. Мак-Нилу, Морозову Е. М., Мяченкову В. И., Никишкову Г. П., Норри Д., Ж. де Фризу, Розину Л. А., Сегерлинду Л. и другим.
В контактной задаче деформированного тела, с использованием МКЭ и других численных методов, известны работы исследователей: Адлуцкого В. Я., Александрова В. М., Алексидзе М. А., Барлама Д. М., Блоха М. В., Божковой Л. В., Власенко Ю. Е., Галанова Б. А., Галкиной Н. С, Гнучий Ю. Б., Можаровского Н. С, Нагиной Е. Л., Потапова С. Д., Рудакова К. Н., Цвика Л. Б., Щеглова Б. А. и других.
При разработке математической модели контактной задачи на основе МКЭ, наиболее успешным следует считать использование вариационно-энергетического подхода, построенного на основе метода перемещений теории упругости и пластичности твердых деформируемых тел. В этом случае, построение функционала потенциальной энергии для сборной конструкции, как деформируемой системы, осуществляется с дополнительным учетом энергии контактного взаимодействия деталей. После минимизации представленного функционала, относительно вектора перемещений, получаем модифицированную глобальную систему алгебраических уравнений МКЭ в ее классическом виде: [K]{S}={F], (1.1) где [к]- глобальная матрица жесткости, которая остается симметричной и положительно определенной для рассматриваемой деформируемой системы; {F}, {S} - глобальные вектор-столбцы сил и неизвестных перемещений.
Первоосновой модификации глобальной системы уравнений (1.1), связанной с реализацией контактной задачи, является использование, в том или ином виде, классического метода Шварца о вариационном неравенстве. Решение этого неравенства осуществляется на различной математической основе: обычными итерационными методами; с использованием оптимизационных методов (множителей Лагранжа, штрафа, релаксационных схем и др.); а также с построением матрицы граничной податливости.
Важным отличительным аспектом выбора методики реализации контактной задачи, является вопрос о способе модификации выражения (1.1). В конечном счете, он сводится к вопросу о принадлежности производной величины потенциальной энергии по перемещению, связанному с контактным взаимодействием между деталями, к правой или левой части системы (1.1).
В работах [32, 85, 138 и др.], величина этой производной выделяется в виде контактных сил, величина которых используется в модификации глобального вектора сил {F} выражения (1.1). Вычисление этих сил производится через итерации относительно характеристик жесткости контактирующих деформируемых тел, а также величин невязки (неравенства) поля перемещений между их контактными поверхностями. Такой подход, с физической точки зрения, вполне оправдан, так как позволяет определить условия контактного взаимодействия деформируемых тел не только при наличии внешней нагрузки, но и при воздействии в них внутренних сил предварительного сопряжения между телами, например, при их взаимных посадках с натягом и других условиях. Внешняя нагрузка, при этом, может полностью отсутствовать. Однако, при использовании представленного подхода в расчете сборных роторных конструкций, где контактное сопряжение происходит по нескольким поверхностям, не зависящим друг от друга и находящимся в различных координатных плоскостях, начинает доминировать фактор неустойчивости решения глобальной системы алгебраических уравнений (1.1).
Применение несовместных функций формы при моделировании изгиба толстых пластин объемными конечными элементами
В любой программе, реализующей МКЭ, ключевой является процедура решения глобальной системы алгебраических уравнений равновесия объекта. Процесс решения занимает большую часть времени счета (до 99% от общего времени), и поэтому, требуется тщательный подбор метода решения и составления алгоритма. В настоящее время в МКЭ наибольшее распространение для решения глобальной системы уравнений получил прямой метод Холецкого. Представленный метод является самым экономичным по числу арифметических операций, позволяет эффективно использовать оперативную память ЭВМ и является абсолютно устойчивым [50].
Метод Холецкого можно рассматривать как симметричный вариант Гауссова исключения для симметричной положительно определенной матрицы.
Рассмотрим глобальную систему алгебраических уравнений равновесия ансамбля конечных элементов: [K]{S} = {F}, (2.139) где [К ] - [п п] симметричная положительно определенная матрица жесткости; {F} - вектор внешней нагрузки размерности п (правых частей); \5} - неизвестный вектор - решение размерности п (узловые перемещения); п - число степеней свободы конечно-элементной модели. Применение к [к] метода Как видно из рисунка 2.17 вычисляются один за другим столбцы [L\, но к той части матрицы, которая на данном шаге еще не разложена, обращения не происходит. После того, как вычислено разложение, остается решить треугольные системы (2.143) и (2.144).
Решение системы (2.143) называется прямой подстановкой. Так как матрица [і]- нижняя треугольная, то первое уравнение имеет следующий вид: ll l=/l- (2.146) Из (2.246) получаем: l= l/ ll- (2.147)
Теперь можно вычесть первый столбец матрицы [/,], умноженный на х\, из вектора {/ }. Система (2.143) преобразуется в треугольную систему порядка (л -1) с треугольной подматрицей матрицы [L ] в качестве матрицы коэффициентов. Теперь можно вычислить Х2, и приведенные соображения могут быть повторены далее, пока не будут вычислены все элементы { х}.
Если известен вектор {х}, то система (2.144) решается обратной подстановкой (по аналогии с описанной выше прямой подстановкой), начиная с последнего уравнения. Рассмотрим процедуру решения системы уравнений, структура ненулевых элементов матрицы коэффициентов которой показана на рис. 2.18, а. Множитель [L] рассматриваемой матрицы представлен нарис. 2.18, б.
Перенумеруем переменные и переупорядочим уравнения так, чтобы последнее стало первым, второе снизу - вторым сверху и так далее, пока бывшее первое уравнение не станет последним. В результате получится эквивалентная система уравнений, для которой матрица коэффициентов [К \ представлена на рис. 2.18, в. Множитель Холецкого переупорядоченной матрицы [I \ приведен на рисунке 2.18, г.
Решая (2.143) и (2.144) с использованием рассмотренных выше матриц, получим решения {$} и \Sj, где решение \S \ является переупорядоченной формой {д}.
Важнейший момент состоит в том, что переупорядочение уравнений и переменных привело к треугольному множителю [Ь\, который разрежен в точности в той мере, что и нижний треугольник [К \.
Приведенный выше пример демонстрирует, что разумное упорядочивание строк и столбцов матрицы коэффициентов может дать огромное сокращение заполнения и, следовательно, экономию времени вычислений и памяти ЭВМ.
Для поиска переупорядочивания используются специальные алгоритмы, основанные на теории графов. В настоящей работе для переупорядочивания системы уравнений использован метод вложенных сечений [50].
Используемая в настоящей работе схема хранения глобальной матрицы предусматривает хранение только логически ненулевых элементов факторизуемой матрицы жесткости [К] и множителя [/_,], рис. 2.20 (см. далее).
Имеется основной вектор KN (LN), содержащий все ненулевые элементы нижнего треугольника [к] (множителя [ь]). Ненулевые элементы, исключая диагональные, хранятся в KN (LN) столбец за столбцом. Диагональные элементы хранятся отдельно в векторе KD (LD). Кроме того, предусмотрены два сопровождающих адресных вектора.
Адресный вектор iK (iL) - указывает начало ненулевых элементов каждого столбца в векторе KN (LN) (последний элемент вектора iK (iL) несет информацию о размере векторов jK (jL)). Адресный вектор jK (jL) - дает строковые индексы ненулевых элементов нижнего треугольника глобальной матрицы (множителя [L ]).
Информация о структуре нулевых и ненулевых элементов множителя глобальной матрицы получается в результате символического разложения. Символическое разложение моделирует численное разложение матрицы. Поскольку в вопросе символического разложения числовые значения элементов матрицы не играют роли, эта задача решается методами теории графов [50].
На рис. 2.21 (см. далее) приведена блок-схема общего алгоритма постановки задачи и решения глобальной системы алгебраических уравнений равновесия.
Алгоритм решения статической контактной задачи
Метод конечных элементов является одним из наиболее эффективных и универсальных численных методов решения инженерных задач. Подмеченное в процессе его применения сходство в общей структуре алгоритмов различных физических приложений, а также дальнейшее развитие вычислительной техники, привело к созданию на его основе мощных универсальных программных комплексов. Фактически, они представляют собой накопленный к настоящему времени объем математического аппарата необходимого для проведения инженерных вычислений. А некоторые из них, в частности такой программный комплекс как MSC.Nastran, определяют собой стандарты конечно-элементного решения физических задач.
В то же время готовые программы не всегда удовлетворяют требованиям конкретной инженерной задачи, например, когда исследователю необходимо рассмотреть некоторую качественно новую сторону изучаемого объекта.
В настоящей работе, одной из таких сторон является применение контактной задачи механики деформируемого тела для динамического анализа сборных конструкций роторов. В этом случае возникает необходимость в создании своих собственных программ, используя накопленный опыт алгоритмизации МКЭ. Вероятно, необходимость в подобной практике сохранится и в дальнейшем, так как в МКЭ и в тех областях, для которых ведется расчет, все время возникают новые подходы к моделированию физических задач.
Осевые жесткости стержней - К\ и К2 соответственно. В узловых точках деформируемой системы реализована одна степень свободы - перемещение вдоль оси х. Граничные условия - запрет перемещений четвертого узла, принадлежащего правому стержню. Между стержнями имеется зазор UQ, моделируемый КЭСКом. В первом узле, принадлежащем левому стержню, приложена внешняя нагрузка - осевая сжимающая сила F. На рис. 3.9 (см. далее) показана КЭ модель сборного узла реального ротора авиационного ГТД (см. главу 4). Модель представляет циклосимметричный сектор полной окружности и предназначена для расчета статического НДС. Моделирование отдельных сборных узлов позволяет анализировать условия сопряжений деталей на различных режимах работы ротора (при различных случаях нагружения). Изменение условий сопряжений конструкции в процессе работы доказывается результатами расчета ротора в состоянии покоя и при максимальной частоте вращения, рис. 3.10,3.11 (см. далее).
Применительно к задаче динамики движения сборной конструкции ротора, уравнение (3.21) приобретает вид: [Mp}+[C]{S}+{FH(S)}+ + {[K]+[KQ(nMKK(S)])lS}={FQ(n t)}+{FK(s)}, где [Kg(S)], {Ffc(S)} - дополнительные слагаемые, описывающие контактное взаимодействие деталей ротора; {- ( )}" вектор сил внутреннего неупругого сопротивления движению; [A"Q(W)] - матрица дополнительной жесткости, обусловленной прецессионным движением ротора в поле центробежных сил; {FQ(«, /)} - вектор внешних сил, зависящих от частоты вращения ротора п и времени t. Вектор {FQ H, t)} может включать: - центробежные силы инерции, обуславливаемые вращением ротора; - динамические силы неуравновешенных вращающихся масс; - эквивалентные силы тепловых расширений от температурной нагрузки; - эквивалентные силы давлений в полостях и на поверхностях ротора; - динамические лопаточные нагрузки от неравномерностей газового потока; - гидро- газодинамические нагрузки уплотнений и опор и т. п..
Матрицы {Fff(c))}, [Кк(б)], {FK(S)} - являются неявно заданными функциями вектора перемещений {S}. В такой постановке уравнение (3.23) не имеет аналитического решения и может быть решено только с использованием численных методов.
Для решения динамического уравнения (3.23) выбран метод прямого численного интегрирования Ньюмарка. Термин "прямое" означает, что перед интегрированием не производится никаких предварительных преобразований.
Численное интегрирование уравнения (3.23) организуется пошагово и осуществляется на основе двух принципов. Первый из них требует выполнение равновесия (3.23) не в любой момент времени t, а только на отдельных коротких его промежутках - At. Решение строится из серии статических решений (шагов интегрирования) на временном интервале динамического воздействия {F( t)}.
Второй применяемый принцип заключается в том, что изменения величин перемещений {д}, скоростей { } и ускорений \S\ учитываются внутри каждого временного интервала At. Предполагается, также, что в начальный момент времени IQ =0 их значения известны. Решение в последующий момент времени t + At находится на основе предыдущих шагов интегрирования.
Выбор метода Ньюмарка продиктован двумя основными его свойствами: - метод является одношаговым; то есть, для выполнения вычислений на следующем временном шаге необходима информация об одном предыдущем шаге, что позволяет легко менять временной шаг в процессе интегрирования; - метод является абсолютно устойчивым, то есть, решение не растет неограниченно при любых начальных условиях и некорректно заданном временном шаге Д/.
Решение задачи о вынужденных колебаниях сборного ротора
С увеличением ускорения , резонансный пик сдвигается вправо, при этом в закритической области возникают биения; - биения объясняются тем, что система не успевает за быстроизменяющейся нагрузкой, и движение здесь складывается из гармоник собственной частоты и нагрузки; - вследствие демпфирования собственная гармоника затухает, и биения прекращаются; - уменьшение ускорения приближает характеристики к установившимся режимам; - величина ускорения , приводящая к существенному сдвигу резонансного пика и вызывающая биения, пропорциональна величине критической частоты.
Результаты исследования служат основой для подбора величины ускорения раскрутки при заданных границах исследуемого диапазона частот вращения и ориентировочных положениях интересующих резонансных областей. Благодаря использованию простейшей схемы, избегаются огромные вычислительные затраты, свойственные сложной модели ротора.
Целью анализа является формирование представлений о точности и особенностях применения выбранного численного метода, что позволит подойти к решению интересующей задачи с известным уровнем достоверности.
Проверка достоверности численного решения разработанной модели и алгоритма динамической контактной задачи осуществлена на основе рассмотрения свободных колебаний виброударной системы с одной степенью свободы (рис. 3.31, справа). Ненапряженному состоянию пружины соответствует нулевое перемещение массы = 0 (нейтральное положение). Преграда, о которую происходит удар точечной массы М, отстоит от нейтрального
Временной шаг численного интегрирования при 20 шагах полного периода: Дг = 0.1047. Начальные условия задачи: »в-в; ,-і:і,-о; в-о. Результаты численного решения задачи представлены на рис. 3.31 и соответствуют классическим представлениям: диаграмма перемещения составлена из верхних фрагментов синусоиды sinyK/M lj, усеченной на величину натяга и .
Особенностью виброударных систем является мгновенная потеря энергии при соударении [8, 42, 43, 94, 96], характеризуемая коэффициентом восстановления относительной скорости: (3.73) о Ы?- ?- 1, ?2+- + гДе д\+, %+; 5j_, &i- - до- и послеударные скорости тел 1 и 2.
На рис. 3.32 показана фазовая траектория системы с ударным демпфированием. интегрирования неточно соответствуют моментам ударов. При увеличении исследуемого временного интервала описанное свойство неточности алгоритма интегрирования проявляется еще более четко, рис. 3.33. Таким образом, если ударное демпфирование исследуемой системы невелико, отпадает необходимость его специального описания при моделировании. Степень численной диссипации можно варьировать изменением штрафной жесткости КЭСК и величиной временного шага интегрирования. представлено исследование достоверности численного решения контактной динамической задачи с использованием двухмассовой виброударной системы, где массы М\ = Mj и жесткости К\ = К.2
Временной интервал между соударениями масс (без учета демпфирования): Т = к М/К =0.993. (3-74 Перемещения каждой массы представляются фрагментами синусоиды: Ssm{jK/M + ip), $.75) амплитуда S и фаза (р которой будут периодически изменяться после каждого соударения. При аналитическом подходе рассматриваемая задача существенно усложнится, если задать К\ФК.2, Ы\ФМ2 и ввести в систему внешнее воздействие. Колебания в таком случае могут принять хаотический характер, и предсказание времени соударений окажется невозможным. При использовании же численного интегрирования, несложно организовать алгоритм адаптации временного шага в окрестностях моментов соударений.
