АНАЛИЗ РЕЗОНАНСНЫХ СВОЙСТВ ТОНКОСТЕННЫХ ЭЛЕМЕНТОВ МУЗЫКАЛЬНЫХ СТРУННЫХ ИНСТРУМЕНТОВ
Шлычков С. В.
Марийский Государственный Технический Университет (Йошкар-Ола)
1.Введение
Качество музыкальных струнных инструментов в значительной мере определяется на стадии проектирования и во многом зависит от степени совершенства расчетных моделей и методов. В настоящее время в практике проектирования широкое распространение получили упрощенные расчетные модели и эмпирические формулы [1,2]. Вместе с тем в последнее время появились разработки, в которых сделаны попытки построения расчетных моделей высокого уровня на базе современных численных методов. В [3] представлена конечноэлементная динамическая модель корпуса скрипки. В наших работах [4,5] разрабатывается уточненная конечноэлементная модель классической гитары (рис.1). В отличие от скрипки для гитары существенным фактором расчетной схемы становятся нелинейности, связанные с начальным натяжением струн.
Известно, что основными звукоизлучающими элементами гитары являются верхняя и нижняя деки. Дека представляет собой тонкостенную, как правило, деревянную панель сложной геометрической формы, нагруженную силами натяжения струн. В качестве расчетной схемы деки применяется предварительно напряженная тонкостенная пластинка. Древесина считается как ортотропный материал. Для дискретизации пластинки используется треугольный конечный элемент (КЭ) [6]. Учитываются сдвиговые деформации и инерция поворота нормали.
В настоящей работе рассматриваются задачи динамики пластинок, выполняется расчет собственных частот и форм. Результаты расчетов методом конечных элементов (МКЭ) сопоставляются с данными известных аналитических решений [7,8]. Исследуется точность и сходимость решений.
2.Основные расчетные соотношения
Задача динамики предварительно напряженной пластинки описывается однородной системой дифференциальных уравнений вида
[м ш+[к- С ы= 0 (2.1)
и системой линейных алгебраических уравнений
[К]^}^}. (2.2)
Здесь, [М], [К] и [О] - матрицы масс, жесткости и начальных напряжений (геометрической жесткости) конструкции; {¿1(0}, {¿(0}, {¿т } - векторы обобщенных ускорений и перемещений соответственно.
Дифференциальное уравнение (2.1) описывает свободные колебания предварительно напряженной динамической системы относительно равновесной формы (конфигурации) (2.2). Напряженно-деформированное состояние (НДС), соответствующее равновесной форме, определяется силами натяжения струн {Рт}. Влияние мембранных усилий на изгибную жесткость деки учитывается при помощи матрицы начальных напряжений [О].
Решение (2.1) находится в виде
ШН"*^}. (2.3)
В этом случае задача динамики сводится к алгебраической проблеме на собственные значения
[К]{^} =^2[М]{^}, (2.4)
где и - соответственно собственные частоты и векторы собственных форм колебаний 0=1, 2,..., р).
Ограничимся расчетом р низших собственных форм и частот, причем p<<m, где m - порядок разрешающих уравнений. В этом случае одним из наиболее эффективных численных методов расчета является метод итераций в подпространстве собственных векторов [9].
Согласно алгоритма метода итераций в подпространстве предварительно строится начальное подпространство с матрицей [x0] размерности (p x m). Для ускорения сходимости при формировании матрицы [x0] используются единичные векторы. Причем, одновременно считается не r, а p = min {2r, r+8} собственных векторов. Это позволяет при вычислении наибольшего собственного значения А обеспечить заданную точность (до 6 значащих цифр) не более, чем за 9 итераций. Решение (2.4) сводится к процедуре обратных итераций для p собственных векторов
[У0] = [M][x0]
K
{x к }={y кчН
{x к }={x к }[0к]
{у к }=[M]{x к }
(2.5)
с одновременным проектированием матриц на подпространство.
K к
{x к }
K
{x к }
M к
= {x к }T № к }.
(2.6)
и решением задачи на собственные значения подпространства
K к
[Q к ] =
M к
[Q к ][л к ],
(2.7)
где
K
= [K-G], [Л] - диагональная матрица с элементами Ai=Wj2, к - номер
итерации.
В результате преобразований (2.6) имеет место редукция (понижение размерности) матриц с (т х т) до (р х р). Для решения (2.7) используется стан-
дартная программа метода Якоби. Критерием сходимости первых собственных частот является условие (А/к)- А/к-1))/ А/к) <е, где е =10-6 - заданная точность решения.
Построение расчетной модели МКЭ начинается с дискретизации. Дискретный аналог тонкой пластинки представляет ансамбль плоских треугольных КЭ (рис.1).
3
Хз
2
X!
Рис.1
Расчетные соотношения получаются на основе смешанной вариационной формулировки [6]. Для аппроксимации перемещений КЭ используется выражение вида
М^Ф]^}. (2.8)
Деформации аппроксимируются зависимостью
{е}=[ю]{а}. (2.9)
Здесь [Ф] и [ю] -матрицы функций формы, {¿} - вектор обобщенных узловых перемещений КЭ, {а} - вектор произвольных постоянных.Вектор перемещений
{^={1 ^2,ю, 01,62}Т, вектор деформаций {е}={еье2,У12,КьК2,%12,¥ь¥2}т Связь деформаций с перемещениями определяется соотношениями типа Коши
{е}=[Ь] {u}, (2.10)
1
где [Ь]- дифференциальный оператор.
Выражение для матрицы жесткости КЭ имеет вид
[К(п)]=[8]т[и]-1[8], (2.11)
где [Б]= Д[со]т МЫф^, [И]= Д[со]т Мю^, п-порядковый номер КЭ.
Матрица масс КЭ определяется выражением
[М(п)]= Д[ф]т [тЦф^Ю,
(2.12)
где [m]=diag[mo, то, то, 11, 12]. Здесь то- погонная масса; 11, 12 - моменты инерции погонной массы относительно координатных осей Х1, Х2. Матрица начальных напряжений определяется формулой
[0(п)]= Ц[Ф]т [о]т [к0 Ь]№, (2.13)
где
N0
N0 к12
_К°2 N2 J
-матрица начальных усилий, которые определяются на
основании решения статической задачи (2.2); [0]-матрица дифференциальных операторов, имеющая следующий вид
Э
[0]=
00 00
Эх 1 Э
дх2
00 00
Отметим, что матрицы [К(п)], [М(п)], [0(п)] являются согласованными, они составляются на основании единых функций формы. В результате матрица масс, матрица жесткости и матрица начальных напряжений КЭ имеют одинаковые размерности (30х30) и блочную структуру заполнения.
З.Анализ точности и сходимости решений МКЭ
На базе КЭ [6] и алгоритма ансамблирования МКЭ разработана программа расчета тонкостенных конструкций ЛБСМ. Для оценки точности расчета ЛБСМ рассмотрим ряд тестовых задач. Используем модели изотропного и ор-
примем следующие упругие постоянные: модуль упругости Е=198 ГПа, коэффициент Пуассона У=0,3. Для ортотропного тела: Е1=19,8 ГПа, Е2=198 ГПа, Vl2=0,03, У21=0,3, 012=7 ГПа, ^3=19,6 ГПа, ©23=19,6 ГПа. Рассмотрим квадратную пластинку с параметрами: сторона a=0,4 м, толщина h=0,01 м. Пластинка разбивается на N=256 КЭ.
• Изотропный материал.
1. Контур пластинки, защемлен (рис.2). В этом случае частота собственных колебаний пластинки определяется формулой [7]:
где Э-цилиндрическая жесткость, а-безразмерный коэффициент.
В табл. 1 представлены коэффициенты а для десяти низших собственных частот, рассчитанных с помощью программы ЛБСМ. Для сравнения приведены
3 3
тотропного тела. Считаем, что плотность р=7,8-10 кг/м . Для изотропного тела
(3.1)
a
Рис. 2
результаты, полученные асимптотическим методом (АМ) [7], и решения [8]. Решение [8] получено в двойных тригонометрических рядах, удовлетворяющих граничным условиям на контуре. Учитывались первые шесть членов ряда. Величины т и п - целые числа, которые определяют форму колебаний.
Таблица 1
т п ЛБСМ АМ [7] Решение [8]
1 1 3,647 3,556 3,646
1 2 7,452 7,386 7,437
2 1 7,452 7,386 7,437
2 2 10,941 10,889 10,965
1 3 13,420 13,337 13,393
3 1 13,485 13,337 13,393
2 3 16,695 16,656 16,717
3 2 16,695 16,656 16,717
3 3 22,146 22,222 —
1 4 21,615 21,313 —
4 1 21,615 21,313 —
2 4 24,621 24,540 24,631
4 2 24,735 24,540 —
3 4 29,838 29,960 —
4 3 29,838 29,960 —
4 4 37,177 37,556 —
Из таблицы видно, что при N=256 решение МКЭ для основной формы
колебаний отличается от аналитического решения [8] менее, чем на 0,03%. Для остальных частот отличие составляет менее 0,1%. Это отвечает требованиям инженерной точности расчетов. Заметим, что для рассмотренной задачи полученное решение МКЭ дает более близкие к [8] значения, чем решение АМ [7].
2. Сторона Х1= а шарнирно оперта, остальные - жестко защемлены по контуру (рис.3).
3. Стороны Х2 = 0 и Х2 =a шарнирно оперты, две другие - жестко защемлены по контуру (рис.4).
В табл.2 приведены значения коэффициентов а, вычисленные с помощью
программы ЛБСМ. Здесь же представлены результаты решений АМ [7] и [8]. А Л
Рис. 3 Рис. 4
Таблица 2
т п Вариант 2 (рис.3) Вариант 3 (рис.4)
ЛБСМ ЛМ[7] Решение [8] ЛБСМ ЛМ[7] Решение [8]
1 1 3,221 3,178 3,238 2,928 2,292 2,937
1 2 6,411 7,174 — 5,526 5,542 —
2 1 7,215 7,174 — 7,037 5,542 —
2 2 10,171 10,166 — 9,534 9,580 —
1 3 11,819 13,191 — 10,345 10,355 —
3 1 13,297 13,191 — 13,171 10,355 —
2 3 15,315 16,114 — 14,095 14,200 —
3 2 16,131 16,114 — 16,556 14,200 —
3 3 21,033 21,160 — 20,054 20,230 —
Из табл.2 видно, что коэффициенты а для основных частот, найденные с помощью МКЭ, отличаются от аналитического решения [8] менее, чем на 0,6% (вариант 2) и менее, чем на 0,3% (вариант 3). По сравнению с АМ [7] решение МКЭ для рассмотренных задач дает более близкие к [8] значения основной частоты (отличие решения [7] для варианта 3 составляет 20%). Для других форм колебаний значения коэффициентов а (МКЭ) незначительно отличаются от значений, полученных АМ [7].
4. Стороны Х1= а и Х2= а шарнирно оперты, две другие - жестко защемлены (рис.5).
5. Сторона Х1=0 жестко защемлена, остальные - шарнирно оперты (рис.6).
Расчетные коэффициенты а представлены в табл.3.
Рис. 5
Рис. 6
Таблица 3
т п Вариант 4 (рис.5) Вариант 5 (рис.6)
ЛБСМ ЛМ [7] Решение [8] ЛБСМ ЛМ [7] Решение [8]
1 1 2,730 2,722 2,746 2,382 2,395 2,401
1 2 6,124 6,133 — 5,209 5,234 —
2 1 6,147 6,133 — 5,930 5,234 —
2 2 9,344 9,389 — 8,653 8,726 —
1 3 11,636 11,600 — 10,146 10,151 —
3 1 11,636 11,600 — 11,503 10,151 —
2 3 14,687 14,778 — 13,432 13,549 —
3 2 14,717 14,778 — 14,184 13,549 —
3 3 19,861 20,055 — 18,827 19,050 —
Для варианта 4 коэффициент а для основной частоты, найденный с помощью МКЭ, отличается от решения [8] менее, чем на 0,6%. Для варианта 5 погрешность решения МКЭ составляет менее 0,8%. Это соответствует инженерной точности расчета. Для других форм колебаний значения коэффициентов а, рассчитанные МКЭ, достаточно близки к значениям, полученным АМ [7].
• Анизотропный материал.
Рассмотрим задачу о собственных колебаниях квадратной пластинки постоянной толщины. Пластинка защемлена по контуру (рис.2). Главные направ-
ления упругости совмещаем с направлениями сторон. Уравнение свободных колебаний пластинки имеет вид [7]:
БцЭ^/Эх4 + 2В12д^/Эх2Эу2 + В22д^%4 +рЬЭ^/Э12 =0 (3.2)
Здесь ^^х, у, ^ - нормальный прогиб; Э11 и О22 - цилиндрические жесткости для главных направлений упругости; О12 - смешанная жесткость. В этом случае формула для частот собственных колебаний пластинки имеет вид [7]:
„2
..........(3.3)
Ю = ¡2
К
Он
рИ
Значения вычисленных коэффициентов а приведены в табл.4.
Таблица 4
т п ЛБОМ АМ [7]
1 1 7,654 7,570
2 1 10,147 9,971
3 1 15,504 14,936
1 2 19,725 20,137
2 2 21,334 21,587
4 1 23,857 22,249
3 2 25,062 25,060
4 2 31,754 30,898
1 3 37,678 39,064
2 3 38,942 40,300
3 3 41,943 42,965
Из таблицы видно, что при N=256 для низшей собственной частоты решение МКЭ дает достаточно близкие к [7] значения. Максимальное отличие высших частот не превышает 4%.
Исследование устойчивости пластинки.
Для оценки влияния мембранных усилий на спектр собственных частот используется матрица Очевидно, что при достижении сжимающими мембранными усилиями определенного значения, матрица [К^] станет вырожденной. Это условие соответствует потере устойчивости типа дивергенции по первой собственной форме. Таким образом, варьируя нагрузку и определяя
> Ч qx2
qxl
X
; к
1
Рис. 7
спектры собственных частот, можно найти значение критической нагрузки, при которой низшая собственная частота обращается в нуль.
Рассмотрим следующие схемы закрепления и нагружения. 1. Квадратная пластинка, жестко защемленная по контуру и сжатая силами qXl, qX2, причем qX1=qX2 (рис.7). Пластинка изготовлена из изотропного материала: модуль упругости Е=198 ГПа, коэффициент Пуассона У=0,3. Огорона a=0,4 м, толщина И=0,01 м. Пластинка разбивается на N=256 КЭ. Критическая нагрузка определяется формулой [10]:
п Ъ
qкр = qxl = qx2 = У—, (3.2)
a2
где у- безразмерный коэффициент. В [10] приводятся следующие значения коэффициента: у = 5,33 (С.П.Тимошенко) и у=5,30 (Тейлор).
Для исследования сходимости к аналитическому решению эта задача была решена при различном числе КЭ. Расчетные значения критической нагрузки приведены в табл.5.
Таблица 5
Число КЭ qкр,•106 Н/м
4 39,438
16 14,524
64 6,647
256 5,938
Тейлор [10] 5,932
Тимошенко С.П. [10] 5,965
На рие.8 результаты расчета представлены графически.
Якр,106 Н/м 40 -
30
20 -
10 5,932
0 50 100 150 200 250
Число КЭ Рис.8
аналил ическое ре шение
Из графика видно, что с увеличением числа КЭ имеет место сходимость решения МКЭ к аналитическому решению [10].
2. Рассмотрим устойчивость квадратной пластинки, у которой стороны Х2=0 и X2=a шарнирно оперты, а две другие - жестко защемлены. Пластина равномерно сжата в одном направлении силами дХ1 (рис.9). Полученные значения критической нагрузки в зависимости от числа КЭ приведены в табл.6.
А
Х2
Рис. 9
Таблица 6
Число КЭ дкр,-106 Н/м
4 345,312
16 13,391
64 6,498
256 6,158
Аналитическое решение [11] 6,205
Электронный журнал «ИССЛЕДОВАНО В РОССИИ» 93 8 http://zhurnal.ape.relarn.ru/articles/2000/0064.pdf
Якр,106 Н/м 40
30
20 -
10 6,205
0
0 50 100 150 200 250
Число КЭ рис.10
На рис.10 приведен график, доказывающий сходимость полученных результатов к аналитическому решению при увеличении числа КЭ.
4.Расчет собственных частот нижней деки гитары
Рассмотрим нижнюю деку семиструнной гитары артикула Н-34Р, изготовленной на Кунгурской музыкальной фабрике. Она изготовлена из
3
древесины. Примем следующие физические параметры: плотность р=500 кг/м , Е1=15,9 ГПа, Е2=0,69 ГПа, v12=0,44, у21=0,038, 012=0,63 ГПа, 013=0,408 ГПа, а23=0,03 ГПа [12]. Геометрические размеры соответствуют указанному артикулу. Дискретный аналог деки представляет ансамбль плоских треугольных КЭ (рис. 11). Дека разбивается на N=732 КЭ.
Определяем собственные частоты деки при различных условиях закрепления по контуру. Значения пяти низших собственных частот приведены в табл.7.
Рис.11
Таблица 7
Защемление Шарнирное опирание
187,8 Гц 86,6 Гц
258,5 Гц 124,5 Гц
269,1 Гц 148,1 Гц
328,7 Гц 202,4 Гц
397,6 Гц 270,3 Гц
Из таблицы видно, что полученные значения спектров собственных частот существенно отличаются друг от друга в зависимости от принимаемой модели
граничных условий. Это обстоятельство следует учитывать при построении расчетных моделей дек музыкальных струнных инструментов.
5.Заключение
В работе представлена расчетная динамическая модель гитарной деки. Задача динамики сформулирована как задача на собственные колебания предварительно напряженной динамической системы с конечным числом степеней свободы. При построении динамической модели применяются смешанная вариационная формулировка и МКЭ с независимой аппроксимацией перемещений и деформаций. Учитываются сдвиговые деформации и инерция поворота нормали. Эти эффекты являются существенными для тонкостенных элементов музыкальных струнных инструментов, изготовляемых, как правило, из древесины, имеющей малую сдвиговую жесткость.
С целью оценки точности модели рассчитаны спектры низших собственных частот пластинок при различных граничных условиях. При этом рассмотрены модели изотропного и анизотропного тел. Сопоставление результатов расчета с данными известных аналитических решений позволяет сделать вывод о достоверности и достаточно высокой точности модели. Сравнительный анализ результатов, полученных при решении задач устойчивости, также дает уверенность в правильной оценке влияния мембранных усилий на динамические свойства пластинок. Такая оценка необходима для учета влияния натяжения струн на динамические свойства конструкции. Характерная сходимость полученного решения МКЭ к точному решению "сверху" дает основание утверждать: конечно-элементная модель являетсяется более жесткой, чем реальная конструкция. Проведенные ранее исследования [4] свойств КЭ [6] при решении статических задач также сведетельствуют о "работоспособности" модели.
Следовательно, используя программу расчета параметров НДС и спектра
собственных частот элементов тонкостенных конструкций ASCM,
составленную на базе КЭ [6], можно получать достоверные результаты.
Литература.
1. Корсаков Г.С. Технология музыкальных инструментов из древесины: Учебное пособие. - Л.: ЛТА, 1986. - 73с.
2. Римский-Корсаков А.В. Дьяконов Н.А. Музыкальные инструменты: Методы исследований и расчеты. - М.: Местная промышленность, 1952.-345с.
3. J. Bretos, C. Santamaria, J. Alonso Moral. Vibrational patterns and frequency responses of the free plates and box of a violin obtained by finite element analysis// Journal Acoustic Society of America. Vol.105 (3), 1999 - P. 1942-1950.
4. Куликов Ю.А., Шлычков С.В. Компьютерная динамическая модель музыкального струнного инструмента как композитной конструкции// Композиционные материалы в авиастроении и народном хозяйстве: Материалы Всероссийской научно-технической конференции, 5-8 октября 1999 г. - Казань, 1999. Ч.11 - С.36.
5. Шлычков С.В. Оценка точности расчетной модели напряженно - деформированного состояния тонкостенных элементов музыкальных струнных инструментов// Электронный журнал "Исследовано в России". 2000. №18. - C.245 - 262. http://zhurnal.mipt.rssi.ru/articles/2000/018.pdf
6. Попов Б.Г. Расчет многослойных конструкций вариационно-матричными методами. - М.: Издательство МГТУ, 1993. - 294с.
7. Болотин В.В., Макаров Б.П., Мишенков Г.В.и др. Асимптотический метод исследования спектра собственных частот упругих пластинок.// Расчеты на прочность: Сб. статей - М.: Машгиз, 1960. Вып.6. - С.231 - 253.
8. Iguchi S., Die Biegungsschwingungen der vierseitig eingespanten rechteckigen Platten. Ing-Arch. Bd. 8. №1. 1937.
9. Бате К., Вилсон Е. Численные методы анализа и метод конечных элементов. - М.: Стройиздат, 1982. - 448с.
10.Тимошенко С.П. Устойчивость упругих систем. - Гостехиздат, 1955.
11.Алфутов Н.А. Основы расчета на устойчивость упругих систем. - М.: Машиностроение, 1 978. - 263с.
12.Ашкенази Е.К., Ганов Э.В. Анизотропия конструкционных материалов: Справочник - Л.: Машиностроение , 1972. - 216 а