Метод расчета колебаний оболочек, содержащих жидкость Текст научной статьи по специальности «Физика»

УДК 533.6: 624.074.4

МЕТОД РАСЧЕТА КОЛЕБАНИЙ ОБОЛОЧЕК, СОДЕРЖАЩИХ ЖИДКОСТЬ © 2002 г. Д.В. Щитов, А.С. Юдин

The numerically-analytical method for solution of a problem about oscillation's of prestressed shells of tanks containing fluid is introduced. With the help of approximations on circumferential and axial coordinates the source problem is resulted in unidimen-

sional boundary value problems for systems of ordinary

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

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

Исходные соотношения. Жидкость. Вводится потенциал скорости Ф, связанный с вектором скорости частиц жидкости: v = -УФ = - gradO. Динамическое давление жидкости р = ржФ,( , где рж - плотность жидкости в состоянии покоя; t=d(...)/dt. В линейном приближении для идеальной жидкости потенциал скорости Ф и давление р должны удовлетворять волновому уравнению АФ - Ф,п/сж2 = 0 (или Ар - р,п/сж2 = 0), где сж - скорость распространения звука в жидкости; Д= V2 - оператор Лапласа (лапласиан); внутренние источники звука отсутствуют. Операторы V и А записывается в криволинейной системе координат, наиболее подходящей для области, занимаемой жидкостью. Так, в цилиндрической системе координат г, 9, z они имеют вид:

V(..) = (.. .),r-er+[(.. ,),o/r]-e0+(.. ,),2-ez,

Д(...) = (...),„+ (...)./r + (...),09+ (•••).гг . где er, е0, ez - орты цилиндрической системы. На границе объема жидкости нормальная скорость выражается через производную от потенциала скорости по нормали граничной поверхности: vn= vn = - Ф,„ = = —п- УФ = w,t, где п — внешняя нормаль; w - нормальная компонента вектора перемещений.

В случае гармонических колебаний зависимость от времени всех функций задачи по методу комплексных амплитуд задается в виде exp(-is0cot ), где i - мнимая единица, s0 = ± 1; со - круговая частота; t- время. Тогда {Ф, р, v, w,...} =

= {Ф,р, v, w, —}ехр(—is0tot), где чертой сверху обозначены амплитуды; p = -s0iupJICO. Волновое уравнение переходит в уравнение Гельмгольца:

equations with independent radial coordinate.

АФ+к2Ф=0, Ар + k2p = 0 , где к = а)/сж. На границе объема жидкости Ф,п = soitow, поэтому амплитуды нормального смещения и давления связаны формулой р,„ =-со2ржш.

Введем характерные параметры Rx, hx - большой и малый размеры; рх - плотность; Ех — модуль упругости; vx - коэффициент Пуассона, а также величины: Ei= hx/Rx - малый параметр, v0=l-vx2, cx=[Ex/(pxv0)]1/2- характерная скорость звука. Перейдем к безразмерным величинам и операторам по формулам:

Ф = Ф/(аЖД А = A-Rx2, V = V Rx, n = n/Rx, w = w/hx, p = pRx2Vo/(Exhx2),

рж = рж./рх, сж = сж./сх, Q = coRx./cx, k = Ш сж. Эти соотношения относятся как к исходным величинам, так и их амплитудам при гармонических колебаниях. Ниже будем работать с амплитудами в безразмерной форме, сохраняя за ними обозначения исходных величин (т.е. для упрощения записи убираются тильды и черточки). Для жидкости основные соотношения в безразмерной форме получают вид:

’ АФ + к2Ф = 0, Ар + к2р = 0, р= - 801рж(й/е1)2Ф,

Ф,п=501е^,р,п=Е1рж(О/Е1)^. (1)

Оболочка. Базовые уравнения малых гармонических колебаний осесимметрично напряженных оболочек вращения включают кинематические, соотношения (2), уравнения движения в усилиях и моментах (3), соотношения упругости (4). Ограничимся вариантом соотношений упругости для ортотропных оболочек. В записи относительно амплитуд уравнения имеют вид:

Еп =u'+k1w + di°^1,

Е22 = V* +tyU + k2W +$2 #2’

Oj = v', fi2=u*-\|/v,

#1 =-w'+kiU, t?2=-w*+k2V,

Ej2 = Q| + £^2 +^2^1 *

К j i K.22 = г?2 + Ф >

Tj =#2, T2 =#* —

Ki2 = У(Ь +к2П1) + (1-у)(т1 +t2)/2; (2)

T,'i + V(TU -T22) + S* +k,(Qn + yH*) +

+ р1ш2и+р1 =0,

S' + T22* + 2v|/(S + yk.H) + k2 (Q22 + 7H') + -

