Научная статья на тему 'Фазовые скорости бегущих волн в цилиндрической оболочке при конечной и бесконечной жесткости поперечного сдвига'

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

CC BY
33
7
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
ОБОЛОЧКА / ФАЗОВАЯ СКОРОСТЬ / ПРОДОЛЬНЫЕ ВОЛНЫ / ВОЛНЫ ИЗГИБА / КРИТИЧЕСКАЯ ЧАСТОТА

Аннотация научной статьи по физике, автор научной работы — Каледин Валерий Олегович, Седова Елена Александровна, Шпакова Юлия Владимировна

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

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по физике , автор научной работы — Каледин Валерий Олегович, Седова Елена Александровна, Шпакова Юлия Владимировна

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

The phase velocity of traveling waves in a cylindrical shell with finite and infinite stiffness transverse shear

An analytical solution to the problem of calculation of phase velocity of traveling waves on the surface of a cylindrical shell on one edge which has a perturbing force per unit length, and the second edge attached. The influence of the choice of the kinematic hypothesis, the results of calculation of phase velocities.

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

УДК 539.3

Фазовые скорости бегущих волн в цилиндрической оболочке при конечной и бесконечной жесткости поперечного сдвига

© В.О. Каледин1, Е.А. Седова1, Ю.В. Шпакова2

1 Кемеровский государственный университет Новокузнецкий институт (филиал),

г. Кемерово, 650043, Россия 2МГТУ им. Н.Э. Баумана, Москва, 105005, Россия

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

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

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

Распространение волн в элементах конструкций типа оболочек и пластин начали рассматривать уже в XIX в. В 1889 г. Рэлей исследовал распространение волн в пластине [12], позже для нее были определены фазовые и групповые скорости. История этого вопроса и соответствующий обзор работ до 1965 г. приведены в [1, 8].

Для цилиндрических оболочек характеристическое уравнение по точной трехмерной теории при осесимметричных деформациях впервые было дано Дж. Гошем в 1923 г. [11]. В 1960 г. Микловиц опубликовал обширный перечень трудов зарубежных авторов о распространении упругих волн в стержнях, плитах, цилиндрических оболочках, полупространстве и неограниченном пространстве [13]. Исследование колебаний оболочек с заполнителем и анализ распространения свободных волн в системе оболочка - инерциальная среда началось значительно позже. Так в 1979 г. А.С. Вольмир в монографии «Оболочки в потоке жидкости и газа: Задача гидроупругости» исследует поведение деформируемых оболочек с протекающей жидкостью. Основное внимание уделено переходным процессам, связанным с резкими изменениями параметров состояния жидкости в том или ином сечении оболочки. Речь идет о цилиндрических оболочках конечной и бесконечной длины [4].

В 1992 г. А.Г. Горшков и В.И. Пожуев в монографии «Стационарные задачи динамики многослойных конструкций» дают уточ-

ненную постановку задач о распространении осесимметричных и неосесимметричных свободных волн в трехслойной цилиндрической оболочке бесконечной длины [5].

В 2000 г. в диссертации В.Г. Григорьева были рассмотрены колебания конструкции, составленной из оболочек вращения (цилиндрические, конические и сферические оболочки), полость которых заполнена жидкостью [7].

В 2004 г. в статье С.М. Белоносова рассмотрена уточненная постановка задачи о пульсирующем течении ньютоновской жидкости в упругой трубке с учетом сдвиговой податливости оболочки и раздельной постановкой условий на внешней и внутренней поверхности [2].

В диссертационной работе Шлаковой Ю.В. рассмотрены волновые процессы, возникающие в цилиндрической оболочке конечной длины при обтекании потоком жидкости [10].

В статье С.А. Бочкарева и В.П. Матвеенко исследовано динамическое поведение нагруженных оболочек вращения, содержащих неподвижную или текущую сжимаемую жидкость. Рассмотрена цилиндрическая оболочка [3].

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

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

Исходные уравнения и граничные условия. Рассмотрим движущуюся цилиндрическую оболочку. Найдем аналитическое решение в следующих допущениях:

все коэффициенты Пуассона ортотропного материала равны нулю;

главные оси анизотропии совпадают с линиями кривизны.

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

В соответствии со сделанными допущениями

r = const, Rs = да,

