CC BY
28
Наука и Образование
МГТУ им. Н.Э. Баумана
Наука и Образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2016. № 05. С. 152-174.
]Э5М 1994-040В
DOI: 10.7463/0516.0839190
Представлена в редакцию: Исправлена:
© МГТУ им. Н.Э. Баумана
УДК 519.711.2
Использование метода конечных элементов совместно с методом возмущений в задаче вычисления расщепления частоты оболочки с дефектом формы срединной поверхности
Вахлярский Д. С.1*, Гуськов А. М.1, Басараб М. А.1, Матвеев В. А.1
10.04.2016 24.04.2016
УакЫу аг&ку gmail.com
1МГТУ им. Н.Э. Баумана, Москва, Россия
В статье рассмотрена задача нахождения расщепления собственных частот оболочки, близ -кой к оболочке вращения, при наличии малых дефектов формы срединной поверхности. Задача решена методом конечных элементов в комбинации с методом возмущений. Приводится описание конечного элемента оболочки вращения с произвольной формой меридиана. Проведена верификация используемого подхода на цилиндрической и полусферической оболочках. Проведено сравнение результатов использованного подхода с результатами, полученными в коммерческом конечно-элементном комплексе Ansys.
Ключевые слова: волновой твердотельный гироскоп, расщепление частот, дефект формы срединной поверхности
Введение
Волновой твердотельный гироскоп (ВТГ) - высокоточный высоконадежный прибор измерения инерциальной информации. ВТГ обладает рядом качеств, наличие которых делает перспективным использование такого типа датчиков в микроэлектромеханических системах. Работая в режиме интегрирующего гироскопа, ВТГ является одним из немногих датчиков, выходной сигнал которого пропорционален непосредственно углу поворота прибора, а не угловой скорости. В настоящее время ведутся интенсивные исследования в области технологии получения чувствительных элементов микроэлектромеханических ВТГ [1,2].
Чувствительный элемент волнового твердотельного гироскопа - резонатор, выполненный в виде тонкостенной оболочки вращения из кварцевого стекла, сапфира, рубина, бронзы [3]. Изготовление резонатора невозможно без погрешностей. Это приводит к расщеплению Д/2 рабочей частоты / свободных колебаний чувствительного элемента (индекс 2 отвечает форме колебаний с 4 узловыми меридианами и без узловых параллелей).
Указанное расщепление, существенным образом влияет на точность показаний прибо-
Влияние геометрических неоднородностей резонатора ВТГ на его динамику ранее рассматривалось в ряде работ. Так, в [4-6] исследуются типы геометрических неоднородностей кольцевого резонатора, наиболее влияющие на расщепление частот и точностные характеристики прибора. При этом как формы колебаний, так и формы неоднородности раскладываются в ряды Фурье с последующим количественным и качественным анализом влияния на динамику искажения формы. Еще в [7,8] был подробно исследован случай второй формы неоднородности (эллиптичность) кольцевого резонатора. Выражения, позволяющие оценить динамику цилиндрической оболочки такого типа (эллиптический цилиндр) приведены в [9]. Вывод уравнений динамики оболочки вращения с ненулевой гауссовой кривизной и произвольным очертанием линии меридиана рассмотрен в [10]. При исследовании колебаний оболочек идеальных в окружном направлении и имеющих неоднородность в меридиональном направлении могут быть использованы численные алгоритмы, в частности, в [11] предложен подход к расчету динамики (собственные формы колебаний и коэффициент прецессии) оболочки вращения с произвольной формой образующей.
Для вычисления / необходимо решить задачу на собственные значения для свободных колебаний оболочки. В операторной форме задача на собственные значения имеет вид [12]:
где Т,0 - линейные операторы; и - собственный вектор;
Ак = (2я/^)2 - собственное значение.
На практике, отклонение геометрических параметров резонатора от идеальных, очень мало. При математическом моделировании, это обстоятельство позволяет применять метод возмущений для вычисления расщепления собственной частоты. Задача нахождения расщепления собственного значения резонатора методом возмущений имеет следующий вид [12]:
ра [3].
1. Метод возмущений
Тик - ЛкСик
(1)
(2)
□
е :
где , - линейные операторы, собственное значение и собственный вектор
невозмущенной задачи, соответствующей идеальной оболочке вращения;
Тп), в и^ - возмущения соответствующих величин порядка п (п > 0).
Окончательная постановка задачи (2) для оболочки с заданной геометрией меридиана и заданными погрешностями изготовления требует записи конкретного вида операторов Т(0), О(0) и возмущений этих операторов Т(п), О(п) до порядка п, где п - искомое возмущение собственного значения.
Дальнейшее математическое описание задачи принципиально зависит от типа дефектов изготовления резонатора. Возможны два варианта, каждому из которых соответствует своя группа дефектов:
• Дефекты, не вызывающие искажения срединной поверхности оболочки резонатора;
• Дефекты, искажающие срединную поверхность оболочки резонатора.
К первой группе дефектов относятся: неоднородность распределения толщины оболочки, плотности, упругих и диссипативных свойств материала в окружном направлении. Влияние данной группы дефектов на расщепление частоты резонатора хорошо изучено [3, 13-16]. Это связано с тем, что срединная поверхность оболочки с дефектами такого типа остается поверхностью вращения. В таком случае, для построения, как невозмущенных операторов Т(0),О(0), так и возмущений Т(п),О(п) в (2) используется известная система уравнений, описывающая колебания оболочек вращения. Эти уравнения имеют достаточно простой вид [17], возмущения операторов задачи (2) могут быть получены аналитически для достаточно широкого класса дефектов. В данной работе, влияние этой группы дефектов на расщепление Д/2 не рассматривается.
Ко второй группе дефектов относится собственно искажение срединной поверхности резонатора. В таких случаях, вопрос вычисления расщепления собственной частоты остается мало исследованным [18]. Это связано с тем, что срединная поверхность оболочки перестает быть поверхностью вращения. В этом случае, задача (1) описывается уравнениями свободных колебаний общей теории оболочек. И хотя в задаче (2) для построения операторов Т(0), О(0) и решения невозмущенной задачи используются уравнения динамики оболочек вращения, для построения возмущений операторов Т(п), О(п) (п>0), уже приходится рассматривать уравнения общей теории оболочек. Эти уравнения имеют достаточно сложный вид. Ситуация осложняется тем, что оболочки с дефектами указанного типа чаще всего невозможно параметризовать ортогональными гауссовыми координатами. Что приводит к еще большему нагромождению аналитических выражений. Ко всему прочему, метод возмущений сам по себе приводит к сильному «разрастанию» аналитических выражений.
В связи с этим, аналитическое вычисление возмущения собственного значения хотя бы первого порядка для конкретной оболочки с простой формой меридиана представляется авторам практически вряд ли возможным, даже в случае простейших возмущений гео-
метрии. Возможно, единственным исключением является цилиндрическая оболочка с граничными условиями Навье [19].
Применение метода конечных элементов (МКЭ) существенно облегчает задачу нахождения расщепления рассматриваемой частоты и устраняет необходимость построения возмущений операторов Т и С. При использовании МКЭ задается реальная геометрия оболочки и вычисляются собственные частоты. Расщепление просто равно разности полученных частот. Но в таком случае, возникает вопрос о точности вычислений, так как дефекты формы срединной поверхности малы.
Чтобы избежать построения возмущений операторов задачи и для оценки точности вычислений частот в коммерческих конечно-элементных комплексах, было разработано решение задачи (1) методом возмущений (2) с использованием специальных конечных элементов (КЭ).
Для задачи свободных колебаний, после дискретизации конечными элементами, от операторной формы (1) приходим к матричной форме задачи на собственные значения:
- А!^ = 0 (3)
где К, М - матрицы жесткости и масс КЭ модели оболочки соответственно;
Чк - собственный вектор;
^ - собственное значение.
После применения метода возмущений получаем следующие уравнения:
г1 : (К(0) -%6>М(0> )<£> = (л^М^ + -К(1) }<£>. (5)
+ + ^0)м<2) -кЯ ' (б)
где индексом 0 отмечены невозмущенные величины, а индексом п (п>0) отмечены возмущения порядка п соответствующих величин.
Из уравнения (4) получаем невозмущенное собственное значение А(0) и соответствующий ему собственный вектор .
Резонатор ВТГ работает на частоте /, соответствующей корню Я(0) задачи (4),
имеющему кратность т = 2. Кратному собственному значению соответствует семейство собственных векторов. На рисунке 1 показана одна из форм колебаний резонатора, соответствующая собственному вектору q(20J) . Вектор q, полученный поворотом волновой
картины, отвечающей q(0|), на произвольный угол, так же является собственным вектором, отвечающим данному собственному значению Я(0). Таким образом, в зависимости от выбранного q(0n , будут получены различные и , определяемые ориентацией де-
фектов резонатора по отношению к . Аналогично для А^) и q2(2) в зависимости от q(02). Согласно свойствам отношения Рэлея [12, 20], необходимо выбрать такие и , которые в первом приближении дают максимальное ( А ^ * j и минимальное
( ^ -'max
( А^1"*) значения l (l. Значения ( А^) и (А^1 отвечают расщепленным собст-
^ -'min ~ ^ -'max ^ 'min
венным частотам (f )max и (f ) оболочки с дефектами.
Рисунок 1 - формы колебаний, отвечающие кратным собственным частотам резонатора без погрешностей
Для определения ( А ^ и ( А ^ необходимо «приспособить» невозмущенные
V / тпу V / rniri
собственные векторы q(0) и q(0) , отвечающие двукратному собственному значению
(0)
к возмущениям К(1) и М(1) [21]. Для этого, представим собственный вектор ч(0) в более
r(i)
(0)
общем виде:
q(0) = C q(0) + C q(0) q2 Clq2(l) C2q2(2).
(7)
После подстановки (7) в (5) для к=2, получаем:
(Кт -^0)М(9)= (К®-Л™М® -^>М<0) )(С1Ч™ + С^). (8)
Неоднородная система уравнений (8) имеет ненулевое решение только в случае, если
,(<>)„ я(0) 12(1) и Чк(2)
её правая часть ортогональна к решениям и q((! однородной системы (4), отвечаю-
щим собственному числу , то есть когда:
(к® - л^м®+ ) = о [ (к® - jjpfcf*-йГм^ )(ciq® + )=о
А( Л® jC = О.
(9)
где С = {С С2 }Т - вектор неизвестных.
Получаем систему уравнений относительно С и С. Условие разрешимости системы (9) - вековое уравнение, представляющее собой равенство нулю её определителя:
с1 е ^ А (41 } ) ) = 0 . (10)
/1
Из уравнения (10) находим Я 2 с i ) и Я 2/2 у Согласно свойствам отношения Рэлея
L2(l)
2(2)
я:
2 с 1 ) V 2 ;
я
(1) = 2 С 2 )
( Я 21 )) , так как (10) является задачей на собственные значения.
____________^ -'min
Из системы (9), используя полученные ( Я ^ ) ) и ( Я ^ )) определяем отвечаю-
^ -'max ^ -'min
щие им C(1) и C(2) нормированные на единицу. Используя (7), находим «приспособленные» невозмущенные собственные векторы q(0(]^ax) и q((im), отвечающие максимальному и минимальному возмущениям собственного значения Я /0 ) кратности m=2.
l2(max) и q 2(min)
решить систему уравнений (5). Поскольку ЯС,0) является собственным значением системы (4), система уравнений (5) является вырожденной. Для её однозначного решения требуются дополнительные уравнения. Рассмотрим разложение собственного вектора по малому параметру в:
Для нахождения возмущения q(1) собственных векторов q^^ и необходимо
Представим все возмущения разложенными на составляющие: сонаправленную не-
возмущенному вектору q(0) и ортогональную ему: После подстановки (12) в (11), имеем:
% - (1+/V+/V+. ОчГ+¿чГ+
(12)
(13)
Из выражения (13) видно, что составляющие, сонаправленные невозмущенному собственному вектору, не приводят к его возмущениям, так как собственный вектор опреде-
лен с точностью до произвольного множителя. Таким образом, все возмущения qj¡n) (п>0)
должны быть ортогональны невозмущенному собственному вектору q['. Поскольку исследуется двукратное собственное значение, имеем две следующие системы для нахождения q(21(p):
В качестве Чад в (14) необходимо подставлять приспособленные к возмущениям
л(0) и л(0) Ч2(шах) и Ч2(шт) •
Для нахождения возмущения Я^2^ второго порядка собственного значения, умножим (6) на Ч(0))Т слева, учитывая (4) и симметрию матриц К(0) и М(0). Получим:
Возмущение ч(22\ вычисляется решением системы (6), при условии его ортогональ-
ности собственным векторам Ч(0)
2(Р)
Получение возмущения более высокого порядка аналогично процедуре полу-
1(2) чения •
Основная часть возмущений собственного значения содержится в первом приближе-1 2 (р ).
С1")
нии Я^А. Однако, как было показано в многочисленных работах, например в [22], что ес-
(1)
ли окружное волновое число равно к, то расщепление собственного значения ДЯд. , отвечающего числу к, вызывается только гармониками дефектов с номером 2к. Таким образом, в случае возникновения необходимости исследования влияния дефектов, распределенных в окружном направлении по гармоникам отличным от 2к, потребуется вычисление
возмущения собственного значения как минимум второго порядка •
2. Описание конечного элемента
Рассматриваемый здесь конечный элемент, имеет два узла, представляет собой тонкий кольцевой элемент оболочки и показан на рисунке 2.
При построении конечного элемента используется реальная геометрия рассматриваемой оболочки:
где <р - гауссовы координаты оболочки;
( 1 ) > (2 ) - координаты сечений оболочки, выделяющий конечный элемент.
Рисунок 2 - Геометрия оболочки и конечного элемента
Перемещения в пределах конечного элемента, в соответствии с методом возмущений, разложены в ряд по малому параметру в:
где - вектор перемещений в пределах элемента;
и (0 ) - вектор перемещений на элементе, соответствующий невозмущенной (симметричной) задаче;
- возмущения вектора перемещений на элементе. Поскольку оболочка замкнута, каждый из векторов и можно разложить в ряд по окружному углу ф.
N.
¡(и)
м,
(19)
(=1
где I - номер гармоники в разложении перемещений по окружному углу ф; N - число удерживаемых в разложении гармоник.
Каждый из векторов и^ (1 = А, С, Б ) в разложении (19) представляется проекциями на оси цилиндрической системы координат:
где - орты глобальной цилиндрической системы координат ( ;
М1("),^"г,и{ь - проекции на орты цилиндрической системы координат(г,<р ,х) .
1/ 1 ' 1/ 2 '"1/ 3
Для аппроксимации проекций м2") используются кубические полиномы:
1/J
где - меридиональная координата в пределах конечного элемента.
В качестве неизвестных в узлах принимаются значения и(П"11, и2^2 и их производные по меридиональной координате
номер в верхнем индексе обозначает номер узла. Таким образом каждая из проекций и представляется следующим образом:
(п) Ч J
где и^ - значение и(П) в первом узле;
и^2 и,
1J - значение 11J во втором узле;
д (п) 1 „ и2П) „ г
ои\ ■ - значение производной 11J по меридиональной координате <; в первом узле;
»
д (п) 2 . и;, „ г
о и> у - значение производной 11■> по меридиональной координате $ во втором уз-
ле.
Учитывая вышеописанные разложения и аппроксимации, имеем следующий вектор узловых перемещений элемента для приближения порядка п:
Перемещения в пределах конечного элемента выражаются через узловые перемещения следующим образом:
где и - арифметический вектор перемещений КЭ; и ( <р) - матрица преобразования;
q(n) - вектор узловых перемещений элемента приближения порядка п. Матрица преобразования и ( <р ) одинакова для всех q(n), так как для каждого и (п) используется одна и та же аппроксимация.
Для получения матрицы жесткости и матрицы масс элемента используются нижеследующие геометрические параметры.
Для вычисления координат первого метрического тензора вычислим первые и вторые производные радиус-вектора оболочки:
дг
д2г
Г
г = ■
где .
Координаты первого метрического тензора:
, (г = 1,2, у = 1.2),
(25)
= Г/ " Г;
(26)
Определитель первого метрического тензора:
а а11а(2 - а1(.
(27)
Орт нормали равен:
й— Г1ХГ2
(28)
Символы Кристоффеля - Шварца:
где Г'к - векторы взаимные к г^(/=1,2, А=1,2).
С учетом (17), разложение координат первого метрического тензора имеет вид:
(29)
Как уже упоминалось, с учетом разложения (17) аналитическая запись выражений (28) - (29) является довольно затруднительной. Аналитическое выражение для матрицы масс и матрицы жесткости, с учетом возмущений геометрии и перемещений, вообще представляется, вряд ли возможным. Поэтому, в данной работе, при численной реализации метода, предлагается использовать для хранения всех возмущенных величин трехмерные массивы, и переопределить все используемые при построении конечного элемента математические операции.
Пояснение данного подхода и удобство его использования продемонстрируем на построении матрицы масс элемента.
Матрица масс конечного элемента оболочки выражается следующим образом [23]:
где I, ] - номера точки интегрирования по схеме Гаусса (рисунок 2); ш¿у ^ - весовой коэффициент интегрирования по схеме Гаусса;
- плотность материала оболочки; - координаты точки интегрирования;
- толщина оболочки;
- определитель первого метрического тензора оболочки.
В выражении (31), при возмущениях геометрии срединной поверхности (17), возмущенным будет только определитель первого метрического тензора . При использовании трехмерного массива для хранения возмущенных параметров и переопределения математических операций вычисление квадратного корня из определителя выполняется следующим образом:
1) Используя выражения (26), находим координаты первого метрического тензора оболочки, вводя операцию скалярного умножения векторов, записанных в трехмерные массивы следующим образом:
2) Далее используя такой же как в (32) метод умножения для скалярных величин находим определитель первого метрического тензора . по формуле (27).
3) Операция извлечения корня / а ( <р ) из определителя, записанного в трехмерный массив, определяется разложением в ряд Тейлора по параметру в:
Для трехмерного массива получаем:
где %(:,:,к)= 1); а(:,:,к)= а(к"1); к = 1,2,3,....
Таким образом, используя трехмерные массивы (для хранения возмущенных величин) и переопределение необходимых математических операций, выражения для матриц жесткости и масс имеют такой же вид, как если бы метод возмущений не использовался. А возмущения этих матриц вычисляются автоматически и записываются по третьему измерению используемых массивов.
Преимущество такого подхода заключается в том, что выражения для матриц жесткости и масс имеют такой же вид, как если бы метод возмущений не использовался. А возмущения этих матриц автоматически записываются по третьему индексу.
Алгоритм вычисления расщепления частоты выглядит следующим образом:
1) Используя разработанный конечный элемент, строятся матрицы К(0),...,К(п), м(0),..., м(п) всей оболочки резонатора вплоть до необходимого порядка п.
2) Из задачи (4) находим невозмущенные собственные числа Я(0 ) и соответствующие им собственные векторы q(0)) и q(20)2).
3) Решая уравнение (10), вычисляются возмущения собственных чисел первого порядка , и отвечающие им векторы констант С и С .
^ -'тах ^ 'тт 2 ) 2 )
4) Используя найденные С и С , по выражению (7) находим приспособленные к
возмущениям невозмущенные собственные вектора q(20^) и q20(im).
<2(1) и q 2(2)
5) Из системы (14) находим возмущения собственных векторов q2;^ и q^ первого
порядка.
6) По выражению (15) вычисляются возмущения собственных чисел второго порядка.
7) Решая две системы (16), находим возмущения собственных векторов второго порядка.
8) Далее, для вычисления приближений более высокого порядка повторяются шаги 6) и 7) с использованием соответствующих формул.
3. Численный эксперимент
Ниже представлены результаты верификации предлагаемого в работе подхода и проведено сравнение с результатами, полученными в коммерческом КЭ комплексе ANSYS. В качестве тестовых примеров, рассмотрены задачи нахождения расщепления
(1)
первого порядка для некруговой цилиндрической оболочки и полусферической обо-
лочки с простейшими возмущениями геометрии. Для цилиндрической оболочки проведено сравнение с теоретическим результатом, полученным в [19].
Г 1)
Параметры конечно-элементных моделей, принятые при вычислении Д/2 приведены в таблице 1.
Таблица 1 - Параметры конечно-элементных моделей.
Параметр Значение
Число элементов вдоль меридиана N е, 10
Число удерживаемых гармоник NN 4
Число точек интегрирования в меридиональном направлении ^ , 4
Число точек интегрирования в окружном направлении N р ^ 20
Число элементов вдоль меридиана в ЛЫ8У8 Nem 100
Число элементов вдоль параллели в ЛЫ8У8 ^ ^ 100
Параметры некруговой цилиндрической оболочки с граничными условиями Навье приведены в таблице 2.
Таблица 2 - Параметры некруговой цилиндрической оболочки.
Параметр цилиндрической оболочки Значение
модуль упругости Е, МПа 73600
плотность р, кг/м3 2210
коэффициент Пуассона 0,17
срединный радиус Я, мм 40
длина Ь, мм 80
толщина к, мм 1
Малый параметр возмущения геометрии е, мм 10-5
Частота Гц 7952,6
Было принято следующее несовершенство срединной поверхности цилиндрической оболочки:
1<р*Щ2и]. (35)
Параметры полусферической оболочки приведены в таблице 3.
Таблица 3 - Параметры полусферической оболочки.
Параметры полусферической оболочки Значение
модуль упругости Е, МПа 75000
плотность р, кг/м3 2600
коэффициент Пуассона 0,1
срединный радиус Я, мм 30
радиус ножки Я, мм 2
толщина к, мм 0,5
Малый параметр возмущения геометрии е, мм 10"6
Частота Гц 654,0
Для полусферической оболочки принято следующее несовершенство срединной поверхности:
= ь^ФЛфФЛк]. (Зб)
На рисунке 3 показана геометрия рассматриваемых оболочек, масштаб несовершенств намеренно увеличен для наглядности.
К *
Рисунок 3 - цилиндрическая и полусферическая оболочки с дефектом срединной поверхности, распределенным в окружном направлении по четвертой гармонике окружного угла.
Как видно из таблицы 4, первое приближение расщепления частоты вычисленное методом возмущений и с помощью АКБУБ практически идентичны, как для цилиндрической оболочки (таблица 2) с возмущениями (35), так и для полусферической оболочки (таблица 3) с возмущениями (36).
Цилиндрическая оболочка Полусферическая оболочка
Д/2( С), Гц Д/2( С)1, Гц С)1, Гц ^2 ^ Гц ^2 ^ Гц
(метод возмущений) (ЛШУ8) (теория [19]) (метод возмущений) (ЛN8У8)
23,854 23,859 23,856 0,16683 0,16672
Г (0) Гп 1 2 су1, Гц f (0) Гц 1 2 су1, f ( 0 ) Гц f 2 су1, Гц /•(0) Гц 1 2 БрЬ, Гц ( 0 ) Гц 1 2 БрЬ, Гц
(метод возмущений) (ЛШУ8) (теория [19]) (метод возмущений) (ЛN8У8)
7952,6 7953,5 - 654,0 652,7
д,(1) ,А0) су1/ ¡2 су1 су1/ ¡2 су1 ,А1) /АО) су1/ ¡2 су1 ,А1) /АО) ,А1) /АО)
(метод возмущений) (ЛШУ8) (теория [19]) (метод возмущений) (ЛШУ8)
0,003 0,003 0,003 0,000255 0,000255
На рисунках 4 и 5 показаны результаты влияния параметров конечно-элементной
(1)
модели на значение расщепления Д^ при использовании метода возмущений. При вычислениях варьируется один параметр, а остальные принимаются согласно таблице 1. Из рисунка 4 видно, что расщепление частоты цилиндрической оболочки в окрестности точки с выбранными по таблице 1 параметрами практически постоянно. Из рисунка 5 видно, что для полусферической оболочки расщепление частоты в окрестности точки с выбранными по таблице 1 параметрами изменяется не более чем на 2%. Таким образом, при использовании предлагаемого в работе конечного элемента совместно с методом возмущений, для рассмотренных тестовых задач сходимость решения обеспечивается.
Рисунок 4 - зависимость расщепления частоты Д/2 с у1 цилиндрической оболочки, вычисленного методом возмущений, в зависимости от: а) количества КЭ вдоль меридиана N е,; б) количества удерживаемых гармоник N1^ в) количества точек интегрирования вдоль меридиана N р г) количества точек
интегрирования вдоль параллели .
Рисунок 5 - зависимость расщепления частоты Д/2 пр ^ полусферической оболочки, вычисленного методом возмущений, в зависимости от: а) количества КЭ вдоль меридиана N е,; б) количества удерживаемых гармоник N1^ в) количества точек интегрирования вдоль меридиана N р г) количества точек
интегрирования вдоль параллели .
Рисунок 6 - зависимость расщепления частоты Д/2 цилиндрической оболочки, вычисленного в А№У8, в зависимости от: а) количества КЭ вдоль меридиана N е,; б) количества КЭ вдоль параллели N е р.
На рисунках 6 и 7 показаны зависимости расщепления собственной частоты цилиндрической и полусферической оболочек от количества элементов в меридиональном и окружном направлениях, вычисленные в программе АКБУБ. При расчете использовался конечный элемент оболочки вЬе11281 с восемью узлами из библиотеки АКБУБ. Из рисунков 6 и 7 видно, что увеличение количества элементов больше 100 как в окружном, так и в меридиональном направлении не приводит к сколько-нибудь значимому уточнению расщепления.
Рисунок 7 - зависимость расщепления частоты Д/2 полусферической оболочки, вычисленного в АКБУБ, в зависимости от: а) количества КЭ вдоль меридиана Ые,с; б) количества КЭ вдоль параллели Ые^.
Заключение
1) Разработан подход совместного использования метода конечных элементов и метода возмущений для вычисления расщепления собственной частоты оболочек вращения с малым дефектом срединной поверхности и произвольной формой меридиана;
2) Результаты решения тестовых задач предложенным методом хорошо согласуются с теоретическим решением и решением в коммерческом программном комплексе АКБУБ;
3) В соответствии с приведенным выше пунктом, для решения задачи оптимизации формы меридиана оболочки с целью уменьшения расщепления частоты возможно использование как представленного в работе метода, так и программного комплекса АКБУБ.
4) Небольшие искажения формы резонатора могут быть связаны с эффектами неоднородности толщины оболочки и поверхностной плотности материала. Вопросы устранения влияния дефектов такого вида (балансировка) применительно к цилиндрическому и полусферическому резонаторам рассмотрены в [24-26]. Интерес в дальнейшем представляет учет совместного влияния как геометрических (профиль срединной поверхности,
толщина), так и физических (плотность, модуль Юнга [27]) неоднородностей на динамику ВТГ, а также исследование возможностей взаимной компенсации этого влияния.
Список литературы
1. Heidari A., Chan M-L., Yang H-A., Jaramillo G., Taheri-Tehrani P., Fonda P., Najar N., Yamazaki K., Lin L., Horsley D. A. Hemispherical wineglass resonators fabricated from the microcrystalline diamond // Journal of Micromechanics and Microengineering. 2013. vol. 23, no. 12, pp. 8. DOI: 10.1088/0960-1317/23/12/125016
2. Pai P., Chowdhury F.K., Mastrangelo C.H., Tabib-Azar M., MEMS-Based hemispherical resonator gyroscopes // Conference: Sensors, 2012 IEEE. DOI: 10.1109/ICSENS.2012.6411346.
3. Лунин Б.С., Матвеев В.А., Басараб М.А. Волновой твердотельный гироскоп. Теория и технология. М.: Радиотехника, 2014. 176 с.
4. Hwang R.S., Fox C.H.J., McWilliam S. The in-plane vibration of thin rings with in-plane profile variations. Part I: General background and theoretical formulation // Journal of Sound and Vibration. 1999. vol. 220, no. 3, pp. 497-516. DOI: 10.1006/jsvi.1998.1963
5. Fox C.H.J., Hwang R.S., and McWilliam S. The in-plane vibration of thin rings with in-plane profile variations. Part II: Application to nominally circular rings // Journal of Sound and Vibration. 1999. vol. 220, no. 3, pp. 517-539. DOI: 10.1006/jsvi.1998.1962
6. Yilmaz E. and Bindel D. Effects of imperfections on solid-wave gyroscope dynamics // in Proc. IEEE Sensors, Baltimore, MD, USA, Nov. 2013, pp. 1331-1334.
7. Sato K. Free flexural vibrations of an elliptical ring in its plane // J. Acoust. Soc. Am. 1975. vol. 57, no. 1, pp.113- 115. DOI: 10.1121/1.380420
8. Brigham G.A. In-plane free vibrations of tapered oval rings // The Journal of the Acoustical Society of America. 1973, vol. 54, no. 2, pp.451-460.
9. Новожилов В.В. Теория тонких оболочек. СПб.: Изд-во С.-Петерб. ун-та, 2010. 380 с.
10. Карачун В.В., Мельник В.Н. Уравнения динамики оболочки вращения с ненулевой гауссовой кривизной и произвольным очертанием линии меридиана // Вестник двигателе-строения, 2009. №3. с.29-36.
11. Басараб М.А., Кравченко В.Ф., Матвеев В.А., Пустовойт В.И. Атомарные функции в задаче определения функций Рэлея и коэффициента прецессии резонатора волнового твердотельного гироскопа // Доклады Академии Наук, 2001, т. 376, № 4, с. 474-479.
12. Коллатц Л. Задачи на собственные значения (с техническими приложениями): пер. с нем. М.: Наука, 1968. 504 с.
13. Меркурьев И.В., Подалков В.В. Динамика микромеханического и волнового твердотельного гироскопов. М.: Физматлит, 2009. 228 c.
14. Астахов С.В. Нелинейные эффекты в динамике волнового твердотельного и микромеханического гироскопов в условиях медленно меняющихся параметров: дис. ... канд. техн. наук. М., 2012. 157 с.
15. Донник А.С. Влияние геометрической неоднородности и упругой анизотропии материала на точностные характеристики волнового твердотельного гироскопа: дис. ... канд. техн. наук. М., 2006. 131 с.
16. Лунин Б.С. Физико-химические основы разработки полусферических резонаторов волновых твердотельных гироскопов. М.: Изд-во МАИ, 2005. 224 c.
17. Бидерман В.Л. Механика тонкостенных конструкций. Статика. М.: Машиностроение, 1977. 488 c.
18. Heidari A., Chan M., Yang H., Jaramillo G., Taheri-Tehrani P., Fonda P., Najar H., Yamazaki K., Lin L., Horsley D. Hemispherical wineglass resonators fabricated from the microcrystalline diamond // Journal of Micromechanics and Microengineering, 2013, Vol. 23, no. 12, pp. 125016-23(8). DOI: 10.1088/0960-1317/23/12/125016
19. Козубняк С.А. Расщепление собственных частот колебаний цилиндрического резонатора волнового твердотельного гироскопа, вызванное возмущением формы // Вестник МГТУ им. Н.Э. Баумана. Сер. Приборостроение. 2015. № 3. C. 39-49. DOI: 10.18698/0236-39332015-3-39-49
20. Ланкастер П. Теория матриц. М.: Наука, 1973. 280 с.
21. Маделунг Э. Математический аппарат физики. М.: ФИЗМАТЛИТ, 1949. 618 с.
22. Нарайкин О.С., Сорокин Ф.Д., Козубняк С.А Расщепление собственных частот кольцевого резонатора твердотельного волнового гироскопа, вызванное возмущением формы / Вестник МГТУ им. Н. Э. Баумана. Сер. Машиностроение. 2012. №6. С. 176-185.
23. Голованов А.И., Тюленева О.Н., Шигабутдинов А.Ф. Метод конечных элементов в статике и динамике тонкостенных конструкций М.: Физматлит, 2006. 392 с.
24. Матвеев В.А., Басараб М.А., Лунин Б.С. Аппроксимация распределения плотности резонатора волнового твердотельного гироскопа по измеренным параметрам дебаланса // Приборы и системы. Управление, контроль, диагностика. 2015. №10. С.9-16.
25. Басараб М.А., Лунин Б.С., Матвеев В.А., Чуманкин Е.А. Балансировка полусферических резонаторов волновых твердотельных гироскопов методом химического травления // Ги-роскопия и навигация. 2015. т.88, №1. С. 61-70.
26. Басараб М.А., Лунин Б.С., Матвеев В.А., Чуманкин Е.А. Статическая балансировка цилиндрических резонаторов волновых твердотельных гироскопов // Гироскопия и навигация. 2014. Т. 85, №2. С.43-51.
27. Chang Chia-Ou, Chang Guo-En, Chou Chan-Shin, et al. In-plane free vibration of a single-crystal silicon ring // Int. Journal of Solids and Structures. 2008, vol. 45, pp.6114-6132. D0I:10.1016/i .ijsolstr.2008.07.033
Science ¿Education
of the Baumail MSTU
Science and Education of the Bauman MSTU, 2016, no. 05, pp. 152-174.
DOI: 10.7463/0516.0839190
Received: 10.04.2016
Revised: 24.04.2016
© Bauman Moscow State Technical Unversity
Using a Combination of FEM and Perturbation Method in Frequency Split Calculation of a Nearly Axisymmetric Shell with Middle Surface Shape Defect
D.S. Vakhlyarskiy1'*, A.M. Guskov1, M.A. Basarab1, V.A. Matveev1
:Bauman Moscow State Technical University, Moscow, Russia
Vakhly arsky DS@ gmail.com
Keywords: hemispherical resonator gyroscope, HRG, frequency split, middle surface shape defect, surface imperfections
This paper proposes a method to calculate the splitting of natural frequency of the shell of hemispherical resonator gyro. (HRG). The paper considers splitting that arises from the small defect of the middle surface, which makes the resonator different from the rotary shell. The presented method is a combination of the perturbation method and the finite element method. The method allows us to find the frequency splitting caused by defects in shape, arbitrary distributed in the circumferential direction. This is achieved by calculating the perturbations of multiple natural frequencies of the second and higher orders. The proposed method allows us to calculate the splitting of multiple frequencies for the shell with the meridian of arbitrary shape.
A developed finite element is an annular element of the shell and has two nodes. Projections of movements are used on the axis of the global cylindrical system of coordinates, as the unknown. To approximate the movements are used polynomials of the second degree. Within the finite element the geometric characteristics are arranged in a series according to the small parameter of perturbations of the middle surface geometry.
Movements on the final element are arranged in series according to the small parameter, and in a series according to circumferential angle. With computer used to implement the method, three-dimensional arrays are used to store the perturbed quantities. This allows the use of regular expressions for the mass and stiffness matrices, when building the finite element, instead of analytic dependencies for each perturbation of these matrices of the required order with desirable mathematical operations redefined in accordance with the perturbation method.
As a test task, is calculated frequency splitting of non-circular cylindrical resonator with Navier boundary conditions. The discrepancy between the results and semi-analytic solution to this problem is less than 1%. For a cylindrical shell is made a comparison of results with solution in ANSYS commercial complex - a difference is less than 1%. For a hemispherical shell was found the frequency splitting. The comparison has shown that a discrepancy between the results
and ANSYS solution is less than 1%. The solution of this problem allows us to estimate further a mutual influence of defects of different nature (shape, thickness, density, modulus of elasticity, etc.) on splitting the frequency of the HRG. This is an urgent problem in terms of balancing the HRG resonators.
References
1. Heidari A., Chan M-L., Yang H-A., Jaramillo G., Taheri-Tehrani P., Fonda P., Najar N., Yamazaki K., Lin L., Horsley D. A. Hemispherical wineglass resonators fabricated from the microcrystalline diamond. Journal of Micromechanics and Microengineering, 2013, vol. 23, no. 12, 8 p. DOI: 10.1088/0960-1317/23/12/125016
2. Pai P., Chowdhury F.K., Mastrangelo C.H., Tabib-Azar M., MEMS-Based hemispherical resonator gyroscopes. Conference: Sensors, 2012 IEEE. DOI: 10.1109/ICSENS.2012.6411346
3. Lunin B.S., Matveyev V.A., Basarab M.A. Volnovoy tverdotel'nyy giroskop. Teoriya i tekhnologiya [Hemispherical resonator gyroscope. Theory and technology]. Moscow, Radiotekhnika Publ., 2014. 176 p. (in Russian).
4. Hwang R.S., Fox C.H.J., and McWilliam S. The in-plane vibration of thin rings with in-plane profile variations. Part I: General background and theoretical formulation. Journal of Sound and Vibration, 1999, vol. 220, no. 3, pp. 497-516. DOI: 10.1006/jsvi.1998.1963
5. Fox C.H.J., Hwang R.S., and McWilliam S. The in-plane vibration of thin rings with in-plane profile variations. Part II: Application to nominally circular rings. Journal of Sound and Vibration, 1999, vol. 220, no. 3, pp. 517-539. DOI: 10.1006/jsvi.1998.1962
6. Yilmaz E. and Bindel D. Effects of imperfections on solid-wave gyroscope dynamics. Proc. IEEE Sensors, Baltimore, MD, USA, Nov. 2013, pp. 1331-1334.
7. Sato K. Free flexural vibrations of an elliptical ring in its plane. J. Acoust. Soc. Am., 1975, vol. 57, no. 1, pp.113-115. DOI: 10.1121/1.380420
8. Brigham G.A. In-plane free vibrations of tapered oval rings. The Journal of the Acoustical Society of America. 1973. vol. 54. no. 2. pp. 451-460.
9. Novozhilov V.V. Teoriya tonkikh obolochek [Thin shells theory]. St Petersburg, SPnU Publ., 2010. 380 p. (in Russian).
10. Karachun V.V., Mel'nik V.N. Dynamic equations for rotational shell with zero Gaussian curvature and spontaneous shape of meridian line. Vestnik dvigatelestroyeniya, 2009, no. 3, pp.2936. (in Russian).
11. Basarab M.A., Kravchenko V.F., Matveyev V.A., Pustovoyt V.I. Atomic functions in the problem of determining Rayleigh function and precession coefficient for hemispherical resonator gyro. Doklady Akademii Nauk, 2001, vol. 376, no. 4. pp. 474-479. (in Russian).
12. Collatz L. Eigenwertaufgaben mit technischen anwendungen. Akademische verlagsgesellschaft Geest & Portig K.-G., 1963. 500 s.
13. Merkur'yev I.V., Podalkov V.V. Dinamika mikromekhanicheskogo i volnovogo tverdo-tel'nogo giroskopov [The dynamics of MEMS and hemispherical resonator gyroscopes]. Moscow, Fizmatlit Publ., 2009. 228 p. (in Russian).
14. Astakhov S.V. Nelineynyye effekty v dinamike volnovogo tverdotel'nogo i mikromekhanicheskogo giroskopov v usloviyakh medlenno menyayushchikhsya parametrov. Diss. kand. tekhn. nauk. [Nonlinear effects in the dynamics of hemispherical resonator and MEMS gyroscopes in case of slowly varying parameters: cand. tech. sci. diss.]. Moscow, 2012. 157 p. (in Russian).
15. Donnik A.S. Vliyaniye geometricheskoy neodnorodnosti i uprugoy anizotropii mate-riala na tochnostnyye kharakteristiki volnovogo tverdotel'nogo giroskopa: diss. kand. tekhn. nauk. [Influence of geometrical inhomogeneity and material elastic anisotropy on accuracy of a hemispherical resonator gyroscope: cand. tech. sci. diss.]. Moscow, 2006. 131 p. (in Russian).
16. Lunin B.S. Fiziko-khimicheskiye osnovy razrabotki polusfericheskikh rezonatorov volnovykh tverdotel'nykh giroskopov [Physical-and-chemical bases of designing hemispherical resonators of solid-state wave gyroscopes]. Moscow, MAI Publ., 2005. 224 p. (in Russian).
17. Biderman V.L. Mekhanika tonkostennykh konstruktsiy. Statika [Mechanics of thin-walled structures. Statics]. Moscow, Mashinostroyeniye Publ., 1977. 488 p. (in Russian).
18. Heidari A., Chan M., Yang H., Jaramillo G., Taheri-Tehrani P., Fonda P., Najar H., Yamazaki K., Lin L., Horsley D. Hemispherical wineglass resonators fabricated from the microcrystalline diamond. Journal of Micromechanics andMicroengineering, 2013, vol. 23, no. 12, pp. 125016-23(8). DOI: 10.1088/0960-1317/23/12/125016
19. Kozubnyak S.A. Splitting of natural frequencies of cylindrical resonator gyro due to non-ideal shape. Vestnik MGTU im. N.E. Baumana. Ser. Priborostroyeniye. = Ser. Instrument Engineering, 2015. no. 3, pp. 39-49. (in Russian). DOI: 10.18698/0236-3933-2015-3-39-49
20. Lankaster P. Teoriya matrits [Theory of matrices]. Moscow, Nauka Publ., 1973. 280 p. (in Russian).
21. Madelung E. Matematicheskiy apparat fiziki [Mathematical apparatus of physics]. Moscow, Fizmatlit Publ., 1949. 618 p. (in Russian).
22. Naraykin O.S., Sorokin F.D., Kozubnyak S.A Natural frequencies splitting of a ring resonator of a solid-state gyroscope, caused by shape imperfection. Vestnik MGTU im. N.E. Baumana. Ser. Mashinostroyeniye = Ser. Mechanical Engineering, 2012, no. 6, pp. 176-185. (in Russian).
23. Golovanov A.I., Tyuleneva O.N., Shigabutdinov A.F. Metod konechnykh elementov v statike i dinamike tonkostennykh konstruktsiy [Finite element method in statics and dynamics of thin-walled structures]. Moscow, Fizmatlit Publ., 2006. 392 p. (in Russian).
24. Matveyev V.A., Basarab M.A., Lunin B.S. Approximation of density distribution of solid-state wave gyro resonator with respect to measured disbalance parameters. Pribory i sistemy. Upravleniye, kontrol', diagnostika = Instruments and Systems: Monitoring, Control, and Diagnostics, 2015, no. 10, pp. 9-16. (in Russian).
25. Basarab M.A., Lunin B.S., Matveyev V.A., Chumankin Ye.A. HRG resonator balancing by chemical etching. Giroskopiya i navigatsiya, 2015, vol. 88, no. 1, pp.61-70. (in Russian).
26. Basarab M.A., Lunin B.S., Matveyev V.A., Chumankin Ye.A. Static balancing of cylindrical resonators of wave solid-state gyroscopes. Giroskopiya i navigatsiya, 2014, vol. 85, no. 2, pp. 43-51. (in Russian).
27. Chang Chia-Ou, Chang Guo-En, Chou Chan-Shin, et al. In-plane free vibration of a single-crystal silicon ring. Int. Journal of Solids and Structures, 2008, vol. 45, pp. 6114-6132. DOI: 10.1016/i.ij solstr.2008.07.033
