Деформация упругой сферической оболочки, закрепленной по экватору, в потоке вязкой несжимаемой жидкости Текст научной статьи по специальности «Физика»

УДК 539.3

Вестник СПбГУ. Сер. 1. 2011. Вып. 3

ДЕФОРМАЦИЯ УПРУГОЙ СФЕРИЧЕСКОЙ ОБОЛОЧКИ, ЗАКРЕПЛЕННОЙ ПО ЭКВАТОРУ,

В ПОТОКЕ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ*

H. В. Наумова1, Б. А. Ершов2, Д. Н. Иванов3

I. С.-Петербургский государственный университет, канд. физ.-мат. наук, доцент, nat_n75@mail.ru

2. С.-Петербургский государственный университет, д-р физ.-мат. наук, профессор, bae_1925@mail.ru

3. С.-Петербургский государственный университет, канд. физ.-мат. наук, ст. науч. сотр., denisiv3@rambler.ru

Введение. Задачу об обтекании абсолютно твердого шара вязкой несжимаемой жидкостью впервые решил Стокс в 1851 г. [1], для этого он использовал приближенный метод, который заключается в пренебрежении в основных уравнениях движения инерционными членами и внешними силами. В работах [2, 3] приведено решение такой задачи для случая, когда сфера является абсолютно жесткой, т. е. деформации поверхности не рассматриваются, и движение жидкости осуществляется при малых числах Рейнольдса. Большая часть опубликованных работ (см., например, [4, 5]) посвящена исследованию наиболее простого случая деформации сферической оболочки под действием равномерного внешнего давления.

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

Постановка задачи. Рассмотрим обтекание сферы радиуса а, центр которой находится в начале координат, потоком вязкой жидкости. Направление потока и сфера показаны на рис. 1. Для определенности будем считать, что на экваторе оболочка жестко закреплена.

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

Уравнения равновесия элемента оболочки. В этом разделе используются обозначения из монографий [6, 7], а именно: и —тангенциальная компонента перемещения, т — нормальная компонента перемещения (прогиб), Н — толщина сферической оболочки, ^ = Н2/12 — малый параметр, Т., Т2 —тангенциальные усилия, N1 —перерезывающее усилие, М1, М2 — моменты, £1, £2 — компоненты деформации, р 1 и рз — тангенциальная и нормальная компоненты поверхностной нагрузки. Если срединная поверхность оболочки отнесена к некоторой конкретной криволинейной системе

* Работа выполнена при финансовой поддержке РФФИ (грант №07-01-00250).

© Н. В. Наумова, Б. А. Ершов, Д. Н. Иванов, 2011

Рис. 1. Сферическая оболочка, обтекаемая потоком вязкой несжимаемой жидкости.

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

cWi dMx

- (Ti + T2)+ N1 ctg в = рз, + ctg в (Mi - M2) = N1.