* v12 = 0 ^ C12 = D12 = 0 (1) Kj = 0.

Уравнения движения ортотропной оболочки вращения в криволинейной ортогональной лагранжевой системе координат (s, 0, n) получены в работе [9]:

д дР

— (rNs) н-----cos фУе +

д^ ^ де е Rs де R

1 дН r r д2и u

+—Qs+rqs+ r J pdn s = 0

д , m дЫе д cos ф ТТ .

— (rP) н---- + cos фР н-(sın фН) н-Н + sın фQе +

дs де дs Rs

t d д\ 0

+ rqe + r J pdn—- = 0, h дї2

д дQе r . r d2w

—(rQs)+— Ns -sın ф^е + rqn+r J Pdn

дs де Rs i

дї2

д дН h3 д2

— (rMs ) + ^ - cos фме - rQs + rms + rP——f

дs де 12 дї2

д ч дМе un h3 52фе

т-(rH)+- cos фн - rQ + rm+r p——i-дs де 12 дї2

= 0,

0,

0,

где и(s,0), v(s,0), w (s,0) - перемещения вдоль дуги, окружности и нормали соответственно; у s (s, 0), ^( s, 0) - углы поворота нормали; Rs, R0 - радиусы главных кривизн; ї - время; h - толщина оболочки; р - плотность материала; ф - угол, образованный нормалью к координатной поверхности и осью вращения; qs, qе, qn - компоненты поверхностной распределенной нагрузки; ms, mе - поверхностные распределенные моменты.

r = R sinф; Р = N0s -MR0 = Ne -^; Н = МЛ= Ма,;

R0 Rs

Ns = C11Ss + C12s0 + C16Ss0 + K11Ks + K12K0 + K162ks0 + ysl11 + у0l21 ; N0 = C12Ss + C22s0 + K12Ks + K22K0 + K262ks0 + уsl12 + у022 ;

Ns0 = C16S 0 + C26s s + C66s s0 + K16K0 + K26k s+K662K0s +