+ р1ш2у + р2 = О,

Сц+УРп+Фгг — к1Т], -к2Т22 +

+ р^иг + рз -р = 0,

+¥(М! —М22) + Н’ -Т,0^ --Тц^-Б0^ =°- (3)

Н'+2\)/\|/ +М22* -<322 -Т221?2 - Б0!?! — Бт?” =0,

Тц = ВцЕц + В12Е22, (1^2)

Мц =ВцКц + (1<->2)

Б =ВззЕ]2, Н = А33Е12. (4)

где (...) = Э(...)/(А1Эа1), (...)* =а(...)/(А2Эа2);

\|/ = А27 А2; аь а2 = 0 - продольная и окружная криволинейные координаты на основной поверхности оболочки; Аь А2 - коэффициенты Лямэ; у = 0, 1 - переключатель между вариантами уравнений, отличающихся выражением для кручения [1]; р - динамическое давление (реакция) жидкости. Величины с ноликами вверху соответствуют компонентам радиально-симметричного статического НДС. Они определяются решениями в общем случае нелинейных краевых задач [2].

В уравнениях (2) - (4) окружная координата отделяется представлением амплитуд компонент динамического НДС тригонометрическими рядами Фурье:

£ [^(1)псО8П0С2 +Ш(_1)п81ППа2] , п=0

Т11 = £ ЕГ11(1)пСОЗпос2 +Тц(_1)п5тпа2] , п=0

и т. д. После перехода к безразмерным величинам формируется каноническая система уравнений относительно основных функций У^ ') = 0,..., 7:

Уо =3(_5) +(1 + У)£1к2Н(_8), У! =МП(8),

^2 =Тп(8)> = СЬод +е15пН(-5)>

у4 = У(-8), У5 = 1?1(8), У6 = и(1), У7 = W(s),

у;= и/аь п, 5, У), У={Уо, У,,..., У7},

] -0,..., 7, где п = п/А2 ,

Р„ = -2 у У0 + 5пТ22(8) +(1-у>|/е1(к1-

-к2)Н(_8) -к2022 —Р]П2У4 — р2,

Р. = У(М22Ы-%)-28пН(.5) + У3/г, +

+Т°У5 + У2#°,

Р2 = ¥(Т22(5) - У2)-бпУ0 -к[У3 + х

х[2укг+ (1 — у) (^ + к2)]-риО!У6 -р„

= — ¥^3 + к] У2 + к2Т22(8) —бп((322 +

+ 281уН(_5))-р1^2У7 -р3 +р,

^4 =Е12(-Ю +¥^4 +^6-СА°&2Ы)’

^5 = Кц(5)/е1 •

Е6=ЕШ5)-к1У7-81г?1°У5(

Р^-У^к^.

Входящие в правые части уравнений величины определяются последовательностью формул:

$2(-5) = к2У4 +8пУ7,

С, =(1-у)¥к1У4 +5п{[(1-у)к1 +к2]У6 -- (2 - у)У5 - 1(/У7 } - !?2(_8) [(1 - у)V+к 2 8,^ ],

1?! = в33 +(1+Зу;е12к2о33, ь2 =(1+У)2 61 к2033, Ь3=ВпОп,

= (У1 ~Ь2С1)/Ь1,

^12(-б) ~ ег (к2Е12(-5) +С1)/(2-у),

Е2ад = &ТТУ4 +\уУб +к2У7>

К 22(8) =е1(^2(-з) +¥у5)-Кц(<,) = (У,-о12к22(5))вп/ьз ,

Ещя) =(^2 ~®12Е22(5))/Вц ,

^22(5) — О.гКщз) + О22К22(5) > т22(в) = в12е11(5) + в 22 е 22(5) .

Н(-з) =2Е>ззК12(_5).

Q22 = —е 1 [впМ 22(в) + Т22#2(-5) +

+ ^в(У0-е1к2(1 + 7)НН))].

Система зависит от номеров окружных гармоник (мод). Индекс п для упрощения записи опущен. Общее решение определяется суперпозицией решений для отдельных мод.

