Практически все современные технические сооружения и аппараты — ракеты и космические станции, самолеты, корабли, автомобили, строительные и гидротехнические сооружения — представляют собой сложные системы, состоящие из совместно функционирующих подсистем. Условия взаимодействия этих подсистем, выделяемых либо пространственно, как часть конструкции, либо в плане выполняемой функции, определяют успешность выполнения главной задачи разрабатываемой системы. Как правило, понятие «сложность» связывается именно с наличием в системе многих компонент, взаимное влияние которых создает проблемы при проведении теоретических исследований, необходимых для ее проектирования.
Физическую основу рассматриваемых систем, несущую все прочие подсистемы, представляет конструкция, скомпонованная из стержневых, тонкостенных или иных элементов, изготовленных из материалов, которые в пределах достаточно малых деформаций могут рассматриваться как упругие. Результатом взаимодействия упругой конструкции с прочими подсистемами и с внешней средой являются ее колебания — периодические или же переходный процесс. Параметры этих колебаний определяют пригодность конструкции к эксплуатации по критериям прочности, амплитудным значениям перемещений, уровням перегрузок или иным конкретным для каждой системы показателям.
Важным этапом исследования динамического поведения разрабатываемой системы является определение динамических характеристик входящей в ее состав упругой конструкции, к числу которых относятся собственные частоты и формы колебаний, амплитудно-фазовые частотные характеристики, динамические коэффициенты влияния (динамические жесткости и динамические податливости) и т. д. Эта информация является исходной для последующего анализа вибраций конструкции.
Обычно упругая конструкция сама представляет собой сложную систему, составленную из относительно более простых подконструкций, механически соединенных между собой и взаимодействующих в процессе совместных колебаний. Это существенно осложняет задачу исследования ее динамических характеристик как экспериментальными, так и расчетными методами. При этом возникающие трудности могут иметь как технический, так и организационный характер:
— размерность математической модели всей конструкции в целом может превышать возможности используемой для расчета вычислительной системы (либо ограничен объем памяти, либо потребное время счета делает задачу невыполнимой);
— конструкция может оказаться слишком велика для проведения вибрационных испытаний (в особенности это относится к летательным и космическим аппаратам, динамические характеристики которых должны определяться при отсутствии какого-либо закрепления);
— многие крупные системы (например, космические станции) обычно формируются из фрагментов, разрабатываемых разными фирмами, находящимися в разных странах на значительном удалении друг от друга, когда сборка всех компонент для проведения испытаний оказывается весьма дорогостоящим и трудновыполнимым мероприятием.
Естественным направлением мысли на пути преодоления указанных проблем является анализ расчлененной на подсистемы конструкции по частям и последующий синтез результатов, полученных для каждой части в отдельности теоретически или экспериментально. Развитие электронной вычислительной техники с середины 1960;х годов придало актуальность разработке универсальных алгоритмов, позволяющих автоматизировать процедуру синтеза при исследовании динамических характеристик сложных механических систем.
Считается, что впервые четко оформленный тензорно-матричный подход к этой проблеме изложен в работах Г. Крона [65]. Предложенная им методология преимущественно ориентирована на анализ электрических сетей и оказалась мало приспособленной к специфике механических задач. Тем не менее, имеются немногочисленные последователи, развивающие это направление [183, 185, 86].
Основные же пути развития теории синтеза динамических характеристик подконструкций определялись с учетом особенностей задач динамики упругих систем. Значительный вклад в этот раздел науки внесли отечественные исследователи, и здесь следует отметить работы Постнова В. А. [70, 89, 90, 85], Вольмира А. С. [28, 29, 30, 79], Шклярчука Ф. Н. [101], Шмакова В. П. [103, 104, 105], Лиходеда А. И. [68, 4, 5], Бурмана З. И. [24]. Среди зарубежных исследователей наиболее заметны работы таких авторов, как Craig R.R. [126, 127, 128, 129, 130], MacNeal R.H. [168].
Среди многообразия подходов к синтезу динамических характеристик выделим, как наиболее физически обоснованный, метод модального синтеза, когда в качестве исходной информации о свойствах подконструкций используются данные об их собственных частотах и формах колебаний. Основой для построения математической модели всей системы в этом случае служит представление колебаний каждой подконструкции в виде ряда, содержащего ее собственные формы (в дальнейшем — модального разложения колебаний). Т. е. собственные формы играют роль координатных функций в описании движения подконструкции.
Заметим, что существуют подходы к решению данной задачи, не основанные на предварительном вычислении динамических характеристик подкон-струкций. Это, например, работа [113], в которой разбиение дискретной модели упругой системы на подконструкции используется, фактически, лишь для более эффективной реализации метода итерирования подпространства при вычислении собственных частот и форм системы. Здесь же упомянем работы [140, 169], в которых предлагается для аппроксимации колебаний подконст-рукций использовать произвольные полные системы базисных функций, а также работы [117, 141, 142], где с помощью специальных итерационных алгоритмов эти базисные функции улучшаются (итерация подпространств на уровне подконструкций).
Тем не менее, наибольшее количество работ посвящено модальному синтезу, поскольку в этом случае удается эффективно сокращать объем исходной информации о подконструкциях и, что самое важное, уменьшать размерность решаемой в процессе синтеза динамических характеристик задачи, основываясь на ограничении исследуемого частотного диапазона.
Важное значение при исследовании сложной системы имеет способ соединения ее компонент (интерфейс системы). Как правило, соединение под-конструкций осуществляется посредством специальных пространственно локализованных узлов, работающих таким образом, что в рамках принимаемой математической модели подконструкции воздействие со стороны этого узла представляет совокупность сосредоточенных обобщенных сил, связанных с соответствующими обобщенными перемещениями в точке. Обычно входящие в узел степени свободы связаны линейными соотношениями, входящими в получаемую при синтезе математическую модель системы непосредственно либо с помощью множителей Лагранжа, как в работе [132].
Сопряжение подконструкций по одномерным и двумерным многообразиям обычно имеет место при искусственном рассечении крупногабаритной конструкции. При использовании в расчетах дискретных моделей подконструкций (как правило, построенных на основе метода конечных элементов) здесь не возникает принципиальных затруднений, поскольку соединение осуществляется посредством коллокации в узлах модели. В работах [140, 169] предлагается метод, основанный на введении специальных весовых функций, связанных с континуальным интерфейсом, что практически означает его дискретизацию (хотя и не пространственную). В работе [155] с этой целью введены граничные обобщенные координаты. Такой подход может быть полезен при использовании аналитических моделей подконструкций. Отметим также работу [174], где в вариационной постановке задачи синтеза используются определенные на границе сопряжения множители Лагранжа, а решение дискре-тизированных по методу Ритца уравнений осуществляется с использованием сингулярного разложения подматриц, соответствующих интерфейсу системы.
Ключевым вопросом при реализации модального синтеза подконструкций является выбор граничных условий, при которых определяются собственные частоты и формы компонент системы (парциальные динамические характеристики). При этом, естественно, не могут варьироваться наложенные на подконструкцию кинематические ограничения, не относящиеся к интерфейсу системы. Свобода выбора существует только для обобщенных перемещений, связанных с соединительными узлами, которые в дальнейшем будем называть внешними степенями свободы подконструкции. Существующие методы модального синтеза можно классифицировать по этому признаку следующим образом:
— методы жестких границ, когда парциальные характеристики подконст-рукций определяются при условии закрепления внешних степеней свободы;
— методы свободных границ, когда парциальные характеристики подкон-струкций определяются при не закрепленных внешних степенях свободы;
— гибридные методы, если возможно частичное закрепление внешних степеней свободы подконструкции при определении ее парциальных характеристик.
В дополнение к перечисленным, существует еще метод, названный в работе [131] методом «нагруженных границ», когда расчет форм колебаний под-конструкции осуществляется не изолированно от остальных компонент системы, а при дополнительных жесткостных и инерционных нагрузках, добавляемых к внешним степеням свободы с целью приблизить эти формы к виду собственных колебаний системы в целом на данной подконструкции. При правильном выборе этих нагрузок решение задачи о собственных значениях при синтезе подконструкций может дать более точный результат. Варианты такого подхода под названием «метода ветвей» описаны в работах [135, 116]. Нагруженные собственные формы используются также в работе [156].
Неудобство этого метода по сравнению с тремя перечисленными выше очевидно ввиду того, что модификация одной из составляющих систему подконструкций при таком подходе вызывает необходимость пересчета парциальных характеристик остальных подконструкций и внесения изменений в соответствующие им базы данных. При этом интуитивность подбора дополнительных нагрузок не гарантирует существенного повышения точности результата.
Варианты метода жестких границ представлены в работах [149, 150, 151, 128, 126, 130, 61]. Отметим, что вариант в форме Крейга-Бэмптона [128] послужил основой распространенного в настоящее время формата обмена данными по динамике подконструкций между кооперированными разработчиками сложных конструкций. В работе [144] описан метод жестких границ применительно к системам с демпфированием, позволяющий учитывать несимметричность, связанную с кориолисовыми силами и взаимодействием системы с внешней средой. В работе [5] описан метод выделения квазистатических составляющих кинематического и силового типов как дополнительных членов модального разложения колебаний подконструции, обеспечивающий существенное повышение точности решения. При этом вычисляемая на первом шаге квазистатическая составляющая кинематического типа представляет собой аналог введенных в методе Крейга-Бэмптона «граничных форм» («constraint modes»), дополняющих модальное разложение.
Методы свободных границ представлены в работах [136, 146, 129, 126, 181]. Отметим, что в работе [129] введено понятие остаточной податливости для приближенного учета влияния в низкочастотном диапазоне не включенных в модальное разложение высших тонов подконструкции, с которым связано понятие «соединительных форм» («attachment modes»), как дополнительных членов этого разложения [126]. Метод учета остаточных эффектов второго порядка предложен в работе [181]. В работе [154] для построения матриц остаточной податливости балочных подсистем используются аналитические выражения для собственных форм и частот.
Гибридный метод описан в работе [168], где также предложены методы учета остаточных эффектов для повышения точности решения.
В отечественной практике получили распространение многоуровневые методы синтеза, в которых допускается поэтапное укрупнение фрагментов сложной системы. При этом синтез группы подконструкций на более низком уровне дает информацию о подконструкции следующего уровня, получаемой посредством их соединения. Метод суперэлементов, представленный в работах [28, 29], ориентирован на использование собственных форм, определяемых с учетом влияния соседних суперэлементов вдоль общих границ, чем сходен с методом работы [116]. Предложенный в работе [57] метод многоуровневой динамической конденсации может быть отнесен к методам жестких границ и на низшем уровне соответствует идеологии работы [128]. При этом, как и в работе [115], подконструкция рассматривается как суперэлемент с внутренними обобщенными неизвестными, соответствующими ее собственным формам.
Отметим как предельный («вырожденный») случай модального синтеза метод, основанный на статической конденсации подконструкции [153, 139, 70], когда вводится жесткая связь внешних степеней свободы с внутренними и последние исключаются из уравнений колебаний. Такой подход весьма прост в реализации, т.к. не требует предварительного определения собственных частот и форм, но имеет весьма ограниченную сферу применения ввиду невысокой точности, поскольку в уравнениях полностью исключается внутренняя динамика подконструкции. Как промежуточный можно рассматривать предложенный в работе [167] вариант упрощенной динамической конденсации, когда вычисленные при фиксированных границах собственные формы используются для приближенного учета внутренней динамики подконструкции посредством линеаризации в окрестности заданного значения частоты.
Заметим, что все упомянутые выше методы направлены на обеспечение удовлетворительной точности результатов синтеза в диапазоне низких частот колебаний, для чего в модальных разложениях удерживаются формы, соответствующие низшим собственным частотам подконструкций. В работах [152, 162] предложены методы учета остаточных эффектов не только высших, но и низших отсекаемых в модальных разложениях тонов, когда интервал исследуемых частот для системы ограничен снизу ненулевым значением.
Подавляющая часть методик синтеза динамических характеристик разрабатывалась таким образом, чтобы в результате синтеза формировались линейные алгебраические системы или системы дифференциальных уравнений. В этом случае на завершающем этапе расчета динамических характеристик системы формулируется линейная задача о собственных значениях, методы решения которой хорошо разработаны. Ради этого при выводе соотношений отбрасываются нелинейные члены и вводятся иные допущения, приводящие к погрешностям, оценка которых представляет собой непростую задачу.
В случае метода жестких границ результирующая система обычно синтезируется аналогично процедуре объединения конечных элементов с внутренними степенями свободы [56] в соответствии с методом перемещений. В случае методов свободных границ процедура синтеза несколько более сложная. В работе [53] описан алгоритм формирования уравнений задачи о собственных значениях для метода остаточных податливостей.
Однако, следует отметить, что для решения многих задач, связанных с разработкой технических систем, включающих сложную упругую конструкцию, определение в чистом виде ее собственных частот и форм колебаний является лишь промежуточной задачей. Непосредственно важными часто оказываются данные, получаемые с использованием этих характеристик, — это амплитудно-фазовые частотные характеристики, передаточные функции, импе-дансы, динамические коэффициенты влияния (динамические жесткости и динамические податливости) и т. д. Характерными примерами здесь могут быть задачи исследования устойчивости систем управления упругими объектами, задачи акустики. Все перечисленные выше величины, являющиеся функциями частоты, относятся в широком смысле к динамическим характеристикам упругой системы. Непосредственное их получение из соотношений синтезированной математической модели упругой системы без промежуточной стадии вычисления ее собственных частот и форм обеспечивает в указанных случаях решение основной задачи.
Подход, представленный в работе Крона [65] и его последователей [183, 185], фактически направлен на вычисление матриц динамических податливо-стей как нелинейных функций частоты с последующим синтезом по методу сил. Получаемая система линейных алгебраических уравнений содержит зависящие от частоты коэффициенты. Это дает возможность построения различного рода амплитудно-фазовых частотных характеристик и прочих перечисленных выше зависимостей. Определение собственных частот и форм колебаний системы здесь также возможно с использованием общих методов поиска решений нелинейных уравнений, разработаны и специализированные методы для рассматриваемого класса задач (см., например, [183, 185]).
Заметим, что достаточно точная информация о динамических коэффициентах влияния упругой системы может эффективно использоваться и для расчета переходных процессов при исследовании динамических нагрузок при помощи численных интегральных преобразований, как, например, в работе [184].
Аналогичный подход представлен в работах [164, 165], но с использованием матриц динамических жесткостей подконструкций. Он проще в реализации, поскольку основан на синтезе систем уравнений по методу перемещений, что ближе исследователям, традиционно работающим с методом конечных элементов. Здесь, как и в работах [65, 183, 185], для построения динамических матриц подконструкций используются данные об их собственных частотах и формах колебаний.
Отметим также работу [2], в которой для формирования нелинейных характеристических уравнений сложных упругих систем предлагается использовать математический аппарат метода факторизованных возмущений, предназначенного для анализа линейных физических систем с взаимодействием. Каждая налагаемая в процессе соединения подсистем связь при этом представляется возмущающим оператором, действующим на исходную совокупность не взаимодействующих подсистем. Существенной составляющей в предлагаемых математических построениях являются выражения для частотных гриновских функций изолированных подсистем. Это весьма жесткое для практического использования метода условие, поскольку построение соответствующих аналитических выражений возможно для ограниченного класса подсистем.
Общей для всех описанных подходов является проблема точности представления динамических свойств подконструкции в синтезированной системе в условиях, когда в модальном разложении может присутствовать ограниченное число собственных форм, — проблема усечения модального разложения. Этот вопрос принципиален, поскольку лежит в основе метода модального синтеза, обеспечивающего снижение размерности решаемых задач и вычислительных затрат за счет ограничения спектра подконструкций.
Очевидно, что если ставится задача исследования динамических свойств системы в диапазоне частот, ограниченных сверху некоторой частотой среза, то в модальных разложениях подконструкций должны учитываться все собственные формы, частоты которых не превосходят эту границу. Отбрасывание форм с более высокими частотами вносит погрешность в математическую модель и приводит к ошибкам в получаемых результатах. Как показывают исследования, в зависимости от способа соединения подконструкций и локальных особенностей их собственных форм влияние высших тонов на низкочастотную динамику системы меняется существенно не монотонно с возрастанием их собственных частот. Вопросам выбора критериев оценки и исследованию величины вносимой погрешности посвящены работы [171, 163, 145, 114, 52]. Предлагаемые критерии предназначены для автоматизации процесса выбора удерживаемых в модальных разложениях собственных форм подконструкций. Однако на практике наиболее употребителен подход, представленный в работах [181, 57] рекомендацией удерживать в модальных разложениях все собственные формы, частоты которых превосходят обусловленную частоту среза в 1,5 — 2 раза.
Подчеркнем, что несмотря на актуальность вопроса, ни один из рассмотренных выше методов модального синтеза не позволяет сформулировать априорную оценку погрешности получаемых результатов. Причина этого в том, что единственным варьируемым параметром, влияющим на точность получаемого результата, остается во всех методах количество удерживаемых собственных форм в модальном разложении колебаний подконструкции. Изначально не предсказуемая в общем случае сложность спектров подконструкций и структура их собственных форм не позволяет делать предварительные оценки погрешности синтеза.
В настоящей работе предлагается принципиально иной подход к оценке погрешности синтеза, основанный на отказе от идеи наращивания числа учитываемых собственных форм сверх минимально необходимого, определяемого количеством тонов с частотами, попадающими в исследуемый интервал. Повышение точности представления динамических свойств подконструкции в ограниченном частотном диапазоне достигается с помощью конструктивного алгоритма формирования вспомогательных членов модального разложения при неизменном наборе сохраняемых собственных форм.
В работе В. П. Шмакова [103] предложено строить решение задачи о гармонических колебаниях одномерной подконструкции в виде разложения по собственным формам базовой задачи, дополненного корректирующей составляющей. Эта корректирующая составляющая строится в виде многочлена относительно квадрата частоты колебаний, коэффициенты которого определяются как решения рекуррентной последовательности статических краевых задач. Поэтому степень многочлена может неограниченно наращиваться. При этом коэффициенты модальных членов разложения изменяются таким образом, что вклад тонов с высокими собственными частотами в области низких частот колебаний уменьшается с ростом порядка корректирующего многочлена. Это является основой для повышения точности получаемых при синтезе подконст-рукций решений частотного уравнения в условиях ограниченного числа учтенных собственных форм.
Аналогичный подход использован в работе А. И. Лиходеда [68] для повышения точности при расчетах нестационарного динамического нагружения упругих систем, а также при синтезе подконструкций методом жестких границ [5], где он трактуется как многократное выделение квазистатической составляющей.
Неограниченное наращивание порядка корректирующего многочлена не приводит к построению сходящегося степенного ряда. Однако, если выделить некоторое количество собственных форм и строить корректирующий многочлен в подпространстве, ортогональном к линейной оболочке этих форм, то соответствующие им (формам) коэффициенты имеют вид, не зависящий от порядка корректирующего многочлена. В таком случае получающийся степенной ряд оказывается сходящимся в ограниченном частотном интервале, если в модальное разложение включены все собственные формы, частоты которых лежат в этом интервале. С учетом этого обстоятельства автором настоящей работы был введен термин «корректирующий ряд» [37, 38, 39]. Коэффициенты этого ряда можно называть корректирующими функциями.
В настоящей работе методика синтеза динамических характеристик строится на основе гибридного подхода, когда собственные формы подконст-рукции определяются при условии частичного закрепления внешних степеней свободы. Метод жестких границ и метод свободных границ рассматриваются как частные случаи. Соединение подконструкций предполагается дискретным, т. е. интерфейс системы конечномерный.
Теоретической основой разработанного здесь метода корректирующих рядов служат две сформулированные автором основные теоремы. Суть первой из них в том, что гармонический отклик подконструкции в ограниченном частотном интервале может быть точно представлен в виде модального разложения, включающего лишь те тона колебаний, собственные частоты которых не превосходят верхней границы этого интервала, и дополненного равномерно сходящимся на этом интервале корректирующим рядом. Вторая теорема утверждает то же самое относительно ограниченного частотного интервала с ненулевой нижней границей, причем в модальном разложении должны присутствовать лишь те собственные формы, частоты которых лежат в этом интервале. В этом случае строятся корректирующие ряды относительно смещенного значения частотного параметра.
Использование этих теорем позволяет принципиально изменить идеологию модального синтеза. Усечение модальных разложений путем отбрасывания высших тонов подконструкций заменяется усечением корректирующего ряда. В процессе доказательства теорем строится алгоритм вычисления коэффициентов корректирующего ряда и, что весьма важно, выводится асимптотическая оценка погрешности его усечения. Использование этой оценки в качестве априорной при проведении пробных расчетов дает возможность реально оценивать точность получаемых результатов.
В настоящей работе принята терминология, в соответствии с которой усеченный корректирующий ряд, содержащий т членов, (фактически, корректирующий многочлен) называется корректирующим рядом т-го порядка. Упомянутая асимптотическая оценка содержит порядок корректирующего ряда в показателе степени, основание которой меньше единицы. Таким образом, наблюдается экспоненциальная сходимость решения с ростом порядка корректирующего ряда.
На основе предложенного вида модального разложения для подконст-рукции формулируются соотношения между внешними обобщенными силами и перемещениями в виде системы уравнений, содержащих набор динамических коэффициентов влияния. Эти коэффициенты могут использоваться для формирования матрицы динамических жесткостей или матрицы динамических податливостей подконструкции. Синтез системы в зависимости от этого можно осуществлять либо по методу перемещений (метод динамических жестко-стей), либо по методу сил (метод динамических податливостей).
Получаемая система уравнений может быть использована для определения собственных частот и форм упругой системы, однако наиболее удобно ее применение для непосредственного вычисления амплитудно-фазовых частотных характеристик, передаточных функций и других характеристик в решении тех задач, где такие данные используются. Следует отметить, что разработанная методика обеспечивает чрезвычайно быстрый расчет частотных характеристик системы при гарантированной высокой точности в заданном частотном диапазоне.
Основные соотношения и теоремы метода сформулированы применительно к обобщенной операторной постановке задачи для континуальных моделей упругих систем. Аналогичные результаты формулируются для дискретных моделей подконструкций, получаемых, например, методом конечных элементов. Векторные коэффициенты корректирующих рядов в этом случае называются корректирующими векторами.
Полученные формулы легко обобщаются на случай пропорционально демпфированных подконструкций. При этом предоставляется возможность задания различных коэффициентов демпфирования в подконструкциях, что существенно расширяет возможности моделирования демпфированных упругих систем.
Для проведения численных исследований и расчетов разработан программный комплекс, содержащий:
— программы расчета динамических характеристик конечноэлементных моделей упругих подконструкций с препроцессором;
— постпроцессорные программы формирования баз данных о подконст-рукциях, содержащих информацию о динамических характеристиках и корректирующих векторах;
— программы исследования спектра составной упругой системы и построения ее частотных характеристик с учетом демпфирования.
Сформированные предварительно базы данных о подконструкциях используются в алгоритме синтеза для определения динамических и частотных характеристик составной системы. Замена или модификация какой-либо компоненты системы отражается лишь в изменении соответствующей ей базы данных. Разработана универсальная структура баз данных для конечноэле-ментных моделей подконструкций.
Программный комплекс строится на основе общих принципов и может дополняться программами расчета необходимых для синтеза данных о под-конструкциях различного типа. В настоящее время разработаны версии для пространственных стержневых конструкций и осесимметричных оболочечных конструкций, содержащих жидкость. Могут использоваться и аналитические соотношения для континуальных моделей подконструкций. В частности, выведены формулы для изгибных колебаний однородных стержней, соединяемых в концевых сечениях.
В процессе исследований автором настоящей работы обнаружено, что рекуррентному алгоритму построения корректирующих векторов присуще свойство неустойчивости, приводящее к быстрому накоплению погрешностей и распаду решения уже при учете 5−6 членов корректирующего ряда [37]. Эффективным средством подавления этой неустойчивости оказалась ортогонали-зация вычисляемых векторов на каждом шаге рекуррентного процесса к учтенным собственным формам подконструкции. Это обстоятельство свидетельствует о том, что формулировки работ [103, 104, 105, 5] не достаточны для построения устойчивого гарантирующего точность результата алгоритма. Построение корректирующих функций неизбежно должно осуществляться в ортогональном к учтенным собственным формам подпространстве (с коррекцией ортогональности на каждом шаге).
Существенные затруднения при реализации модального синтеза вызывает наличие среди собственных частот подконструкции нулевых значений, соответствующих формам движения твердого тела. Предложенный автором способ использования частотного сдвига при вычислении корректирующих векторов позволяет избежать трудностей и не вносить в алгоритм принципиальных изменений [39].
Помимо гарантированной оценки точности результата синтеза для метода корректирующих рядов важна оценка дополнительных затрат, связанных с вычислением корректирующих векторов, в сравнении с расширением набора собственных форм подконструкции. Учитывая то обстоятельство, что факторизация матрицы для решения последовательности статических задач выполняется однократно, можно приблизительно оценить затраты на вычисление одного дополнительного корректирующего вектора как на порядок меньшие по сравнению с затратами на вычисление дополнительной собственной формы. При этом проведенные исследования показали, что при использовании 5−10 членов корректирующего ряда относительная погрешность результатов синтеза не превышает 10−5-10−6. Это свидетельствует о неоспоримом преимуществе предложенного метода, особенно существенном для подконструкций с высокой плотностью спектра (что характерно для сложных пространственных конструкций). Его эффективность демонстрируется в настоящей работе на примере сложной составной конструкции орбитальной космической станции, модифицируемой в процессе сборки и изменяющей свою конфигурацию в процессе функционирования (например, в связи со стыковкой с транспортными кораблями или добавлением новых модулей).
Наличие среди компонентов сложной упругой системы гидроупругих звеньев, представляющих собой упругие конструкции, взаимодействующие с ограниченными объемами жидкости (в частности, содержащие жидкость во внутренних полостях), существенно усложняет задачу исследования ее динамических свойств. Такие задачи актуальны для исследования динамики летательных аппаратов, гидротехнических сооружений, проблем транспортировки жидких грузов. Успех их решения непосредственно зависит от правильности математической постановки задачи и эффективности методов определения динамических характеристик упругих конструкций с жидкостью, входящих в состав системы.
Как правило, конструкции, содержащие жидкость, представляют собой тонкостенные сосуды, динамическое поведение которых хорошо описывается с помощью соотношений теории оболочек.
Первые результаты в области динамики колебаний упругих оболочек, взаимодействующих с жидкостью, связаны с именами Рэлея, Жуковского Н. Е. и других исследователей конца XIX — начала XX веков. Однако период наиболее интенсивной разработки данной проблемы относится ко второй половине XX века в связи с развитием ракетной и авиационной техники, а также электронных вычислительных систем.
В основе современного этапа развития методов решения данного класса задач лежат работы таких отечественных исследователей как Моисеев Н. Н. [73], Рабинович Б. И. [92], Шмаков В .П. [102], Рапопорт И. М. [94]. Среди зарубежных специалистов здесь можно отметить работы Abramson H.N., Kana D.D., Lindholm U.S., Bauer H.F. Значительный вклад в решение проблемы внесли исследования Григолюка Э. И., Шклярчука Ф. Н., Горшкова А. Г., Балабуха Л. И., Балакирева Ю. Г., Лампера Р. Е., Пожалостина А. А. и других.
Вообще, число опубликованных к настоящему времени работ, посвященных задаче расчета динамических характеристик оболочек, содержащих жидкость, весьма велико. Основную часть из них составляют работы, в которых исследуются оболочки определенной формы и решения получаются, как правило, с помощью какого-либо вариационного метода, причем выбор координатных функций определяется формой полости. В этих работах получены приближенные или точные формулы для ряда оболочек простой геометрической формы. Однако для практики, где сложность конструкторских решений часто затрудняет получение аналитических оценок и не всегда позволяет использовать простые модели, наибольшую ценность представляют универсальные численные методы, не налагающие жестких ограничений на форму и параметры исследуемых конструкций. Здесь следует отметить разработанный под руководством В. П. Шмакова метод расчета динамических характеристик оболочек вращения с жидкостью [17, 18, 23], основанный на разложении потенциала смещений жидкости в ряд по собственным функциям гидродинамической задачи (решаемой методом Ритца) и использовании при решении уравнений теории оболочек метода ортогональной прогонки. Отметим также предложенный Р. Е. Лампером метод Ритца с варьируемым параметром [3, 67], позволяющий рассчитывать динамические характеристики широкого класса осе-симметричных баков с жидкостью. Применялись при решении указанной задачи и прямые численные методы: метод конечных разностей [6, 7] и метод суммарных представлений [8, 9], также основанный на конечно-разностной аппроксимации дифференциальных операторов.
В настоящее время известны решения задач для весьма широкого класса оболочек с жидкостью: исследовались цилиндрические оболочки с различными днищами (плоскими, сферическими, коническими), сферические оболочки, конические оболочки, соосные цилиндрические оболочки. Мы не указываем здесь эти работы ввиду большого их числа — библиография по этому вопросу содержит несколько сотен названий. Заметим однако, что, несмотря на это, ряд вопросов долгое время оставался не решенным. Неизвестны, например, работы, в которых изучались бы колебания тороидальной оболочки с жидкостью. Лишь единичные работы касаются задачи о колебаниях систем из двух баков с промежуточными разделительными днищами, причем рассмотрены весьма частные случаи. При этом на практике, например, в изделиях ракетной техники часто встречаются баки весьма сложной конфигурации, с двусвязными полостями, с разделительными днищами, когда невозможно рассматривать отдельно колебания жидкости в баках горючего и окислителя. Значительно усложняют разработку простых моделей для аналитических оценок и приближенных расчетов такие факторы, как переменность толщины оболочки, наличие шпангоутов, подкрепляющих ребер.
Все это вызвало потребность разработки простого и универсального алгоритма, позволяющего быстро выполнить расчет с учетом как можно большего числа конструктивных особенностей (возможно, ценой увеличения затрат времени работы вычислительной системы). В наибольшей степени таким условиям при расчете упругих конструкций отвечает метод конечных элементов [56]. В настоящей работе развиваются теоретические и методические основы его применения для расчета динамических характеристик оболочечных конструкций с жидкостью, функционирующих самостоятельно либо входящих в состав сложной упругой системы в качестве подконструкций.
Применению метода конечных элементов к расчету динамики упругих конструкций, взаимодействующих с жидкостью, посвящен целый ряд работ преимущественно зарубежных авторов. Если попытаться их классифицировать, то можно выделить две группы работ.
К первой группе отнесем те работы, в которых движение частиц жидкости описывается непосредственно векторами их смещений. В [147, 148] разработаны дискретные элементы, моделирующие свойства частиц сжимаемой или несжимаемой жидкости, которые могут быть включены в конечноэлементную модель упругой конструкции, содержащей жидкость. В комментарии к этим работам [124] отмечается, что эти элементы построены лишь на основе физических соображений и интуиции, и предлагается выражение для лагранжиана, на основе которого можно строить любые конечные элементы жидкости, причем несжимаемость может быть обеспечена наложением некоторых ограничений на степени свободы элемента. Фактически такое выражение для лагранжиана и использовано в работе [161], где описаны конечные элементы для аппроксимации сжимаемой жидкости. К этой же группе работ относится и работа [177], где метод конечных элементов разработан для расчета неосесиммет-ричных колебаний оболочек вращения с несжимаемой жидкостью. Отметим, что подход, использованный в [177], отличается чрезмерной усложненностью и громоздкостью построений, не свойственными, вообще говоря, методу конечных элементов.
В работе [97] обсуждается применение конечных элементов, предложенных в [161], к расчету осесимметричных колебаний оболочек вращения с жидкостью. Результаты сопоставляются с результатами расчета при помощи элемента с матрицей жесткости, сформированной на основе принципа минимума дополнительной энергии. При этом точность, полученная в первом случае, оказывается значительно хуже, чем во втором.
Если пренебречь влиянием поверхностных волн, то эквивалентным подходом является интерпретация сжимаемой жидкости в виде упругой (акустической) среды с нулевым модулем сдвига [110, 111]. Такой подход позволяет использовать с небольшими дополнениями имеющиеся конечноэлементные программы общего назначения.
Методы этой первой группы содержат ряд общих сложностей. Во-первых, это обеспечение условий на поверхности контакта упругого тела с жидкостью, для выполнения которых необходимо либо налагать дополнительные ограничения, как, например, в [177], либо вводить специальные конечные элементы, как предлагается в [110]. Во-вторых, как отмечается в [161], появляется набор собственных форм с нулевой частотой, соответствующих циркуляционным течениям жидкости, что усложняет решение задачи о собственных значениях.
Такие же циркуляционные формы получаются в работе [143], где для описания движения жидкости использованы уравнения в перемещениях в ла-гранжевой форме и сформулирован вариационный принцип, используемый для построения конечных элементов. Для исключения этих форм в этой работе предложено использовать метод штрафных функций.
Ко второй группе отнесем работы, в которых используется предположение о безвихревом движении жидкости, что позволяет описывать его при помощи одной скалярной функции (давление, потенциал скоростей, потенциал смещений). Это не только исключает формы циркуляционного течения жидкости, но и снижает размерность дискретизированной задачи.
В работах [187, 189, 160] поведение несжимаемой жидкости описывается потенциалом скоростей, который с помощью интегральных соотношений теории потенциала исключается из уравнений движения конструкции. Такой подход (часто используемый при решении задач гидроупругости) позволяет ограничиться описанием движения лишь оболочки и свободной поверхности жидкости, и соответственно лишь их дискретизацией при использовании метода конечных элементов, что особенно удобно при решении трехмерных задач. Однако при этом матрица масс конечноэлементной модели теряет весьма существенное свойство, а именно, ленточную структуру. Это резко увеличивает требуемые размеры памяти, используемой при расчетах вычислительной системы. Такой же результат получается в [56, 123, 112, 54], а также в [125] для несжимаемой жидкости, где сначала выполняется конечноэлементное разбиение жидкости, а затем из дискретной системы уравнений исключаются переменные, соответствующие жидкостным потенциалам (или давлению). Разбиение жидкости на конечные элементы с сохранением ленточной структуры матриц, как показано в [161], позволяет более эффективно использовать машинную память.
Вопрос о формировании матриц дискретизированной системы в случае, когда разбиение осуществляется как для упругой конструкции, так и для жидкости, может решаться несколькими способами. В работах [56, 193, 112, 125], где поведение несжимаемой жидкости описывается давлением, отдельно строятся системы дискретных уравнений для упругой конструкции и для жидкости. При построении же матриц, определяющих взаимосвязь этих двух систем, давление жидкости и перемещения смоченной поверхности рассматриваются как внешние воздействия по отношению соответственно к первой и ко второй из них. Получающиеся системы уравнений несимметричны, а методы их симметризации, как и при исключении соответствующих жидкости переменных, нарушают ленточную структуру матриц. В работах [138, 191, 192] предложены методы решения получаемых при таком подходе задач о собственных значениях, экономно использующие память вычислительной системы.
В работах [118, 119, 120, 121, 188] движение несжимаемой жидкости описывается при помощи потенциала скоростей. Для построения конечноэле-ментных систем уравнений здесь все уравнения (в том числе и соотношения на поверхности контакта конструкции с жидкостью) записываются в вариационной форме. Дискретизация задачи при таком подходе соответствует методу Бубнова-Галеркина. Для учета условий контакта с жидкостью при этом, вообще говоря, необходимо вводить специальные конечные элементы.
Однако наиболее простым и естественным является подход, основанный на дискретизации функционала действия, условие стационарности которого удовлетворяется решением задачи. Такой подход использован в работе [54], где решается задача о колебаниях пластинки на поверхности несжимаемой жидкости, заполняющей жесткий резервуар, а также в работе [170], являющейся развитием работы [118] в плане учета сжимаемости жидкости. При этом из условия стационарности функционала следуют все соотношения, в том числе и граничные условия контакта упругого тела и жидкости, что приводит к их автоматическому учету в дискретизированных системах уравнений. Это обеспечивается надлежащим выбором вида функционала. В работах [32, 33] такой подход называется смешанным вариационным принципом. Фактически он эквивалентен используемому в [118] в смысле построения матриц конечноэле-ментной модели.
В работах [170, 175, 176] для описания поведения сжимаемой жидкости используются две функции: потенциал смещений и давление.
Введение
второй независимой переменной при учете сжимаемости оказывается необходимым для получения симметричных матричных уравнений. В противном случае, как, например, в [125, 133], где движение жидкости описывается лишь давлением, приходится осуществлять дополнительные преобразования с целью симметризации системы.
Для получения правильного результата при использовании численного метода, в частности, метода конечных элементов, первоочередное значение имеет адекватность исходной математической модели (системы уравнений, вариационного принципа) исследуемому физическому явлению. Возможность пренебрежения в математической модели влиянием отдельных факторов с целью упрощения процедуры решения должна быть не только обоснованной, но и соответствующие преобразования должны быть выполнены корректно, не приводя к противоречиям с физическими законами. Поэтому целая глава настоящей работы посвящена формулировке математической модели колебаний упругой конструкции с жидкостью в условиях однородного гравитационного поля.
В настоящей работе рассматриваются малые колебания упругой конструкции, содержащей жидкость во внутренних полостях (либо взаимодействующей с ограниченным объемом жидкости). Колебания жидкости также малы, и поэтому задача может рассматриваться в линейной постановке. Линеаризация исходных, вообще говоря, нелинейных соотношений должна осуществляться в окрестности статического равновесного состояния системы. Тем самым при определении ее динамических характеристик учитывается влияние статического нагружения. Материал конструкции считается линейно упругим, жидкость — идеальной и несжимаемой. Конструкция находится под действием однородного гравитационного поля, определяющего положение недеформиро-ванной свободной поверхности жидкости.
Традиционно при формулировке замкнутой системы уравнений, описывающих колебания упругой конструкции, содержащей во внутренней полости жидкость, исследователи основывались на совмещении двух различных способов описания движения сплошной среды. Если уравнения малых колебаний упругого тела представляют лагранжев подход, то для описания движения жидкости наиболее распространен метод Эйлера. При этом используется предположение о потенциальности течения жидкости, позволяющее существенно упростить соответствующие уравнения. Такое совмещение разнородных методов повышает вероятность появления в постановке задачи некорректностей, одна из которых рассмотрена ниже. (Хотя, конечно, корректное построение определяющих соотношений возможно и на основе смешанного эйлерово-лагранжева подхода.).
В настоящей работе вывод уравнений осуществлен на основе лагранжева подхода как для упругой конструкции, так и для заполняющей ее полости жидкости. Техника последовательного единого лагранжева подхода в проблемах динамического взаимодействия разнородных сред изложена в работах [58, 96, 59], где основное внимание уделено вопросам нелинейности получаемых соотношений. Применительно к случаю малых колебаний контактные соотношения линеаризованы в окрестности стационарных значений переменных параметров. В предположении о безвихревом движении жидкости для его описания введен потенциал смещений. Выведена формула для лагранжиана консервативной системы «упругая конструкция — жидкость». Сформулирован вариационный принцип смешанного типа, среди условий стационарности которого содержатся контактные соотношения на смоченной поверхности упругой конструкции, что важно для эффективной реализации метода конечных элементов.
Важным моментом в описанных построениях является корректный учет влияния гравитационного поля на динамику системы.
Влияние интенсивности однородного гравитационного поля на собственные частоты колебаний упругих конструкций, содержащих жидкость, может осуществляться через два различных механизма взаимодействия этого поля с конструкцией. Первый из них реализуется опосредованно через статическую составляющую напряженно-деформированного состояния конструкции. Математически учет этого эффекта сводится к линеаризации уравнений малых колебаний относительно этой статической составляющей. Если материал конструкции линейно упругий, то определяющую роль здесь играет геометрическая нелинейность в выражении для тензора деформаций.
Второй канал влияния связан с работой, выполняемой гравитационными силами, действующими на жидкость, в процессе колебаний гидроупругой системы. Одно из проявлений этой работы связано с волнообразованием на свободной поверхности жидкости. Однако этим механизм явления не исчерпывается.
Тесно связанным с волнообразованием оказывается эффект колебаний гидростатического давления жидкости на стенки полости вследствие их движения. Локальное вертикальное смещение контактирующей с жидкостью стенки сосуда приводит к изменению давления на нее вследствие изменения глубины. Появляющийся при учете этого фактора член в уравнениях малых колебаний имеет первый порядок малости, т. е. по такому формальному признаку отброшен быть не может и должен присутствовать в динамических соотношениях на контактной поверхности.
Впервые корректное выражение для динамических условий на контактной поверхности приведено применительно к оболочкам вращения в работе Э. И. Григолюка и Ф. Н. Шклярчука [33]. При их выводе принципиален учет того обстоятельства, что гидростатическая сила является «следящей». Это выражается в том, что, во-первых, направление действия силы меняется в соответствии с поворотом нормали к деформированной смоченной поверхности, а во-вторых, изменение площади элемента поверхности вследствие растяжения или сжатия в процессе деформации влияет на величину совершаемой этой силой работы. Эти два обстоятельства отмечены в работе [20] применительно к проблемам устойчивости упругих оболочек. Только при учете этих факторов удается записать выражение для потенциальной энергии гравитационых сил, действующих на жидкость, и следовательно, выражение для лагранжиана гидроупругой системы.
В процессе исследований выяснилось, что выражение для потенциальной энергии гравитационных сил жидкости, приведенное в работе [33], верно только для гладких поверхностей контакта жидкости с конструкцией. В настоящей работе получено выражение, применимое для произвольного случая кусочно-гладких смоченных поверхностей.
Заметим, что в большинстве из рассмотренных выше работ гравитационными эффектами полностью пренебрегается. Это допустимо в тех случаях, когда для исследования представляют интерес лишь колебания, обусловленные преимущественно упругостью стенок сосуда, а их спектр существенно отделен от собственных частот колебаний свободной поверхности жидкости. В этом случае на свободной поверхности задается условие равенства нулю давления или потенциала смещений (или скоростей). Естественно, отсутствует гравитационная составляющая и в контактных соотношениях. Такой подход можно считать предельным случаем упрощения гидроупругой задачи. Он весьма распространен, поскольку его результаты дают точность, достаточную для широкого круга приложений.
Однако, поскольку вносимая гравитационным членом поправка в эффективную жесткость стенок сосуда, как правило, весьма мала, возникает вопрос о возможности пренебрежения ее величиной в уравнениях математической модели при сохранении членов, описывающих волнообразование на свободной поверхности.
В настоящей работе показано, что простое вычеркивание соответствующего члена в динамических контактных соотношениях недопустимо, поскольку приводит к противоречивой с точки зрения механики математической модели. А именно, система получает дополнительную связь, препятствующую ее смещению как жесткого тела в направлении вектора гравитации. Оценки показывают, что вносимая этой связью погрешность в расчетный спектр системы может быть значительной.
Это свидетельствует о взаимной связи колебаний свободной поверхности и деформации контактной поверхности как единой системы, состояние которой определяет потенциальную энергию гравитационных сил жидкости.
В работе показано также, что, варьируя способ определения аддитивной константы в потенциале смещений, можно получить эквивалентную систему уравнений, в которой вычеркивание гравитационного члена в динамических контактных соотношениях не приводит к противоречиям. Фактически, в результате такого преобразования свободная поверхность жидкости выделяется в подсистему с собственной потенциальной энергией гравитационных сил. Вследствие указанного выбора константы в уравнениях колебаний свободной поверхности жидкости и в лагранжиане системы появляется интегральная величина — среднее смещение свободной поверхности.
Введение
такой интегральной величины вызывает затруднения при реализации метода конечных элементов, обесценивающие это «упрощение».
Необходимо отметить, что в перечисленных выше основополагающих работах отечественных исследователей [73, 92, 102, 62, 94], где влияние гравитации учитывалось только в форме поверхностных волн, в уравнениях также присутствуют указанные интегральные члены, обеспечивающие непротиворечивость математической модели.
Если же говорить о работах зарубежных специалистов, то библиографические исследования показывают, что все они (за исключением тех, где гравитационные эффекты не учитываются совсем) содержат ту или иную форму некорректности. Единственным обнаруженным исключением являются работы [108, 109], где исследуются математические свойства краевых задач, сформулированных на основе результатов работ [73, 94].
Несмотря на относительную давность результата работы [33], он остался, по-видимому, полностью вне сферы внимания зарубежных специалистов. Некорректный подход к формулировке динамических условий контакта упругого тела с жидкостью, когда не учитываются две отмеченные выше особенности действия гидростатической силы в условиях колебаний, приводит к появлению несимметричных операторов.
В качестве примера можно привести работу [118], где с целью достижения симметричности введена совершенно не оправданная с точки зрения механики гипотеза, состоящая в пренебрежении касательными к поверхности контакта смещениями. Результатом явились уравнения, не инвариантные относительно сдвига системы и, как следствие, дополнительно наложенные на систему связи.
Связи эти действуют не только в вертикальном направлении, как в рассмотренном выше случае неверных уравнений для свободной поверхности при отбрасывании гравитационной компоненты в динамических контактных соотношениях, но и в горизонтальном направлении. При этом жесткости этих связей могут быть как положительными, так и отрицательными. В работе [32] отмечено, что такой подход допустим, если нормальные перемещения смоченной поверхности значительно больше тангенциальных. Однако в случае систем, колебания которых характеризуются значительными средними смещениями содержащихся в них масс жидкости, указанные связи должны оказывать существенное влияние на результаты расчетов. (Отметим единственный геометрический вариант, когда эти связи не возникают, — полость, образованная строго вертикальными стенками и горизонтальным плоским днищем.).
Часто такие же формулы получаются из-за небрежности при выводе динамических контактных условий, выполняемом с помощью наглядных геометрических представлений, как это сделано в работе [133]. На таких же формулах основано выражение для энергии сил тяжести конечного элемента в работах [161, 124]. Такое же выражение использовано в работе [160] для потенциальной энергии объема жидкости. На основе таких же формул построены конечные элементы жидкости в распространенном программном комплексе АКБУБ версии 5.4. Ссылки на эти же формулы содержатся в работах отечественных авторов [55, 74].
Возможность получить удобные формулировки, связав гравитационную потенциальную энергию жидкого объема (и конечного элемента жидкости) лишь с параметрами его собственной деформации, отвлекает разработчиков от механической сущности проблемы. Форма занимаемого жидкостью объема определяется удерживающими его стенками сосуда, а работа гравитационных сил жидкости осуществляется на перемещениях этих стенок (а не контактирующих с ними частиц жидкости в условиях взаимного скольжения). Отсюда следует невозможность конструирования жидкостного конечного элемента, в котором учитывался бы этот фактор.
Иного рода усилия были приложены разработчиками известного конеч-ноэлементного комплекса иА1/ЫА8ТКАК [186], где в результате применения метода конечных элементов к системе с несимметричными операторами получаются уравнения с несимметричными матрицами. Для решения проблемы собственных значений в состав комплекса включены специально разработанные алгоритмы, требующие нетривиальных ресурсов вычислительных систем. Вопрос о точности получаемых результатов требует отдельного исследования.
Достоин удивления факт значительного количества работ, в которых при не учете влияния гравитации на контактные соотношения уравнения поверхностных волн сформулированы неверно — без учета среднего смещения свободной поверхности. В качестве примеров укажем работы [193, 175, 176, 122]. С попытками устранить возникающие при этом противоречия связано, по-видимому, появление работ [134, 172, 173], где введены усложненные формулировки, основанные на труднообъяснимом в рамках линейной задачи введении понятия «отсчетного» состояния для учета влияния движения стенок сосуда.
Широкий класс практических задач связан с определением динамических характеристик осесимметричных оболочечных конструкций, полости которых частично заполнены жидкостью. Вектор гравитационных сил в этом случае коллинеарен продольной оси конструкции. В настоящей работе задача расчета динамических характеристик конструкций такого типа рассмотрена подробно и разработан алгоритм метода конечных элементов для ее решения.
Осесимметричность конструкции позволяет в цилиндрической системе координат свести пространственную задачу определения динамических характеристик к двумерной путем разложения неизвестных в ряд Фурье по окружной координате. Расщепление системы уравнений для гармоник на независимые подсистемы приводит к тому, что собственные формы представляют собой колебания с целым числом волн по окружности (в случае нулевого значения это продольно-радиальные и крутильные колебания).
Рассмотренная выше постановка задачи о колебаниях конструкций, содержащих жидкость, переработана с учетом гипотез теории тонких упругих оболочек и нелинейности геометрических соотношений [98, 77]. Линеаризация в окрестности статического напряженно-деформированного состояния, возникающего в результате действия давления газа в полостях конструкции и сил тяжести, позволяет учесть влияние этих факторов на собственные частоты и формы колебаний. Рассмотрен общий случай произвольного количества несвязанных между собой полостей, частично заполненных жидкостями разной плотности. Допускается наличие промежуточных между этими полостями стенок, касающихся жидкости с обеих сторон. Учитываются такие конструктивные особенности, как наличие упругих силовых шпангоутов, а также несимметричного относительно срединной поверхности оболочки подкрепления в виде часто расположенных продольных и кольцевых ребер.
Жидкость считается идеальной и несжимаемой, движение ее безвихревое и описывается при помощи потенциала смещений. Для применения метода конечных элементов использован описанный выше вариационный принцип в виде условия стационарности функционала, из которого следуют уравнение несжимаемости, граничные условия на свободной поверхности жидкости, а также условия на контактной поверхности жидкости и оболочки с учетом влияния однородного гравитационного поля.
В работе описаны конечные элементы для дискретизации жидкости (включая элементы свободной поверхности) и оболочечные конечные элементы, описаны способы формирования их матриц и способ учета влияния статического напряженно-деформированного состояния на матрицу жесткостей оболочечного элемента. Заметим, что оболочечные элементы сформулированы таким образом, что учтена возможность контакта с жидкостью, и это освободило от необходимости введения специальных контактных элементов.
На основе этих конечных элементов разработан программный комплекс, предназначенный для расчета динамических характеристик оболочечных конструкций с жидкостью при осесимметричных колебаниях, при колебаниях с заданным числом волн по окружности, с учетом влияния внутреннего давления и собственного веса конструкции.
Этот программный комплекс расширяет возможности разработанных ранее программ [34, 35], эксплуатировавшихся на вычислительных машинах БЭСМ-6, единой серии ЕС-1050, ЕС-1060, на компьютерах типа IBM PC. Программа расчета динамических характеристик оболочечных конструкций с жидкостью при осесимметричных колебаниях включена в ФАП по ракетно-космической технике РКА [36], использовалась во всех конструкторских бюро ракетно-космической отрасли в расчетах динамических характеристик топливных баков и корпусов жидкостных ракет.
Благодаря тщательной отладке и тестированию, а также использованию эффективных численных методов достигнута высокая надежность алгоритма, не имевшего отказов при проведении тысяч расчетов. Решающую роль в распространении программ сыграла разработка удобных препроцессорных и постпроцессорных процедур, а также возможность вести расчет практически от чертежа, без трудно обосновываемых упрощений и введения промежуточных расчетных моделей с использованием механических аналогов.
Новый программный комплекс позволяет (в дополнение к прежним возможностям) в полном объеме учитывать гравитационные эффекты. Это важно ввиду того, что на практике далеко не всегда реализуются условия, позволяющие пренебрегать влиянием волн на поверхности жидкости. В работе приведен также пример конструкции, расчет динамических характеристик которой в принципе невозможен без учета гравитационной составляющей соотношений на смоченной поверхности.
Применительно к проблемам динамики жидкостных ракет важное преимущество программного комплекса связано с возможностью определять динамические характеристики корпуса целиком, без промежуточной схематизации его стержнем с упруго подвешенными массами (что по сути является способом синтеза динамических характеристик). Такая схематизация затруднительна в упомянутом выше случае вложенных друг в друга баков с промежуточным разделительным днищем, когда невозможно рассматривать колебания жидкости в баках по отдельности. Кроме того, заметное влияние на собственные формы корпуса оказывает связь между радиальными и продольными колебаниями стенок несущих баков за счет отличия от нуля коэффициента Пуассона (не учитываемая в стержневой модели корпуса), что показано в работе [22] на примере продольных колебаний конструкции, составленной из отрезков стержня и цилиндрических баков, частично заполненных жидкостью.
В программном комплексе предусмотрена также возможность модального синтеза подконструкций (отдельных отсеков) с использованием метода корректирующих рядов, обеспечивающего гарантированную высокую точность результатов, в случае когда расчет полной модели корпуса затруднен из-за высокой размерности задачи. В работе проведено сравнительное исследование точности результатов, получаемых различными методами синтеза.
Возбуждение динамических процессов в упругой или гидроупругой системе в составе сложного технического объекта связано, как правило, с работой источника энергии. Интенсивность его воздействия на конструкцию связана с изменением параметров ее деформации, причем в этом процессе наряду с собственными динамическими характеристиками источника энергии значительную роль могут играть и специальные регулирующие устройства. Такая система образует замкнутый контур, в котором могут развиваться автоколебания, отрицательно влияющие на условия функционирования объекта.
Одна из актуальных проблем такого рода связана с исследованием продольных автоколебаний жидкостной ракеты в полете, возникающих в результате взаимодействия колебаний корпуса ракеты и вызываемых ими колебаний жидкости в топливных магистралях с колебаниями давления в камере сгорания и тяги двигателя.
Это явление достаточно хорошо изучено с точки зрения исследования условий возникновения автоколебаний при неустойчивости в замкнутом контуре «корпус — топливная магистраль — двигатель». Здесь значительный вклад внесли отечественные ученые Микишев Г. Н. [71], Рабинович Б. И. [91], Шмаков В. П. [93], Колесников К. С. [63, 64], Натанзон М. С. [75, 76], Пилипенко В. В. [88], Балакирев Ю. Г. [13, 16], а также ряд американских исследователей [178, 180, 190].
Тем не менее, в значительной мере дискуссионным остается вопрос о достигаемых при продольной неустойчивости амплитудах колебаний, которые могут оказаться невелики, не требуя поэтому специальных конструктивных мероприятий по их подавлению.
В любом из звеньев указанного замкнутого контура могут проявляться нелинейные зависимости, влияющие на предельные значения амплитуд автоколебаний. Влияние нелинейных характеристик двигательной установки, обусловленных кавитационными явлениями в насосах, на амплитуды продольных колебаний жидкостных ракет рассматривалось в работах [88, 87]. Динамика корпуса описывалась при этом линейными соотношениями. В результате расчетов получено качественное соответствие форм колебаний параметров системы типичным данным летных испытаний.
Тем не менее, упругая конструкция корпуса также может играть значительную роль в ограничении амплитуд. Предпосылкой для проведения исследований в этом направлении послужило экспериментально установленное явление возбуждения неосесимметричных форм колебаний осесимметричных оболочек с жидкостью при продольном воздействии, существенно снижающее порог, при котором проявляются их нелинейные свойства [158, 1]. Механизм этого явления имеет характер параметрического возбуждения при действии вызванных продольным воздействием периодических колебаний осесиммет-ричного напряженного состояния оболочки. Поэтому при исследовании влияния нелинейности корпуса на динамику продольных автоколебаний нельзя ограничиться рассмотрением лишь осесимметричных форм колебаний оболочек, как это делается при исследовании устойчивости.
Мы рассмотрим в данной работе ракеты тандемной компоновки, корпуса которых с хорошей степенью точности можно рассматривать как осесиммет-ричные оболочечные конструкции с жидкостью.
Параметрические колебания оболочек с жидкостью привлекают внимание исследователей с начала шестидесятых годов. Важность изучения этого эффекта особо отмечена в обзорном докладе Э. И. Григолюка [32] на VII Всесоюзной конференции по теории оболочек и пластинок в 1969 году.
По-видимому, статья [26] - первая из работ, где теоретически рассмотрено это явление. В ней показано, что периодические осесимметричные усилия в оболочке играют роль параметрических нагрузок для неосесимметричных форм колебаний и в качестве примера рассмотрена цилиндрическая оболочка с шарнирно опертыми краями, заполненная жидкостью и находящаяся под действием переменных во времени продольной сжимающей силы, распределенной продольной нагрузки и гидростатической силы. Задача о динамической неустойчивости методом Бубнова-Галеркина сводится к системе уравнений типа Хилла. Однако усилия, характеризующие невозмущенное напряженно-деформированное состояние, определяются без учета инерции оболочки и динамического взаимодействия оболочки с содержащейся в ней жидкостью, что ограничивает область применимости полученных соотношений.
Аналогичный подход использован в работах [10, 11] при исследовании параметрических колебаний цилиндрической оболочки, заполненной жидкостью медленно меняющейся глубины, под действием периодической продольной сжимающей силы, а также в работе [49] при исследовании устойчивости колебаний цилиндрической оболочки с пологим сферическим дном и массами на торцах, частично заполненной жидкостью, под действием приложенных к торцам периодических продольных сил.
В работах [99, 100], где исследовались параметрические колебания цилиндрической оболочки при продольных воздействиях, осесимметричные усилия в оболочке определяются из решения задачи об осесимметричных колебаниях оболочки совместно с частично заполняющей ее жидкостью. Из соотношений нелинейной теории пологих цилиндрических оболочек в работах [80, 81] выведены уравнения нелинейных параметрических осесимметричных колебаний цилиндрической оболочки, содержащей жидкость. Получены условия возбуждения неосесимметричного тона и выражения для амплитуды колебаний в окрестности главного параметрического резонанса. Построение областей неустойчивости осуществлялось в данных работах приближенно с помощью метода гармонического баланса.
Аналогичные результаты, но в менее точной постановке, были получены в работе [158]. Исследование различного вида параметрических резонансов цилиндрической оболочки с жидкостью при гармонических колебаниях ее основания выполнено в работах [82, 83, 84].
Необходимо отметить большое значение экспериментальных исследований параметрических колебаний цилиндрической оболочки с жидкостью, описанных в работах [158, 1]. Экспериментальные данные содержатся также в работе [69].
В настоящей работе описана разработанная автором методика оценки амплитуды автоколебаний системы, состоящей из осесимметричной оболочеч-ной конструкции с жидкостью, гидравлически связанного с ней трубопровода и регулятора, управляющего величиной действующей на конструкцию продольной силы зависимости от давления на выходе из трубопровода. Эта методика позволяет учитывать геометрическую нелинейность поведения оболочек и связанный с ней эффект параметрического возбуждения неосесимметричных форм колебаний.
Первым этапом методики является расчет динамических характеристик оболочечной конструкции с жидкостью с помощью упомянутого выше программного комплекса, основанного на методе конечных элементов. Вычисляются как осесимметричные, так и неосесимметричные собственные формы, частоты которых лежат в представляющем интерес диапазоне.
Второй этап заключается в вычислении коэффициентов при нелинейных членах (квадратичных и кубичных) в уравнениях колебаний конструкции, записанных в нормальных координатах. Эти коэффициенты представлены интегралами по поверхности оболочки и вычисляются при помощи той же конеч-ноэлементной сетки, которая использовалась на первом этапе. В работе проанализирована структура определяемых этими коэффициентами нелинейных связей между осесимметричными и неосесимметричными тонами в зависимости от числа волн и сдвига по окружности (в случае кратных частот). Сформулированные условия неравенства нулю нелинейных коэффициентов значительно снижают объем вычислительной работы.
На третьем этапе осуществляется расчет областей параметрического возбуждения неосесимметричных форм колебаний конструкции при гармоническом продольном силовом воздействии для оценки степени возбудимости этих форм (и выделения наиболее важных для учета на следующем этапе). Границы областей строятся на плоскости «частота — амплитуда воздействия». Для этого с использованием данных первых двух этапов формируются уравнения в вариациях относительно стационарного осесимметричного отклика конструкции, образующие систему линейных однородных дифференциальных уравнений с периодическими коэффициентами. В предположении, что амплитуда воздействия достаточно мала, параметры этого осесимметричного отклика определяются аналитически из линейных уравнений (этому условию соответствуют нижние границы областей параметрического возбуждения неосесимметричных тонов). Для построения областей используется метод, основанный на вычислении мультипликаторов системы, предложенный в работе [21] и дающий точное положение границ.
И последний, четвертый, этап состоит в исследовании амплитуд колебаний, развивающихся в автоколебательной системе. Для его реализации в работе выведена система уравнений и сформулирована задача Коши, описывающая переходный процесс. В уравнениях учтено взаимное влияние колебаний жидкости в баках и в присоединенных к ним трубопроводах. Динамика источника энергии (двигателя) описывается интегро-дифференциальным уравнением типа Вольтерра, ядро которого вычисляется на основе данных о его амплитудно-фазовых частотных характеристиках.
Сложный вид нелинейности, содержащейся в этой системе, ограничивает возможность ее аналитического исследования. Поэтому исследование проводится методом численного интегрирования при различных начальных условиях до выхода колебательного процесса на максимальные значения амплитуд. Для решения систем интегро-дифференциальных уравнений общего вида была разработана и отлажена подпрограмма, на основе которой составлена программа численного решения уравнений продольных колебаний.
Проведенное численное исследование показало, что, во-первых, учет геометрической нелинейности оболочек в рамках предположения об осесим-метричности колебаний не дает существенного эффекта в плане ограничения амплитуд продольных автоколебаний, и во-вторых, учет эффекта параметрического возбуждения неосесимметричных тонов оказывает существенное влияние на характер переходного процесса и приводит к значительному ограничению амплитуд.
Следовательно, благодаря указанному эффекту нелинейность корпуса может быть определяющим фактором в ограничении амплитуд продольных автоколебаний жидкостной ракеты наряду с нелинейными кавитационными эффектами на входе в насосы и нелинейным демпфированием.
Основные результаты представленных в данной диссертации исследований опубликованы в работах [14, 15, 34 — 48, 106, 107, 182].
Диссертация состоит из введения, пяти глав, разбитых на 26 разделов, заключения и списка литературы.
Основные результаты, полученные в процессе выполненных исследований, можно сформулировать следующим образом.
1. Сформулированы и доказаны основные теоремы метода корректирующих рядов, составляющие принципиально новую идеологическую основу модального синтеза подконструкций при исследовании динамических свойств сложных систем в ограниченном частотном интервале. Операция усечения ряда из собственных форм подконструкции заменяется усечением степенного корректирующего ряда при конечном числе собственных форм в модальном разложении колебаний. Коэффициенты степенного ряда вычисляются рекур-рентно с помощью последовательности статических задач.
2. Получена асимптотическая оценка погрешности усечения модального разложения при увеличении порядка корректирующего ряда, дающая априорную оценку точности результатов синтеза подконструкций.
3. Выведены соотношения метода корректирующих рядов как для дискретных моделей подконструкций, так и для континуальных моделей. Рассмотрены различные варианты синтеза в зависимости от условий закрепления внешних степеней свободы подконструкций при определении их собственных частот и форм: методы жестких и свободных границ, а также гибридный метод, когда часть внешних степеней свободы закреплена, а часть свободна. Исследованы различные методы формирования матриц динамических коэффициентов влияния подконструкций.
4. Разработан численно устойчивый алгоритм вычисления корректирующих векторов (или функций). Показано, что в ходе рекуррентного процесса они должны вычисляться в подпространстве, ортогональном к учтенным в разложении собственным формам, при этом на каждом шаге должна выполняться дополнительная ортогонализация решения.
5. Разработан программный комплекс для исследования динамических характеристик сложных пространственных стержневых систем, включающий программу расчета собственных частот и форм колебаний подконструкций, препроцессор, постпроцессор для формирования баз данных, содержащих информацию о динамических характеристиках подконструкций и применяемых при синтезе корректирующих векторов, а также программы для исследования динамических характеристик составной конструкции с использованием метода корректирующих рядов. На основе аналогичных принципов разработан программный комплекс для исследования составных осесимметричных оболочеч-ных конструкций, содержащих жидкость. Проведены численные исследования сходимости метода корректирующих рядов на ряде простых примеров, а также на примерах конструкций типа корпусов жидкостных ракет и орбитальной космической станции (стержневая модель). Полученные результаты демонстрируют высокие показатели точности метода корректирующих рядов и скорости вычислений при невысокой требовательности к параметрам вычислительных систем.
6. Сформулирована непротиворечивая постановка краевой задачи и вариационный принцип для описания динамического поведения упругой конструкции, взаимодействующей с ограниченным объемом жидкости в условиях однородного гравитационного поля. Предложенная формулировка применима к произвольным регулярным (кусочно-гладким) поверхностям контакта конструкции с жидкостью.
7. Разработан алгоритм решения задачи определения динамических характеристик осесимметричных оболочечных конструкций, взаимодействующих с ограниченными объемами жидкости, основанный на методе конечных элементов. Вычисляются как осесимметричные, так и неосесимметричные формы колебаний с учетом влияния начального напряженно-деформированного состояния, обусловленного внутренним давлением в полостях конструкции (включая гидростатическое давление) и ее собственным весом. Программная реализация алгоритма на различных электронных вычислительных системах показала его высокую надежность в процессе многолетней эксплуатации в условиях конструкторских бюро ракетно-космической отрасли. Он включен в фонд алгоритмов РКА и использовался в процессе про-ектно-конструкторских работ при создании ракет Зенит, Энергия-Буран, Космос, Рокот, Прибой и других. Новый программный комплекс позволяет в полном объеме учитывать гравитационные эффекты, связанные с образованием поверхностных волн и деформацией поверхности контакта жидкости со стенками сосуда.
8. Разработана методика исследования амплитуд продольных автоколебаний жидкостной ракеты с учетом нелинейности деформаций оболочек корпуса и эффекта параметрического возбуждения неосесимметричных форм колебаний. Проведенное исследование показало, что, во-первых, учет геометрической нелинейности в рамках предположения об осесимметричности колебаний не дает существенного эффекта в плане ограничения амплитуд продольных автоколебаний, и во-вторых, учет эффекта параметрического возбуждения неосесимметричных тонов оказывает существенное влияние на характер переходного процесса и приводит к значительному ограничению амплитуд продольных колебаний. Этот результат показывает необходимость учета нелинейности поведения корпуса при оценке максимальных амплитуд продольных колебаний, развивающихся в связи с неустойчивостью в контуре «корпус — топливная магистраль — двигатель», а методика открывает путь исследования данного вопроса в ходе проектно-конструкторских разработок данного класса изделий.
— 305 -ЗАКЛЮЧЕНИЕ.
В диссертации представлена методология исследования динамических свойств сложных упругих и гидроупругих систем, основанная на корректном и непротиворечивом подходе к задаче определения динамических характеристик входящих в систему конструкций, содержащих жидкость, и высокоэффективном и надежном методе модального синтеза подконструкций, обеспечивающем оценку точности получаемых результатов.