+ R (K16s0 + K26s s +(66s0s + D16K0 + D26k s +D66 2k s0^y sl13 + уЄ0 R0

N05 - C168e + C2685 + C66s05 + K16k0 + K26Ks + K662k0s +

N0s - C168e

66’ 0 s

Ms - K118s + K128e + K168s0 + D11Ks + D12k0 + D162ks0 + ysl14 + yel24;

M0 - K128s + K228e + K268se + D12Ks + D22k0 + D262ks0 + ysl15 + y0l25;

Mse - M0s - K168s + K268e + K668s0 + D16K5 +

+D26k0 + D662k s0 + y sl16 +yel26;

6s - Фцу s +^12ye + l118s + l12se + l138s0 + l14Ks + l15Ke + ^16Ks0 ;

Qe - ^21y s + Ф22ye + l218s + l228e + l238s0 + l24Ks + l25Ke + ^26Ks •

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

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

При использовании гипотезы Кирхгофа - Лява тензор деформаций срединной поверхности содержит шесть компонент:

(3)

8 -

T

[8s ([ e) 8e ([ e) 8s0 (^ e) Ks (s, e) k0 (^ e) K5e (s e)] ,

s s (s, 0) = ^0> + w(s e»

ds

Rs

cos Ф

00 (s, 0) =-------^ u (s, 0) + -

1 dv(s, 0) sin Ф

r d0

w(s, 0),

0 (s, 0) = I+Mi0) _ С0ЇФ. V(S, 0),

r d0

ds

к s (л Q)=d-

ds

f

dw(s, 0) u (s, 0)

Л

ds

R

s J

(4)

к0 (s, 0) =

cos Ф

dw(s, 0) u (s, 0)

Л 1 f

ds

R.

s J

Sin Ф

dv(s, 0) d2 w(s, 0)

d

d0

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

2к„(s,0) = “2f®Ş®> - £“2 v(s, 0)) +

r I ds r J

I_d_

r d0

u (s, 0) + w(s, 0)cos Ф dw(s, 0)

R

V s

ds

J

Тогда в системе уравнений движения (2) первое уравнение примет вид

dNs h d2u 0

r —- _ rph —— = 0,

ds dt2

(5)

второе и пятое уравнения выполняются тождественно; остальные преобразуем, выразив из четвертого Qs и подставив в третье:

d 2Ms

ds 2

d 2 w

r _ N0 _ rph —— = 0.

(6)

Подставим в полученные уравнения выражения для продольной

й N С du й й N C w

погонной силы Ns =Сіі , поперечной погонной силы N0 = C22 —

ds r

и продольного момента Ms = D

f d 2 w ^

11

ds2

V us J

с учетом принятых допу-

щений и кинематической гипотезы. Получим систему уравнений движения в перемещениях в случае гипотезы Кирхгофа - Лява:

rC _ roh îl - 0

" Ss2 ° St2 - 0'

^ S4w 'l ^ w

rD,

11

Ss4

S 2 w

_ C22--roh—— - 0

22 r St2

(7)

При использовании гипотезы Тимошенко деформации координатной поверхности выразим через линейные перемещения координатной поверхности: u(s, 9), v(s, 9), w(s, 9) - вдоль дуги, окружности и нормали соответственно - и углы поворота нормали ys(s, 9) и y9(s, 9) [6].

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

rC1

\rK

S 2и 1 Ss2

Ss

_ roh

S 2 и St2

- 0,

- + rK

S 2 w C

22

Ss2

w _ roh

S 2 w St2

- 0,

S2W Sw h3 S2y

rD11—2Г _ rK1y + rK1--ъ ro---2- 0.

11 Ss2 1 s 1 Ss *12 S2

(8)

Отличие системы (7) от (8) состоит в том, что (8) содержит три уравнения вместо двух в связи с большим числом степеней свободы.

Расчет фазовых скоростей. Первое уравнение в обеих системах, описывающее продольные колебания, одинаково и после сокращения на r и деления на С11 представляет собой хорошо известное волновое уравнение

S2и ph S2и 0

~SF ~ C~Иё ~ ,

описывающее продольную бегущую волну

и(s, t) - J(s) cos(01t + Фи )

(9)

(10)

с фазовой скоростью, равной скорости звука в материале оболочки:

апр -

C1l, где C11 - E1h.

oh

(11)

Волна изгиба описывается решением системы (7) или (8) в зависимости от выбранной кинематической гипотезы. В случае гипотезы

Кирхгофа - Лява найдем решение второго уравнения в виде гармонической функции w (5, t) = B (s) cos (ш 2t + фw). Из второго уравнения системы (7), разделив на w, получим следующее уравнение:

-rDn ^W)4 - С22- + ш2rph = 0. (12)

r

Для нахождения скорости волны изгиба выразим из уравнения

(12) ФW:

ф = 4 YW *

r 2 ro2ph - С22

r 2 d

11

Разделим частоту на фW, получим

аК-Л

шл/r ■ tfD— 4r2pho2 - C22

(13)

(14)

D Eh C E

где D11 = 12 ; C22 = E2h.

Выражение (14) дает точное значение фазовой скорости волны изгиба в принятых допущениях в случае гипотезы Кирхгофа - Лява.

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

д_

ds

(rQs )- N0- rPh

d 2 w dt2

0,

<

д

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

(rMs )- rQs + rP

h3 д 2v s

12 dt2

0,

(15)

с граничными условиями w (0, t) = 0, уs (0, t) = 0, Ms (L, t) = 0, Qs (L, t) = J ■ cosQt. Найдем решение системы в виде

(w ) = ( C1 )

V^ s у

V C2 у

(16)

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

Таким образом, исходная система будет 2

pho2 - % -°-K1

a

- io K

aT

- İO K

aT

P^ o2 + K1 + Dn 12 1 11

2 o

V aT У

1 1

_1 0

_1

1_ 1

0

_1

(17)

Условием существования нетривиального решения данной системы является вырожденность матрицы коэффициентов, откуда получим характеристическое уравнение относительно фазовой скорости:

-ш4 A.K.i+

aT

2

+ 1 рИш - 2

ш4D11ph - ш2 C22D11 - ш4K рИ

з У

,2 C

22

V

Г рИ3

12

1

a

(18)

12

ш2 + K

= 0.

Введем замену -1 = f и получим квадратное уравнение отно-

сительно f которое имеет два корня: f =

-b -Vb2 - 4dc

2d

и

-b + Vb2 - 4dc

f2 =------2d------' где

d = -ш4 D11K1,

4

b = o D11ph -o'

C22 D1

c = I pho2 -

C

r 2

, Ґ .3

-o4 K

ph3 12

22

ph3 2 ~

----o + K

12 1

(19)

Оценим знак каждого из корней, заметив, что d < 0. Величина 4dc будет положительной, если c < 0, т. е. необходимо чтобы вы- C2

полнилось условие o <—£ 2 w т.

phr 2

С Es

если o 2 <—22----------= o 2рит. Анализ всех возможных комбина-

phr2 Es - 5 Gsn ^

2 C22 (1)

' ' = oyL.. Условие b < 0 выполняется,

a

ций знаков b и dc показал, что в случае, если ю2 больше минималь-

,(1)

>(2)

ного из ю крит, ю крит, то f1 > 0, а f2 < 0. Если частота меньше минимальной из двух критических частот, то в случае, когда выражение b2 - 4dc положительное, оба корня отрицательные, а в случае, когда оно отрицательное - комплексные.

Таким образом, можно сделать вывод, что уравнение (18) имеет четыре корня, два из которых, в случае, если ю2 квадрат частоты

превышает минимальное из ю

(1)

ю

(2)

крит ’ крит

, чисто мнимые, а два - дей-

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

Для пары действительных фазовых скоростей константы С1 и

С2 могут быть найдены как нетривиальные решения вырожденной системы алгебраических линейных уравнений (17):

Ci = ±/yC2,

(20)

где у =

2

C22 Ю Тг 7 2

—22 + —ў K - Р^ю 2 r аТ

Кл

аТ

В случае пары мнимых скоростей аТ = ± аТ i:

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

С = ±РС2,

(21)

C22 ю

+ _2 -К + р^ю ^

где Р =

а

К1 ^

аТ

Окончательно общее решение системы (17) для случая двух действительных и двух чисто мнимых корней примет вид

Y (s, t ) = C1

f1 1 İ03S f 1 1 ms f 1 ^ ras f 11

e аТ + C2 еаТ + C3 e аТ + C4

1ІУ) V-i У) VP) v-Py