Каноническая система уравнений нелинейной ’ осесимметричной деформации в квадратичном приближении после перехода к безразмерным величинам И введения ОСНОВНЫХ функций Уо = Тц, У1= <3и, У2= Мц, Уз = и, у4 = \У, у5 = -§1 имеет вид (нолики для упрощения записи опущены;] = 0,..., 5):' У]=^(«1.У)» У = {Уо.-»У5}» (5)

*о =У(Т22-Уо)-ЦУ1-Ч1.

^ =-¥У1+к1У0 + к2Т22-Чз>

{2 =¥(М22-у2)+у,/е1+у0у5,

^ =Еп-к1у4-е1у^/2,

Ь =к1Уз-У5. ^5 = Кп/е1 ;

Е22=уу3+ку4, К22=е,\|/у5,

Кп = (у2 —В12К22)ВП/Ь3,

Ец = (Уо ~®12Е22)/В11,

Т22 = В12ЕИ + В22Е22,

М22 = 012Кц + В22К22. (6)

Решение краевых задач для уравнений (5), (6) и исследование статического НДС выполнялось ранее [3,4].

Объект и модель. Рассмотрим жидкость в цилиндрической емкости, стоящей на горизонтальной плоскости (рис.1).

Считаем, что емкость изготовлена из изотропного листового материала толщиной Ь, с модулем Юнга Е, коэффициентом Пуассона V. При транспортировке могут возбуждаться вертикальные, горизонтальные и поворотные колебания. Если известны данные для усредненных вибрационных характеристик разных типов дорог [5], то можно считать заданными амплитуды ускорений (скоростей, смещений) на опорном кольцевом контуре, т.е. кинематические воздействия. Поскольку днище, имея составную конфигурацию, по существу является пологой оболочкой [3], то условия на нижней границе для жидкости можно перенести на плоскость. Таким образом, жидкость можно рассматривать находящейся в цилиндрическом объеме высоты Нж. Обычно заполнение объема всей емкости составляет 95 %, так что Нж близко к высоте бочки.

Формирование разрешающей системы. Уравнение Гельмгольца (1), записанное относительно безразмерной амплитуды динамического давления, в цилиндрической системе координат имеет вид:

Р.гг + Р.,/г + Р.ее/г2 + Р.гг + к2р = 0 (7)

Решение ищется в виде двойного тригонометрического ряда Фурье по координатам 0 и г с коэффициентами, зависящими от г:

N М __

р(г,0,г)=^П рпт( 1,1 )(г>С05 т г+

л=0т=0

+Рпт(1,-1)(г)зт (тг)]со5(п0)+

+[Р„.п(-1.1)0')СО8(тг)+

+Рпт(-1.-1)(г)5т(т2)]5т(п0)]}, (8)

где ш=тл/НЖ) т - число полуволн по оси ъ\ п -число волн в окружном направлении; Нж - высота столба жидкости.

Поскольку в уравнения Гельмгольца входят четные (вторые) производные по г и 0, то после подстановки (8) в (7) получаются однотипные системы для определения рпш(± 1, ± 1)*

рпт(± 1, ± 1)>г г "+■ Рпт(± 1, ± +

"1"Кшп(г)Рпт(± 1, ± 1) — 0, --(9)

где Кпга(г) =Кт- (п/г)2, Кт=к2-т2.

Дальнейшие преобразования и выделение соответствующих решений из (8) связаны с использованием гранично-контактных условий и учетом вариантов вызывающих колебания воздействий.

Условия, соответствующие свободной верхней границе (г = Нж) и совместным колебаниям жидкости и днища (при г = 0), имеют вид: р(г,Нж) = 0, т.е.