Моменты Mi и M2 связаны с углом деформации 71 соотношениями Mi = М4 (■^ТГ + v ctg 07Л , М2 = м4 (у^г + ctg 071

¿в

de

dw

где 7i = —— + и. ¿в

(1) (2)

(3)

(4)

(5)

Тангенциальные усилия Т и Т2 связаны с компонентами деформации Е\, £2 соотношениями:

Ti = £ 1 + V£2, T2 = £2 + V£i, du

где si = — + го, £2 = и ctg 0 + го. ¿в

(6) (7)

Тангенциальная и нормальная компоненты поверхностной нагрузки р\ и рз получены в работах [2, 3] и имеют следующий вид:

pi = Ä'sinö, рз = —К cos в, —.

2 Ena

(8)

где Д — коэффициент динамической вязкости, V — коэффициент Пуассона, Е — модуль Юнга. Безразмерные величины связаны с размерными следующими соотношениями:

1-v2 - -

{Ti,N1}=-=—{Ti,N1}, Mi

Ena

1 - v2 ^ 1 - _ .

———тгMi, {/?, и, го} = —{/г, и, го}, i = 1,2. Ena2 a

Рассмотрим решения на промежутке 0 ^ в ^ п/2, а для построения полного решения будем продолжать значения найденных функций с учетом их четности (или нечетности) на оставшуюся часть. Обоснованность этого факта будет рассмотрена позже.

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

Для определенности будем считать, что на экваторе оболочка закреплена таким образом, что

х(1)(п/2) = и(2)(п/2) = 0, У(1)(п/2) = ш(2)(п/2) =0,

М(1)(п/2) = (п/2) = 0, ^ (п/2) = (п/2).

г(2)

„(2),

(9) (10) (11)

Кроме того, из физических соображений осесимметрии рассматриваемой задачи (предполагается, что нет особенностей в полюсах) получаем

(1)

(0)

(2)

(п) = 0,

7 (1)(0) = 7(2)(п) =0,

М(1)(0) = (п) =0.

(2)

(12)

(13)

(14)

Пренебрегая членами второго порядка малости по относительной толщине Н, т.е., полагая ц = 0 в (4), получаем безмоментную систему (1)—(2), из которой можно найти функции Т (в) и Т2 (в):

ТЛв) = -^-(-Ксовв + Сг),

БШ в

Т2 (в) = -Тг -рз = —^-(-Ксовв + С^ + Ксовв. бш2 в

(15)

(16)

Выразив из (6) компоненты деформации £1 и £2 через усилия Т1 (в) и Т2 (в), получим систему двух уравнений для определения и и ш. В частности, для определения перемещения и получим следующее уравнение:

¿и 1

2

1 - V V Бт2 в

(-К СОБ в + С1) - К СОБ в .

(17)

Решение линейного дифференциального уравнения (17) будет иметь вид

и(к)(в)

эт в

1 - V

К в - К 1п яш в + С1

(к)

1п

в tg -й 2

СОБ в

+ С.

(к)

(18)

Неизвестные постоянные С(к) и С2к) определяются из граничных условий (9), (12). Определим тангенциальную компоненту перемещения:

,(1)(в) =

К в т в

1 - V

— СОБ в

1 + соэ в

— 1п [1 + СОБ (

0 < в <

(19)

п

Нормальные перемещения ш точек срединной поверхности находятся с помощью формул (7) и имеют следующий вид:

(0) =

K

1 - v2

[1 + v - cos 0 [2 + v + (1 + v) ln [1 + cos <

0 < 0 <

2

(20)

Отметим, что ш(к), найденное по формуле (20), не удовлетворяет условию (10), т.к. порядки нулевого и первого приближений для функциий ш(к) одинаковые. В связи с этим обстоятельством данная задача может быть корректно разрешена только с учетом краевого эффекта [8] вблизи закрепленного экватора.

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

У = (и, ш, Т1, 7Ь М1, М1)т; тогда систему (1)—(5) можно представить в матричном виде:

где P = (pi, рз, 0, 0, 0, 0)T а матрицы А(0,м) и B(0, м) имеют вид

А(0,м) =

(21)

0(1 - v) 0 1 0 0 0

(1 + v) 00 0 0 -1

0 00 ctg 0/(1 - v) 1 0

1 00 0 0 0

0 00 м4 0 0

0 10 0 0 0

/611 0 0 0 0 1 \

621 2(1 + v) 0 0 0 626

0 0 0 634 0 -1

641 (1 + v) -1 0 0 0

0 0 0 654 -1 0

-1 0 0 1 0 0

B(0, м)

Ниже выписаны значения некоторых коэффициентов матрицы В(0,м): bii = - ctg2 0(1 - v), 621 = (1+ v) ctg 0, 626 = - ctg 0,

634 = - ctg2 0^4(1 - v), 641 = v ctg 0,

654 = v^4 ctg 0v.

Приближенное решение краевой задачи (21), (9)—(13) можно представить в виде (см., например, [8])

У = Уа + Уь. (22)

Вектор Уа описывает главное напряженно-деформированное состояние и является решением безмоментной системы, получающейся из (21) при ^ = 0. В нашем случае

V ( (к) (к) гр(к) (к)

Уа = («4 % % Т1а , 71а

Ta, Yi?, 0, , где wifc) и wifc) находятся по формулам (19), (20),

а Ta и 7ia могут быть вычислены по формулам (6), (7) и (5).

п

0 0 0 0 0

0 м0 0 0 0 0

0 0 м1 0 0 0

0 0 0 м-1 0 0

0 0 0 0 м2 0

V0 0 0 0 0 м1

Вектор Yь описывает простой краевой эффект [8] вблизи закрепленного экватора (в = п/2). Следуя теории краевого эффекта, вектор Yb будем искать в виде

Yb = Y(0)exp - [ qdt , (23)

W J

где "Y = (T, w, T\, 7i, M\, N\)T. После подстановки (22) в систему (21) для нахождения компонент вектора Y получим новую систему:

Л(в, м)^® + (в(в, м) + 1А(в, ^ Y (в) = 0. (24)

Вектор YT (в) можно представить в виде

Y(в) = D • Y(в), (25)

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

D=

Пренебрегая более высокими степенями малого параметра д, т. е. отбрасывая слагаемое A(в,д)дY(в)/дв в системе (24) и учитывая (25), получаем алгебраическую систему уравнений

Т + ^ в(1 - V)« = 0, т(1 + V) + 7ц = 0, (26)

N + 7(1 - V) = 0, М1 - ^71 = 0, (27)

N - цМ1 = 0, 71 + цт = 0. (28)

В силу однородности системы (26)—(28) ее определитель Д = ц4 + (1 - V2) должен равняться нулю, следовательно,

q = —-=— (1 — г/2)1^4 , где г — мнимая единица. (29)

а/2

Аналитически разрешая систему (26)—(28), находим вектор Y(в), а затем Y (в):

/ 2 3 4 \ Т

= • (30)

Таким образом, составляющая тъ функции прогиба (вторая компонента вектора Yb) определяется по следующей формуле:

т=~ и' ехр ^ У4 ^ I •

В частности, с учетом (29) перепишем выражение (31):

шь = -

1 + V

¿=1

• ехр -М

J ^ ^ I .

о /

Неизвестные постоянные О находятся путем подстановки ш = ша + в граничные условия (10). Процедура подстановки и определение значений постоянных О осуществлялась в пакете прикладных программ МаШеша^са 5.1. Окончательно, для нахождения функции шь получим

шь (в) = -

К

1 - V

exp

СОБ

(33)

1

1

Анализ полученных результатов. Сравнение численных и асимптотических расчетов, полученных по формулам (19), (20), (22), (33), показывает достаточно хорошее совпадение результатов. Максимальная относительная погрешность составляет 20%. Численные расчеты проведены с помощью метода конечных элементов (МКЭ) в пакете А№УЯ 5.6.

В качестве примера была рассмотрена сферическая оболочка радиуса а = 0.2 м, толщиной Н = 0.002 м с коэффициентом Пуассона V = 0.3, модулем Юнга Е = 2.07 • 1011 Н/м2. Скорость потока и = 0.01 м/с, коэффициент динамической вязкости М =1 Пас.

На рис. 2 представлено меридианальное сечение поверхности сферической оболочки. Сплошной линией показана деформированная поверхность, штриховой линией — поверхность до деформации.

Рис. 2. Меридианальное сечение поверхности сферической оболочки.

Рис. 3. График функции и (в). Сплошная линия — асимптотические формулы, пунктирная линия — МКЭ.

Графики функций и (в) и ш (в) представлены на рис. 3 и рис. 4, соответственно. В направлении горизонтальной оси откладывается угол в, по вертикали — значения функций и (в) и ш (в) в метрах. Аналогичную картину для деформаций сферической оболочки, находящейся в потоке вязкой несжимаемой жидкости, можно наблюдать и в пакете А№УЯ. На графиках, представленных на рис. 3 и рис. 4, сплошная линия

w(6)

0.002

-0.00]

0.00]

в

Рис. 4• График функции ю (в). Сплошная линия — асимптотические формулы, пунктирная линия — МКЭ.

-0.002

соответствует асимптотическим расчетам по формулам (19), (20), (22), (33), а пунктирная линия — численным расчетам, полученных методом конечных элементов.

Как нетрудно заметить, w является нечетной функцией, а и — четной функцией относительно прямой 0 = п/2. В нашем случае условия нечетности функции w, w(0) = — w(-0), и четности функции и, и(0) = и(-0), выполняются для симметричных (относительно экватора) точек верхнего (северного) и нижнего (южного) полушарий. Это свойство функций использовалось при аналитическом интегрировании системы (1)—(5). Так, например, для четных функций в данной системе f (0) = — f (п — 0), для нечетных f (0) = f (п — 0). После введения новой переменной £ = п — 0 для точек южного полушария получается система, аналогичная (1)—(5), относительно £. Поэтому вместо меридиана можно рассматривать его половину, а для построения полного решения продолжать значения найденных функций с учетом их четности (нечетности) на оставшуюся часть.

Литература

1. Stokes G. G. On the effect of the internal friction of fluid on the motion of pendulums // Math. and Phys. Papers. Vol. 3. 1851.

2. Кочин Н. Е., Кибель И. А., Розе Н. Е. Теоретическая гидромеханика. Т. 2. М., 1963.

3. Ландау Л. Д., Лифшиц Е.М. Гидродинамика. Т. 4. М., 1988.

4. Гурьянов Н. Г. Сферическая оболочка, находящаяся под действием нагрузки, равномерно распределенной по площадке // Исслед. по теор. пластин и оболочек. Вып. 5. Казань: Изд-во Казанского ун-та, 1967. C. 136-147.

5. Слепов Б. И. Устойчивость незамкнутой сферической оболочки под действием равномерного внешнего давления // Исслед. по теор. пластин и оболочек. Вып. 4. Казань: Изд-во Казанского ун-та, 1966. C. 139-151.

6. Гольденвейзер А. Л., Лидский В. Б., Товстик П.Е. Свободные колебания тонких упругих оболочек. М., 1979.

7. Филиппов С. Б. Теория сопряженных и подкрепленных оболочек. СПб., 1999.

8. Новожилов В. В. Теория тонких оболочек. СПб.: Изд-во С.-Петерб. ун-та, 2010. 380 с.

Статья поступила в редакцию 22 апреля 2011 г.