®s

âj

(22)

Найдем из граничных условий значения констант, решив следующую систему уравнений:

1

іу

-iraL

ye a

1 1

_1 О

_1

Q 0

Q 0

Q _ 1

О

_1

1

1

-іУ

iraL

iraL

P

-raL

-pe ат

1

-P

raL

-pe

(23)

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

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

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

Рис. 1. Зависимость фазовой скорости волны изгиба от жесткости поперечного сдвига

Из рисунка ясно, что независимо от частоты, начиная с определенного значения жесткости поперечного сдвига, скорости волн изгиба стремятся к предельному значению.

Найдем предельное значение фазовой скорости изгибной волны из уравнения (18) при К ^ а>;

а = -

юУг • фРЦ

2 3

ш r ph

ш4 r 2 fph3 ]

D11 1 24 J

(24)

(r(ш2 - С22 )

Сравним полученный результат со скоростью изгибной волны при использовании гипотезы Кирхгофа - Лява (14), для чего найдем разность предельной фазовой скорости в оболочке Тимошенко (24) и фазовой скорости в оболочке Кирхгофа - Лява (14);

а аК-Л = аК-Л

W - С22 )

-у±^у2 +(r2pho2 -С22)

— 1

(25)

где у

2 3

ш r ph

Прежде всего отметим, что действительная фазовая скорость получается при неотрицательном r phш - С22, т. е. когда круговая частота вынужденных колебаний больше отношения скорости продольной волны в окружном направлении к радиусу оболочки. При этом условии корень в знаменателе имеет одно действительное и одно мнимое значение. Рассматривая только действительное значение, проведем замену;

^r 2phш'

Р = ^ у

Таким образом

s(p) = Р

•J-1 + -/

(26)

(27)

где s (р) - убывающая функция с вертикальной асимптотой в нуле, т. е. s(p) в окрестности нуля стремится к бесконечности (рис. 2).

Рис. 2. Относительная разность фазовых скоростей, найденных по моделям Кирхгофа - Лява и Тимошенко