1(-1)тРпт(.,(г) = 0; (10)

т-О

М

Р(Г,0) = р(г) = £ Рпт(1)(г); (И)

т=0 4

Р(Г,0), г = - рж(02/£,) \У(Г), ИЛИ

м

4 £т Рпт(-1)(Г)= - Рж (£2 /8!МГ), (12)

ш=0

где w - нормальные прогибы днища. Соотношение (11) определяет динамическую реакцию жидкости, действующую на оболочку.

Рассмотрим вариант вертикальных колебаний, для которых имеются данные по амплитудам вынуждающих ускорений [5]. После перехода к амплитудам смещений усредненные данные для разных типов дорог имеют вид, представленный на рис.2 в логарифмических координатах в диапазоне от 4 до 200 Гц.

Поскольку наиболее податливым и нагруженным элементом емкости является днище, то приемлемой может быть модель сосуда с деформируемым дном и жесткими стенками. Этому варианту соответствует нулевая окружная мода п = 0 и К„т(г) = Кт. Тогда на жесткой цилиндрической стенке (г=а) смещения равны нулю, т.е.

м _ _ 1 Р(а>2),г = X [Рт(1)д(а)со8т2 + Рт(-1),г(а)51пт2.1=0>

т=0

откуда Рт(1) г(а) =0, ргаМ),г(а) = 0, т = 0.......М.

Рис. 2

9і

Нж 2: 10

ф 5 71 2: :0

Г _

Г

Рис.1

Цилиндрическая (полярная) система координат имеет особенность в полюсе. Поэтому в центре ставятся одельные условия на границе малого выреза г = г0. Для варианта жесткой стенки рт(п, г = 0. Рт(-1).г=0 при г = г0. В варианте мягкой стенки Рт<1) = 0, Рт(-1)= 0 при Г = Г0.

Обозначим: рт<±1)= Рот(±1,1). (•••)' = (•••),.•• Уравнения (9) получают вид:

Рш(I) Рш(1) ^ + К-тРш(1) = Рт(-1) Рт(-1)^Г +

+кт рт(-1)= о, ш = 0,..„М.

Обозначив р'1|П) = цга([), запишем уравнения гармонических колебаний жидкости в виде системы стандартной формы:

Рш(1) — Чт(1) = 0, Ят(1) + Ят(1/Г +КтРт(1)= 0|

Рт(-1) — Чш(-1)= 0, Чт(_!)+ Чт(-1/г +

+ ктрт(.1)=0, ш = 0,М. (13)

Используем условия на верхней и нижней границах жидкости. Из (10) следует:

SH) mPm(l)(r) + (-1)М рМ(1) (г) = о, ИЛИ

ш=0

М‘1

Рм<1)(г)= - S(~l) +mpm(l)(r).

(14)

На основе(12)

1

Рмс-1) (г) = ~К^(г)/М- — ЕМРт(-1)(г) (15)

М т=0

где Кп=Нжрж02/(л(,).

С учетом (14), (15) система уравнений (13) преобразуется в следующую:

Рт(1) = Чт(1)> Рм(1)“ ЯМ(1).

Чт(1) ~ — Чт(1/Г — Ктрт(1),

М-1

Чм(1) = ~Чм(1) ^Г + ^мХ(-1) Рт(1)(Г)’

т=0

Рш(-|) = Ят(-1)> Рм(-1) = М(-1>;

Ят(-1) “■ ~ Чт(-1)^Г “ КшРт(-1)> — 0,.. .,М—1.

Чм(1)= “ Чми/г +(КмКп/М)\у +

+(Км/М) £тр■

ш=0

Краевые условия: в центре (г = г0): я, = 0 для жесткой стенки или р, = 0 для мягкой стенки; на краевом контуре (г = а) я, = 0; ]= 0,..., М.

При осесимметричных колебаниях оболочки Уо и У4 равны нулю. В системе уравнений связанных колебаний дна емкости и жидкости основные функции в унифицированных обозначениях 2нумеруются от нуля в следующей последовательности: 2о=Уг, Ъ\ = Уз, Ъг=Х\, 24=У7, =

^6+гп Рт(1)> ^7+М+т— Чт(!)» ^8+2М+ш— РтС-1)? ^9+ЗМ+т —

=Ят(-о> т = 0,...,М. Общее число неизвестных функций (порядок системы) N = 6 + 4(М + 1).

Рассмотрим систему при М=1 (N=14). В этом случае набор основных функций для жидкости будет следующим: 2б = Роа)> = Ркп> — Чо<п> 29 = яко, 2ю=Ро(-1), 2ц=р1(_1), ^!2=Яо(-1)-. 2^1з=Як—и- Соответствующая часть системы для жидкости имеет вид:

7' = 7 7'—7 7'-7'

^6 8’ 7 " 9 ’ 8 10 ’

Гд = -29/г-ВД,

2ю=212, Zl2=

7;з=-21з/г+К1Кш24+К1. (16)

Аналогично для варианта

М=2 (N=18): Z6= р0(1), Ъ~1- рщ), Z8 = Яго)!

