CC BY
38
Известия Тульского государственного университета Естественные науки. 2013. Вып. 1. С. 56-66
Механика
УДК 519.632.4
Численное исследование свободных колебаний упругого тела вращения
С. Д. Алгазин
Аннотация. Рассматривается трёхмерная задача о вычислении свободных колебаний первой краевой задачи теории упругости в теле вращения. На доступной для вычислений сетке из 900 узлов найдены собственные частоты, совпадающие с одномерным тестом с 3-7 знаками после запятой.
Ключевые слова: свободные колебания, численный алгоритм без насыщения, трёхмерная теория упругости.
Введение. Рассмотрим векторное уравнение свободных колебаний теории упругости для однородной изотропной среды:
где ш = (1 — 2а)-1 и а — постоянная Пуассона, р — плотность, от — частота колебания, ц — модуль сдвига. Ставится задача об исследовании спектра оператора в левой части уравнения (1) при краевых условиях первой задачи:
Здесь Q — тело вращения вокруг оси (O,x3), x = (xi,x2,x3),0 — его меридиональное сечение.
При всех конечных значениях ш, кроме ш = —1, оператор А* = = А + + ш grad div эллиптичен. Следовательно, решения первой краевой задачи теории упругости достаточно гладкие. Чтобы воспользоваться этой гладкостью, ниже построен метод дискретизации первой краевой задачи теории упругости, не имеющий насыщения [1, 2]. Наиболее распространённым в настоящее время методом решения задач механики деформируемого твёрдого тела является метод конечных элементов. Его недостатки общеизвестны: аппроксимируя перемещение кусочно-линейной функцией, мы получаем, что напряжения разрывные. Вместе с тем, следует заметить, что большинство задач механики деформируемого твёрдого тела описывается уравнениями эллиптического типа, которые имеют гладкие решения. Представляется актуальным разработать алгоритмы, которые учитывали бы эту гладкость.
Аи + и grad div и + Л2и = 0, Л2 = ,х Є Q,
(1)
(2)
Идея таких алгоритмов принадлежит К.И.Бабенко [2]. Эта идея высказана им в начале 70-х годов прошлого века. Многолетнее применение этой методики в эллиптических задачах на собственные значения автором настоящей работы доказало их высокую эффективность. Например, рассматривалась задача на собственные значения для нулевого уравнения Бесселя, на сетке из 23 узлов первое собственное значение этой задачи определено с 28 знаками после запятой. В отличие от классических разностных методов и метода конечных элементов, где зависимость скорости сходимости от числа узлов сетки степенная, здесь имеем экспоненциальное убывание погрешности.
1. Дискретизация. Введем систему криволинейных координат (гД ф), связанную с декартовыми координатами (хі, Х2, хз) соотношениями
Пусть G — область, получаемая меридиональным сечением тела Q. Выберем функции и и v следующим образом. Пусть ф = ^(z), ф = и + iv, z = r ■ ехр^в) — конформное отображение круга \z\ ^ І на внутренность области G. Удобно считать (r, в, ф) сферическими координатами, тогда соотношения (3) задают отображение шара единичного радиуса на внутренность тела Q. Поверхность шара единичного радиуса переходит при отображении (3) в поверхность тела Q. Тогда краевые условия, заданные на дQ, переносятся на поверхность шара.
Обычно при использовании криволинейных координат уравнения для векторных величин записываются в проекциях на оси собственного базиса, координатные векторы которого направлены по касательным к координатным линиям. Этот базис зависит от координат точки пространства. В данном случае такой подход неудобен, так как отображение (3) теряет однозначность на оси хз (если v = 0, то ф любое). Это вызывает появление особенностей в решении, которые вызваны не существом дела, а «плохой» системой координат. Отметим, что сферическая система координат обладает аналогичным «недостатком». Выход из этого положения следующий: оставим в качестве искомых функций проекции вектора скорости иг (i = І, 2, 3) на оси декартовой системы координат, а независимые переменные xl, x2, хз заменим подстановкой (1.1) на г^,ф. Тогда частные производные по декартовым координатам xi (i = І, 2, З) выразятся через производные по r, в и ф: U(xl,x2,x3) = U(v ■ cosф,v ■ sinф,и):
xl = v(r, в) cos ф, x2 = v(r, в) sin ф, x3 = и^,в).
(3)
<
dU rv'd dU rv'r dU
„ dx3 w2 dr w2 дв ’
где a(r, в) = —m'g/w2, (w2 = и'2в + v'2e), в(r, в) = (І + rdgv'r/w2)/v/e.
^ Оу 1 ди ди 1 Оу
Так как выполняются условия Коши-Римана: — = — —, — = — —,
дг г дв дг г дв
то система координат (г, в, ф) ортогональна и в этой системе координат
лапласиан скалярной функции имеет вид
ДФ
vw2
д_
дг
/ дФ \ д / у дФ \
\У дг ) + дв \г дв )
w2 = (ду/дв)2 + (ди/дв)2.
1 д2Ф
+ у2 дф2
(4)
Тогда вместо задачи (1)—(2) имеем внутреннюю задачу в шаре единичного радиуса. Причем на его границе ставятся нулевые граничные условия. Далее будем считать, что конформное отображение круга единичного радиуса на внутренность области С известно. Заметим, что для численного построения конформного отображения имеются надежные алгоритмы [3].
Для дискретизации лапласиана (4) с однородным краевым условием применим методику, описанную в [1].
Таким образом, получаем дискретный лапласиан в виде ^-матрицы:
2
к=0
Здесь штрих означает, что слагаемое при к = 0 берется с коэффициентом 1/2; знак 0 — кронекерово произведение матриц; Н — матрица размера Ь х Ь с элементами
Нщ =ео8 к 2П^гь ^ , (г,І = 1, 2,..., Ь);
Лк — матрица дискретного оператора, соответствующего дифференциальному оператору
vw2
д_
дг
дФ д
ГОТІ + дві
д ( у дФ \
г дв
к2
- Ф, к = 0,1,...,1
с краевым условием
Ф|
Г=1
0.
(5)
(6)
Для дискретизации дифференциального оператора (5), (6) выберем по в сетку, состоящую из п узлов:
г
г
п, ч (2^ — 1)п
ви = ^(уи + 1), Уи = еой є», Єи =--------2П-, V = 1, 2,--,п,
а также применим интерполяционную формулу
«(«) = . У = > - п),
^ (у - УV)’ ' (7)
5^ = §(0и), V = 1, 2,...,п; Tra(y)=еos(n агеео8(у)).
Первую и вторую производные по 0, входящие в соотношения (5), получим дифференцированием интерполяционной формулы (7).
По г выберем сетку, состоящую из т узлов:
К ^ - 1)п
Ги = ~(г„ + 1), ху =ео8 XV, XV = -------- ---, V = 1, 2,...,т,
2 2т
а также применим интерполяционную формулу
?(Г) = Е , ,Т"(Г)(Г - ^ „, <ь = Ф'»), * = 2г - 1. (8)
*=1 т ~1^ (ГV - 1)(* - ^)
Первую и вторую производные по г, входящие в выражение (5), найдем дифференцированием интерполяционной формулы (8).
Исходная система примет вид:
Д-и(1) + ш^— + )?u(1 = 0, Ди(2) + ш-^— + \2u(2') = 0,
oxj 0x2 (q)
(3) d- 2 (3) - du(1) du(2) du(1)
Ди(3) + ш — + \2u(3) =0, - = ——+ —— + —— ^
ихз dxi 0x2 0x3
, du(1 du(1) 1 , du(1) 0u(2') 0u(2')
^ a cos ф —-------I p cos ф^т--------------------sin ф^-—I a sin ф —-I p sin ф +
dr d- v дф dr d-
1 , du(2) rv'a du(3) rv'r du(3) -
+ -cos ф^~; I 2 2 яа = -■
v дф w2 dr w2 d-
В дискретном виде на сетке (-v, rM, фк), v = 1, 2, ■ ■■, n; л = 1, 2,..., m; k = = 0,1, ■■■, 2l получаем:
(m \ / n \
I] D$1 ^ifc + p»v COs фь [^2 DvX uu^k) -
Ml=1 ) Vi=i )
sin фк( £ Dkt1 + +QW sin фк( £ DSiuSifc) +
Vfci=0 / \^i = 1 /
+p,vsin фк( £ DVVi uViu) + v~cos ф^ £ Dkki uiiO +
\vi = 1 ) \hi =0 )
+
(тр В{т) и(3) \ _ (Г^) [г2) г = Г) \ 2^ а^1к I и2 в — в„ К)1=1 '
(£ »!% Оь)
2 I Г — Г ' і
в — ви
)
Е(в) п(3 I — в ь ^ууг %) I — °")ь ■
В криволинейной системе координат (3) уравнения (9) выглядят так:
+ Л2п(1) — 0,
Ап,(1 + ш Ап(2 + ш
.дв п .дв 1 . дв
а 008 ф — + в 008 ф —----------------81П ф —
дг дв V дф
.дв п .дв 1 дв а 81П ф — + в 81П ф— + - 008 ф — дг дв V дф
+ Л2п(2) — о,
А43 + ш
( го'в\ дв ( п)'г \
\и2 ) дг \ и2 )
дв_ и2 ) дв
+ Л2п(3) — 0.
В дискретном виде они запишутся следующим образом:
Нп(1 + ш
а)и008 <
1
---------8ІП фь
Нп(2) + ш
а) 8іп (
(Е ^в^ь) + в»„ 008 < \)2 = 1 )
( 21
\Ь2=0
ьь2 ви)ь2
+ Л2п(1) — о,
1
+------- 008 фь
(т _ \ /и \
^Ь ) + вV 81п фь Е В(и2 в1*2)ь \ +
)2 = 1 / \^2 = 1 )
( 21 Е
ь2 =0
ЬЬ2 °^)ь2
+ Л2п(2) — о,
Нп(3) + ш
(го^) г — г^ ( Е
в — ви К)2 = 1 '
(—) (^ п(0) в )
\и2І г — г^ 1 ^^2 в^2)Ь
в — в^ ^2 = 7
+ Л2п(3) — 0-
X
С учётом (9) имеем
Hu(1) + ш [b(11u(1) + b(12)u(2) + b(1’3)u(3) ] + \2u(1) = 0,
Hu(2) + ш [b(21)u(1) + b(2,2)u(2) + b(2,3)u(3)] + A2u(1) = 0,
Hu(3) + ш [b(3,1)u(1) + b(32)u(2) + b(3,3)u(3)] + \2u(3) = 0^
Распишем формулы для b(%,j\ i = 1, 2, 3; j = 1, 2, 3. Имеем
b(11) :
m(
.V cos фк I V D(®’r) (a.. V cos фкS.... Skk D(r)
b(1,2)
(m __
E D(M2) ^v cos фк5vviSkkiDMbi) +
M2 = 1
ID^Ip,iv cos фк6kkiD^i - -^Sin фк5VViD^i)
\ vMiV J
(n
E DVV2 (e.v2 cos фкS.Mi Skki D^)
V2 = 1
+DVti (
1 sin фк( £ Dkk( - -j— sin фк2 bvVi b.Mi D^ki UMV \k2=0 V v.V
MMiu kki^ v2vi
V2 = 1
aM Vi cos фк Skki D(;M i - sin фк d.Mi Dk4>l}i
vM Vi
+Dk>ki (aM V cos фkiS v Vi DMr) i + ^ cos фki SMMi D
(m
E DM>M,r) {aM2V sin фкsvVi Skki dSm2)
M2 = 1
ID^I p,i v sin фк Skki D(Vi + -L cos фк Svvi D>)
\ vMiV J
(n
E DVX (pMv2 sinфкsmmiSkkiDV%i)
V2 = 1
+DieJi( aMVi sin фк Skki DMrM)i + cos фк SMM i Dk4ki)
\ vMVi J
- -1 sin фк( £ D<jk>i( vL COS фк2 Svvi Smm i D(£k i) vMV \k2=0 VvMV >
+
+
+Dkki (aMV sin фк i Svvi DMM i + ^ sin фк i Smm i DlVi)
Ъ(1,3):
(
а,іcos фk
V D{e’r) Z_> Г,,2
2=1
(rv\
\w2 )
\
2 ) r = r,, 6vv 1 6kk 1 D(rl 1
2
+
в = ві
/
+DX2)
( rv'r \
\w^J r = r
\
№ 1 в = ві
§kk i DV0l)1
+
+в,„ cos 1
5^ — ( g)
V2 = 1 \
w2 I r = r, ^ 1 kk l V2V1
в = ві
+Dll
іі
( ЧЛ
\w2 J
2 ) r = r, 6kkl DWl
в = в„_
))
§,, i §kki DH \\
J/
+
___— sin фkd(ф) < ( ) § D(r) — ( —2 ) D(0) 6
v,V фDkki]\ w2 ) r = r, §VV1 D,, 1 ^ w2 ) r = r, Dvvi6,,i
в = ві
в = ві
Ъ(2,1):
а,і sin
(т __
5-/ Г\ч,2 (уа,2і cos <Pk§W1 §kki D,2,1) + ,2 = 1
в
+Df,rM вп v cos фk6kki D^ — — sin фk§wi Dkk ) ) +
(ф)
J,1і
wi^ kki
(n
Yj Г(і2 (ві cos ^ §,,1 §kki d(2ii) + і2=1
+ DH1 ^ а,і1 cos ^ §kki D<£, 1 —
(Е Гф (
\k2=0 V
+-------cos (
D
,і
(ф)
DZH — — sin фk2 §ni §,,i D%ki ) +
v,y
)
kk1 (а,і cos фki §іі1 D(jl,i + в,„ cos фki §,,i Dffi) )
2
Ъ,2>2)
Ъ,2>3)
1
+-------
У,„
а,„ 81П
(т _
5-/ В№2 {^,2^ 81П фк5ии1 5к*1 +
,2 = 1
+ВЙ? (V В1Пфк6кк1 ВУ) + -±- С08фк5ищВ^Л) + \ 7,1^ / /
(и
Е В^1^2 (Р,»2 Й1П фк5,,1 5кк1 +
V2 = 1
+ВИ\( а,»1 ®1п фк 5кк1 в,г,1 + сое фк 5,,1 Вкк0) +
\ у,^1 / /
+ С08 ф^ ^ Вкфк2 ^ 7^ 008 фк2 5™15,,1 +
\к2=0 У
+Вк11 (а,^ зт фк15^10$1 + ^ ®1П фк15,,1 ,
(
а,и ®1П фк
V в{в’г) /_> В,,2
,2 = 1
( Гув\ \ы2 )
\
2 ) г = Г,, 5ии15кк1 В,2,1
0 = 0„
+ВХ:)
/ гь^ )
\
5кк1 В^1
У}2 } Г = Г,1 "кк1~ ™1
0 = 0у
+
))
+
/
+в^ 8Ш<
и
В,в) — (5 5 в,в)
/ , ВУУ2. I ш2 I Г = Г, 5,,15кк1 В^2^1
'2 = 1 \ 0 = 0и.
+в^1
( гул
\ы2 )
\
5кк10,
ы2 ) г = г, кк1 ^
+
))
+
)
гое фкБ,Ф) < () 5 Б,г) — () в,в) 5
ШЬ фкВкк1\\ ы2 ) Г = Г, 0ии1 В,,1 I ы2 } Г = Г, В™15
0 = 0и
щ2 у Г = Г, ии1 ^1
0 = 0„
1
Ъ(зл):
(f \ f т _
w?) r = rjE D,2 (<W cos Фk Sn,, Skki Г£п)
в = ві ',2 = 1
+
+Dtr)(в,г v cos фkUkki D^ — ^ sin фkUni Dkki
,1і
(w?) r = r, (E Гіі2 {в,і, cos <kU,Ukk,Гі2і1) +
в = ві ^і2=1
+Гі% (а,і,cos ^Ukki — -^sin ^U,,i Dk4>k1
\ и,і1
Ъ(3.2):
(f \ / т _
wvi) г = гЛ Е ГІЙ2) (а,2іsinФ^іА* <>„) +
в = ві ^,2=1
+Г(,:)(в,1 v sin ФkUkki Di% + ^ cos ФkUni D^i
,i і
( w2) r=(e гіі2 {в,„ sin a U,,i Ukki D>,%i)+
в = ві ^і2=1
+Гіі\ (а,иsln ^ Ukki cos ^ U,,i Dkkk1j),
\ и,і1 / /
Ъ(3.3):
/ rv0\ \w2 )
О І гу» ---------- /V»
w2 r = r, в = ві
т
V D(0’r) Z_> ,,2
,2 = 1
( 4 )
\w2 )
\
U Uuu D(r) Uiil Ukki dU2u
w2 ) r = r,2 іі1 kkl ,2,i
в = ві
+
/
/
+D(0’r)
+D,,i
/ rv2 \ VW у
О І гу» ---------- /у»
w2 r = r i в = ві
Ukki Гіі1
и
( rv2 \
\w )
0 I fy* — fy*
w2 J ' =',
в = ві
n
V D0
іі2
і2 =1
( rv’r \
\w)
\
w2 I Г = Г,, U,,l Ukkl Гі2і1 w / ,
в = ві„
+
+Б(А)
1
[ ( )
\иі2 )
\
5ккі
р(г)
„1
))
Здесь Б^), Б(г), В(в) и — матрицы численного дифференцирования по г без удовлетворения краевым условиям, по г с нулевым краевым условием при г — 1, по 9 и ф; 5 — дельта Кронекера.
2. Результаты численных расчётов. Расчёты проводились для шара и области близкой к шару с эпитрохоидой в меридиональном сечении (для задач 2 и 3 конформное отображение задаётся формулой ф(х) — г(1 + єгПр), где пр — 4, є — 1/6 для задачи 2? и пр — 12, є — 1/16 для задачи 3), на сетке из 900 — 10 х 10 х 9 узлов. Таким образом, размер решаемой спектральной задачи 2700 х 2700. Результаты расчётов представлены в таблице. В первой колонке таблицы указан номер действительного собственного значения для шара, во второй колонке приведены простые собственные значения для шара, удержано количество знаков, совпавшее с одномерным тестом [4, стр. 720], в 3-й и 4-й колонках приведены простые собственные значения для возмущённого шара, а — 0.25.
Результаты расчёта простых частот X2 для шара и областей, близких к шару, а — 0.25
№ Шар Эпитрохоида є — 1/6, пр — 4 Эпитрохоида є — 1/16, пр — 12
4 20,19072854 18,4696406693752 20,5137012323477
30 59,6795 58,6037971034038 54,5883426026319
79 118,902 118,350831138150 109,895709114932
136 197,852 192,040762581846 194,761657499711
Выводы. Таким образом, на сетке из 900 узлов возможно определить первые собственные частоты шара и области близкой к шару с 3-7 знаками после запятой. Такой эффект достигнут применением метода без насыщения для решения спектральной задачи.
1
Список литературы
1. Алгазин С.Д. Численные алгоритмы классической математической физики. М.: Диалог-МИФИ, 2010. 240 с.
2. Бабенко К.И. Основы численного анализа. М.: Наука, 1986. 744 с.; Издание второе, исправленное и дополненное, под редакцией А.Д. Брюно. Москва-Ижевск, РХД, 2002. 847 с.
3. Казанджан Э.П. Об одном численном методе конформного отображения односвязных областей. М.: Ин-т прикладной математики АН СССР, 1977. Препр. № 82. 60 с.
4. Новацкий В. Теория упругости. М.: Мир, 1975. 872 с.
Алгазин Сергей Дмитриевич (algazinsd@mail.ru), д.ф.-м.н., ведущий научный сотрудник, Институт проблем механики им. А.Ю. Ишлинского РАН, Москва.
Numerical examination of the free oscillations of an elastic body
gyrations
S. D. Algazin
Abstract. The three-dimensional problem about an evaluation of free oscillations of the first boundary value problem of the theory of elasticity in a solid of revolution is considered. On a grid accessible to evaluations from 900 knots outcomes are received: the discovered fundamental frequencies coincide with the one-dimensional test with 3-7 signs after a comma.
Keywords: free oscillations, numerical algorithm without saturation, the three-dimensional theory of elasticity.
Algazin Sergey (algazinsd@mail.ru), doctor of physical and mathematical sciences, leading researcher, Institute for Problems in Mechanics of RAS, Moscow.
Поступила 08.11.2012