Заключение. Таким образом, при малых частотах (близких к

критической величине ш

модель Тимошенко в пределе

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

ЛИТЕРАТУРА

[1] Айнола Л. Волновые процессы деформаций упругих плит и оболочек. Изв. АН ЭССР. Сер. физ.-мат. и техн. наук, 1965, № 1, с. 3-62.

[2] Белоносов С.М. Пульсирующее осесимметричное возмущение ньютоновой жидкости в длинной цилиндрической упругой трубке. Динамика сплошной среды, Новосибирск, 2004, вып. 122, с. 3-39.

[3] Бочкарев С.А., Матвеенко В.П. Численное моделирование устойчивости нагруженных оболочек вращения при внутреннем течении жидкости. Прикладная механика и техническая физика, Новосибирск, 2008, т. 49, № 2, с. 185-195.

[4] Вольмир А.С. Оболочки в потоке жидкости и газа: Задача гидроупругости, Москва, Наука, 1979, 320 с.

[5] Горшков А.Г. Стационарные задачи динамики многослойных конструкций, Москва, Машиностроение, 1992, 224 с.

[6] Григоренко Я.М. Изотропные и анизотропные оболочки вращения переменной жесткости, Киев, Наукова думка, 1973, 228 с.

[7] Григорьев В.Г. Методология исследования динамических свойств сложных упругих и гидроупругих систем. Дис. ... д-ра техн. наук. Москва, 2000, 328 с.

[8] Нигул У.К. Волновые процессы деформации оболочек и пластин. Тр. VII Всесоюзной конференции по теории оболочек и пластин. Днепропетровск, 1969, Москва, 1970, с. 846-883.

[9] Седова Е.А. Решение связанной задачи гидроупругости. Сб. тр. IXМежрегиональной научно-практической конференции студентов и аспирантов, в 3 т. 10 апреля 2009 г., Новокузнецк, 2009, т. 1. с. 8-11.

[10] Шпакова Ю.В. Статическая прочность и колебания подкрепленных оболочек вращения из слоистых композиционных материалов. Автореф. дис. . канд. техн. наук. Томск, 2007, 16 с.

[11] Ghosh J. Longitudinal Vibrations of a Hollow Cylinder. Bull. Calc. Math. Soc, 1923-1924, vol. 14, рр. 31-40.

[12] Lord Rayleigh. On the Free Vibrations of an Infinite Plate of Homogeneous Isotropic Elastic Matter. Proc. Lond. Maht. Soc, 1889, vol. 20, рр. 225-234.

[13] Miklowitz J. Recent Developments in Elastic Wave Propagation. Appl. Mech. Revs, 1960, vol. 13 (12), рр. 865-878.

Статья поступила в редакцию 27.06.2013

Ссылку на эту статью просим оформлять следующим образом:

Каледин В.О., Седова Е.А., Шпакова Ю.В. Фазовые скорости бегущих волн в цилиндрической оболочке при конечной и бесконечной жесткости поперечного сдвига. Инженерный журнал: наука и инновации, 2013, вып. 9. URL: http://engjoumal.ru/catalog/mathmodel/technic/958.html

Каледин Валерий Олегович родился в 1955 г., окончил Харьковский авиационный институт в 1978 г. Д-р техн. наук, профессор, декан факультета информационных технологий Новокузнецкого института (филиала) федерального государственного бюджетного образовательного учреждения высшего профессионального образования «Кемеровский государственный университет». Область научных интересов: механика деформированного твердого тела и механика композиционных материалов. e-mail: vkaled@nkfi.ru

Седова Елена Александровна родилась в 1985 г., окончила Кемеровский государственный университет в 2007 г. Канд. физ.-мат. наук, доцент Новокузнецкого института (филиала) федерального государственного бюджетного образовательного учреждения высшего профессионального образования «Кемеровский государственный университет». Область научных интересов: механика деформированного твердого тела, механика жидкости и газа. e-mail: sedovaea@yandex.ru

Шпакова Юлия Владимировна родилась в 1982 г., окончила Кемеровский государственный университет в 2004 г. Канд. техн. наук, доцент МГТУ им. Н.Э. Баумана. Область научных интересов: механика деформированного твердого тела, механика композиционных материалов и конструкций. e-mail: shpakovayuliya@bmstu.ru

i Надоели баннеры? Вы всегда можете отключить рекламу.