29 = Чо(1), ящ), Zu= ц2(1); 2^= Рои).

213= Рк-п> 2(4= р2(-1); 215= Чо(-1). 2^= Чк-1).

2п — Чг<-1) > 26= 29, 27=, 28 = 2П,

29 = -29/г—K0Z6, Z10 = — 210/г— К,27,

2ц = -2,,/г-к0(г6 -27) г12=г15,

213=216, 7<и=Ъ\т, 215 =-215/г-Ко212,

216 = — Z16/r—К^13 ЪХ1 — Z17/r +

+К2К2йг4+к2г11/2. (17)

Динамическое давление жидкости действует на оболочку через правую часть уравнения, содержащего нормальную компоненту нагрузки:

м м

Р4=...-р,£2 \У + £рД1),или Р4=...-р1£22г4 + £г^6.

1=0 j=0

Проанализируем систему (16) первого приближения. Все ее правые части не зависят от 2ц, а два уравнения с правыми частями под номерами 10 и 12 образуют независимую подсистему. Исключая эти уравнения и уравнение Z^, = 2ц, .получим основную систему, определяющую колебания оболочки с учетом динамической реакции жидкости. Аналогичный анализ и формирование основной системы необходимо выполнять и в более высоких приближениях. Так, для системы (17) второго приближения автономную подсистему составляют уравнения 2;2=215, г^б, ^-г.з/г-ВДи,

216= -^б/г-К^з, а Ъ14 отсутствует в правых частях уравнений и может быть определено после решения основной задачи. Отделение этих уравнений от главной системы обеспечивает эффективность и устойчивость счета последней.

После перенумерации оставшихся основных функций и подключения уравнений (5), (6) разрешающая система в случае М = 1 имеет вид: у'о= [(V-1 )у0 + В (1 -У2)(у3+Ф0у4)/г]/г- Р1^2у3, у\= (к|+\|к2)уо+к2В(1-\'2)(уз/г+к2у4)—

- Р1^2У4+Уб+У7>

у,2=0О1Уо+^-1)у2/г+[О(1-у2)е,/г2+Т°11]у5, у'з = уо/В-уу3/г-(к, +ук2)у4-е, 0° 1 у5,

У'4=У5. у'5 = У2/(£1В)-Уу3/г, у'б — У8> у'т=у9,

у'ъ = КоУб -у8/г, у'д = К,у7-у9/г,

у'ю= К^Кд^М—ую/г, (18)

где у3=г)^=0...9;уш=213.

Краевые задачи решались численными методами двусторонней пристрелки, а также дифференциальной прогонки с ортогонализацией по С.К. Годунову [2]. Оказалось, что в актуальном диапазоне частот

(до к=0,12, или f < 200Гц) при возбуждении колебаний вертикальными смещениями опорного контура сходимость обеспечивается уже системой (18).

На рис. 3 показана динамическая интенсивность напряжений для частоты Г = 100 Гц на срединной поверхности и лицевых поверхностях днища при единичной безразмерной амплитуде вертикальных

Рис. 3

В силу линейности задачи уровни, соответствующие рис. 2, для всех характеристик получаются пропорциональным пересчетом. Аналогично трансформируются амплитудно-частотные характеристики.

Работа выполнена при финансовой поддержке гранта РФФИ 00-15-96087, 01-01-00681.

Литература

1 .Юдин А.С., Рукина Т.И., Шевченко В.И. // Изв. СКНЦ ВШ. Естеств. науки. 1981. № 3. С.32-36.

2. Кармишин А.В. и др. Статика и динамика тонкостенных оболочечных конструкций. М., 1975.

Ъ.Юдин А.С., Сафроненко В Г., Щитов Д.В. // Современные проблемы механики сплошной среды. Тр. VI Между-нар. конф. Ростов н/Д, 2001. Т.2. С.158-161.

4. Щитов Д В. Н Тр. аспирантов и соискателей РГУ. Т. 7. Ростов н/Д, 2001. С.17-19.

5. Сафроненко В.Г. и др. И Современные проблемы механики сплошной среды. Тр. V Междунар. конф. Ростов н/Д, 2000. Т.2. С.161-163.

НИИ механики и прикладной математики РГУ

20 июня 2002 г.