Математическая модель динамического поведения конструкций на основе низкомодульных материалов Текст научной статьи по специальности «Физика»

УДК 624.04

В.В. Павлоградский, Р.В. Бульбович

Пермский государственный технический университет

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ДИНАМИЧЕСКОГО ПОВЕДЕНИЯ КОНСТРУКЦИЙ НА ОСНОВЕ НИЗКОМОДУЛЬНЫХ МАТЕРИАЛОВ

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

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

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

Для решения задачи используется следующий вариационный принцип: среди всех допустимых значений перемещений фактически существующие определяются вариационным уравнением [1], которое записывается в следующем виде:

ЕЖ1 - и)

V У \^У У V ^ У

+2ц

дих диу диу ди„ дих ди„

+ ■

- + ■

ч дх ду ду дг дх дг у

Г ди диу У (диу ди л2

+

+

[ду дх ) [ & ду)

+

( дих ди ^ 2

+ х + —-

N д дх )

+

где п - количество конечных элементов; V - объем конечного элемента;

5 - площадь конечного элемента;

1 - 2ц

X - постоянная Ляме, X =

2

и2 - перемещения вдоль осей координат;

х, у, z - координаты;

Е - модуль упругости материала;

ц - коэффициент Пуассона материала;

Рх, Ру, Pz - поверхностные силы вдоль осей координат.

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

Уравнение равновесия системы конечных элементов, находящейся в состоянии движения, имеет вид

[ М ] {5} + [ К ]{5} = ДО)}, (2)

где [ К ] - матрица жесткости конструкции;

[М ] - матрица масс конструкции;

{5} - вектор неизвестных (перемещений и их производных);

{5} - вектор узловых ускорений;

)} - вектор нагрузки в момент времени t.

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

При прямом интегрировании уравнения (2) используется численная пошаговая процедура. Термин «прямое» означает, что перед интегрированием не производится никаких преобразований уравнений. Прямое интегрирование основано на двух идеях. Во-первых, удовлетворение условий равновесия (2) требуется не в любой момент времени t, а только на отдельных коротких отрезках времени At. Это означает, что равновесие с учетом сил инерции рассматривается в дискретных точках временного интервала. Следовательно, становится возможным эффективное использование в методах прямого интегрирования всего вычислительного аппарата статистического анализа. Во-вторых, учитывается изменение перемещений, скоростей и ускорений внутри каждого временного интервала At. Именно способ учета этих изменений определяет точность, устойчивость и экономичность процедуры решения.

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

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

Ф(. - вектор г-й собственной формы; юг - г-я собственная частота колебаний.

Соотношение (3) - уравнение равновесия системы с одной степенью свободы с единичной массой и жесткостью ю2.

Решение каждого уравнения (3) можно получить с помощью интеграла Дюамеля [3]:

Интеграл Дюамеля должен в общем случае вычисляться численно. Для получения полной реакции системы необходимо найти решения всех п уравнений (3). Перемещения узловых точек получаются суперпозицией реакций системы по всем формам:

Таким образом, для получения реакции системы методом разложения по собственным формам требуется, во-первых, вычислить собственные значения и собственные векторы системы (2), затем решить уравнения равновесия (3) и, наконец, сложить реакции по каждой собственной форме в соответствии с (4).

Из методов прямого интегрирования наибольшее распространение нашли методы центральных разностей, метод Хаболта, 0-метод Вилсона и метод Ньюмарка [2]. В данной статье используется метод Ньюмарка как метод, обладающий лучшими характеристиками точности и устойчивости.

хі (ї) + ю2 хі (ї) = г (ї),

(3)

где г (ї) = ф*Я(ї);

і= 1, 2, ..., п;

п

(4)

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

E(гю) = E*(cos фЕ + i sin фЕ ) (5)

где Е* - динамический модуль упругости;

фЕ - угол сдвига фаз между продольной деформацией и напряжением в образце при гармоническом нагружении частотой со.

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

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

полученными экспериментально методами идентификации [5]. Искомые параметры (перемещения узлов, деформации и напряжения) также становятся комплексными величинами.

Для анализа динамического поведения конструкций необходимо также знать резонансные частоты, которые можно определить из амплитудно-частотной характеристики (АЧХ) в исследуемом диапазоне частот при проведении гармонического анализа. Однако построение АЧХ - это трудоемкий процесс, поэтому для определения резонансных частот необходимо решать собственную задачу. Задачи о собственных колебаниях имеют важное значение для анализа поведения конструкций при динамических нагрузках. Процедура МКЭ сводит задачу об определении собственных значений к алгебраической проблеме нахождения ненулевых решений следующей системы линейных алгебраических уравнений:

([ К]-X [ М ]){«)= 0, (6)

где X = ю2;

ю - собственная частота конструкции.

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

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

ю — Ю1 + т2,

где ю1 - действительная часть комплексной собственной частоты;

ю2 - мнимая часть комплексной собственной частоты.

Коэффициент демпфирования колебаний, являющийся наряду с собственной частотой важнейшей динамической характеристикой конструкции, определяется из соотношения [3]:

ю2

у — —. ю1

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

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

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

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

к

"її

, 1 ,

Рис. 1. Двухслойный стержень

На рис. 2 представлена зависимость первой собственной частоты продольных колебаний двухслойного стержня (/ = ю / 2л) от модуля упругости низкомодульного материала. Из рисунка видно, что при изменении модуля упругости от 10 до 100 МПа собственная частота практически не изменяется и в основном определяется продольной собственной частотой металлической оболочки, равной 312 Гц. При увеличении модуля упругости от 100 до 1000 МПа первая частота продольных колебаний двухслойного стержня увеличивается от 318 до 362 Гц.

На рис. 3 представлена зависимость коэффициента демпфирования колебаний двухслойного стержня от модуля упругости низкомодульного материала. Угол сдвига фаз между напряжением и деформацией низкомодульного материала принимался равным 0,2. Модуль упругости низкомодульного материала изменялся от 10 до 1000 МПа. Из рисунка видно, что при небольших значениях модуля упругости демпфирование в конструкции будет практически отсутствовать. При Е = 1000 МПа коэффициент демпфирования составляет 0,026.

У

Рис. 2. Зависимость собственной частоты колебаний от модуля упругости

Рис. 3. Зависимость коэффициента демпфирования колебаний от модуля упругости

На рис. 4 представлена зависимость коэффициента демпфирования колебаний двухслойного стержня от угла сдвига фаз между напряжением и деформацией для двух значений динамического модуля упругости низкомодульного материала 100 и 1000 МПа. Из рисунка видно, что эта зависимость в обоих случаях является линейной.

На рис. 5 показано изменение коэффициента усиления колебаний Р в области первой собственной частоты вязкоупругой конструкции для различных значений коэффициента демпфирования. Из этих кривых видно, что в области, когда частота возмущающей силы мала по сравнению с первой собственной частотой, коэффициент усиления незначительно отличается от единицы, а перемещения являются такими же, как в случае статического действия силы P sin ю?. После прохождения первого резонанса величина коэффициента усиления приближа-

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

0,1

ОЛ 0.1 фЕ

Рис. 4. Зависимость коэффициента демпфирования колебаний от угла сдвига фЕ : 1 - Е = 100 МПа;

2 - Е* = 100 МПа

Рис. 5. Изменение коэффициента усиления колебаний в области первой собственной частоты при различных значениях коэффициента демпфирования конструкции: 1 - у = 0,02; 2 - у = 0,05; 3 - у = 0,10

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

Работа выполнена при финансовой поддержке ФЦП «Исследования и разработки по приоритетным направлениям развития научнотехнического комплекса России на 2007-2012 годы» (госконтракт № 02.518.11.7135).

Библиографический список

1. Прочность, устойчивость, колебания: справочник. Т. 1 / под ред. И.А. Биргера и Я.Г. Пановко. - М.:Машиностроение, 1968. - 831 с.

2. Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов. - М.: Стройиздат, 1982. - 448 с.

3. Тимошенко С.П., Янг Д.Х., Уивер У. Колебания в инженерном деле. - М.: Машиностроение, 1985. - 472 с.

4. Бленд Д. Теория линейной вязкоупругости. - М.: Мир, 1965. -

230 с.

5. Бульбович Р.В., Пальчиковский В.Г., Павлоградский В.В. Метод определения динамических деформационных свойств мягких вязкоупругих материалов //Наука-производству. - М.: Вираж-Центр, 1999. - № 12 (25). - С. 14-18.

6. Корн Г., Корн Т. Справочник по математике. - М.: Наука, 1973. - 832 с.

Получено 6.12.2010